# -*- 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 Union, Optional, Dict
import copy
import warnings
import numpy as np
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Mesh import Mesh, ElementsContainer, AllElements
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshTools.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh, GetMonoIndex
from Muscat.MeshTools.MeshModificationTools import ComputeSkin
from Muscat.MeshTools.MeshModificationTools import CleanLonelyNodes
[docs]
def CreateMeshOfTriangles(points: ArrayLike, triangles: ArrayLike) -> Mesh:
"""Create a Unstructured mesh using only points
and the connectivity matrix for the triangles.
Parameters
----------
points : ArrayLike
The positions of the point (nb Points,3)
triangles : ArrayLike
The connectivity of the triangles (nb Triangles, 3)
Returns
-------
Mesh
A instance of Mesh with the triangles
"""
return CreateMeshOf(points, triangles, elemName=ED.Triangle_3)
[docs]
def CreateMeshOf(points: ArrayLike, connectivity: ArrayLike, elemName: ED.ElementType) -> Mesh:
"""Helper function to create a mesh of homogeneous elements.
Parameters
----------
points : ArrayLike
The positions of the point (nb Points,3)
connectivity : ArrayLike
The connectivity of the elements (nb of elements, number of point per element)
elemName : str, optional
the name of the elements to create. , by default None
Returns
-------
Mesh
An instance of Mesh with the given point and elements
"""
res = Mesh()
with res.WithModification():
res.nodes = np.array(points, dtype=MuscatFloat)
res.originalIDNodes = np.arange(
0, res.GetNumberOfNodes(), dtype=MuscatIndex)
elements = res.GetElementsOfType(elemName)
elements.connectivity = np.asarray(connectivity, dtype=MuscatIndex)
elements.originalIds = np.arange(
0, elements.connectivity.shape[0], dtype=MuscatIndex)
elements.cpt = elements.connectivity.shape[0]
elements.tags.CreateTag(str(
ED.dimensionality[elemName]) + "D").SetIds(np.arange(elements.GetNumberOfElements()))
return res
[docs]
def CreateSquare(dimensions: ArrayLike = [2, 2], origin: ArrayLike = [-1.0, -1.0], spacing: ArrayLike = [1.0, 1.0], ofTriangles: bool = False) -> Mesh:
"""Create a mesh of a square
Parameters
----------
dimensions : ArrayLike, optional
Number of points in every direction, by default [2,2]
origin : ArrayLike, optional
origin of the mesh , by default [-1.0,-1.0]
spacing : ArrayLike, optional
x and y size of every element, by default [1.,1.]
ofTriangles : bool, optional
, by default False
Returns
-------
Mesh
_description_
"""
dimensions = np.asarray(dimensions, dtype=MuscatIndex)
spacing = np.asarray(spacing, dtype=float)
origin = np.asarray(origin, dtype=float)
myMesh = CreateConstantRectilinearMesh(
dimensions=dimensions, origin=origin, spacing=spacing)
with myMesh.WithModification(verifyIntegrity=False):
# corners
d = dimensions - 1
s = spacing
indices = [[0, 0, 0], [d[0], 0, 0], [0, d[1], 0], [d[0], d[1], 0]]
for n in indices:
idx = GetMonoIndex(n, dimensions)
name = "x" + ("0" if n[0] == 0 else "1")
name += "y" + ("0" if n[1] == 0 else "1")
myMesh.nodesTags.CreateTag(name, False).SetIds([idx])
if ofTriangles:
from Muscat.MeshTools.MeshTetrahedrization import Tetrahedrization
myMesh = Tetrahedrization(myMesh)
ComputeSkin(myMesh, inPlace=True)
for selection in ElementFilter(dimensionality=2)(myMesh):
selection.elements.GetTag("2D").SetIds(selection.indices)
mMin, mMax = myMesh.ComputeBoundingBox()
tol = np.min(spacing) / 10
for coordinateI in range(2):
for suffix in [0, 1]:
name = ["X", "Y"][coordinateI] + str(suffix)
val = mMin[coordinateI] * \
(1 - suffix) + mMax[coordinateI] * suffix
for selection in ElementFilter(dimensionality=1, zone=lambda xyz: abs(xyz[:, coordinateI] - val) - tol)(myMesh):
selection.elements.GetTag(name).SetIds(selection.indices)
return myMesh
[docs]
def CreateDisk(
nr: MuscatIndex = 10, nTheta: MuscatIndex = 10, r0: MuscatFloat = 0.5, r1: MuscatFloat = 1.0, theta0: MuscatFloat = 0, theta1: MuscatFloat = np.pi / 2, ofTriangles: bool = False
) -> Mesh:
"""Function to create a disk section
Parameters
----------
nr : MuscatIndex, optional
number of points in the radial direction, by default 10
nTheta : MuscatIndex, optional
number of point in the angular direction, by default 10
r0 : MuscatFloat, optional
internal radius, by default 0.5
r1 : MuscatFloat, optional
external radius, by default 1.
theta0 : MuscatFloat, optional
start angle, by default 0
theta1 : MuscatFloat, optional
end angle, by default np.pi/2
ofTriangles : bool, optional
if false, quand are used for the mesh , by default False
if true, triangles are used for the mesh
Returns
-------
Mesh
An Mesh of a disk sector
"""
myMesh = CreateSquare(dimensions=[nr, nTheta], origin=[r0, theta0], spacing=[(
r1 - r0) / (nr - 1), (theta1 - theta0) / (nTheta - 1)], ofTriangles=ofTriangles)
r = myMesh.nodes[:, 0].copy()
theta = myMesh.nodes[:, 1].copy()
myMesh.nodes[:, 0] = r * np.cos(theta)
myMesh.nodes[:, 1] = r * np.sin(theta)
myMesh.elements[ED.Bar_2].tags["X0"].name = "R0"
myMesh.elements[ED.Bar_2].tags["X1"].name = "R1"
myMesh.elements[ED.Bar_2].tags["Y0"].name = "Theta0"
myMesh.elements[ED.Bar_2].tags["Y1"].name = "Theta1"
return myMesh
[docs]
def CreateCube(dimensions: ArrayLike = [2, 2, 2], origin: ArrayLike = [-1.0, -1.0, -1.0], spacing: ArrayLike = [1.0, 1.0, 1.0], ofTetras: bool = False) -> Mesh:
"""Create a Mesh of a cube
Parameters
----------
dimensions : ArrayLike, optional
Number of point in every direction, by default [2,2,2]
origin : ArrayLike, optional
origin of the lowest point, by default [-1.0,-1.0,-1.0]
spacing : ArrayLike, optional
size in the 3 direction of every element, by default [1.,1.,1.]
ofTetras : bool, optional
if False: mesh of hexahedrons
if True: mesh of tetrahedrons
by default False
Returns
-------
Mesh
A mesh of the cube
"""
dimensions = np.asarray(dimensions, dtype=MuscatIndex)
spacing = np.asarray(spacing, dtype=MuscatFloat)
origin = np.asarray(origin, dtype=MuscatFloat)
myMesh = CreateConstantRectilinearMesh(
dimensions=dimensions, origin=origin, spacing=spacing)
with myMesh.WithModification(verifyIntegrity=False):
# corners
d = dimensions - 1
s = spacing
indices = [[0, 0, 0], [d[0], 0, 0], [0, d[1], 0], [d[0], d[1], 0], [
0, 0, d[2]], [d[0], 0, d[2]], [0, d[1], d[2]], [d[0], d[1], d[2]]]
for n in indices:
idx = GetMonoIndex(n, dimensions)
name = "x" + ("0" if n[0] == 0 else "1")
name += "y" + ("0" if n[1] == 0 else "1")
name += "z" + ("0" if n[2] == 0 else "1")
myMesh.nodesTags.CreateTag(name, False).SetIds([idx], unique=True)
if ofTetras:
from Muscat.MeshTools.MeshTetrahedrization import Tetrahedrization
myMesh = Tetrahedrization(myMesh)
ComputeSkin(myMesh, inPlace=True)
for selection in ElementFilter(dimensionality=3)(myMesh):
selection.elements.GetTag("3D").SetIds(
selection.indices, unique=True)
mMin, mMax = myMesh.ComputeBoundingBox()
tol = np.min(spacing) / 10
for coordinateI in range(3):
for suffix in [0, 1]:
name = ["X", "Y", "Z"][coordinateI] + str(suffix)
val = mMin[coordinateI] * \
(1 - suffix) + mMax[coordinateI] * suffix
for selection in ElementFilter(dimensionality=2, zone=lambda xyz: abs(xyz[:, coordinateI] - val) - tol)(myMesh):
selection.elements.GetTag(name).SetIds(
selection.indices, unique=True)
return myMesh
[docs]
def CreateMeshFromCellsDict(
points: ArrayLike, cellsDict: Dict[str, ArrayLike], pointFields: Optional[Dict[str, ArrayLike]] = None, cellFields: Optional[Dict[str, ArrayLike]] = None, cellFieldInDictOrder: bool = True, originalIds: Optional[ArrayLike] = None, cellsOriginalIds: Dict[str, ArrayLike] = {}
) -> Mesh:
"""Function to create a Mesh from points and a dict of element
Data is converted of the native format (MuscatFloat, MuscatIndex) if needed
Parameters
----------
points : ArrayLike
Nodes coordinates of size [nb_points x [2 or3]]
cellsDict : dict[Union[ElementType,str],ArrayLike]
A dict with (key= elementType, value=connectivity), connectivity start from 0 pointing to the node indices
pointFields : dict[str,ArrayLike], optional
fields at nodes each fields must be of size [nb_points, nb_components], by default None
cellFields : dict[str,ArrayLike], optional
fields at elements each fields must be of size [total_number_of_elements, nb_components], by default None
_description_, by default None
cellFieldInDictOrder : bool, by default True
if True the cell fields are interpreted in the cells_dict insertion order (insertion order is preserved in python)
if False the cell field are interpreted in the canonic Muscat ("alphanumeric") order.
originalIds : ArraLike, optional
Nodes original IDs of size [nb_points x 1], generated if not specified
cellsOriginalIds : dict[Union[ElementType,str], ArraLike], default {}
Cells original IDs. If any element type is not specified in the dictionary, the original IDs will be generated
Returns
-------
Mesh
A Mesh with the field inside
"""
mesh = Mesh()
if originalIds is None:
mesh.SetNodes(points, generateOriginalIDs=True)
else:
mesh.SetNodes(points, originalIDNodes=originalIds,
generateOriginalIDs=False)
cpt = 0
cellsDict = dict(zip(
[ED.ElementType[k] if k is str else k for k in cellsDict.keys()], cellsDict.values()))
cellsOriginalIds = dict(zip(
[ED.ElementType[k] if k is str else k for k in cellsOriginalIds.keys()], cellsOriginalIds.values()))
for cell_type, cell_conn in cellsDict.items():
cell_conn = np.asarray(cell_conn, dtype=MuscatIndex)
mesh.elements.GetElementsOfType(cell_type).AddNewElements(
cell_conn, originalIds=cellsOriginalIds.get(cell_type, np.arange(cpt, cpt + cell_conn.shape[0])))
if pointFields is not None:
for pointFieldName, pointField in pointFields.items():
mesh.nodeFields[pointFieldName] = pointField
if cellFields is not None:
if cellFieldInDictOrder:
offset = mesh.ComputeGlobalOffset()
for cellFieldName, cellField in cellFields.items():
cellField = np.asarray(cellField)
data = np.empty_like(cellField)
cpt = 0
for cell_type, cell_conn in cellsDict.items():
elements = mesh.elements.GetElementsOfType(cell_type)
nbElems = elements.GetNumberOfElements()
data[offset[elements.elementType]: offset[elements.elementType] +
nbElems] = cellField[cpt: cpt + nbElems]
cpt += nbElems
mesh.elemFields[cellFieldName] = data
else:
for cellFieldName, cellField in cellFields.items():
mesh.elemFields[cellFieldName] = cellFields[cellFieldName]
return mesh
[docs]
def MeshToSimplex(mesh: Mesh, inPlace=True) -> Mesh:
"""(EXPERIMENTAL) Convert mesh to only tetrahedron/triangle/bars/points
Warning!!! This function is dangerous, we don't check the compatibility at each interface between
the elements.
Parameters
----------
mesh : Mesh
A mesh
inPlace: bool, optional, by default True
Returns
-------
Mesh
A mesh composed only by simplex (this is the mesh inplace is True)
"""
warnings.warn("MeshToSimplex will be removed in future version Muscat 2.4, please use Muscat.MeshTools.MeshTetrahedrization.Tetrahedrization function", DeprecationWarning)
from Muscat.MeshContainers.Mesh import ElementsContainer, AllElements
ae = AllElements()
dataBase = {}
dataBase[ED.Hexahedron_8] = (ED.Tetrahedron_4, [[0, 6, 2, 3], [0, 6, 3, 7], [
0, 6, 7, 4], [0, 6, 4, 5], [0, 6, 5, 1], [0, 6, 1, 2]])
dataBase[ED.Quadrangle_4] = (ED.Triangle_3, [[0, 1, 2], [0, 2, 3]])
for elementName, data in mesh.elements.items():
res = data
if elementName in dataBase:
nbElements = data.GetNumberOfElements()
if nbElements == 0:
continue
newElementType, extractor = dataBase[elementName]
nbNE = len(extractor) # number Of New Elements Per Element
res = ElementsContainer(newElementType)
nbElements = data.GetNumberOfElements()
res.Allocate(nbElements * nbNE)
conn = data.connectivity
for i in range(nbNE):
res.connectivity[i: nbElements *
nbNE: nbNE] = conn[:, extractor[i]]
res.originalIds = np.repeat(data.originalIds, nbNE)
for tagName in data.tags.keys():
ids = data.tags[tagName].GetIds()
res.tags.CreateTag(tagName).SetIds(
np.repeat(ids, nbNE) * nbNE + np.tile(range(nbNE), len(ids)))
elif elementName in [ED.Triangle_3, ED.Triangle_6, ED.Tetrahedron_4, ED.Tetrahedron_10, ED.Bar_2, ED.Bar_3, ED.Point_1]:
pass
else:
raise (Exception("Don't know how to convert {} to simplices".format(
elementName))) # pragma: no cover
if res.elementType in ae:
ae[res.elementType].Merge(res)
else:
ae.AddContainer(res)
if inPlace:
mesh.elements = ae
mesh.props.pop("IsConstantRectilinear", None)
return mesh
else:
res = Mesh()
res.elements = ae
res.nodes = mesh.nodes.view()
res.nodes.flags.writeable = False
res.originalIDNodes = mesh.originalIDNodes.view()
res.originalIDNodes.flags.writeable = False
res.nodesTags = res.nodesTags.Copy()
res.nodeFields = mesh.nodeFields
res.props.pop("IsConstantRectilinear", None)
return res
[docs]
def ToQuadraticMesh(inputMesh: Mesh) -> Mesh:
"""Function to convert any linear mesh to a quadratic mesh.
Nodes fields and element fields are lost.
Parameters
----------
inputMesh : Mesh
the in put mesh
Returns
-------
Mesh
A mesh composed only by quadratic elements
"""
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP2
from Muscat.FE.Fields.FEField import FEField
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.FE.Fields.FieldTools import FillFEField
dim = inputMesh.GetPointsDimensionality()
numbering = ComputeDofNumbering(inputMesh, LagrangeSpaceP2)
res = Mesh()
res.nodes = np.empty((numbering.size, dim), dtype=MuscatFloat)
def GetPos(j):
return lambda x: x[j]
for i in range(dim):
newPos_ = FEField("newPos_", inputMesh, LagrangeSpaceP2, numbering)
newPos_.data = res.nodes[:, i]
FillFEField(newPos_, [(ElementFilter(), GetPos(i))])
res.originalIDNodes = np.zeros(numbering.size, dtype=MuscatIndex) - 1
newIds = np.empty(inputMesh.GetNumberOfNodes(), dtype=MuscatIndex)
for i in range(inputMesh.GetNumberOfNodes()):
dof = numbering.GetDofOfPoint(i)
newIds[i] = dof
res.originalIDNodes[newIds] = range(inputMesh.GetNumberOfNodes())
for tag in inputMesh.nodesTags:
res.nodesTags.CreateTag(tag.name).SetIds(newIds[tag.GetIds()])
import Muscat.MeshContainers.ElementsDescription as ED
toQuad = {}
toQuad[ED.Point_1] = ED.Point_1
toQuad[ED.Bar_2] = ED.Bar_3
toQuad[ED.Bar_3] = ED.Bar_3
toQuad[ED.Triangle_3] = ED.Triangle_6
toQuad[ED.Triangle_6] = ED.Triangle_6
toQuad[ED.Tetrahedron_4] = ED.Tetrahedron_10
toQuad[ED.Tetrahedron_10] = ED.Tetrahedron_10
toQuad[ED.Quadrangle_4] = ED.Quadrangle_9
toQuad[ED.Quadrangle_8] = ED.Quadrangle_9
toQuad[ED.Quadrangle_9] = ED.Quadrangle_9
toQuad[ED.Hexahedron_8] = ED.Hexahedron_27
toQuad[ED.Hexahedron_20] = ED.Hexahedron_27
toQuad[ED.Hexahedron_27] = ED.Hexahedron_27
toQuad[ED.Wedge_6] = ED.Wedge_18
toQuad[ED.Wedge_15] = ED.Wedge_18
toQuad[ED.Pyramid_5] = ED.Pyramid_14
toQuad[ED.Pyramid_13] = ED.Pyramid_14
toQuad[ED.Pyramid_14] = ED.Pyramid_14
for name, data in inputMesh.elements.items():
if data.GetNumberOfElements() == 0:
continue
elements = res.elements.GetElementsOfType(toQuad[name])
elements.AddNewElements(numbering[name], data.originalIds)
for tag in data.tags:
elements.tags.CreateTag(tag.name, False).SetIds(tag.GetIds())
res.PrepareForOutput()
return res
[docs]
def QuadToLin(inputMesh: Mesh, divideQuadElements: bool = True, linearizedMiddlePoints: bool = False) -> Mesh:
"""Convert a quadratic mesh to a linear mesh
Nodes fields and element fields are lost.
Parameters
----------
inputMesh : Mesh
the input mesh
divideQuadElements : bool, optional
if the quadratic element must be divided, by default True
linearizedMiddlePoints : bool, optional
if the middle point of the quadratic element must line in the segment formed by the extreme values , by default False
Returns
-------
Mesh
a linear mesh
"""
res = Mesh()
res.CopyProperties(inputMesh)
res.nodes = np.copy(inputMesh.GetPosOfNodes())
res.originalIDNodes = np.arange(
0, res.GetNumberOfNodes(), dtype=MuscatIndex)
res.nodesTags = inputMesh.nodesTags.Copy()
dataBaseOneElement = {}
dataBaseOneElement[ED.Tetrahedron_10] = (ED.Tetrahedron_4, [[0, 1, 2, 3]])
dataBaseOneElement[ED.Triangle_6] = (ED.Triangle_3, [[0, 1, 2]])
dataBaseOneElement[ED.Quadrangle_8] = (ED.Quadrangle_4, [[0, 1, 2, 3]])
dataBaseOneElement[ED.Hexahedron_20] = (
ED.Hexahedron_8, [[0, 1, 2, 3, 4, 5, 6, 7]])
dataBaseOneElement[ED.Bar_3] = (ED.Bar_2, [[0, 1]])
dataBaseSubDived = {}
dataBaseSubDived[ED.Tetrahedron_10] = (ED.Tetrahedron_4, [[0, 4, 6, 7], [1, 5, 4, 8], [
2, 6, 5, 9], [7, 8, 9, 3], [4, 5, 6, 7], [4, 5, 7, 8], [5, 6, 7, 9], [5, 7, 8, 9]])
dataBaseSubDived[ED.Triangle_6] = (
ED.Triangle_3, [[0, 3, 5], [1, 4, 3], [2, 5, 4], [3, 4, 5]])
dataBaseSubDived[ED.Quadrangle_8] = (ED.Quadrangle_4, [[0, 1, 2, 3]])
dataBaseSubDived[ED.Hexahedron_20] = (
ED.Hexahedron_8, [[0, 1, 2, 3, 4, 5, 6, 7]])
dataBaseSubDived[ED.Bar_3] = (ED.Bar_2, [[0, 2], [2, 1]])
# the average of the first 2 columns is affected to the last
linearizedMiddlePoints = {}
linearizedMiddlePoints[ED.Tetrahedron_10] = np.array(
[[0, 1, 4], [1, 2, 5], [2, 0, 6], [0, 3, 7], [1, 3, 8], [2, 3, 9]])
linearizedMiddlePoints[ED.Triangle_6] = np.array(
[[0, 1, 3], [1, 2, 4], [2, 0, 5]])
linearizedMiddlePoints[ED.Quadrangle_8] = np.empty((0, 3))
linearizedMiddlePoints[ED.Hexahedron_20] = np.array([[0, 1, 8], [1, 2, 9], [2, 3, 10], [3, 0, 11], [
4, 5, 12], [5, 6, 13], [6, 7, 14], [7, 4, 15], [0, 4, 16], [1, 5, 17], [2, 6, 18], [3, 7, 19]])
linearizedMiddlePoints[ED.Bar_3] = np.array([[0, 1, 2]])
for elementName in inputMesh.elements.keys():
quadElement = inputMesh.elements[elementName]
if elementName in dataBaseSubDived:
if divideQuadElements:
# new Element Type
nET, extractor = dataBaseSubDived[elementName]
else:
# new Element Type
nET, extractor = dataBaseOneElement[elementName]
lmp = linearizedMiddlePoints[elementName]
linearElements = res.GetElementsOfType(nET)
initNbElem = linearElements.GetNumberOfElements()
nbOfNewElements = len(extractor)
linearElements.Reserve(
initNbElem + quadElement.GetNumberOfElements() * nbOfNewElements)
linearElements.cpt = initNbElem + quadElement.GetNumberOfElements() * \
nbOfNewElements
for i in range(quadElement.GetNumberOfElements()):
quadConn = quadElement.connectivity[i, :]
for j in range(nbOfNewElements):
linearElements.connectivity[initNbElem + i *
nbOfNewElements + j, :] = quadConn[extractor[j]]
if linearizedMiddlePoints:
for k in range(lmp.shape[0]):
res.nodes[quadConn[lmp[k, 2]], :] = (
res.nodes[quadConn[lmp[k, 0]], :] + res.nodes[quadConn[lmp[k, 1]], :]) / 2
elif ED.linear[elementName]:
linearElements = res.GetElementsOfType(elementName)
initNbElem = linearElements.GetNumberOfElements()
linearElements.Reserve(
initNbElem + quadElement.GetNumberOfElements())
nbOfNewElements = 1
linearElements.connectivity[initNbElem: initNbElem +
quadElement.GetNumberOfElements(), :] = quadElement.connectivity
linearElements.cpt = initNbElem + quadElement.GetNumberOfElements()
else:
raise Exception(
f"Error : Quadratic to linear QuadToLin not available for this type of element ({elementName}) ") # pragma: no cover
for originalTag in quadElement.tags:
destinationTag = linearElements.GetTag(originalTag.name)
ids = originalTag.GetIds()
for i in range(originalTag.cpt):
for t in range(nbOfNewElements):
destinationTag.AddToTag(
initNbElem + ids[i] * nbOfNewElements + t)
destinationTag.Tighten()
res.ComputeGlobalOffset()
res.PrepareForOutput()
if divideQuadElements == False:
CleanLonelyNodes(res)
return res
[docs]
def MirrorMesh(inmesh: Mesh, x: Optional[MuscatFloat] = None, y: Optional[MuscatFloat] = None, z: Optional[MuscatFloat] = None, propagateToNodeFields: bool=False) -> Mesh:
"""Compute a new symmetric mesh from the inmesh
All element of inMesh are copied to the output
The user is responsible to call RemoveDoubles to eliminate the possible double nodes on the symmetry plane
Parameters
----------
inmesh : Mesh
input mesh
x : Optional[MuscatFloat], optional
the position of the yz plane to build the symmetric part, if Non, no symmetry to the plane yz
, by default None
y : Optional[MuscatFloat], optional
the position of the xz plane to build the symmetric part, if Non, no symmetry to the plane xz
, by default None
z : Optional[MuscatFloat], optional
the position of the xy plane to build the symmetric part, if Non, no symmetry to the plane xy
_description_, by default None
propagateToNodeFields : bool
True if node fields should be symmetrized together with the mesh
Returns
-------
Mesh
the output mesh
"""
nbPoints = inmesh.GetNumberOfNodes()
outMesh = Mesh()
with outMesh.WithModification(verifyIntegrity=False):
outMesh.CopyProperties(inmesh)
d = 0
if x is not None:
d += 1
if y is not None:
d += 1
if z is not None:
d += 1
outMesh.nodes = np.empty(
(nbPoints * (2**d), inmesh.GetPointsDimensionality()), dtype=MuscatFloat)
outMesh.originalIDNodes = np.empty(
(nbPoints * (2**d),), dtype=MuscatIndex)
if propagateToNodeFields:
for fieldName in inmesh.nodeFields:
outMeshField = np.empty( (nbPoints * (2**d),), dtype=MuscatFloat)
outMeshField[0:nbPoints] = inmesh.nodeFields[fieldName]
outMesh.nodeFields[fieldName] = outMeshField
# copy of points:
outMesh.nodes[0:nbPoints, :] = inmesh.nodes
# copy of points:
outMesh.originalIDNodes[0:nbPoints] = inmesh.originalIDNodes
outMesh.nodesTags = inmesh.nodesTags.Copy()
cpt = nbPoints
def increaseTags(tags, oldSize):
for tag in tags:
ids = tag.GetIds()[:] # make a copy
tag.SetIds(np.hstack((ids, ids + oldSize)))
def symmetrizeNodeFields(outMesh, cpt):
if propagateToNodeFields:
for fieldName in outMesh.nodeFields:
outMesh.nodeFields[fieldName][cpt: (2 * cpt)] = outMesh.nodeFields[fieldName][0:cpt]
if x is not None:
vec = np.array(
[
[-1, 1, 1],
],
dtype=float,
)
outMesh.nodes[cpt: (
2 * cpt), :] = (outMesh.nodes[0:cpt, :] - [x, 0, 0]) * vec + [x, 0, 0]
outMesh.originalIDNodes[cpt: (
2 * cpt)] = outMesh.originalIDNodes[0:cpt]
increaseTags(outMesh.nodesTags, cpt)
symmetrizeNodeFields(outMesh, cpt)
cpt = cpt * 2
if y is not None:
vec = np.array(
[
[1, -1, 1],
],
dtype=float,
)
outMesh.nodes[cpt: (
2 * cpt), :] = (outMesh.nodes[0:cpt, :] - [0, y, 0]) * vec + [0, y, 0]
outMesh.originalIDNodes[cpt: (
2 * cpt)] = outMesh.originalIDNodes[0:cpt]
increaseTags(outMesh.nodesTags, cpt)
symmetrizeNodeFields(outMesh, cpt)
cpt = cpt * 2
if z is not None:
vec = np.array(
[
[1, 1, -1],
],
dtype=float,
)
outMesh.nodes[cpt: (
2 * cpt), :] = (outMesh.nodes[0:cpt, :] - [0, 0, z]) * vec + [0, 0, z]
outMesh.originalIDNodes[cpt: (
2 * cpt)] = outMesh.originalIDNodes[0:cpt]
increaseTags(outMesh.nodesTags, cpt)
symmetrizeNodeFields(outMesh, cpt)
cpt = cpt * 2
for name, vals in inmesh.elements.items():
nbElements = vals.GetNumberOfElements()
outElements = outMesh.GetElementsOfType(name)
outElements.Reserve(nbElements * (2**d))
outElements.connectivity[0:nbElements, :] = vals.connectivity
outElements.tags = copy.deepcopy(vals.tags)
cpt = nbElements
pointCpt = nbPoints
permutation = ED.mirrorPermutation[name]
if x is not None:
outElements.connectivity[cpt: (
2 * cpt), :] = (outElements.connectivity[0:cpt, :] + pointCpt)[:, permutation]
pointCpt = pointCpt * 2
increaseTags(outElements.tags, cpt)
cpt = cpt * 2
if y is not None:
outElements.connectivity[cpt: (
2 * cpt), :] = (outElements.connectivity[0:cpt, :] + pointCpt)[:, permutation]
pointCpt = pointCpt * 2
increaseTags(outElements.tags, cpt)
cpt = cpt * 2
if z is not None:
outElements.connectivity[cpt: (
2 * cpt), :] = (outElements.connectivity[0:cpt, :] + pointCpt)[:, permutation]
pointCpt = pointCpt * 2
increaseTags(outElements.tags, cpt)
cpt = cpt * 2
outElements.cpt = cpt
return outMesh
[docs]
def Create0DElementContainerForEveryPoint(mesh: Mesh) -> ElementsContainer:
"""Create a ElementContainer with 0D elements based on the point of the mesh
The user is responsible to put this ElementContainer back to a mesh (using merge for example)
the Nodal tags are transferred to the new 0D elements
the originalIds is constructed based on the nodes (range(nb Nodes))
Parameters
----------
mesh : Mesh
the input mesh
Returns
-------
ElementsContainer
An ElementsContainer with the new 0D elements
"""
nbNodes = mesh.GetNumberOfNodes()
elements1D = ElementsContainer(ED.Point_1)
elements1D.tags = mesh.nodesTags.Copy()
elements1D.connectivity = np.arange(nbNodes)
elements1D.connectivity.shape = (nbNodes, 1)
elements1D.originalIds = np.arange(nbNodes)
elements1D.cpt = mesh.GetNumberOfNodes()
return elements1D
[docs]
def SubDivideMesh(mesh: Mesh, level: MuscatIndex = 1) -> Mesh:
"""Subdivide a mesh
Parameters
----------
mesh : Mesh
The input mesh
level : MuscatIndex, optional
Number of times the mesh must be divided, by default 1
Returns
-------
Mesh
A new mesh divided level times
"""
if level == 0:
return mesh
res = Mesh()
subdivisionAlmanac = {}
subdivisionAlmanac[ED.Point_1] = [(ED.Point_1, [0])]
subdivisionAlmanac[ED.Bar_2] = [(ED.Bar_2, [0, 2]), (ED.Bar_2, [2, 1])]
subdivisionAlmanac[ED.Quadrangle_4] = [(ED.Quadrangle_4, [0, 4, 8, 7]), (ED.Quadrangle_4, [
4, 1, 5, 8]), (ED.Quadrangle_4, [8, 5, 2, 6]), (ED.Quadrangle_4, [7, 8, 6, 3])]
subdivisionAlmanac[ED.Hexahedron_8] = [
(ED.Hexahedron_8, [0, 8, 24, 11, 16, 22, 26, 20]),
(ED.Hexahedron_8, [8, 1, 9, 24, 22, 17, 21, 26]),
(ED.Hexahedron_8, [24, 9, 2, 10, 26, 21, 18, 23]),
(ED.Hexahedron_8, [11, 24, 10, 3, 20, 26, 23, 19]),
(ED.Hexahedron_8, [16, 22, 26, 20, 4, 12, 25, 15]),
(ED.Hexahedron_8, [22, 17, 21, 26, 12, 5, 13, 25]),
(ED.Hexahedron_8, [26, 21, 18, 23, 25, 13, 6, 14]),
(ED.Hexahedron_8, [20, 26, 23, 19, 15, 25, 14, 7]),
]
subdivisionAlmanac[ED.Triangle_3] = [(ED.Triangle_3, [0, 3, 5]), (ED.Triangle_3, [
3, 1, 4]), (ED.Triangle_3, [3, 4, 5]), (ED.Triangle_3, [5, 4, 2])]
subdivisionAlmanac[ED.Tetrahedron_4] = [
(ED.Tetrahedron_4, [0, 4, 6, 7]),
(ED.Tetrahedron_4, [1, 5, 4, 8]),
(ED.Tetrahedron_4, [2, 6, 5, 9]),
(ED.Tetrahedron_4, [7, 8, 9, 3]),
(ED.Tetrahedron_4, [5, 6, 7, 9]),
(ED.Tetrahedron_4, [5, 6, 4, 7]),
(ED.Tetrahedron_4, [5, 7, 4, 8]),
(ED.Tetrahedron_4, [5, 9, 7, 8]),
]
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo, LagrangeSpaceP2
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.IntegrationRules import NodalEvaluationP2Quadrature
numberingGeo = ComputeDofNumbering(
mesh, LagrangeSpaceGeo, fromConnectivity=True)
numberingP2 = ComputeDofNumbering(mesh, LagrangeSpaceP2)
# Generation of nodes
res.nodes = np.empty((numberingP2.size, mesh.nodes.shape[1]), dtype=MuscatFloat)
res.originalIDNodes = np.zeros(res.nodes.shape[0], dtype=MuscatIndex) - 1
oldToNewDofs = np.zeros(mesh.GetNumberOfNodes(), dtype=MuscatIndex)
oldToNewDofs[numberingP2.doftopointLeft] = numberingP2.doftopointRight
res.originalIDNodes[oldToNewDofs] = mesh.originalIDNodes
for tag in mesh.nodesTags.keys():
name = mesh.nodesTags[tag].name
ids = mesh.nodesTags[tag].GetIds()
res.nodesTags.CreateTag(name).SetIds(oldToNewDofs[ids])
for elemType, data in mesh.elements.items():
spaceGeo = LagrangeSpaceGeo[elemType]
spaceGeo.Create()
elementaryQuadrature = NodalEvaluationP2Quadrature[elemType]
sGeoAtIp = spaceGeo.GetSpaceEvaluatedAt(
elementaryQuadrature.points, elementaryQuadrature.weights)
nGeo = numberingGeo[elemType]
nP2 = numberingP2[elemType]
# generation of nodes
for sf in range(elementaryQuadrature.GetNumberOfPoints()):
geoNs = sGeoAtIp.valN[sf]
for c in range(mesh.nodes.shape[1]):
res.nodes[nP2[:, sf], c] = np.sum(mesh.nodes[:, c][nGeo] * geoNs, axis=1)
# generation of elements
for t, nn in subdivisionAlmanac[elemType]:
newElements = res.GetElementsOfType(t)
offset = newElements.GetNumberOfElements()
tne = newElements.AddNewElements(nP2[:, nn], data.originalIds)
for tag in data.tags.keys():
name = data.tags[tag].name
ids = data.tags[tag].GetIds()
if len(ids) == 0:
continue
newElements.tags.CreateTag(name, False).AddToTag(ids + offset)
mesh.PrepareForOutput()
return SubDivideMesh(res, level - 1)
[docs]
def CreateEdgeMesh(mesh: Mesh, ef: Optional[ElementFilter] = None) -> Mesh:
"""Create a new mesh from the incoming mesh composed only by the edges off all the elements (Bar_2)
all element tags are lost.
Parameters
----------
mesh : Mesh
input mesh
ef : Optional[ElementFilter], optional
element Filter to select the elements to work on, by default None
Returns
-------
Mesh
A mesh composed only with Bar_2
Raises
------
RuntimeError
In the case the faces of the current elements are not Bar_2
"""
if ef is None:
ef = ElementFilter(dimensionality=mesh.GetElementsDimensionality())
res = mesh.View()
res.elements = AllElements()
edges = {}
nex = mesh.GetNumberOfElements()
elements = res.elements.GetElementsOfType(ED.Bar_2)
for selection in ef(mesh):
dim = ED.dimensionality[selection.elementType]
if dim == 3:
faces = ED.faces2[selection.elementType]
elif dim == 2:
faces = ED.faces1[selection.elementType]
else:
faces = [(ED.Bar_2, [0, 1])]
for faceType, localFaceConnectivity in faces:
if faceType != ED.Bar_2: # pragma: no cover
raise RuntimeError(
f"can only treat faces of type ED.Bar_2 nor ({faceType})")
globalFaceConnectivity = selection.elements.connectivity[
selection.indices, :][:, localFaceConnectivity]
key = np.array([nex**x for x in range(ED.numberOfNodes[faceType])])
globalFaceConnectivitySorted = np.sort(
globalFaceConnectivity, axis=1)
unique_bars = np.unique(globalFaceConnectivitySorted, axis=0)
elements.AddNewElements(unique_bars)
unique_bars_final = np.unique(elements.connectivity, axis=0)
elements.connectivity = unique_bars_final
elements.cpt = elements.connectivity.shape[0]
elements.originalIds = np.arange(
elements.connectivity.shape[0], dtype=MuscatIndex)
return res
# ------------------------- CheckIntegrity ------------------------
[docs]
def CheckIntegrity_CreateEdgeMesh(GUI: bool = False) -> str:
"""CheckIntegrity_CreateEdgeMesh"""
a = CreateUniformMeshOfBars(startPoint=[0, 0, 0], stopPoint=[0, 0, 0])
b = CreateEdgeMesh(a)
aa = CreateDisk()
bb = CreateEdgeMesh(aa)
aaa = CreateCube()
bbb = CreateEdgeMesh(aaa)
if GUI: # pragma: no cover
from Muscat.Actions.OpenInParaView import OpenInParaView
OpenInParaView(b)
return "ok"
[docs]
def CheckIntegrity_CreateMeshFromCellsDict(GUI: bool = False) -> str:
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 0], [0.5, 1.5, 0.5]]
cellsDict = {ED.Triangle_3: [[2, 4, 3]], ED.Quadrangle_4: [[0, 1, 2, 3]]}
pointFields = {"my_pt_field": np.random.randn(len(points))}
cellFields = {"my_cell_field": np.asarray([1, 2])}
res = CreateMeshFromCellsDict(
points, cellsDict, pointFields=pointFields, cellFields=cellFields, cellFieldInDictOrder=False)
np.testing.assert_equal(res.elemFields["my_cell_field"], [1, 2])
res.VerifyIntegrity()
assert len(points) == res.GetNumberOfNodes()
res = CreateMeshFromCellsDict(
points, cellsDict, pointFields=pointFields, cellFields=cellFields, cellFieldInDictOrder=True)
print(res.elemFields["my_cell_field"])
np.testing.assert_equal(res.elemFields["my_cell_field"], [2, 1])
res.VerifyIntegrity()
print(res)
return "OK"
[docs]
def CheckIntegrity_CreateDisk(GUI: bool == False) -> str:
"""CheckIntegrity_CreateDisk"""
a = CreateDisk()
return "ok"
[docs]
def CheckIntegrity_SubDivideMesh(GUI: bool == False) -> str:
points = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
[0, 0, 1], [1, 0, 1.5], [1, 1, 1], [0, 1, 1.5]]
hexahedrons = [
[0, 1, 2, 3, 4, 5, 6, 7],
]
mesh = CreateMeshOf(points, hexahedrons, ED.Hexahedron_8)
mesh.nodesTags.CreateTag("FirstPoint").AddToTag(0)
mesh.GetElementsOfType(ED.Hexahedron_8)
mesh.GetElementsOfType(ED.Hexahedron_8).tags
mesh.GetElementsOfType(ED.Hexahedron_8).tags.CreateTag("OnlyHex")
mesh.GetElementsOfType(ED.Hexahedron_8).tags.CreateTag(
"OnlyHex", False).AddToTag(0)
mesh.GetElementsOfType(ED.Hexahedron_8).tags.CreateTag("OnlyHexEmpty")
outMesh = SubDivideMesh(mesh, 1)
print(mesh)
print(outMesh)
if GUI: # pragma: no cover
from Muscat.Bridges.vtkBridge import PlotMesh
PlotMesh(outMesh)
from Muscat.IO.XdmfWriter import WriteMeshToXdmf
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tempdir = TemporaryDirectory.GetTempPath()
WriteMeshToXdmf(tempdir + "CheckIntegrity_SubDivideMesh.xdmf", outMesh,
PointFields=[outMesh.originalIDNodes], PointFieldsNames=["originalIDNodes"])
print(tempdir)
np.testing.assert_equal(outMesh.GetNumberOfNodes(), 27)
if outMesh.GetNumberOfElements() != 8:
raise # pragma: no cover
return "ok"
[docs]
def CheckIntegrity_Create0DElementContainerForEveryPoint(GUI: bool == False) -> str:
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tets = [
[0, 1, 2, 3],
]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("FirstPoint").AddToTag(0)
print(mesh)
zeroDElements = Create0DElementContainerForEveryPoint(mesh)
print(zeroDElements)
np.testing.assert_equal(
zeroDElements.GetNumberOfElements(), mesh.GetNumberOfNodes())
return "ok"
[docs]
def CheckIntegrity_MirrorMesh(GUI: bool == False) -> str:
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tets = [
[0, 1, 2, 3],
]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
mesh.nodesTags.CreateTag("FirstPoint").AddToTag(0)
mesh.GetElementsOfType(ED.Tetrahedron_4).tags.CreateTag(
"OnlyTet").AddToTag(0)
outMesh = MirrorMesh(mesh, x=0)
np.testing.assert_equal(outMesh.GetNumberOfNodes(), 8)
np.testing.assert_equal(outMesh.GetNumberOfElements(), 2)
outMesh = MirrorMesh(mesh, y=0, z=0)
np.testing.assert_equal(outMesh.GetNumberOfNodes(), 16)
np.testing.assert_equal(outMesh.GetNumberOfElements(), 4)
mesh.nodeFields["nodeField"] = np.arange(mesh.GetNumberOfNodes(), dtype=MuscatFloat)
outMesh = MirrorMesh(mesh, x=0, y=0, z=0, propagateToNodeFields=True)
np.testing.assert_equal(len(outMesh.nodeFields["nodeField"]), 32)
return "ok"
[docs]
def CheckIntegrity_QuadToLin(GUI: bool == False) -> str:
myMesh = Mesh()
myMesh.nodes = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [0.5, 0, 0], [
0.5, 0.5, 0], [0, 0.5, 0], [0, 0, 0.5], [0.5, 0, 0.5], [0, 0.5, 0.5]], dtype=float)
tag = myMesh.GetNodalTag("linPoints")
tag.AddToTag(0)
tag.AddToTag(1)
tag.AddToTag(2)
tag.AddToTag(3)
import Muscat.MeshContainers.ElementsDescription as ED
elements = myMesh.GetElementsOfType(ED.Tetrahedron_10)
elements.AddNewElement([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 0)
elements = myMesh.GetElementsOfType(ED.Triangle_6)
elements.AddNewElement([0, 1, 2, 4, 5, 6], 1)
elements = myMesh.GetElementsOfType(ED.Bar_3)
elements.AddNewElement([0, 1, 4], 2)
elements = myMesh.GetElementsOfType(ED.Bar_2)
elements.AddNewElement([0, 1], 3)
myMesh.AddElementToTagUsingOriginalId(3, "LinElements")
print(myMesh)
linMesh = QuadToLin(myMesh, divideQuadElements=False)
print(linMesh)
print(QuadToLin(myMesh, divideQuadElements=True))
print(QuadToLin(myMesh, divideQuadElements=True, linearizedMiddlePoints=True))
from Muscat.MeshTools.MeshFieldOperations import QuadFieldToLinField
QuadFieldToLinField(myMesh, np.arange(myMesh.GetNumberOfNodes()))
QuadFieldToLinField(myMesh, np.arange(myMesh.GetNumberOfNodes()), linMesh)
return "ok"
[docs]
def CheckIntegrity_CreateMeshOfTriangles(GUI: bool == False) -> str:
res = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [
0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
print(res)
return "OK"
[docs]
def CheckIntegrity_CreateCube(GUI: bool = False) -> str:
mesh = CreateCube(dimensions=[20, 21, 22], spacing=[
2.0, 2.0, 2.0], ofTetras=False)
mesh = CreateCube(dimensions=[20, 21, 22], spacing=[
2.0, 2.0, 2.0], ofTetras=True)
return "ok"
[docs]
def CheckIntegrity_CreateSquare(GUI: bool = False) -> str:
mesh = CreateSquare(dimensions=[20, 21], spacing=[
2.0, 2.0], ofTriangles=False)
mesh = CreateSquare(dimensions=[20, 21], spacing=[
2.0, 2.0], ofTriangles=True)
return "ok"
[docs]
def CheckIntegrity_ToQuadraticMesh(GUI: bool = False) -> str:
inMesh = CreateCube()
inMesh.GetElementsOfType(ED.Triangle_3) # add empty element container
outMesh = ToQuadraticMesh(inMesh)
return "ok"
[docs]
def CheckIntegrity_MeshToSimplex(GUI: bool = False) -> str:
inMesh = CreateCube()
inMesh.GetElementsOfType(ED.Hexahedron_8).GetTag(
"First Hexahedron_8").SetIds([0])
MeshToSimplex(inMesh)
inMesh = CreateSquare()
inMesh.GetElementsOfType(ED.Quadrangle_4).GetTag(
"First Quadrangle_4").SetIds([0])
inMesh.GetElementsOfType(ED.Triangle_3)
inMesh.GetElementsOfType(ED.Hexahedron_8)
MeshToSimplex(inMesh)
MeshToSimplex(inMesh, inPlace=False)
return "ok"
[docs]
def CheckIntegrity(GUI: bool = False) -> str:
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
toTest = [
CheckIntegrity_CreateEdgeMesh,
CheckIntegrity_ToQuadraticMesh,
CheckIntegrity_CreateUniformMeshOfBars,
CheckIntegrity_CreateCube,
CheckIntegrity_CreateSquare,
CheckIntegrity_CreateMeshOfTriangles,
CheckIntegrity_MeshToSimplex,
CheckIntegrity_QuadToLin,
CheckIntegrity_MirrorMesh,
CheckIntegrity_Create0DElementContainerForEveryPoint,
CheckIntegrity_SubDivideMesh,
CheckIntegrity_CreateDisk,
CheckIntegrity_CreateMeshFromCellsDict,
]
return RunListOfCheckIntegrities(toTest, GUI=GUI)
if __name__ == "__main__":
print(CheckIntegrity(False)) # pragma: no cover