Thanks a lot
About :
I believe I finally managed to build a 1D elements mesh in 3D space without gmesh nor meshio with the following :
# 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+1, i+2] 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)
says :
geometrical dimension : 3
topological dimension : 1
Which is what I need for future use
But then determine_point_ownership in compute_cell_contributions fails to return the cells :
Cells: []
Basis values: []
I tried to get some help from this implementation that seems to be 3D
But the code fails at the same place.
I am still working on this but any hint would be welcome.
Here bellow is my code without (?) the Colab stuff:
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, CellType
import basix
import matplotlib.pyplot as plt
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-6)
# 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.geometry.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
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+1, i+2] 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)
# 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
(gdim,))) # function space dimension
# Define the variational problem
u = ufl.TrialFunction(V) # Displacement (unknown)
v = ufl.TestFunction(V) # Test function
E = 1e7 # Young's modulus
I = 1e-4 # Moment of inertia
A = 1e-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 * 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)
b.x.array[:] = 0
# Apply the point load at the free end of the beam
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, 3), 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)
##############################
# not tested after this point
##############################
for cell, basis_value in zip(cells, basis_values):
dofs = V.dofmap.cell_dofs(cell)
b.x.array[dofs] += basis_value * -1.0 # Applying a negative load
# Apply the boundary conditions to the right-hand side vector
dolfinx.fem.petsc.apply_lifting(b.vector, [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
uh.x.array