# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
from typing import Optional, Union
import numpy as np
from Muscat.Types import ArrayLike, MuscatFloat
from Muscat.ImplicitGeometry.ImplicitGeometryFactory import RegisterClass
from Muscat.Helpers.Logger import Debug
[docs]
def dsin(x):
r = np.mod(x + np.pi / 2, 2 * np.pi) - np.pi
return -(abs(r) - np.pi / 2)
[docs]
def dcos(x):
return dsin(x + np.pi / 2)
[docs]
def smoothmin(a, b, k=1):
h = np.maximum(k - np.abs(a - b), 0.0) / k
return np.minimum(a, b) - h * h * k * (1.0 / 4.0)
[docs]
def smoothmax(a, b, k=1):
return -smoothmin(-a, -b, k=k)
[docs]
class ImplicitGeometryBase:
def __init__(self):
super().__init__()
self.insideOut = False
[docs]
def ApplyVector(self, support, cellCenter=False, dim=None):
if type(support).__module__ == np.__name__:
if len(support.shape) == 2 and support.shape[1] == 2:
support = np.hstack((support, np.zeros((support.shape[0], 1), dtype=support.dtype)))
d = self.GetDistanceToPoint(support)
else:
if cellCenter:
from Muscat.Containers.MeshTools import GetElementsCenters
pos = GetElementsCenters(support, dim=dim)
else:
pos = support.GetPosOfNodes()
if pos.shape[1] == 2:
pos = np.hstack((pos, np.zeros((pos.shape[0], 1), dtype=pos.dtype)))
d = self.GetDistanceToPoint(pos)
return d
[docs]
def GetGradientDistanceToPoint(self, pos: np.ndarray, dx: Optional[Union[MuscatFloat, ArrayLike]] = None) -> np.ndarray:
"""Compute the numerical gradient of the implicit geometry
the child classes can overwrite it for a more efficient or exact result.
Parameters
----------
pos : np.ndarray
the position to evaluate the gradient
dx : Optional[MuscatFloat,Arraylike], optional
the step (per coordinate) to be used to compute the numerical gradient, by default dx is 1e-6
Returns
-------
np.ndarray
a matrix of size (nb point, nb coordinates) with the gradient for each coordinate (d_distance/dx)
"""
if len(pos.shape) == 1:
pos = np.asarray([pos], dtype=MuscatFloat)
ndim = pos.shape[1]
if dx is None:
dx = [1.0e-6] * ndim
elif not hasattr(dx, "__len__"):
dx = [dx] * ndim
if pos.shape[1] == 2:
pos = np.hstack((pos, np.zeros((pos.shape[0], 1), dtype=pos.dtype)))
res = np.zeros_like(pos)
for i in range(ndim):
if dx[i] != 0:
pos[:, i] -= dx[i]
res[:, i] = -self.GetDistanceToPoint(pos)
pos[:, i] += 2 * dx[i]
res[:, i] += self.GetDistanceToPoint(pos)
pos[:, i] -= dx[i]
res[:, i] /= 2 * dx[i]
return res[:,0:ndim]
[docs]
def ApplyInsideOut(self, res):
if self.insideOut:
return -res
return res
def __call__(self, support, cellCenter=False):
return self.ApplyVector(support, cellCenter=cellCenter)
[docs]
def GetDistanceToPoint(self, pos):
raise (Exception("Please overload this function"))
# def __str__(self):
# return "+-+-"
####################### cache object ################################
[docs]
class ImplicitGeometryCachedData(ImplicitGeometryBase):
def __init__(self, internalImplicitGeometry):
super().__init__()
self.internalImplicitGeometry = internalImplicitGeometry
self.cacheData = None
self.cachedInputGetDistanceToPoint = None
self.cachedInputApplyVector = None
self.cachedInputcellCenter = False
self.cachedInputcelldim = 0
[docs]
def GetDistanceToPoint(self, pos):
if self.cacheData is None:
Debug("building Cache")
self.cacheData = self.internalImplicitGeometry.GetDistanceToPoint(pos)
elif pos is not self.cachedInputGetDistanceToPoint:
Debug("rebuilding Cache")
self.cacheData = self.internalImplicitGeometry.GetDistanceToPoint(pos)
else:
pass
self.cachedInputGetDistanceToPoint = pos
self.cachedInputApplyVector = None
return self.cacheData
[docs]
def ApplyVector(self, support, cellCenter=False, dim=None):
if self.cacheData is not None:
if support is self.cachedInputApplyVector and cellCenter == self.cachedInputcellCenter and dim == self.cachedInputcelldim:
return self.cacheData
res = super().ApplyVector(support, cellCenter=cellCenter, dim=dim)
self.cachedInputApplyVector = support
self.cachedInputcellCenter = cellCenter
self.cachedInputcelldim = dim
return res
####################### cache object ################################
[docs]
class ImplicitGeometryDelayedInit(ImplicitGeometryBase):
def __init__(self, name, ops={}):
super().__init__()
self.internalImplicitGeometry = None
self.__name = name
self.__ops = ops
[docs]
def Init(self):
if self.internalImplicitGeometry is None:
from Muscat.ImplicitGeometry.ImplicitGeometryFactory import Create
self.internalImplicitGeometry = Create(self.__name, self.__ops)
## release memory
self.__name = None
self.__ops = None
if self.internalImplicitGeometry is None:
raise (Exception("Error in the initialisation of : " + str(self.__name))) # pragma: no cover
[docs]
def ApplyVector(self, support, cellCenter=False, dim=None):
self.Init()
return self.internalImplicitGeometry.ApplyVector(support, cellCenter=cellCenter)
[docs]
def GetDistanceToPoint(self, pos):
self.Init()
return self.internalImplicitGeometry.GetDistanceToPoint(pos)
def __str__(self):
res = "ImplicitGeometryDelayedInit :\n"
if self.__name is None:
res += self.internalImplicitGeometry.__str__()
else:
res = str(self.__name) + "\n"
res += str(self.__ops) + "\n"
return res
[docs]
def TryToCreate(name, ops=None):
from Muscat.ImplicitGeometry.ImplicitGeometryFactory import Create
res = Create(name, ops)
if res is None:
raise Exception("Unable to Create a " + str(name)) # pragma: no cover
return res
[docs]
def CheckIntegrity(GUI: bool = False):
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
myMesh3D = CreateConstantRectilinearMesh(dimensions=[30, 30, 30], spacing=[2 / 29] * 3, origin=[-1] * 3)
myMesh2D = CreateConstantRectilinearMesh(dimensions=[30, 30], spacing=[2 / 29] * 2, origin=[-1] * 2)
class DummyImplicitGeometry(ImplicitGeometryBase):
def __init__(self):
super().__init__()
def GetDistanceToPoint(self, pos):
return self.ApplyInsideOut(np.zeros(pos.shape[0]))
RegisterClass("Dummy", DummyImplicitGeometry, withError=False)
TryToCreate("Dummy")(myMesh3D)
TryToCreate("Dummy")(myMesh2D)
TryToCreate("Dummy").ApplyVector(myMesh3D, cellCenter=True)
TwoPoints3D = np.array([[0.0, 0.0, 0.0], [1.0, 2.0, 3.0]], dtype=float)
TryToCreate("Dummy")(TwoPoints3D)
try:
obj = ImplicitGeometryBase()
obj.GetDistanceToPoint(TwoPoints3D)
raise # pragma: no cover
except:
pass
TwoPoints2D = np.array([[0.0, 0.0], [1.0, 2.0]], dtype=float)
obj = TryToCreate("Dummy")
obj.insideOut = True
obj(TwoPoints2D)
####################################################################
myObj6 = ImplicitGeometryCachedData(obj)
myObj6(myMesh3D)
myObj6(myMesh3D)
myObj6.GetDistanceToPoint(TwoPoints3D)
myObj6.GetDistanceToPoint(TwoPoints3D)
myObj6.GetGradientDistanceToPoint(np.array([0, 0]))
myObj6.GetGradientDistanceToPoint(TwoPoints3D)
myObj6.GetGradientDistanceToPoint(TwoPoints3D, 1)
#######################################################################
res = ImplicitGeometryDelayedInit("Dummy")
print(res)
res(myMesh3D)
res.GetDistanceToPoint(TwoPoints3D)
print(res)
return "ok"
if __name__ == "__main__": # pragma: no cover
print(CheckIntegrity(GUI=False))