Linear elasticity using a disk

variational formulation
1

Hai,

I was running a linear elastic problem but encountered an error

raceback (most recent call last):
  File "/lfs1/usrscratch/phd/am19d027/job527879/test3d.py", line 76, in <module>
    sigmaT = project(sigma(u),WV,solver_type="cg")
  File "/lfs/sware/anaconda3_2019/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/projection.py", line 132, in project
    A, b = assemble_system(a, L, bcs=bcs,
  File "/lfs/sware/anaconda3_2019/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/assembling.py", line 382, in assemble_system
    assembler.assemble(A_tensor, b_tensor)
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 successfully call PETSc function 'MatXIJSetPreallocation'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1608660238187/work/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  640ef70247ae212e16ae9f24fc9bc603b506f78a

The code was

mesh = Mesh('disK3d.xml')
V = VectorFunctionSpace(mesh, "CG", 2)
WV=TensorFunctionSpace(mesh, "DG", 1)

u = Function(V)
solve(a == L, u, bc,solver_parameters={'linear_solver':'cg'})

sigmaT = project(sigma(u),WV,solver_type="cg")
2

Have you tried to run the same problem on a built in mesh (such as a unit square)?
Please also note that to be able to help you we need a minimal code that reproduces the error message; see for instance: Read before posting: How do I get my question answered? - #4 by nate

3

Hai,

I have not tried on a square mesh. But, when I am using a coarse mesh, the code is working fine. It is not working for finer mesh. My fine mesh was having 17e6 elements.
If in that case which solver can be used?

Also, If want to run parallel, how can I modify the solver. I need to extract the say if find the the sigma 1, at location say center (0,0) ,
sigma_1 = 0.5*(sigmaT[0,0] + sigmaT[1,1]) + sqrt( (0.5*(sigmaT[0,0] - sigmaT[1,1]))**2 + (sigmaT[0,1])**2 )
Thanks

4

You need to

  1. specify what variational problem you are solving
  2. provide a code that can reproduce the issue

See for instance: MPI - Gather and distribute a function to all processes (That I can evaluate at a point) - #4 by dokken

5

I was solving a disk under diametral compression.

mesh = Mesh('disk3dnew.xml')
V = VectorFunctionSpace(mesh, "CG", 2)
WV=TensorFunctionSpace(mesh, "DG",0)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[1]<-(Radius-(5E-3)) and on_boundary  
    
class top(SubDomain):
    def inside(self, x, on_boundary):
        return  x[1]>(Radius-(5E-3))  and on_boundary 

tp=top()
bt=bottom()
boundary = MeshFunction("size_t", mesh,mesh.topology().dim()-1)
boundary.set_all(0)
tp.mark(boundary, 1)
bt.mark(boundary, 2)

bc = DirichletBC(V, Constant((0.0,0.0,0.0)), boundary,2)
ds = Measure("ds",domain=mesh, subdomain_data=boundary)
n = FacetNormal(mesh)

Area= assemble(1*ds(1))
Traction= (-500)/Area

def sigma(u):
    return 2.0*mu*sym(nabla_grad(u)) + lmbda*tr(sym(nabla_grad(u)))*Identity(len(u))

u = TrialFunction(V)
v = TestFunction(V)

a = inner(nabla_grad(v),sigma(u))*dx
L = dot(v, Traction*n)*ds(1)


u = Function(V)
solve(a == L, u, bc,solver_parameters={'linear_solver':'cg'})

plot_ = File ("./stress.pvd")
sigmaT = project(sigma(u),WV,solver_type="cg")
plot_<<sigmaT
6

First of all, in your code, your DirichletBC is 3D, i.e.

does that mean that your mesh is three-dimensional (i.e. a manifold)?

If not, you should check the geometrical dimension of your mesh, i.e.

print(mesh.geometry().dim())

and make sure it is 2, not 3

As you have not supplied the mesh, I cannot reproduce your problem.
Here is a minimal working code using a built in mesh:

from dolfin import *
N = 100
mesh = UnitSquareMesh(N, N)
mu = 0.4
lmbda = 1e5
V = VectorFunctionSpace(mesh, "CG", 2)
WV = TensorFunctionSpace(mesh, "DG", 0)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary


class top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1) and on_boundary


tp = top()
bt = bottom()
boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary.set_all(0)
tp.mark(boundary, 1)
bt.mark(boundary, 2)

bc = DirichletBC(V, Constant((0.0, 0.0)), boundary, 2)
ds = Measure("ds", domain=mesh, subdomain_data=boundary)
n = FacetNormal(mesh)

Area = assemble(1*ds(1))
Traction = (-300)/Area


def sigma(u):
    return 2.0*mu*sym(nabla_grad(u)) + lmbda*tr(sym(nabla_grad(u)))*Identity(len(u))


u = TrialFunction(V)
v = TestFunction(V)

a = inner(nabla_grad(v), sigma(u))*dx
L = dot(v, Traction*n)*ds(1)


u = Function(V)
solve(a == L, u, bc, solver_parameters={
      'linear_solver': 'cg', "preconditioner": "amg"})

plot_ = File("./stress.pvd")
sigmaT = project(sigma(u), WV, solver_type="cg")
plot_ << sigmaT
7

Ya, it was 3d. It is a simple disk (or cylinder), say 20 mm dia with 3 mm thick

8

Start by making the problem reproducable by using s unit cube as opposed to a box.

9

Thanks,
May I know why did you add a precondition to the solver parameters

solve(a == L, u, bc, solver_parameters={
      'linear_solver': 'cg', "preconditioner": "amg"})
10

Because any variational problem should be properly preconditioned. Using default settings/the linear algebra solver as a black box is never recommended.

For linear elasticity, I usually use the following:

You are also advices to set the nullspace of your problem; Bitbucket

1 Like
11

Ok, I got it. Thanks Dokken