# Solver effor for linear elastic problem

**URL:** <https://fenicsproject.discourse.group/t/solver-effor-for-linear-elastic-problem/13640>\
**Category:** General\
**Created:** [February 16, 2024, 7:12pm UTC](https://fenicsproject.discourse.group/t/solver-effor-for-linear-elastic-problem/13640 "2024-02-16T19:12:53Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Devin](https://avatars.discourse-cdn.com/v4/letter/d/278dde/32.png) [@Devin](https://fenicsproject.discourse.group/u/Devin)\
**Post date:** [February 16, 2024, 7:12pm UTC](https://fenicsproject.discourse.group/t/solver-effor-for-linear-elastic-problem/13640/1 "2024-02-16T19:12:53Z")

</div>

I was running a linear elastic problem. The python script are shown as follows:

```auto
from __future__ import print_function
from dolfin import *
from ufl import nabla_div
from mpi4py import MPI
import numpy as np
import meshio
import gmsh
import sys

path = 'sample0_300_8973410.stl'

#Scaled variables

E = 2.5e11
nv = 0.27

lambda_ = E * nv / ((1 + nv) * (1 - 2 * nv))
mu = 0.5 * E / (1 + nv)

comm = MPI.COMM_WORLD

if comm.rank == 0:

    gmsh.initialize()
    gmsh.merge(path)

    s = gmsh.model.getEntities(2)
    l = gmsh.model.geo.addSurfaceLoop([e[1] for e in s])
    V0 = gmsh.model.geo.addVolume([l])

    gmsh.model.addPhysicalGroup(3, [V0], 1)
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(3)

    #obtain the model dimension
    dimensionBox = gmsh.model.getBoundingBox(-1, -1)

    gmsh.write("sample0_100.msh")
    gmsh.finalize()

    def create_mesh(mesh, cell_type, prune_z=False):
        cells = mesh.get_cells_type(cell_type)
        cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
        points = mesh.points[:, :2] if prune_z else mesh.points
        out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
        return out_mesh

    mesh_from_file = meshio.read("sample0_100.msh")

    tetra_mesh = create_mesh(mesh_from_file, "tetra", prune_z=False)
    meshio.write("tetra_mesh.xdmf", tetra_mesh)

    print(mesh_from_file)

else:
    dimensionBox = None
comm.Barrier()

dimensionBox = comm.bcast(dimensionBox, root=0)

xmin = dimensionBox[0]
ymin = dimensionBox[1]
zmin = dimensionBox[2]
xmax = dimensionBox[3]
ymax = dimensionBox[4]
zmax = dimensionBox[5]

mesh = Mesh()
with XDMFFile("tetra_mesh.xdmf") as infile:
     infile.read(mesh)

V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1.5

def boundary_Left(x, on_boundary):
    return on_boundary and near(x[0], xmin, tol)

def boundary_Right(x, on_boundary):
    return on_boundary and near(x[0], xmax, tol)

def boundary_Front(x, on_boundary):
    return on_boundary and near(x[1], ymin, tol)

def boundary_Back(x, on_boundary):
    return on_boundary and near(x[1], ymax, tol)

def boundary_Top(x, on_boundary):
    return on_boundary and near(x[2], zmax, tol)

def boundary_Bottom(x, on_boundary):
    return on_boundary and near(x[2], zmin, tol)

applied_strain = 1.0e-4
applied_disp = applied_strain * (xmax - xmin)

bc_l = DirichletBC(V.sub(0), Constant(0), boundary_Left)
bc_r = DirichletBC(V.sub(0), Constant(applied_disp), boundary_Right)
bc_f = DirichletBC(V.sub(1), Constant(0), boundary_Front)
bc_ba = DirichletBC(V.sub(1), Constant(0), boundary_Back)
bc_bo = DirichletBC(V.sub(2), Constant(0), boundary_Bottom)
bc_t = DirichletBC(V.sub(2), Constant(0), boundary_Top)

bcs = [bc_l, bc_r, bc_f, bc_ba, bc_bo, bc_t]

# Define strain and stress
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))

a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs, solver_parameters = {'linear_solver':'mumps'})

```

When I run the job, I got the error message:

Error: Unable to solve linear system using PETSc Krylov solver.  
\*\*\* Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED\_PC\_FAILED, residual norm ||r|| = 0.000000e+00).  
\*\*\* Where: This error was encountered inside PETScKrylovSolver.cpp.

Does anyone have idea about what was wrong in my model? Thanks.

---

<div class="post-metadata">

**Author:** ![Devin](https://avatars.discourse-cdn.com/v4/letter/d/278dde/32.png) [@Devin](https://fenicsproject.discourse.group/u/Devin)\
**Post date:** [February 16, 2024, 7:13pm UTC](https://fenicsproject.discourse.group/t/solver-effor-for-linear-elastic-problem/13640/2 "2024-02-16T19:13:51Z")

</div>

The input file can be found through the link:[https://1drv.ms/u/s!AjzJPj4bEblbgR-bVpb2CdQxUOu6?e=XrNhGB](https://1drv.ms/u/s!AjzJPj4bEblbgR-bVpb2CdQxUOu6?e=XrNhGB)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 20, 2024, 8:01pm UTC](https://fenicsproject.discourse.group/t/solver-effor-for-linear-elastic-problem/13640/3 "2024-03-20T20:01:15Z")

</div>

As you are providing a very big model (4 million nodes), which I cannot even perform the meshing part of on my 16 GB laptop. You also havent told us how many processes you are using, I find it hard to give you guidance.