# -*- 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 numpy as np
import math
import BasicTools.Helpers.ParserHelper as PH
from BasicTools.ImplicitGeometry.ImplicitGeometryFactory import RegisterClass
from BasicTools.ImplicitGeometry.ImplicitGeometryBase import ImplicitGeometryBase,dsin,dcos
from BasicTools.ImplicitGeometry.ImplicitGeometryOperators import ImplicitGeometryIntersection,ImplicitGeometryUnion
from BasicTools.Containers.Filters import ElementFilter
[docs]class ImplicitGeometryWrapped(ImplicitGeometryBase):
"""Wrapper to add the ImplicitGeometry API arround a numpy array
field : numpy array
"""
def __init__(self, field=None ):
super(ImplicitGeometryWrapped,self).__init__()
self.field = field
[docs] def ApplyVector(self, support,cellCenter=False):
if type(support).__module__ == np.__name__:
if support.shape[0] != len(self.field):
raise(Exception("suport and field not compatible"))
else:
if cellCenter :
if support.GetNumberOfElements() != len(self.field):
raise(Exception("suport and field not compatible"))
elif support.GetNumberOfNodes() != len(self.field):
raise(Exception("suport and field not compatible"))
return self.ApplyInsideOut(self.field)
[docs] def GetDistanceToPoint(self,_pos):
raise Exception("Cant use this function") # pragma: no cover
RegisterClass("Wrapped",ImplicitGeometryWrapped)
[docs]def CreateImplicitGeometryExternalSurface(ops):
"""
ImplicitGeometry of the external surface of a mesh (support)
support : Grid to compute the surface
<All id="x" [support=""| ls=""] >
"""
res = ImplicitGeometryExternalSurface()
if ("support" in ops or "ls" in ops ) :
if "ls" in ops:
sup = ops["ls"].support
else:
sup = ops["support"]
res.SetSupport(sup)
else:
raise(Exception('Need a (ls or support) '))
return res
[docs]class ImplicitGeometryExternalSurface(ImplicitGeometryBase):
"""ImplicitGeometry of the external surface of a mesh (support)
support : mesh to compute the surface
"""
def __init__(self,support = None):
super(ImplicitGeometryExternalSurface,self).__init__()
self.offset = 1E-3
if support is not None:
self.SetSupport(support)
[docs] def SetSupport(self, support):
from BasicTools.Containers.NativeTransfer import NativeTransfer
self.__nativeTransfer = NativeTransfer()
self.__nativeTransfer.SetVerbose(False)
from BasicTools.FE.Fields.FEField import FEField
from BasicTools.FE.FETools import PrepareFEComputation
space, numberings, offset, NGauss = PrepareFEComputation(support,numberOfComponents=1)
self.__field = FEField("", mesh=support, space=space, numbering=numberings[0])
self.__nativeTransfer.SetSourceFEField(self.__field, ElementFilter(support, dimensionality=support.GetElementsDimensionality() ) )
self.__nativeTransfer.SetTransferMethod("Interp/Clamp")
from BasicTools.Containers.UnstructuredMeshModificationTools import ComputeSkin, CleanLonelyNodes
skinmesh = ComputeSkin(support,inPlace=False)
CleanLonelyNodes(skinmesh,inPlace=True)
space, numberings, offset, NGauss = PrepareFEComputation(skinmesh,numberOfComponents=1)
self.__Sfield = FEField("", mesh=skinmesh, space=space, numbering=numberings[0])
self.__SnativeTransfer = NativeTransfer()
self.__SnativeTransfer.SetVerbose(False)
self.__SnativeTransfer.SetSourceFEField(self.__Sfield )
self.__SnativeTransfer.SetTransferMethod("Interp/Clamp")
[docs] def GetDistanceToPoint(self,pos):
self.__nativeTransfer.SetTargetPoints(pos)
self.__nativeTransfer.Compute()
# status must be 1 inside or 3 Clamp not other
status = self.__nativeTransfer.GetStatus()
self.__SnativeTransfer.SetTargetPoints(pos)
self.__SnativeTransfer.Compute()
op = self.__SnativeTransfer.GetOperator()
admissible = {1, 3}
if not admissible.issuperset(np.unique(status)):
raise RuntimeError("Error in field transfer.")
# the expression status -2 generate a -1 for the inside and 1 for the outside
vectorDistance = op.dot(self.__Sfield.mesh.nodes)- pos
distance = np.sqrt(np.sum((vectorDistance )**2,axis=1))
return distance*(status.flatten()-2) + self.offset
RegisterClass("All",ImplicitGeometryExternalSurface,CreateImplicitGeometryExternalSurface)
####################### objects ################################
[docs]def CreateImplicitGeometryByETag(ops):
res = ImplicitGeometryByETag()
if "eTags" not in ops:
raise(Exception('Need a eTags'))
if ("support" in ops or "ls" in ops ) :
if "ls" in ops:
sup = ops["ls"].support
else:
sup = ops["support"]
else:
raise(Exception('Need a (ls or support) '))
res.offset = (PH.ReadFloat(ops.get("offset",res.offset)))
res.eps = (PH.ReadFloat(ops.get("eps",res.eps)))
if "dim" in ops:
res.dim = PH.ReadInt(ops["dim"])
res.SetSupportAndZones(sup,PH.ReadStrings(ops["eTags"] ))
return res
[docs]class ImplicitGeometryByETag(ImplicitGeometryBase):
"""ImplicitGeometry based in a eTags from a mesh (support) or levelset
support : mesh from where to extract the eTags
ls : levelset from where to extract the support to extract the eTags
etags : list of etags
offset : offset from the eTag (negative will grow the zone )
Note: For the moment works only for volume and surface eTag (no lines)
the implementation compute the distance on the the skin of the etag
"""
def __init__(self ):
super(ImplicitGeometryByETag,self).__init__()
self.op = None
# we add an epsiloin to treat cases where all 4 nodes on one tetra are
# in the iso (case of a tetra in a corrner/edge of a zone)
self.eps = -1.e-4
self.offset = 0.
[docs] def SetSupportAndZones(self,support, etags):
from BasicTools.Containers.UnstructuredMeshInspectionTools import ExtractElementsByElementFilter
from BasicTools.Containers.Filters import ElementFilter
sup = ExtractElementsByElementFilter(support, elementFilter=ElementFilter(mesh= support, tags= etags) )
if sup.GetElementsDimensionality() == support.GetElementsDimensionality():
self.op = ImplicitGeometryExternalSurface()
self.op.offset = 0
self.op.SetSupport(sup)
else:
self.op = ImplicitGeometryStl()
self.op.SetSurface(sup)
def __str__(self):
res = "ImplicitGeometryByETag\n"
#res += " support: " + str(self.support)+"\n"
return res
[docs] def ApplyVector(self, support,cellCenter=False):
vals = self.op.ApplyVector(support,cellCenter=cellCenter)
return self.ApplyInsideOut(vals+(self.eps+self.offset))
[docs] def GetDistanceToPoint(self,_pos):
vals = self.op.GetDistanceToPoint(_pos)
return self.ApplyInsideOut(vals+(self.eps+self.offset))
RegisterClass("ETag",ImplicitGeometryByETag,CreateImplicitGeometryByETag)
[docs]def CreateImplicitGeometryByETagII(ops):
res = ImplicitGeometryByETagII()
if "eTags" not in ops:
raise(Exception('Need a eTags'))
if ("support" in ops or "ls" in ops ) :
if "ls" in ops:
sup = ops["ls"].support
else:
sup = ops["support"]
if "dim" in ops:
dim = PH.ReadInt(ops["dim"])
else:
dim = None
res.SetSupportAndZones(sup,PH.ReadStrings(ops["eTags"] ), dim = dim)
else:
raise(Exception('Need a (ls or support) '))
res.offset = (PH.ReadFloat(ops.get("offset",res.offset)))
return res
[docs]class ImplicitGeometryByETagII(ImplicitGeometryBase):
"""ImplicitGeometry based in a eTags from a mesh (support) or levelset
support : mesh from where to extract the eTags
ls : levelset from where to extract the support to extract the eTags
etags : list of etags
offset : offset from the eTag (negative will grow the zone )
Note: for surface and line the resulting fields is possitive (no inside)
the implementation puts zeros in the interface beetwen the elements in the
etag and the other elements (1 outside, 0 in the interface, -1 inside)
Some special treatement is done for elements with all the nodes with zero values
(presence of values -0.5 and 0.5 in the resulting field)
"""
def __init__(self ):
super(ImplicitGeometryByETagII,self).__init__()
self.offset = 0.
from BasicTools.Containers.Filters import ElementFilter
self.elementFilter = ElementFilter()
self.offset = 0.
self.numberOfElementPseudoDistance = 10
[docs] def SetSupportAndZones(self,support, etags, dim=None):
self.elementFilter.mesh = support
self.elementFilter.tags = etags
self.elementFilter.dimensionality = dim
def __str__(self):
res = "ImplicitGeometryByETagII\n"
res += " " + str(self.elementFilter)
return res
[docs] def ApplyVector(self, support, cellCenter=False):
if cellCenter:
raise (Exception("Not implemented"))
from BasicTools.Containers.Filters import ElementFilterToImplicitField
self.elementFilter.mesh = support
vals = ElementFilterToImplicitField(self.elementFilter,self.numberOfElementPseudoDistance)
res = self.ApplyInsideOut(vals+(self.offset))
return res
[docs] def GetDistanceToPoint(self,_pos):
if _pos is self.elementFilter.mesh.nodes:
from BasicTools.Containers.Filters import ElementFilterToImplicitField
vals = ElementFilterToImplicitField(self.elementFilter,self.numberOfElementPseudoDistance)
return self.ApplyInsideOut(vals+(self.offset))
raise
RegisterClass("ETagII",ImplicitGeometryByETagII,CreateImplicitGeometryByETagII)
[docs]class ImplicitGeometryAxisAlignBox(ImplicitGeometryBase):
"""ImplicitGeometry based ona Axis Alinged Box
origin : origin in the mox (Xmin Ymin Zmin)
size : size of the box (Xl, Yl Zl)
"""
def __init__(self, origin=None, size=None ):
super(ImplicitGeometryAxisAlignBox,self).__init__()
if origin is None:
self.origin = np.array([0.,0.,0.], dtype= float)
else:
self.origin = np.array(origin, dtype= float)
if size is None:
self.size = np.array([1.,1.,1.], dtype= float)
else:
if type(size) is float:
self.size = np.array([1.,1.,1.], dtype= float)*size
else:
self.size = size
def __str__(self):
res = "ImplicitGeometryAxisAlignBox \n"
res += " origin: " + str(self.origin)+"\n"
res += " size: " + str(self.size)+"\n"
return res
[docs] def GetBoundingMin(self):
return self.origin
[docs] def GetBoundingMax(self):
return self.origin + self.size
[docs] def SetSupport(self,support):
support.ComputeBoundingBox()
self.origin = support.boundingMin
self.size = support.boundingMax - support.boundingMin
[docs] def GetDistanceToPoint(self,_pos):
walls = []
data = [[-1, 0, 0],
[ 0,-1, 0],
[ 0, 0,-1]]
for normal in data:
Obj = ImplicitGeometryPlane()
Obj.point=self.GetBoundingMin()
Obj.normal=np.array(normal,dtype=float)
walls.append(Obj)
data = [[ 1, 0, 0],
[ 0, 1, 0],
[ 0, 0, 1]]
for normal in data:
Obj = ImplicitGeometryPlane()
Obj.point=self.GetBoundingMax()
Obj.normal=np.array(normal,dtype=float)
walls.append(Obj)
return self.ApplyInsideOut(ImplicitGeometryIntersection(walls).ApplyVector(_pos))
RegisterClass("AABox",ImplicitGeometryAxisAlignBox)
[docs]def CreateImplicitGeometrySphereFromNTag(ops):
res = ImplicitGeometrySphereFromNTag()
res.SetRadius(PH.ReadFloat(ops.get("radius",1.)))
res.nTag = PH.ReadString(ops["nTag"])
res.SetSupport(ops["ls"].support)
return res
[docs]class ImplicitGeometrySphereFromNTag(ImplicitGeometryBase):
"""ImplicitGeometry based on a ntags
ls : levelset to extrat ntag
nTag : name of the nTag to retrieve the centers
radius : radius of the sphere (r)
"""
def __init__(self,radius= 1. ):
super(ImplicitGeometrySphereFromNTag,self).__init__()
self.SetRadius(radius)
self.centers = None
self.nTag = None
[docs] def SetLs(self,ls):
self.SetSupport(ls.support)
[docs] def SetSupport(self,support):
self.centers = support.nodes[support.nodesTags[self.nTag].GetIds(),:]
[docs] def SetRadius(self,val):
self.radius = val
def __str__(self):
res = "ImplicitGeometrySphereFromNTag \n"
res += " radius: " + str(self.radius) + "\n"
res += " nTags:" + str(self.nTag) + "\n"
return res
[docs] def GetDistanceToPoint(self,_pos):
if len(_pos.shape) == 1:
pos = np.array([_pos])
else:
pos = _pos
res = np.sqrt(np.sum((pos-self.centers[0,:])**2,axis=1))-self.radius
for i in range(1,self.centers.shape[0]):
res2 = np.sqrt(np.sum((pos-self.centers[i,:])**2,axis=1))-self.radius
res = np.minimum(res,res2)
return self.ApplyInsideOut(res)
RegisterClass("SphereFromNTag",ImplicitGeometrySphereFromNTag,CreateImplicitGeometrySphereFromNTag )
[docs]class ImplicitGeometrySphere(ImplicitGeometryBase):
"""ImplicitGeometry based on a sphere
center : center of the sphere (X Y Z)
radius : radius of the sphere (r)
"""
def __init__(self,radius= 1.,center= None ):
super(ImplicitGeometrySphere,self).__init__()
if center is not None:
self.SetCenter(PH.ReadFloats(center))
else:
self.SetCenter(PH.ReadFloats([0., 0., 0.]) )
self.SetRadius(PH.ReadFloat(radius))
def __str__(self):
res = "ImplicitGeometrySphere \n"
res += " center: " + str(self.center)+"\n"
res += " radius: "+str(self.radius)
return res
[docs] def SetCenter(self,center):
self.center = PH.ReadFloats(center)
[docs] def SetRadius(self,val):
self.radius = PH.ReadFloat(val)
[docs] def GetDistanceToPoint(self,_pos):
if len(_pos.shape) == 1:
pos = np.array([_pos])
else:
pos = _pos
res = np.sqrt(np.sum((pos-self.center)**2,axis=1))-self.radius
return self.ApplyInsideOut(res)
[docs] def GetBoundingMin(self):
res = self.center - self.radius
return res
[docs] def GetBoundingMax(self):
res = self.center + self.radius
return res
RegisterClass("Sphere",ImplicitGeometrySphere )
[docs]class ImplicitGeometryCylinder(ImplicitGeometryBase):
"""ImplicitGeometry based on a cylinder
center1 : first point on the axis (X Y Z)
center2 : second point on the axis (X Y Z)
radius : radius of the cylinder (r)
wcups : if we compute the distance to the cups of the cylinder
"""
def __init__(self,center1=None,center2=None,radius=1.,wcups=True):
super(ImplicitGeometryCylinder,self).__init__()
if center1 is None:
self.center1 = np.zeros((3,))
else:
self.center1 = np.array(center1,copy=True,dtype=float)
if center2 is None:
self.center2 = np.array([1.,0.,0.])
else:
self.center2 = np.array(center2,copy=True,dtype=float)
self.radius = radius
self.wcups = wcups
#self.infinit = infinit
def __str__(self):
res = "ImplicitGeometryCylinder \n"
res += " center1: " + str(self.center1)+"\n"
res += " center2: " + str(self.center2)+"\n"
res += " radius: "+str(self.radius)
res += " width cups: "+str(self.wcups)
#res += " infinit: "+str(self.infinit)
return res
[docs] def GetDistanceToPoint(self,_pos):
if len(_pos.shape) == 1:
pos = np.array([_pos])
else:
pos = _pos
u = (self.center2 -self.center1).astype(float) # director vector
nu = np.linalg.norm(u)
u /= nu
# distance to the center of the cylinder
d = np.linalg.norm(np.cross(pos-self.center1,u),axis=1) - self.radius
if self.wcups:
# distance to the planes (the cups)
pA1 = ImplicitGeometryPlane(point=self.center1, normal=-u)
pA2 = ImplicitGeometryPlane(point=self.center2, normal=u)
pA = ImplicitGeometryIntersection([pA1,pA2]).GetDistanceToPoint(pos)
res = np.maximum(d,pA)
else:
res = d
return self.ApplyInsideOut(res)
[docs] def GetBoundingMin(self):
res = np.minimum(self.center1, self.center2)
diff = self.center2 -self.center1
l = math.sqrt(np.linalg.norm(diff))
if diff[0] > 0:
res[0] = res[0]-l*self.radius/abs(diff[0])
if diff[1] > 0:
res[1] = res[1]-l*self.radius/abs(diff[1])
if diff[2] > 0:
res[2] = res[2]-l*self.radius/abs(diff[2])
return res
[docs] def GetBoundingMax(self):
res = np.maximum(self.center1, self.center2)
diff = self.center2 - self.center1
l = math.sqrt(np.linalg.norm(diff))
if diff[0] > 0:
res[0] = res[0]+l*self.radius/abs(diff[0])
if diff[1] > 0:
res[1] = res[1]+l*self.radius/abs(diff[1])
if diff[2] > 0:
res[2] = res[2]+l*self.radius/abs(diff[2])
return res
RegisterClass("Cylinder",ImplicitGeometryCylinder )
[docs]def CreateImplicitGeometryPlane(ops):
obj = ImplicitGeometryPlane()
PH.ReadProperties(ops,["point","offset"], obj)
obj.SetNormal(PH.ReadFloats(ops["normal"]))
# _setProps = ["normal":np.array([1.,0,0])]
return obj
[docs]class ImplicitGeometryPlane(ImplicitGeometryBase):
"""ImplicitGeometry based on a plane
point : point on the plane (X Y Z)
normal : normale of the plane (X Y Z)
offset : offset to the plane (r)
"""
def __init__(self, point=None, normal=None, offset = 0.0):
super(ImplicitGeometryPlane,self).__init__()
if point is None:
self.point = np.array([0.,0.,0.],dtype=float,copy=True)
else:
self.point = np.array(point,copy=True)
if normal is None:
self.SetNormal(np.array([1,0,0], dtype =float ) )
else :
self.SetNormal(normal)
self.offset = float(offset)
[docs] def SetNormal(self,invec):
self.normal = np.array(invec,copy=True)
self.normal = self.normal/np.linalg.norm(self.normal)
#normal = property(fset=SetNormal)
def __str__(self):
res = "ImplicitGeometryPlane\n"
res += " point: " + str(self.point) + "\n"
res += " normal: " + str(self.normal) + "\n"
res += " offset: " + str(self.offset) + "\n"
return res
[docs] def GetDistanceToPoint(self,pos):
d = self.normal.dot(self.point)
res = np.sum(pos*self.normal,axis=1) - d - self.offset
#res = np.sum(self.normal*(pos-self.point),axis=1)
return self.ApplyInsideOut(res)
RegisterClass("Plane",ImplicitGeometryPlane,CreateImplicitGeometryPlane)
[docs]class ImplicitGeometryGyroid(ImplicitGeometryBase):
"""ImplicitGeometry based on a gyroid
scale : scale of the gyroid (1. defaul)
offset : offset of the gyroid (X Y Z)
type : [0-1] version of gyroid (0 default)
wall : [True-False] wall or body (false default)
wallThickness : parameter to control the wall thickness
"""
def __init__(self, scale=1., offset=[0.,0.,0.] ):
super(ImplicitGeometryGyroid,self).__init__()
self.scale = scale
self.offset = np.array(offset,dtype=float)
self.type = 0
self.wall = False
self.wallThickness = 0.5
def __str__(self):
res = "ImplicitGeometryGyroid \n"
res += " scale: " + str(self.scale)+"\n"
res += " offset: "+str(self.offset)+"\n"
return res
[docs] def GetDistanceToPoint(self,pos):
npos = pos/self.scale +self.offset
x = npos[:,0]*np.pi
y = npos[:,1]*np.pi
z = npos[:,2]*np.pi
if self.type == 0:
res = (dsin(x)*dcos(y) + dsin(y)*dcos(z) + dsin(z)*dcos(x))
res *= self.scale/(1.55*np.pi)
elif self.type == 1:
res = np.sin(x)*np.cos(y) + np.sin(y)*np.cos(z) + np.sin(z)*np.cos(x)
#sin = np.sin
#cos = np.cos
#ng = 0.01+np.sqrt(abs(sin(x)*sin(y) - cos(y)*cos(z))**2 + abs(sin(x)*sin(z) - cos(x)*cos(y))**2 + abs(sin(y)*sin(z) - cos(x)*cos(z))**2)
#res /= ng
# aproximation of the signed distance dist = f/norm(grad(f))
else:
res = 0.5*(+np.sin(2.*x)*np.cos(y)*np.sin(z)
+np.sin(2.*y)*np.cos(z)*np.sin(x)
+np.sin(2.*z)*np.cos(x)*np.sin(y)
)-0.5*(
+np.cos(2.*x)*np.cos(2.*y)
+np.cos(2.*y)*np.cos(2.*z)
+np.cos(2.*z)*np.cos(2.*x)
)+0.15
if self.wall:
res = np.abs(res)-self.wallThickness/2.
return self.ApplyInsideOut(res)
RegisterClass("Gyroid",ImplicitGeometryGyroid)
[docs]class ImplicitGeometry60D(ImplicitGeometryBase):
"""ImplicitGeometry using walls at 60 degrees in the (x,y) plane
lx : scale
w : wall thickness
"""
def __init__(self, lx=1., w=0.1 ):
super(ImplicitGeometry60D,self).__init__()
self.lx = lx
self.w = w
def __str__(self):
res = "ImplicitGeometry60D \n"
res += " lx: " + str(self.lx)+"\n"
res += " width: "+str(self.w)+"\n"
return res
[docs] def GetDistanceToPoint(self,pos):
v = math.sqrt(3.)
llx = self.lx/2.
lw =self.w/2.
aw = lw/llx
mpos = abs((np.mod(pos/llx,[2.,2*v,1.])-[1,v,0]) )
pA1 = ImplicitGeometryPlane(point=np.array([0.,0,0.]), normal=np.array([0,1.,0,]), offset = aw )
normal = np.array([-v,1.,0])
pB1 =ImplicitGeometryPlane(point=[0.,0,0], normal=-normal, offset= aw)
pB2 =ImplicitGeometryPlane(point=[0.,0,0], normal= normal, offset= aw)
pB = ImplicitGeometryIntersection([pB1,pB2])
pC1 =ImplicitGeometryPlane(point=[0,v,0], normal=np.array([0,-1,0,]), offset = aw)
res = ImplicitGeometryUnion([pA1,pB,pC1]).GetDistanceToPoint(mpos)
return self.ApplyInsideOut(res)
RegisterClass("60D",ImplicitGeometry60D)
[docs]class ImplicitGeometryHoneycomb(ImplicitGeometryBase):
"""ImplicitGeometry Honycomb (in the xy plane)
lx : scale (1)
w : wall thickness (0.1)
translate : offset of (0 0 0)
"""
def __init__(self, lx=1., w=0.1 ):
super(ImplicitGeometryHoneycomb,self).__init__()
self.lx = float(lx)
self.w = float(w)
self.translate = np.array([0.,0.,0.], dtype= float)
def __str__(self):
res = "ImplicitGeometryHoneycomb \n"
res += " lx: " + str(self.lx)+"\n"
res += " width: "+str(self.w)+"\n"
res += " translate: "+str(self.translate)+"\n"
return res
[docs] def GetDistanceToPoint(self,pos):
v = math.sqrt(3.)
v2 = math.sqrt(1./3.)
llx = self.lx/2.
lw =self.w/2.
mpos = abs((np.mod((pos+self.translate)/llx,[2.,2*v,1.])-[1.,v,0.]) )
pA1 = ImplicitGeometryPlane(point=np.array([(lw/llx),0.,0.]), normal=np.array([1.,0.,0.,]))
pA2 = ImplicitGeometryPlane(point=[0.,v2,0.], normal=[0.,1.,0.])
pA = ImplicitGeometryIntersection([pA1,pA2])
normal = np.array([-v2,1.,0.])
pB1 =ImplicitGeometryPlane(point=[0.,v2,0.], normal=-normal, offset= lw/llx)
pB2 =ImplicitGeometryPlane(point=[0.,v2,0.], normal=normal, offset= lw/llx)
pB = ImplicitGeometryIntersection([pB1,pB2])
pC1 =ImplicitGeometryPlane(point=[1-(lw/llx),0.,0.], normal=np.array([-1.,0.,0.,]))
pC2 =ImplicitGeometryPlane(point=[1.,v-v2,0], normal=[0.,-1.,0.,])
pC = ImplicitGeometryIntersection([pC1,pC2])
#res = ImplicitGeometryUnion([pA,pC]).GetDistanceToPoint(mpos)
res = ImplicitGeometryUnion([pA,pB,pC]).GetDistanceToPoint(mpos)
#res = ImplicitGeometryUnion([pA,pB,pC]).GetDistanceToPoint(mpos)
#res = ImplicitGeometryUnion([pA,pB,pC]).GetDistanceToPoint(mpos)
return self.ApplyInsideOut(res)
RegisterClass("Honeycomb",ImplicitGeometryHoneycomb)
import os
## hack to make work loading stl using point as decimal separator (,/.)
os.environ["LANG"] = "en_UK"
[docs]class ImplicitGeometryStl(ImplicitGeometryBase):
"""ImplicitGeometry based on a external stlfile
filename : stl filename to be loaded
option 'asVolume' is on as default this means the surface represent a close volume.
The algorithm will try to generate a close representation of the levelset.
if asVolume is False the a positive levelset is generated around the surface.
Please use an offset to add 'volume to the surface'
"""
def __init__(self):
super(ImplicitGeometryStl,self).__init__()
self.implicitFunction = None
self.filename = ''
self.boundingMin = [0,0,0]
self.boundingMax = [0,0,0]
self.surface = None
self.asVolume = None
[docs] def GetBoundingMin(self):
return self.boundingMin
[docs] def GetBoundingMax(self):
return self.boundingMax
[docs] def SetFileName(self,filenameSTL):
self.LoadFromFile(filenameSTL)
[docs] def LoadFromFile(self,filenameSTL):
from BasicTools.IO.IOFactory import InitAllReaders
InitAllReaders()
from BasicTools.IO.UniversalReader import ReadMesh
mesh = ReadMesh(filename=filenameSTL)
mesh.ConvertDataForNativeTreatment()
if mesh.GetNumberOfNodes() == 0:# pragma: no cover
raise ValueError( "No point data could be loaded from '" + filenameSTL)
return self.SetMesh(mesh)
[docs] def SetMesh(self,mesh):
# check if we have only element of dimensionality 2 or less
from BasicTools.Containers.Filters import ElementFilter
if mesh.GetElementsDimensionality() == 2:
return self.SetSurface(mesh)
from BasicTools.Containers.UnstructuredMeshInspectionTools import ExtractElementsByElementFilter
bulkmesh = ExtractElementsByElementFilter(mesh,ElementFilter(mesh,dimensionality=mesh.GetElementsDimensionality()),copy=True)
from BasicTools.Containers.UnstructuredMeshModificationTools import ComputeSkin
ComputeSkin(bulkmesh, inPlace= True)
skin = ExtractElementsByElementFilter(bulkmesh,ElementFilter(bulkmesh, dimensionality=bulkmesh.GetElementsDimensionality()-1),copy=False)
return self.SetSurface(skin, asVolume=True)
[docs] def SetSurface(self, mesh, asVolume = None):
# check if we have element of dimensionality 3
if mesh.GetElementsDimensionality() == 3:
return self.SetMesh(mesh)
if asVolume is None:
from BasicTools.Containers.UnstructuredMeshModificationTools import ComputeSkin
if ComputeSkin(mesh, md=2, inPlace=False).GetNumberOfElements() == 0 :
self.asVolume = True
else:
self.asVolume = False
else:
self.asVolume = asVolume
self.surface = mesh
from BasicTools.Containers.NativeTransfer import NativeTransfer
from BasicTools.Containers.Filters import ElementFilter
self.__nativeTransfer = NativeTransfer()
self.__nativeTransfer.SetVerbose(False)
from BasicTools.FE.Fields.FEField import FEField
from BasicTools.FE.FETools import PrepareFEComputation
space, numberings, offset, NGauss = PrepareFEComputation(mesh,numberOfComponents=1)
self.__field = FEField("", mesh=mesh, space=space, numbering=numberings[0])
self.__nativeTransfer.SetSourceFEField(self.__field, ElementFilter(mesh, dimensionality=mesh.GetElementsDimensionality() ) )
from BasicTools.NumpyDefs import PBasicFloatType
self.normals = np.zeros((mesh.GetNumberOfNodes(),3), dtype=PBasicFloatType)
nodes = mesh.nodes
for name, data, ids in ElementFilter(mesh, dimensionality=mesh.GetElementsDimensionality() ):
if data.connectivity.shape[1] >=3 :
for i in range(data.GetNumberOfElements()):
p0 = nodes[data.connectivity[i,0],:]
p1 = nodes[data.connectivity[i,1],:]
p2 = nodes[data.connectivity[i,2],:]
normal = np.cross(p1-p0,p2-p0)
normal = normal/np.linalg.norm(normal)
self.normals[data.connectivity[i,:],:] += [normal]
else:
for i in range(data.GetNumberOfElements()):
dp = nodes[data.connectivity[i,1],:] - nodes[data.connectivity[i,0],:]
normal = [dp[1], - dp[0]]
normal = normal/np.linalg.norm(normal)
self.normals[data.connectivity[i,:],:2] += [normal]
self.normals /= np.linalg.norm(self.normals, axis=1)[:,None]
[docs] def GetDistanceToPoint(self, pos):
if self.asVolume:
self.__nativeTransfer.SetTransferMethod("Interp/Extrap")
else:
self.__nativeTransfer.SetTransferMethod("Interp/Clamp")
if len(pos.shape) == 1:
targetPoints = np.empty((1,3),dtype=float)
targetPoints[0,0:pos.shape[0]] = pos
elif pos.shape[1] != 3:
targetPoints = np.empty((pos.shape[1],3),dtype=float)
targetPoints[:,0:pos.shape[0]] = pos
else:
targetPoints = pos
self.__nativeTransfer.SetTargetPoints(targetPoints)
self.__nativeTransfer.Compute()
op = self.__nativeTransfer.GetOperator()
vdist = op.dot(self.__field.mesh.nodes)- pos
dist = np.sqrt(np.sum((vdist )**2,axis=1))
if self.asVolume:
tnormals = op.dot(self.normals)
dist *= np.sign( np.sum(vdist*tnormals,axis=1))
dist *= -1.0
return self.ApplyInsideOut(dist)
def __str__(self):
res = "ImplicitGeometryStl \n"
res += " fileName : " + str(self.filename)+"\n"
res += " InsideOut : " + str(self.insideOut)+"\n"
res += " surface : " + str(self.surface)+"\n"
res += " asVolume : " + str(self.asVolume)+"\n"
return res
[docs]def CreateImplicitGeometryStl(ops):
res = ImplicitGeometryStl()
if "filename" in ops:
res.LoadFromFile(ops["filename"])
del ops["filename"]
PH.ReadProperties(ops,ops.keys(),res)
return res
RegisterClass("StlFile",ImplicitGeometryStl,CreateImplicitGeometryStl)
[docs]class ImplicitGeometryAnalytical(ImplicitGeometryBase):
"""ImplicitGeometry based on a function f(x,y,z) or f(pos)
expr : function expression
"""
def __init__(self):
super(ImplicitGeometryAnalytical,self).__init__()
self.expression = "0"
[docs] def SetExpression(self,string):
self.expression = string
[docs] def GetDistanceToPoint(self,pos):
if len(pos.shape) == 1:
res = np.zeros(1,dtype=float)
x, y, z = pos
else:
res = np.zeros(pos.shape[0],dtype=float)
x = pos[:,0]
y = pos[:,1]
z = pos[:,2]
from sympy import sympify
ls = {"x":x,"y":y,"z":z,"pos":pos}
res = np.array(sympify(self.expression,ls),dtype=float)
return self.ApplyInsideOut(res)
def __str__(self):
res = "ImplicitGeometryAnalytical \n"
res += " expression : " + str(self.expression)+"\n"
res += " InsideOut : " + str(self.insideOut)+"\n"
return res
[docs]def CreateImplicitGeometryAnalytical(ops):
res = ImplicitGeometryAnalytical()
if "expr" in ops:
res.SetExpression(ops["expr"])
return res
RegisterClass("Analytical",ImplicitGeometryAnalytical,CreateImplicitGeometryAnalytical)
[docs]def InitHoles(ls, nx, ny, nz, r):
nNodes = ls.support.GetDimensions()
grids = np.ix_(
np.linspace(0.0, 2 * nx * np.pi, nNodes[0]),
np.linspace(0.0, 2 * ny * np.pi, nNodes[1]),
np.linspace(0.0, 2 * nz * np.pi, nNodes[2]))
ls.phi = -np.cos(grids[0]) * np.cos(grids[1]) * np.cos(grids[2]) + r - 1.0
#self.phi = 0.2 * np.ceil(np.maximum(self.phi, 0.0)) - 0.1
ls.phi.shape = (np.prod(ls.phi.shape),)
[docs]class ImplicitGeometryHoles(ImplicitGeometryBase):
"""ImplicitGeometry to create holes
n : number of hole in each direction (3 3 3)
r : radius of holes (0.5)
offset : offset (0 0 0 )
support : (grid) to compute the bounding box to generate the holes
"""
def __init__(self):
super(ImplicitGeometryHoles,self).__init__()
self.r = 0.5
self.n = np.array([3.,3.,3.])
self.offset = np.array([0.,0.,0.])
self.type= ''
self.boundingMin = np.array([0.,0.,0.])
self.boundingMax = np.array([1.,1.,1.])
[docs] def SetSupport(self,support):
support.ComputeBoundingBox()
self.boundingMin = support.boundingMin[:]
self.boundingMax = support.boundingMax[:]
[docs] def GetBoundingMin(self):
return self.boundingMin
[docs] def GetBoundingMax(self):
return self.boundingMax
[docs] def SetNumberOfHoles(self,data):
if len(data) != 3:
raise
self.n = np.array(data,dtype=int)
[docs] def GetDistanceToPoint(self,pos):
l = self.boundingMax - self.boundingMin
l[l ==0] = 1.
l = list(l)
if len(l) == 2:
l.append(1.)
if self.type == "Original":
if len(pos.shape) == 1:
res = np.zeros(1,dtype=float)
res[0] = -np.cos(pos[0]) * np.cos(pos[1]) * np.cos(pos[2]) + self.r - 1.0
else:
res = -np.cos(self.offset[0] + (self.n[0]*np.pi/l[0])*pos[:,0]) * \
np.cos(self.offset[1] + (self.n[1]*np.pi/l[1])*pos[:,1]) * \
np.cos(self.offset[2] + (self.n[2]*np.pi/l[2])*pos[:,2]) + self.r - 1.0
elif self.type == "Aligned":
balls = [ ImplicitGeometrySphere(radius=self.r,center=[x,y,z])
for x in np.linspace(self.boundingMin[0],self.boundingMax[0], np.int(self.n[0]))
for y in np.linspace(self.boundingMin[1],self.boundingMax[1], np.int(self.n[1]))
for z in np.linspace(self.boundingMin[2],self.boundingMax[2], np.int(self.n[2]))]
Ores = ImplicitGeometryUnion(balls)
Ores.insideOut=True
res = Ores.GetDistanceToPoint(pos)
else:
dl = l/self.n
mpos = abs(np.mod((pos+self.offset)+dl,2*dl))-dl
sA = ImplicitGeometrySphere(radius=self.r,center=[0,0,0])
sB0 = ImplicitGeometrySphere(radius=self.r,center=[dl[0],dl[1],0])
sB1 = ImplicitGeometrySphere(radius=self.r,center=[-dl[0],-dl[1],0])
sB2 = ImplicitGeometrySphere(radius=self.r,center=[dl[0],-dl[1],0])
sB3 = ImplicitGeometrySphere(radius=self.r,center=[-dl[0],dl[1],0])
sC0 = ImplicitGeometrySphere(radius=self.r,center=[dl[0],0,dl[2]])
sC1 = ImplicitGeometrySphere(radius=self.r,center=[-dl[0],0,-dl[2]])
sC2 = ImplicitGeometrySphere(radius=self.r,center=[-dl[0],0,dl[2]])
sC3 = ImplicitGeometrySphere(radius=self.r,center=[dl[0],0,-dl[2]])
sD0 = ImplicitGeometrySphere(radius=self.r,center=[0,dl[1],dl[2]])
sD1 = ImplicitGeometrySphere(radius=self.r,center=[0,-dl[1],-dl[2]])
sD2 = ImplicitGeometrySphere(radius=self.r,center=[0,-dl[1],dl[2]])
sD3 = ImplicitGeometrySphere(radius=self.r,center=[0,dl[1],-dl[2]])
Ores = ImplicitGeometryUnion([sA,sB0,sB1,sB2,sB3,sC0,sC1,sC2,sC3,sD0,sD1,sD2,sD3])
Ores.insideOut=True
res = Ores.GetDistanceToPoint(mpos)
return self.ApplyInsideOut(res)
def __str__(self):
res = "ImplicitHoles \n"
res += " self.r : " + str(self.r)+"\n"
res += " InsideOut : " + str(self.insideOut)+"\n"
return res
RegisterClass("Holes",ImplicitGeometryHoles)
[docs]def CheckIntegrity(GUI=False):
def MustFail(func):
try:
func()
raise #pragma: no cover
except:
pass
from BasicTools.Containers.ConstantRectilinearMesh import ConstantRectilinearMesh
myMesh = ConstantRectilinearMesh()
myMesh.SetDimensions([5,6,7])
myMesh.SetSpacing(1./(myMesh.GetDimensions()-1)*2)
myMesh.SetOrigin([-1.,-1.,-1.])
myMesh.elements["hex8"].tags.CreateTag("2elems").SetIds([0,1])
myMesh.nodesTags.CreateTag("3points").SetIds([2,3,4])
print(myMesh)
OnePoint3D = np.array([1,2,3])
TwoPoints3D = np.array([[0.,0.,0.],[1.,2.,3.]], dtype=float)
import BasicTools.TestData as TestData
from BasicTools.Helpers.Tests import TestTempDir
testDataPath = TestData.GetTestDataPath()
##################### ImplicitGeometrySphere ################################
IGSphere = ImplicitGeometrySphere(center=[0,0,0],radius=0.1)
IGSphere = ImplicitGeometrySphere()
IGSphere.center = np.array([.5,-0.3,0.5])
IGSphere.radius = 0.7
IGSphere_Data = IGSphere.ApplyVector(myMesh)
print(IGSphere)
print(IGSphere.GetBoundingMin() )
print(IGSphere.GetBoundingMax() )
IGSphere(np.array([1,2,3]))
###########################################################################
IGWrapped = ImplicitGeometryWrapped(IGSphere_Data)
IGWrapped(myMesh)
ImplicitGeometryWrapped(TwoPoints3D[:,0])(TwoPoints3D)
from functools import partial
MustFail(partial(ImplicitGeometryWrapped(TwoPoints3D[0:1,0]),TwoPoints3D))
MustFail(partial(ImplicitGeometryWrapped(TwoPoints3D[0:1,0]),myMesh,cellCenter=True))
MustFail(partial(ImplicitGeometryWrapped(TwoPoints3D[0:1,0]),myMesh,cellCenter=False))
#########################ImplicitGeometryExternalSurface#########################
IGExternalSurface = CreateImplicitGeometryExternalSurface({"support":myMesh})
IGExternalSurface(TwoPoints3D)
############################ImplicitGeometryByETag###############################
IGByETag = CreateImplicitGeometryByETag({"support":myMesh, "eTags":["2elems"]})
print(IGByETag)
IGByETag(TwoPoints3D)
IGByETag.GetDistanceToPoint(TwoPoints3D)
class LS():
pass
LS.support = myMesh
IGByETag = CreateImplicitGeometryByETag({"ls":LS, "eTags":["2elems"]})
MustFail(partial(CreateImplicitGeometryByETag,{"support":myMesh} ))
MustFail(partial(CreateImplicitGeometryByETag,{"eTags":["2elems"]} ))
######################### ImplicitGeometryAxisAlignBox ###################
IGAxisAlignBox = ImplicitGeometryAxisAlignBox()
IGAxisAlignBox = ImplicitGeometryAxisAlignBox(origin=[0,1,2],size=[0.1,0.2,0.3])
IGAxisAlignBox = ImplicitGeometryAxisAlignBox(origin=[0,1,2],size=3.)
IGAxisAlignBox.GetBoundingMin()
IGAxisAlignBox.GetBoundingMax()
IGAxisAlignBox.SetSupport(myMesh)
IGAxisAlignBox(myMesh)
print(IGAxisAlignBox)
###################### ImplicitGeometrySphereFromNTag ###################
IGSphereFromNTag= CreateImplicitGeometrySphereFromNTag({"radius":0.1,"nTag":"3points","ls":LS})
IGSphereFromNTag.SetLs(LS)
IGSphereFromNTag(myMesh)
IGSphereFromNTag(OnePoint3D)
print(IGSphereFromNTag)
###################### ImplicitGeometryCylinder ###################
IGCylinder = ImplicitGeometryCylinder(wcups=False)
IGCylinder(OnePoint3D)
IGCylinder = ImplicitGeometryCylinder(center1=[0,0,0], center2=[1,2,3],radius=0.1)
IGCylinder(TwoPoints3D)
print(IGCylinder)
IGCylinder.GetBoundingMin()
IGCylinder.GetBoundingMax()
###################### ImplicitGeometryPlane ###################
IGGeometryPlane = CreateImplicitGeometryPlane({"point":[0,0,0],"normal":[1,1,1]})
print(IGGeometryPlane)
###################### ImplicitGeometryGyroid ###################
IGGyroid = ImplicitGeometryGyroid()
IGGyroid = ImplicitGeometryGyroid(scale=0.1,offset=[0.,0.2,0.3])
print(IGGyroid)
IGGyroid(myMesh)
IGGyroid.type = 1
IGGyroid(myMesh)
IGGyroid.type = 2
IGGyroid(myMesh)
IGGyroid.wall = True
IGGyroid(myMesh)
###################### ImplicitGeometry60D ###################
IG60D = ImplicitGeometry60D()
print(IG60D)
IG60D(myMesh)
###################### ImplicitGeometryHoneycomb ###################
IGHoneycomb = ImplicitGeometryHoneycomb()
print(IGHoneycomb)
IGHoneycomb(myMesh)
###################### ImplicitGeometryStl ###################
IGStl = CreateImplicitGeometryStl({"filename":testDataPath+"stlsphere.stl"})
IGStl.SetFileName(testDataPath+"stlsphere.stl")
IGStl.GetBoundingMin()
IGStl.GetBoundingMax()
IGStl(myMesh)
print(IGStl)
######################## ImplicitGeometryStl with open stl #######################
IGStl = CreateImplicitGeometryStl({"filename":testDataPath+"stlOpenSphere.stl"})
dist = IGStl(myMesh)
myMesh.nodeFields["dist"] = dist
if GUI:
from BasicTools.Actions.OpenInParaView import OpenInParaView
OpenInParaView(myMesh,filename=TestTempDir().GetTempPath()+ "LsOpenSphere.xdmf",run=False)
print(IGStl)
###################### ImplicitGeometryAnalytical ###################
IGAnalytical = CreateImplicitGeometryAnalytical({"expr":"2*y+1+z"})
IGAnalytical(myMesh)
IGAnalytical(OnePoint3D)
print(IGAnalytical)
###################### ImplicitGeometryHoles ###################
IGHoles = ImplicitGeometryHoles()
IGHoles.SetSupport(myMesh)
IGHoles.GetBoundingMin()
IGHoles.GetBoundingMax()
IGHoles.SetNumberOfHoles([2,3,4])
IGHoles(myMesh)
print(IGHoles)
IGHoles.type = "Original"
IGHoles(myMesh)
IGHoles(OnePoint3D)
MustFail(partial(IGHoles.SetNumberOfHoles,[1,2]))
############################ImplicitGeometryByETagII###############################
IGByETagII = CreateImplicitGeometryByETagII({"support":myMesh, "eTags":["2elems"]})
print(IGByETagII)
ETagII = IGByETagII(myMesh)
return "ok"
if __name__ == '__main__':# pragma: no cover
print(CheckIntegrity(GUI=True))