# -*- 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, List, Tuple
import numpy as np
import os
import tempfile
import platform
import getpass
import sympy
from sympy import Symbol, DiracDelta, Matrix, Expr
from sympy.utilities.lambdify import lambdify
from Muscat.Types import NDArray
from Muscat.Helpers.Cache import CachedResultDecorator
from Muscat.FE.Spaces.SpaceBase import SpaceBase, SpaceAtIntegrationPoints
[docs]
def DiracDeltaNumeric(data, der=None):
if abs(data) > 1e-15:
return 0
else:
return 1
cachePath = tempfile.gettempdir() + os.sep + f"{getpass.getuser()}_Muscat_SymSpaceBase_Create_{str(platform.python_version())}_{str(sympy.__version__)}safe_to_delete" + os.sep
@CachedResultDecorator(name="Create", path=cachePath, needDill=True)
def ProcessShapeFunctions(nbDim: int, symN: Matrix, NExceptions: List[Tuple[Expr, NDArray]], dNExceptions: List[Tuple[Expr, NDArray]]):
"""Function to process shape function and to compute derivatives and lambdas for the evaluation
Parameters
----------
nbDim : int
number of local coordinates of the shape functions
symN : Matrix
A sympy matrix containing the shape function for every dof
NExceptions : :List[Tuple[Expr, NDArray]]
the Value exceptions for the shape function (see PyrSpaces.py for a use of it)
dNExceptions : :List[Tuple[Expr, NDArray]]
the Value exceptions for the first derivative of shape function (see PyrSpaces.py for a use of it)
Returns
-------
lcoords : List[Symbols]
List of the coordinate of the elements
symN : Matrix
return of the argument symN
fct_N_Matrix : Callable
Function to evaluate the shape functions on a local position
symdNdxi : Matrix
Matrix of the derivatives of the shape functions
fct_dNdxi_Matrix : Callable
Function to evaluate the shape functions on a local position
symdNdxidxi : Matrix
Matrix of the second derivatives of the shape functions
fct_dNdxidxi_Matrix : Callable
Function to evaluate the shape functions on a local position
"""
allcoords = (Symbol("xi"), Symbol("eta"), Symbol("phi"))
lcoords = tuple(allcoords[x] for x in range(nbDim))
subsList = [(DiracDelta(0), 1.0), (DiracDelta(0, 1), 1.0), (DiracDelta(0, 2), 1.0)]
lambdifyList = [{"DiracDelta": DiracDeltaNumeric}, "numpy"]
nbSF = len(symN)
############# shape function treatment ########################
clean_N = symN.subs(subsList)
fct_N_Matrix = lambdify(allcoords, [clean_N[i] for i in range(nbSF)], lambdifyList)
if NExceptions is not None:
def GetModifiedShapeFunctions(_fct_N_Matrix, NExceptions):
def ShapeFunctions(xi_n, eta_n, phi_n):
for condition, values in NExceptions:
if abs(condition.subs({allcoords[0]: xi_n, allcoords[1]: eta_n, allcoords[2]: phi_n})) < 1.0e-6:
return values
else:
return _fct_N_Matrix(xi_n, eta_n, phi_n)
return ShapeFunctions
fct_N_Matrix = GetModifiedShapeFunctions(fct_N_Matrix, NExceptions)
############# shape functions first derivative #################
symdNdxi = [[None] * nbSF for i in range(nbDim)]
for i in range(nbDim):
for j in range(nbSF):
symdNdxi[i][j] = symN[j].diff(lcoords[i])
symdNdxi = Matrix(symdNdxi)
if symdNdxi.shape == (0, 0):
fct_dNdxi_Matrix = lambda xi, chi, phi: np.empty((0, 0))
else:
fct_dNdxi_Matrix = lambdify(allcoords, symdNdxi.subs(subsList), lambdifyList)
if dNExceptions is not None:
def GetModifiedShapeFunctionsDef(fct_dNdxi_Matrix, dNExceptions):
def ShapeFunctionsDer(xi_n, eta_n, phi_n):
for condition, values in dNExceptions:
if abs(condition.subs({allcoords[0]: xi_n, allcoords[1]: eta_n, allcoords[2]: phi_n})) < 1.0e-6:
return values
else:
return fct_dNdxi_Matrix(xi_n, eta_n, phi_n)
return ShapeFunctionsDer
fct_dNdxi_Matrix = GetModifiedShapeFunctionsDef(fct_dNdxi_Matrix, dNExceptions)
############ shape functions second derivative ################
symdNdxidxi = [None] * nbSF
fct_dNdxidxi_Matrix = [None] * nbSF
for i in range(nbSF):
symdNdxidxi[i] = [[0] * nbDim for j in range(nbDim)]
for j in range(nbDim):
for k in range(nbDim):
func = symN[i].diff(lcoords[j]).diff(lcoords[k])
symdNdxidxi[i][j][k] = func
symdNdxidxi[i] = Matrix(symdNdxidxi[i])
fct_dNdxidxi_Matrix[i] = lambdify(allcoords, symdNdxidxi[i].subs(subsList), lambdifyList)
return lcoords, symN, fct_N_Matrix, symdNdxi, fct_dNdxi_Matrix, symdNdxidxi, fct_dNdxidxi_Matrix
[docs]
class SymSpaceBase(SpaceBase):
def __init__(self):
super().__init__()
# ----symbolic part use Create() to pass to the next step ------------
self.xi = Symbol("xi")
self.eta = Symbol("eta")
self.phi = Symbol("phi")
self.symN = None
self.symdNdxi = None
self.__created__ = False
def __getstate__(self):
return None
def __setstat__(self, state):
self.__init__()
[docs]
def GetNumberOfShapeFunctions(self):
return len(self.symN)
[docs]
def Created(self):
return self.__created__
[docs]
def Create(self, force=False):
if self.__created__ and not force:
return
self.__created__ = True
nbDim = self.GetDimensionality()
self.lcoords, self.symN, self.fct_N_Matrix, self.symdNdxi, self.fct_dNdxi_Matrix, self.symdNdxidxi, self.fct_dNdxidxi_Matrix = ProcessShapeFunctions(
nbDim, self.symN, self.__dict__.get("NExceptions", None), self.__dict__.get("dNExceptions", None)
)
# This is not perfect because for P2 spaces we check only the number of shape functions
if np.all([dof[0] == "P" and nb == dof[1] for nb, dof in enumerate(self.dofAttachments)]):
self.fromConnectivityCompatible = True
[docs]
def GetSpaceEvaluatedAt(self, points, weights):
self.Create()
return SpaceAtIntegrationPoints(space=self, points=points, weights=weights)
[docs]
def GetPosOfShapeFunction(self, i, Xi):
valN = self.GetShapeFunc(self.posN[i, :])
return np.dot(valN, Xi).T
[docs]
def GetShapeFunc_default(self, xi=0, chi=0, phi=0):
return self.fct_N_Matrix(xi, chi, phi)
[docs]
def GetShapeFunc(self, qcoor):
return np.asarray(self.GetShapeFunc_default(*qcoor), dtype=float)
[docs]
def GetShapeFuncDer_default(self, xi=0, chi=0, phi=0):
return self.fct_dNdxi_Matrix(xi, chi, phi)
[docs]
def GetShapeFuncDer(self, qcoor):
return np.asarray(self.GetShapeFuncDer_default(*qcoor), dtype=float)
[docs]
def GetShapeFuncDerDer_default(self, xi=0, chi=0, phi=0):
return [np.asarray(x(xi, chi, phi), dtype=float) for x in self.fct_dNdxidxi_Matrix]
[docs]
def GetShapeFuncDerDer(self, qcoor):
return self.GetShapeFuncDerDer_default(*qcoor)
[docs]
def CheckIntegrity():
np.testing.assert_equal(DiracDeltaNumeric(0), 1.0)
np.testing.assert_equal(DiracDeltaNumeric(1), 0.0)
a = SymSpaceBase()
a.Created()
return "ok"