Mesh size error to calculate stress in 3d elasticity

1

In a 3d linear elasticity model i hope to calculate Displacement and Stress tensor and i had a error for the model mesh size. I found the solution in UMFPACK error: out of memory despite system having free memory - #12 by ncarmona, but it only works to calculate the solution for u, when i try to calculate the Stress Tensor i have an error again. The code is:

from dolfin import*
from mshr import*
from fenics import*
from ufl import nabla_div
import matplotlib.pyplot as plt
import numpy as np

mesh = Mesh("Caseron_actual.xml")
#plot(mesh);
#Scaled Variables
E = Constant(20e3) #Módulo de Young
nu = Constant(0.3) #Radio de Poisson
mu = E/(2*(1+nu)) #Módulo de cizalle
lambda_ = E*nu/((1+nu)*(1-2*nu)) #Lamé
rho = Constant(0.027) #densidad

# Strain function
def epsilon(u):
    return sym(nabla_grad(u))

# Stress function
def sigma(u):
    return lambda_*nabla_div(u)*Identity(3) + 2*mu*epsilon(u)
# lambda is a reserved python keyword, naming convention recommends using a single trailing underscore for such cases

#Load
g_z = 0.0
b_z = -rho
g = Constant((0.0,g_z,0.0))
b = Constant((0.0,0.0,b_z))

#Function Spaces
V = VectorFunctionSpace(mesh, "P", 2)
u_tr = TrialFunction(V)
u_test = TestFunction(V)

a = inner(sigma(u_tr),epsilon(u_test))*dx #bilineal form
l = dot(b,u_test)*dx + dot(g, u_test)*ds #lineal form

x_max = 1216.0
y_max = 929.
tol=1e-8
def bottom(x, on_boundary):
    return (on_boundary and near(x[2],0,tol))
bc1 =DirichletBC(V.sub(2), Constant(0.),bottom)

def left(x, on_boundary):
    return (on_boundary and near(x[0], 0,tol))
bc2 =DirichletBC(V.sub(0), Constant(0.),left)

def right(x, on_boundary):
    return (on_boundary and near(x[0], x_max,tol))
bc3 =DirichletBC(V.sub(0), Constant(0.),right)

def north(x, on_boundary):
    return (on_boundary and near(x[1], y_max,tol))
bc4 =DirichletBC(V.sub(1), Constant(0.),north)

def south(x, on_boundary):
    return (on_boundary and near(x[1], 0,tol))
bc5 =DirichletBC(V.sub(1), Constant(0.),south)
bc=[bc1,bc2,bc3,bc4,bc5]

u = Function(V)
#problem = LinearVariationalProblem(a, l, u, bc)
#solver = LinearVariationalSolver(problem)
#solver.solve();
solve(a==l, u, bc, solver_parameters={'linear_solver' : 'mumps'})

#Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
#sig = Function(Vsig, name="Cauchy_Stress")
#sig.assign(project(-sigma(u), Vsig));

file_results = XDMFFile("Verde_CA.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
#file_results.write(sig, 0.)

When i remove the # in the lines to calculate the stress the model dont work. I hope you can help me.

2

This is the error

WhatsApp Image 2022-09-08 at 11.54.10 AM
WhatsApp Image 2022-09-08 at 11.54.10 AM1001×457 189 KB

3

You can change the solver in project, as shown in: Out of memory error - #4 by dokken

1 Like
4

Thanks Dokken, its work perfect!

5

Hi dokken, im working whit a bigger 3d mesh and i cant run my model. I am using this solver:

u = Function(V,name="Displacement")
solve(a==l, u, bc, solver_parameters={"linear_solver" : "gmres", "preconditioner" : "ilu"},form_compiler_parameters={"optimize":True})

And the error is:

imagen
imagen1466×515 78.3 KB

Hope you can help me

6

In general I would not use ilu to solve linear elasticity problems.

I would use cg and gamg (and attaching the near nullspace).

Anyhow, without an example reproducing the error, I cant give any more concrete suggestions

7

Im sorry, my code for the error is:

from dolfin import*
from mshr import*
from fenics import*
from ufl import nabla_div
import numpy as np

mesh = Mesh("FASE2.xml")
E = Constant(20e3) #Módulo de Young
nu = Constant(0.3) #Radio de Poisson
mu = E/(2*(1+nu)) #Módulo de cizalle
lambda_ = E*nu/((1+nu)*(1-2*nu)) #Lamé
rho = Constant(0.027) #densidad

# Strain function
def epsilon(u):
    return sym(nabla_grad(u))

# Stress function
def sigma(u):
    return lambda_*nabla_div(u)*Identity(3) + 2*mu*epsilon(u)
# lambda is a reserved python keyword, naming convention recommends using a single trailing underscore for such cases

#Load
g_z = 0.0
b_z = -rho
g = Constant((0.0,g_z,0.0))
b = Constant((0.0,0.0,b_z))

#Function Spaces
V = VectorFunctionSpace(mesh, "P", 2)
u_tr = TrialFunction(V)
u_test = TestFunction(V)

a = inner(sigma(u_tr),epsilon(u_test))*dx #bilineal form
l = dot(b,u_test)*dx + dot(g, u_test)*ds #lineal form

x_max = 1216.0
y_max = 929.0
tol=1e-8
def bottom(x, on_boundary):
    return (on_boundary and near(x[2],0,tol))
bc1 =DirichletBC(V.sub(2), Constant(0.),bottom)

def left(x, on_boundary):
    return (on_boundary and near(x[0], 0,tol))
bc2 =DirichletBC(V.sub(0), Constant(0.),left)

def right(x, on_boundary):
    return (on_boundary and near(x[0], x_max,tol))
bc3 =DirichletBC(V.sub(0), Constant(0.),right)

def north(x, on_boundary):
    return (on_boundary and near(x[1], y_max,tol))
bc4 =DirichletBC(V.sub(1), Constant(0.),north)

def south(x, on_boundary):
    return (on_boundary and near(x[1], 0,tol))
bc5 =DirichletBC(V.sub(1), Constant(0.),south)
bc=[bc1,bc2,bc3,bc4,bc5]

u = Function(V,name="Displacement")
solve(a==l, u, bc, solver_parameters={"linear_solver" : "mumps", "preconditioner" : "ilu"},form_compiler_parameters={"optimize":True})


file_results = XDMFFile("Verde_CA.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
8

The point is that this does not help, as you have not supplied the mesh.xml file, or any way of generating the mesh file.

As I said above, I don’t think you are using a good solver for elasticity problems.

9

How can i share the xml file??

10

i can give you the file in one drive link:

11

As stated now three times, this is most likely due to your solver setup.
See for instance Bitbucket for a sensible configuration of your solver.

1 Like
12

Thanks Dokken it works perfect! And sorry for your lost time and my insistance.