# -*- 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 __future__ import annotations
from typing import Any, NewType, Dict, Tuple, List, Optional, Union, Generator
import numpy as np
from Muscat.Types import MuscatFloat, ArrayLike
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.Helpers.Factory import Factory
import Muscat.FE.Spaces.QuadSpaces as QS
import Muscat.FE.Spaces.BarSpaces as BaS
import Muscat.FE.Spaces.WedgeSpaces as WS
import Muscat.FE.Spaces.HexaSpaces as HS
import Muscat.FE.Spaces.TetSpaces as TS
import Muscat.FE.Spaces.TriSpaces as TrS
import Muscat.FE.Spaces.PyrSpaces as PyS
from Muscat.LinAlg.MatVecOperations import ImmutableView
[docs]
class ElementaryQuadrature:
__slots__ = ("elementType", "weights", "points")
def __init__(self, elementType: ED.ElementType, points: ArrayLike, weights: ArrayLike):
self.weights: np.ndarray = ImmutableView(np.asarray(weights, dtype=MuscatFloat))
points = np.asarray(points, dtype=MuscatFloat, order='C')
if points.shape[1] < 3:
points = np.hstack((points, np.zeros((points.shape[0], 3 - points.shape[1]), dtype=MuscatFloat)))
elif points.shape[1] > 3: # pragma: no cover
raise RuntimeError(f"integration points must have at most 3 coordinates {points.shape}")
self.points: np.ndarray = ImmutableView(np.asarray(points, dtype=MuscatFloat))
self.elementType: ED.ElementType = ED.ElementType(elementType)
if len(self.weights) != len(self.points): # pragma: no cover
raise RuntimeError(f"Incompatible integration rule for elements type {elementType}: ({self.weights},{self.points})")
[docs]
def GetNumberOfPoints(self):
return len(self.weights)
def __str__(self):
return f"Integration Rule for {self.elementType}: points {self.points}, weights {self.weights}"
[docs]
def CopyFor(self, elementType: ED.ElementType) -> ElementaryQuadrature:
return ElementaryQuadrature(elementType, points=self.points, weights=self.weights)
[docs]
class MeshQuadrature:
__slots__ = ("name", "elementQuadrature", "__weakref__")
def __init__(self, name: str):
self.name: str = str(name)
self.elementQuadrature: Dict[ED.ElementType, ElementaryQuadrature] = dict()
[docs]
def AddElementaryQuadrature(self, elementaryQuadrature: ElementaryQuadrature):
if elementaryQuadrature.elementType in self.elementQuadrature: #pragma: no cover
raise RuntimeError(f"Quadrature for this ({elementaryQuadrature.elementType}) element already set)")
self.elementQuadrature[elementaryQuadrature.elementType] = elementaryQuadrature
def __getitem__(self, elementType: ED.ElementType) -> ElementaryQuadrature:
return self.elementQuadrature[elementType]
[docs]
def keys(self) -> ED.ElementType:
return self.elementQuadrature.keys()
[docs]
def values(self) -> Generator[ElementaryQuadrature]:
return self.elementQuadrature.values()
def __call__(self) -> MeshQuadrature:
return self
[docs]
def items(self) -> Generator[ED.ElementType, ElementaryQuadrature]:
return self.elementQuadrature.items()
def __str__(self):
res = f"MeshQuadrature(name='{self.name}')\n"
for v in self.elementQuadrature.values():
res += " " + str(v)
return res
[docs]
class IntegrationRulesFactory(Factory):
_Catalog = {}
_SetCatalog = set()
def __init__(self):
super().__init__()
[docs]
def RegisterIntegrationRule(ir: MeshQuadrature):
return IntegrationRulesFactory.RegisterClass(ir.name, ir, withError=True)
[docs]
def GetIntegrationRuleByName(ruleName: str) -> MeshQuadrature:
return IntegrationRulesFactory.Create(ruleName)
[docs]
def GetIntegrationRule(ruleName: Optional[str] = None, rule: Optional[MeshQuadrature] = None) -> MeshQuadrature:
"""Return the integration rule for the mesh using the name or an integration rule. In the case where
both are None then LagrangeIsoParamQuadrature is returned
Parameters
----------
ruleName : str, optional
name of the integration rule, by default None
rule : MeshQuadrature, optional
integration rule in the case ruleName is None, by default None
Returns
-------
MeshQuadrature
the integration rule
"""
if ruleName is None:
if rule is None:
return GetIntegrationRuleByName("LagrangeIsoParamQuadrature")
else:
return rule
else:
if rule is None:
return GetIntegrationRuleByName(ruleName)
else:
if ruleName != rule.name:
raise RuntimeError("must give ruleName or rule not both")
return rule
[docs]
def TensorProductGauss(dim: int, nbPoints: int = 2) -> Tuple[np.ndarray, np.ndarray]:
"""Tensor product of gauss point in dimension (dim) using nbPoints (number of points) per dimension
Parameters
----------
dim : int
the dimensionality of the tensor product
nbPoints : int, optional
Number of points per dimension , by default 2
Returns
-------
IntegrationRuleType
Integration Rule (points and weights)
"""
import math
if nbPoints == 2:
p = [-math.sqrt(1.0 / 3) / 2 + 0.5, math.sqrt(1.0 / 3) / 2 + 0.5]
w = [1.0 / 2.0, 1.0 / 2.0]
elif nbPoints == 3:
# https://fr.wikipedia.org/wiki/M%C3%A9thodes_de_quadrature_de_Gauss
p = [-math.sqrt(3.0 / 5.0) / 2 + 0.5, 0.5, math.sqrt(3.0 / 5.0) / 2 + 0.5]
w = [5.0 / 18.0, 8.0 / 18.0, 5.0 / 18.0]
else:
raise # pragma: no cover
if dim == 1:
return (
np.array(
[
p,
]
).T,
np.array(w),
)
return TensorProdHomogeneous(dim, np.array([p]).T, np.array(w))
# Utilities to construct 2D and 3D quadratures
[docs]
def TensorProdHomogeneous(dim: int, p: ArrayLike, w: ArrayLike):
if dim == 2:
return TensorProd(p, w, p, w)
elif dim == 3:
return TensorProd(p, w, p, w, p, w)
else:
raise # pragma: no cover
[docs]
def TensorProd(p1: ArrayLike, w1: ArrayLike, p2: ArrayLike, w2: ArrayLike, p3: Optional[ArrayLike] = None, w3: Optional[ArrayLike] = None):
PRes = []
WRes = []
p1 = np.asarray(p1, dtype=MuscatFloat)
w1 = np.asarray(w1, dtype=MuscatFloat)
p2 = np.asarray(p2, dtype=MuscatFloat)
w2 = np.asarray(w2, dtype=MuscatFloat)
if p3 is None:
for i in range(len(p1)):
for j in range(len(p2)):
res: List[MuscatFloat] = []
res.extend(p1[i, :])
res.extend(p2[j, :])
PRes.append(res)
WRes.append(w1[i] * w2[j])
return (np.array(PRes), np.array(WRes))
else:
p3 = np.asarray(p3, dtype=MuscatFloat)
w3 = np.asarray(w3, dtype=MuscatFloat)
for k in range(len(p3)):
for j in range(len(p2)):
for i in range(len(p1)):
res = []
res.extend(p1[i, :])
res.extend(p2[j, :])
res.extend(p3[k, :])
PRes.append(res)
# Pres.append([p1[i,0], p2[j,0]])
WRes.append(w1[i] * w2[j] * w3[k])
return (np.array(PRes), np.array(WRes))
# p23,w23 = TensorProd(p2,w2,p3,w3)
# return TensorProd(p1,w1,p23,w23)
##########################################################################################################################
elementCenterQuadrature = MeshQuadrature("ElementCenter")
# 0D elements
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0.0]], weights=[1.0]))
# 1D elements
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=[[1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=[[1.0 / 2]], weights=[1.0]))
# 2D elements
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=[[1.0 / 3, 1.0 / 3]], weights=[0.5]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=[[1.0 / 3, 1.0 / 3]], weights=[0.5]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=[[1.0 / 2, 1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=[[1.0 / 2, 1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=[[1.0 / 2, 1.0 / 2]], weights=[1.0]))
# 3D elements
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_4, points=[[1.0 / 4, 1.0 / 4, 1.0 / 4]], weights=[1.0 / 6]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_10, points=[[1.0 / 4, 1.0 / 4, 1.0 / 4]], weights=[1.0 / 6]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_8, points=[[1.0 / 2, 1.0 / 2, 1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_20, points=[[1.0 / 2, 1.0 / 2, 1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=[[1.0 / 2, 1.0 / 2, 1.0 / 2]], weights=[1.0]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=[[1.0 / 3, 1.0 / 3, 1.0 / 2]], weights=[0.5]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_15, points=[[1.0 / 3, 1.0 / 3, 1.0 / 2]], weights=[0.5]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_18, points=[[1.0 / 3, 1.0 / 3, 1.0 / 2]], weights=[0.5]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_5, points=[[3.0 / 8, 3 / 8, 1.0 / 4]], weights=[1.0 / 3]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_13, points=[[3.0 / 8, 3 / 8, 1.0 / 4]], weights=[1.0 / 3]))
elementCenterQuadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_14, points=[[3.0 / 8, 3 / 8, 1.0 / 4]], weights=[1.0 / 3]))
RegisterIntegrationRule(elementCenterQuadrature)
################# Lagrange P1 ##################################
LagrangeP1Quadrature = MeshQuadrature("LagrangeP1Quadrature")
# 0D elements
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0]], weights=[1.0]))
# 1D elements
__p, __w = TensorProductGauss(dim=1, nbPoints=2)
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=__p, weights=__w))
# 2D elements
__p = 1.0 / 6.0 * np.array([[1.0, 1.0], [4.0, 1.0], [1.0, 4.0]])
__w = 1.0 / 6.0 * np.array([1.0, 1.0, 1])
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=__p, weights=__w))
__p, __w = TensorProductGauss(dim=2, nbPoints=2)
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=__p, weights=__w))
# 3D elements
__p = np.array([[1.0 / 4, 1.0 / 4, 1.0 / 4]], dtype=MuscatFloat)
__w = np.array([1.0 / 6])
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_4, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_10, points=__p, weights=__w))
__p, __w = TensorProductGauss(dim=3, nbPoints=2)
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_8, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_20, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=__p, weights=__w))
__triLP1: ElementaryQuadrature = LagrangeP1Quadrature[ED.Triangle_3]
__barLP1: ElementaryQuadrature = LagrangeP1Quadrature[ED.Bar_2]
__p, __w = TensorProd(__triLP1.points[:, 0:2], __triLP1.weights, __barLP1.points[:, 0:1], __barLP1.weights)
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_15, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_18, points=__p, weights=__w))
__p = np.array([[3.0 / 8, 3.0 / 8, 1.0 / 4]], dtype=MuscatFloat)
__w = np.array([1.0 / 3])
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_5, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_13, points=__p, weights=__w))
LagrangeP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_14, points=__p, weights=__w))
RegisterIntegrationRule(LagrangeP1Quadrature)
################# Lagrange P2 ##################################
LagrangeP2Quadrature = MeshQuadrature("LagrangeP2Quadrature")
# 0D elements
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0]], weights=[1.0]))
# 1D elements
__p, __w = TensorProductGauss(dim=1, nbPoints=3)
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=__p, weights=__w))
# 2D elements
__p = np.array(
[
[0.445948490915, 0.445948490915, 0],
[0.108103018168, 0.445948490915, 0],
[0.445948490915, 0.108103018168, 0],
[0.091576213509, 0.091576213509, 0],
[0.816847572980, 0.091576213509, 0],
[0.091576213509, 0.816847572980, 0],
],
dtype=MuscatFloat,
)
__w = np.array([0.1116907948390, 0.1116907948390, 0.1116907948390, 0.0549758718276, 0.0549758718276, 0.0549758718276], dtype=MuscatFloat)
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=__p, weights=__w))
__p, __w = TensorProductGauss(dim=2, nbPoints=3)
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=__p, weights=__w))
# 3D elements
# from https://www.code-aster.org/V2/doc/v13/en/man_r/r3/r3.01.01.pdf
p0, p1, p2 = (1.0 / 4.0, 1.0 / 6.0, 1.0 / 2.0)
__p = np.array([[p1, p1, p1], [p2, p1, p1], [p1, p2, p1], [p1, p1, p2], [p0, p0, p0]])
w0, w1 = (-2.0 / 15, 3.0 / 40)
__w = np.array([w1, w1, w1, w1, w0])
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_4, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_10, points=__p, weights=__w))
__p, __w = TensorProductGauss(dim=3, nbPoints=3)
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_8, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_20, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=__p, weights=__w))
__triLP2: ElementaryQuadrature = LagrangeP2Quadrature[ED.Triangle_3]
__barLP2: ElementaryQuadrature = LagrangeP2Quadrature[ED.Bar_2]
__p, __w = TensorProd(__triLP2.points[:, 0:2], __triLP2.weights, __barLP2.points[:, 0:1], __barLP2.weights)
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_15, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_18, points=__p, weights=__w))
# article of R.S. Browning, K.T. Danielson and D.L. Littlefield Computer Methods in Applied Mechanics and Engineering 405 (2023)
# 5 points quadrature
p0, p1, p2 = (0.4868644955601477, 0.1666666666666667, 0.7000000000000000)
p3 = np.array([[-p0, -p0, p1], [p0, -p0, p1], [p0, p0, p1], [-p0, p0, p1], [0, 0, p2]])
p4 = np.array([[-1.0 / 2, 0, -1.0 / 2], [-1.0 / 2, 0, -1.0 / 2], [0, 0, 1]]) @ p3.T + np.array([[1.0 / 2], [1.0 / 2], [0]]) # change of coordinates to have a reference pyrmamid between 0 and 1
__p = p4.T
w0, w1 = (0.2812500000000000, 0.208333333333333)
w = np.array([w0, w0, w0, w0, w1])
__w = w / 4 # renormalisation of the volume
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_5, points=__p, weights=__w))
# article of R.S. Browning, K.T. Danielson and D.L. Littlefield Computer Methods in Applied Mechanics and Engineering 405 (2023)
# 13 points quadrature
p0, p1, p2, p3, p4, p5, p6 = (0.3851039921187038, 0.4034583196072821, 0.5315787743696198, 0.4285714285714285, 0.3392857142857143, 0.0849673202614379, 0.7621970180376850)
p9 = np.array(
[[-p0, -p0, p3], [p0, -p0, p3], [p0, p0, p3], [-p0, p0, p3], [-p1, 0, p4], [p1, 0, p4], [0, -p1, p4], [0, p1, p4], [-p2, -p2, p5], [p2, -p2, p5], [p2, p2, p5], [-p2, p2, p5], [0, 0, p6]]
)
p10 = np.array([[-1.0 / 2, 0, -1.0 / 2], [-1.0 / 2, 0, -1.0 / 2], [0, 0, 1]]) @ p9.T + np.array([[1.0 / 2], [1.0 / 2], [0]]) # change of coordinates to have a reference pyrmamid between 0 and 1
__p = p10.T
w0, w1, w2, w3 = (0.0840821256038648, 0.0561359290874341, 0.1756270761022237, 0.069952810159243)
w4 = np.array([w0, w0, w0, w0, w1, w1, w1, w1, w2, w2, w2, w2, w3])
__w = w4 / 4 # renormalisation of the volume
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_13, points=__p, weights=__w))
LagrangeP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_14, points=__p, weights=__w))
RegisterIntegrationRule(LagrangeP2Quadrature)
################# Lagrange LagrangeIsoParamQuadrature ##################################
LagrangeIsoParamQuadrature = MeshQuadrature("LagrangeIsoParamQuadrature")
for name in LagrangeP1Quadrature.keys():
if ED.degree[name] == 1:
LagrangeIsoParamQuadrature.AddElementaryQuadrature(LagrangeP1Quadrature[name])
else:
LagrangeIsoParamQuadrature.AddElementaryQuadrature(LagrangeP2Quadrature[name])
RegisterIntegrationRule(LagrangeIsoParamQuadrature)
##### Nodal P1 Integration points for the evaluation of post quantities at nodes ######
NodalEvaluationP1Quadrature = MeshQuadrature("NodalEvaluationP1")
# 0D elements
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0.0]], weights=[1.0]))
# 1D elements
__p = BaS.Bar_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=__p, weights=np.ones(__p.shape[0])))
# 2D elements
__p = TrS.Tri_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=__p, weights=np.ones(__p.shape[0])))
__p = QS.Quad_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=__p, weights=np.ones(__p.shape[0])))
# 3D elements
__p = TS.Tet_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_4, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_10, points=__p, weights=np.ones(__p.shape[0])))
__p = HS.Hexa_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_8, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_20, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=__p, weights=np.ones(__p.shape[0])))
__p = WS.Wedge_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_15, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_18, points=__p, weights=np.ones(__p.shape[0])))
__p = PyS.Pyr_P1_Lagrange().posN
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_5, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_13, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP1Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_14, points=__p, weights=np.ones(__p.shape[0])))
RegisterIntegrationRule(NodalEvaluationP1Quadrature)
##### Nodal P2 Integration points for the evaluation of post quantities at nodes ######
NodalEvaluationP2Quadrature = MeshQuadrature("NodalEvaluationP2")
# 0D elements
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0.0]], weights=[1.0]))
# 1D elements
__p = BaS.Bar_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=__p, weights=np.ones(__p.shape[0])))
# 2D elements
__p = TrS.Tri_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=__p, weights=np.ones(__p.shape[0])))
__p = QS.Quad_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=__p, weights=np.ones(__p.shape[0])))
# 3D elements
__p = TS.Tet_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_4, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Tetrahedron_10, points=__p, weights=np.ones(__p.shape[0])))
__p = HS.Hexa_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_8, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_20, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=__p, weights=np.ones(__p.shape[0])))
__p = WS.Wedge_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_15, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_18, points=__p, weights=np.ones(__p.shape[0])))
__p = PyS.Pyr13_P2_Lagrange().posN
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_5, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_13, points=__p, weights=np.ones(__p.shape[0])))
NodalEvaluationP2Quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Pyramid_14, points=__p, weights=np.ones(__p.shape[0])))
RegisterIntegrationRule(NodalEvaluationP2Quadrature)
##### Nodal IsoParametric Integration points for the evaluation of post quantities at nodes ######
NodalEvaluationIsoParamQuadrature: MeshQuadrature = MeshQuadrature("NodalEvaluationIsoParam")
for name in NodalEvaluationP1Quadrature.keys():
if ED.degree[name] == 1:
NodalEvaluationIsoParamQuadrature.AddElementaryQuadrature(NodalEvaluationP1Quadrature[name])
else:
NodalEvaluationIsoParamQuadrature.AddElementaryQuadrature(NodalEvaluationP2Quadrature[name])
RegisterIntegrationRule(NodalEvaluationIsoParamQuadrature)
################# trapezoidalOrderGenerated Quadratures ##################################
[docs]
def GenerateTrapezoidalOrderRuleOfOrder(i):
name = f"TrapezoidalOrder{i}"
ir = IntegrationRulesFactory._Catalog.get(name, (None, None))
if ir[0] is not None:
return ir[0]
__quadrature: MeshQuadrature = MeshQuadrature(name)
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Point_1, points=[[0]], weights=[1.0]))
if i == 1:
w = np.ones(1)
p = np.array([[0.5]])
elif i == 2:
w = np.ones(2) / 2
p = np.array([[0, 1]]).T
else:
w = np.full(i, 1 / (2 * i - 2))
w[1:-1] *= 2
p = np.arange(i)[:, np.newaxis] / (i - 1)
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_2, points=p, weights=w))
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Bar_3, points=p, weights=w))
__p, __w = TensorProdHomogeneous(2, p, w)
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_4, points=__p, weights=__w))
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_8, points=__p, weights=__w))
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Quadrangle_9, points=__p, weights=__w))
__p, __w = TensorProdHomogeneous(3, p, w)
ec = ElementaryQuadrature(elementType=ED.Hexahedron_8, points=__p, weights=__w)
__quadrature.AddElementaryQuadrature(ec)
__quadrature.AddElementaryQuadrature(ec.CopyFor(elementType=ED.Hexahedron_20))
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Hexahedron_27, points=__p, weights=__w))
if i == 1:
p = 1.0 / 2.0 * np.array([[1, 1.0]], dtype=MuscatFloat)
w = 1.0 / 2.0 * np.array([1.0], dtype=MuscatFloat)
else:
p = []
w = np.ones(i * (i + 1) // 2)
cpt = 0
for j in range(i):
phi0 = 1.0 / (i - 1) * j
for k in range(i - j):
phi1 = 1.0 / (i - 1) * (k)
p.append([phi0, phi1])
# halp of the contribution of the border points
if j == 0 or j == i - 1 or k == 0 or k == (i - j - 1):
w[cpt] *= 0.5
# 1/6 for the corner points
if (j == 0 and k == 0) or (j == 0 and k == i - j - 1) or (j == i - 1 and k == 0):
w[cpt] = 1.0 / 6.0
cpt += 1
w = np.array(w)
# normalization of the weight to sum = 1/2
w *= 1 / (np.sum(w) * 2)
p = np.array(p)
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_3, points=p, weights=w))
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Triangle_6, points=p, weights=w))
bar2rule = __quadrature[ED.Bar_2]
__p, __w = TensorProd(p[:, 0:2], w, bar2rule.points[:, 0:1], bar2rule.weights)
__quadrature.AddElementaryQuadrature(ElementaryQuadrature(elementType=ED.Wedge_6, points=__p, weights=__w))
IntegrationRulesFactory.RegisterClass(__quadrature.name, __quadrature, ir[1], withError=False)
return __quadrature
trapezoidalOrderGenerated = 6
# trapezoidal Rules
[docs]
def GetConstructorTrapezoidalOrderRuleOfOrder(i):
def WithOptions(opt=None):
return GenerateTrapezoidalOrderRuleOfOrder(i)
return WithOptions
for i in range(1, trapezoidalOrderGenerated):
IntegrationRulesFactory.RegisterClass(name=f"TrapezoidalOrder{i}", classType=None, constructor=GetConstructorTrapezoidalOrderRuleOfOrder(i), withError=True)
[docs]
def CheckIntegrity(GUI: bool = False):
IntegrationRulesFactory()
from Muscat.Helpers.CheckTools import MustFailFunctionWith
GetIntegrationRule("LagrangeIsoParamQuadrature")
MustFailFunctionWith(RuntimeError, GetIntegrationRule, ruleName="LagrangeP1Quadrature", rule=elementCenterQuadrature)
Volume: Dict[ED.ElementType, Union[MuscatFloat, float]] = {
ED.Point_1: MuscatFloat(1.0),
ED.Bar_2: MuscatFloat(1.0),
ED.Bar_3: MuscatFloat(1.0),
ED.Quadrangle_4: 1.0,
ED.Quadrangle_8: 1.0,
ED.Quadrangle_9: MuscatFloat(1.0),
ED.Hexahedron_8: 1.0,
ED.Hexahedron_20: 1.0,
ED.Hexahedron_27: 1.0,
ED.Triangle_3: 0.5,
ED.Triangle_6: 0.5,
ED.Wedge_6: 0.5,
ED.Wedge_15: 0.5,
ED.Wedge_18: 0.5,
ED.Tetrahedron_4: 1.0 / 6.0,
ED.Tetrahedron_10: 1.0 / 6.0,
ED.Pyramid_5: 1.0 / 3.0,
ED.Pyramid_13: 1.0 / 3.0,
ED.Pyramid_14: 1.0 / 3.0,
}
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.Helpers.TextFormatHelper import TFormat
for quadratureName in IntegrationRulesFactory.keys():
meshQuadrature = IntegrationRulesFactory.Create(quadratureName)
print(meshQuadrature)
GetIntegrationRule()
GetIntegrationRule(ruleName=quadratureName)
GetIntegrationRule(rule=meshQuadrature)
GetIntegrationRule(ruleName=quadratureName, rule=meshQuadrature)
meshQuadrature.items()
print(TFormat.Center(quadratureName, "*"))
for elementType, elemQuadrature in zip(meshQuadrature.keys(), meshQuadrature.values()):
lw = len(elemQuadrature.weights)
lp = len(elemQuadrature.points)
lc = elemQuadrature.GetNumberOfPoints()
print(f"{quadratureName} {elementType} has {lw} integration points")
print(elemQuadrature)
# print(str(elemQuadrature))
np.testing.assert_equal(lp, lw)
np.testing.assert_equal(lp, lc)
# Nodal... not designed for integration
if quadratureName in ["NodalEvaluationP1", "NodalEvaluationP2", "NodalEvaluationIsoParam"]:
continue
referenceVolume = Volume[elementType]
if np.abs(np.sum(elemQuadrature.weights) - referenceVolume) > 1e-10: # pragma: no cover
print(elemQuadrature)
print("Measure : ", np.sum(elemQuadrature.weights))
print("Measure Error: ", np.sum(elemQuadrature.weights) - referenceVolume)
raise RuntimeError(f"{quadratureName} {elementType} does not integrate constant function")
# End check of constant function ##################
# Checking linear function ########################
# can't integrate over a point
if elementType in [ED.Point_1]:
continue
# TrapezoidalP1 cant integrate linear functions
if quadratureName in ["TrapezoidalOrder1", "ElementCenter"]:
continue
def f1(x):
return x[0] - 0.5
if referenceVolume == 1.0:
intRef = 0.0
elif referenceVolume == 0.5:
intRef = (1 / 3 - 0.5) / 2
elif referenceVolume == 1 / 6:
intRef = -0.25 / 6
elif referenceVolume == 1 / 3:
intRef = -1.0 / 24
else: # pragma: no cover
raise
space_ipValues = LagrangeSpaceGeo[elementType].GetSpaceEvaluatedAt(elemQuadrature.points, elemQuadrature.weights)
integral = 0
for ip in range(len(elemQuadrature.weights)):
Jacobian, JDet, JInv = space_ipValues.GetJacobianDeterminantAndInverseAtIP(ip, LagrangeSpaceGeo[elementType].posN)
integral += JDet * elemQuadrature.weights[ip] * f1(elemQuadrature.points[ip])
# integral = data[1].dot([f1(x) for x in data[0]])
if np.abs(integral - intRef) > 1e-10: # pragma: no cover
print("function : f(x) = x-0.5")
print("int S f(x)dx = {}".format(integral))
print("Integral Exact: ", intRef)
print("Integral Error: ", np.abs(integral - intRef))
print(elemQuadrature)
raise RuntimeError(f"{quadratureName} {elementType} does not integrate f(x) = x-0.5")
# ElementCenterEval and LagrangeP1Quadrature cant integrate quadratic functions
if quadratureName in ["ElementCenter", "LagrangeP1Quadrature"] or quadratureName.startswith("Trapezoidal"):
continue
# this is a degree 1 integration
if quadratureName == "LagrangeIsoParamQuadrature" and elementType == ED.Tetrahedron_4:
continue
if quadratureName == "LagrangeIsoParamQuadrature" and elementType == ED.Pyramid_5:
continue
# Checking quadratic integration
def f2(x):
return (x[0] - 0.5) ** 2
if referenceVolume == 1.0:
intRef = 1.0 / 12.0
elif referenceVolume == 0.5:
intRef = 1 / 24.0
elif referenceVolume == 1 / 6:
intRef = 1.0 / 60.0
elif referenceVolume == 1 / 3:
intRef = 3.0 / 120
else: # pragma: no cover
raise
space_ipValues = LagrangeSpaceGeo[elementType].GetSpaceEvaluatedAt(elemQuadrature.points, elemQuadrature.weights)
integral = 0
for ip in range(len(elemQuadrature.weights)):
Jacobian, JDet, JInv = space_ipValues.GetJacobianDeterminantAndInverseAtIP(ip, LagrangeSpaceGeo[elementType].posN)
integral += JDet * elemQuadrature.weights[ip] * f2(elemQuadrature.points[ip])
# integral = data[1].dot([f1(x) for x in data[0]])
if np.abs(integral - intRef) > 1e-10: # pragma: no cover
print("function : f(x) = (x-0.5)**2")
print("int S f(x)dx = {}".format(integral))
print("Integral Exact: ", intRef)
print("Integral Error: ", np.abs(integral - intRef))
# print(data)
raise Exception(f"{quadratureName} {elementType} does not integrate f(x) = (x-0.5)**2")
return "ok"
if __name__ == "__main__":
print(CheckIntegrity(True)) # pragma: no cover