# -*- 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 Callable, Union, Optional, Iterable, Tuple
import numpy as np
try: # pragma: no cover
from paraview.vtk import vtkPolyData, vtkUnstructuredGrid, vtkPoints, vtkIdList, vtkImageData, vtkCellArray, vtkMultiBlockDataSet, vtkRectilinearGrid
from paraview.vtk.util import numpy_support
except:
try:
# faster import only needed classes
from vtkmodules.vtkCommonCore import vtkPoints, vtkIdList
from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkUnstructuredGrid, vtkImageData, vtkCellArray, vtkMultiBlockDataSet, vtkRectilinearGrid
from vtkmodules.util import numpy_support
except: # pragma: no cover
from Muscat.Helpers.Logger import Info
Info("vtk not found, some functionalities will not be available")
from Muscat.Types import ArrayLike, MuscatFloat, MuscatIndex
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.Mesh import Mesh, ElementsContainer
from Muscat.Containers.MeshCreationTools import CreateMeshOfTriangles
from Muscat.TestData import GetTestDataPath
nbPointsTo2DCells = np.array([0, 1, 3, 5, 9], dtype=MuscatIndex)
# from file vtkCellType.h of the vtk sources
vtkNameByNumber = {}
vtkNameByNumber[0] = "VTK_EMPTY_CELL"
vtkNameByNumber[1] = "VTK_VERTEX"
vtkNameByNumber[2] = "VTK_POLY_VERTEX"
vtkNameByNumber[3] = "VTK_LINE"
vtkNameByNumber[4] = "VTK_POLY_LINE"
vtkNameByNumber[5] = "VTK_TRIANGLE"
vtkNameByNumber[6] = "VTK_TRIANGLE_STRIP"
vtkNameByNumber[7] = "VTK_POLYGON"
vtkNameByNumber[8] = "VTK_PIXEL"
vtkNameByNumber[9] = "VTK_QUAD"
vtkNameByNumber[10] = "VTK_TETRA"
vtkNameByNumber[11] = "VTK_VOXEL"
vtkNameByNumber[12] = "VTK_HEXAHEDRON"
vtkNameByNumber[13] = "VTK_WEDGE"
vtkNameByNumber[14] = "VTK_PYRAMID"
vtkNameByNumber[15] = "VTK_PENTAGONAL_PRISM"
vtkNameByNumber[16] = "VTK_HEXAGONAL_PRISM"
vtkNameByNumber[21] = "VTK_QUADRATIC_EDGE"
vtkNameByNumber[22] = "VTK_QUADRATIC_TRIANGLE"
vtkNameByNumber[23] = "VTK_QUADRATIC_QUAD"
vtkNameByNumber[36] = "VTK_QUADRATIC_POLYGON"
vtkNameByNumber[24] = "VTK_QUADRATIC_TETRA"
vtkNameByNumber[25] = "VTK_QUADRATIC_HEXAHEDRON"
vtkNameByNumber[26] = "VTK_QUADRATIC_WEDGE"
vtkNameByNumber[27] = "VTK_QUADRATIC_PYRAMID"
vtkNameByNumber[28] = "VTK_BIQUADRATIC_QUAD"
vtkNameByNumber[29] = "VTK_TRIQUADRATIC_HEXAHEDRON"
vtkNameByNumber[30] = "VTK_QUADRATIC_LINEAR_QUAD"
vtkNameByNumber[31] = "VTK_QUADRATIC_LINEAR_WEDGE"
vtkNameByNumber[32] = "VTK_BIQUADRATIC_QUADRATIC_WEDGE"
vtkNameByNumber[33] = "VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON"
vtkNameByNumber[34] = "VTK_BIQUADRATIC_TRIANGLE"
vtkNameByNumber[35] = "VTK_CUBIC_LINE"
vtkNameByNumber[41] = "VTK_CONVEX_POINT_SET"
vtkNameByNumber[42] = "VTK_POLYHEDRON"
vtkNameByNumber[51] = "VTK_PARAMETRIC_CURVE"
vtkNameByNumber[52] = "VTK_PARAMETRIC_SURFACE"
vtkNameByNumber[53] = "VTK_PARAMETRIC_TRI_SURFACE"
vtkNameByNumber[54] = "VTK_PARAMETRIC_QUAD_SURFACE"
vtkNameByNumber[55] = "VTK_PARAMETRIC_TETRA_REGION"
vtkNameByNumber[56] = "VTK_PARAMETRIC_HEX_REGION"
vtkNameByNumber[60] = "VTK_HIGHER_ORDER_EDGE"
vtkNameByNumber[61] = "VTK_HIGHER_ORDER_TRIANGLE"
vtkNameByNumber[62] = "VTK_HIGHER_ORDER_QUAD"
vtkNameByNumber[63] = "VTK_HIGHER_ORDER_POLYGON"
vtkNameByNumber[64] = "VTK_HIGHER_ORDER_TETRAHEDRON"
vtkNameByNumber[65] = "VTK_HIGHER_ORDER_WEDGE"
vtkNameByNumber[66] = "VTK_HIGHER_ORDER_PYRAMID"
vtkNameByNumber[67] = "VTK_HIGHER_ORDER_HEXAHEDRON"
vtkNameByNumber[68] = "VTK_LAGRANGE_CURVE"
vtkNameByNumber[69] = "VTK_LAGRANGE_TRIANGLE"
vtkNameByNumber[70] = "VTK_LAGRANGE_QUADRILATERAL"
vtkNameByNumber[71] = "VTK_LAGRANGE_TETRAHEDRON"
vtkNameByNumber[72] = "VTK_LAGRANGE_HEXAHEDRON"
vtkNameByNumber[73] = "VTK_LAGRANGE_WEDGE"
vtkNameByNumber[74] = "VTK_LAGRANGE_PYRAMID"
# ---------------------------------------------------------------------------
vtkNumberByElementName = {}
vtkNumberByElementName[ED.Point_1] = 1
vtkNumberByElementName[ED.Bar_2] = 3
vtkNumberByElementName[ED.Triangle_3] = 5
vtkNumberByElementName[ED.Quadrangle_4] = 9
vtkNumberByElementName[ED.Tetrahedron_4] = 10
# vtkNumberByElementName[ED.Hexahedron_8] = 11 # voxel
vtkNumberByElementName[ED.Hexahedron_8] = 12
vtkNumberByElementName[ED.Wedge_6] = 13
vtkNumberByElementName[ED.Pyramid_5] = 14
vtkNumberByElementName[ED.Bar_3] = 21
vtkNumberByElementName[ED.Triangle_6] = 22
vtkNumberByElementName[ED.Quadrangle_8] = 23
vtkNumberByElementName[ED.Tetrahedron_10] = 24
vtkNumberByElementName[ED.Hexahedron_20] = 25
vtkNumberByElementName[ED.Wedge_15] = 26
vtkNumberByElementName[ED.Pyramid_13] = 27
vtkNumberByElementName[ED.Quadrangle_9] = 28
vtkNumberByElementName[ED.Hexahedron_27] = 29
elementNameByVtkNumber = {}
elementPermutationByVtkNumber = {}
for key, vtkNumber in vtkNumberByElementName.items():
elementNameByVtkNumber[vtkNumber] = key
elementNameByVtkNumber[4] = ED.Bar_2
elementNameByVtkNumber[8] = ED.Quadrangle_4 # voxel
elementPermutationByVtkNumber[8] = [0, 1, 3, 2]
elementNameByVtkNumber[11] = ED.Hexahedron_8 # voxel
elementPermutationByVtkNumber[11] = [0, 1, 3, 2, 4, 5, 7, 6]
elementNameByVtkNumber[2] = ED.Point_1 # voxel
# if a field is of type [..] and the min max are 0 and 1 then the field is
# converted to a tag. the first type is used to encode tags in vtk
tagsTypes = [np.int8, np.uint8, int]
[docs]def VtkFieldToNumpyField(vtkField) -> Tuple[str, np.ndarray]:
from vtkmodules.util import numpy_support
from vtkmodules.vtkCommonCore import vtkStringArray
name = vtkField.GetName()
if isinstance(vtkField, vtkStringArray):
return (name, np.array([vtkField.GetValue(x) for x in range(vtkField.GetNumberOfValues())], dtype=np.str_))
data = numpy_support.vtk_to_numpy(vtkField)
return (name, data)
[docs]def NumpyFieldToVtkField(fieldData: np.ndarray, fieldName: str):
from vtkmodules.util import numpy_support
from vtkmodules.vtkCommonCore import vtkStringArray
if type(fieldData[0]) in [str, object, np.str_]: # working on list of string or numpy of objects
# for the moment only work for scalar (1 components fields )
VTK_data = vtkStringArray()
VTK_data.SetNumberOfComponents(1)
VTK_data.SetNumberOfTuples(len(fieldData))
for i, v in enumerate(fieldData):
VTK_data.InsertValue(i, str(v))
VTK_data.SetName(fieldName)
return VTK_data
else:
outputtype = numpy_support.get_vtk_array_type(fieldData.dtype)
VTK_data = numpy_support.numpy_to_vtk(num_array=fieldData, deep=True, array_type=outputtype)
if len(fieldData.shape) > 1:
VTK_data.SetNumberOfComponents(fieldData.shape[1])
VTK_data.SetName(fieldName)
return VTK_data
[docs]def ApplyVtkPipeline(mesh: Mesh, op: Callable) -> Mesh:
vtkMesh = MeshToVtk(mesh)
vtkOutputMesh = op(vtkMesh)
return VtkToMesh(vtkOutputMesh)
[docs]def PlotMesh(mesh: Union[Mesh, "vtkPolyData", "vtkUnstructuredGrid"], TagsAsFields: bool = False): # pragma: no cover
"""Plot the input mesh using vtk,
Parameters
----------
mesh : Union[Mesh, "vtkPolyData", "vtkUnstructuredGrid"]
mesh to be plotted
TagsAsFields : bool, optional
if the mesh is a Muscat mesh (Mesh) then this option can be used to control if the tags are converted to field before plotting , by default False
"""
from vtkmodules.vtkFiltersGeometry import vtkGeometryFilter
from vtkmodules.vtkRenderingCore import vtkPolyDataMapper, vtkActor, vtkRenderer, vtkRenderWindow, vtkRenderWindowInteractor, vtkColorTransferFunction
from vtkmodules.vtkFiltersCore import vtkAssignAttribute
from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
from vtkmodules.vtkInteractionWidgets import vtkButtonWidget, vtkTexturedButtonRepresentation2D, vtkOrientationMarkerWidget
from vtkmodules.vtkIOImage import vtkPNGReader
from vtkmodules.vtkRenderingAnnotation import vtkAxesActor
import vtkmodules.vtkRenderingOpenGL2
import vtkmodules.vtkRenderingFreeType
from Muscat.Containers.Mesh import Mesh
if isinstance(mesh, Mesh):
vtkMesh = MeshToVtk(mesh, TagsAsFields=TagsAsFields)
else:
vtkMesh = mesh
vGF = vtkGeometryFilter()
vGF.SetInputData(vtkMesh)
vGF.Update()
if vtkMesh.GetPointData():
nbArrays = vtkMesh.GetPointData().GetNumberOfArrays()
else:
nbArrays = 0
cylinderMapper = vtkPolyDataMapper()
if nbArrays > 0:
out2 = vtkAssignAttribute()
out2.SetInputConnection(vGF.GetOutputPort())
array = vtkMesh.GetPointData().GetArray(0)
out2.Assign(array.GetName(), "SCALARS", "POINT_DATA")
cylinderMapper.SetInputConnection(out2.GetOutputPort())
else:
cylinderMapper.SetInputConnection(vGF.GetOutputPort())
cylinderActor = vtkActor()
cylinderActor.SetMapper(cylinderMapper)
ren = vtkRenderer()
renWin = vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
style = vtkInteractorStyleTrackballCamera()
style.SetDefaultRenderer(ren)
iren.SetInteractorStyle(style)
ren.AddActor(cylinderActor)
ren.GradientBackgroundOn()
ren.SetBackground(0.3176, 0.3412, 0.4314)
ren.SetBackground2(0, 0, 0.1647)
renWin.SetSize(800, 600)
buttonWidget = vtkButtonWidget()
buttonWidget.SetInteractor(iren)
class Listener:
def __init__(self):
self.fildcpt = 0
self.minmaxs = dict()
def processStateChangeEvent(self, obj, ev):
self.fildcpt += 1
if vtkMesh.GetPointData().GetNumberOfArrays() == 0:
return
nb = self.fildcpt % vtkMesh.GetPointData().GetNumberOfArrays()
array = vtkMesh.GetPointData().GetArray(nb)
arrayName = array.GetName()
out2.Assign(arrayName, "SCALARS", "POINT_DATA")
if arrayName in self.minmaxs:
lo, hi = self.minmaxs[arrayName]
else:
lo, hi = array.GetRange()
self.minmaxs[arrayName] = (lo, hi)
lut = vtkColorTransferFunction()
lut.SetColorSpaceToHSV()
lut.SetColorSpaceToDiverging()
lut.AddRGBPoint(lo, 0.23137254902000001, 0.298039215686, 0.75294117647100001)
lut.AddRGBPoint((lo + hi) / 2, 0.86499999999999999, 0.86499999999999999, 0.86499999999999999)
lut.AddRGBPoint(hi, 0.70588235294099999, 0.015686274509800001, 0.149019607843)
lut.Build()
cylinderMapper.SetInterpolateScalarsBeforeMapping(True)
cylinderMapper.SetScalarRange(lo, hi)
cylinderMapper.SetUseLookupTableScalarRange(True)
cylinderMapper.SetLookupTable(lut)
cylinderMapper.SetScalarModeToUsePointData()
cylinderMapper.ScalarVisibilityOn()
cylinderMapper.SelectColorArray(array.GetName())
print("Plot of field {} min/max, {}/{}".format(arrayName, lo, hi))
renWin.Render()
pass
listener = Listener()
listener.processStateChangeEvent(None, None)
buttonWidget.AddObserver("StateChangedEvent", listener.processStateChangeEvent)
buttonRepresentation = vtkTexturedButtonRepresentation2D()
buttonRepresentation.SetNumberOfStates(1)
r = vtkPNGReader()
fileName = GetTestDataPath() + "Next.png"
r.SetFileName(fileName)
r.Update()
image = r.GetOutput()
buttonRepresentation.SetButtonTexture(0, image)
buttonRepresentation.SetPlaceFactor(1)
buttonRepresentation.PlaceWidget([0, 50, 0, 50, 0, 0])
buttonWidget.SetRepresentation(buttonRepresentation)
buttonWidget.On()
iren.Initialize()
axesActor = vtkAxesActor()
axes = vtkOrientationMarkerWidget()
axes.SetOrientationMarker(axesActor)
axes.SetInteractor(iren)
axes.EnabledOn()
axes.InteractiveOn()
ren.ResetCamera()
ren.GetActiveCamera().Zoom(1.5)
renWin.Render()
# Start the event loop.
iren.Start()
renWin.Finalize()
iren.TerminateApp()
# del renWin, iren
[docs]def MeshToVtk(mesh: Mesh, vtkobject=None, TagsAsFields=False) -> Union["vtkPolyData", "vtkUnstructuredGrid"]:
# From www.vtk;org/wp-content/updloads/2015/04/file-formats.pdf
if vtkobject is None:
usePoly = True
for elementsName, elementContainer in mesh.elements.items():
if ED.dimensionality[elementsName] == 3:
usePoly = False
break
if usePoly:
output = vtkPolyData()
else:
output = vtkUnstructuredGrid()
else:
output = vtkobject # pragma: no cover
if isinstance(output, vtkPolyData):
usePoly = True
else:
usePoly = False
if mesh.GetNumberOfNodes() == 0:
return output
if mesh.GetNumberOfElements() == 0:
return output
output.Allocate(mesh.GetNumberOfElements())
VTK_originalIDNodes = NumpyFieldToVtkField(mesh.originalIDNodes, "originalIds")
output.GetPointData().AddArray(VTK_originalIDNodes)
pts = vtkPoints()
if mesh.nodes.shape[1] == 3:
pts.SetData(numpy_support.numpy_to_vtk(num_array=mesh.nodes, deep=False))
else:
p = np.zeros((mesh.GetNumberOfNodes(), 3), dtype=MuscatFloat)
p[:, 0:2] = mesh.nodes
pts.SetData(numpy_support.numpy_to_vtk(num_array=p, deep=False))
output.SetPoints(pts)
VTK_originalIDsEl = NumpyFieldToVtkField(mesh.GetElementsOriginalIDs(), "originalIds")
output.GetCellData().AddArray(VTK_originalIDsEl)
if usePoly == False:
cellTypes = np.empty(mesh.GetNumberOfElements(), dtype=MuscatIndex)
offsets = np.empty(mesh.GetNumberOfElements() + 1, dtype=np.int64)
cpt = MuscatIndex(0)
offsetcpt = MuscatIndex(0)
for elementsName, elementContainer in mesh.elements.items():
nbElement = elementContainer.GetNumberOfElements()
cellTypes[cpt : cpt + nbElement] = vtkNumberByElementName[elementsName]
offsets[cpt : cpt + nbElement] = offsetcpt + np.arange(nbElement) * elementContainer.GetNumberOfNodesPerElement()
offsetcpt += elementContainer.GetNumberOfNodesPerElement() * nbElement
cpt += nbElement
offsets[cpt] = offsetcpt
connectivity = np.empty(offsetcpt, dtype=np.int64)
offsetcpt = MuscatIndex(0)
for elementsName, elementContainer in mesh.elements.items():
nbElement = elementContainer.GetNumberOfElements()
nbObjects = elementContainer.GetNumberOfNodesPerElement() * nbElement
connectivity[offsetcpt : offsetcpt + nbObjects] = elementContainer.connectivity.flatten()
offsetcpt += nbObjects
offsets = numpy_support.numpy_to_vtkIdTypeArray(offsets, deep=True)
connectivity = numpy_support.numpy_to_vtkIdTypeArray(connectivity, deep=True)
cellArray = vtkCellArray()
cellArray.SetData(offsets, connectivity)
output.SetCells(cellTypes.tolist(), cellArray)
else:
for elementsName, elementContainer in mesh.elements.items():
pointIds = vtkIdList()
npe = elementContainer.GetNumberOfNodesPerElement()
pointIds.SetNumberOfIds(npe)
vtkNumber = vtkNumberByElementName[elementsName]
for e in range(elementContainer.GetNumberOfElements()):
for i in range(npe):
pointIds.SetId(i, elementContainer.connectivity[e, i])
output.InsertNextCell(vtkNumber, pointIds)
if hasattr(mesh, "nodeFields"):
for name, data in mesh.nodeFields.items():
if data is None: # pragma: no cover
continue
if np.size(data) // mesh.GetNumberOfNodes() != np.size(data) / mesh.GetNumberOfNodes(): # pragma: no cover
print("field (" + str(name) + ") is not consistent : it has " + str(np.size(data)) + " values and the mesh has " + str(mesh.GetNumberOfNodes()) + " nodes")
raise Exception("field (" + str(name) + ") is not consistent : it has " + str(np.size(data)) + " values and the mesh has " + str(mesh.GetNumberOfNodes()) + " nodes")
VTK_data = NumpyFieldToVtkField(data, name)
output.GetPointData().AddArray(VTK_data)
continue
if TagsAsFields:
tagMask: np.ndarray = np.empty(mesh.GetNumberOfNodes(), dtype=tagsTypes[0])
for tag in mesh.nodesTags:
tag.GetIdsAsMask(output=tagMask)
VTK_data = NumpyFieldToVtkField(tagMask, tag.name)
output.GetPointData().AddArray(VTK_data)
continue
nbElements = mesh.GetNumberOfElements()
if nbElements > 0:
if hasattr(mesh, "elemFields"):
for name, data in mesh.elemFields.items():
if data is None: # pragma: no cover
continue
if np.size(data) / nbElements != np.size(data) // nbElements: # pragma: no cover
print("field (" + str(name) + ") is not consistent : it has " + str(np.size(data)) + " values and the mesh has " + str(nbElements) + " elements")
continue
VTK_data = NumpyFieldToVtkField(data, name)
output.GetCellData().AddArray(VTK_data)
continue
if TagsAsFields:
elementTags = mesh.GetNamesOfElementTags()
for tagname in elementTags:
ids = mesh.GetElementsInTag(tagname)
tagMask = np.zeros(nbElements, dtype=tagsTypes[0])
tagMask[ids] = True
VTK_data = NumpyFieldToVtkField(tagMask, tagname)
output.GetCellData().AddArray(VTK_data)
continue
return output
[docs]def AddCellsFromVtkCellArrayToMesh(cells, cellTypesArray: Union[ArrayLike, None], out: Mesh, originalIdsOffset=0):
"""Function to migrate cells from a vtkCellArray to a Muscattured Mesh
Parameters
----------
cells : vtkCellArray
the cell to be created in the mesh
cellTypesArray : ArrayLike
the vtk cell type to use for in the vtkCellArray. In the case is none we infer
the type of cell from the number of nodes (for the polys in a vtkPolyData)
out : Mesh
the output mesh
originalIdsOffset: offset to apply to the original ids
"""
offsets = numpy_support.vtk_to_numpy(cells.GetOffsetsArray())
if cellTypesArray is None:
cellTypesArray = nbPointsTo2DCells[offsets[1:] - offsets[:-1]]
types, typeNB = np.unique(cellTypesArray, return_counts=True)
elemtype_index = dict()
for t, tnb in zip(types, typeNB):
et = elementNameByVtkNumber[t]
elements = out.GetElementsOfType(et)
elements.Reserve(elements.GetNumberOfElements() + tnb)
elemtype_index[t] = (elements, np.where(cellTypesArray == t)[0])
connectivity = numpy_support.vtk_to_numpy(cells.GetConnectivityArray())
nbNewElements = 0
for vtktype, (elements, index) in elemtype_index.items():
cpt0 = elements.cpt
numberOfNodesPerElement = elements.GetNumberOfNodesPerElement()
nbNewElements = len(index)
cpt1 = cpt0 + nbNewElements
mask = np.arange(numberOfNodesPerElement, dtype=MuscatIndex) * np.ones((len(index), numberOfNodesPerElement), dtype=MuscatIndex) + offsets[index, None]
localConnectivity = connectivity[mask]
localConnectivity.shape = (nbNewElements, numberOfNodesPerElement)
elements.connectivity[cpt0:cpt1, :] = localConnectivity
elements.originalIds[cpt0:cpt1] = index + originalIdsOffset
if vtktype in elementPermutationByVtkNumber:
elements.connectivity[cpt0:cpt1, :] = elements.connectivity[cpt0:cpt1, elementPermutationByVtkNumber[vtktype]]
elements.cpt += nbNewElements
nbNewElements += nbNewElements
return nbNewElements
[docs]def VtkToMesh(vtkmesh: Union["vtkPolyData", "vtkUnstructuredGrid"], meshobject=None, FieldsAsTags=True):
if meshobject is None:
out = Mesh()
else:
out = meshobject
from vtkmodules.util import numpy_support
if vtkmesh.IsA("vtkImageData"):
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
out = CreateConstantRectilinearMesh(spacing=vtkmesh.GetSpacing(), origin=vtkmesh.GetOrigin(), dimensions=vtkmesh.GetDimensions())
originalIDNodesReady = True
elif vtkmesh.IsA("vtkPolyData"):
data = vtkmesh.GetPoints().GetData()
out.nodes = np.ascontiguousarray(numpy_support.vtk_to_numpy(data), dtype=MuscatFloat)
originalIDNodesReady = False
nc = vtkmesh.GetNumberOfCells()
cpt = MuscatIndex(0)
verts = vtkmesh.GetVerts()
if verts.GetNumberOfCells() != 0:
cpt += AddCellsFromVtkCellArrayToMesh(verts, np.full(verts.GetNumberOfCells(), 1), out, originalIdsOffset=cpt)
lines = vtkmesh.GetLines()
if lines.GetNumberOfCells() != 0:
cpt += AddCellsFromVtkCellArrayToMesh(lines, np.full(lines.GetNumberOfCells(), 3), out, originalIdsOffset=cpt)
polys = vtkmesh.GetPolys()
if polys.GetNumberOfCells() != 0:
cpt += AddCellsFromVtkCellArrayToMesh(polys, None, out, originalIdsOffset=cpt)
strips = vtkmesh.GetStrips()
if strips.GetNumberOfCells() != 0:
offsets = numpy_support.vtk_to_numpy(strips.GetOffsetsArray())
nbtrianglesPerCell = offsets[1:] - offsets[:-1] - 2
print(nbtrianglesPerCell)
newTris = np.sum(nbtrianglesPerCell)
elements = out.GetElementsOfType(ED.Triangle_3)
cpt0 = elements.cpt
elements.Reserve(elements.GetNumberOfElements() + newTris)
connectivity = numpy_support.vtk_to_numpy(strips.GetConnectivityArray())
for connCpt, tris in enumerate(nbtrianglesPerCell):
for tri in range(tris):
elements.connectivity[elements.cpt, :] = connectivity[offsets[connCpt] + tri : offsets[connCpt] + tri + 3]
elements.originalIds[elements.cpt] = cpt0 + cpt
elements.cpt += 1
cpt += 1
cpt0 += 1
elif vtkmesh.IsA("vtkStructuredGrid"):
data = vtkmesh.GetPoints().GetData()
out.nodes = np.ascontiguousarray(numpy_support.vtk_to_numpy(data), dtype=MuscatFloat)
originalIDNodesReady = False
dimensions = [0, 1, 2]
vtkmesh.GetCellDims(dimensions)
from Muscat.Containers.ElementsContainers import StructuredElementsContainer
if dimensions[2] > 1:
ec = StructuredElementsContainer(ED.Hexahedron_8)
ec.SetDimensions(np.asarray(dimensions) + 1)
else:
ec = StructuredElementsContainer(ED.Quadrangle_4)
ec.SetDimensions(np.asarray(dimensions[0:2]) + 1)
out.elements.AddContainer(ec)
elif vtkmesh.IsA("vtkRectilinearGrid"):
points = vtkPoints()
vtkmesh.GetPoints(points)
data = points.GetData()
out.nodes = np.ascontiguousarray(numpy_support.vtk_to_numpy(data), dtype=MuscatFloat)
originalIDNodesReady = False
dimensions = [0, 1, 2]
dimensions = vtkmesh.GetDimensions()
from Muscat.Containers.ElementsContainers import StructuredElementsContainer
if dimensions[2] > 1:
ec = StructuredElementsContainer(ED.Hexahedron_8)
ec.SetDimensions(dimensions)
else:
ec = StructuredElementsContainer(ED.Quadrangle_4)
ec.SetDimensions(np.asarray(dimensions[0:2]))
out.elements.AddContainer(ec)
else:
data = vtkmesh.GetPoints().GetData()
out.nodes = np.ascontiguousarray(numpy_support.vtk_to_numpy(data), dtype=MuscatFloat)
originalIDNodesReady = False
AddCellsFromVtkCellArrayToMesh(vtkmesh.GetCells(), vtkmesh.GetCellTypesArray(), out)
if vtkmesh.GetPointData().GetNumberOfArrays():
for f in range(vtkmesh.GetPointData().GetNumberOfArrays()):
data = vtkmesh.GetPointData().GetArray(f)
(name, field) = VtkFieldToNumpyField(data)
if name == "originalIds":
out.originalIDNodes = np.asarray(field, dtype=MuscatIndex)
originalIDNodesReady = True
else:
if FieldsAsTags and field.size == out.GetNumberOfNodes() and field.dtype in tagsTypes and np.min(field) >= 0 and np.max(field) <= 1:
out.nodesTags.CreateTag(name).SetIds(np.where(field)[0])
else:
out.nodeFields[name] = field
if not originalIDNodesReady:
out.originalIDNodes = np.arange(out.GetNumberOfNodes(), dtype=MuscatIndex)
out.PrepareForOutput()
EOIds = out.GetElementsOriginalIDs()
EOIds = np.argsort(EOIds)
if vtkmesh.GetCellData().GetNumberOfArrays():
for f in range(vtkmesh.GetCellData().GetNumberOfArrays()):
data = vtkmesh.GetCellData().GetAbstractArray(f)
(name, field) = VtkFieldToNumpyField(data)
Elfield = np.empty(field.shape, dtype=field.dtype)
if len(field.shape) > 1:
Elfield[EOIds, :] = field
else:
Elfield[EOIds] = field
if name == "originalIds":
out.SetElementsOriginalIDs(Elfield)
else:
if FieldsAsTags and field.size == out.GetNumberOfElements() and field.dtype in tagsTypes and np.min(field) >= 0 and np.max(field) <= 1:
cpt = MuscatIndex(0)
for elname, data in out.elements.items():
nn = data.GetNumberOfElements()
data.tags.CreateTag(name).SetIds(np.where(Elfield[cpt : cpt + nn])[0])
cpt += nn
else:
out.elemFields[name] = Elfield
continue
return out
[docs]def GetOutputVtk(request, inInfoVec, outInfoVec, copyAttr: bool = True, outputNumber=0) -> Union["vtkPolyData", "vtkUnstructuredGrid"]: # pragma: no cover
from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid
output = vtkUnstructuredGrid.GetData(outInfoVec, outputNumber)
if copyAttr:
input0 = GetInputVtk(request, inInfoVec, outInfoVec)
output.CopyAttributes(input0)
return output
[docs]def SetOutputMuscat(request, inInfoVec, outInfoVec, outMesh: Mesh, TagsAsFields: bool = False) -> None: # pragma: no cover
output = GetOutputVtk(request, inInfoVec, outInfoVec, False)
MeshToVtk(outMesh, output, TagsAsFields=TagsAsFields)
[docs]def VtkToPartitionedMesh(vtkObject, OP: Callable[[Union["vtkPolyData", "vtkUnstructuredGrid"]], Mesh] = VtkToMesh):
"""Experimental, Vtk To PartitionedMesh
Parameters
----------
vtkObject : vtkMultiBlockDataSet
_description_
OP : Callable[[Union[vtkPolyData, vtkUnstructuredGrid]], Mesh], default VtkToMesh() function
call operator to convert vtk mesh to mesh
Returns
-------
PartitionedMesh
A partitioned mesh with every block as a Mesh
"""
from Muscat.Experimental.Containers.PartitionedMesh import PartitionedMesh
if vtkObject.IsA("vtkMultiBlockDataSet"):
res = PartitionedMesh()
nb = vtkObject.GetNumberOfBlocks()
for i in range(nb):
block = vtkObject.GetBlock(i)
res.storage.append(OP(block))
return res
else:
res = PartitionedMesh()
res.storage.append(OP(vtkObject))
return res
[docs]def CellDataToPoint(mesh: Mesh, cellfields: np.ndarray) -> np.ndarray:
"""Applies the CellDataToPointData from vtk.
Supported only for the dimensionality of the mesh (no mix of elements of different
dimensions)
Parameters
----------
mesh : Mesh
Mesh containing the cells and vertices concerned by the conversion.
cellfield : np.ndarray
of size (number of fields, number of elements). Cell fields to convert to Point field.
Returns
-------
np.ndarray
of size (number of points, number of fields). Field converted at the vertices of the mesh.
"""
from Muscat.Bridges import vtkBridge as vB
from vtkmodules.util import numpy_support
from vtkmodules.vtkFiltersCore import vtkCellDataToPointData
vtkMesh = vB.MeshToVtk(mesh)
nbFields = cellfields.shape[0]
for i in range(nbFields):
vtkMesh.GetCellData().AddArray(vB.NumpyFieldToVtkField(cellfields[i, :], "field_" + str(i)))
cellToPoint = vtkCellDataToPointData()
cellToPoint.SetInputData(vtkMesh)
cellToPoint.Update()
nbArrays = vtkMesh.GetCellData().GetNumberOfArrays()
cellData = cellToPoint.GetOutput().GetPointData()
res = [numpy_support.vtk_to_numpy(cellData.GetArray(i)) for i in range(nbArrays - nbFields, nbArrays)]
return np.array(res)
[docs]def ReadMeshAndPopulateVtkObject(filename, vtkobject=None, TagsAsFields=False): # pragma: no cover
from Muscat.IO.UniversalReader import ReadMesh as UReadMesh
mesh = UReadMesh(filename)
from Muscat.Bridges.vtkBridge import MeshToVtk
return MeshToVtk(mesh, vtkobject, TagsAsFields=TagsAsFields)
[docs]def CheckIntegrity_VtkToMesh(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]
tet = [[0, 1, 2, 3]]
res = CreateMeshOf(points, tet, elemName=ED.Tetrahedron_4)
res.nodeFields = {"x": res.nodes[:, 0].flatten(), "Pos": res.nodes}
res.nodesTags.CreateTag("FirstPoint").AddToTag(0)
res.elemFields = {"SecondPoint": res.GetElementsOfType(ED.Tetrahedron_4).connectivity[:, 1].flatten().astype(float), "conn": res.GetElementsOfType(ED.Tetrahedron_4).connectivity}
res.GetElementsOfType(ED.Tetrahedron_4).tags.CreateTag("FirstTetrahedron").AddToTag(0)
sol = MeshToVtk(res, TagsAsFields=True)
print("CheckIntegrity_VtkToMesh :")
print(res)
print(VtkToMeshOnlyMeta(sol, FieldsAsTags=True).elemTags)
resII = VtkToMesh(sol, FieldsAsTags=True)
print(resII)
from Muscat.Containers.MeshTools import IsClose
print(res)
print(resII)
if not IsClose(res, resII): # pragma: no cover
raise (Exception("The meshes are not equal"))
return "ok"
[docs]def CheckIntegrity_VtkToMesh2D(GUI: bool = False) -> str:
res = CreateMeshOfTriangles([[0, 0], [1, 0], [1, 1], [1, 1]], [[0, 1, 2], [0, 2, 3]])
res.nodeFields = {"x": res.nodes[:, 0].flatten(), "Pos": res.nodes}
res.nodesTags.CreateTag("FirstPoint").AddToTag(0)
res.elemFields = {"SecondPoint": res.GetElementsOfType(ED.Triangle_3).connectivity[:, 1].flatten().astype(float), "conn": res.GetElementsOfType(ED.Triangle_3).connectivity}
res.GetElementsOfType(ED.Triangle_3).tags.CreateTag("FirstTriangle").AddToTag(0)
sol = MeshToVtk(res, TagsAsFields=True)
print("CheckIntegrity_VtkToMesh :")
print(res)
print(VtkToMeshOnlyMeta(sol))
resII = VtkToMesh(sol, FieldsAsTags=True)
resII.nodes = resII.nodes[:, 0:2]
print(resII)
from Muscat.Containers.MeshTools import IsClose
print(res)
print(resII)
if not IsClose(res, resII): # pragma: no cover
raise (Exception("The meshes are not equal"))
return "ok"
[docs]def CheckIntegrity_MeshToVtk(GUI: bool = False) -> str:
res = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
res.nodeFields = {"x": res.nodes[:, 0].flatten(), "Pos": res.nodes}
res.nodesTags.CreateTag("FirstPoint").AddToTag(0)
res.elemFields = {
"SecondPoint": res.GetElementsOfType(ED.Triangle_3).connectivity[:, 1].flatten(),
"conn": res.GetElementsOfType(ED.Triangle_3).connectivity,
"FE Names": np.array(["c2d3", "c2d3"], dtype=str),
}
res.GetElementsOfType(ED.Triangle_3).tags.CreateTag("FirstTriangle").AddToTag(0)
sol = MeshToVtk(res, TagsAsFields=True)
print(sol)
resII = VtkToMesh(sol, meshobject=Mesh())
from Muscat.Containers.MeshTools import IsClose
if not IsClose(res, resII): # pragma: no cover
raise (Exception("The meshes are not equal"))
# test a 2D mesh
res = CreateMeshOfTriangles([[0, 0], [1, 0], [0, 1], [1, 1]], [[0, 1, 2], [2, 1, 3]])
if GUI: # pragma: no cover
res.nodeFields["Field1"] = np.array([30, 20, 30, 1])
res.nodeFields["Field2"] = np.array([0, 1, 0, 1]) + 0.1
PlotMesh(res)
sol = MeshToVtk(res)
PlotMesh(sol)
return "OK"
[docs]def CheckIntegrity_RectilinearMesh(GUI: bool = False) -> str:
print("*** CheckIntegrity_RectilinearMesh 3D **** ")
rgrid = vtkRectilinearGrid()
xCoords = np.linspace(0, 1, num=3)
yCoords = np.linspace(0.1, 1.2, num=4)
zCoords = np.linspace(0.2, 1.4, num=5)
rgrid.SetDimensions(len(xCoords), len(yCoords), len(zCoords))
rgrid.SetXCoordinates(NumpyFieldToVtkField(xCoords, "xCoords"))
rgrid.SetYCoordinates(NumpyFieldToVtkField(yCoords, "yCoords"))
rgrid.SetZCoordinates(NumpyFieldToVtkField(zCoords, "zCoords"))
crmII = VtkToMesh(rgrid, FieldsAsTags=True)
np.testing.assert_equal(crmII.GetNumberOfNodes(), rgrid.GetNumberOfPoints())
np.testing.assert_equal(crmII.GetNumberOfElements(), rgrid.GetNumberOfCells())
print("*** CheckIntegrity_RectilinearMesh 2D **** ")
rgrid = vtkRectilinearGrid()
rgrid.SetDimensions(len(xCoords), len(yCoords), 1)
rgrid.SetXCoordinates(NumpyFieldToVtkField(xCoords, "xCoords"))
rgrid.SetYCoordinates(NumpyFieldToVtkField(yCoords, "yCoords"))
zCoords = np.linspace(0, 0, num=1)
rgrid.SetZCoordinates(NumpyFieldToVtkField(zCoords, "zCoords"))
crmII = VtkToMesh(rgrid, FieldsAsTags=True)
np.testing.assert_equal(crmII.GetNumberOfNodes(), rgrid.GetNumberOfPoints())
np.testing.assert_equal(crmII.GetNumberOfElements(), rgrid.GetNumberOfCells())
return "ok"
[docs]def CheckIntegrity_ConstantRectilinearMesh(GUI: bool = False) -> str:
print("*** CheckIntegrity_ConstantRectilinearMesh 3D **** ")
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
crm = CreateConstantRectilinearMesh(dimensions=[2, 3, 4], spacing=[0.1, 0.2, 0.3], origin=[-1, -1, -1])
crm.nodesTags.CreateTag("FirstPoint").SetIds([0])
crm.nodeFields["nodeData"] = np.arange(crm.GetNumberOfNodes()) * 0.3
crm.nodeFields["nodeData3coom"] = np.zeros((crm.GetNumberOfNodes(), 3)) * 0.3
vtkcrm = MeshToVtk(crm, TagsAsFields=True)
crmII = VtkToMesh(vtkcrm, FieldsAsTags=True)
from Muscat.Containers.MeshTools import IsClose
print(crm)
print(vtkcrm)
print(crmII)
print(crm.nodeFields["nodeData"])
crmII.nodeFields["nodeData"] = crmII.nodeFields["nodeData"].flatten()
crmII.nodeFields["nodeData3coom"] = crmII.nodeFields["nodeData3coom"].reshape(-1, 3)
if not IsClose(crm, crmII): # pragma: no cover
raise (Exception("The meshes are not equal"))
print("*** CheckIntegrity_ConstantRectilinearMesh 2D**** ")
from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
crm = CreateConstantRectilinearMesh(dimensions=[2, 3], spacing=[0.1, 0.2], origin=[-1, -1])
nodes = np.zeros((crm.GetNumberOfNodes(), 3), dtype=MuscatFloat)
nodes[:, :2] = crm.nodes
crm.nodes = nodes
crm.nodesTags.CreateTag("FirstPoint").SetIds([0])
crm.nodeFields["nodeData"] = np.arange(crm.GetNumberOfNodes()) * 0.3
crm.nodeFields["nodeData3coom"] = np.zeros((crm.GetNumberOfNodes(), 3)) * 0.3
vtkcrm = MeshToVtk(crm, TagsAsFields=True)
crmII = VtkToMesh(vtkcrm, FieldsAsTags=True)
from Muscat.Containers.MeshTools import IsClose
print(crm)
print(vtkcrm)
print(crmII)
print(crm.nodeFields["nodeData"])
crmII.nodeFields["nodeData"] = crmII.nodeFields["nodeData"].flatten()
crmII.nodeFields["nodeData3coom"] = crmII.nodeFields["nodeData3coom"].reshape(-1, 3)
if not IsClose(crm, crmII): # pragma: no cover
raise (Exception("The meshes are not equal"))
return "ok"
[docs]def checkIntegrity_ApplyVtkPipeline(GUI) -> str:
res = CreateMeshOfTriangles([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], [[0, 1, 2], [0, 2, 3]])
class op:
def __call__(self, vtkinput):
from vtkmodules.vtkFiltersModeling import vtkLinearExtrusionFilter
from vtkmodules.vtkFiltersCore import vtkTriangleFilter
extrude = vtkLinearExtrusionFilter()
extrude.SetInputData(vtkinput)
extrude.SetExtrusionTypeToNormalExtrusion()
extrude.SetVector(2.0, 2.0, 2.0)
extrude.Update()
cleaner = vtkTriangleFilter()
cleaner.SetInputData(extrude.GetOutput())
cleaner.Update()
return cleaner.GetOutput()
res = ApplyVtkPipeline(res, op())
if GUI: # pragma: no cover
PlotMesh(res)
return "ok"
[docs]def CheckIntegrity_CellDataToPoint(GUI: bool = False) -> str:
from Muscat.Containers.MeshCreationTools import CreateMeshOf
points = [[-0.5, -0.5, -0.5], [2.5, -0.5, -0.5], [-0.5, 2.5, -0.5], [-0.5, -0.5, 2.5], [2.5, 2.5, 2.5]]
tets = [[0, 1, 2, 3]]
mesh = CreateMeshOf(points, tets, ED.Tetrahedron_4)
cellField = np.array([[1.0], [2.0]])
CellDataToPoint(mesh, cellField)
return "ok"
[docs]def CheckIntegrity_Empty(GUI: bool = False) -> str:
emptyMesh = Mesh()
vtkMesh = MeshToVtk(emptyMesh)
vtkMesh = MeshToVtk(emptyMesh, vtkPolyData())
emptyMesh.SetNodes([[0, 0, 0]])
vtkMesh = MeshToVtk(emptyMesh, vtkUnstructuredGrid())
VtkToMeshOnlyMeta(vtkMesh)
VtkToMeshOnlyMeta(None)
return "ok"
[docs]def CheckIntegrity_Pixel(GUI: bool = False) -> str:
vtkMesh = vtkUnstructuredGrid()
vtkMesh.Allocate(1)
pts = vtkPoints()
pts.SetData(numpy_support.numpy_to_vtk(num_array=np.asarray([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]), deep=False))
vtkMesh.SetPoints(pts)
pointIds = vtkIdList()
pointIds.SetNumberOfIds(4)
for i in range(4):
pointIds.SetId(i, [0, 1, 3, 2][i])
vtkMesh.InsertNextCell(8, pointIds)
mesh = VtkToMesh(vtkMesh)
assert mesh.GetNumberOfElements() == 1
assert mesh.GetNumberOfNodes() == 4
vtkMesh = vtkImageData()
vtkMesh.SetDimensions([2, 2, 2])
vtkMesh.SetOrigin([0, 0, 0])
vtkMesh.SetSpacing([1, 1, 1])
mesh = VtkToMesh(vtkMesh)
return "ok"
[docs]def CheckIntegrity_Poly(GUI: bool = False) -> str:
vtkMesh = vtkPolyData()
vtkMesh.Allocate(3)
pts = vtkPoints()
pts.SetData(numpy_support.numpy_to_vtk(num_array=np.asarray([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]]), deep=False))
vtkMesh.SetPoints(pts)
pointIds = vtkIdList()
for vtkElemNb, indices in zip([1, 3, 5, 6], [[0], [0, 1], [0, 1, 2], [0, 1, 3, 2]]):
pointIds.SetNumberOfIds(len(indices))
for i in range(len(indices)):
pointIds.SetId(i, indices[i])
vtkMesh.InsertNextCell(vtkElemNb, pointIds)
myMesh = VtkToMesh(vtkMesh)
print(vtkMesh)
print(myMesh)
assert myMesh.GetNumberOfElements() == 5
return "ok"
[docs]def CheckIntegrity_VtkToPartitionedMesh(GUI: bool = False) -> str:
vtkMB = vtkMultiBlockDataSet()
vtkMB.SetNumberOfBlocks(2)
res1 = CreateMeshOfTriangles([[0, 0, 2], [1, 0, 2], [0, 1, 2], [0, 0, 3]], [[0, 1, 2], [0, 2, 3]])
res2 = CreateMeshOfTriangles([[0, 0, 2], [1, 0, 2], [0, 1, 2], [0, 0, 3]], [[0, 1, 2], [0, 2, 3]])
vtkMB.SetBlock(0, MeshToVtk(res1))
vtkMB.SetBlock(1, MeshToVtk(res1))
VtkToPartitionedMesh(vtkMB)
VtkToPartitionedMesh(MeshToVtk(res1))
return "ok"
[docs]def CheckIntegrity(GUI: bool = False) -> str:
try:
numpy_support
except: # pragma: no cover
return "skip"
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities
toTest: list[Callable[[bool], str]] = [
CheckIntegrity_MeshToVtk,
CheckIntegrity_VtkToMesh2D,
CheckIntegrity_VtkToMesh,
CheckIntegrity_ConstantRectilinearMesh,
CheckIntegrity_RectilinearMesh,
checkIntegrity_ApplyVtkPipeline,
CheckIntegrity_CellDataToPoint,
CheckIntegrity_Empty,
CheckIntegrity_Pixel,
CheckIntegrity_Poly,
CheckIntegrity_VtkToPartitionedMesh,
]
return RunListOfCheckIntegrities(toTest)
if __name__ == "__main__":
print(CheckIntegrity(True)) # pragma: no cover