Hyperelasticity example issues

General
1 
#add code here

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl
import basix

from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 10.0
                                #############For 3D case#######################
# domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)
# B = fem.Constant(domain, default_scalar_type((0, 0,0)))
# T = fem.Constant(domain, default_scalar_type((0, 0,0)))

                                #############For 2D case#######################
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0.0, 0.0), (10.0,1.0)), n=(20,5), cell_type=mesh.CellType.quadrilateral)
B = fem.Constant(domain, default_scalar_type((0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0)))

#element = ufl.vectorelement("Lagrange", domain.ufl_cell(),1)
#V=fem.FunctionSpace(domain, element)
#V = fem.VectorFunctionSpace(mesh, "Lagrange", 1)
element = basix.ufl.element("Lagrange", basix.CellType.quadrilateral, 1)
mesh = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0.0, 0.0), (10.0,1.0)), n=(20,5), cell_type=mesh.CellType.quadrilateral) #dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, 1, 1, 1, cell_type=dolfinx.mesh.CellType.hexahedron)
V = dolfinx.fem.functionspace(domain, element)


def left(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
left_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, left)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets])
marked_values = np.hstack([np.full_like(left_facets, 1)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))

bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

v = ufl.TestFunction(V)
u = fem.Function(V)

# Identity tensor
I = ufl.variable(ufl.Identity(3))
# Deformation gradient
def grad_3D(u):
    return ufl.as_matrix([[u[0].dx(0), u[0].dx(1), 0], 
                          [u[1].dx(0), u[1].dx(1), 0], 
                          [0, 0, 0]])
F = ufl.variable(I + grad_3D(u))

# Spatial dimension
d = len(u)
# Identity tensor
#I = ufl.variable(ufl.Identity(d))
# Deformation gradient
#F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = ((mu/2) * (((J**(-2./3.)) * Ic) - 3)) + ((lmbda/2) * (((J**2.) /2) - (1/2) - ufl.ln(J)) )
ee=-1.5*mu
dd=mu
#psi = (mu/2*fem.inner(F, F) + ee - dd*fem.ln(J))*fem.dx - rho_0*J*fem.dot(b, u)*fem.dx + surface_def*fe.inner(g, u)*ds(1) #psi = (aa*fe.inner(F, F) + ee - dd*fe.ln(J))*fe.dx - rho_0*J*fe.dot(b, u)*fe.dx + surface_def*fe.inner(g, u)*ds(1)
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(grad_3D(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds

problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.solve(u)

I’m trying to solve the hyperelasticity example for a 2d case, this is my code, but I’m having several issues, starting from the dirichlet BCs where I receive this error message:
RuntimeError: Rank mis-match between Constant and function space in DirichletBC
could you please check my code.
Thank you

2

Start by changing this to

u_bc = np.full(domain.geometry.dim, 0., dtype =default_scalar_type)

Next time, please post the full stack trace of your error message.

3

Thank you very much dokken I’ll try it

4

Hi dokken, it seems that nothing changes. This is the error:

TypeError Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/dolfinx/fem/bcs.py in dirichletbc(value, dofs, V)
186 try:
→ 187 bc = bctype(_value, dofs, V)
188 except TypeError:

TypeError: init(): incompatible function arguments. The following argument types are supported:
1. init(self, g: ndarray[dtype=float64, writable=False, order=‘C’], dofs: ndarray[dtype=int32, writable=False, shape=(), order=‘C’], V: dolfinx.cpp.fem.FunctionSpace_float64) → None
2. init(self, g: dolfinx.cpp.fem.Constant_float64, dofs: ndarray[dtype=int32, writable=False, shape=(), order=‘C’], V: dolfinx.cpp.fem.FunctionSpace_float64) → None
3. init(self, g: dolfinx.cpp.fem.Function_float64, dofs: ndarray[dtype=int32, writable=False, shape=(), order=‘C’]) → None
4. init(self, g: dolfinx.cpp.fem.Function_float64, dofs: collections.abc.Sequence[ndarray[dtype=int32, writable=False, shape=(), order=‘C’]], V: dolfinx.cpp.fem.FunctionSpace_float64) → None

Invoked with types: dolfinx.cpp.fem.DirichletBC_float64, ndarray, ndarray, dolfinx.fem.function.FunctionSpace

During handling of the above exception, another exception occurred:

RuntimeError Traceback (most recent call last)
1 frames
/usr/local/lib/python3.10/dist-packages/dolfinx/fem/bcs.py in dirichletbc(value, dofs, V)
187 bc = bctype(_value, dofs, V)
188 except TypeError:
→ 189 bc = bctype(_value, dofs, V._cpp_object)
190 else:
191 bc = bctype(_value, dofs)

RuntimeError: Rank mis-match between Constant and function space in DirichletBC

5

If you want an element for a vector space you need to add

element = basix.ufl.element("Lagrange", basix.CellType.quadrilateral, 1, shape=(mesh.geometry.dim, )

7

Thank you dokken, now it works, however I’m getting a displacement plot that is strange.
This is the code I’m using:

import dolfinx.plot as plot
import pyvista
#pyvista.start_xvfb()
#print(pyvista.global_theme.jupyter_backend)
pyvista.global_theme.jupyter_backend = "html"
print(pyvista.global_theme.jupyter_backend)
#pyvista.set_jupyter_backend("html")
pyvista.global_theme.notebook=True

p = pyvista.Plotter()
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

u_3D = np.zeros((u_geometry.shape[0], 3))
u_3D[:, :2] = u.x.array.reshape(-1, 2)
grid.point_data["Displacement"] = u_3D
grid.set_active_vectors("Displacement")
warped = grid.warp_by_vector("Displacement", factor=1000)

#p = pyvista.Plotter()
p.window_size = (800, 300)
p.add_mesh(warped)
#edges = warped.extract_all_edges()
#p.add_mesh(edges, color="k", line_width=1, opacity=0.5)
p.view_xy()
p.zoom_camera(2.5)

p.show()

DISPLbeam
DISPLbeam978×362 11.5 KB

(The undeformed is horizontal).
Do you know how to fix it? It’s ok also to have only the contour of the beam in the undeformed shape.
Thank you

8

This is unrelated to the initial question.
Please also note that the code you give in this reply cannot reproduce the issue.

9

Thank you dokken

#add code here

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl
import basix

from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 10.0
                                #############For 3D case#######################
# domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)
# B = fem.Constant(domain, default_scalar_type((0, 0,0)))
# T = fem.Constant(domain, default_scalar_type((0, 0,0)))

                                #############For 2D case#######################
domain = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0.0, 0.0), (10.0,1.0)), n=(20,5), cell_type=mesh.CellType.quadrilateral)
B = fem.Constant(domain, default_scalar_type((0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0)))

#element = ufl.vectorelement("Lagrange", domain.ufl_cell(),1)
#V=fem.FunctionSpace(domain, element)
#V = fem.VectorFunctionSpace(mesh, "Lagrange", 1)
element = basix.ufl.element("Lagrange", basix.CellType.quadrilateral, 1, shape=(mesh.geometry.dim, )
mesh = mesh.create_rectangle(comm=MPI.COMM_WORLD, points=((0.0, 0.0), (10.0,1.0)), n=(20,5), cell_type=mesh.CellType.quadrilateral) #dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, 1, 1, 1, cell_type=dolfinx.mesh.CellType.hexahedron)
V = dolfinx.fem.functionspace(domain, element)


def left(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
left_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, left)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets])
marked_values = np.hstack([np.full_like(left_facets, 1)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.full(domain.geometry.dim, 0., dtype =default_scalar_type)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))

bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

v = ufl.TestFunction(V)
u = fem.Function(V)

# Identity tensor
I = ufl.variable(ufl.Identity(3))
# Deformation gradient
def grad_3D(u):
    return ufl.as_matrix([[u[0].dx(0), u[0].dx(1), 0], 
                          [u[1].dx(0), u[1].dx(1), 0], 
                          [0, 0, 0]])
F = ufl.variable(I + grad_3D(u))

# Spatial dimension
d = len(u)
# Identity tensor
#I = ufl.variable(ufl.Identity(d))
# Deformation gradient
#F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = ((mu/2) * (((J**(-2./3.)) * Ic) - 3)) + ((lmbda/2) * (((J**2.) /2) - (1/2) - ufl.ln(J)) )
ee=-1.5*mu
dd=mu
#psi = (mu/2*fem.inner(F, F) + ee - dd*fem.ln(J))*fem.dx - rho_0*J*fem.dot(b, u)*fem.dx + surface_def*fe.inner(g, u)*ds(1) #psi = (aa*fe.inner(F, F) + ee - dd*fe.ln(J))*fe.dx - rho_0*J*fe.dot(b, u)*fe.dx + surface_def*fe.inner(g, u)*ds(1)
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(grad_3D(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds

problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.solve(u)

That’s the code that works and the lines for the plot are the ones I have written before.
If you prefer, I can open a new discussion.

10

Please open a new post, and remove all code not required to reproduce the result. In your code you have several lines of code included that are commented out.

11

Ok, then I’ll open a new one