### Approach and split the tensor component

I was trying to split the tensor into 4 components since it is the second-order tensor. I have seen other posts like

```
W= TensorFunctionSpace(mesh,'Lagrange',2)
sigma_w=project(E/(1+nu)*(sym(grad(u))+tr(grad(u))*nu/(1-2*nu)*Identity(d)),W)
```

or

```
W= VectorFunctionSpace(mesh,'Lagrange',2,dim=2)
sigma_w=project(E/(1+nu)*(sym(grad(u))+tr(grad(u))*nu/(1-2*nu)*Identity(d)),W)
```

The first one gave the error

Traceback (most recent call last):

File "/home/qiyaopeng/Desktop/PhD/StudyNotes/Literature/Jiao's paper/FEniCS Demo/elasticity_Smalltri.py", line 135, in <module>

sigma_w=project(E/(1+nu)*(sym(grad(u))+tr(grad(u))*nu/(1-2*nu)*Identity(d)),W)

File "/usr/lib/python3/dist-packages/dolfin/fem/projection.py", line 142, in project

form_compiler_parameters=form_compiler_parameters)

File "/usr/lib/python3/dist-packages/dolfin/fem/assembling.py", line 375, in assemble_system

assembler = cpp.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)

File "/usr/lib/python3/dist-packages/dolfin/cpp/fem.py", line 2698, in __init__

_fem.SystemAssembler_swiginit(self, _fem.new_SystemAssembler(a, L, bcs))

RuntimeError:

*** -------------------------------------------------------------------------

*** DOLFIN encountered an error. If you are not able to resolve this issue

*** using the information listed below, you can ask for help at

***

*** fenics-support@googlegroups.com

***

*** Remember to include the error message listed below and, if possible,

*** include a *minimal* running example to reproduce the error.

***

*** -------------------------------------------------------------------------

*** Error: Unable to assemble system.

*** Reason: expected a linear form for L.

*** Where: This error was encountered inside SystemAssembler.cpp.

*** Process: 0

***

*** DOLFIN version: 2017.2.0

*** Git changeset: unknown

*** -------------------------------------------------------------------------

and the second one gave the error of mismatched shapes. Then I had no idea how to fix it.

Here comes my code

```
from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import numpy as np
x0=20
y0=15
nx=40
ny=20
mesh=RectangleMesh(Point(-x0,-y0),Point(x0,y0),nx,ny)
### FEM to solve the PDE
### parameters value
E=1
nu=0.49
V = VectorFunctionSpace(mesh, 'P', 2)
# def epsilon(u):
# return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# #return sym(nabla_grad(u))
# def sigma(u):
# return E/(1+nu)*(epsilon(u)+tr(epsilon(u))*Identity(d)*nu/(1-2*nu))
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension() # space dimension
W= TensorFunctionSpace(mesh,'Lagrange',2)
# W=VectorFuctionSpace(mesh,'Lagrange',2,dim=2)
sigma_w=project(E/(1+nu)*(sym(grad(u))+tr(grad(u))*nu/(1-2*nu)*Identity(d)),W)
```

Thanks in advance for the kind help!

Best wishes,

Alice

### 1 Answer

`u`

as a `Function`

instead of a `TrialFunction`

and the code works.
`TrialFunction`

is only used as argument in a bilinear form. As such it is more of an abstract UFL object rather than something filled with content like nodal values. Writing

`u = Function(V)`

initializes zero values everywhere associated with the `Function`

object which then can actually be used in the computation.
```
from __future__ import print_function
from fenics import *
from mshr import *
from dolfin import *
import numpy as np
x0=20
y0=15
nx=40
ny=20
mesh=RectangleMesh(Point(-x0,-y0),Point(x0,y0),nx,ny)
subdomain_vert=np.array([[0,0,1],[0,1.5,1.5]])
center=np.mean(subdomain_vert,axis=1)
### FEM to solve the PDE
### parameters value
E=1
nu=0.49
P=2.08
kappa=3
V = VectorFunctionSpace(mesh, 'P', 1)
# def epsilon(u):
# return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# #return sym(nabla_grad(u))
# def sigma(u):
# return E/(1+nu)*(epsilon(u)+tr(epsilon(u))*Identity(d)*nu/(1-2*nu))
u = Function(V)
v = TestFunction(V)
d = v.geometric_dimension() # space dimension
W= TensorFunctionSpace(mesh,'CG',1)
sigma_w=project(E/(1+nu)*(sym(grad(u))+tr(grad(u))*nu/(1-2*nu)*Identity(d)),W)
[s11,s12,s21,s22]= sigma_w.split(True)
[u1,u2]=u.split(True)
a1=kappa*u1*ds+(s11*grad(v)[0,0]+s12*grad(v)[0,1])*dx
a2=kappa*u2*ds+(s12*grad(v)[1,0]+s22*grad(v)[1,1])*dx
L1=Constant(0)*v[0]*dx
L2=Constant(0)*v[1]*dx
###########################################
###the following code doesnt work at all
###########################################
#A1=assemble(a1) #this one already cannot work
#b=assemble(L1)
#delta=PointSource( V, Point ( 0.0, 0.0 ), 100.0)
#delta.apply(b)
##or even
#solve(a1==L1,u1) #cannot work either
```

The error was in the following:

*** Error: Unable to perform just-in-time compilation of form.

*** Reason: ffc.jit failed with message:

Traceback (most recent call last):

File "/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py", line 142, in jit

result = ffc.jit(ufl_object, parameters=p)

File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 218, in jit

module = jit_build(ufl_object, module_name, parameters)

File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 134, in jit_build

generate=jit_generate)

File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 167, in jit

header, source, dependencies = generate(jitable, name, signature, params["generator"])

File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 67, in jit_generate

prefix=module_name, parameters=parameters, jit=True)

File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 143, in compile_form

prefix, parameters, jit)

File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects

analysis = analyze_ufl_objects(ufl_objects, kind, parameters)

File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 94, in analyze_ufl_objects

for form in forms)

File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 94, in <genexpr>

for form in forms)

File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 178, in _analyze_form

do_apply_restrictions=True)

File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 388, in compute_form_data

check_form_arity(preprocessed_form, self.original_form.arguments()) # Currently testing how fast this is

File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 152, in check_form_arity

check_integrand_arity(itg.integrand(), arguments)

File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 147, in check_integrand_arity

raise ArityMismatch("Integrand arguments {0} differ from form arguments {1}.".format(args, arguments))

ufl.algorithms.check_arities.ArityMismatch: Integrand arguments () differ from form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 0, None),).

.

*** Where: This error was encountered inside jit.py.

*** Process: 0

***

*** DOLFIN version: 2017.2.0

*** Git changeset: unknown

*** -------------------------------------------------------------------------

What else can I do? Thanks in advance

`solve`

instruction.If you project an object like your stress onto a function space you lose its generality as an abstract object. In your case,

`s11, s12`

... are simply zero everywhere since you use the zero initialized values of `u`

for calculation.Why don't you use a vectorial notation and write sigma as combination of UFL objects? If I translated your forms correctly you have

\[

a(u,v) =\int_{\Omega}\mathbf\sigma \cdot\!\cdot \operatorname{grad} v\,\mathrm dx+\int_{\partial\Omega}\kappa \mathbf u\,\mathrm d s

\]

You can write this as (just an excerpt from your code)

```
u = TrialFunction(V)
v = TestFunction(V)
epsilon = 0.5*(nabla_grad(u) + nabla_grad(u).T)
sigma = E/(1+nu)*(epsilon+tr(epsilon)*Identity(2)*nu/(1-2*nu))
a = kappa*u*ds + inner(sigma, grad(v))*dx
```

Then analogously define your linear form in a vectorial manner. For solution you define a new function in which the solution values get stored:

```
u = Function(V)
solve(a==L,u)
```

After this, the function `u`

contains the solution. If you're now interested in values of stress, you can project.

As a side note, you can also directly access the components of a tensor. If sigma is defined as above, you can for example write `sigma[0,0]`

do get the \(\sigma_{11}\) component.

$$ $-\nabla\cdot\sigma=\sum_{\left\{i,j\right\}}P\textbf{n(x)}\delta\left(x-x_{\left\{i,j\right\}}\right)\bigtriangleup\Gamma_{\left\{i,j\right\}}$−∇·

`σ`=∑

_{{i,j}}

`P`

**n(x)**

`δ`(

`x`−

`x`

_{{i,j}})△Γ

_{{i,j}}

with boundary condintion

$\frac{\partial\sigma}{\partial\textbf{n}}+\kappa\textbf{u}=0$∂

`σ`∂

**n**+

`κ`

**u**=0

Here,

**n(x)**is the inward direction (vector).

Since the PointSource magnitude can only be constant, I have no idea how to assemble with the inward direction and I decided to piecewise the equation in x direction and y direction respectively.

Therefore, all the problem comes from delta function.

Thanks a lot for your kind help and explanation.