Merge branch '546-support-for-more-vtk-elements' into 'development'

use general Lagrange elements instead of P1 elements

Closes #546

See merge request damask/DAMASK!1128

Author: Javier Velo

python/damask/_geomgrid.py

@@ -1382,4 +1382,4 @@ def get_grain_boundaries(self,
             connectivity.append(np.block([base_nodes + o[i][k] for k in range(4)]))
 
         coords = grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
-        return VTK.from_unstructured_grid(coords,np.vstack(connectivity),'QUAD')
+        return VTK.from_unstructured_grid(coords,np.vstack(connectivity),'QUADRILATERAL')

python/damask/_vtk.py

@@ -30,10 +30,10 @@
     vtkRectilinearGrid,
     vtkUnstructuredGrid,
     vtkPolyData,
-    VTK_TRIANGLE,
-    VTK_QUAD,
-    VTK_TETRA,
-    VTK_HEXAHEDRON,
+    VTK_LAGRANGE_TRIANGLE,
+    VTK_LAGRANGE_QUADRILATERAL,
+    VTK_LAGRANGE_TETRAHEDRON,
+    VTK_LAGRANGE_HEXAHEDRON,
 )
 from vtkmodules.vtkIOLegacy import (
     vtkGenericDataObjectReader,
@@ -242,11 +242,14 @@ def from_image_data(cells: IntSequence,
     @staticmethod
     def from_unstructured_grid(nodes: np.ndarray,
                                connectivity: np.ndarray,
-                               cell_type: Literal['TRIANGLE', 'TETRA', 'QUAD', 'HEXAHEDRON']) -> 'VTK':
+                               cell_type: Literal['TRIANGLE', 'TETRAHEDRON', 'QUADRILATERAL', 'HEXAHEDRON']
+                              ) -> 'VTK':
         """
         Create VTK of type vtkUnstructuredGrid.
 
         This is the common type for mesh solver results.
+        Elements are of the Lagrange type and all elements
+        should have the same type and order.
 
         Parameters
         ----------
@@ -255,7 +258,7 @@ def from_unstructured_grid(nodes: np.ndarray,
         connectivity : numpy.ndarray of np.dtype = np.int64
             Cell connectivity (0-based), first dimension determines #Cells,
             second dimension determines #Nodes/Cell.
-        cell_type : {'TRIANGLE', 'QUAD', 'TETRA', 'HEXAHEDRON'}
+        cell_type : {'TRIANGLE', 'QUADRILATERAL', 'TETRAHEDRON', 'HEXAHEDRON'}
             Name of the vtkCell subclass.
 
         Returns
@@ -273,9 +276,10 @@ def from_unstructured_grid(nodes: np.ndarray,
 
         vtk_data = vtkUnstructuredGrid()
         vtk_data.SetPoints(vtk_nodes)
-        cell_types = {'TRIANGLE':VTK_TRIANGLE, 'QUAD':VTK_QUAD,
-                      'TETRA'   :VTK_TETRA,    'HEXAHEDRON':VTK_HEXAHEDRON}
-        vtk_data.SetCells(cell_types[cell_type.split("_",1)[-1].upper()],cells)
+        # gracefully accept 'VTK_QUAD', 'vtk_lagrange_quadrilateral', etc.
+        cell_types = {'TRIA':VTK_LAGRANGE_TRIANGLE,    'QUAD':VTK_LAGRANGE_QUADRILATERAL,
+                      'TETR':VTK_LAGRANGE_TETRAHEDRON, 'HEXA':VTK_LAGRANGE_HEXAHEDRON}
+        vtk_data.SetCells(cell_types[cell_type.split('_')[-1].upper()[:4]],cells)
 
         return VTK(vtk_data)
 

python/tests/test_VTK.py

@@ -83,14 +83,17 @@ def test_polyData_VTKHDF(np_rng,tmp_path):
     vtkhdf = VTK.load(tmp_path/'polyData.vtkhdf')
     assert (v.as_ASCII() == vtkhdf.as_ASCII())
 
-@pytest.mark.parametrize('cell_type,n',[
-                                        ('VTK_hexahedron',8),
-                                        ('TETRA',4),
-                                        ('quad',4),
-                                        ('VTK_TRIANGLE',3)
-                                       ]
+# https://defelement.org/elements/lagrange.html
+@pytest.mark.parametrize('cell_type,n_nodes',[
+                                              ('VTK_hexahedron',lambda k: (k+1)**3),
+                                              ('TETRA',lambda k: ((k+1)*(k+2)*(k+3))//6),
+                                              ('quad',lambda k: (k+1)**2),
+                                              ('VTK_TRIANGLE',lambda k: ((k+1)*(k+2))//2),
+                                             ]
                         )
-def test_unstructuredGrid(np_rng,tmp_path,cell_type,n):
+@pytest.mark.parametrize('order',[1,2,3,4,5])
+def test_unstructuredGrid(np_rng,tmp_path,cell_type,n_nodes,order):
+    n = n_nodes(order)
     nodes = np_rng.random((n,3))
     connectivity = np_rng.choice(np.arange(n),n,False).reshape(-1,n)
     v = VTK.from_unstructured_grid(nodes,connectivity,cell_type)
@@ -104,15 +107,17 @@ def test_unstructuredGrid(np_rng,tmp_path,cell_type,n):
     assert (string == vtu.as_ASCII() == vtk.as_ASCII())
 
 @pytest.mark.xfail(np.lib.NumpyVersion(vtkVersion.GetVTKVersion()) < '9.4.0',
-                    reason = 'not available in VTK < 9.4')
-@pytest.mark.parametrize('cell_type,n',[
-                                        ('VTK_hexahedron',8),
-                                        ('TETRA',4),
-                                        ('quad',4),
-                                        ('VTK_TRIANGLE',3)
-                                       ]
+                   reason = 'not available in VTK < 9.4')
+@pytest.mark.parametrize('cell_type,n_nodes',[
+                                              ('VTK_lagrange_hexa',lambda k: (k+1)**3),
+                                              ('VTK_LAGRANGE_TETRAHEDRON',lambda k: ((k+1)*(k+2)*(k+3))//6),
+                                              ('quadrilateral',lambda k: (k+1)**2),
+                                              ('TRIANGLE',lambda k: ((k+1)*(k+2))//2),
+                                             ]
                         )
-def test_unstructuredGrid_VTKHDF(np_rng,tmp_path,cell_type,n):
+@pytest.mark.parametrize('order',[1,2,3,4,5])
+def test_unstructuredGrid_VTKHDF(np_rng,tmp_path,cell_type,n_nodes,order):
+    n = n_nodes(order)
     nodes = np_rng.random((n,3))
     connectivity = np_rng.choice(np.arange(n),n,False).reshape(-1,n)
     v = VTK.from_unstructured_grid(nodes,connectivity,cell_type)

src/Marc/element.f90

@@ -763,10 +763,10 @@ subroutine tElement_init(self,elemType)
       self%vtkType  = 'TRIANGLE'
     case(2)
       self%cellFace = CELLFACE2
-      self%vtkType  = 'QUAD'
+      self%vtkType  = 'QUADRILATERAL'
     case(3)
       self%cellFace = CELLFACE3
-      self%vtkType  = 'TETRA'
+      self%vtkType  = 'TETRAHEDRON'
     case(4)
       self%cellFace = CELLFACE4
       self%vtkType  = 'HEXAHEDRON'

src/result.f90

@@ -99,7 +99,7 @@ subroutine result_init(restart)
   if (.not. file_exists .or. .not. restart) then
     call result_createJobFile()
     call result_addAttribute('DADF5_version_major',1)
-    call result_addAttribute('DADF5_version_minor',2)
+    call result_addAttribute('DADF5_version_minor',3)
     call get_command_argument(0,commandLine)
     call result_addAttribute('creator',trim(commandLine)//' '//DAMASK_VERSION)
     call result_addAttribute('created',now())