Approach and split the tensor component


184
views
0
3 months ago by
Hello, everyone!

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
Community: FEniCS Project

1 Answer


3
9 weeks ago by
Declare u as a Function instead of a TrialFunction and the code works.
that does work! Thanks a lot but could you please explain further why I should use Function rather than TrialFunction?
written 9 weeks ago by Alice Peng  
As far as I know, a 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.
written 9 weeks ago by klunkean  
Then does that mean, I cannot write the bilinear form with u and sigma now? Because my complete codes are 
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

written 9 weeks ago by Alice Peng  
I think your approach does not really fit into the FEniCS line of thought. Preferably you would work with ufl objects up to the 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.

written 9 weeks ago by klunkean  
Okay, thanks for the answer. Now all my problem gets back to the PointSource again... My original PDE is 
$$      $-\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}Pn(x)δ(xx{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.
written 9 weeks ago by Alice Peng  
Maybe open a new question thread for this, since it really doesn't fit this question anymore (Approach and split the tensor component)
written 9 weeks ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »