Consider the following MWE:
import ufl
import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import dolfinx.fem.petsc
import dolfinx.io
import ufl
from dolfinx import fem
from dolfinx.mesh import create_mesh
import basix
def compute_cell_contributions(V, points):
"""
Determine the cells that contain the specified points and compute the basis
function values at those points.
Args:
V (FunctionSpace): The finite element function space.
points (np.array): Array of points to be evaluated.
Returns:
cells (np.array): Array of cell indices containing the points.
basis_values (np.array): Array of basis function values at the points.
"""
# Extract the mesh from the function space
mesh = V.mesh
# Determine point ownership using the mesh's CPP object and points.
# `1e-6` is a tolerance value to handle numerical precision issues.
ownership_data = dolfinx.cpp.geometry.determine_point_ownership(
mesh._cpp_object, points, 1e-8
)
# Get the points and cells that own the specified points
owning_points = ownership_data.dest_points
cells = ownership_data.dest_cells
# Convert owning points to a numpy array and reshape to (num_points, 3)
owning_points = np.asarray(owning_points).reshape(-1, 3)
# Get mesh node coordinates
mesh_nodes = mesh.geometry.x
# Get the coordinate mapping (cmap) for the mesh
cmap = mesh.geometry.cmap
# Initialize an array to hold reference coordinates
ref_x = np.zeros((len(cells), mesh.topology.dim), dtype=mesh.geometry.x.dtype)
# Loop over each point and corresponding cell
for i, (point, cell) in enumerate(zip(owning_points, cells)):
# Get the degrees of freedom (dofs) for the cell
geom_dofs = mesh.geometry.dofmap[cell]
# Pull back the point to reference coordinates using the coordinate map (cmap)
ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
# Define a trial function in the function space V
u = ufl.TrialFunction(V)
# Get the number of degrees of freedom per cell and multiply by block size
num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
# If there are any cells that own the points
if len(cells) > 0:
# Create an expression for evaluating the basis functions
expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
# Evaluate the expression on the mesh for the given cells
values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))
# Extract the basis values for the degrees of freedom
basis_values = values[: num_dofs : num_dofs * len(cells)]
else:
# If no cells own the points, create an empty basis_values array
basis_values = np.zeros((0, num_dofs), dtype=dolfinx.default_scalar_type)
# Return the cells and corresponding basis function values
return cells, basis_values
# Geometral
length = 1.0
# Define the number of intervals and the end points
n_intervals = 10
# Geometral
length = 1.0
# Define the number of intervals and the end points
n_intervals = 10
end_points = [0.0, length] # along x[0]
# Create points and cells for the interval mesh
points = np.linspace(end_points[0], end_points[1], n_intervals * 2 + 1)
cells = np.array(
[[i, i + 2, i + 1] for i in range(0, 2 * n_intervals, 2)], dtype=np.int64
)
# Initialize a new 3D coordinate array
coordinates = np.zeros((points.shape[0], 3))
coordinates[:, 0] = points # Set the x-coordinates
# Create the UFL element
element = basix.ufl.element(
"Lagrange",
"interval", # 1D cells
2, # degree of the finite element
shape=(3,),
) # 3D space
# Create the mesh
domain = create_mesh(MPI.COMM_WORLD, cells, coordinates, element)
# from dolfinx.mesh import create_interval
# domain = create_interval(MPI.COMM_WORLD, 10, [0, length])
# Get the geometrical and topological dimensions of the new mesh
gdim = domain.geometry.dim
tdim = domain.topology.dim
# from dolfinx.mesh import create_interval
# domain = create_interval(MPI.COMM_WORLD, n_intervals, [0, length])
# Get the geometrical and topological dimensions of the new mesh
gdim = domain.geometry.dim
tdim = domain.topology.dim
print("geometrical dimension : ", gdim)
print("topological dimension : ", tdim)
# Define the function space
V = fem.functionspace(
domain,
(
"Lagrange",
2, # polynomial degree
(3,),
),
) # function space dimension
# Define the variational problem
u = ufl.TrialFunction(V) # Displacement (unknown)
v = ufl.TestFunction(V) # Test function
E = dolfinx.fem.Constant(domain, 1.0e7) # Young's modulus
I = dolfinx.fem.Constant(domain, 1.0e-4) # Moment of inertia
A_c = dolfinx.fem.Constant(domain, 1.0e-2) # Cross-sectional area
# Define the stiffness matrix for a beam in 3D
a = (
E * I * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
+ A_c * ufl.inner(u, v) * ufl.dx
)
a_compiled = dolfinx.fem.form(a)
# Apply boundary condition (fixed at the left end)
dofs = dolfinx.fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
zero = dolfinx.fem.Function(V)
zero.x.array[:] = 0
bc = dolfinx.fem.dirichletbc(zero, dofs)
# Initialize the right-hand side vector
b = dolfinx.fem.Function(V)
if domain.comm.rank == 0: # If this is the root process ...
print("root process")
# define the point load position as a 3D with dtype match the mesh
point_load_position = np.array([[length, 0, 0]], dtype=domain.geometry.x.dtype)
else:
point_load_position = np.zeros(
(0, domain.geometry.dim), dtype=domain.geometry.x.dtype
)
# Compute the contributions to the right-hand side from the point source
# points = point_load_position
cells, basis_values = compute_cell_contributions(V, point_load_position)
print("Cells:", cells)
print("Basis values:", basis_values)
for cell, basis_value in zip(cells, basis_values):
dofs = V.dofmap.cell_dofs(cell)
for i, dof in enumerate(dofs):
for bs in range(V.dofmap.bs):
b.x.array[dof * V.dofmap.bs + i] += basis_value[i * V.dofmap.bs + bs] * -1.0
# Apply the boundary conditions to the right-hand side vector
b_petsc = b.x.petsc_vec
dolfinx.fem.petsc.apply_lifting(b_petsc, [a_compiled], [[bc]])
b.x.scatter_reverse(dolfinx.la.InsertMode.add)
dolfinx.fem.petsc.set_bc(b.vector, [bc])
b.x.scatter_forward()
# Assemble the system matrix
A = dolfinx.fem.petsc.assemble_matrix(a_compiled, bcs=[bc])
A.assemble()
# Set up the linear solver
ksp = PETSc.KSP().create(domain.comm)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.PREONLY)
ksp.getPC().setType(PETSc.PC.Type.LU)
ksp.getPC().setFactorSolverType("mumps")
# Solve the linear system
uh = dolfinx.fem.Function(V)
ksp.solve(b.vector, uh.vector)
uh.x.scatter_forward()
# extract coordinates
DoFs_x = V.tabulate_dof_coordinates()[:, 0]
# extract displacement
print(uh.x.array)
A.destroy()
ksp.destroy()
with dolfinx.io.VTXWriter(domain.comm, "u.bp", [uh], engine="BP4") as bp:
bp.write(0.0)
Note that I changed
to
cells = np.array(
[[i, i + 2, i + 1] for i in range(0, 2 * n_intervals, 2)], dtype=np.int64
)
since the nodes has to following the ordering defined in at DefElement
I also changed several other things in your code.
Note that it is unclear to me why you need the interval to be embedded in 3D. Why cant the interval be a 1D interval with a 3D vector function space defined on it? Does it really need to be a 1D manifold in 3D space?