Applying displacement BCs without geometric function

dolfinx
1

Hello,

I am performing FEA on human long bones. The mesh is imported to dolfinx from a .xdmf file. I have to apply a displacement on the femur’s head however I cannot define a geometric function to find the right dofs or facets with dolfinx.fem.locate_dofs_geometrical or dolfinx.fem.locate_dofs_topological.

To solve this I opened the mesh in Paraview and selected the nodes. I used this in dolfinx to find the closest nodes and their corresponding dof number. However I have to set the BCs for the three components of the displacement and I need only for u_x and u_z. I am trying to build the equivalent to dolfinx.fem.locate_dofs_geometrical((V_u.sub(0), V_u.sub(0).collapse()), top) manually - without geometric function, as I did imposing boundary conditions to the three displacement components (shown in the code below).

I would appreciate any help, thank you :slight_smile:

## Example for illustration

# Mesh construction
filename = "meshes/"+boneID+"_LinMsh_hRef_prox1.xdmf"  

with XDMFFile(MPI.COMM_WORLD, filename, "r") as file:
    mesh = file.read_mesh(name="Grid")

element_u = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree=1, dim=3)
V_u = dolfinx.FunctionSpace(mesh, element_u)
u_ = dolfinx.Function(V_u, name="Boundary Displacement") # Full vector BC
  
def find_ClosestPts(coor_E,coor):
    # Create KD-Tree from the E-coordinates
    tree = KDTree(coor)
    # Initialize an array to store DOF values associated with mesh points
    nearest_values = np.full(coor_E.shape[0], 0)  # Initialize with NaN
    
    for i in range(len(coor_E)):
        # Query the KD-Tree to find the closest point for each mesh node
        dist, index = tree.query(coor_E[i])
        # Check if the distance is within the maximum threshold
        max_distance = 2.0  # Set the maximum allowed distance (optional)
        if dist <= max_distance:
            nearest_values[i] = int(index)
    return nearest_values

def top(coor):
    file_dofs_topBC = "meshes/top_BCs_dofs_" + boneID + "_LinMsh_selRef_prox1.csv"
    # Load the coordinates from the CSV file into a numpy array (x, y, z)
    coor_csv = np.genfromtxt(file_dofs_topBC, delimiter=",", skip_header=1, usecols=(1, 2, 3))
    
    return find_ClosestPts(coor_csv,coor)

coords = mesh.geometry.x  # Get all coordinates of the mesh
dofs_uImp_top = top(coords)
dolfinx.fem.DirichletBC(u_, dofs_uImp_top.astype(np.int32))
2

Could you illustrate this with a built in mesh, (like a unit cube), and specify some nodes, say ((1, 0.345, 0), (0.82, 0.32, 1)) and set up the exact same code? This would make the example reproducible for others and easier for anyone to give some meaningful advice.

So your question is, how to set up a DirichletBC for a sub-space, where you want to find the dofs closest to the aforementioned set of points?

A followup question is then, will you always use first order Lagrange elements, or do you want a solution that works for arbitrary order spaces?

3

Thank your for your quick answer and support :slight_smile:

Indeed you got my question right - I managed to impose a full displacement (ux,uy,uz) with this method but I need only to impose on (ux,uz). Also I would appreciate solutions working for higher order Lagrange elements although something working for linear elements is the most urgent.

# Import required libraries
import numpy as np
import pandas as pd
import dolfinx
import ufl
from mpi4py import MPI
import os

# Find closest point- inspiration from the technique to import E hetero 
from scipy.spatial import KDTree

def find_ClosestPts(coor_E,coor):
    # Create KD-Tree from the E-coordinates
    tree = KDTree(coor)
    # Initialize an array to store DOF values associated with mesh points
    nearest_values = np.full(coor_E.shape[0], 0)  # Initialize with NaN
    
    for i in range(len(coor_E)):
        # Query the KD-Tree to find the closest point in coor_E for each mesh node
        dist, index = tree.query(coor_E[i])
        # Check if the distance is within the maximum threshold
        max_distance = 2.0  # Set the maximum allowed distance (optional)
        if dist <= max_distance:
            # Assign the E value from coor_E to the closest mesh point
            nearest_values[i] = int(index)
    return nearest_values
    
# Parameters
tdim= 3 # Dimensions of the problem


# The top function is specific since we will direclty provide the degrees of freedom 
def top(coor):
    coor_csv = np.array([[1,0.345,0],[0.82,0.32,1]])
    return find_ClosestPts(coor_csv,coor)


from dolfinx import BoxMesh, DirichletBC, Function, VectorFunctionSpace, cpp  
from dolfinx.cpp.mesh import CellType             
# Mesh construction
mesh = BoxMesh(
    MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                     np.array([1.0, 1.0, 1.0])], [12, 12, 12],
    CellType.tetrahedron, dolfinx.cpp.mesh.GhostMode.shared_facet)
   
   
# Function space
deg_pol_sol = 1
    
print(mesh.ufl_cell())


element_u = ufl.VectorElement("Lagrange", mesh.ufl_cell(), degree=deg_pol_sol, dim=tdim)
V_u = dolfinx.FunctionSpace(mesh, element_u)
V_ux = V_u.sub(0).collapse()
V_uz = V_u.sub(2).collapse()

# Define the state
u = dolfinx.Function(V_u, name="Displacement")
u_ = dolfinx.Function(V_u, name="Boundary Displacement") # Full vector BC
## ux and uz components of uImp, however we don't know about uy
# u_x = dolfinx.Function(V_ux,name="Boundary Displacement Ux") 
# u_z = dolfinx.Function(V_uz,name="Boundary Displacement Uz") 

# Boundary definitions
coords = mesh.geometry.x  # Get all coordinates of the mesh
dofs_uImp_top = top(coords)

print(dofs_uImp_top)

u_.interpolate(lambda x: (1*np.ones_like(x[0]), 0*np.ones_like(x[1]), np.ones_like(x[2]))) # Full vector as BC
# u_x.interpolate(lambda x: (np.ones_like(x[0])))
# u_z.interpolate(lambda x: (np.ones_like(x[2])))


bcs_u = [dolfinx.fem.DirichletBC(u_, dofs_uImp_top.astype(np.int32))] # Full vector BC
4

You seem to be using a very old version of DOLFINx (v0.3?). Please note that we strongly advice to upgrade to later versions of DOLFINx, as there has been tons of bug fixes, new features etc since that release.

I’ve upgraded the code to 0.9, and added the necessary tooling to solve the problem

# Import required libraries
import numpy as np
import dolfinx.fem.petsc
from mpi4py import MPI
import basix.ufl

# Find closest point- inspiration from the technique to import E hetero
from scipy.spatial import KDTree


def find_ClosestPts(coor_E, coor):
    # Create KD-Tree from the E-coordinates
    tree = KDTree(coor)
    # Initialize an array to store DOF values associated with mesh points
    nearest_values = np.full(coor_E.shape[0], 0, dtype=np.int32)  # Initialize with NaN

    for i in range(len(coor_E)):
        # Query the KD-Tree to find the closest point in coor_E for each mesh node
        dist, index = tree.query(coor_E[i])
        # Check if the distance is within the maximum threshold
        max_distance = 2.0  # Set the maximum allowed distance (optional)
        if dist <= max_distance:
            # Assign the E value from coor_E to the closest mesh point
            nearest_values[i] = int(index)
    return nearest_values


# Parameters
tdim = 3  # Dimensions of the problem


# The top function is specific since we will direclty provide the degrees of freedom
def top(coor):
    coor_csv = np.array([[1, 0.345, 0], [0.82, 0.32, 1]])
    return find_ClosestPts(coor_csv, coor)


import dolfinx

# Mesh construction
mesh = dolfinx.mesh.create_box(
    MPI.COMM_WORLD,
    [np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
    [12, 12, 12],
    dolfinx.mesh.CellType.tetrahedron,
    ghost_mode=dolfinx.mesh.GhostMode.shared_facet,
)


# Function space
deg_pol_sol = 1


element_u = basix.ufl.element(
    "Lagrange", mesh.basix_cell(), degree=deg_pol_sol, shape=(mesh.geometry.dim,)
)
V_u = dolfinx.fem.functionspace(mesh, element_u)
V_ux, x_to_V = V_u.sub(0).collapse()
V_uz, z_to_V = V_u.sub(2).collapse()

# Define the state
u = dolfinx.fem.Function(V_u, name="Displacement")
u_ = dolfinx.fem.Function(V_u, name="Boundary Displacement")  # Full vector BC
## ux and uz components of uImp, however we don't know about uy
# u_x = dolfinx.Function(V_ux,name="Boundary Displacement Ux")
# u_z = dolfinx.Function(V_uz,name="Boundary Displacement Uz")

# Boundary definitions
coords_x = V_ux.tabulate_dof_coordinates()
coords_z = V_uz.tabulate_dof_coordinates()
dofs_uImp_top_x = top(coords_x)
dofs_uImp_top_z = top(coords_z)

dofs_x = (np.asarray(x_to_V, dtype=np.int32)[dofs_uImp_top_x], dofs_uImp_top_x)
dofs_z = (np.asarray(z_to_V, dtype=np.int32)[dofs_uImp_top_z], dofs_uImp_top_z)


ux_ = dolfinx.fem.Function(V_ux)
ux_.x.array[:] = 1

uz_ = dolfinx.fem.Function(V_uz)
uz_.x.array[:] = 3


bcs_u = [
    dolfinx.fem.dirichletbc(ux_, dofs_x, V_u.sub(0)),
    dolfinx.fem.dirichletbc(uz_, dofs_z, V_u.sub(2)),
]  # Full vector BC

test = dolfinx.fem.Function(V_u)
dolfinx.fem.petsc.set_bc(test.x.petsc_vec, bcs_u)
with dolfinx.io.VTXWriter(mesh.comm, "test.bp", [test]) as bp:
    bp.write(0.0)

yielding

image
image1920×1491 26.1 KB