Tetrahedral mesh deformation case process

General
1

I would like to analyze the process of deformation of an object of elastic material when it is subjected to a squeezing force using a tetrahedral mesh that I have defined. I would like to know why the code I wrote didn’t work. I hope someone can solve my problem. Using dolfinx version 0.8.0 libraries

import numpy as np
import dolfinx
import ufl
from mpi4py import MPI
from scipy.spatial import Delaunay, ConvexHull
import dolfinx.mesh
import dolfinx.io
import basix
import open3d as o3d
import dolfinx.fem.petsc

def read_point_cloud(filename):
    points = np.loadtxt(filename, delimiter=" ")
    points = points[:, 0:3]
    return points

def create_tetrahedral_mesh(points, simplices):
    
    tetrahedra = simplices
    element =  basix.ufl.element("CG", 'tetrahedron', 1, shape=(3,))
    mesh_tetrahedral = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, tetrahedra,points, element)
    return mesh_tetrahedral



def tetrahedron_edge_lengths(vertices):
    edges = [
        np.linalg.norm(vertices[i] - vertices[j])
        for i in range(4) for j in range(i + 1, 4)
    ]
    return edges
def dense_point(points):
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)
    sparse_pcd = pcd

    points = np.asarray(sparse_pcd.points)
    hull = ConvexHull(points)

    mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(sparse_pcd, 0.03)
    mesh.compute_vertex_normals()

    voxel_size = 0.006  
    voxel_grid = o3d.geometry.VoxelGrid.create_from_triangle_mesh_within_bounds(
    mesh, voxel_size, mesh.get_min_bound(), mesh.get_max_bound())

    def is_point_in_hull(point, hull):
        return np.all(np.dot(hull.equations[:, :-1], point) + hull.equations[:, -1] <= 1e-12)

    dense_points = []
    for voxel in voxel_grid.get_voxels():
        voxel_center = voxel.grid_index * voxel_grid.voxel_size + voxel_grid.origin
        if is_point_in_hull(voxel_center, hull):
            dense_points.append(voxel_center)

    dense_pcd = o3d.geometry.PointCloud()
    dense_pcd.points = o3d.utility.Vector3dVector(np.array(dense_points))

    merged_pcd = sparse_pcd + dense_pcd

    points = np.asarray(merged_pcd.points)
    return points

def exclud_tetrahedron(points, length_threshold):
    tri = Delaunay(points)

    used_points_indices = set()
    filtered_simplices = []
    for simplex in tri.simplices:
        edge_lengths = tetrahedron_edge_lengths(points[simplex])

        if all(length < length_threshold for length in edge_lengths):
            filtered_simplices.append(simplex)
            used_points_indices.update(simplex)  
    return filtered_simplices, used_points_indices
            
def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u):
    return lmbda * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu * epsilon(u)

def boundary(x):
    return np.isclose(x[2], 0)
    # return np.full(x.shape[0], np.isclose(x[:, 2], 0))


def normalize_point_cloud(points):
    centroid = np.mean(points, axis=0)
    normalized_points = points - centroid
    return normalized_points
#################################

if __name__ == "__main__":

    points = read_point_cloud("/home/jiao/chen_test/code/2f347d95d4a2dae3e762cf5917cef4ef.txt")

    points = dense_point(points)
    
    one = Delaunay(points).simplices

    # print(type(one))

    simplices, used_points_indices = exclud_tetrahedron(points, 0.05)

    filtered_dense_points_all = points[list(used_points_indices)]
 
    simplices_dense, used_points_indices_dense = exclud_tetrahedron(filtered_dense_points_all, 0.05)
    finally_simplices = np.array(simplices_dense)

    tetrahedral_mesh = create_tetrahedral_mesh(filtered_dense_points_all, finally_simplices)
    tetrahedral_mesh_deformed = create_tetrahedral_mesh(filtered_dense_points_all, finally_simplices)

    #################################
    E = 1000  
    nu = 0.3  
    mu = E / (2 * (1 + nu))  
    lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))  
    #################################

    V = dolfinx.fem.functionspace(tetrahedral_mesh, ("CG", 1 ,(tetrahedral_mesh.geometry.dim,)))

    #################################
    u_bc = dolfinx.fem.Function(V)
    u_bc.interpolate(lambda x: np.zeros((3, x.shape[1])))

    bc = dolfinx.fem.dirichletbc(u_bc, dolfinx.fem.locate_dofs_geometrical(V, boundary))
    #################################

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
    L = ufl.inner(dolfinx.fem.Constant(tetrahedral_mesh, np.array([0.0, 0.0, -10000.0], dtype=np.float64)), v) * ufl.dx

    problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
    u = problem.solve()

    u_values = u.x.array
######################################
    with dolfinx.io.VTKFile(MPI.COMM_WORLD, "original_mesh1.vtk", "w") as vtkfile:
        vtkfile.write_mesh(tetrahedral_mesh)

    original_coords = tetrahedral_mesh.geometry.x.copy()
    deformed_coords = original_coords + u.vector.array.reshape(-1, 3)

    tetrahedral_mesh_deformed.geometry.x[:] = deformed_coords

    with dolfinx.io.VTKFile(MPI.COMM_WORLD, "deformed_mesh1.vtk", "w") as vtkfile:
        vtkfile.write_mesh(tetrahedral_mesh_deformed)
######################################
2

Please explain why it does not work. Provide an error message and try to strip down your problem to a MWE for making it easier for people to help you.

3

Thank you for your answer. I would like to use point cloud to generate tetrahedral meshes that simulates an object undergoing deformation. When the above code is tested with a cubic mesh there is no problem and the corresponding deformation occurs. However, when changing to my own data the object will not change. I don’t know why?

4

Are you sure the mesh is correctly generated? You can try interpolating a known analytical function to check. Then be sure to apply correctly boundary conditions.

5

Thank you for your answer. It should be a problem with my mesh generation. I generated tetrahedra from the original point cloud using Delaunay triangulation and then used it directly in FEA to be able to deform. However, I want to refine the object, so after removing the larger tetrahedra and the point cloud that is not divided into tetrahedra, no deformation occurs, and I don’t know what the problem is?