# -*- 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 Dict, Union, List, Any, Tuple, Optional
import bisect
import copy
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
from Muscat.Helpers.Decorators import froze_it
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.ElementsContainers import ElementsContainer, AllElements
from Muscat.Containers.Tags import Tag, Tags
from Muscat.LinAlg.MatVecOperations import ImmutableView
[docs]@froze_it
class Mesh():
"""
Class to store an unstructured (i.e. general) mesh:
* self.nodes : the nodes positions
* self.originalIdNodes : the ids of the previous mesh/file
* self.elements : all the element in the mesh (of type AllElements)
The manual construction of this class must always end with a call to the
function mesh.PrepareForOutput() or with the context manager mesh.WithModification().
This will ensure the correctness of the internal data
"""
def __init__(self) -> None:
super().__init__()
self.nodes: np.ndarray = np.empty((0, 3), dtype=MuscatFloat)
"""Nodes of the mesh (this are point in 2D or 3D always of format MuscatFloat)"""
self.nodesTags = Tags()
"""Tags for the nodes """
self.nodeFields: Dict[str, np.ndarray] = {}
""" This is a dictionary containing fields defined over the nodes (all the nodes).
The keys are the names of the fields
the values are the actual data of size (nb nodes, nb of components)"""
self.elemFields: Dict[str, np.ndarray] = {}
""" This is a dictionary containing fields defined over the elements (all the elements).
The keys are the names of the fields
the values are the actual data of size (nb elements, nb of components)"""
self.props: Dict[str, Any] = {}
"""Metadata this is just a dictionary that can be used to transport
information with the mesh, please use the class name as key of the
object using/setting the information.
Special keys are:
- "IsConstantRectilinear" : True if the mesh is constant Rectilinear
- "spacing" : a np.array with the spacing in each direction
- "origin" : a np.array with the origin of the mesh
- "dimensions" : a np.array with the number of cells in each direction
"""
self.originalIDNodes: np.ndarray = np.empty((0,), dtype=MuscatIndex)
"""The role of the originalIDNodes is double.
In the case the origin of the mesh is a file from disk, it is used to store the "node number" (a concept present in many file formats).
In the case the mesh is a transformation of an initial mesh (extraction for example) the originalIDNodes is used to track the origin (the index) of each entity.
"""
self.elements: AllElements = AllElements()
"""Attribute to store for each element type a instance of a element container (:py:class:`Muscat.Containers.ElementsContainers.ElementsContainer` or :py:class:`Muscat.Containers.ElementsContainers.StructuredElementsContainer` ) containing the connectivity, originalIds and the tags.
"""
[docs] def View(self) -> Mesh:
"""return a copy of this object recursively, numpy arrays are view to the original ones
but with the flags.writeable = False
Returns
-------
Mesh
a view of Mesh
"""
res = Mesh()
res.nodes = ImmutableView(self.nodes)
res.originalIDNodes = ImmutableView(self.originalIDNodes)
res.nodesTags = self.nodesTags.View()
res.elements = self.elements.View()
for k, v in self.nodeFields.items():
res.nodeFields[k] = ImmutableView(v)
for k, v in self.elemFields.items():
res.elemFields[k] = ImmutableView(v)
res.props = self.props.copy()
return res
def __eq__(self, other: Union[Mesh, Any]) -> bool:
if not isinstance(other, Mesh):
return False
if self is other:
return True
if self.nodesTags != other.nodesTags:
return False
if len(self.nodeFields) != len(other.nodeFields):
return False
nodeFieldsNames = set()
nodeFieldsNames.update(self.nodeFields.keys())
nodeFieldsNames.update(other.nodeFields.keys())
for name in nodeFieldsNames:
v1 = self.nodeFields.get(name, None)
v2 = other.nodeFields.get(name, None)
if v1 is None or v2 is None:
return False
if not np.array_equal(v1, v2):
return False
if len(self.elemFields) != len(other.elemFields):
return False
for k, v in self.elemFields.items():
if k in other.elemFields:
if not np.array_equal(v, other.elemFields[k]):
return False
else:
return False
if self.props != other.props:
return False
if not np.array_equal(self.nodes, other.nodes):
return False
if not np.array_equal(self.originalIDNodes, other.originalIDNodes):
return False
if self.elements != other.elements:
return False
return True
[docs] def CopyProperties(self, other: Mesh) -> None:
self.props = copy.deepcopy(other.props)
[docs] def GetElementsOfType(self, typename: Union[ED.ElementType, str]) -> ElementsContainer:
"""return the element container for the element name (typename)
Parameters
----------
typename : str
the name of the elements to extract
Returns
-------
ElementContainer
The Element Container for element
"""
return self.elements.GetElementsOfType(typename)
[docs] def GetNodalTag(self, tagName: str) -> Tag:
"""return the Tag (instance of the class) with name tagName.
If the tag does not exist a new is created """
return self.nodesTags.CreateTag(tagName, False)
[docs] def GetElementsDimensionality(self) -> int:
"""Return the maximal dimension of the elements
Returns
-------
int
the max of all elements dimensionality
"""
return int(np.max([ED.dimensionality[elementType] for elementType in self.elements.keys()]))
[docs] def GenerateManufacturedOriginalIDs(self, nodesOffset: int = 0, elementsOffset: int = 0) -> None:
"""function to generate a valid original ids for the nodes and the elements
Parameters
----------
offset : int, optional
offset to use to start counting nodes and elements.
This is useful for code needing numbering stating from 1, by default 0
"""
self.originalIDNodes = np.arange(self.GetNumberOfNodes(), dtype=MuscatIndex)
self.originalIDNodes += nodesOffset
counter = 0
for key, value in self.elements.items():
value.originalIds = np.arange(counter, counter+value.GetNumberOfElements(), dtype=MuscatIndex)
value.originalIds += elementsOffset
counter += value.GetNumberOfElements()
[docs] def ConvertDataForNativeTreatment(self) -> None:
"""Ensure all the data is compatible with the cpp treatment (continuous in C order)
"""
self.originalIDNodes = np.asarray(self.originalIDNodes, dtype=MuscatIndex, order="C")
self.nodes = np.asarray(self.nodes, dtype=MuscatFloat, order="C")
for data in self.elements.values():
data.ConvertDataForNativeTreatment()
[docs] def WithModification(self) -> ClosingMeshAutomatically:
"""Context manager to release all the extra memory of the mesh after modification.
this context manager call mesh.PrepareForOutput() automatically at the end of the block
Example
-------
with mesh.WithModification():
# mesh.PrepareForOutput() is called at the end of this block,
# even if code in the block raises an exception
"""
class ClosingMeshAutomatically():
def __init__(self, mesh: Mesh) -> None:
self.mesh: Mesh = mesh
def __enter__(self) -> None:
pass
def __exit__(self, type, value, traceback) -> None:
self.mesh.PrepareForOutput()
return ClosingMeshAutomatically(self)
[docs] def GetNumberOfNodes(self) -> MuscatIndex:
"""return the total number of nodes in the mesh
Returns
-------
int
Number of node in the mesh
"""
return self.nodes.shape[0]
[docs] def SetNodes(self, nodes: ArrayLike, originalIDNodes: Optional[ArrayLike] = None, generateOriginalIDs: bool = False) -> None:
"""Set nodes and original Ids in the correct internal format.
Parameters
----------
nodes : ArrayLike
Nodes positions of size (number of node, space dimensionality )
originalIDNodes : ArrayLike, optional
original Ids, by default None
generateOriginalIDs : bool, optional
if True to generate original ids with numpy.arange, by default False
"""
self.nodes = np.require(nodes, dtype=MuscatFloat, requirements=['C', 'A'])
if originalIDNodes is not None:
self.originalIDNodes = np.require(originalIDNodes, dtype=MuscatIndex, requirements=['C', 'A'])
elif generateOriginalIDs:
self.originalIDNodes = np.arange(self.GetNumberOfNodes(), dtype=MuscatIndex)
[docs] def GetPointsDimensionality(self) -> int:
"""Return the number of coordinates of the points
Returns
-------
int
number of columns in the point array
"""
return self.nodes.shape[1]
[docs] def GetNumberOfElements(self, dim: Optional[int] = None) -> MuscatIndex:
"""Compute and return the total number of elements in the mesh
Parameters
----------
dim : int, optional
dimensionality filter, by default None
Returns
-------
int
number of element in the mesh (filtered by dim if dim != None)
"""
numberOfElements:MuscatIndex = 0
if dim == None:
for elementType, data in self.elements.items():
numberOfElements += data.GetNumberOfElements()
else:
for elementType, data in self.elements.items():
if ED.dimensionality[elementType] == dim:
numberOfElements += data.GetNumberOfElements()
return numberOfElements
[docs] def Merge(self, other: Mesh) -> None:
"""Merge the other mesh into this. no treatment is done to remove double
nodes or elements.
Parameters
----------
other : Mesh
the mesh to be merge into this (self)
"""
def MergeFields(dest, origin, dSize, oSize):
names = set()
names.update(dest.keys())
names.update(origin.keys())
for name in names:
dv = dest.get(name, None)
ov = origin.get(name, None)
# if field absent in the origin (other)
if dv is not None and ov is None:
# take the shape of destination
# this can have more than one column
shape = list(dv.shape)
# the size of the temporary empty field
shape[0] = oSize
ov = np.zeros(shape, dtype=dv.dtype)
elif dv is None and ov is not None:
# take the shape of the origin
shape = list(ov.shape)
shape[0] = dSize
dv = np.zeros(shape, dtype=ov.dtype)
shape = list(ov.shape)
if len(shape) == 1:
dest[name] = np.hstack((dv, ov))
else:
dest[name] = np.vstack((dv, ov))
# if name == "ZeroPointField":
# exit(0)
MergeFields(self.nodeFields, other.nodeFields, self.GetNumberOfNodes(), other.GetNumberOfNodes())
MergeFields(self.elemFields, other.elemFields, self.GetNumberOfElements(), other.GetNumberOfElements())
nbNodes = self.GetNumberOfNodes()
self.nodes = np.vstack((self.nodes, other.nodes))
self.originalIDNodes = np.hstack((self.originalIDNodes, -other.originalIDNodes))
for tagName in other.nodesTags.keys():
tag = other.nodesTags[tagName]
self.nodesTags.CreateTag(tagName, False).AddToTag(tag.GetIds()+nbNodes)
self.MergeElements(other, force=True, nodesOffset=nbNodes)
[docs] def MergeElements(self, other: Mesh, force: bool = False, nodesOffset: MuscatIndex = 0) -> None:
"""Merge the elements for a second mesh into this
the nodes array must be the same (not only equal)
the user can force the merge if needed (force argument)
Parameters
----------
other : Mesh
the other mesh
force : bool, optional
True to force the merge even if the nodes are not the same, by default False
"""
if (self.nodes is not other.nodes) and (not force):
raise (RuntimeError("the two meshes does not share the same nodes (potentially dangerous)"))
self.elements.MergeElements(other.elements, nodesOffset)
[docs] def ComputeGlobalOffset(self) -> Dict[ED.ElementType, MuscatIndex]:
""" Return a dict with with key the elementType and the value
the number of elements in the mesh before this type
Recompute the Global Offset,
"""
offsets = dict()
cpt = 0
for data in self.elements.values():
offsets[data.elementType] = cpt
cpt += data.GetNumberOfElements()
return offsets
[docs] def ComputeBoundingBox(self) -> Tuple[np.ndarray, np.ndarray]:
"""Compute the bounding box (min and max) of the mesh
Returns
-------
Tuple[np.ndarray, np.ndarray]
the min, max
"""
return np.amin(self.nodes, axis=0), np.amax(self.nodes, axis=0)
[docs] def AddNodeToTagUsingOriginalId(self, oid: int, tagname: str) -> None:
"""add a node (using the original id ) to a tag (tagname)
Parameters
----------
oid : int
Original id node
tagname : str
Tag name
Raises
------
Exception
if the original id is not found
"""
w = np.where(self.originalIDNodes == oid)
if len(w[0]) > 0:
self.GetNodalTag(tagname).AddToTag(w[0])
else:
raise Exception("No node with id " + str(oid)) # pragma: no cover
[docs] def AddElementToTagUsingOriginalId(self, oid: int, tagname: str) -> None:
"""Add an element (using the originalId) to a tag (tagname)
Parameters
----------
oid : int
the original id of the element
tagname : str
Tag name
Raises
------
Exception
if the element with the original id is not found
"""
for data in self.elements.values():
w = np.where(data.originalIds[:data.cpt] == oid)
if len(w[0]) > 0:
data.tags.CreateTag(tagname, False).AddToTag(w[0])
break
else:
raise Exception("No element with id " + str(oid)) # pragma: no cover
[docs] def AddElementsToTag(self, globalElemNumbers: ArrayLike, tagname: str) -> None:
"""Add elements (using the global element number) to a tag (tagname)
Parameters
----------
globalElemNumbers : ArrayLike
the list of the global ids of the element
tagname : str
Tag name
Raises
------
Exception
if some element indices are greater than the number of elements.
"""
if len(globalElemNumbers) == 0:
return
elementNotTreated = np.unique(globalElemNumbers)
cpt = MuscatIndex(0)
cpt2 = MuscatIndex(0)
oldIndex = 0
for data in self.elements.values():
cpt2 = data.GetNumberOfElements() + cpt
index = bisect.bisect_left(elementNotTreated, cpt2)
dataToTreat = elementNotTreated[oldIndex:index]
if len(dataToTreat):
tag = data.tags.CreateTag(tagname, False)
tag.AddToTag(dataToTreat-cpt)
tag.RemoveDoubles()
cpt += data.GetNumberOfElements()
oldIndex = index
if elementNotTreated[-1] >= cpt2:
raise RuntimeError(f"No element found {elementNotTreated}") # pragma: no cover
[docs] def AddElementToTag(self, globalElemNumber: int, tagname: str) -> None:
"""Add an element (using the global element number) to a tag (tagname)
Parameters
----------
globalElemNumber : int
the element number
tagname : str
Tag name
Raises
------
Exception
if the the element is not found
"""
return self.AddElementsToTag([globalElemNumber], tagname)
[docs] def GetPosOfNodes(self) -> np.ndarray:
"""return the position of all the nodes
Returns
-------
np.ndarray
The position for all the nodes in the mesh
"""
return self.nodes
[docs] def GetElementsOriginalIDs(self, ef=None) -> np.ndarray:
"""return a single list with all the original Ids concatenated
Parameters
----------
ef : ElementFilter, optional
filter to select the information, by default None
Returns
-------
np.ndarray
array with all the originalId for the selected elements
"""
if ef is None:
from Muscat.Containers.Filters.FilterObjects import ElementFilter
ef = ElementFilter()
from Muscat.Containers.Filters.FilterTools import GetFrozenFilter
fef = GetFrozenFilter(ef, self)
res = np.empty(fef.GetElementSelectionSize(), dtype=MuscatIndex)
for selection in fef(self):
res[selection.GetSelectionSlice()] = selection.elements.originalIds[selection.indices]
return res
[docs] def SetElementsOriginalIDs(self, originalIDs: ArrayLike) -> None:
"""Set from a single list all the element original Ids
Parameters
----------
originalIDs : ArrayLike
the list of all the original ids in the order of the mesh
"""
originalIDs = np.asarray(originalIDs, dtype=MuscatIndex)
cpt = 0
for data in self.elements.values():
data.originalIds = originalIDs[cpt:data.GetNumberOfElements()+cpt]
cpt += data.GetNumberOfElements()
[docs] def GetElementsInTag(self, tagname: str, useOriginalId: bool = False) -> np.ndarray:
"""return a list with the ids of the elements in a tag
Parameters
----------
tagname : str
the tag name
useOriginalId : bool, optional
to return the list of original ids and not the global numbers, by default False
Returns
-------
np.ndarray
the list off all the element in the tag named "tagname"
"""
res = np.zeros((self.GetNumberOfElements(),), dtype=MuscatIndex)
cpt = 0
offsets = self.ComputeGlobalOffset()
for name, data in self.elements.items():
if tagname in data.tags:
tag = data.tags[tagname].GetIds()
if useOriginalId:
res[cpt:cpt+len(tag)] = data.originalIds[tag]
else:
res[cpt:cpt+len(tag)] = offsets[name]+tag
cpt += len(tag)
return res[0:cpt]
[docs] def PrepareForOutput(self) -> None:
"""Function to free the extra memory used after the incremental creation of a mesh
and final treatment (offset computation).
"""
self.nodesTags.Tighten()
for data in self.elements.values():
data.tags.RemoveDoubles()
data.Tighten()
self.VerifyIntegrity()
def __str__(self) -> str:
"""Return a string with a summary of the mesh
Returns
-------
str
summary of the mesh
"""
res = "Mesh \n"
res += " Number Of Nodes : {} \n".format(self.GetNumberOfNodes())
res += " Tags : " + " ".join(["("+x.name+":"+str(len(x))+")" for x in self.nodesTags]) + "\n"
res += " Number Of Elements : {} \n".format(self.GetNumberOfElements())
for data in self.elements.values():
res += str(data)
if len(self.nodeFields.keys()) > 0:
res += " nodeFields : " + str(list(self.nodeFields.keys())) + "\n"
if len(self.elemFields.keys()) > 0:
res += " elemFields : " + str(list(self.elemFields.keys())) + "\n"
return res
[docs] def Clean(self) -> None:
"""Remove unnecessary data :
1) empty tags
2) empty element containers
"""
self.nodesTags.RemoveEmptyTags()
for k in self.elements.keys():
if self.elements[k].GetNumberOfElements() == 0:
del self.elements[k]
continue
self.elements[k].tags.RemoveEmptyTags()
[docs] def VerifyIntegrity(self) -> None:
"""Integrity verification of the mesh.
Raises
------
RuntimeError
if the integrity of the mesh is compromised
"""
def CheckType(array, type, name):
if array.dtype != type: # pragma: no cover
raise RuntimeError(f"{name} with the wrong type ={array.dtype} : expected = {type}")
# verification nodes an originalIdNodes are compatible
if len(self.nodes.shape) != 2:
raise RuntimeError(f"Error in the shape of nodes: {self.nodes.shape=}")
if self.nodes.flags['C_CONTIGUOUS'] is False:
raise RuntimeError("Error in the order of nodes")
if len(self.originalIDNodes.shape) != 1:
print(f"{self.originalIDNodes.shape=}")
raise RuntimeError("Error in the shape of originalIDNodes")
if self.originalIDNodes.flags['C_CONTIGUOUS'] is False:
raise RuntimeError("Error in the order of originalIDNodes")
if self.originalIDNodes.shape[0] != self.nodes.shape[0]:
print(f"{self.originalIDNodes.shape=}")
print(f"{self.nodes.shape=}")
raise RuntimeError("nodes and originalIDNodes are incompatible")
CheckType(self.nodes,MuscatFloat,"Nodes")
CheckType(self.originalIDNodes, MuscatIndex, "originalIDNodes")
nbNodes = self.nodes.shape[0]
# verification of min max in nodes tags
for nodeTagName, data in self.nodesTags.items():
ids = data.GetIds()
CheckType(ids, MuscatIndex, f"nodeTagName {nodeTagName}")
if len(ids) == 0:
continue
if ids[0] < 0:
raise Exception("Ids of '" + nodeTagName + "' tag out of bound (<0)")
if ids[-1] >= nbNodes:
print(f"{ids=}")
print(f"{nbNodes=}")
raise RuntimeError("Ids of '" + nodeTagName + "' tag out of bound > nbNodes")
# verification of elements
for elementType, data in self.elements.items():
if data.connectivity.flags['C_CONTIGUOUS'] is False:
raise Exception("Error : connectivity not C_CONTIGUOUS")
CheckType(data.connectivity, MuscatIndex, f"connectivity for {elementType}")
if len(data.connectivity.shape) != 2:
raise Exception(f"Wrong shape of connectivity of elements '{elementType}'")
if data.originalIds.shape[0] != data.connectivity.shape[0]:
print(f"{data.originalIds.shape[0]=}")
print(f"{data.connectivity.shape[0]=}")
raise Exception(f"connectivity and originalIds are incompatible '{elementType}'")
CheckType(data.originalIds, MuscatIndex, f"originalIds for {elementType}")
if data.originalIds.flags['C_CONTIGUOUS'] is False:
raise Exception("Error in the order of originalIds")
if ED.numberOfNodes[elementType] != data.connectivity.shape[1]:
print(f"{elementType=}")
print(f"{ED.numberOfNodes[elementType]=}")
print(f"{data.connectivity.shape[1]=}")
raise Exception("Incompatible number of columns of the connectivity array")
if len(data.connectivity) and np.amin(data.connectivity) < 0:
raise Exception(f"connectivity of '{elementType}' out of bound (<0)")
if len(data.connectivity) and np.amax(data.connectivity) >= nbNodes:
print(f"{data.connectivity=}")
print(f"{nbNodes=}")
raise Exception(f"connectivity of '{elementType}' out of bound > nbNodes")
nbElements = data.connectivity.shape[0]
# verification of min max in nodes tags
for tagName, tagData in data.tags.items():
ids = tagData.GetIds()
CheckType(ids, MuscatIndex, f"Element tag for {elementType} TagName {tagName}")
if len(ids) == 0:
continue
if ids[0] < 0:
print(f"{nbElements=}")
print(f"{ids=}")
raise Exception("Ids of '" + tagName + "' tag out of bound (<0)")
if ids[-1] >= nbElements:
print(f"{nbElements=}")
print(f"{ids=}")
raise Exception("Ids of '" + tagName + "' tag out of bound > nbElements")
[docs]def CheckIntegrity() -> str:
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles
from Muscat.Helpers.CheckTools import MustFailFunction
mesh1 = CreateMeshOfTriangles([[0, 0, 0], [1, 2, 3], [0, 1, 0]], [[0, 1, 2]])
mesh2 = CreateMeshOfTriangles([[0, 0, 0], [1, 2, 3], [0, 1, 0]], [[0, 1, 2]])
copy.copy(mesh1)
assert (mesh1 == mesh1)
assert (mesh1 != " ")
assert (mesh1 == mesh2)
mesh2.GetElementsOfType(ED.Triangle_6)
assert (mesh1 != mesh2)
mesh2.originalIDNodes[0] = 10
assert (mesh1 != mesh2)
mesh2.nodes[0, 0] = 0.1
assert (mesh1 != mesh2)
mesh2.props["extraProp"] = 1
assert (mesh1 != mesh2)
mesh1.nodeFields["1"] = np.array(1)
assert (mesh1 != mesh2)
mesh2.nodeFields["2"] = np.array(2)
assert (mesh1 != mesh2)
del mesh2.nodeFields["2"]
mesh2.nodeFields["1"] = np.array(2)
assert (mesh1 != mesh2)
del mesh1.nodeFields["1"]
del mesh2.nodeFields["1"]
mesh1.elemFields["1"] = np.array(1)
assert (mesh1 != mesh2)
mesh2.elemFields["2"] = np.array(2)
assert (mesh1 != mesh2)
del mesh2.elemFields["2"]
mesh2.elemFields["1"] = np.array(2)
assert (mesh1 != mesh2)
del mesh1.elemFields["1"]
del mesh2.elemFields["1"]
mesh1.nodesTags.CreateTag("test")
assert (mesh1 != mesh2)
mesh1.CopyProperties(mesh2)
np.testing.assert_equal(mesh1.GetNamesOfElementTags(), ['2D'])
mesh1.ConvertDataForNativeTreatment()
mesh1.SetNodes(mesh2.nodes, mesh2.originalIDNodes)
mesh1.SetNodes(mesh2.nodes, generateOriginalIDs=True)
np.testing.assert_equal(mesh1.GetElementsOriginalIDs(), [0])
mesh1.SetElementsOriginalIDs([1])
print(mesh1.elements.storage[0])
elements = mesh1.GetElementsOfType(ED.Triangle_3)
print(elements)
print("------------")
elements = mesh1.GetElementsOfType(ED.Bar_2)
elements.AddNewElement([1, 2], 1)
mesh1.nodeFields["ZeroPointField"] = np.zeros((mesh1.GetNumberOfNodes(), 1))
mesh1.elemFields["ZeroElementField"] = np.zeros((mesh1.GetNumberOfElements(), 1))
elements.GetNumberOfNodesPerElement()
print(elements)
print(mesh1)
mesh1.AddElementToTag(1, "SecondElement")
mesh1.AddElementToTag(1, "SecondElement(2)")
np.testing.assert_equal(mesh1.GetNumberOfElements(), 2)
boundingMin, boundingMax = mesh1.ComputeBoundingBox()
print(boundingMin)
print(boundingMax)
mesh1.nodesTags.CreateTag('toto')
print(mesh1.GetNodalTag('toto'))
print(mesh1.GetNodalTag('toto2'))
mesh1.AddElementToTagUsingOriginalId(1, "bars")
np.testing.assert_equal(mesh1.GetPosOfNodes()[1, 1], 2)
print(mesh1.PrepareForOutput())
print(mesh1.GetElementsInTag("bars"))
print(mesh1.GetElementsInTag("bars", useOriginalId=True))
print(mesh1.GetElementsOfType(ED.Bar_2).GetTag("toto"))
print(mesh1)
mesh1.DeleteElemTags(["SecondElement"])
mesh1.DeleteElemTags("SecondElement(2)")
print(mesh1)
resII = CreateMeshOfTriangles([[0, 0, 0], [1, 2, 3], [0, 1, 0]], [[0, 1, 2]])
resII.AddNodeToTagUsingOriginalId(0, "First Point")
resII.GenerateManufacturedOriginalIDs()
if resII.GetPointsDimensionality() != 3:
raise # pragma: no cover
if resII.GetElementsDimensionality() != 2:
raise # pragma: no cover
resII.Clean()
try:
resII.MergeElements(mesh1)
raise # pragma: no cover
except:
pass
resII.MergeElements(mesh1, force=True)
print("----")
print(resII)
print(resII.GetNumberOfElements(dim=2))
print(resII.elements.GetTagsNames())
del resII.elements[ED.Triangle_3]
resII.nodeFields["scalarNF"] = np.zeros(resII.GetNumberOfNodes(), dtype=MuscatFloat)
mesh1.nodeFields["scalarNF"] = np.arange(mesh1.GetNumberOfNodes(), dtype=MuscatFloat)
resII.nodeFields["vectorNF"] = np.zeros((resII.GetNumberOfNodes(), 2), dtype=MuscatFloat)
mesh1.nodeFields["vectorNF"] = np.ones((mesh1.GetNumberOfNodes(), 2), dtype=MuscatFloat)
resII.elemFields["scalarEF"] = np.zeros(resII.GetNumberOfElements(), dtype=MuscatFloat)
mesh1.elemFields["scalarEF"] = np.zeros(mesh1.GetNumberOfElements(), dtype=MuscatFloat)
resII.elemFields["vectorEF"] = np.zeros((resII.GetNumberOfElements(), 2), dtype=MuscatIndex)
mesh1.elemFields["vectorEF"] = np.ones((mesh1.GetNumberOfElements(), 2), dtype=MuscatIndex)
resII.elemFields["vectorEF_X"] = np.zeros((resII.GetNumberOfElements(), 2), dtype=MuscatFloat)
mesh1.elemFields["vectorEF_Y"] = np.ones((mesh1.GetNumberOfElements(), 2), dtype=MuscatFloat)
nSize1 = resII.GetNumberOfNodes()
eSize1 = resII.GetNumberOfElements()
nSize2 = mesh1.GetNumberOfNodes()
eSize2 = mesh1.GetNumberOfElements()
resII.Merge(mesh1)
for name, value in resII.elemFields.items():
np.testing.assert_equal(value.shape[0], eSize1 + eSize2)
for name, value in resII.nodeFields.items():
np.testing.assert_equal(value.shape[0], nSize1 + nSize2)
mesh1.GetElementsOfType(ED.Triangle_3).tags.CreateTag("invalid").SetIds([1000])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Triangle_3).tags["invalid"].SetIds([-1])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).connectivity = np.array([[100, 0]])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).connectivity = np.array([[-10, 0]])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).connectivity = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
mesh1.GetElementsOfType(ED.Bar_2).originalIds = np.arange(3)
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).originalIds = np.arange(10)[slice(0, 5, 2)]
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).originalIds = np.arange(2)
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).connectivity = np.array([0])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Bar_2).connectivity = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])[slice(0, 3, 2), :]
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.nodesTags.CreateTag("invalid").SetIds([100])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.nodesTags.CreateTag("invalid", False).SetIds([-1])
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.originalIDNodes = np.arange(100)
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.originalIDNodes = np.arange(100)[slice(0, 50, 2)]
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.originalIDNodes = np.zeros((2, 3))
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.nodes = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]])[slice(0, 3, 2), :]
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.nodes = np.zeros((2, 3, 4))
MustFailFunction(mesh1.VerifyIntegrity)
mesh1.GetElementsOfType(ED.Triangle_3).tags.CreateTag("invalid", False).SetIds([])
mesh1.GetElementsOfType(ED.Triangle_6).tags.CreateTag("empty", False)
mesh1.Clean()
# check View.
mesh1.View()
return "ok"
if __name__ == '__main__':
print(CheckIntegrity()) # pragma: no cover