# -*- 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.
#
from typing import Callable, Union, Tuple
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.Filters.FilterObjects import ElementFilter, NodeFilter, NM_AT_LEAST_ONE
from Muscat.Containers.Filters.FilterOperators import IntersectionFilter, DifferenceFilter
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.MeshFieldOperations import TransportPosToPoints
from Muscat.Containers.MeshTetrahedrization import Tetrahedrization, TestMeshValidity, CopyElements
from Muscat.Containers.MeshModificationTools import CleanDoubleElements
[docs]def DistanceToSurface(mesh, surfMesh, out = None, method = "Interp/Clamp"):
from Muscat.Containers.MeshFieldOperations import TransportPos
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.FE.Fields.FEField import FEField
tspace, tnumbering, _, _ = PrepareFEComputation(mesh, numberOfComponents = 1)
tnumbering = tnumbering[0]
names = ["x", "y", "z"][0: mesh.GetPointsDimensionality()]
InterfacePosFields = TransportPos(surfMesh, mesh, tspace, tnumbering, method = method)
MeshPosFields = np.array([FEField("pos_"+names[x], mesh, space = tspace, numbering = tnumbering, data = mesh.nodes[:, x]) for x in range(len(names))])
if out is None:
res = np.sqrt(np.sum((InterfacePosFields - MeshPosFields)**2)).data
return res
else:
out[:] = np.sqrt(np.sum((InterfacePosFields - MeshPosFields)**2)).data
return out
[docs]def Redistance(mesh, phi, method = "Interp/Extrap"):
sign = np.sign(phi)
res = ComputeDistanceToIsoZero(mesh, phi, mesh.nodes, method = method)
res *= sign
return res
[docs]def ComputeDistanceToIsoZero(mesh, phi, targetPoints, method = "Interp/Extrap"):
IGTM = IGToMesh(mesh, phi)
newMesh = IGTM.ComputeInterfaceMesh()
pos = TransportPosToPoints(newMesh, targetPoints, method = method)
return np.sqrt(np.sum((targetPoints - pos)**2, axis = 1))
[docs]class IGToMesh:
def __init__(self, imesh = None, phi = None):
self.inputMesh = imesh
self.SetPhi(phi)
self.volMesh = None
[docs] def SetPhi(self, phi, snapValues = True, snapTol = 1e-5):
if phi is None:
self.phi = None
return
self.phi = np.copy(phi)
phimin = min(self.phi)
phimax = max(self.phi)
self.phi[abs(self.phi) < abs(phimax-phimin)*snapTol] = 0.
[docs] def ComputeInterfaceMesh(self):
from scipy.spatial import Delaunay
inputNodes = self.inputMesh.nodes
self.volMesh = Mesh()
self.volMesh.CopyProperties(self.inputMesh)
phi = self.phi
if np.max(phi) < 0 or np.min(phi) > 0:
print("Warning: non iso zero on phi")
pid = dict()
elems = set()
cpt = 0
vcpt = 0
xyz = []
vxyz = []
omesh = Mesh()
ef1 = ElementFilter(zone = [lambda x: phi <= 0])
ef2 = ElementFilter(zone = [lambda x: -phi <= 0])
ef1.zonesTreatment = "AT_LEAST_ONE"
ef2.zonesTreatment = "AT_LEAST_ONE"
ef = IntersectionFilter([ef1, ef2])
def AddElement(pos_id, neg_id, td, p):
if pos_id is not None:
if np.cross(xyz[p[2]]-xyz[p[1]], xyz[p[0]]-xyz[p[1]]).dot(inputNodes[pos_id, :]-pcontrol_point) > 0:
td.AddNewElement(p, 0)
else:
td.AddNewElement(np.flip(p), 0)
elif neg_id is not None:
if np.cross(xyz[p[2]]-xyz[p[1]], xyz[p[0]]-xyz[p[1]]).dot(inputNodes[neg_id, :]-ncontrol_point) > 0:
td.AddNewElement(np.flip(p), 0)
else:
td.AddNewElement(p, 0)
else:
raise
for selection in ef(self.inputMesh):
if ED.dimensionality[selection.elementType] == 1:
# continue
pointElements = omesh.elements.GetElementsOfType(ED.Point_1)
for eid in selection.indices:
lconn = selection.elements.connectivity[eid, :]
phi0 = self.phi[lconn[0]]
phi1 = self.phi[lconn[1]]
x = phi0/(phi0-phi1)
xpos = self.inputMesh.nodes[lconn[0], :]*(1-x) + self.inputMesh.nodes[lconn[1], :]*(x)
xyz.append(xpos)
cpt += 1
pointElements.AddNewElement(len(xyz)-1, 0)
continue
barElements = omesh.elements.GetElementsOfType(ED.Bar_2)
trisElements = omesh.elements.GetElementsOfType(ED.Triangle_3)
quadElements = omesh.elements.GetElementsOfType(ED.Quadrangle_4)
if ED.dimensionality[selection.elementType] == 2:
faces = ED.faces1[selection.elementType]
else:
faces = ED.faces2[selection.elementType]
for eid in selection.indices:
cutpoints = []
neg_pointsO = []
pos_pointsO = []
neg_pointsG = []
pos_pointsG = []
lconn = selection.elements.connectivity[eid, :]
# if zero on the a point
neg_id = None
pos_id = None
for n in lconn:
if phi[n] == 0:
if n not in pid:
nid = cpt
xyz.append(inputNodes[n, :])
pid[n] = cpt
cpt += 1
else:
nid = pid[n]
cutpoints.append(nid)
neg_pointsO.append(n)
pos_pointsO.append(n)
elif phi[n] < 0:
neg_id = n
neg_pointsO.append(n)
elif phi[n] > 0:
pos_id = n
pos_pointsO.append(n)
# if zero on a bars
for n, d in faces:
llcon = lconn[d[0:2]]
lphi = phi[llcon]
pphi = lphi > 0
nphi = lphi < 0
if np.any(pphi) and np.any(nphi):
t = tuple(np.sort(llcon))
if t in pid:
nid, vnid = pid[t]
else:
a0 = lphi[1]/(lphi[1] - lphi[0])
a1 = 1-a0
nid = cpt
vnid = vcpt
newpoint = inputNodes[llcon[1], :]*a1 + inputNodes[llcon[0], :]*a0
xyz.append(newpoint)
vxyz.append(newpoint)
pid[t] = (cpt, vcpt)
cpt += 1
vcpt += 1
cutpoints.append(nid)
neg_pointsG.append(vnid)
pos_pointsG.append(vnid)
if ED.dimensionality[selection.elementType] == 2:
if len(cutpoints) > 2:
raise (Exception("More than 2 point of intersection"))
if len(cutpoints) < 2:
continue
barElements.AddNewElement(cutpoints, 0)
continue
if len(cutpoints) < 3:
continue
scp = tuple(np.sort(cutpoints))
if scp in elems:
continue
elems.update(scp)
X = np.array([xyz[p] for p in cutpoints])
# ---# X is [n_pts x 3]
X_mean = X.mean(axis=0, keepdims=True) # [1 x 3]
X_centered = X - X_mean # [n_pts x 3]
valp, vecp = np.linalg.eigh(np.matmul(X_centered.T, X_centered))
assert (valp[0] <= valp[1] <= valp[2])
plane_base = vecp[:, 1:] # [3 x 2]
XX = np.matmul(X_centered, plane_base) # [n_pts x 2]
def pca_transform(x):
return np.matmul(x - X_mean, plane_base)
def pca_inverse_transform(x):
return np.matmul(x, plane_base.T) + X_mean
# use this to check the orientation
if pos_id is not None:
projected_point = pca_transform(inputNodes[[pos_id], :]) # [1 x 2]
pcontrol_point = pca_inverse_transform(projected_point)[0] # [1 x 3]
elif neg_id is not None:
projected_point = pca_transform(inputNodes[[neg_id], :]) # [1 x 2]
ncontrol_point = pca_inverse_transform(projected_point)[0] # [1 x 3]
if len(cutpoints) == 3 or len(cutpoints) == 4:
ang = [np.arctan2(x, y) for x, y in XX]
s = np.argsort(ang)
cutpoints2 = [cutpoints[x] for x in s]
if len(cutpoints) == 3:
AddElement(pos_id, neg_id, trisElements, cutpoints2)
else:
AddElement(pos_id, neg_id, quadElements, cutpoints2)
else:
tri = Delaunay(XX)
for simple in tri.simplices:
p = [cutpoints[x] for x in simple]
AddElement(pos_id, neg_id, trisElements, p)
omesh.nodes = np.array(xyz, dtype = MuscatFloat, ndmin = 2)
omesh.GenerateManufacturedOriginalIDs()
omesh.PrepareForOutput()
return omesh
[docs]def CutMeshIsoZero(mesh: Mesh, isoZeroFunction :Union[Callable[[ArrayLike], ArrayLike], np.ndarray], snapTolerance = None, computePointFieldTransferOperator = False) -> Tuple[Mesh,np.ndarray ]:
"""Cut a mesh using the iso zero computed by the function isoZeroFunction.
This function works for meshes composed only by Point_1, Bar_2, Triangle_3 and Tetrahedron_4
Parameters
----------
mesh : Mesh
Mesh to cut by the iso zero
isoZeroFunction : Callable[[ArrayLike], np.ndarray] or ArrayLike
if Callable Function to compute the level-set field.
if arrayLike iso at the nodes
snapTolerance : float,
tolerance to apply the snap of the izoZero to the nodes in absolute values (zero == no snap, default 1e-6 of the bounding box )
computePointFieldTransferOperator : bool, default False
if True a transfer operator for a point field is constructed
Returns
-------
Mesh
a new simplex mesh cut by the iso zeros value of the isoZeroFunction
pointFieldTransferOperator
None if computePointFieldTransferOperator is False
Raises
------
RuntimeError
if iso zero not found
RuntimeError
if elements of non supported types are present
"""
ptf_values = []
ptf_i = []
ptf_j = []
if callable(isoZeroFunction):
izoAtPoints = isoZeroFunction(mesh.nodes)
else:
izoAtPoints = np.asarray(isoZeroFunction, dtype = MuscatFloat)
if snapTolerance is None:
minbb, maxbb = mesh.ComputeBoundingBox()
snapTolerance = np.linalg.norm(maxbb-minbb)*1e-6
izoAtPoints = np.copy(izoAtPoints)
if snapTolerance > 0:
izoAtPoints[abs(izoAtPoints) < snapTolerance] = 0.0
filterP = ElementFilter(nMask = (izoAtPoints > 0), nMaskTreatment = NM_AT_LEAST_ONE)
filterN = ElementFilter(nMask = (izoAtPoints < 0), nMaskTreatment = NM_AT_LEAST_ONE)
filterInterface1 = IntersectionFilter([filterP, filterN]) # elements that will be cut by the isozero
from Muscat.Containers.Filters.FilterTools import GetComplementaryFilter
compFilter = GetComplementaryFilter(filterInterface1)
filter0 = ElementFilter(nMask = (izoAtPoints == 0), nMaskTreatment = NM_AT_LEAST_ONE)
filterInterface2 = IntersectionFilter([compFilter, filter0]) # elements that touch isozero but remain in the negative or positive side
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
interfaceMesh = ExtractElementsByElementFilter(inmesh = mesh, elementFilter = filterInterface2, copy = True)
# list of points to be added to the mesh
nbNodes = mesh.GetNumberOfNodes()
xyz = []
xyzIndices = {}
def AddNewPoint(pId0: MuscatIndex, pId1: MuscatIndex, phi0: MuscatFloat, phi1: MuscatFloat) -> MuscatIndex:
"""Function to add a new point to the mesh and recover the id
if the point exists, no double insertion is done.
Parameters
----------
pId0 : MuscatIndex
index of the first point on the mesh
pId1 : MuscatIndex
index of the second point on the mesh
phi0 : MuscatFloat
level set value at the first point
phi1 : MuscatFloat
level set value at the second point
Returns
-------
MuscatIndex
the id of the point
"""
opair = (pId0, pId1) if pId0 < pId1 else (pId1, pId0)
if opair in xyzIndices:
return xyzIndices[opair]
r = phi0/(phi0-phi1)
xpos = mesh.nodes[pId0, :]*(1-r) + mesh.nodes[pId1, :]*r
xyz.append(xpos)
newPointId = len(xyz) + nbNodes - 1
xyzIndices[opair] = newPointId
if computePointFieldTransferOperator:
ptf_i.extend([newPointId]*2)
ptf_j.extend([pId0, pId1])
ptf_values.extend([(1-r), r])
return newPointId
# we start with the elements that touch the isozero but remain in the positive or negative side
# we have to construct new elements (point, bar or triangle) when all but one node of the element are on the interface
for selection in ElementFilter()(interfaceMesh):
localConnectivity = selection.elements.connectivity
if selection.elementType == ED.Point_1:
continue
elif selection.elementType == ED.Bar_2:
pointsElements = interfaceMesh.elements.GetElementsOfType(ED.Point_1)
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
for n in range(2):
if phi[n] != 0: # n is the node that is not on the isozero
pointsElements.AddNewElement(conn[(n+1)%2], -1)
break
elif selection.elementType == ED.Triangle_3:
barsElements = interfaceMesh.elements.GetElementsOfType(ED.Bar_2)
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
for n in range(3):
if phi[n] != 0 and phi[(n+1)%3] == 0 and phi[(n+2)%3] == 0: # n is the node that is not on the isozero
if phi[n] > 0:
barsElements.AddNewElement([conn[(n+2)%3], conn[(n+1)%3]], -1)
else:
barsElements.AddNewElement([conn[(n+1)%3], conn[(n+2)%3]], -1)
break
elif selection.elementType == ED.Tetrahedron_4:
trisElements = interfaceMesh.elements.GetElementsOfType(ED.Triangle_3)
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
for n in range(4):
if phi[n] != 0 and phi[(n+1)%4] == 0 and phi[(n+2)%4] == 0 and phi[(n+3)%4] == 0: # n is the node that is not on the isozero
if (n%2 == 0 and phi[n] > 0) or (n%2 == 1 and phi[n] < 0):
trisElements.AddNewElement([conn[(n+3)%4], conn[(n+2)%4], conn[(n+1)%4]], -1)
else :
trisElements.AddNewElement([conn[(n+1)%4], conn[(n+2)%4], conn[(n+3)%4]], -1)
break
else:
raise RuntimeError(f"Cant Treat this type of elements {selection.elementType}, use FromQuadToLinear and/or MeshTetrahedrization functions in Containers first")
CleanDoubleElements(interfaceMesh)
pointsElements = interfaceMesh.elements.GetElementsOfType(ED.Point_1)
barsElements = interfaceMesh.elements.GetElementsOfType(ED.Bar_2)
trisElements = interfaceMesh.elements.GetElementsOfType(ED.Triangle_3)
quadsElements = interfaceMesh.elements.GetElementsOfType(ED.Quadrangle_4)
tetsElements = interfaceMesh.elements.GetElementsOfType(ED.Tetrahedron_4)
pyrsElements = interfaceMesh.elements.GetElementsOfType(ED.Pyramid_5)
wedgesElements = interfaceMesh.elements.GetElementsOfType(ED.Wedge_6)
# we treat the elements that cross the isozero, in this case we must cut them into new elements by adding points and also create an element (point, bar or triangle) on the interface
for selection in filterInterface1(mesh):
tags = selection.elements.tags
localConnectivity = selection.elements.connectivity
if selection.elementType == ED.Bar_2:
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
nodeId = AddNewPoint(conn[0], conn[1], phi[0], phi[1])
pointsElements.AddNewElement(nodeId, -1)
id2 = barsElements.AddNewElements([[conn[0], nodeId], [nodeId, conn[1]]], [eid, eid])-1
id1 = id2-1
for tag in tags:
if eid in tag.GetIds():
barsElements.tags.CreateTag(tag.name, False).AddToTag([id1, id2])
elif selection.elementType == ED.Triangle_3:
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
for n in range(3):
if phi[n] == 0: # here n is the node on the isozero
nodeId = AddNewPoint(conn[(n+1)%3], conn[(n+2)%3], phi[(n+1)%3], phi[(n+2)%3])
if phi[(n+1)%3] > 0:
barsElements.AddNewElement([conn[n], nodeId], -1)
else:
barsElements.AddNewElement([nodeId, conn[n]], -1)
id2 = trisElements.AddNewElements([[conn[n], nodeId, conn[(n+2)%3]], [conn[n], conn[(n+1)%3], nodeId]], [eid, eid])-1
id1 = id2-1
for tag in tags:
if eid in tag.GetIds():
trisElements.tags.CreateTag(tag.name, False).AddToTag([id1, id2])
break
elif phi[(n+1)%3]*phi[(n+2)%3] > 0: # here n is the node isolated whereas (n+1)%3 and (n+2)%3 are on the same side of the isozero
nodeId1 = AddNewPoint(conn[n], conn[(n+1)%3], phi[n], phi[(n+1)%3])
nodeId2 = AddNewPoint(conn[n], conn[(n+2)%3], phi[n], phi[(n+2)%3])
if phi[n] > 0: # the orientation depends on the sign of phi[n] to have a normal from negative to positive side
barsElements.AddNewElement([nodeId2, nodeId1], -1)
else:
barsElements.AddNewElement([nodeId1, nodeId2], -1)
id1 = quadsElements.AddNewElement([nodeId1, conn[(n+1)%3], conn[(n+2)%3], nodeId2], eid)-1
id2 = trisElements.AddNewElement([conn[n], nodeId1, nodeId2], eid)-1
for tag in tags:
if eid in tag.GetIds():
quadsElements.tags.CreateTag(tag.name, False).AddToTag(id1)
trisElements.tags.CreateTag(tag.name, False).AddToTag(id2)
break
elif selection.elementType == ED.Tetrahedron_4:
for eid in selection.indices:
conn = localConnectivity[eid, :]
phi = izoAtPoints[conn]
m = np.nonzero(phi)[0]
l = np.nonzero(phi==0)[0]
if l.size == 2: # the tet has two nodes on the isozero
nodeId = AddNewPoint(conn[m[0]], conn[m[1]], phi[m[0]], phi[m[1]])
if (m[0]%2 == 0 and phi[m[0]] > 0): # the orientation depends on the sign of phi[m[n]] but also on the rotation of the tet with respect to the reference tet
trisElements.AddNewElement([conn[l[1]], conn[l[0]], nodeId], -1)
else:
trisElements.AddNewElement([conn[l[0]], conn[l[1]], nodeId], -1)
conn1 = np.copy(conn)
conn1[m[0]] = nodeId
conn2 = np.copy(conn)
conn2[m[1]] = nodeId
id2 = tetsElements.AddNewElements([conn1, conn2], [eid, eid])-1
id1 = id2-1
for tag in tags:
if eid in tag.GetIds():
tetsElements.tags.CreateTag(tag.name, False).AddToTag([id1, id2])
elif l.size == 1: # the tet has one node on the isozero
for n in range(3):
if phi[m[(n+1)%3]]*phi[m[(n+2)%3]] > 0: # n is the node isolated in one side of the isozero
nodeId1 = AddNewPoint(conn[m[n]], conn[m[(n+1)%3]], phi[m[n]], phi[m[(n+1)%3]])
nodeId2 = AddNewPoint(conn[m[n]], conn[m[(n+2)%3]], phi[m[n]], phi[m[(n+2)%3]])
if (l[0]%2 == 0 and phi[m[n]] < 0) or (l[0]%2 == 1 and phi[m[n]] > 0):
trisElements.AddNewElement([conn[l[0]], nodeId2, nodeId1], -1)
else:
trisElements.AddNewElement([conn[l[0]], nodeId1, nodeId2], -1)
if l[0]%2 == 0:
id1 = tetsElements.AddNewElement([nodeId2, nodeId1, conn[m[n]], conn[l[0]]], eid)-1
id2 = pyrsElements.AddNewElement([nodeId2, conn[m[(n+2)%3]], conn[m[(n+1)%3]], nodeId1, conn[l[0]]], eid)-1
else:
id1 = tetsElements.AddNewElement([nodeId1, nodeId2, conn[m[n]], conn[l[0]]], eid)-1
id2 = pyrsElements.AddNewElement([nodeId1, conn[m[(n+1)%3]], conn[m[(n+2)%3]], nodeId2, conn[l[0]]], eid)-1
for tag in tags:
if eid in tag.GetIds():
tetsElements.tags.CreateTag(tag.name, False).AddToTag(id1)
pyrsElements.tags.CreateTag(tag.name, False).AddToTag(id2)
break
elif l.size == 0: # the tet has no node on the isozero
for n in range(4):
if phi[n]*phi[(n+1)%4] < 0 and phi[n]*phi[(n+2)%4] < 0 and phi[n]*phi[(n+3)%4] < 0: # n is the node isolated in one side of the isozero surface
nodeId1 = AddNewPoint(conn[n], conn[(n+1)%4], phi[n], phi[(n+1)%4])
nodeId2 = AddNewPoint(conn[n], conn[(n+2)%4], phi[n], phi[(n+2)%4])
nodeId3 = AddNewPoint(conn[n], conn[(n+3)%4], phi[n], phi[(n+3)%4])
if (n%2 == 0 and phi[n] > 0) or (n%2 == 1 and phi[n] < 0):
trisElements.AddNewElement([nodeId3, nodeId2, nodeId1], -1)
else:
trisElements.AddNewElement([nodeId1, nodeId2, nodeId3], -1)
if n%2 == 0:
id1 = tetsElements.AddNewElement([nodeId3, nodeId2, nodeId1, conn[n]], eid)-1
id2 = wedgesElements.AddNewElement([conn[(n+3)%4], conn[(n+2)%4], conn[(n+1)%4], nodeId3, nodeId2, nodeId1], eid)-1
else:
id1 = tetsElements.AddNewElement([nodeId1, nodeId2, nodeId3, conn[n]], eid)-1
id2 = wedgesElements.AddNewElement([conn[(n+1)%4], conn[(n+2)%4], conn[(n+3)%4], nodeId1, nodeId2, nodeId3], eid)-1
for tag in tags:
if eid in tag.GetIds():
tetsElements.tags.CreateTag(tag.name, False).AddToTag(id1)
wedgesElements.tags.CreateTag(tag.name, False).AddToTag(id2)
break
else:
for n in range(3):
if phi[3]*phi[n] > 0: # n is the node in the same side of the isozero as node 3
nodeId1 = AddNewPoint(conn[n], conn[(n+1)%3], phi[n], phi[(n+1)%3])
nodeId2 = AddNewPoint(conn[n], conn[(n+2)%3], phi[n], phi[(n+2)%3])
nodeId3 = AddNewPoint(conn[3], conn[(n+2)%3], phi[3], phi[(n+2)%3])
nodeId4 = AddNewPoint(conn[3], conn[(n+1)%3], phi[3], phi[(n+1)%3])
if phi[3] < 0:
quadsElements.AddNewElement([nodeId4, nodeId3, nodeId2, nodeId1], -1)
else:
quadsElements.AddNewElement([nodeId1, nodeId2, nodeId3, nodeId4], -1)
id2 = wedgesElements.AddNewElements([[nodeId1, nodeId2, conn[n], nodeId4, nodeId3, conn[3]], [nodeId1, nodeId4, conn[(n+1)%3], nodeId2, nodeId3, conn[(n+2)%3]]], [eid, eid])-1
id1 = id2-1
for tag in tags:
if eid in tag.GetIds():
wedgesElements.tags.CreateTag(tag.name, False).AddToTag([id1, id2])
break
else:
raise RuntimeError(f"Can t Treat this type of elements {selection.elementType}, use FromQuadToLinear and/or MeshTetrahedrization functions first")
if len(xyz) :
interfaceMesh.nodes = np.vstack((mesh.nodes, np.asarray(xyz)))
interfaceMesh.originalIDNodes = np.append(np.arange(nbNodes), -np.ones(len(xyz), dtype = MuscatIndex))
izoAtPoints = np.append(izoAtPoints, np.zeros(len(xyz), dtype = MuscatFloat))
outputMesh = Tetrahedrization(interfaceMesh)
from Muscat.Containers.Filters.FilterTools import GetFrozenFilter
frozenFilter = GetFrozenFilter(ElementFilter(), interfaceMesh)
# we contruct originalIds of the outputMesh with respect to the input mesh
for selection in ElementFilter()(outputMesh):
globalIds = selection.elements.originalIds
for eid in selection.indices:
idGlobal = globalIds[eid] # global number interfaceMesh of element eid
for selection2 in frozenFilter(interfaceMesh):
idLocal = idGlobal - selection2.meshOffset
if idLocal >= 0 and idLocal < selection2.elements.GetNumberOfElements():
selection.elements.originalIds[eid] = selection2.elements.originalIds[idLocal]
outputFilter = DifferenceFilter([compFilter, filter0])
outputMesh.MergeElements(ExtractElementsByElementFilter(mesh, outputFilter), force = True)
TestMeshValidity(outputMesh)
outputMesh.PrepareForOutput()
# create new tags
for selection in ElementFilter(nMask = (izoAtPoints > 0), nMaskTreatment = NM_AT_LEAST_ONE)(outputMesh):
selection.elements.tags.CreateTag("isoPositive", False).AddToTag(selection.indices)
for selection in ElementFilter(nMask = (izoAtPoints < 0), nMaskTreatment = NM_AT_LEAST_ONE)(outputMesh):
selection.elements.tags.CreateTag("isoNegative", False).AddToTag(selection.indices)
for selection in ElementFilter(nMask = (izoAtPoints == 0), nMaskTreatment = "ALL_NODES")(outputMesh):
selection.elements.tags.CreateTag("isoZero", False).AddToTag(selection.indices)
if computePointFieldTransferOperator:
from scipy import sparse
ptf_values = np.hstack((np.ones(nbNodes, dtype = MuscatFloat).ravel(), np.asarray(ptf_values, dtype = MuscatFloat).ravel()))
ptf_i = np.hstack((np.arange(nbNodes, dtype = MuscatIndex), np.asarray(ptf_i, dtype = MuscatIndex)))
ptf_j = np.hstack((np.arange(nbNodes, dtype = MuscatIndex), np.asarray(ptf_j, dtype = MuscatIndex)))
to = sparse.coo_matrix((ptf_values, (ptf_i, ptf_j)), shape = (outputMesh.GetNumberOfNodes(), nbNodes), copy = False)
return outputMesh, to
else:
return outputMesh, None
[docs]def CheckIntegrity(GUI: bool = False):
import Muscat.Containers.MeshCreationTools as UMCT
from Muscat.ImplicitGeometry.ImplicitGeometryObjects import ImplicitGeometrySphere
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
mesh2D = UMCT.CreateSquare(dimensions = [41, 41], origin = [-1, -1], spacing = [2/40, 2/40], ofTriangles = True)
def F1(x):
return np.linalg.norm(x, axis = 1) - 0.75
def F2(x):
return np.linalg.norm(x+[1,1], axis = 1) - 0.75
def F3(x):
return np.cos(2.1*np.pi*np.linalg.norm(x+[1,1], axis = 1))
splitMeshI, transfertOP = CutMeshIsoZero(mesh2D, F1, computePointFieldTransferOperator=True)
for elements in splitMeshI.elements:
elements.tags.RenameTag("isoPositive", "plastic", noError = True)
elements.tags.RenameTag("isoNegative", "fiber", noError = True)
elements.tags.RenameTag("isoZero", "interface", noError = True)
splitMeshI.nodeFields["isoI"] = transfertOP.dot(F1(mesh2D.nodes))
splitMeshI.nodeFields["isoII"] = transfertOP.dot(F2(mesh2D.nodes))
splitMeshI.nodeFields["isoIII"] = transfertOP.dot(F3(mesh2D.nodes))
splitMeshII, transfertOPII = CutMeshIsoZero(splitMeshI, splitMeshI.nodeFields["isoII"], computePointFieldTransferOperator = True, snapTolerance = 0)
for elements in splitMeshII.elements:
elements.tags.RenameTag("isoPositive", "plasticII", noError = True)
elements.tags.RenameTag("isoNegative", "fiberII", noError = True)
elements.tags.RenameTag("isoZero", "interfaceII", noError = True)
splitMesh = CutMeshIsoZero(splitMeshII, transfertOPII.dot(splitMeshI.nodeFields["isoIII"]), snapTolerance=1e-5)[0]
if GUI:
from Muscat.Actions.OpenInParaView import OpenInParaView
mesh2D.nodeFields["isoI"] = F1(mesh2D.nodes)
mesh2D.nodeFields["isoII"] = F2(mesh2D.nodes)
mesh2D.nodeFields["isoIII"] = F3(mesh2D.nodes)
OpenInParaView(mesh2D)
OpenInParaView(splitMesh)
print("Creating mesh")
n = 15
nn = n-1
myMesh = UMCT.CreateCube(dimensions = [n, n, n], origin = [0, 0, 0], spacing = [1./nn, 1./nn, 1./nn], ofTetras = True)
print("Creating phi")
sB0 = ImplicitGeometrySphere(radius = 0.5, center = [0, 0, 0])
phi = sB0(myMesh)
oo = IGToMesh(myMesh, phi)
omesh = oo.ComputeInterfaceMesh()
print(omesh)
ooII = IGToMesh(myMesh, abs(phi)+0.1)
omeshII = ooII.ComputeInterfaceMesh()
print(omeshII)
tempdir = TemporaryDirectory.GetTempPath()
WriteMeshToXdmf(tempdir + 'SphereIsoZero.xdmf', omesh)
print("Iso zero mesh in " + tempdir +'SphereIsoZero.xdmf')
phi2 = phi*2
NewPhi = Redistance(myMesh, phi2, method="Interp/Clamp")
WriteMeshToXdmf(tempdir + 'SpherePhiNew.xdmf', myMesh, PointFieldsNames = ["Phi", "PhiX2", "NewPhi"], PointFields = [phi, phi2, NewPhi])
print("phis on original mesh in " + tempdir + 'SpherePhiNew.xdmf')
error = max(np.abs(phi-NewPhi))
if error > 1e-2:
return "KO"
from Muscat.IO.GeofReader import ReadGeof
from Muscat.TestData import GetTestDataPath
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
tempdir = TemporaryDirectory.GetTempPath()
filename = GetTestDataPath()+"cube2.geof"
mesh = ReadGeof(filename)
print(mesh)
CleanDoubleElements(mesh)
mesh = Tetrahedrization(mesh)
print(mesh)
filename = tempdir + "inputMesh.xdmf"
PointFieldsNames = list(mesh.nodeFields.keys())
PointFields = list(mesh.nodeFields.values())
WriteMeshToXdmf(filename, mesh, PointFieldsNames = PointFieldsNames, PointFields = PointFields)
def F4(x):
return np.cos(1.3*np.pi*np.linalg.norm(x, axis = 1))
outputMesh, transfertOP = CutMeshIsoZero(mesh, F4, computePointFieldTransferOperator = True)
print(outputMesh)
filename = tempdir + "outputMesh.xdmf"
PointFieldsNames = list(outputMesh.nodeFields.keys())
PointFields = list(outputMesh.nodeFields.values())
WriteMeshToXdmf(filename, outputMesh, PointFieldsNames = PointFieldsNames, PointFields = PointFields)
print(tempdir)
return 'OK'
if __name__ == '__main__': # pragma: no cover
print(CheckIntegrity(GUI=False))