# -*- 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.MeshContainers.Mesh import Mesh
from Muscat.Types import ArrayLike
from Muscat.MeshTools.MeshCreationTools import CreateMeshOfTriangles
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.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.bytes_ 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.MeshContainers.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 ComputeSignedDistance(mesh:Mesh,pos: np.ndarray)-> np.ndarray :
"""Function to compute the signed distance
Parameters
----------
mesh : Mesh
Mesh from which the signed distance field needs to be computed
pos : np.ndarray
positions where to compute the signed distance
Returns
-------
np.ndarray
signed distance field evaluated at pos
"""
from Muscat.ImplicitGeometry.ImplicitGeometryObjects import ImplicitGeometryExternalSurface
IGExternalSurface = ImplicitGeometryExternalSurface(mesh,method="normal") # legacy method is "bulk"
return -1*IGExternalSurface(pos) # inversing sign due to reverse convention from level set
[docs]
def CheckIntegrity_ComputeSignedDistance():
from Muscat.MeshTools.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
mesh1 = CreateConstantRectilinearMesh(dimensions=[10,10,10], spacing=[1, 1, 1], origin=[0, 0, 0])
mesh2 = CreateConstantRectilinearMesh(dimensions=[20,20,20], spacing=[1, 1, 1], origin=[-5, -5, -5])
ComputeSignedDistance(mesh1,mesh2.nodes)
[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.MeshTools.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():
CheckIntegrity_GetCellCenters()
CheckIntegrity_TagsToRefs()
CheckIntegrity_ComputeSignedDistance()
return "OK"
if __name__ == "__main__":
print(CheckIntegrity()) # pragma: no cover