# -*- 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
from Muscat.Types import MuscatFloat
[docs]class SpaceBase():
"""Base class for a finite element interpolation space for one type of element
"""
def __init__(self):
super().__init__()
self.name: str = ""
self.fromConnectivityCompatible: bool = False
self.posN = None
@property
def geoSupport(self):
return self._geoSupport
@geoSupport.setter
def geoSupport(self, gs):
self._geoSupport = gs
[docs] def GetDimensionality(self):
return self.geoSupport.dimensionality # pragma: no cover
[docs] def GetSpaceEvaluatedAt(self, points, weights):
raise NotImplementedError("GetSpaceEvaluatedAt not implemented")
[docs] def ClampParamCoordinates(self, xietaphi):
""" Clamps the xi eta and phi to othe intervals [0,1] for the first dim coordinates,
the others are set to zero """
res = xietaphi.copy()
for cpt, d in enumerate(xietaphi):
if cpt < self.GetDimensionality():
res[cpt] = max(0., d)
res[cpt] = min(1., res[cpt])
else:
res[cpt] = 0
return res
[docs] def GetNormal(self, Jacobian):
# Edge in 2D
if Jacobian.shape[0] == 1 and Jacobian.shape[1] == 2:
res = np.array([Jacobian[0, 1], -Jacobian[0, 0]], dtype=MuscatFloat)
# res = np.array([Jacobian[1,:] -Jacobian[0,:]],dtype =MuscatFloat) #ANCIENNE VERSION
# Edge in 3D, we return the xy projection of the normal
elif Jacobian.shape[0] == 1 and Jacobian.shape[1] == 3:
res = np.array([Jacobian[0, 1], -Jacobian[0, 0]], dtype=MuscatFloat)
# surface in 3D
elif Jacobian.shape[0] == 2 and Jacobian.shape[1] == 3:
res = np.cross(Jacobian[0, :], Jacobian[1, :])
else:
raise RuntimeError(f"cant compute the normal of an element with a jacobian of shape {Jacobian.shape}")
# normalization
res /= np.linalg.norm(res)
return res
[docs] def GetJackAndDet(self, Nfder, xcoor):
"""Deprecated Please use GetJacobianDeterminantAndInverse """
return self.GetJacobianDeterminantAndInverse(Nfder, xcoor)
[docs] def GetJacobianDeterminantAndInverse(self, Nfder: np.ndarray, xcoor: np.ndarray):
dim = self.GetDimensionality()
if dim == 0:
return np.array([[1.]]), 1, lambda x: x
Jacobian = np.dot(Nfder, xcoor)
s = xcoor.shape[1]
if dim == s:
Jdet = np.linalg.det(Jacobian)
if dim == 3:
def jinv(vec, jacobian=Jacobian):
m1, m2, m3, m4, m5, m6, m7, m8, m9 = jacobian.flatten()
determinant = m1*m5*m9 + m4*m8*m3 + m7*m2*m6 - m1*m6*m8 - m3*m5*m7 - m2*m4*m9
return np.dot(np.array([[m5*m9-m6*m8, m3*m8-m2*m9, m2*m6-m3*m5],
[m6*m7-m4*m9, m1*m9-m3*m7, m3*m4-m1*m6],
[m4*m8-m5*m7, m2*m7-m1*m8, m1*m5-m2*m4]]), vec) / determinant
return Jacobian, Jdet, jinv
elif dim == 2:
return Jacobian, Jdet, np.linalg.inv(Jacobian).dot
elif dim == 0:
Jdet = 1
elif dim == 1:
Jdet = np.linalg.norm(Jacobian)
elif dim == 2:
Jdet = np.linalg.norm(np.cross(Jacobian[0, :], Jacobian[1, :]))
q, r = np.linalg.qr(Jacobian)
qt = q.T
def jinv(vec, qt=qt, r=r): return np.linalg.lstsq(r, np.dot(qt, vec), rcond=None)[0]
return Jacobian, Jdet, jinv
[docs]class SpaceAtIntegrationPoints():
def __init__(self, space=None, points=None, weights=None):
if space is not None and points is not None and weights is not None:
self.SetIntegrationRule(space, points, weights)
else:
self.nbPoints = 0
self.space: SpaceBase = SpaceBase()
self.weights = np.empty((0), dtype=MuscatFloat)
self.points = np.empty((0, 3), dtype=MuscatFloat)
self.valN = [None]*self.nbPoints
self.valdphidxi = [None]*self.nbPoints
[docs] def SetIntegrationRule(self, space, points, weights):
self.space: SpaceBase = space
self.weights = weights
self.points = points
self.nbPoints = len(weights)
self.valN = [None]*self.nbPoints
self.valdphidxi = [None]*self.nbPoints
for pp in range(self.nbPoints):
point = points[pp]
self.valN[pp] = space.GetShapeFunc(point)
self.valdphidxi[pp] = space.GetShapeFuncDer(point)
[docs] def GetJacobianDeterminantAndInverseAtIP(self, pp, xcoor):
return self.space.GetJacobianDeterminantAndInverse(self.valdphidxi[pp], xcoor)
[docs] def GetJacobianDeterminantAndInverse(self, Nfder, xcoor):
return self.space.GetJacobianDeterminantAndInverse(Nfder, xcoor)
[docs] def GetJackAndDet(self, Nfder, xcoor):
raise Exception('Deprecated Please use GetJacobianDeterminantAndInverse ')
[docs] def GetNormal(self, Jacobian):
return self.space.GetNormal(Jacobian)
[docs]def CheckIntegrity(GUI=False):
SpaceBase()
SpaceAtIntegrationPoints()
return "OK"
if __name__ == '__main__':
print(CheckIntegrity()) # pragma: no cover