Can you help me apply compute_cell_contributions method to 1D cantilever beam example?

dolfinx
1

Following bending-of-the-simply-supported-beam I renewed efforts to adapt compute_cell_contributions to my FEniCSx 1D cantilever beam under point load.
Thank you for this already.
Now there is this error I don’t manage to get through :

cells, basis_values = compute_cell_contributions(V, point_load_position)

says :

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-26-03019a981a3c> in <cell line: 13>()
     11 
     12 # Compute the contributions to the right-hand side from the point source
---> 13 cells, basis_values = compute_cell_contributions(V, point_load_position)
     14 for cell, basis_value in zip(cells, basis_values):
     15     dofs = V.dofmap.cell_dofs(cell)

<ipython-input-6-eb0dc5d80503> in compute_cell_contributions(V, points)
     13     """
     14     mesh = V.mesh
---> 15     _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
     16         mesh._cpp_object, points, 1e-6)
     17     owning_points = np.asarray(owning_points).reshape(-1, 3)

TypeError: cannot unpack non-iterable PointOwnershipData_float64 object

Could you please help me through ?

The full code goes like :

# From :

# Create a point source for Poisson problem
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

import os
arch = os.getenv("ARGS", "real")

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl
    import dolfinx
else:
    try:
        import ufl
        import dolfinx
    except ImportError:
        if arch != "complex":
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        else:
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-complex.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        import ufl
        import dolfinx

# Install meshio for mesh handling
!pip install meshio

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import dolfinx.fem.petsc
import dolfinx.io
import ufl
import meshio
from dolfinx import fem

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.
    """
    mesh = V.mesh
    _, _, owning_points, cells = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6)
    owning_points = np.asarray(owning_points).reshape(-1, 3)
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmaps[0]
    ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    u = ufl.TrialFunction(V)
    num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
    if len(cells) > 0:
        expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))
        basis_values = values[:num_dofs:num_dofs*len(cells)]
    else:
        basis_values = np.zeros((0, num_dofs), dtype=dolfinx.default_scalar_type)
    return cells, basis_values

# Create a 1D mesh for the cantilever beam in 3D space using meshio
length = 10.0
N = 10  # Number of elements
points = np.zeros((N+1, 3))  # N+1 array of zeros in 3D space
points[:, 0] = np.linspace(0, length, N+1)  # Distribute points along x from 0 to 'length'
cells = np.array([[i, i+1] for i in range(N)])  # Line elements connectivity
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))  # Write XDMF mesh file

# Reading the mesh from the XDMF file
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid")

# Get the geometrical and topological dimensions of the mesh from the 'domain' object.
gdim = domain.geometry.dim    # = 3
tdim = domain.topology.dim    # = 1

# Define the function space
# V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
V = fem.functionspace(domain,
                      ("Lagrange",
                       1,                            # 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:
    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
cells, basis_values = compute_cell_contributions(V, point_load_position)

#####################
# 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()

# Write the solution to file in VTK format for visualization
with dolfinx.io.VTXWriter(domain.comm, "uh.pvd", [uh], engine="BP4") as bp:
    bp.write(0.0)
2

Edited your post: it’s not polite to tag people in the very first message of a topic. It would have been fine to do if you were replying to their message, but you decided to create a new topic, hence you aren’t replying to anyone.

3

I am really sorry it felt that way.
I was replied in this thread : bending-of-the-simply-supported-beam with the advise of having a look to compute_cell_contributions.
But since I didn’t want to deviate Kumar’s thread from his initial issue, I created this new one and and tagged from here to make the link.
I thought it was the right thing to do, sorry again.

4

Consider the following MWE:

# From :

# Create a point source for Poisson problem
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

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

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.
    """
    mesh = V.mesh
    ownership_data = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6)
    owning_points = ownership_data.dest_points
    cells = ownership_data.dest_cells
    owning_points = np.asarray(owning_points).reshape(-1, 3)
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    u = ufl.TrialFunction(V)
    num_dofs = V.dofmap.dof_layout.num_dofs * V.dofmap.bs
    if len(cells) > 0:
        expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF)
        values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))
        basis_values = values[:num_dofs:num_dofs*len(cells)]
    else:
        basis_values = np.zeros((0, num_dofs), dtype=dolfinx.default_scalar_type)
    return cells, basis_values

domain = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)

# Get the geometrical and topological dimensions of the mesh from the 'domain' object.
gdim = domain.geometry.dim    # = 3
tdim = domain.topology.dim    # = 1

# Define the function space
# V = dolfinx.fem.FunctionSpace(domain, ("Lagrange", 1))
V = fem.functionspace(domain,
                      ("Lagrange",
                       1,                            # 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
length=1
# Apply the point load at the free end of the beam
if domain.comm.rank == 0:
    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
cells, basis_values = compute_cell_contributions(V, point_load_position)

#####################
# not tested after this point

print(cells, basis_values)
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()

# Write the solution to file in VTK format for visualization
with dolfinx.io.VTXWriter(domain.comm, "uh.pvd", [uh], engine="BP4") as bp:
    bp.write(0.0)

Please note that for further questions on the subject, please:

  1. Tell us how you installed DOLFINx (you have several parts of your code indicating that you are using colab and you cannot execute the code you have posted outside of colab/jupyter without modifying it.
  2. Use a built in mesh, to avoid that people have to install additional dependencies
5

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
6

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?

7

Thanks a lot for your time.

My final goal is to calculate a sailboat rig. Mast, boom, spreaders, running and standing lines they can all be represented by 1D elements but they interact with each other in a 3D space. That is why in my mind … is this the right way to do ?

Your corrected code runs, but I can’t get it to match beam theory :

maximal deflection =  -0.000999
beam theory =  -0.000333

and deformation is strange :
Untitled
I would be so glad if you could have a look.
This is how I get displacements :

DoFs_x = V.tabulate_dof_coordinates()[:, 0]
sort = np.argsort(DoFs_x)
u_y = uh.x.array[1::3]
plt.plot(DoFs_x[sort], u_y[sort], 'o-')

My full code is here bellow :

# ===================================
# Import
# ===================================

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

# ================================
# Constants
# ================================

length = 1.0       # beam length
n_intervals = 10   # beam number of intervals
P = -1.            # point load magnitude
E_val = 1.0e7      # Young's modulus
I_val = 1.0e-4     # Moment of inertia

# ================================
# Mesh
# ================================

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
    )

coordinates = np.zeros((points.shape[0], 3))  # init 3D coordinate array

coordinates[:, 0] = points                    # set x-coordinates

element = basix.ufl.element(                  # define element
    "Lagrange",                               # type
    "interval",                               # 1D cells
    2,                                        # degree
    shape=(3,),                               # 3D space
    )

domain = create_mesh(                         # Create the mesh
    MPI.COMM_WORLD,
    cells,
    coordinates,
    element
    )

gdim = domain.geometry.dim                    # Get geometrical dimension
tdim = domain.topology.dim                    # Get topological dimension

# Define the function space
V = fem.functionspace(
    domain,             # mesh
     (
         "Lagrange",    # type
         2,             # degree
         (3,),          # 3D space
         ),
    )

# Define the variational problem
u = ufl.TrialFunction(V)                     # Displacement
v = ufl.TestFunction(V)                      # Test function
E = dolfinx.fem.Constant(domain, E_val)      # Young's modulus
I = dolfinx.fem.Constant(domain, I_val)     # 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
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] * P

# 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 x coordinates
DoFs_x = V.tabulate_dof_coordinates()[:, 0]
sort = np.argsort(DoFs_x)
# extract y displacements
u_y = uh.x.array[1::3]
# Plot the y-deformation along the x-axis
import matplotlib.pyplot as plt
plt.plot(DoFs_x[sort], u_y[sort], 'o-')
plt.xlabel("x")
plt.ylabel("y-deformation")
plt.show()

# maximal deflection at the free end of the beam
max_deflection = u_y[sort][-1]

# Calculate value with beam theory formula
beam_theory_deflection = (P * length**3) / (3 * E_val *I_val)

# Print result
print("Point load only")
print("maximal deflection = ", max_deflection)           #
print("beam theory = ",        beam_theory_deflection)   #

A.destroy()
ksp.destroy()

with dolfinx.io.VTXWriter(domain.comm, "u.bp", [uh], engine="BP4") as bp:
    bp.write(0.0)
8

Is the point load supposed to be in only one of the three coordinate directions (i.e. the y direction?)

And when you say they interact in 3D, does that mean that the lines are going to be curved in 3D space, or will they always be in the plane, but having unknowns corresponding to x,y,z displacement?

9

Here is a slightly modified version that puts a point source only in y-direction:

# ===================================
# Import
# ===================================

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 evalating 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

# ================================
# Constants
# ================================

length = 1.0       # beam length
n_intervals = 15   # beam number of intervals
P = -1.            # point load magnitude
E_val = 1.0e7      # Young's modulus
I_val = 1.0e-4     # Moment of inertia

# ================================
# Mesh
# ================================

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
    )

coordinates = np.zeros((points.shape[0], 3))  # init 3D coordinate array

coordinates[:, 0] = points                    # set x-coordinates

element = basix.ufl.element(                  # define element
    "Lagrange",                               # type
    "interval",                               # 1D cells
    2,                                        # degree
    shape=(3,),                               # 3D space
    )

domain = create_mesh(                         # Create the mesh
    MPI.COMM_WORLD,
    cells,
    coordinates,
    element
    )

gdim = domain.geometry.dim                    # Get geometrical dimension
tdim = domain.topology.dim                    # Get topological dimension

# Define the function space
V = fem.functionspace(
    domain,             # mesh
     (
         "Lagrange",    # type
         2,             # degree
         (3,),          # 3D space
         ),
    )

# Define the variational problem
u = ufl.TrialFunction(V)                     # Displacement
v = ufl.TestFunction(V)                      # Test function
E = dolfinx.fem.Constant(domain, E_val)      # Young's modulus
I = dolfinx.fem.Constant(domain, I_val)     # 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 ...
    # 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
cells, basis_values = compute_cell_contributions(V.sub(1), point_load_position)

for cell, basis_value in zip(cells, basis_values):
    dofs = V.sub(1).dofmap.cell_dofs(cell)
    for i, dof in enumerate(dofs):
        for bs in range(V.sub(1).dofmap.bs):
            b.x.array[dof * V.sub(1).dofmap.bs + bs] += basis_value[i * V.sub(1).dofmap.bs + bs] * P

# 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 y displacements
uy = uh.sub(1).collapse()
coords = uy.function_space.tabulate_dof_coordinates()[:, 0]
x_sort = np.argsort(coords)
# Plot the y-deformation along the x-axis
import matplotlib.pyplot as plt
plt.plot(coords[x_sort], uy.x.array[x_sort], 'o-')
plt.xlabel("x")
plt.ylabel("y-deformation")
plt.savefig("deflection.png")

# maximal deflection at the free end of the beam
max_deflection = uy.x.array[x_sort][-1]

# Calculate value with beam theory formula
beam_theory_deflection = (P * length**3) / (3 * E_val *I_val)

# Print result
print("Point load only")
print("maximal deflection = ", max_deflection)           #
print("beam theory = ",        beam_theory_deflection)   #

A.destroy()
ksp.destroy()

with dolfinx.io.VTXWriter(domain.comm, "u.bp", [uh], engine="BP4") as bp:
    bp.write(0.0)

This gives the plot:

image
image1282×568 19.5 KB

As I do not have a background from solid mechanics, I can’t say if your beam theory estimates are correct. I can only comment that if you change P, the displacement does consistently scale with a factor of 3 compared to your exact solution.

10

Hi and thanks for answering

In this testing example yes
But for futur use I’ll need point load to be any combination of Fx Fy Fz

lines are going to be curved in 3D space

here is a link to the basic “Cantilever beam with end load” I am trying to calculate in FEniCS :

it gives w_c = PL^3/(3EI) for the maximum deflection
But furthermore the deflection is not linear as the above codes suggest

700×723 12.6 KB

11

As I showed above you can select what direction the point source should be put in.

I don’t know how your variational formulation:

links to elasticity theory, as this scaled system has a scaling of 1e3 for the stiffness component, 1e-2 for the mass matrix, which makes this problem very close to a Poisson equation, whose analytical solution in a linear solution.

You would have to provide further context to justify this variational formulation.

12

Hi, this is not the correct variational formulation for beams. Please check
https://bleyerj.github.io/comet-fenicsx/tours/beams/beams_3D/beams_3D.html

1 Like