Distance computation from a surface mesh

Imports

import numpy as np
import pyvista as pv
pv.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots
from Muscat.Simple import *
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.Bridges.PyVistaBridge import MeshToPyVista
from Muscat.Containers.MeshCreationTools import CreateDisk
from Muscat.Containers.MeshFieldOperations import GetFieldTransferOp
from Muscat.Containers.MeshTools import GetElementsCenters

Mesh Creation and Target Points

# create a mesh and be add and extra zeros to obtain a mesh in 3D
diskMesh = CreateDisk(nr = 10, nTheta= 10, r0=0.5, r1= 1., theta0= 0, theta1=np.pi/2, ofTriangles=True)
diskMesh.nodes = np.hstack((diskMesh.nodes,np.zeros((diskMesh.GetNumberOfNodes(),1), dtype=diskMesh.nodes.dtype ) ) )

# define some target points
TargetPoints = np.array(
                        [[1.0,1.0,0],
                         [0.2,0.6,0.6],
                         [-.2,0.8,0],
                         [0.6,0.2,0.0],
                         ])

Transfer Operator

Compute the transfer operator from an isoparametric field mesh to the target points.

• By default the FEField is isoparametric (number of nodes == number of dofs)

• By default the method is “Interp/Clamp”

P, status, entities =  GetFieldTransferOp(FEField('input',diskMesh),TargetPoints)

Extract Used Elements and Compute Elements Centers

# compute the position of the information extraction
dataExtractionPoints = P.dot(diskMesh.nodes)
# generate a mask to plot elements used for the ext
mask =np.zeros(diskMesh.GetNumberOfElements(), dtype=bool)
# be sure to put true only in the case of 1: "Interp"=1 or "Clamp"=3
mask[entities] = np.logical_or(status ==1,status == 3)
# extract elements using the mask
usedElements = ExtractElementsByElementFilter(diskMesh, ElementFilter(eMask=mask))
#compute the center of the used elements
elementCenters = GetElementsCenters(usedElements)

Console Output

np.set_printoptions(precision=3,suppress=True)
print(f"Number of target points: {TargetPoints.shape[0]}")
print("Position if the location to extract the information ")
print(dataExtractionPoints)
print(f"Status : {status.flatten()}")
print(f"Entities : {entities.flatten()}")

Number of target points: 4
Position if the location to extract the information
[[0.704 0.704 0.   ]
 [0.2   0.6   0.   ]
 [0.    0.8   0.   ]
 [0.6   0.2   0.   ]]
Status : [3 3 3 1]
Entities : [188  86 143  75]

Graphical Output

p = pv.Plotter()
p.add_mesh(MeshToPyVista(diskMesh),style="surface",show_edges=True)
p.add_mesh(MeshToPyVista(usedElements), line_width=5,color="RED")
for tp in range(TargetPoints.shape[0]):
    p.add_ruler(dataExtractionPoints[tp,:], TargetPoints[tp,:], number_labels=2, show_labels=False, title="")
    #p.add_mesh(pv.Sphere(radius=0.001, center= TargetPoints[tp,:]))
    p.add_points(TargetPoints[tp,:],style="points")
    p.add_points(TargetPoints, render_points_as_spheres=True, point_size=1)
    p.add_point_labels(TargetPoints, labels= list(map(lambda x: "Target Point:"+str(x),range(TargetPoints.shape[0]))),shape_opacity=0.5, bold=False)
    p.add_point_labels(elementCenters, labels= list(map(lambda x : "Source Element:" +str(x) ,entities.flatten())),point_size=5, shape_opacity=0.5,justification_vertical="top",justification_horizontal="right",fill_shape=False,shape=None, bold=False)
p.show()