# -*- 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 os
import numpy as np
import math
from Muscat.Types import MuscatIndex
import Muscat.Helpers.ParserHelper as PH
from Muscat.ImplicitGeometry.ImplicitGeometryFactory import RegisterClass
from Muscat.ImplicitGeometry.ImplicitGeometryBase import ImplicitGeometryBase, dsin, dcos
from Muscat.ImplicitGeometry.ImplicitGeometryOperators import ImplicitGeometryIntersection, ImplicitGeometryUnion
from Muscat.Containers.Filters.FilterObjects 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().__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().__init__()
self.offset = 1E-3
if support is not None:
self.SetSupport(support)
[docs] def SetSupport(self, support):
from Muscat.Containers.NativeTransfer import NativeTransfer
self.__nativeTransfer = NativeTransfer()
self.__nativeTransfer.SetVerbose(False)
from Muscat.FE.Fields.FEField import FEField
from Muscat.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(dimensionality=support.GetElementsDimensionality() ) )
self.__nativeTransfer.SetTransferMethod("Interp/Clamp")
from Muscat.Containers.MeshModificationTools 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().__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.
self.dim = None
[docs] def SetSupportAndZones(self, support, etags):
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.Containers.Filters.FilterObjects import ElementFilter
sup = ExtractElementsByElementFilter(support, elementFilter=ElementFilter(eTag= 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 = self.__class__.__name__ + "\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 = []
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().__init__()
self.offset = 0.
from Muscat.Containers.Filters.FilterObjects import ElementFilter
self.elementFilter = ElementFilter()
self.offset = 0.
self.numberOfElementPseudoDistance = 10
self.support = None
[docs] def SetSupportAndZones(self, support, etags, dim=None):
self.support = support
self.elementFilter.SetETags(etags)
self.elementFilter.dimensionality = dim
def __str__(self):
res = self.__class__.__name__ + "\n"
res += " " + str(self.elementFilter)
return res
[docs] def ApplyVector(self, support, cellCenter=False):
if cellCenter:
raise (Exception("Not implemented"))
from Muscat.Containers.Filters.FilterTools import ElementFilterToImplicitField
vals = ElementFilterToImplicitField(self.elementFilter, self.support,self.numberOfElementPseudoDistance)
res = self.ApplyInsideOut(vals+(self.offset))
return res
[docs] def GetDistanceToPoint(self, _pos):
if _pos is self.support.nodes:
from Muscat.Containers.Filters.FilterTools import ElementFilterToImplicitField
vals = ElementFilterToImplicitField(self.elementFilter, self.support, 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().__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 = self.__class__.__name__ + "\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):
boundingMin, boundingMax = support.ComputeBoundingBox()
self.origin = boundingMin
self.size = boundingMax - 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().__init__()
self.SetRadius(radius)
self.centers = None
self.nTag = None
[docs] def SetLs(self, ls):
self.SetSupport(ls.support)
[docs] def SetSupport(self, support):
if len(support.nodes.shape) == 2 and support.nodes.shape[1] ==2:
support3D = np.hstack((support.nodes,np.zeros((support.nodes.shape[0],1) ,dtype=support.nodes.dtype ) ) )
self.centers = support3D[support.nodesTags[self.nTag].GetIds(), :]
else:
self.centers = support.nodes[support.nodesTags[self.nTag].GetIds(), :]
[docs] def SetRadius(self, val):
self.radius = val
def __str__(self):
res = self.__class__.__name__ + "\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().__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 = self.__class__.__name__ + "\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 ImplicitGeometryEllipse(ImplicitGeometrySphere):
"""ImplicitGeometry based on a Ellipse
You can set one of the radius to be infinite so as to have an infinite ellipse on one dimension
- center : center of the ellipse (X Y Z)
- radius : radius of the ellipse (a b c)
"""
def __init__(self, radius=[1.,1.,1.], center=None):
super().__init__(radius,center)
[docs] def GetDistanceToPoint(self, _pos):
if len(_pos.shape) == 1:
pos = np.array([_pos])
else:
pos = _pos
# the scaling factor here is taken as the max radius that is not infinite
# otherwise to have the real distance to the ellipse one should solve a convex problem
res = (np.sqrt(np.sum(((pos-self.center)/self.radius)**2, axis=1))-1)*np.min(self.radius)
return self.ApplyInsideOut(res)
RegisterClass("Ellipse", ImplicitGeometryEllipse)
[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().__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 = self.__class__.__name__ + "\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().__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 = self.__class__.__name__ + "\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().__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 = self.__class__.__name__ + "\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().__init__()
self.lx = lx
self.w = w
def __str__(self):
res = self.__class__.__name__ + "\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().__init__()
self.lx = float(lx)
self.w = float(w)
self.translate = np.array([0., 0., 0.], dtype=float)
def __str__(self):
res = self.__class__.__name__ + "\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)
# 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().__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 Muscat.IO.IOFactory import InitAllReaders
InitAllReaders()
from Muscat.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 Muscat.Containers.Filters.FilterObjects import ElementFilter
if mesh.GetElementsDimensionality() == 2:
return self.SetSurface(mesh)
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
bulkmesh = ExtractElementsByElementFilter(mesh, ElementFilter(dimensionality=mesh.GetElementsDimensionality()), copy=True)
from Muscat.Containers.MeshModificationTools import ComputeSkin
ComputeSkin(bulkmesh, inPlace=True)
skin = ExtractElementsByElementFilter(bulkmesh, ElementFilter(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 Muscat.Containers.MeshModificationTools 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 Muscat.Containers.NativeTransfer import NativeTransfer
from Muscat.Containers.Filters.FilterObjects import ElementFilter
self.__nativeTransfer = NativeTransfer()
self.__nativeTransfer.SetVerbose(False)
from Muscat.FE.Fields.FEField import FEField
from Muscat.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(dimensionality=mesh.GetElementsDimensionality()))
from Muscat.Types import MuscatFloat
self.normals = np.zeros((mesh.GetNumberOfNodes(), 3), dtype=MuscatFloat)
usedNodes = np.zeros(mesh.GetNumberOfNodes(), dtype=bool)
nodes = mesh.nodes
for selection in ElementFilter(dimensionality=mesh.GetElementsDimensionality())(mesh):
if selection.elements.connectivity.shape[1] >= 3:
for i in range(selection.elements.GetNumberOfElements()):
p0 = nodes[selection.elements.connectivity[i, 0], :]
p1 = nodes[selection.elements.connectivity[i, 1], :]
p2 = nodes[selection.elements.connectivity[i, 2], :]
normal = np.cross(p1-p0, p2-p0)
normal = normal/np.linalg.norm(normal)
self.normals[selection.elements.connectivity[i, :], :] += [normal]
selection.elements.GetNodesIndexFor(selection.indices)
usedNodes[selection.elements.GetNodesIndexFor(selection.indices)] = True
else:
for i in range(selection.elements.GetNumberOfElements()):
dp = nodes[selection.elements.connectivity[i, 1], :] - nodes[selection.elements.connectivity[i, 0], :]
normal = [dp[1], -dp[0]]
normal = normal/np.linalg.norm(normal)
self.normals[selection.elements.connectivity[i, :], :2] += [normal]
usedNodes[selection.elements.GetNodesIndexFor(selection.indices)] = True
norm = np.linalg.norm(self.normals, axis=1)[:, None]
norm[usedNodes==0] = 1
self.normals /= norm
[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 = self.__class__.__name__ + "\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().__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 = self.__class__.__name__ + "\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().__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):
self.boundingMin, self.boundingMax = support.ComputeBoundingBox()
[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], MuscatIndex(self.n[0]))
for y in np.linspace(self.boundingMin[1], self.boundingMax[1], MuscatIndex(self.n[1]))
for z in np.linspace(self.boundingMin[2], self.boundingMax[2], MuscatIndex(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: bool = False):
from Muscat.Helpers.CheckTools import MustFailFunction
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
myMesh = CreateConstantRectilinearMesh(dimensions=[5, 6, 7], spacing=[2/.4, 2/.5, 2/.6], origin=[-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 Muscat.TestData as TestData
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
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]))
##################### ImplicitGeometryEllipse ################################
IGEllipse = ImplicitGeometryEllipse(center=[0, 0, 0], radius=[0.1,0.5,0.7])
IGEllipse = ImplicitGeometryEllipse()
IGEllipse.center = np.array([.5, -0.3, 0.5])
IGEllipse.radius = np.array([0.1,0.5,np.inf])
IGEllipse_Data = IGEllipse.ApplyVector(myMesh)
print(IGEllipse)
print(IGEllipse.GetBoundingMin())
print(IGEllipse.GetBoundingMax())
IGEllipse(np.array([1, 2, 3]))
###########################################################################
IGWrapped = ImplicitGeometryWrapped(IGSphere_Data)
IGWrapped(myMesh)
ImplicitGeometryWrapped(TwoPoints3D[:, 0])(TwoPoints3D)
from functools import partial
MustFailFunction(partial(ImplicitGeometryWrapped(TwoPoints3D[0:1, 0]), TwoPoints3D))
MustFailFunction(partial(ImplicitGeometryWrapped(TwoPoints3D[0:1, 0]), myMesh, cellCenter=True))
MustFailFunction(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"]})
MustFailFunction(partial(CreateImplicitGeometryByETag, {"support": myMesh}))
MustFailFunction(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 Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(myMesh, filename=TemporaryDirectory().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)
MustFailFunction(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))