# -*- 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.
#
import numpy as np
from scipy.sparse import coo_matrix
from Muscat.Types import MuscatIndex, MuscatFloat, ArrayLike
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.MeshContainers.ElementsContainers import StructuredElementsContainer, ElementsContainer
[docs]
def GetConstantRectilinearPoints(dimensions, origin, spacing) -> np.ndarray:
"""
Create a point array in the space
Returns
-------
numpy.array
A 2-dimensional array, the first axis corresponds to the node
index, the second axis corresponds to space dimension index.
"""
x = np.arange(dimensions[0])*spacing[0]+origin[0]
y = np.arange(dimensions[1])*spacing[1]+origin[1]
dim = len(dimensions)
nbNodes = np.prod(dimensions)
if dim == 2:
xv, yv = np.meshgrid(x, y, indexing='ij', copy=False)
nodes = np.empty((nbNodes, 2), dtype=MuscatFloat)
nodes[:, 0] = xv.ravel()
nodes[:, 1] = yv.ravel()
else:
z = np.arange(dimensions[2])*spacing[2]+origin[2]
xv, yv, zv = np.meshgrid(x, y, z, indexing='ij', copy=False)
nodes = np.empty((nbNodes, 3), dtype=MuscatFloat)
nodes[:, 0] = xv.ravel()
nodes[:, 1] = yv.ravel()
nodes[:, 2] = zv.ravel()
return nodes
[docs]
def GetMonoIndex(indices: ArrayLike, dimensions: ArrayLike) -> np.ndarray:
"""return the mono indices from multi indices
Parameters
----------
indices : ArrayLike
the ij or ijk indices
dimensions: ArrayLike
Number of element in each dimension
Returns
-------
np.ndarray
the mono index
"""
indices = np.asarray(indices,dtype=MuscatIndex)
dimensions = np.asarray(dimensions,dtype=MuscatIndex)
if len(indices.shape) == 1:
indexs = indices[np.newaxis]
else:
indexs = indices
if len(dimensions) == 3:
planeSize = dimensions[1] * dimensions[2]
res = planeSize*indexs[:, 0]
res += indexs[:, 1]*dimensions[2]
res += indexs[:, 2]
return res
else:
planeSize = dimensions[1]
return planeSize*indexs[:, 0]+indexs[:, 1]
[docs]
def GetMultiIndex(indices:ArrayLike, dimensions: ArrayLike) -> np.ndarray:
"""Return the multi-index (ijk) from linear indices
Parameters
----------
indices : ArrayLike
the indices of the nodes to treat
dimensions: ArrayLike
Number of element in each dimension
Returns
-------
np.ndarray
an array with the ijk indices for every point in indices
size (nb points, 2 (ij) or 3 (ijk) )
"""
indices = np.asarray(indices,dtype=MuscatIndex)
if len(dimensions) == 3:
planesize = dimensions[1] *dimensions[2]
res = np.empty((len(indices),3),dtype=MuscatIndex)
res[:,0] = indices // planesize
resyz = indices - res[:,0]*(planesize)
res[:,1] = resyz // dimensions[2]
res[:,2] = resyz - res[:,1]*dimensions[2]
return res
else:
res = np.empty((len(indices),2),dtype=MuscatIndex)
res[:,0] = indices // dimensions[1]
res[:,1] = indices - res[:,0]*(dimensions[1])
return res
[docs]
def CreateConstantRectilinearMesh(dimensions, origin, spacing) -> Mesh:
res = Mesh()
with res.WithModification(verifyIntegrity=False):
if len(dimensions) == 3:
ec = StructuredElementsContainer(ED.Hexahedron_8)
else:
ec = StructuredElementsContainer(ED.Quadrangle_4)
ec.SetDimensions(dimensions)
res.elements.AddContainer(ec)
res.nodes = GetConstantRectilinearPoints(dimensions=dimensions, origin=origin, spacing=spacing)
res.originalIDNodes = np.arange(res.GetNumberOfNodes(), dtype=MuscatIndex)
res.props["spacing"] = np.asarray(spacing,dtype=MuscatFloat)
res.props["origin"] = np.asarray(origin,dtype=MuscatFloat)
res.props["dimensions"] = np.asarray(dimensions, dtype=MuscatIndex)
res.props["IsConstantRectilinear"] = True
return res
[docs]
def ToUnstructuredMesh(mesh:Mesh, inPlace=True) -> Mesh:
if inPlace:
for et, data in mesh.elements.items():
if isinstance(data, StructuredElementsContainer):
ec = ElementsContainer(et)
ec.connectivity = data.connectivity.copy()
ec.originalIds = data.originalIds
ec.cpt = data.GetNumberOfElements()
ec.tags = data.tags
mesh.elements.AddContainer(ec)
return mesh
else:
res = Mesh()
res.nodes = mesh.nodes.copy()
res.originalIDNodes = np.arange(mesh.GetNumberOfNodes(),dtype=MuscatIndex)
res.props = {k:v for k,v in mesh.props.items() if k not in ["spacing", "origin", "dimensions", "IsConstantRectilinear"]}
for et, data in mesh.elements.items():
if isinstance(data, StructuredElementsContainer ):
ec = ElementsContainer(et)
ec.connectivity = data.connectivity.copy()
ec.originalIds = np.arange(ec.GetNumberOfElements(), dtype=MuscatIndex)
ec.cpt = data.GetNumberOfElements()
ec.tags = data.tags.Copy()
res.elements.AddContainer(ec)
else:
res.elements.AddContainer(data)
return res
[docs]
def CheckIntegrity(GUI: bool = False) -> str:
cmesh = CreateConstantRectilinearMesh(dimensions=[2,2,2],origin=[0,0,0], spacing=[1,1,1])
mesh = ToUnstructuredMesh(cmesh, inPlace = False)
for et, data in mesh.elements.items():
if isinstance(data, StructuredElementsContainer):
raise RuntimeError("The mesh still have a StructuredElementsContainer as elements")
mesh = ToUnstructuredMesh(cmesh, inPlace = True)
if mesh is not cmesh:
raise RuntimeError("Inplace error")
for et, data in mesh.elements.items():
if isinstance(data, StructuredElementsContainer):
raise RuntimeError("The mesh still have a StructuredElementsContainer as elements")
np.testing.assert_equal(GetMonoIndex([[0,0,0]],[2,3,4])[0], 0)
np.testing.assert_equal(GetMonoIndex([[1,2,3]],[2,3,4])[0], 23)
return "ok"
if __name__ == '__main__': # pragma: no cover
print(CheckIntegrity(True))
print("done")