from typing import Dict, List, Tuple
import numpy as np
from numpy.typing import NDArray
from Muscat.Containers.ElementsContainers import AllElements, ElementsContainer
from Muscat.Containers.Mesh import Mesh
from Muscat.FE.Spaces.QuadSpaces import Quad8_P2_Lagrange
from Muscat.FE.Spaces.HexaSpaces import Hexa20_P2_Lagrange
import Muscat.Containers.ElementsDescription as ED
from Muscat.Types import MuscatIndex, MuscatFloat
from Muscat.Containers.MeshInspectionTools import GetVolumePerElement
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.FE.Fields.FieldTools import FillFEField
[docs]def Tetrahedrization(mesh: Mesh) -> Mesh:
"""Convert all elements in a mesh in simplex elements
Warning : originalIds of the new mesh are related to the global number of
the elements in the input mesh
Parameters
----------
mesh : Mesh
Input mesh
Returns
-------
Mesh
A new mesh containing simplicial elements
References
----------
Julien Dompierre, Paul Labbé, Marie-Gabrielle Vallet, Ricardo Camarero,
"How to subdivide pyramids, prisms and hexahedra into tetrahedra",
International Meshing, 1999
"""
# Elements that will be preserved
simplex = (
ED.Point_1,
ED.Bar_2,
ED.Bar_3,
ED.Triangle_3,
ED.Triangle_6,
ED.Tetrahedron_4,
ED.Tetrahedron_10,
)
# Elements that will be transformed
change = {
ED.Quadrangle_4: ED.Triangle_3,
ED.Quadrangle_9: ED.Triangle_6,
ED.Pyramid_5: ED.Tetrahedron_4,
ED.Pyramid_14: ED.Tetrahedron_10,
ED.Wedge_6: ED.Tetrahedron_4,
ED.Wedge_18: ED.Tetrahedron_10,
ED.Hexahedron_8: ED.Tetrahedron_4,
ED.Hexahedron_27: ED.Tetrahedron_10,
}
# Connectivity of the different simplexes created from each element type
# Nodes ordered according to VTK rules (https://vtk.org/doc/nightly/html)
II: Dict[ED.ElementType, List[List[Tuple[int, ...]]]] = {
ED.Quadrangle_4: [[(0, 1, 2), (0, 2, 3)], [(0, 1, 3), (1, 2, 3)]],
ED.Quadrangle_9: [
[(0, 1, 2, 4, 5, 8), (0, 2, 3, 8, 6, 7)],
[(0, 1, 3, 4, 8, 7), (1, 2, 3, 5, 6, 8)],
],
ED.Pyramid_5: [[(0, 1, 2, 4), (0, 2, 3, 4)], [(1, 2, 3, 4), (1, 3, 0, 4)]],
ED.Pyramid_14: [
[(0, 1, 2, 4, 5, 6, 13, 9, 10, 11), (0, 2, 3, 4, 13, 7, 8, 9, 11, 12)],
[(1, 2, 3, 4, 6, 7, 13, 10, 11, 12), (1, 3, 0, 4, 13, 8, 5, 10, 12, 9)],
],
ED.Wedge_6: [
[(0, 1, 2, 5), (0, 1, 5, 4), (0, 4, 5, 3)],
[(0, 1, 2, 4), (0, 4, 2, 5), (0, 4, 5, 3)],
],
ED.Wedge_18: [
[
(0, 1, 2, 5, 6, 7, 8, 17, 16, 14),
(0, 1, 5, 4, 6, 16, 17, 15, 13, 10),
(0, 4, 5, 3, 15, 10, 17, 12, 9, 11)
],
[
(0, 1, 2, 4, 6, 7, 8, 15, 13, 16),
(0, 4, 2, 5, 15, 16, 8, 17, 10, 14),
(0, 4, 5, 3, 15, 10, 17, 12, 9, 11)
],
],
ED.Hexahedron_8: [
[
(0, 1, 2, 5),
(0, 2, 7, 5),
(0, 2, 3, 7),
(0, 5, 7, 4),
(2, 7, 5, 6),
(-1, -1, -1, -1) # Dummy values
], # There are only 5 tetrahedra created in that case
[
(0, 5, 7, 4),
(0, 1, 7, 5),
(1, 6, 7, 5),
(0, 7, 2, 3),
(0, 7, 1, 2),
(1, 7, 6, 2)
],
[
(0, 4, 5, 6),
(0, 3, 7, 6),
(0, 7, 4, 6),
(0, 1, 2, 5),
(0, 3, 6, 2),
(0, 6, 5, 2)
],
[
(0, 2, 3, 6),
(0, 3, 7, 6),
(0, 7, 4, 6),
(0, 5, 6, 4),
(1, 5, 6, 0),
(1, 6, 2, 0)
],
],
ED.Hexahedron_27: [
[
(0, 1, 2, 5, 8, 9, 24, 22, 17, 21),
(0, 2, 7, 5, 24, 23, 20, 22, 21, 25),
(0, 2, 3, 7, 24, 10, 11, 20, 23, 19),
(0, 5, 7, 4, 22, 25, 20, 16, 12, 15),
(2, 7, 5, 6, 24, 21, 23, 18, 14, 13),
(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1) # Dummy values
], # There are only 5 tetrahedra created in that case
[
(0, 5, 7, 4, 22, 25, 20, 16, 12, 15),
(0, 1, 7, 5, 8, 26, 20, 22, 17, 25),
(1, 6, 7, 5, 21, 14, 26, 17, 13, 25),
(0, 7, 2, 3, 20, 23, 24, 11, 19, 10),
(0, 7, 1, 2, 20, 26, 8, 24, 23, 9),
(1, 7, 6, 2, 26, 14, 21, 9, 23, 18)
],
[
(0, 4, 5, 6, 16, 12, 22, 26, 25, 13),
(0, 3, 7, 6, 11, 19, 20, 26, 23, 14),
(0, 7, 4, 6, 20, 15, 16, 26, 14, 25),
(0, 1, 2, 5, 8, 9, 24, 22, 17, 21),
(0, 3, 6, 2, 11, 23, 26, 24, 10, 18),
(0, 6, 5, 2, 26, 13, 22, 24, 18, 21),
],
[
(0, 2, 3, 6, 24, 10, 11, 26, 18, 23),
(0, 3, 7, 6, 11, 19, 20, 26, 23, 14),
(0, 7, 4, 6, 20, 15, 16, 26, 14, 25),
(0, 5, 6, 4, 22, 13, 26, 16, 12, 25),
(1, 5, 6, 0, 17, 13, 21, 8, 22, 26),
(1, 6, 2, 0, 21, 18, 9, 8, 26, 24),
],
],
}
# Rotate the elements so that the node with lowest global index
# is at the bottom-left corner
JJ: Dict[ED.ElementType, List[Tuple[int, ...]]] = {
ED.Wedge_6: [
(0, 1, 2, 3, 4, 5),
(1, 2, 0, 4, 5, 3),
(2, 0, 1, 5, 3, 4),
(3, 5, 4, 0, 2, 1),
(4, 3, 5, 1, 0, 2),
(5, 4, 3, 2, 1, 0),
],
ED.Wedge_18: [
(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17),
(1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9, 13, 14, 12, 16, 17, 15),
(2, 0, 1, 5, 3, 4, 8, 6, 7, 11, 9, 10, 14, 12, 13, 17, 15, 16),
(3, 5, 4, 0, 2, 1, 11, 10, 9, 8, 7, 6, 12, 14, 13, 17, 16, 15),
(4, 3, 5, 1, 0, 2, 9, 11, 10, 6, 8, 7, 13, 12, 14, 15, 17, 16),
(5, 4, 3, 2, 1, 0, 10, 9, 11, 7, 6, 8, 14, 13, 12, 16, 15, 17),
],
ED.Hexahedron_8: [
(0, 1, 2, 3, 4, 5, 6, 7),
(1, 2, 3, 0, 5, 6, 7, 4),
(2, 3, 0, 1, 6, 7, 4, 5),
(3, 0, 1, 2, 7, 4, 5, 6),
(4, 7, 6, 5, 0, 3, 2, 1),
(5, 4, 7, 6, 1, 0, 3, 2),
(6, 5, 4, 7, 2, 1, 0, 3),
(7, 6, 5, 4, 3, 2, 1, 0),
],
ED.Hexahedron_27: [
(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26),
(1, 2, 3, 0, 5, 6, 7, 4, 9, 10, 11, 8, 13, 14, 15, 12, 17, 18, 19, 16, 22, 23, 21, 20, 24, 25, 26),
(2, 3, 0, 1, 6, 7, 4, 5, 10, 11, 8, 9, 14, 15, 12, 13, 18, 19, 16, 17, 21, 20, 23, 22, 24, 25, 26),
(3, 0, 1, 2, 7, 4, 5, 6, 11, 8, 9, 10, 15, 12, 13, 14, 19, 16, 17, 18, 23, 22, 20, 21, 24, 25, 26),
(4, 7, 6, 5, 0, 3, 2, 1, 15, 14, 13, 12, 11, 10, 9, 8, 16, 19, 18, 17, 22, 23, 20, 21, 25, 24, 26),
(5, 4, 7, 6, 1, 0, 3, 2, 12, 15, 14, 13, 8, 11, 10, 9, 17, 16, 19, 18, 21, 20, 22, 23, 25, 24, 26),
(6, 5, 4, 7, 2, 1, 0, 3, 13, 12, 15, 14, 9, 8, 11, 10, 18, 17, 16, 19, 23, 22, 21, 20, 25, 24, 26),
(7, 6, 5, 4, 3, 2, 1, 0, 14, 13, 12, 15, 10, 9, 8, 11, 19, 18, 17, 16, 20, 21, 23, 22, 25, 24, 26),
],
}
# Rotation of 120 or 240 degrees around the line between nodes 0 and 6
KK: Dict[ED.ElementType, List[Tuple[Tuple[int, int, int], Tuple[int, ...]]]] = {
ED.Hexahedron_8: [
((0, 0, 1), (0, 4, 5, 1, 3, 7, 6, 2)),
((1, 1, 0), (0, 4, 5, 1, 3, 7, 6, 2)),
((0, 1, 0), (0, 3, 7, 4, 1, 2, 6, 5)),
((1, 0, 1), (0, 3, 7, 4, 1, 2, 6, 5)),
],
ED.Hexahedron_27: [
(
(0, 0, 1),
(0, 4, 5, 1, 3, 7, 6, 2, 16, 12, 17, 8, 19, 14, 18, 10, 11, 15, 13, 9, 24, 25, 20, 21, 22, 23, 26)
),
(
(1, 1, 0),
(0, 4, 5, 1, 3, 7, 6, 2, 16, 12, 17, 8, 19, 14, 18, 10, 11, 15, 13, 9, 24, 25, 20, 21, 22, 23, 26)
),
(
(0, 1, 0),
(0, 3, 7, 4, 1, 2, 6, 5, 11, 19, 15, 16, 9, 18, 13, 17, 8, 10, 14, 12, 22, 23, 24, 25, 20, 21, 26)
),
(
(1, 0, 1),
(0, 3, 7, 4, 1, 2, 6, 5, 11, 19, 15, 16, 9, 18, 13, 17, 8, 10, 14, 12, 22, 23, 24, 25, 20, 21, 26)
),
],
}
newMesh = Mesh()
newMesh.elements = AllElements()
offset = mesh.ComputeGlobalOffset()
for eltype in mesh.elements.keys():
eloffset = offset[eltype]
if eltype in simplex:
# Copy simplexes from the input mesh into the newMesh
CopyElements(
newMesh.GetElementsOfType(eltype),
mesh.GetElementsOfType(eltype),
eloffset,
)
elif eltype in (ED.Quadrangle_4, ED.Quadrangle_8, ED.Pyramid_5 , ED.Pyramid_14):
newElt = newMesh.GetElementsOfType(change[eltype])
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts]
# Find the proper way to cut the element thanks to condI
I = np.array(II[eltype], dtype = MuscatIndex)
condI = np.argmin(V[:, :4], axis = 1) % 2
VI = np.empty((numOfElts, I.shape[1], I.shape[2]), dtype = MuscatIndex)
for k in range(I.shape[0]):
ind = (condI == k)
VI[ind] = V[ind][:, I[k]]
# Create simplexes with connectivity VI and the same id as the original element
VI = VI.reshape(I.shape[1] * numOfElts, I.shape[2])
ids = np.outer(
np.arange(eloffset, eloffset + numOfElts, dtype = MuscatIndex), np.ones(I.shape[1], dtype = MuscatIndex)
).flatten()
newElt.AddNewElements(VI, ids)
elif eltype in (ED.Wedge_6, ED.Wedge_18):
newElt = newMesh.GetElementsOfType(change[eltype])
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts]
# Apply a rotation depending on condJ
J = np.array(JJ[eltype], dtype = MuscatIndex)
condJ = np.argmin(V[:, :6], axis = 1)
VJ = np.empty_like(V, dtype = MuscatIndex)
for k in range(J.shape[0]):
ind = (condJ == k)
VJ[ind] = V[ind][:, J[k]]
# Find the proper way to cut the element thanks to condI
I = np.array(II[eltype], dtype = MuscatIndex)
condI = np.argmin(VJ[:, [1, 2, 5, 4]], axis = 1) % 2
VIJ = np.empty((numOfElts, I.shape[1], I.shape[2]), dtype = MuscatIndex)
for k in range(I.shape[0]):
ind = (condI == k)
VIJ[ind] = VJ[ind][:, I[k]]
# Create the tetrahedra with connectivity VIJ and the same id as the original element
VIJ = VIJ.reshape(I.shape[1] * numOfElts, I.shape[2])
ids = np.outer(
np.arange(eloffset, eloffset + numOfElts, dtype = MuscatIndex), np.ones(I.shape[1], dtype = MuscatIndex)
).flatten()
newElt.AddNewElements(VIJ, ids)
elif eltype in (ED.Hexahedron_8, ED.Hexahedron_27):
newElt = newMesh.GetElementsOfType(change[eltype])
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts]
# Apply a first rotation depending on condJ
J = np.array(JJ[eltype], dtype = MuscatIndex)
condJ = np.argmin(V[:, :8], axis = 1)
VJ = np.empty_like(V, dtype = MuscatIndex)
for k in range(J.shape[0]):
ind = (condJ == k)
VJ[ind] = V[ind][:, J[k]]
# Apply a second rotation depending on condK
K = KK[eltype]
faces: List[Tuple[int, int, int, int]] = [(2, 6, 5, 1), (7, 6, 2, 3), (5, 6, 7, 4)]
condK = np.argmin(VJ[:, faces], axis = 2) % 2
VJK = VJ
for k in range(len(K)):
ind = (condK == K[k][0]).all(1)
VJK[ind] = VJ[ind][:, K[k][1]]
# Find the proper way to cut the element thanks to condI
I = np.array(II[eltype], dtype = MuscatIndex)
condI = np.sum(condK, axis = 1)
VIJK = np.empty((numOfElts, I.shape[1], I.shape[2]), dtype = MuscatIndex)
for k in range(I.shape[0]):
ind = (condI == k)
VIJK[ind] = VJK[ind][:, I[k]]
# Create the tetrahedra with connectivity VIJ and the same id as the original element
VIJK = VIJK.reshape((I.shape[1] * numOfElts, I.shape[2]))
ids = np.outer(
np.arange(eloffset, eloffset + numOfElts, dtype = MuscatIndex), np.ones(I.shape[1], dtype = MuscatIndex)
).flatten()
# Handle the case where only 5 tetrahedra are created by removing the dummy tet
cond0 = (VIJK[:, 0] != VIJK[:, 1])
VIJK0 = VIJK[cond0]
ids0 = ids[cond0]
newElt.AddNewElements(VIJK0, ids0)
else:
print(f"Elements of type {eltype} cannot be converted currently")
print("Use function CompleteQuadraticMesh or LinearizeQuadraticMesh first")
CopyElements(
newMesh.GetElementsOfType(eltype),
mesh.GetElementsOfType(eltype),
eloffset,
)
# Copy nodes from the input mesh in the newMesh
newMesh.nodes = mesh.nodes
newMesh.originalIDNodes = np.arange(mesh.GetNumberOfNodes())
# Transfer node tags
for tag in mesh.nodesTags:
newMesh.nodesTags.CreateTag(tag.name, False).SetIds(tag.GetIds())
TransferElementTags(mesh, newMesh)
# Transfer node Fields
for key, val in mesh.nodeFields.items():
newMesh.nodeFields[key] = val
# Transfer element Fields
TransferElementFields(mesh, newMesh)
newMesh.PrepareForOutput()
return newMesh
[docs]def FromLinearToIncompleteQuad(mesh: Mesh) -> Mesh:
"""Convert linear elements to incomplete quadratic elements
Warning : originalIds of the new mesh are related to the global number of the elements in the input mesh
Parameters
----------
mesh : Mesh
the input mesh
Returns
-------
Mesh
a new mesh
"""
complete: Dict[ED.ElementType, ED.ElementType] = {
ED.Bar_2: ED.Bar_3,
ED.Triangle_3: ED.Triangle_6,
ED.Tetrahedron_4: ED.Tetrahedron_10,
ED.Quadrangle_4: ED.Quadrangle_8,
ED.Hexahedron_8: ED.Hexahedron_20,
ED.Wedge_6: ED.Wedge_15,
ED.Pyramid_5: ED.Pyramid_13
}
# edges of the elements where to add nodes
# Nodes ordered according to VTK rules (https://vtk.org/doc/nightly/html)
faces: Dict[ED.ElementType, List[Tuple[int, ...]]] = {
ED.Bar_2: [(0, 1)],
ED.Triangle_3: [(0, 1), (1, 2), (2, 0)],
ED.Quadrangle_4: [(0, 1), (1, 2), (2, 3), (3, 0)],
ED.Tetrahedron_4: [(0, 1), (1, 2), (2, 0), (0, 3), (1, 3), (2, 3)],
ED.Hexahedron_8: [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)],
ED.Wedge_6: [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (0, 3), (1, 4), (2, 5)],
ED.Pyramid_5: [(0, 1), (1, 2), (2, 3), (3, 0), (0, 4), (1, 4), (2, 4), (3, 4)]
}
offset = mesh.ComputeGlobalOffset()
numOfNodes = mesh.GetNumberOfNodes()
newMesh = Mesh()
newNodesNum: List[int] = [] # Global numbers of the created nodes
newNodesFace: List[NDArray] = [] # Faces where a new node has been added
newNodesCoord: List[MuscatIndex] = [] # Coordinates of new created nodes
newMesh.elements = AllElements()
for eltype in mesh.elements.keys():
eloffset = offset[eltype]
if eltype in complete.keys():
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts] # global numbering of the nodes
F = faces[eltype]
point = np.zeros((numOfElts, len(F)), dtype=MuscatIndex)
newElt = newMesh.GetElementsOfType(complete[eltype])
for i, Fi in enumerate(F):
VF = np.sort(V[:, Fi], axis=1) # Sort the global numbers of the nodes in edge i
# to identify if we have already added this node
sumCond = np.zeros(numOfElts, dtype=np.bool_)
for num, face in zip(newNodesNum, newNodesFace):
cond = np.empty_like(sumCond, dtype=np.bool_)
np.all(VF == face, axis=1, out=cond)
point[cond, i] = num
np.logical_or(sumCond, cond, out=sumCond)
np.logical_not(sumCond, out=sumCond)
VFsum = VF[sumCond]
face, inv = np.unique(VFsum, return_inverse=True, axis=0)
indCount = face.shape[0]
rg = np.arange(numOfNodes, numOfNodes + indCount)
newNodesNum.extend(rg)
newNodesFace.extend(face)
newNodesCoord.extend(0.5*np.sum(mesh.nodes[face], axis=1))
numOfNodes += indCount
point[sumCond, i] = rg[inv]
# Create new elements with additional nodes
newElt.AddNewElements(
np.append(V, point, axis=1), range(eloffset, eloffset + numOfElts)
)
else:
# Copy elements of the input mesh in the newMesh
CopyElements(
newMesh.GetElementsOfType(eltype),
mesh.GetElementsOfType(eltype),
eloffset,
)
# Copy nodes from input mesh and add new created nodes
newMesh.nodes = np.zeros((numOfNodes, mesh.GetPointsDimensionality()), dtype=MuscatFloat)
newMesh.nodes[:mesh.GetNumberOfNodes()] = mesh.nodes
if newNodesNum:
newMesh.nodes[mesh.GetNumberOfNodes():] = newNodesCoord
newMesh.originalIDNodes = np.arange(numOfNodes)
# Transfer nodeTags
for tag in mesh.nodesTags:
newMesh.nodesTags.CreateTag(tag.name).SetIds(tag.GetIds())
# if all the nodes of a face share the same tag, the new node created in the middle of the face will also have this tag
for num, face in zip(newNodesNum, newNodesFace):
if np.isin(face, tag.GetIds()).all():
newMesh.nodesTags[tag.name].AddToTag(num)
# Transfer elementTags
TransferElementTags(mesh, newMesh)
# Transfer node fields
for key, val in mesh.nodeFields.items():
field = np.zeros((numOfNodes, val.shape[1]))
field[:mesh.GetNumberOfNodes()] = val
for num, face in zip(
newNodesNum, newNodesFace
):
field[num] = 0.5*np.sum(val[face], axis=0)
newMesh.nodeFields[key] = field
# Transfer element fields
TransferElementFields(mesh, newMesh)
newMesh.PrepareForOutput()
return newMesh
[docs]def FromIncompleteQuadToQuad(mesh: Mesh) -> Mesh:
"""Convert incomplete quadratic elements to full quadratic elements
Warning : originalIds of the new mesh are related to the global number of the elements in the input mesh
Parameters
----------
mesh : Mesh
the input mesh
Returns
-------
Mesh
a new mesh
"""
complete: Dict[ED.ElementType, ED.ElementType] = {
ED.Quadrangle_8: ED.Quadrangle_9,
ED.Hexahedron_20: ED.Hexahedron_27,
ED.Wedge_15: ED.Wedge_18,
ED.Pyramid_13: ED.Pyramid_14,
}
# Faces of the elements where to add nodes
# Nodes ordered according to VTK rules (https://vtk.org/doc/nightly/html)
faces: Dict[ED.ElementType, List[Tuple[int, ...]]] = {
ED.Quadrangle_8: [(0, 1, 2, 3, 4, 5, 6, 7)],
ED.Hexahedron_20: [
(3, 0, 4, 7, 11, 16, 15, 19),
(1, 2, 6, 5, 9, 18, 13, 17),
(0, 1, 5, 4, 8, 17, 12, 16),
(2, 3, 7, 6, 10, 19, 14, 18),
(0, 3, 2, 1, 11, 10, 9, 8),
(4, 5, 6, 7, 12, 13, 14, 15),
(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19),
],
ED.Wedge_15: [
(0, 1, 4, 3, 6, 13, 9, 12),
(1, 2, 5, 4, 7, 14, 10, 13),
(2, 0, 3, 5, 8, 12, 11, 14),
],
ED.Pyramid_13: [(0, 1, 2, 3, 5, 6, 7, 8)],
}
spaceQuad8 = Quad8_P2_Lagrange()
spaceQuad8.Create()
spaceHex20 = Hexa20_P2_Lagrange()
spaceHex20.Create()
shapeFun = {
"8": spaceQuad8.GetShapeFunc([0.5, 0.5]),
"20": spaceHex20.GetShapeFunc([0.5, 0.5, 0.5]),
}
offset = mesh.ComputeGlobalOffset()
numOfNodes = mesh.GetNumberOfNodes()
newMesh = Mesh()
newNodesNum: List[int] = [] # Global numbers of the created nodes
newNodesFace: List[NDArray] = [] # Faces where a new node has been added
newNodesCoord: List[MuscatIndex] = [] # Coordinates of new created nodes
newMesh.elements = AllElements()
for eltype in mesh.elements.keys():
eloffset = offset[eltype]
if eltype in complete.keys():
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts] # global numbering of the nodes
F = faces[eltype]
point = np.zeros((numOfElts, len(F)), dtype=MuscatIndex)
newElt = newMesh.GetElementsOfType(complete[eltype])
for i, Fi in enumerate(F):
VF = V[:, Fi]
sumCond = np.zeros(numOfElts, dtype=np.bool_)
for num, face in zip(newNodesNum, newNodesFace):
# Sort the global numbers of the nodes in face i
# to identify if we have already added this node
cond = np.empty_like(sumCond, dtype=np.bool_)
np.all((np.sort(VF[:, :8], axis=1) == np.sort(face[:8])), axis=1, out=cond)
point[cond, i] = num
np.logical_or(sumCond, cond, out=sumCond)
np.logical_not(sumCond, out=sumCond)
indCount = np.count_nonzero(sumCond)
face = VF[sumCond]
rg = range(numOfNodes, numOfNodes + indCount)
newNodesNum.extend(rg)
newNodesFace.extend(face)
newNodesCoord.extend(
np.tensordot(mesh.nodes[face], shapeFun[str(len(Fi))], axes=(1, 0))
)
numOfNodes += indCount
point[sumCond, i] = rg
# Create new elements with additional nodes
newElt.AddNewElements(
np.append(V, point, axis=1), range(eloffset, eloffset + numOfElts)
)
elif eltype in complete.keys():
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
V = elt.connectivity[:numOfElts] # global numbering of the nodes
F = faces[eltype]
point = np.zeros((numOfElts, len(F)), dtype=MuscatIndex)
newElt = newMesh.GetElementsOfType(complete[eltype])
for i, Fi in enumerate(F):
VF = V[:, Fi]
sumCond = np.zeros(numOfElts, dtype=np.bool_)
for num, face in zip(newNodesNum, newNodesFace):
# Sort the global numbers of the nodes in face i
# to identify if we have already added this node
cond = np.empty_like(sumCond, dtype=np.bool_)
np.all((np.sort(VF[:, :8], axis=1) == np.sort(face[:8])), axis=1, out=cond)
point[cond, i] = num
np.logical_or(sumCond, cond, out=sumCond)
np.logical_not(sumCond, out=sumCond)
indCount = np.count_nonzero(sumCond)
face = VF[sumCond]
rg = range(numOfNodes, numOfNodes + indCount)
newNodesNum.extend(rg)
newNodesFace.extend(face)
newNodesCoord.extend(
np.tensordot(mesh.nodes[face], shapeFun[str(len(Fi))], axes=(1, 0))
)
numOfNodes += indCount
point[sumCond, i] = rg
# Create new elements with additional nodes
newElt.AddNewElements(
np.append(V, point, axis=1), range(eloffset, eloffset + numOfElts)
)
else:
# Copy elements of the input mesh in the newMesh
CopyElements(
newMesh.GetElementsOfType(eltype),
mesh.GetElementsOfType(eltype),
eloffset,
)
# Copy nodes from input mesh and add new created nodes
newMesh.nodes = np.zeros((numOfNodes, mesh.GetPointsDimensionality()), dtype=MuscatFloat)
newMesh.nodes[:mesh.GetNumberOfNodes()] = mesh.nodes
if newNodesNum:
newMesh.nodes[mesh.GetNumberOfNodes():] = newNodesCoord
newMesh.originalIDNodes = np.arange(numOfNodes)
# Transfer nodeTags
for tag in mesh.nodesTags:
newMesh.nodesTags.CreateTag(tag.name).SetIds(tag.GetIds())
# if all the nodes of a face share the same tag, the new node created in the middle of the face will also have this tag
for num, face in zip(
newNodesNum, newNodesFace
):
if np.isin(face, tag.GetIds()).all():
newMesh.nodesTags[tag.name].AddToTag(num)
# Transfer elementTags
TransferElementTags(mesh, newMesh)
# Transfer node fields
for key, val in mesh.nodeFields.items():
field = np.zeros((numOfNodes, val.shape[1]))
field[:mesh.GetNumberOfNodes()] = val
for num, face in zip(
newNodesNum, newNodesFace
):
field[num] = np.tensordot(val[face], shapeFun[str(len(face))], axes=[0, 0])
newMesh.nodeFields[key] = field
# Transfer element fields
TransferElementFields(mesh, newMesh)
newMesh.PrepareForOutput()
return newMesh
[docs]def FromQuadToLinear(mesh: Mesh) -> Mesh:
"""Convert incomplete or complete quadratic elements in a mesh into linear elements
Warning : originalIds of the new mesh are related to the global number of the elements in the input mesh
Parameters
----------
mesh : Mesh
the input mesh
Returns
-------
Mesh
a new mesh
"""
dataBase: Dict[ED.ElementType, Tuple[ED.ElementType, int]] = {}
dataBase[ED.Tetrahedron_10] = (ED.Tetrahedron_4, 4)
dataBase[ED.Triangle_6] = (ED.Triangle_3, 3)
dataBase[ED.Quadrangle_8] = (ED.Quadrangle_4, 4)
dataBase[ED.Hexahedron_20] = (ED.Hexahedron_8, 8)
dataBase[ED.Bar_3] = (ED.Bar_2, 2)
dataBase[ED.Pyramid_13] = (ED.Pyramid_5, 5)
dataBase[ED.Wedge_15] = (ED.Wedge_6, 6)
dataBase[ED.Quadrangle_9] = (ED.Quadrangle_4, 4)
dataBase[ED.Hexahedron_27] = (ED.Hexahedron_8, 8)
dataBase[ED.Pyramid_14] = (ED.Pyramid_5, 5)
dataBase[ED.Wedge_18] = (ED.Wedge_6, 6)
offset = mesh.ComputeGlobalOffset()
newMesh = Mesh()
newMesh.elements = AllElements()
keepNodes = np.ones(mesh.GetNumberOfNodes(), dtype=np.bool_)
for eltype in mesh.elements.keys():
if eltype in dataBase.keys():
elt = mesh.GetElementsOfType(eltype)
remove = elt.connectivity[:, dataBase[eltype][1]:]
keepNodes[remove.flatten()] = False
# Copy nodes from input mesh and remove useless nodes
rang = np.arange(mesh.GetNumberOfNodes(), dtype=MuscatIndex)
saveIds = rang[keepNodes]
newMesh.nodes = mesh.nodes[saveIds, :]
newMesh.originalIDNodes = saveIds
perm = np.ones(mesh.GetNumberOfNodes(), dtype=MuscatIndex)*mesh.GetNumberOfNodes()
perm[saveIds] = np.arange(saveIds.shape[0])
for eltype in mesh.elements.keys():
elt = mesh.GetElementsOfType(eltype)
numOfElts = elt.GetNumberOfElements()
eltOffset = offset[eltype]
if eltype in dataBase.keys():
newElt = newMesh.GetElementsOfType(dataBase[eltype][0])
V = elt.connectivity[:numOfElts, :dataBase[eltype][1]]
else:
newElt = newMesh.GetElementsOfType(eltype)
V = elt.connectivity[:numOfElts, :]
newElt.Reserve(newElt.cpt + numOfElts)
newElt.connectivity[newElt.cpt:, :] = perm[V]
newElt.originalIds[newElt.cpt:] = np.arange(eltOffset, numOfElts + eltOffset)
newElt.cpt += numOfElts
# Transfer nodeTags
for tag in mesh.nodesTags:
ids = perm[tag.GetIds()[keepNodes[tag.GetIds()]]]
newMesh.nodesTags.CreateTag(tag.name, errorIfAlreadyCreated=False).SetIds(ids)
# Transfer elementTags
TransferElementTags(mesh, newMesh)
# Transfer node fields
for key, val in mesh.nodeFields.items():
newMesh.nodeFields[key] = val[saveIds]
# Transfer element fields
TransferElementFields(mesh, newMesh)
newMesh.PrepareForOutput()
return newMesh
[docs]def CopyElements(
this: ElementsContainer, other: ElementsContainer, gidOffset=MuscatIndex(0)
) -> None:
"""
Copy the elements from the container `other` into container `this`.
Warning : the ids of the nodes must be the same in the two containers
Warning : here the `originalIds` of the elements relate to the global
numbering of the elements in other whereas in `ElementsContainer.Merge` they relate
to the local numbering.
"""
if this is other:
raise RuntimeError("this and other are the same object")
if other.cpt == 0:
return
this.Reserve(this.cpt + other.cpt)
this.connectivity[this.cpt:, :] = other.connectivity[:other.cpt, :]
this.originalIds[this.cpt:] = np.arange(gidOffset, other.cpt + gidOffset)
this.cpt += other.cpt
[docs]def TransferElementFields(inMesh: Mesh, outMesh: Mesh) -> None:
"""Function to Transfer elemFields from the original mesh to the derived mesh
Warning : originalIds of outmesh are related to the global number of the elements in inmesh
Parameters
----------
inMesh : Mesh
the source mesh, we extract the fields from inMesh.elemFields.
outMesh : Mesh
The target mesh, we push the new fields into outMesh.elemFields.
"""
# Compute the transfer array for the elemFields
cpt2 = 0
oid = np.empty(outMesh.GetNumberOfElements(), dtype=int)
for name in outMesh.elements.keys():
elements = outMesh.elements[name]
oid[cpt2:cpt2 + elements.GetNumberOfElements()] = (
elements.originalIds
) # global number of the elements in the inMesh
cpt2 += elements.GetNumberOfElements()
for key, val in inMesh.elemFields.items():
outMesh.elemFields[key] = val[oid]
[docs]def TestMeshValidity(mesh):
for eltype in mesh.elements.keys():
ff = ElementFilter(elementType=eltype)
vol = GetVolumePerElement(mesh, elementFilter=ff)
if (vol<=0).any():
nb = np.sum(vol<=0)
print("Warning:", nb, "elements of type", eltype, "have non positiv volume")
return
[docs]def CheckIntegrity():
import time
from Muscat.IO.GeofReader import ReadGeof
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
from Muscat.TestData import GetTestDataPath
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tic = time.time()
filename = GetTestDataPath() + "multielements3D.geof"
mesh = ReadGeof(filename)
mesh.nodeFields["coor"] = mesh.nodes
mesh.ConvertDataForNativeTreatment()
print(mesh)
TestMeshValidity(mesh)
filename = TemporaryDirectory().GetTempPath() + "tetrahedrizationIn.xdmf"
PointFieldsNames = list(mesh.nodeFields.keys())
PointFields = list(mesh.nodeFields.values())
WriteMeshToXdmf(
filename, mesh, PointFieldsNames=PointFieldsNames, PointFields=PointFields
)
nmesh = FromLinearToIncompleteQuad(mesh)
print(nmesh)
TestMeshValidity(nmesh)
nmesh = FromIncompleteQuadToQuad(nmesh)
print(nmesh)
TestMeshValidity(nmesh)
newMesh = Tetrahedrization(nmesh)
print(newMesh)
TestMeshValidity(newMesh)
filename = TemporaryDirectory().GetTempPath() + "tetrahedrizationQuadOut.xdmf"
PointFieldsNames = list(newMesh.nodeFields.keys())
PointFields = list(newMesh.nodeFields.values())
WriteMeshToXdmf(
filename, newMesh, PointFieldsNames = PointFieldsNames, PointFields = PointFields
)
nmesh = FromQuadToLinear(mesh)
print(nmesh)
TestMeshValidity(nmesh)
newMesh = Tetrahedrization(nmesh)
print(newMesh)
TestMeshValidity(newMesh)
filename = TemporaryDirectory().GetTempPath() + "tetrahedrizationLinOut.xdmf"
PointFieldsNames = list(newMesh.nodeFields.keys())
PointFields = list(newMesh.nodeFields.values())
WriteMeshToXdmf(
filename, newMesh, PointFieldsNames = PointFieldsNames, PointFields = PointFields
)
print(f"Execution time: {int(time.time() - tic)} seconds")
return "ok"
if __name__ == "__main__":
print(CheckIntegrity()) # pragma: no cover