Generating mesh and applying bc for 3D linear elasticity problem using known nodes and tetrahedron connectivity

General
1

Hi, I am trying to run a modified example of the linear elasticity example in:
https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html

Using an imported 3D tetrahedral mesh. I am getting the following error:
RuntimeError: Entity-to-cell connectivity has not been computed.
even when using:
mesh.topology.create_connectivity(mesh.topology.dim, 0)

Why is this happening?

How can I extract the facets to nodes connectivity matrix from the computed mesh?

Thank you,
Alex

points = mesh2.points
cells  = mesh2.cells[1].data

finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))
domain = ufl.Mesh(finiteElement)
mesh   = mesh.create_mesh(MPI.COMM_WORLD, cells, points, domain)
mesh.topology.create_connectivity(mesh.topology.dim, 0)

V = fem.functionspace(mesh, finiteElement)

fdim = mesh.topology.dim - 1
boundary_facets = np.array([1])

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
2

You probably need

mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)

or

mesh.topology.create_connectivity(mesh.topology.dim - 1, 0)

(or the other way around with the arguments in case I remember wrongly).

Having a fully reproducible code would have jogged my memory, but the code you posted is not runnable by anyone. For the future, please read Read before posting: How do I get my question answered? and adhere to those guidelines

3

Hi Francesco, thank you for the prompt reply:

I am trying to write a fully reproducible code by replacing points and cells with the following:

external_mesh = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)

points = external_mesh.geometry.x

Do you know how to extract the tetrahedron connectivity matrix?

How could you extract the surface of the mesh (i.e nodes and triangular face connectivity?)

4

I was able to make the code reproducible. However I get the following error at the end:
AttributeError: ‘Mesh’ object has no attribute ‘topology’

# Common libraries
import numpy as np


# Fenicsx FEA libraries
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import basix


mu = 1
rho = 1
gamma = 10
beta = 1.25
lambda_ = beta
g = gamma


# Make nodes and tetrahedral connectivity matrix
external_mesh = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)

points = external_mesh.geometry.x
cells  = external_mesh.topology.connectivity(3,0).array.reshape((-1,4))


# Import nodes and connectivity matrix into domain and mesh objects
finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))
domain = ufl.Mesh(finiteElement)
mesh   = mesh.create_mesh(MPI.COMM_WORLD, cells, points, domain)
mesh.topology.create_connectivity(mesh.topology.dim, 0)


V = fem.functionspace(mesh, finiteElement)


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


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds


problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
5

You are overwriting the mesh module here. Please be careful when declaring variable names to ensure that you are not reusing the name of an import.

Also I do not understand why you need:

as you already have a dolfinx mesh when defining the external_mesh variable.

Ref this error. If you had posted the full error message, we would be able to see what line of your code triggered this.

Changing that code to tetrahedral elements should be as easy as
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])], [20, 6, 6], cell_type=mesh.CellType.tetrahedron)
and I surprised me if this would trigger your aforementioned error

6

Hi Jorgen,

Thank you for the prompt reply; the reason I am making external_mesh is to have a minimal running code. In the actual code the nodes and connectivity come from a mesh from another library.

I changed the code to not use mesh as a variable name (below). Now the error I get is:
Traceback (most recent call last):

File ~/.local/lib/python3.10/site-packages/spyder_kernels/customize/utils.py:209 in exec_encapsulate_locals
exec_fun(compile(code_ast, filename, “exec”), globals)

File ~/robotDepowdering/simple_fenics_example.py:41
part_mesh.topology.create_connectivity(mesh.topology.dim, 0)

AttributeError: module ‘dolfinx.mesh’ has no attribute ‘topology’

# Common libraries
import numpy as np


# Fenicsx FEA libraries
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import basix


mu = 1
rho = 1
gamma = 10
beta = 1.25
lambda_ = beta
g = gamma


# Make nodes and tetrahedral connectivity matrix
external_mesh = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)

points = external_mesh.geometry.x
cells  = external_mesh.topology.connectivity(3,0).array.reshape((-1,4))


# Import nodes and connectivity matrix into domain and mesh objects
finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))
domain = ufl.Mesh(finiteElement)
part_mesh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, domain)
part_mesh.topology.create_connectivity(mesh.topology.dim, 0)


V = fem.functionspace(part_mesh, finiteElement)


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


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds


problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
7

See the previous reply

8

Hi Francesco, I got it. Here is the solution:

# Common libraries
import numpy as np


# Fenicsx FEA libraries
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import basix


mu = 1
rho = 1
gamma = 10
beta = 1.25
lambda_ = beta
g = gamma


# Make nodes and tetrahedral connectivity matrix
external_mesh = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)

points = external_mesh.geometry.x
cells  = external_mesh.topology.connectivity(3,0).array.reshape((-1,4))


# Import nodes and connectivity matrix into domain and mesh objects
finiteElement = basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3, ))

domain = mesh.create_mesh(MPI.COMM_WORLD, cells, points, finiteElement)

V = fem.functionspace(domain, finiteElement)


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


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

print(domain)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)


T = fem.Constant(domain, default_scalar_type((0, 0, 0)))


ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds



problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)