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

**URL:** <https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225>\
**Category:** General\
**Created:** [November 6, 2024, 7:28pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225 "2024-11-06T19:28:13Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [November 6, 2024, 7:28pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/1 "2024-11-06T19:28:13Z")

</div>

Hi, I am trying to run a modified example of the linear elasticity example in:  
[https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity\_code.html](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

```auto
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)

```

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [November 6, 2024, 7:45pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/2 "2024-11-06T19:45:57Z")

</div>

> [@alex24](#):
>
> ```auto
> mesh.topology.create_connectivity(mesh.topology.dim, 0)
> 
> ```

You probably need

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

```

or

```auto
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?](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21) and adhere to those guidelines

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [November 6, 2024, 8:31pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/3 "2024-11-06T20:31:29Z")

</div>

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:

```auto
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?)

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [November 6, 2024, 11:46pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/4 "2024-11-06T23:46:01Z")

</div>

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

```auto
# 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)

```

---

<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:** [November 7, 2024, 6:01am UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/5 "2024-11-07T06:01:06Z")

</div>

> [@alex24](#):
>
> ```auto
> mesh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, domain)
> mesh.topology.create_connectivity(mesh.topology.dim, 0)
> 
> ```

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:

> [@alex24](#):
>
> ```auto
> # 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)
> 
> ```

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

> [@alex24](#):
>
> 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)

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

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [November 7, 2024, 2:41pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/6 "2024-11-07T14:41:16Z")

</div>

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’

```auto
# 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)

```

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [November 7, 2024, 2:49pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/7 "2024-11-07T14:49:55Z")

</div>

See the previous reply

> [@dokken](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [November 7, 2024, 5:23pm UTC](https://fenicsproject.discourse.group/t/generating-mesh-and-applying-bc-for-3d-linear-elasticity-problem-using-known-nodes-and-tetrahedron-connectivity/16225/8 "2024-11-07T17:23:44Z")

</div>

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

```auto
# 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)

```