# -*- 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 sympy.matrices import Matrix
import Muscat.Containers.ElementsDescription as ED
from Muscat.Types import MuscatFloat
from Muscat.FE.Spaces.SymSpace import SymSpaceBase
[docs]class PyrSpaceBase(SymSpaceBase):
__oneOverSqrt = np.sqrt(3)
def __init__(self):
super(PyrSpaceBase,self).__init__()
self.dimensionality = 3
self.geoSupport = ED.GeoSupport.GeoPyr
[docs] def ClampParamCoordinates(self,xietaphi):
res = np.core.umath.maximum(xietaphi, 0)
#
# _._
# __-_/\_ -___
# / / \ /
# / / \ /
# / / \ /
# //_________\/
#"
#diag = np.array([0.,0.,1.]) - np.array([1.,1.,0.])
#diag /= np.linalg.norm(diag)
#xv1 = np.array([1.,0.,1.])# - np.array([1.,1.,0.])
#xv1 /= np.linalg.norm(xv1)
#perPlanxz = np.cross(diag,xv1)
perPlanXZ = np.array([-0.40824829, 0.81649658, 0.40824829])
#yv1 = np.array([0.,1.,1.])# - np.array([1.,1.,0.])
#yv1 /= np.linalg.norm(yv1)
#perPlanyz = np.cross(yv1,diag)
perPlanYZ =np.array([ 0.81649658, -0.40824829, 0.40824829])
d1 = (res-[1,1,0]).dot(perPlanXZ)
d2 = (res-[1,1,0]).dot(perPlanYZ)
#d3 = np.sqrt(eta**2/3 - 2*eta*phi/3 + 2*eta*xi/3 + phi**2/3 - 2*phi*xi/3 + xi**2/3 + 2/3)
#if d1 > 0 and d2 > 0 :
# res[2] = 0
#if res.dot([us2,us2,0]) >= us2 :#and res[1]+res[2] <= 1:
if res[0]+res[2] >= 1 and d1 < 0 :
dif = res[0]-res[2]
res[0] = (1+dif)/2.
res[2] = 1.-res[0]
if res[1]+res[2] > 1 and d2 < 0 :
dif = res[1]-res[2]
res[1] = (1+dif)/2.
res[2] = 1.-res[1]
if d1 >=0 and d2 >= 0 :
res[:] = res[1]/3 + res[0]/3 - res[2]/3 + 1/3, res[1]/3 + res[0]/3 - res[2]/3 + 1/3, -res[1]/3 + res[2]/3 - res[0]/3 + 2/3
return np.core.umath.maximum(np.core.umath.minimum(res, 1), 0)
[docs]class Pyr_P0_Global(PyrSpaceBase):
def __init__(self):
super(Pyr_P0_Global,self).__init__()
self.symN = Matrix([1])
self.posN = np.array([[None,None,None]])
self.dofAttachments = [("G",None,None)]
[docs]class Pyr_P0_Lagrange(PyrSpaceBase):
def __init__(self):
super(Pyr_P0_Lagrange,self).__init__()
self.symN = Matrix([1])
self.posN = np.array([[3/8,3/8,1/4]])
self.dofAttachments = [("C",0,None)]
[docs]class Pyr_P1_Lagrange(PyrSpaceBase):
def __init__(self):
super(Pyr_P1_Lagrange,self).__init__()
xi = self.xi
eta = self.eta
phi = self.phi
omz_inv = 1. / (1.-phi)
self.symN = Matrix([(1-xi-phi)*(1-eta-phi)*omz_inv,
(xi)*(1-eta-phi)*omz_inv,
(xi)*(eta)*omz_inv,
(1-xi-phi)*(eta)*omz_inv,
phi])
self.posN = np.array([[0,0,0],
[1,0,0],
[1,1,0],
[0,1,0],
[0,0,1],
])
self.dofAttachments = [("P",0,None),
("P",1,None),
("P",2,None),
("P",3,None),
("P",4,None) ]
self.NExceptions = [(1-phi, np.array([0,0,0,0,1], dtype=MuscatFloat)) ]
self.dNExceptions = [(1-phi, np.array([[ -1.0, -1.0, -1.0],
[ 1.0, 0.0, 0.0],
[ 0.0, 0.0, 0.0],
[ 0.0, 1.0, 0.0],
[ 0.0, 0.0, 1.0]],dtype=MuscatFloat).T) ]
[docs]class Pyr19_P2_Lagrange(PyrSpaceBase):
def __init__(self):
super(Pyr19_P2_Lagrange,self).__init__()
# article of R.S. Browning, K.T. Danielson and D.L. Littlefield Computer Methods in Applied Mechanics and Engineering 405 (2023)
# change of variable : translation and rotation of the pyramid to be between 0 and 1
eta1 = -2*self.xi-self.phi+1
eta2 = -2*self.eta-self.phi+1
eta3 = self.phi
omz_inv = 1./(eta3-1.)
eta1p = eta1+eta3-1
eta1m = eta1-eta3+1
eta2p = eta2+eta3-1
eta2m = eta2-eta3+1
# shape functions in Table 3
phi19 = -eta3*eta1m*eta1p*eta2m*eta2p*omz_inv**3*16/3 # volume
phi18 = -eta1*eta3*eta1m*eta2m*eta2p*omz_inv**2*27/8 # triangles
phi17 = -eta2*eta3*eta2p*eta1m*eta1p*omz_inv**2*27/8 # triangles
phi16 = -eta1*eta3*eta1p*eta2m*eta2p*omz_inv**2*27/8 # triangles
phi15 = -eta2*eta3*eta2m*eta1m*eta1p*omz_inv**2*27/8 # triangles
phi14tilde = eta1m*eta1p*eta2m*eta2p*omz_inv**2
phi14 = phi14tilde-phi19*9/16 # quadrangle
phi13 = eta3*eta1m*eta2p*omz_inv-(phi17+phi18)*4/9-phi19*3/16 # vertical edges
phi12 = -eta3*eta1p*eta2p*omz_inv-(phi16+phi17)*4/9-phi19*3/16 # vertical edges
phi11 = eta3*eta1p*eta2m*omz_inv-(phi15+phi16)*4/9-phi19*3/16 # vertical edges
phi10 = -eta3*eta1m*eta2m*omz_inv-(phi15+phi18)*4/9-phi19*3/16 # vertical edges
phi9 = eta1m*eta2m*eta2p*omz_inv/2-phi14tilde/2-phi18*4/9 # base edges # error in the article, it is phi18 and not phi15
phi8 = -eta1m*eta1p*eta2p*omz_inv/2-phi14tilde/2-phi17*4/9 # base edges # error in the article, it is phi17 and not phi15
phi7 = -eta1p*eta2m*eta2p*omz_inv/2-phi14tilde/2-phi16*4/9 # base edges # error in the article, it is phi16 and not phi15
phi6 = eta1m*eta1p*eta2m*omz_inv/2-phi14tilde/2-phi15*4/9 # base edges
phi5 = eta3*(2*eta3-1)+(phi15+phi16+phi17+phi18)/9+phi19/8 # vertices
phi4 = (eta1-eta2-1)*eta1m*eta2p*omz_inv/4+phi14tilde/4+(phi17+phi18)/9+phi19*3/64 # vertices # error in the article, it is 3/64 and not 1/64
phi3 = (eta1+eta2+1)*eta1p*eta2p*omz_inv/4+phi14tilde/4+(phi16+phi17)/9+phi19*3/64 # vertices # error in the article, it is 3/64 and not 1/64
phi2 = -(eta1-eta2+1)*eta1p*eta2m*omz_inv/4+phi14tilde/4+(phi15+phi16)/9+phi19*3/64 # vertices
phi1 = -(eta1+eta2-1)*eta1m*eta2m*omz_inv/4+phi14tilde/4+(phi15+phi18)/9+phi19*3/64 # vertices
self.symN = Matrix([phi1, phi2, phi3, phi4, phi5, # vertices
phi6, phi7, phi8, phi9, # base edges
phi10, phi11, phi12, phi13, # vertical edges
phi14, # quadrangle
phi15, phi16, phi17, phi18, # triangles
phi19 # volume
])
self.posN = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], # vertices
[0.5,0,0], [1,0.5,0], [0.5,1,0], [0,0.5,0], # base edges
[0,0,0.5], [0.5,0,0.5], [0.5,0.5,0.5], [0,0.5,0.5], # vertical edges
[0.5,0.5,0], # quadrangle
[1/3,0,1/3], [2/3,1/3,1/3], [1/3,2/3,1/3], [0,1/3,1/3], # triangles
[3/8,3/8,1/4] # volume
])
self.dofAttachments = [("P",0,None), ("P",1,None), ("P",2,None), ("P",3,None), ("P",4,None),
("F2",0,None), ("F2",1,None), ("F2",2,None), ("F2",3,None),
("F2",4,None), ("F2",5,None), ("F2",6,None), ("F2",7,None),
("F",0,None),
("F",1,None), ("F",2,None), ("F",3,None), ("F",4,None),
("C",0,None)
]
self.NExceptions = [(1-eta3, np.array([0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],dtype=MuscatFloat)) ]
self.dNExceptions = [(1-eta3, np.array([[ 1 , 1 , 1 ], [ -1 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , -1 , 0 ], [ 0 , 0 , 3 ],
[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ],
[ -4 , -4 , -4 ], [ 4 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 4 , 0 ],
[ 0 , 0 , 0 ],
[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ],
[ 0 , 0 , 0 ]],dtype=MuscatFloat).T) ]
[docs]class Pyr14_P2_Lagrange(PyrSpaceBase):
def __init__(self):
super(Pyr14_P2_Lagrange,self).__init__()
# change of variable : translation and rotation of the pyramid to be between 0 and 1
eta1 = -2*self.xi-self.phi+1
eta2 = -2*self.eta-self.phi+1
eta3 = self.phi
omz_inv = 1./(eta3-1.)
eta1p = eta1+eta3-1
eta1m = eta1-eta3+1
eta2p = eta2+eta3-1
eta2m = eta2-eta3+1
# shape functions
phi14 = eta1m*eta1p*eta2m*eta2p*omz_inv**2 # quadrangle
phi13 = eta3*eta1m*eta2p*omz_inv # vertical edges
phi12 = -eta3*eta1p*eta2p*omz_inv # vertical edges
phi11 = eta3*eta1p*eta2m*omz_inv # vertical edges
phi10 = -eta3*eta1m*eta2m*omz_inv # vertical edges
phi9 = eta1m*eta2m*eta2p*omz_inv/2-phi14/2 # base edges
phi8 = -eta1m*eta1p*eta2p*omz_inv/2-phi14/2 # base edges
phi7 = -eta1p*eta2m*eta2p*omz_inv/2-phi14/2 # base edges
phi6 = eta1m*eta1p*eta2m*omz_inv/2-phi14/2 # base edges
phi5 = eta3*(2*eta3-1) # vertices
phi4 = (eta1-eta2-1)*eta1m*eta2p*omz_inv/4+phi14/4 # vertices
phi3 = (eta1+eta2+1)*eta1p*eta2p*omz_inv/4+phi14/4 # vertices
phi2 = -(eta1-eta2+1)*eta1p*eta2m*omz_inv/4+phi14/4 # vertices
phi1 = -(eta1+eta2-1)*eta1m*eta2m*omz_inv/4+phi14/4 # vertices
self.symN = Matrix([phi1, phi2, phi3, phi4, phi5, # vertices
phi6, phi7, phi8, phi9, # base edges
phi10, phi11, phi12, phi13, # vertical edges
phi14, # quadrangle
])
self.posN = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], # vertices
[0.5,0,0], [1,0.5,0], [0.5,1,0], [0,0.5,0], # base edges
[0,0,0.5], [0.5,0,0.5], [0.5,0.5,0.5], [0,0.5,0.5], # vertical edges
[0.5,0.5,0], # quadrangle
])
self.dofAttachments = [("P",0,None), ("P",1,None), ("P",2,None), ("P",3,None), ("P",4,None),
("F2",0,None), ("F2",1,None), ("F2",2,None), ("F2",3,None),
("F2",4,None), ("F2",5,None), ("F2",6,None), ("F2",7,None),
("F",0,None),
]
self.NExceptions = [(1-eta3, np.array([0,0,0,0,1,0,0,0,0,0,0,0,0,0],dtype=MuscatFloat)) ]
self.dNExceptions = [(1-eta3, np.array([[ 1 , 1 , 1 ], [ -1 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , -1 , 0 ], [ 0 , 0 , 3 ],
[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ],
[ -4 , -4 , -4 ], [ 4 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 4 , 0 ],
[ 0 , 0 , 0 ]],dtype=MuscatFloat).T) ]
[docs]class Pyr13_P2_Lagrange(PyrSpaceBase):
def __init__(self):
super(Pyr13_P2_Lagrange,self).__init__()
# change of variable : translation and rotation of the pyramid to be between 0 and 1
eta1 = -2*self.xi-self.phi+1
eta2 = -2*self.eta-self.phi+1
eta3 = self.phi
omz_inv = 1./(eta3-1.)
eta1p = eta1+eta3-1
eta1m = eta1-eta3+1
eta2p = eta2+eta3-1
eta2m = eta2-eta3+1
# shape functions
phi13 = eta3*eta1m*eta2p*omz_inv # vertical edges
phi12 = -eta3*eta1p*eta2p*omz_inv # vertical edges
phi11 = eta3*eta1p*eta2m*omz_inv # vertical edges
phi10 = -eta3*eta1m*eta2m*omz_inv # vertical edges
phi9 = eta1m*eta2m*eta2p*omz_inv/2 # base edges
phi8 = -eta1m*eta1p*eta2p*omz_inv/2 # base edges
phi7 = -eta1p*eta2m*eta2p*omz_inv/2 # base edges
phi6 = eta1m*eta1p*eta2m*omz_inv/2 # base edges
phi5 = eta3*(2*eta3-1) # vertices
phi4 = (eta1-eta2-1)*eta1m*eta2p*omz_inv/4 # vertices
phi3 = (eta1+eta2+1)*eta1p*eta2p*omz_inv/4 # vertices
phi2 = -(eta1-eta2+1)*eta1p*eta2m*omz_inv/4 # vertices
phi1 = -(eta1+eta2-1)*eta1m*eta2m*omz_inv/4 # vertices
self.symN = Matrix([phi1, phi2, phi3, phi4, phi5, # vertices
phi6, phi7, phi8, phi9, # base edges
phi10, phi11, phi12, phi13, # vertical edges
])
self.posN = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], # vertices
[0.5,0,0], [1,0.5,0], [0.5,1,0], [0,0.5,0], # base edges
[0,0,0.5], [0.5,0,0.5], [0.5,0.5,0.5], [0,0.5,0.5], # vertical edges
])
self.dofAttachments = [("P",0,None), ("P",1,None), ("P",2,None), ("P",3,None), ("P",4,None),
("F2",0,None), ("F2",1,None), ("F2",2,None), ("F2",3,None),
("F2",4,None), ("F2",5,None), ("F2",6,None), ("F2",7,None),
]
self.NExceptions = [(1-eta3, np.array([0,0,0,0,1,0,0,0,0,0,0,0,0],dtype=MuscatFloat)) ]
self.dNExceptions = [(1-eta3, np.array([[ 1 , 1 , 1 ], [ -1 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , -1 , 0 ], [ 0 , 0 , 3 ],
[ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 0 , 0 ],
[ -4 , -4 , -4 ], [ 4 , 0 , 0 ], [ 0 , 0 , 0 ], [ 0 , 4 , 0 ],
],dtype=MuscatFloat).T) ]
[docs]def CheckIntegrity(GUI=False):
from Muscat.Containers.MeshCreationTools import CreateCube
n = 5
n2 = n-1
mesh = CreateCube(dimensions=[n,n,n],origin=[0,0,0],spacing=[1/(n2),1/(n2),1/(n2)])
points = mesh.nodes
mypyr=Pyr13_P2_Lagrange()
mypyr.Create()
cp = np.empty_like(points)
sfc = np.empty((mesh.GetNumberOfNodes(),mypyr.GetNumberOfShapeFunctions()))
sf = np.empty((mesh.GetNumberOfNodes(),mypyr.GetNumberOfShapeFunctions()))
for cpt, p in enumerate(points):
cp[cpt,:] = mypyr.ClampParamCoordinates(p)
sfc[cpt,:] = mypyr.GetShapeFunc(cp[cpt,:])
sf[cpt,:] = mypyr.GetShapeFunc(p)
mesh.nodeFields["cp"] = cp
for cpt in range(mypyr.GetNumberOfShapeFunctions()):
mesh.nodeFields[f"sfc{cpt}"] = sfc[:,cpt]
mesh.nodeFields[f"sf{cpt}"] = sf[:,cpt]
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(mesh,run=GUI)
print(mypyr.GetShapeFunc([0,0,0.5]) )
print(mypyr.GetShapeFunc([0.5,0,0.5]))
print(mypyr.GetShapeFunc([0.5,0.5,0.5]))
print(mypyr.GetShapeFunc([0,0.5,0.5]))
print(mypyr.GetShapeFunc([0,0,1]))
print("---------------------------")
print(mypyr.GetShapeFuncDer([0,0,1]))
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(True))# pragma: no cover