Hi all,
I’m currently using DOLFINx version 0.8.0 and attempting to apply a fixed boundary condition on one edge of a 3D model. However, I am encountering an error.
untimeError: Entity-to-cell connectivity has not been computed.
Here is the example, a simple elasticity problem with a cube under uniform pressure, aiming to apply a fixed boundary condition to one bottom edge
import dolfinx
from mpi4py import MPI
import numpy as np
from dolfinx.fem import (Constant, Function, functionspace, dirichletbc, locate_dofs_topological, Expression)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import VTXWriter
from dolfinx.mesh import locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction, dx, indices, as_tensor, Measure)
from basix.ufl import element
# -------------------------Material Properties
E = 1 * 130e9
nu = 0.28
laa = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# ......................................... mesh file
H= 0.76e-3
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([H, H, H])],[15, 15, 15], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200
Area2 = H**2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
# ......................................... Define element
EL = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)) # Element
Space = functionspace(mesh, EL)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-10
def Edge(x):
return np.logical_and(np.isclose(x[0], 0, atol=tol), np.isclose(x[2], 0, atol=tol))
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]
###############
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# # ...................................... BCs
clamp_facets_edge_x = locate_entities_boundary(mesh, mesh.topology.dim - 2, Edge)
clamp_dofs_x_edge = locate_dofs_topological((Space.sub(0)), mesh.topology.dim - 2, clamp_facets_edge_x)
bc0_edge_x = dirichletbc(Constant(mesh, ScalarType(0.0)), clamp_dofs_x_edge, Space.sub(0))
bc = [bc0_edge_x]
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l], (i, j, k, l)) # Elasticity tensor
def StressT(u):
return as_tensor((Elas(laa, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])
# ......................................... Define Variational Formulation
a = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
L = tr_S2[i] * del_u[i] * ds(2)
################
Result = Function(Space)
problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()
I applied this setup without issues in an older version. I’ve searched for solutions but haven’t identified the problem in this version.
Thanks