# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
import copy
from typing import Optional, Tuple
import numpy as np
from scipy.spatial import KDTree, distance
import time
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.ElementsContainers import AllElements, ElementsContainer
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.Tags import Tags
import Muscat.Helpers.Logger as Logger
[docs]def FindDuplicates(nodes:ArrayLike, tol:Optional[MuscatFloat]=1e-16) -> np.ndarray:
"""Find duplicate pairs of points. Two points are coincident if the distance between is lower than tol
Parameters
----------
nodes : ArrayLike
array of with the position of the nodes (number of nodes x number of coordinates)
tol : Optional[MuscatFloat], optional
tolerance with respect to the bounding box to detect duplicated nodes, by default 1e-16
Returns
-------
np.ndarray
nx2 array of coincident pairs (sorted such that res[:,1]>res[:,0] and res[:,0] is sorted)
"""
index = KDTree(nodes)
inds = index.query_pairs(r=tol)
if len(inds) == 0:
return np.empty((0,2), dtype=MuscatIndex)
duplicate_inds = np.array(list(inds), dtype=MuscatIndex)
duplicate_inds = np.sort(duplicate_inds, axis=1) # such that duplicate_inds[:,1]>duplicate_inds[:,0]
idx = np.argsort(duplicate_inds[:,0]) # such that duplicate_inds[:,0] are sorted (masters first)
duplicate_inds = duplicate_inds[idx,:]
return duplicate_inds
[docs]def CleanDoubleNodes(mesh: Mesh, tol: Optional[MuscatFloat]=None, nodesToTestMask: ArrayLike=None):
"""Remove double nodes for the input mesh
Parameters
----------
mesh : Mesh
the in put mesh
tol : Optional[MuscatFloat], optional
the tolerance, by default value is = np.linalg.norm(boundingMax - boundingMin)*1e-7
if tol is zero a faster algorithm is used
nodesToTestMask : ArrayLike, optional
a mask of len number of nodes with the value True for nodes to be tested (this can increase the speed of the algorithm),
tests all pairs of nodes present in the mask
by default None
"""
boundingMin,boundingMax= mesh.ComputeBoundingBox()
if tol is None:
tol = np.linalg.norm(boundingMax -boundingMin)*1e-7
nbNodes = mesh.GetNumberOfNodes()
toKeep = np.zeros(nbNodes, dtype=bool)
newIndex = np.zeros(nbNodes, dtype=MuscatIndex)
#---# find duplicates
if nodesToTestMask is None:
duplicate_inds = FindDuplicates(mesh.nodes, tol=tol)
else:
duplicate_inds = FindDuplicates(mesh.nodes[nodesToTestMask], tol=tol)
duplicate_inds = np.arange(mesh.GetNumberOfNodes(), dtype=MuscatIndex)[nodesToTestMask][duplicate_inds]
if len(duplicate_inds) == 0:
return
master_slave = np.arange(nbNodes)
for master, slave in duplicate_inds:
if master_slave[master] != master:
master_slave[slave] = master_slave[master]
else:
master_slave[slave] = master
#---# fill newIndex
maskNodes = (master_slave == np.arange(nbNodes))
newIndex[maskNodes] = np.arange(np.sum(maskNodes))
newIndex[~maskNodes] = newIndex[master_slave][~maskNodes]
toKeep[maskNodes] = True
mesh.nodes = mesh.nodes[toKeep,:]
mesh.originalIDNodes = np.ascontiguousarray(np.where(toKeep)[0],dtype=MuscatIndex)
for tag in mesh.nodesTags :
tag.SetIds(np.unique(newIndex[tag.GetIds()]) )
for elementName in mesh.elements.keys():
elements = mesh.elements[elementName]
elements.connectivity = newIndex[elements.connectivity]
mesh.PrepareForOutput()
[docs]def CleanLonelyNodes(mesh:Mesh, inPlace:bool=True)-> Tuple[np.ndarray, Mesh]:
"""Remove nodes not used by the elements
Parameters
----------
mesh : Mesh
the input mesh
inPlace : bool, optional
if true the nodes are removed in place, by default True
Returns
-------
(usedNodes, cleanMesh)
usedNodes: a mask of type ndarray (size the number of nodes in mesh) with the nodes present on the cleanMesh
cleanMesh: Mesh the mesh without lonely nodes
"""
usedNodes = np.zeros(mesh.GetNumberOfNodes(),dtype=bool )
for data in mesh.elements.values():
usedNodes[data.connectivity.ravel()] = True
cpt = 0
NewIndex = np.zeros(mesh.GetNumberOfNodes(),dtype=MuscatIndex )-1
for n in range(mesh.GetNumberOfNodes()):
if usedNodes[n]:
NewIndex[n] = cpt
cpt += 1
originalIDNodes = np.where(usedNodes)[0]
#filter the nodes
if inPlace:
mesh.nodes = mesh.nodes[usedNodes ,:]
mesh.originalIDNodes = mesh.originalIDNodes[originalIDNodes]
newTags = Tags()
#node tags
for tag in mesh.nodesTags :
newTags.CreateTag(tag.name).SetIds(NewIndex[np.extract(usedNodes[tag.GetIds()],tag.GetIds() )])
mesh.nodesTags = newTags
#renumbering the connectivity matrix
for elementName in mesh.elements.keys():
elements = mesh.elements[elementName]
elements.connectivity = NewIndex[elements.connectivity]
for name,data in mesh.nodeFields.items():
mesh.nodeFields[name] = mesh.nodeFields[name][usedNodes]
return usedNodes, mesh
else:
res = Mesh()
res.nodes = mesh.nodes[usedNodes ,:]
res.originalIDNodes = originalIDNodes
#node tags
for tag in mesh.nodesTags :
outTag = res.GetNodalTag(tag.name)
outTag.SetIds(NewIndex[np.extract(usedNodes[tag.GetIds()],tag.GetIds() )])
for elementName in mesh.elements.keys():
elements = mesh.elements[elementName]
outElements = res.GetElementsOfType(elementName)
outElements.connectivity = NewIndex[elements.connectivity]
outElements.tags = elements.tags.Copy()
for name,data in mesh.nodeFields.items():
mesh.nodeFields[name] = mesh.nodeFields[name][usedNodes]
return usedNodes, res
[docs]def CleanDoubleElements(mesh: Mesh):
"""Remove double elements on the mesh (inPlace), even if a permutation exists
Tags of duplicated elements are merged into a single element.
Only element fields of the first element are kept.
Parameters
----------
mesh : Mesh
the input mesh
"""
newElementContainer = AllElements()
elemFieldMask = np.zeros((mesh.GetNumberOfElements(), ), dtype = MuscatIndex)
maskCpt = 0 # mask counter
elementCpt = 0 # element counter
for elemName, data in mesh.elements.items():
nbe = data.GetNumberOfElements()
if nbe == 0:
continue
#in the case of a StructuredElementsContainer we are sure we dont have double elements
if not isinstance(data, ElementsContainer):
newElementContainer.AddContainer(data)
elemFieldMask[maskCpt : maskCpt + nbe] = np.arange(nbe) + elementCpt
maskCpt += nbe
elementCpt += nbe
continue
data.Tighten()
_, indices, inverse = np.unique(np.sort(data.connectivity, axis = 1), return_index = True, return_inverse = True, axis = 0)
newElements = ElementsContainer(elemName)
newNbe = len(indices)
newElements.cpt = newNbe
newElements.connectivity = data.connectivity[indices, :]
newElements.originalIds = data.originalIds[indices]
elemFieldMask[maskCpt : maskCpt + newNbe] = indices + elementCpt
maskCpt += newNbe
elementCpt += nbe
for tag in data.tags:
newElements.tags.CreateTag(tag.name).SetIds(inverse[tag.GetIds()])
newElementContainer.AddContainer(newElements)
mesh.elements = newElementContainer
elemFieldMask = elemFieldMask[: maskCpt]
for name, field in mesh.elemFields.items():
if len(field.shape) == 1:
mesh.elemFields[name] = field[elemFieldMask]
else:
mesh.elemFields[name] = field[elemFieldMask, :]
[docs]def DeleteElements(mesh: Mesh, mask: ArrayLike, updateElementFields: bool= False):
"""Delete elements on the input mesh (inplace)
Parameters
----------
mesh : Mesh
the input mesh
mask : ArrayLike
the mash of elements to delete
updateElementFields: bool
True to update Element fields
"""
from Muscat.Containers.MeshInspectionTools import ExtractElementsByMask
OriginalNumberOfElements = mesh.GetNumberOfElements()
dataToChange = {}
offsets = mesh.ComputeGlobalOffset()
for name, data in mesh.elements.items():
offset = offsets[name]
localMask = mask[offset:offset+data.GetNumberOfElements()]
dataToChange[name] = ExtractElementsByMask(data, np.logical_not(localMask))
temp = Mesh()
temp.elements.storage = mesh.elements.storage[:]
temp.elemFields = dict(mesh.elemFields)
for name, data in dataToChange.items():
mesh.elements.AddContainer(data)
if OriginalNumberOfElements != mesh.GetNumberOfElements():
if updateElementFields:
from Muscat.Containers.MeshFieldOperations import CopyFieldsFromOriginalMeshToTargetMesh
CopyFieldsFromOriginalMeshToTargetMesh(temp, mesh)
else:
print("Number Of Elements Changed: Dropping elemFields")
mesh.elemFields = {}
for name, data in mesh.elements.items():
data.originalIds = temp.elements.GetElementsOfType(name).originalIds[data.originalIds]
mesh.PrepareForOutput()
[docs]def DeleteInternalFaces(mesh:Mesh):
"""Delete faces not present on the skin (internal faces) (inplace)
Parameters
----------
mesh : Mesh
the mesh
"""
from Muscat.Containers.MeshInspectionTools import ExtractElementsByMask
OriginalNumberOfElements = mesh.GetNumberOfElements()
skin = ComputeSkin(mesh)
dataToChange = {}
for name,data in mesh.elements.items():
if ED.dimensionality[name] != 2:
continue
ne = data.GetNumberOfElements()
mask = np.zeros(ne, dtype=bool)
data2 = skin.elements[name]
ne2 =data2.GetNumberOfElements()
surf2 = {}
key = np.array([ne**x for x in range(ED.numberOfNodes[name]) ])
for i in range(ne2):
cc = data2.connectivity[i,:]
lc = np.sort(cc)
elementHash = np.sum(lc*key)
surf2[elementHash] = [1,cc]
for i in range(ne):
cc = data.connectivity[i,:]
lc = np.sort(cc)
elementHash = np.sum(lc*key)
if elementHash in surf2:
mask[i] = True
if np.any(mask):
dataToChange[name] = ExtractElementsByMask(data,mask)
for name, data in dataToChange.items():
mesh.elements.AddContainer(data)
if OriginalNumberOfElements != mesh.GetNumberOfElements():
print("Number Of Elements Changed: Dropping elemFields")
mesh.elemFields = {}
mesh.PrepareForOutput()
[docs]def ComputeSkin(mesh:Mesh, md:Optional[int]=None, inPlace:bool=False, skinTagName: str="Skin")->Mesh:
"""Compute the skin of a mesh (mesh), if md (mesh dimensionality) is None
the mesh.GetPointsDimensionality() is used to filter the element to compute
the skin.
Warning if some elements are duplicated the behavior of computeSkin with
the option inPlace=True is not defined (Use CleanDoubleElements before
ComputeSkin)
Warning: element fields are not handled by this function and will be deleted
Parameters
----------
mesh : Mesh
the input mesh
md : Optional[int], optional
mesh dimensionality to extract element to compute the skin, by default None
inPlace : bool, optional, by default False
if True, the skin is added to the original mesh, a tag "Skin" is created.
if False a new mesh is returned with the elements of the skin (even if some element of the skin are already present on the original mesh).
If the user merged the two meshes, he must call CleanDoubleElements to clean the mesh after.
skinTagName: str, default "Skin"
name of the created eTag
Returns
-------
Mesh
The mesh containing the skin, (the mesh if inPlace == True, or only the skin if inPlace == False)
"""
if md is None:
md = mesh.GetPointsDimensionality()
# first we count the total number of potential individual skin elements
totalCpt = {}
DFilter = ElementFilter(dimensionality = md)
for selection in DFilter(mesh):
ids = selection.indices
faces = ED.faces1[selection.elementType]
nbElements = selection.elements.GetNumberOfElements()
for faceType,localFaceConnectivity in faces:
cpt = totalCpt.setdefault(faceType,0) + nbElements
totalCpt[faceType] = cpt
if inPlace :
# we add the element already present on the mesh
D1filter = ElementFilter(dimensionality = md-1)
for selection in D1filter(mesh):
nbElements = selection.elements.GetNumberOfElements()
cpt = totalCpt.setdefault(selection.elementType, 0) + nbElements
totalCpt[selection.elementType] = cpt
# Allocation of matrices to store all skin elements
storage = {}
for k,v in totalCpt.items():
nn = ED.numberOfNodes[k]
storage[k] = np.empty((v, nn), dtype = int)
totalCpt[k] = 0
if inPlace :
# fill storage with skin element already on the mesh
for selection in D1filter(mesh):
cpt = totalCpt[selection.elementType]
nbElements = selection.elements.GetNumberOfElements()
storage[selection.elementType][cpt : cpt + nbElements, :] = selection.elements.connectivity
totalCpt[selection.elementType] += nbElements
# fill storage with skin of elements of dimensionality D
for selection in DFilter(mesh):
elemName = selection.elements.elementType
data = selection.elements
ids = selection.indices
faces = ED.faces1[elemName]
nbElements = data.GetNumberOfElements()
for faceType,localFaceConnectivity in faces:
globalFaceConnectivity = data.connectivity[:, localFaceConnectivity]
cpt = totalCpt[faceType]
storage[faceType][cpt : cpt + nbElements, :] = globalFaceConnectivity
totalCpt[faceType] += nbElements
if inPlace:
res = mesh
if mesh.elemFields.keys():
Logger.Warning("Element fields are not handled by this function and will be deleted")
mesh.elemFields = {}
else:
res = Mesh()
res.nodesTags = mesh.nodesTags
res.originalIDNodes = mesh.originalIDNodes
res.nodes = mesh.nodes
res.elements = type(mesh.elements)()
# recover only the unique elements
for elemName, cpt in totalCpt.items():
if cpt == 0:
continue
store = storage[elemName]
_, index, counts = np.unique(np.sort(store, axis = 1), return_index = True, return_counts = True, axis = 0)
# the elements
if inPlace:
uniqueElements = index[counts == 2]
if elemName in mesh.elements:
meshSkinElements = mesh.elements.GetElementsOfType(elemName)
nbMElements = meshSkinElements.GetNumberOfElements()
ids = uniqueElements[uniqueElements < nbMElements]
# tag element already present on the original mesh
meshSkinElements.tags.CreateTag(skinTagName, False).SetIds(ids)
else:
nbMElements = 0
else:
nbMElements = 0
# all the elements present only 1 time not in the original mesh
uniqueElements = index[counts == 1]
ids = uniqueElements[uniqueElements >= nbMElements]
if len(ids) == 0 :
continue
if inPlace:
elementContainer = mesh.elements.GetElementsOfType(elemName)
else:
elementContainer = res.GetElementsOfType(elemName)
newElements = store[ids, :]
# add and tag new elements into the output
elementContainer.tags.CreateTag(skinTagName, False).AddToTag(np.arange(newElements.shape[0]) + elementContainer.GetNumberOfElements())
elementContainer.AddNewElements(newElements)
return res
[docs]def ComputeFeatures(inputMesh: Mesh, featureAngle: MuscatFloat = 90., skin: Optional[Mesh] = None) -> Tuple[Mesh, Mesh]:
"""Compute features, element of dimensionality dim-2 (edges) for a given angle
Parameters
----------
inputMesh : Mesh
the input mesh
featureAngle : MuscatFloat, optional
the detection angle: edges with faces with an angle smaller or equal than featureAngle will be consider as ridges.
We use the internal angle, by default 90.
skin : Optional[Mesh], optional
if the user can provide the skin (to save time), by default None
Returns
-------
Tuple[Mesh,Mesh]
edgeMesh: a mesh containing the edges
skinMesh: a mesh containing the skin
"""
if skin is None:
skinMesh = ComputeSkin(inputMesh)
# we have to merge all the 2D elements from the original mesh to the skin mesh
for selection in ElementFilter(dimensionality = 2)(inputMesh):
skinMesh.GetElementsOfType(selection.elementType).Merge(selection.elements)
CleanDoubleElements(skinMesh)
else:
skinMesh = skin
skinMeshSave = copy.deepcopy(skinMesh)
md = skinMesh.GetPointsDimensionality()
nex = skinMesh.GetNumberOfElements()
# we use the original id to count the number of time the faces is used
surf = {}
for name, data in skinMesh.elements.items():
if ED.dimensionality[name] != md-1:
continue
faces = ED.faces1[name]
numberOfNodes = ED.numberOfNodes[name]
ne = data.GetNumberOfElements()
for faceType, localFaceConnectivity in faces:
globalFaceConnectivity = data.connectivity[:, localFaceConnectivity]
if not faceType in surf:
surf[faceType] = {}
surf2 = surf[faceType]
key = np.array([nex**x for x in range(ED.numberOfNodes[faceType])])
for i in range(ne):
bariCentre = np.sum(skinMesh.nodes[data.connectivity[i, :], :], axis = 0)/numberOfNodes
cc = globalFaceConnectivity[i, :]
lc = np.sort(cc)
elementHash = np.sum(lc*key)
edgeVector = skinMesh.nodes[lc[0], :] - skinMesh.nodes[lc[1], :]
planeVector = bariCentre - skinMesh.nodes[lc[1], :]
normal = np.cross(edgeVector, planeVector)
normal /= np.linalg.norm(normal)
if elementHash in surf2:
surf2[elementHash][0] +=1
normal1 = surf2[elementHash][1]
dot = normal.dot(normal1)
dot = dot if dot < 1 else 1
dot = dot if dot > -1 else -1
angle = np.arccos(dot)
surf2[elementHash][2] = 180*angle/np.pi
else:
surf2[elementHash] = [1, normal, None, cc]
edgeMesh = Mesh()
edgeMesh.nodes = inputMesh.nodes
edgeMesh.nodesTags = inputMesh.nodesTags
edgeMesh.originalIDNodes = inputMesh.originalIDNodes
numberOfNodes = inputMesh.GetNumberOfNodes()
for name, surf2 in surf.items():
if len(surf2) == 0:
continue
data = edgeMesh.GetElementsOfType(name)
for data2 in surf2.values():
if data2[0] == 1 or data2[0] > 2:
data.AddNewElement(data2[3], -1)
elif data2[0] == 2 and data2[2] <= featureAngle:
data.AddNewElement(data2[3], -1)
for elType in [ED.Bar_2, ED.Bar_3]:
if elType not in edgeMesh.elements:
continue
bars = edgeMesh.GetElementsOfType(elType)
bars.tags.CreateTag("Ridges").SetIds(np.arange(bars.GetNumberOfElements()))
# Now we compute the corners
#The corner are the points:
# 1) where 3 Ridges meet
# 2) or touched by only one Ridge
# 3) or the angle of 2 ridges is bigger than the featureAngle
extractedBars = edgeMesh.GetElementsOfType(ED.Bar_2)
extractedBars.Tighten()
originalBars = inputMesh.GetElementsOfType(ED.Bar_2)
originalBars.Tighten()
#corners
mask = np.zeros((numberOfNodes), dtype = bool)
almanac = {}
for bars in [originalBars, extractedBars]:
intMask = np.zeros((numberOfNodes), dtype = np.int8)
np.add.at(intMask, bars.connectivity.ravel(), 1)
mask[intMask > 2 ] = True
mask[intMask == 1 ] = True
for bar in range(bars.GetNumberOfElements()):
id0, id1 = bars.connectivity[bar, :]
p0 = inputMesh.nodes[id0, :]
p1 = inputMesh.nodes[id1, :]
vec1 = (p1-p0)
vec1 /= np.linalg.norm(vec1)
for p, sign in zip([id0, id1], [1, -1]):
if mask[p]:
#already in the mask no need to treat this point
continue
if p in almanac:
vec2 = almanac[p][1]
norm = np.dot(vec1, vec2*sign)
norm = norm if norm < 1 else 1
norm = norm if norm > -1 else -1
angle = (180/np.pi)*np.arccos(norm)
almanac[p][2] = angle
almanac[p][0] +=1
else:
almanac[p] = [1, vec1*sign, id0]
for p,data in almanac.items():
if (180-data[2]) >= featureAngle:
mask[p] = True
edgeMesh.nodesTags.CreateTag("Corners", False).AddToTag(np.where(mask)[0])
edgeMesh.PrepareForOutput()
skinMesh.PrepareForOutput()
return (edgeMesh,skinMeshSave)
[docs]def NodesPermutation(mesh: Mesh, per:ArrayLike):
"""Function to do a permutation of the nodes in a mesh (in place)
Parameters
----------
mesh : Mesh
mesh to be modified
per : ArrayLike
the permutation vector ([1,0,2,3] to permute first an second node)
"""
newNodes = mesh.nodes[per,:]
perII =np.argsort(per).astype(MuscatIndex,casting='same_kind',copy=False)
mesh.nodes = newNodes
for tag in mesh.nodesTags:
ids = tag.GetIds()
newIds = perII[ids]
tag.SetIds(newIds)
for data in mesh.elements.values():
data.connectivity = perII[data.connectivity]
for name,data in mesh.nodeFields.items():
mesh.nodeFields[name] = mesh.nodeFields[name][per]
[docs]def AddTagPerBody(inmesh:Mesh)-> np.ndarray:
"""Generate nodal tag (in the form of `Body_+int`) per body
a body is defined by all the nodes connected by the elements (connectivity filter in vtk )
Parameters
----------
inmesh : Mesh
the input mesh
Returns
-------
np.ndarray
pointsPerBody : a vector with the number of point in every body.
len(pointsPerBody) is the number of unconnected bodies in the mesh
Raises
------
Exception
in the case of an internal error
"""
from Muscat.Containers.MeshInspectionTools import ComputeNodeToNodeConnectivity
dualGraph,usedPoints = ComputeNodeToNodeConnectivity(inmesh)
# Connectivity walk
nbOfNodes = inmesh.GetNumberOfNodes()
treated = np.zeros(inmesh.GetNumberOfNodes(), dtype=bool)
nextPoint = np.zeros(inmesh.GetNumberOfNodes(), dtype=MuscatIndex)
# we start from the first point and body number 0
nextPointCpt = 0
bodyCpt = 0
cpt = 0
pointsPerBody = []
while True:
# we finish all the points
if cpt == nbOfNodes :
break
# we already explored all the point in this body
if cpt == nextPointCpt:
initialPoint = np.argmax(treated == False)
#(edge case) in the case we have only Trues in the treated
if treated[initialPoint]:# pragma: no cover
break
treated[initialPoint] = True
nextPoint[nextPointCpt] = initialPoint
nextPointCpt +=1
tagName = "Body_"+str(bodyCpt)
tag = inmesh.GetNodalTag(tagName)
tag.AddToTag(initialPoint)
bodyCpt += 1
pointsInThisBody = 1
else:
raise Exception ("Error in this function")# pragma: no cover
while cpt < nextPointCpt:
workingIndex = nextPoint[cpt]
indexes = dualGraph[workingIndex,0:usedPoints[workingIndex]]
for index in indexes:
if not treated[index] :
treated[index] = True
nextPoint[nextPointCpt] = index
nextPointCpt +=1
tag.AddToTag(index)
pointsInThisBody +=1
cpt += 1
pointsPerBody.append(pointsInThisBody)
return pointsPerBody
[docs]def Morphing(mesh: Mesh, targetDisplacement: np.ndarray, targetDisplacementMask: np.ndarray, radius:Optional[MuscatFloat]=None, max_steps:int = 1, verbose:bool = False)-> np.ndarray:
"""method for computing the deform mesh knowing displacement of some nodes.
the user can push the morphed point back to the mesh by doing : mesh.node = morphedPoints
note: https://www.researchgate.net/publication/288624175_Mesh_deformation_based_on_radial_basis_function_interpolation_computers_and_structure
Parameters
----------
mesh : Mesh
the input mesh, the use only the point position information
targetDisplacement : numpy array
numpy array: the known displacement of control nodes (shape [number_of_of_known_nodes,3])
targetDisplacementMask : numpy array
numpy array: the ids of control nodes (list of ids or boolean array) in the same order as targetDisplacement
radius : Optional[MuscatFloat], optional
you can choose a radius by setting radius to a value, by default np.linalg.norm(boundingMax-boundingMin)/2
max_steps : int, optional
default: 1; number of max steps for decomposing the morphing
verbose : bool, optional
default: False; print additional information when True
Returns
-------
np.ndarray
morphedPoints: the morphed points
"""
if verbose:
print(f"RBF morphing starts with {targetDisplacementMask.shape[0]} control points ...")
grad_max=0.8
step=0
nb_steps=1
new_nodes = np.copy(mesh.GetPosOfNodes())
if radius is None:
boundingMin, boundingMax = mesh.ComputeBoundingBox()
r=np.linalg.norm(boundingMax-boundingMin)/2
else:
assert radius > 0, "radius cannot be negative or zero"
r=radius
while step<nb_steps:
border_nodes=new_nodes[targetDisplacementMask,:]
rhs=targetDisplacement/nb_steps
d_r = distance.cdist(border_nodes, border_nodes)/r
M = (1-d_r)**4*(4*d_r+1)
M[d_r>1]=0.
np.fill_diagonal(M,1.)
if verbose:
start = time.time()
try:
alpha = np.linalg.solve(M,rhs)
except np.linalg.LinAlgError:
print("ill-conditionned operator, trying a least squares")
alpha = np.linalg.lstsq(M,rhs, rcond=10**(9))[0]
if verbose:
print(f"time to solve to solve linear system: {round(time.time()-start, 2)} s")
start = time.time()
d_r = distance.cdist(new_nodes, border_nodes)/r
y=(1-d_r)**4*(4*d_r+1)
inds = d_r>=1.
y[inds] = 0.
s = np.dot(y, alpha)
if step==0 and max_steps>1:
y=-20*d_r*(1-d_r)**3
y[inds] = 0.
ds = np.dot(y, alpha)/r
nb_steps=min(max_steps, int(np.floor(np.max(np.linalg.norm(ds,axis=1)/grad_max)))+1)
s=s/nb_steps
new_nodes += s
step += 1
if verbose:
print(f"time to update node positions: {round(time.time()-start, 2)} s")
if verbose:
print("RBF morphing ended")
return new_nodes
[docs]def LowerNodesDimension(mesh:Mesh) -> Mesh:
"""Remove the last columns of the mesh nodes coordinate.
This is done in place
Parameters
----------
mesh : Mesh
the input mesh
Returns
-------
Mesh
the mesh
"""
newDim = mesh.GetPointsDimensionality() - 1
mesh.nodes = np.ascontiguousarray(mesh.nodes[:,0:newDim])
return mesh
[docs]def RigidBodyTransformation(mesh: Mesh, rotationMatrix: np.ndarray, translationVector: np.ndarray):
"""In place rigid body transformation.
new pos = Q.pos + translation
the rotation matrix Q should verify:
QtQ = I = QQt et det Q = 1
Parameters
----------
mesh : Mesh
The input mesh
rotationMatrix : np.ndArray
a 3x3 rotation matrix
translationVector : np.ndArray
a vector of size 3 with the translation
"""
assert np.linalg.norm(np.dot(rotationMatrix.T, rotationMatrix) - np.eye(3)) < 1.e-12
assert np.linalg.norm(np.dot(rotationMatrix, rotationMatrix.T) - np.eye(3)) < 1.e-12
assert (np.linalg.det(rotationMatrix) - 1) < 1.e-12
mesh.nodes = np.dot(rotationMatrix, mesh.nodes.T).T
for i in range(3):
mesh.nodes[:,i] += translationVector[i]
[docs]def ComputeRigidBodyTransformationBetweenTwoSetOfPoints(setPoints1: np.ndarray, setPoints2: np.ndarray)-> Tuple[np.ndarray, np.ndarray]:
""" Compute the rotation and the translation operator from two sets of points
setPoints1 and setPoints2 have the same dimension (nbeOfPoints,dimension)
Parameters
----------
setPoints1 : np.ndarray
First set of points
setPoints2 : np.ndarray
Second set of points
Returns
-------
Tuple[np.ndarray, np.ndarray]
A, b such that setPoints1.T approx A*setPoints2.T + b
"""
assert setPoints1.shape == setPoints2.shape
dim = setPoints1.shape[1]
setPoints1 = np.hstack((setPoints1, np.ones((setPoints1.shape[0],1))))
res = np.linalg.lstsq(setPoints1, setPoints2, rcond=None)[0]
A = res[:dim,:dim].T
b = res[dim,:].T
return A, b
#------------------------- CheckIntegrity ------------------------
[docs]def CheckIntegrity_AddTagPerBody(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles
from Muscat.Containers.MeshCreationTools import MirrorMesh
res = CreateMeshOfTriangles([[0,0,0],[1,0,0],[0,1,0],[0,0,1] ], [[0,1,2],[0,2,3]])
resII = MirrorMesh(res,x=0,y=0,z=0)
print( resII.nodes)
print( resII.GetElementsOfType(ED.Triangle_3))
print(AddTagPerBody(resII))
if len(resII.nodesTags) != 8: # pragma: no cover
raise # pragma: no cover
print(resII.nodesTags)
return "ok"
[docs]def CheckIntegrity_CleanDoubleNodes(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0] ]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
CleanDoubleNodes(mesh)
if mesh.GetNumberOfNodes() != 4: # pragma: no cover
raise Exception(f"{mesh.GetNumberOfNodes()} != 4")# pragma: no cover
#test non double nodes
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1]]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
CleanDoubleNodes(mesh)
if mesh.GetNumberOfNodes() != 4: # pragma: no cover
raise Exception(f"{mesh.GetNumberOfNodes()} != 4")# pragma: no cover
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,0,0] ]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("OnePoint").SetIds([0,1])
CleanDoubleNodes(mesh,tol =0)
print(mesh.nodes)
if mesh.GetNumberOfNodes() != 4: # pragma: no cover
raise Exception(f"{mesh.GetNumberOfNodes()} != 4")# pragma: no cover
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[0,0,0] ]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
CleanDoubleNodes(mesh, nodesToTestMask= np.array([True,False,True,False,True], dtype=bool) )
if mesh.GetNumberOfNodes() != 4: # pragma: no cover
raise Exception(f"{mesh.GetNumberOfNodes()} != 4")# pragma: no cover
return "ok"
[docs]def CheckIntegrity_CleanLonelyNodes(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("OnePoint").SetIds([0,1])
mesh.nodeFields["Scalar Field"] = np.arange(mesh.GetNumberOfNodes())
CleanLonelyNodes(mesh)
if mesh.GetNumberOfNodes() != 4:
raise# pragma: no cover
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("OnePoint").SetIds([0,1])
mesh.nodeFields["Vector Field"] = np.zeros((mesh.GetNumberOfNodes(),3))
usedNodes, out = CleanLonelyNodes(mesh,inPlace=False)
if out.GetNumberOfNodes() != 4: # pragma: no cover
print(out)
raise# pragma: no cover
return "ok"
[docs]def CheckIntegrity_ComputeFeatures(GUI:bool =False):
from Muscat.Containers.MeshCreationTools import CreateConstantRectilinearMesh
res2 = CreateConstantRectilinearMesh(dimensions = [2, 3, 4], origin = [-1.0, -1.0, -1.0], spacing = [2./2, 2./3, 2/4])
from Muscat.Containers.MeshTetrahedrization import Tetrahedrization
res2 = Tetrahedrization(res2)
print("Mesh")
print(res2)
edges, skin = ComputeFeatures(res2, featureAngle = 80)
print("edges mesh")
print(edges)
print("SkinMesh")
print(skin)
if GUI : # pragma: no cover
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh = edges, filename = "edges.xmf")
OpenInParaView(mesh = skin, filename = "skin.xmf")
for name,data in edges.elements.items():
res2.GetElementsOfType(name).Merge(data)
OpenInParaView(res2, filename = "all+edges.xmf")
print(res2)
edges, skin = ComputeFeatures(res2, featureAngle = 80, skin = skin)
print(edges)
return "ok"
[docs]def CheckIntegrity_ComputeSkin(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateConstantRectilinearMesh
from Muscat.Containers.MeshTetrahedrization import Tetrahedrization
res2 = CreateConstantRectilinearMesh(dimensions = [2, 3, 4], origin = [-1.0, -1.0, -1.0], spacing = [2/2, 2/3, 2/4])
res2 = Tetrahedrization(res2)
skin = ComputeSkin(res2, inPlace = False)
print(skin)
if GUI : # pragma: no cover
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(skin,filename="CheckIntegrity_ComputeSkin_InPlace_False.xmf")
print(res2)
res2.MergeElements(skin)
ComputeSkin(res2, inPlace = True)
print(res2)
if GUI : # pragma: no cover
OpenInParaView(res2,filename="CheckIntegrity_ComputeSkin_InPlace_True.xmf")
return "ok"
[docs]def CheckIntegrity_Morphing(GUI:bool = False):
from Muscat.Helpers.CheckTools import MustFailFunction
from Muscat.Containers.MeshCreationTools import CreateCube
mesh = CreateCube(dimensions=[20,21,22],spacing=[2.,2.,2.],ofTetras=True)
targetDisplacement = np.empty((mesh.GetNumberOfNodes(),3),dtype=float)
targetDisplacementMask = np.empty(mesh.GetNumberOfNodes(),dtype=int)
cpt = 0
print(mesh)
for name,data in mesh.elements.items():
if ED.dimensionality[name] != 2:
continue
ids = data.GetNodesIndexFor(data.GetTag("X0").GetIds())
print(ids)
targetDisplacementMask[cpt:cpt+len(ids)] = ids
targetDisplacement[cpt:cpt+len(ids),:] = 0
cpt += len(ids)
ids = data.GetNodesIndexFor(data.GetTag("X1").GetIds())
print(ids)
targetDisplacementMask[cpt:cpt+len(ids)] = ids
targetDisplacement[cpt:cpt+len(ids),:] = [[0,0,10]]
cpt += len(ids)
targetDisplacement = targetDisplacement[0:cpt,:]
targetDisplacementMask = targetDisplacementMask[0:cpt]
new_p1 = Morphing(mesh, targetDisplacement,targetDisplacementMask)
new_p2 = Morphing(mesh, targetDisplacement,targetDisplacementMask, radius = 20.)
MustFailFunction(Morphing, mesh, targetDisplacement,targetDisplacementMask, radius = 0.)
new_p2 = Morphing(mesh, targetDisplacement,targetDisplacementMask, radius = 1., max_steps = 1, verbose = True)
new_p0 = np.copy(mesh.nodes)
new_p0[targetDisplacementMask,:] += targetDisplacement
mesh.nodeFields["morph0"] = new_p0
mesh.nodeFields["morph1"] = new_p1
mesh.nodeFields["morph2"] = new_p2
if GUI : # pragma: no cover
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh=mesh)
return "ok"
[docs]def CheckIntegrity_NodesPermutation(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[1,2,3],[4,5,6],[0,1,0],[0,0,1] ]
tets = [[0,1,2,3],[3,1,0,2]]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("firstPoint").AddToTag(0)
mesh.nodeFields["nodal field"] = np.arange(4)+10
mesh.nodeFields["nodal field tensor"] = np.zeros((4,3),dtype=MuscatIndex)
mesh.nodeFields["nodal field tensor"][0,:] = [2,3,4]
NodesPermutation(mesh, [1,0,2,3])
print( mesh.nodes[0,:] )
if np.any( mesh.nodes[0,:]-[4,5,6] ): # pragma: no cover
raise(Exception("error in the permutation"))
np.testing.assert_equal(mesh.nodesTags["firstPoint"].GetIds() ,[1])
np.testing.assert_equal(mesh.nodeFields["nodal field"],[11,10,12,13])
np.testing.assert_equal(mesh.nodeFields["nodal field tensor"][1,1],3)
return "ok"
[docs]def CheckIntegrity_DeleteInternalFaces(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],[3,0,2,4]]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
## add 2 triangles
triangles = mesh.GetElementsOfType(ED.Triangle_3)
triangles.AddNewElement([0,1,2],0)
triangles.AddNewElement([3,0,2],1)
triangles.AddNewElement([3,0,2],2)
#print(mesh)
DeleteInternalFaces(mesh)
print(mesh)
return "ok"
[docs]def CheckIntegrity_RigidBodyTransformation(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],[3,0,2,4]]
mesh = CreateMeshOf(points,tets,ED.Tetrahedron_4)
## add 2 triangles
triangles = mesh.GetElementsOfType(ED.Triangle_3)
triangles.AddNewElement([0,1,2],0)
triangles.AddNewElement([3,0,2],1)
triangles.AddNewElement([3,0,2],2)
#print(mesh)
theta = np.pi/2.
A = np.array([[1, 0, 0],
[0, np.cos(theta), -np.sin(theta)],
[0, np.sin(theta), np.cos(theta)]])
b = np.array([1,0,0])
RigidBodyTransformation(mesh, A, b)
return "ok"
[docs]def CheckIntegrity_ComputeRigidBodyTransformationBetweenTwoSetOfPoints(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points1 = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],[3,0,2,4]]
mesh1 = CreateMeshOf(points1,tets,ED.Tetrahedron_4)
mesh2 = CreateMeshOf(points1,tets,ED.Tetrahedron_4)
theta = np.pi/3.
A = np.array([[1, 0, 0],
[0, np.cos(theta), -np.sin(theta)],
[0, np.sin(theta), np.cos(theta)]])
b = np.array([1,0,0])
RigidBodyTransformation(mesh2, A, b)
A, b = ComputeRigidBodyTransformationBetweenTwoSetOfPoints(mesh1.nodes, mesh1.nodes)
return "ok"
[docs]def CheckIntegrity_DeleteElements(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points1 = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2],[1,2,3],[2,3,4]]
mesh1 = CreateMeshOf(points1,tets,ED.Triangle_3)
mesh1.elements[ED.Triangle_3].originalIds += 10
mask = np.zeros(3)
mask[[0,2] ] = 1
mesh1.elemFields["testField"] = np.arange(mesh1.GetNumberOfElements() , dtype=MuscatFloat) +0.1
DeleteElements(mesh1, mask , True)
assert( len(mesh1.elemFields["testField"]) ==1 )
assert( mesh1.elemFields["testField"][0] == 1.1 )
print("toto",mesh1.elemFields["testField"] )
assert( len(mesh1.elements[ED.Triangle_3].originalIds) == 1 )
assert( mesh1.elements[ED.Triangle_3].originalIds[0] == 11 )
return "ok"
[docs]def CheckIntegrity_CleanDoubleElements(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points1 = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 1]]
tets = [[0, 1, 2, 3], [3, 0, 2, 4], [3, 0, 2, 4], [1, 0, 2, 4]]
mesh1 = CreateMeshOf(points1,tets,ED.Tetrahedron_4)
mesh1.GetElementsOfType(ED.Tetrahedron_4).tags.CreateTag("tag1", False).SetIds([0, 1])
mesh1.GetElementsOfType(ED.Tetrahedron_4).tags.CreateTag("tag2", False).SetIds([2, 3])
mesh1.GetElementsOfType(ED.Tetrahedron_4).originalIds = np.array([5, 6, 7, 8])
mesh1.GetElementsOfType(ED.Triangle_3)
mesh1.elemFields["testField"] = np.arange(4)
mesh1.elemFields["testField4x3"] = np.zeros((4, 3))
CleanDoubleElements(mesh1)
if mesh1.GetNumberOfElements() != 3: # pragma: no cover
raise Exception("ko")
print(mesh1.GetElementsOfType(ED.Tetrahedron_4).originalIds)
print(mesh1)
np.testing.assert_equal(mesh1.GetElementsOfType(ED.Tetrahedron_4).originalIds, [5, 8, 6])
assert mesh1.GetElementsOfType(ED.Tetrahedron_4).tags["tag1"].GetIds()[1] == 2
assert mesh1.GetElementsOfType(ED.Tetrahedron_4).tags["tag2"].GetIds()[1] == 2
return "ok"
[docs]def CheckIntegrity_LowerNodesDimension(GUI:bool=False):
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points1 = [[0,0,0],[1,0,0],[0,1,0],[0,0,1],[1,1,1] ]
tets = [[0,1,2,3],[3,0,2,4],[3,0,2,4], [1,0,2,4] ]
mesh1 = CreateMeshOf(points1,tets,ED.Tetrahedron_4)
LowerNodesDimension(mesh1)
if mesh1.nodes.shape[1] != 2: # pragma: no cover
raise Exception("ko")
return "ok"
[docs]def CheckIntegrity(GUI:bool=False):
functionsToTest= [
CheckIntegrity_NodesPermutation,
CheckIntegrity_Morphing,
CheckIntegrity_ComputeFeatures,
CheckIntegrity_ComputeSkin,
CheckIntegrity_DeleteInternalFaces,
CheckIntegrity_CleanDoubleNodes,
CheckIntegrity_CleanLonelyNodes,
CheckIntegrity_AddTagPerBody,
CheckIntegrity_RigidBodyTransformation,
CheckIntegrity_ComputeRigidBodyTransformationBetweenTwoSetOfPoints,
CheckIntegrity_CleanDoubleElements,
CheckIntegrity_CopyElementTags,
CheckIntegrity_DeleteElements,
CheckIntegrity_ConvertNTagsToETags,
CheckIntegrity_LowerNodesDimension
]
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
return RunListOfCheckIntegrities(functionsToTest)
if __name__ == '__main__':
print(CheckIntegrity(False))# pragma: no cover