# -*- 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, UnionFilter
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.MeshFieldOperations import TransportPosToPoints
[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)
interfaceMesh = IGTM.ComputeInterfaceMesh()
pos = TransportPosToPoints(interfaceMesh, 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:
lcoon = selection.elements.connectivity[eid, :]
phi0 = self.phi[lcoon[0]]
phi1 = self.phi[lcoon[1]]
x = phi0/(phi0-phi1)
xpos = self.inputMesh.nodes[lcoon[0], :]*(1-x) + self.inputMesh.nodes[lcoon[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 = []
lcoon = selection.elements.connectivity[eid, :]
# if zero on the a point
neg_id = None
pos_id = None
for n in lcoon:
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 = lcoon[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, createPointElements = True ) -> Tuple[Mesh,np.ndarray ]:
"""Cut a mesh using the iso zero computed by the function isoZeroFunction.
This function works only for 2D meshes composed only by Point_1, Bar_2 and Triangle_3
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 is 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 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 np.sign(np.min(izoAtPoints)) == np.sign(np.max(izoAtPoints)) :
print(izoAtPoints)
raise ValueError("Iso zero not found")
if snapTolerance is None:
##here we code the span function
minbb, maxbb = mesh.ComputeBoundingBox()
snapTolerance = np.linalg.norm(maxbb-minbb)*1e-6
if snapTolerance > 0:
izoAtPoints = np.copy(izoAtPoints)
izoAtPoints[abs(izoAtPoints)< snapTolerance ] = 0.0
zeroElements = IntersectionFilter([
ElementFilter(nMask=izoAtPoints <= 0 , nMaskTreatment = NM_AT_LEAST_ONE),
ElementFilter(nMask=izoAtPoints >= 0 , nMaskTreatment = NM_AT_LEAST_ONE)
])
interfaceMesh = Mesh()
interfaceMesh.nodes = mesh.nodes
interfaceMesh.originalIDNodes = mesh.originalIDNodes
#list of point to be added to the mesh
xyz = []
xyzIndices = {}
pointElements = interfaceMesh.elements.GetElementsOfType(ED.Point_1)
pointIsoZeroTag = pointElements.tags.CreateTag("isoZero")
barsElements = interfaceMesh.elements.GetElementsOfType(ED.Bar_2)
barsIsoZeroTag = barsElements.tags.CreateTag("isoZero")
barsIsoPTag = barsElements.tags.CreateTag("isoPositive")
barsIsoNTag = barsElements.tags.CreateTag("isoNegative")
trisElements = interfaceMesh.elements.GetElementsOfType(ED.Triangle_3)
trisIsoZeroTag = trisElements.tags.CreateTag("isoZero")
trisIsoPTag = trisElements.tags.CreateTag("isoPositive")
trisIsoNTag = trisElements.tags.CreateTag("isoNegative")
# we use the negative part of the connectivity int to store the point origin
# (positive = original mesh, negative = added point)
def AddNewPoint(pId0:MuscatIndex, pId1:MuscatIndex, phi0:MuscatFloat, phi1:MuscatFloat) -> MuscatIndex:
"""Function to add a new point to the mesh and recover the id+1
if the point are exist, 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
th eid of the t
"""
opair = (pId0,pId1) if pId0 < pId1 else (pId1,pId0)
if opair in xyzIndices:
return xyzIndices[opair]
x= phi0/(phi0-phi1)
xpos = mesh.nodes[pId0, :]*(1-x) + mesh.nodes[pId1, :]*(x)
xyz.append(xpos)
newPointId = len(xyz)
xyzIndices[opair] = newPointId
if computePointFieldTransferOperator:
ptf_i.extend([newPointId-1]*2)
ptf_j.extend([pId0, pId1])
ptf_values.extend([(1-x), x])
return newPointId
def PropagateTags(tags, eid, telements, ids):
for tag in tags:
intag = np.intersect1d(tag.GetIds(),[eid])
if len(intag):
telements.tags.CreateTag(tag.name,False).AddToTag(ids)
for selection in zeroElements(mesh):
if ED.dimensionality[selection.elementType] == 0:
for eid in selection.indices:
lcoon = selection.elements.connectivity[eid, :]
if createPointElements:
pid = pointElements.AddNewElement(lcoon[0], 0) - 1
pointIsoZeroTag.AddToTag(pid)
elif ED.dimensionality[selection.elementType] == 1:
# intersection elements
for eid in selection.indices:
lcoon = selection.elements.connectivity[eid, :]
phi0 = izoAtPoints[lcoon[0]]
phi1 = izoAtPoints[lcoon[1]]
opair = (lcoon[0],lcoon[1]) if lcoon[0] < lcoon[1] else (lcoon[1],lcoon[0])
if phi0 == 0 and phi1 == 0:
if createPointElements:
pid = pointElements.AddNewElement([lcoon[0]], 0) - 1
pointIsoZeroTag.AddToTag(pid)
pid = pointElements.AddNewElement([lcoon[1]], 0) - 1
pointIsoZeroTag.AddToTag(pid)
bar_I_Id = barsElements.AddNewElement(lcoon, 0 ) - 1
PropagateTags(selection.elements.tags, eid, barsElements, [bar_I_Id] )
barsIsoZeroTag.AddToTag(bar_I_Id)
elif phi0 == 0:
if createPointElements:
pid = pointElements.AddNewElement([lcoon[0]], 0) - 1
pointIsoZeroTag.AddToTag(pid)
bar_I_Id = barsElements.AddNewElement(lcoon, 0 ) - 1
PropagateTags(selection.elements.tags, eid, barsElements, [bar_I_Id] )
if phi1 > 0:
barsIsoPTag.AddToTag(bar_I_Id)
else:
barsIsoNTag.AddToTag(bar_I_Id)
elif phi1 == 0:
if createPointElements:
pid = pointElements.AddNewElement(lcoon[1], 0) - 1
pointIsoZeroTag.AddToTag(pid)
bar_I_Id = barsElements.AddNewElement(lcoon, 0 ) - 1
PropagateTags(selection.elements.tags, eid, barsElements, [bar_I_Id] )
if phi0 > 0:
barsIsoPTag.AddToTag(bar_I_Id)
else:
barsIsoNTag.AddToTag(bar_I_Id)
else:
nodeId = AddNewPoint(lcoon[0], lcoon[1], phi0, phi1)
if createPointElements:
pid = pointElements.AddNewElement(-1*(nodeId), 0) - 1
pointIsoZeroTag.AddToTag(pid)
bar_I_Id = barsElements.AddNewElement([lcoon[0], -nodeId], 0 ) - 1
bar_II_Id = barsElements.AddNewElement([ -nodeId, lcoon[1]], 0 ) - 1
PropagateTags(selection.elements.tags, eid, barsElements, [bar_I_Id, bar_II_Id] )
if phi0 < 0:
barsIsoNTag.AddToTag(bar_I_Id)
barsIsoPTag.AddToTag(bar_II_Id)
else:
barsIsoPTag.AddToTag(bar_I_Id)
barsIsoNTag.AddToTag(bar_II_Id)
pointElements.Tighten()
barsElements.Tighten()
elif ED.dimensionality[selection.elementType] == 2:
if selection.elementType == ED.Triangle_3:
#
# 2 2 (-)
# /| /|
# / | / |
# / | / |
# -3 / | / |
# ./________| -2 / _| -2
# / \ / | / _-¨/ |
# / \ / | / _-¨ ¨ / |
# /_______\/_____| /_-¨_____/_____|
# 0 -1 1 0 -1 1
# (-) (+)
# example for (-1,1,-1)
database = {}
database[(-1, 1, -1)] = ((1,2), # nodes to be created (positives )
( np.array([0,-1,-2]),np.array([0,-2,2]),), # tris of the negative part (negative numbers are given with respect the first tuple)
(np.array([-1,-2]),) , # bar of the iso zero
(np.array([-1,1,-2]),) # tris of the positive part
)
database[(-1, 1, 1)] = ((1,3), (np.array([ 0,-1,-2]), ), (np.array([-1,-2]),), (np.array([-1, 1,-2]),np.array([1, 2,-2]),))
database[( 1, 1,-1)] = ((2,3), (np.array([ 2,-2,-1]), ), (np.array([-2,-1]),), (np.array([ 0, 1,-2]),np.array([1,-1,-2]),))
database[( 1,-1, 1)] = ((1,2), (np.array([ 1,-2,-1]), ), (np.array([-2,-1]),), (np.array([ 0,-1,-2]),np.array([0,-2, 2]),))
database[(-1,-1, 1)] = ((2,3), (np.array([ 0, 1,-2]), np.array([1, -1,-2])), (np.array([-1,-2]),), (np.array([-1, 2,-2]), ))
database[( 1,-1,-1)] = ((1,3), (np.array([ 1,-2,-1]), np.array([1, 2, -2])), (np.array([-2,-1]),), (np.array([ 0,-1,-2]), ))
database[(-1, 0,+1)] = ((3,) ,(np.array([ 0, 1,-1]),), (np.array([ 1,-1]),), (np.array([-1, 1, 2]),))
database[(-1,+1, 0)] = ((1,) ,(np.array([ 0,-1, 2]),), (np.array([-1, 2]),), (np.array([ 1, 2,-1]),))
database[( 0,-1,+1)] = ((2,) ,(np.array([ 0, 1,-1]),), (np.array([-1, 0]),), (np.array([ 0,-1, 2]),))
database[( 0,+1,-1)] = ((2,) ,(np.array([ 0,-1, 2]),), (np.array([ 0,-1]),), (np.array([ 0, 1,-1]),))
database[(+1,-1, 0)] = ((1,) ,(np.array([ 1, 2,-1]),), (np.array([ 2,-1]),), (np.array([ 0,-1, 2]),))
database[(+1, 0,-1)] = ((3,) ,(np.array([-1, 1, 2]),), (np.array([-1, 1]),), (np.array([ 0, 1,-1]),))
database[( 0, 1, 1)] = ((), (), (), (np.array([0, 1, 2]),))
database[( 1, 0, 1)] = ((), (), (), (np.array([0, 1, 2]),))
database[( 1, 1, 0)] = ((), (), (), (np.array([0, 1, 2]),))
database[( 0, -1, -1)] = ((), (np.array([0, 1, 2]),), (), ())
database[( -1, 0, -1)] = ((), (np.array([0, 1, 2]),), (), ())
database[( -1,-1, 0)] = ((), (np.array([0, 1, 2]),), (), ())
database[( 1, 0, 0)] = ((), (), (np.array([2,1,0]),), (np.array([0,1, 2]),) )
database[( 0, 1, 0)] = ((), (), (np.array([2,0,0]),), (np.array([0,1, 2]),) )
database[( 0, 0, 1)] = ((), (), (np.array([1,0,0]),), (np.array([0,1, 2]),) )
database[(-1, 0, 0)] = ((), (np.array([0,1, 2]),), (np.array([1,2,0]),), ())
database[( 0,-1, 0)] = ((), (np.array([0,1, 2]),), (np.array([0,2,0]),), ())
database[( 0, 0,-1)] = ((), (np.array([0,1, 2]),), (np.array([0,1,0]),), ())
edges = [v[1] for v in ED.faces1[selection.elementType]]
for eid in selection.indices:
lcoon = selection.elements.connectivity[eid, :]
lphi = izoAtPoints[lcoon]
signs = tuple(np.sign(lphi))
if signs not in database:
print(signs)
continue
newPoints, negativeElems, isoZeroBar, positiveElems = database.get(signs)
nodeId = np.empty(2,dtype=int)
for cccc, i in enumerate(newPoints):
I0, I1 = edges[i-1]
nodeId[cccc] = AddNewPoint(lcoon[I0], lcoon[I1],lphi[I0], lphi[I1] )
newElems = []
for barConn in isoZeroBar:
gconn = lcoon[abs(barConn)*(barConn>=0)]*(barConn>=0) - nodeId[abs(barConn)*(barConn<0)-1]*(barConn<0)
newElems.append( barsElements.AddNewElement(gconn[:2], 0 ) - 1 )
barsIsoZeroTag.AddToTag(newElems)
#PropagateTags(selection.elements.tags, eid, barsElements, newElems )
newElems = []
for ne in negativeElems:
gconn = lcoon[abs(ne)*(ne>=0)]*(ne>=0) - nodeId[abs(ne)*(ne<0)-1]*(ne<0)
tri_I_ID = trisElements.AddNewElement(gconn, -100 ) - 1
trisIsoNTag.AddToTag(tri_I_ID)
newElems.append(tri_I_ID)
for pe in positiveElems:
gconn = lcoon[abs(pe)*(pe>=0)]*(pe>=0) - nodeId[abs(pe)*(pe<0)-1]*(pe<0)
tri_I_ID = trisElements.AddNewElement(gconn, -200) - 1
trisIsoPTag.AddToTag(tri_I_ID)
newElems.append(tri_I_ID)
PropagateTags(selection.elements.tags, eid, trisElements, newElems )
trisElements.Tighten()
barsElements.Tighten()
else:
raise RuntimeError(f"Cant Treat this type of elements {selection.elementType}")
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.Containers.Filters.FilterTools import GetComplementaryFilter
comFilter = GetComplementaryFilter(zeroElements)
output = ExtractElementsByElementFilter(mesh,comFilter)
for selection in IntersectionFilter( [comFilter,ElementFilter(nMask=izoAtPoints <= 0 , nMaskTreatment = NM_AT_LEAST_ONE)] )(output):
selection.elements.tags.CreateTag("isoNegative",False).AddToTag(selection.indices)
for selection in IntersectionFilter( [comFilter,ElementFilter(nMask=izoAtPoints >= 0 , nMaskTreatment = NM_AT_LEAST_ONE)] )(output):
selection.elements.tags.CreateTag("isoPositive",False).AddToTag(selection.indices)
for elems in interfaceMesh.elements:
elems.connectivity[elems.connectivity<0] = abs(elems.connectivity[elems.connectivity<0]) + output.GetNumberOfNodes()-1
if len(xyz) :
output.nodes = np.vstack((output.nodes, np.asarray(xyz)))
output.originalIDNodes = np.hstack((output.originalIDNodes, -(np.arange(len(xyz))+1)))
output.MergeElements(interfaceMesh, force=True)
#propagate tags from the old mesh to the new isoZero
for name in output.GetNamesOfElementTags():
if name in ["isoNegative","isoPositive", "isoZero"]:
continue
for selection in ElementFilter(eTag="isoZero",nFilter=NodeFilter(eTag=name))(output):
selection.elements.tags.CreateTag(name,False).AddToTag(selection.indices)
if computePointFieldTransferOperator:
from scipy import sparse
ptf_values = np.hstack( (np.ones(mesh.nodes.shape[0],dtype=MuscatFloat).ravel(), np.asarray(ptf_values,dtype=MuscatFloat).ravel() ) )
ptf_i = np.hstack( (np.arange(mesh.nodes.shape[0],dtype=MuscatIndex), np.asarray(ptf_i,dtype=MuscatIndex)+mesh.nodes.shape[0]) )
ptf_j = np.hstack( (np.arange(mesh.nodes.shape[0],dtype=MuscatIndex), np.asarray(ptf_j,dtype=MuscatIndex)) )
to = sparse.coo_matrix((ptf_values, (ptf_i, ptf_j)), shape=(output.GetNumberOfNodes(), mesh.GetNumberOfNodes()), copy=False)
return output, to
else:
return output, 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")
# print("NewPhi")
# print(NewPhi)
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"
return 'OK'
if __name__ == '__main__': # pragma: no cover
print(CheckIntegrity(GUI=True))