Source code for Muscat.Containers.MeshTools
# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
from typing import Optional
import numpy as np
from Muscat.Containers.Mesh import Mesh
from Muscat.Types import ArrayLike
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.ElementsContainers import ElementContainerLike
[docs]def IsClose(mesh1, mesh2) -> bool:
"""Verified if two meshes are close (in the sense of numpy.isclose)
meshes must have the same :
- nodes
- nodes tags
- elements
- element tags
Parameters
----------
mesh1 : _type_
first mesh to be compare
mesh2 : _type_
second mesh to be compare
Returns
-------
bool
True if mesh1 and mesh2 are close enough
"""
if type(mesh1) != type(mesh2):
print("types not equal")
return False
if not np.all(np.isclose(mesh1.nodes, mesh2.nodes)):
print(mesh1.nodes)
print(mesh2.nodes)
print("nodes not equal")
return False
for tag1 in mesh1.nodesTags:
tag2 = mesh2.nodesTags[tag1.name]
if not np.all(np.isclose(tag1.GetIds(), tag2.GetIds())):
print("Nodal tag " + str(tag1.name) + " not equal")
return False
for name, data1 in mesh1.elements.items():
data2 = mesh2.elements[name]
if not np.all(np.isclose(data1.connectivity, data2.connectivity)):
print("Connectivity for " + str(name) + " not equal")
return False
for tag1 in data1.tags:
tag2 = data2.tags[tag1.name]
if not np.all(np.isclose(tag1.GetIds(), tag2.GetIds())):
print("Tag " + str(tag1.name) + " is not equal for element" + str(name))
return False
def CompareFields(fields1, fields2):
for name, data1 in fields1.items():
data2 = fields2[name]
if len(data1) != len(data2):
print("Field " + str(name) + " has different size")
return False
if data1.dtype == data2.dtype and data1.dtype == object:
if not np.all([d1 == d2 for d1, d2 in zip(data1, data2)]):
print("Field " + str(name) + " not equal")
return False
continue
if data1.dtype.type is np.string_ or data1.dtype.char == "U":
if not np.all(data1 == data2):
print("Field " + str(name) + " not equal")
return False
else:
if not np.all(np.isclose(data1, data2)):
print("Field " + str(name) + " not equal")
return False
if CompareFields(mesh1.nodeFields, mesh2.nodeFields) == False:
return False
if CompareFields(mesh1.elemFields, mesh2.elemFields) == False:
return False
return True
[docs]def GetElementsCenters(mesh: Optional[Mesh] = None, nodes: Optional[ArrayLike] = None, elements: Optional[ElementContainerLike] = None, dim: Optional[int] = None) -> np.ndarray:
"""Internal function to compute the element centers.
Waring!!!! This function is used in the filters implementation
no Filter can appear in this implementation
Parameters
----------
mesh : Mesh, optional
if mesh is not none the element center for all the element is calculated, by default None
nodes : ArrayLike, optional
if mesh is non, nodes and elements must be supplied to compute the element center only for
the ElementContainer, by default None
elements : ElementContainerLike, optional
if mesh is non, nodes and elements must be supplied to compute the element center only for
the ElementContainer, by default None
dim : int, optional
the dimensionality (int) to filter element to be treated, by default None
Returns
-------
np.ndarray
the center for each element
"""
#
if mesh is not None and elements is not None:
raise (Exception("Cannot treat mesh and element at the same time"))
def ProcessElements(nod, els):
connectivity = els.connectivity[0 : els.cpt, :]
localRes = np.zeros((els.GetNumberOfElements(), nod.shape[1]))
for i in range(nod.shape[1]):
localRes[:, i] += np.sum(nod[connectivity, i], axis=1)
localRes /= connectivity.shape[1]
return localRes
if mesh is not None:
pos = mesh.GetPosOfNodes()
res = np.empty((mesh.GetNumberOfElements(dim=dim), pos.shape[1]))
cpt = 0
from Muscat.Containers.Filters.FilterObjects import ElementFilter
ff = ElementFilter(dimensionality=dim)
for selection in ff(mesh):
res[cpt : cpt + selection.elements.GetNumberOfElements()] = ProcessElements(mesh.nodes, selection.elements)
cpt += selection.elements.GetNumberOfElements()
else:
return ProcessElements(nodes, elements)
return res
[docs]def TagsToRefs(mesh: Mesh, offset: int = 1):
"""Function to convert the tags (strings) of the mesh (for nodes and elements) in one reference number (int).
Needed for some mesh format as .mesh
Parameters
----------
mesh : a mesh with tags
offset : int greater than one
Returns
-------
refByNode : a np.array giving at position i the reference number (int) for the node i in the mesh
dictNodesRefsToTags : a dictionnary giving for each reference number the corresponding node tags
refByElement : a np.array giving at position i the reference number (int) for the element i in the mesh
dictElemRefsToTags : a dictionnary giving for each reference number the corresponding element tags
"""
nodesCorrTab = np.zeros((mesh.GetNumberOfNodes(), len(mesh.nodesTags)), dtype = bool)
for i, tag in enumerate(mesh.nodesTags):
nodesCorrTab[tag.GetIds(), i] = True
nodesTags, nodesInverse = np.unique(nodesCorrTab, axis = 0, return_inverse = True)
dictNodesRefsToTags = {}
nodeTagNames = mesh.nodesTags.keys()
refByNode = np.zeros(mesh.GetNumberOfNodes(), dtype = int)
k = offset
for i, nodesTags_i in enumerate(nodesTags):
tags = [nodeTagNames[j] for j in np.where(nodesTags_i == True)[0]]
if tags:
dictNodesRefsToTags[k] = tags
refByNode[nodesInverse == i] = k
k += 1
elementsCorrTab = np.zeros((mesh.GetNumberOfElements(), len(mesh.GetNamesOfElementTags())), dtype = bool)
for i, tag in enumerate(mesh.GetNamesOfElementTags()):
elementsCorrTab[mesh.GetElementsInTag(tag), i] = 1
elementsTags, elementsInverse = np.unique(elementsCorrTab, axis = 0, return_inverse = True)
dictElemRefsToTags = {}
elementTagNames = mesh.GetNamesOfElementTags()
refByElement = np.zeros(mesh.GetNumberOfElements(), dtype = int)
for i, elementsTags_i in enumerate(elementsTags):
tags = [elementTagNames[j] for j in np.where(elementsTags_i == True)[0]]
if tags:
dictElemRefsToTags[k] = tags
refByElement[elementsInverse == i] = k
k += 1
return refByNode, dictNodesRefsToTags, refByElement, dictElemRefsToTags
[docs]def RefsToTags(mesh: Mesh, dictNodesRefsToTags: dict, dictElemRefsToTags: dict, prefix: str = "new"):
"""Function to retrieve the tags (strings) of the mesh (for nodes and elements)
from the elementFields and nodeFields "ref" which are deleted after conversion to Tags.
Needed for some mesh format as .mesh
Parameters
----------
mesh : a mesh with an element field and a node field named "ref"
dictNodesRefsToTags : a dictionnary giving for each reference number the corresponding node tags
dictElemRefsToTags : a dictionnary giving for each reference number the corresponding element tags
prefix : the prefix to add if new reference numbers have been created
Returns
-------
mesh : the mesh with the tags
"""
refByNode = mesh.nodeFields["refs"]
for ref, tags in dictNodesRefsToTags.items():
for tag in tags:
mesh.nodesTags.CreateTag(tag, False).AddToTag(np.where(refByNode == ref)[0])
newRefs = np.setdiff1d(np.unique(refByNode), list(dictNodesRefsToTags.keys()))
if prefix:
for ref in newRefs:
mesh.nodesTags.CreateTag("NTag_{0}_{1}".format(prefix, ref), False).AddToTag(np.where(refByNode == ref)[0])
del mesh.nodeFields["refs"]
refByElement = mesh.elemFields["refs"]
for ref, tags in dictElemRefsToTags.items():
for tag in tags:
mesh.AddElementsToTag(np.where(refByElement == ref)[0], tag)
del mesh.elemFields["refs"]
mesh.nodesTags.RemoveDoubles()
return mesh
[docs]def CheckIntegrity_GetCellCenters():
mesh1 = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
res = GetElementsCenters(mesh1)
print(res)
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
mesh2 = CreateConstantRectilinearMesh(dimensions=[2, 3, 2], spacing=[1, 1, 1], origin=[0, 0, 0])
res = GetElementsCenters(mesh=mesh2)
print(res)
return "ok"
[docs]def CheckIntegrity_TagsToRefs():
from Muscat.IO.GmshReader import ReadGmsh
from Muscat.TestData import GetTestDataPath
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
tempdir = TemporaryDirectory.GetTempPath()
filename = GetTestDataPath() + "dent3D.msh"
mesh0 = ReadGmsh(filename)
mesh0.ConvertDataForNativeTreatment()
print(mesh0)
refByNode, dictNodesRefsToTags, refByElement, dictElemRefsToTags = TagsToRefs(mesh0)
from Muscat.IO.MeshWriter import WriteMesh
WriteMesh(tempdir + "TagsToRefs", mesh0, binary=True, nodalRefNumber=refByNode, elemRefNumber=refByElement)
from Muscat.IO.MeshReader import ReadMesh
nmesh = ReadMesh(tempdir + "TagsToRefs", ReadRefsAsField=True)
print(nmesh)
mesh1 = RefsToTags(nmesh, dictNodesRefsToTags, dictElemRefsToTags)
print(mesh1)
if mesh0 == mesh1:
return "OK"
if __name__ == "__main__":
print(CheckIntegrity()) # pragma: no cover