Source code for Muscat.ImplicitGeometry.ImplicitGeometryObjects

# -*- 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, MuscatFloat
import Muscat.MeshContainers.ElementsDescription as ED
import Muscat.Helpers.ParserHelper as PH
from Muscat.MeshContainers.Mesh import Mesh

from Muscat.ImplicitGeometry.ImplicitGeometryFactory import RegisterClass
from Muscat.ImplicitGeometry.ImplicitGeometryBase import ImplicitGeometryBase, dsin, dcos
from Muscat.ImplicitGeometry.ImplicitGeometryOperators import ImplicitGeometryIntersection, ImplicitGeometryUnion
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.FE.FETools import ComputeNormalsAtPoints


[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=""] > """ method = ops.get("method","bulk") res = ImplicitGeometryExternalSurface(method=method) 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): def __init__(self, support: Mesh = None,method: str ="bulk"): """ImplicitGeometry of the external surface of a mesh (support) Parameters ---------- support : Mesh, optional mesh to compute the surface, by default None method : str, optional method to use in order to compute distance, either "bulk" either "normal", by default "bulk" the "bulk" method is the legacy method using two field transfers whereas the "normal" method uses one field transfer and one skin integral that should be more efficient """ super().__init__() assert method in ["bulk","normal"] self._method=method self.offset = 1E-3 if support is not None: self.SetSupport(support)
[docs] def SetSupport(self, support: Mesh): from Muscat.MeshTools.TransferBackEnds.NativeTransfer import NativeTransfer from Muscat.FE.Fields.FEField import FEField from Muscat.FE.FETools import PrepareFEComputation from Muscat.MeshTools.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") if self.method=="bulk": self.__nativeTransfer = NativeTransfer() self.__nativeTransfer.SetVerbose(False) self.__nativeTransfer.SetUseEdgeSearch(True) 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")
@property def method(self): return self._method
[docs] def GetDistanceToPoint(self, pos: np.ndarray) -> np.ndarray: self.__SnativeTransfer.SetTargetPoints(pos) self.__SnativeTransfer.Compute() op = self.__SnativeTransfer.GetOperator() vectorDistance = np.zeros( (pos.shape[0], max(pos.shape[1], self.__Sfield.mesh.nodes.shape[1] )) , dtype=MuscatFloat) transfered_nodes = op.dot(self.__Sfield.mesh.nodes) vectorDistance[:,0:self.__Sfield.mesh.nodes.shape[1]] = transfered_nodes vectorDistance[:,0:pos.shape[1]] -= pos distance = np.sqrt(np.sum((vectorDistance )**2,axis=1)) if self.method == "bulk": self.__nativeTransfer.SetTargetPoints(pos) self.__nativeTransfer.Compute() # status must be 1 inside or 3 Clamp not other status = self.__nativeTransfer.GetStatus() 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 return distance*(status.flatten()-2) + self.offset elif self.method == "normal": vectorNormal = np.zeros( (pos.shape[0], max(pos.shape[1], self.__Sfield.mesh.nodes.shape[1] )) , dtype=MuscatFloat) normals = ComputeNormalsAtPoints(self.__Sfield.mesh) vectorNormal[:,0:self.__Sfield.mesh.nodes.shape[1]] = op.dot(normals) signDistance = np.sign(np.sum(vectorDistance*vectorNormal,axis=-1)) return -1*signDistance*distance + 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.MeshTools.MeshInspectionTools import ExtractElementsByElementFilter from Muscat.MeshContainers.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.MeshContainers.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.MeshContainers.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.MeshContainers.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.MeshContainers.Filters.FilterObjects import ElementFilter if mesh.GetElementsDimensionality() == 2: return self.SetSurface(mesh) from Muscat.MeshTools.MeshInspectionTools import ExtractElementsByElementFilter bulkmesh = ExtractElementsByElementFilter(mesh, ElementFilter(dimensionality=mesh.GetElementsDimensionality()), copy=True) from Muscat.MeshTools.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.MeshTools.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.MeshTools.TransferBackEnds.NativeTransfer import NativeTransfer from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter self.__nativeTransfer = NativeTransfer() self.__nativeTransfer.SetVerbose(False) self.__nativeTransfer.SetUseEdgeSearch(True) 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())) 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 = np.zeros( (pos.shape[0], max(pos.shape[1], self.__field.mesh.nodes.shape[1] )) , dtype=MuscatFloat) vdist[:,0:self.__field.mesh.nodes.shape[1]] = op.dot(self.__field.mesh.nodes) vdist[:,0:pos.shape[1]] -= 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.MeshTools.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh myMesh = CreateConstantRectilinearMesh(dimensions=[5, 6, 7], spacing=[2/.4, 2/.5, 2/.6], origin=[-1., -1., -1.]) myMesh.elements[ED.Hexahedron_8].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) myMesh2D = CreateConstantRectilinearMesh(dimensions=[4, 5], spacing=[2/.3, 2/.4], origin=[-1., -1.]) OnePoint2D = np.array([-0.1, 1.]) TwoPoints2D = np.array([[1., -1.], [1., 4.]], 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(TwoPoints2D) IGExternalSurface2D = CreateImplicitGeometryExternalSurface({"support": myMesh2D}) IGExternalSurface2D(TwoPoints3D) ################ ImplicitGeometryExternalSurface Normal method ################### IGExternalSurfaceNormalMethod = CreateImplicitGeometryExternalSurface({"support": myMesh,"method":"normal"}) IGExternalSurfaceNormalMethod(TwoPoints2D) IGExternalSurfaceNormalMethod2D = CreateImplicitGeometryExternalSurface({"support": myMesh2D,"method":"normal"}) IGExternalSurfaceNormalMethod2D(TwoPoints3D) ## Testing if methods make the same result assert np.allclose(IGExternalSurfaceNormalMethod(TwoPoints3D),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))