Merge topic 'FixLagrangePointParametricCoordinates'

26d3405 add multi degree support to vtkUpdateCellsV8toV9
c610e5b fix: overwrite PointParametricCoordinates instead of writting at its tail
1736482 fix PointParametricCoordinates not updated when Lagrange cell degrees change

Acked-by: Kitware Robot <[email protected]>
Acked-by: T.J. Corona <[email protected]>
Merge-request: !6929

Author: tjcorona

Common/DataModel/vtkHigherOrderHexahedron.cxx

@@ -669,17 +669,21 @@ int vtkHigherOrderHexahedron::PointIndexFromIJK(int i, int j, int k, const int*
 }
 
 vtkIdType vtkHigherOrderHexahedron::NodeNumberingMappingFromVTK8To9(
-  const vtkIdType order, const vtkIdType node_id_vtk8)
+  const int order[3], const vtkIdType node_id_vtk8)
 {
-  const int numPtsPerEdgeWithoutCorners = static_cast<int>(order) - 1;
-
-  int offset = 8 + 10 * numPtsPerEdgeWithoutCorners;
-  if ((node_id_vtk8 < offset) || (node_id_vtk8 >= offset + 2 * numPtsPerEdgeWithoutCorners))
+  int numPtsPerEdgeWithoutCorners[3];
+  numPtsPerEdgeWithoutCorners[0] = order[0] - 1;
+  numPtsPerEdgeWithoutCorners[1] = order[1] - 1;
+  numPtsPerEdgeWithoutCorners[2] = order[2] - 1;
+
+  int offset = 8 + 4 * (numPtsPerEdgeWithoutCorners[0] + numPtsPerEdgeWithoutCorners[1]) +
+    2 * numPtsPerEdgeWithoutCorners[2];
+  if ((node_id_vtk8 < offset) || (node_id_vtk8 >= offset + 2 * numPtsPerEdgeWithoutCorners[2]))
     return node_id_vtk8;
-  else if (node_id_vtk8 < offset + numPtsPerEdgeWithoutCorners)
-    return node_id_vtk8 + numPtsPerEdgeWithoutCorners;
+  else if (node_id_vtk8 < offset + numPtsPerEdgeWithoutCorners[2])
+    return node_id_vtk8 + numPtsPerEdgeWithoutCorners[2];
   else
-    return node_id_vtk8 - numPtsPerEdgeWithoutCorners;
+    return node_id_vtk8 - numPtsPerEdgeWithoutCorners[2];
 }
 
 /**\brief Given the index, \a subCell, of a linear approximating-hex, translate pcoords from that
@@ -764,6 +768,8 @@ void vtkHigherOrderHexahedron::SetUniformOrderFromNumPoints(vtkIdType numPts)
 
 void vtkHigherOrderHexahedron::SetOrder(int s, int t, int u)
 {
+  if (this->PointParametricCoordinates && (Order[0] != s || Order[1] != t || Order[2] != u))
+    this->PointParametricCoordinates->Reset();
   Order[0] = s;
   Order[1] = t;
   Order[2] = u;

Common/DataModel/vtkHigherOrderHexahedron.h

@@ -104,7 +104,7 @@ class VTKCOMMONDATAMODEL_EXPORT vtkHigherOrderHexahedron : public vtkNonLinearCe
   virtual vtkHigherOrderInterpolation* getInterp() = 0;
 
   static vtkIdType NodeNumberingMappingFromVTK8To9(
-    const vtkIdType order, const vtkIdType node_id_vtk8);
+    const int order[3], const vtkIdType node_id_vtk8);
 
 protected:
   vtkHigherOrderHexahedron();

Common/DataModel/vtkHigherOrderInterpolation.cxx

@@ -1044,14 +1044,13 @@ void vtkHigherOrderInterpolation::AppendCurveCollocationPoints(
     pts = vtkSmartPointer<vtkPoints>::New();
   }
 
-  vtkIdType existing = pts->GetNumberOfPoints();
   vtkIdType np = order[0] + 1;
-  pts->SetNumberOfPoints(existing + np);
+  pts->SetNumberOfPoints(np);
   vtkVector3d e0(0., 0., 0.);
   vtkVector3d e1(+1., 0., 0.);
 
   // Insert corner points
-  vtkIdType sn = existing;
+  vtkIdType sn = 0;
   pts->SetPoint(sn++, e0.GetData());
   pts->SetPoint(sn++, e1.GetData());
 
@@ -1072,11 +1071,10 @@ void vtkHigherOrderInterpolation::AppendQuadrilateralCollocationPoints(
     pts = vtkSmartPointer<vtkPoints>::New();
   }
 
-  vtkIdType existing = pts->GetNumberOfPoints();
   vtkIdType np = (order[0] + 1) * (order[1] + 1);
-  pts->SetNumberOfPoints(existing + np);
+  pts->SetNumberOfPoints(np);
   // Insert corner points
-  vtkIdType sn = existing;
+  vtkIdType sn = 0;
   for (int ii = 0; ii < 4; ++ii)
   {
     vtkVector3d cc(hexCorner[ii]);
@@ -1118,11 +1116,10 @@ void vtkHigherOrderInterpolation::AppendHexahedronCollocationPoints(
     pts = vtkSmartPointer<vtkPoints>::New();
   }
 
-  vtkIdType existing = pts->GetNumberOfPoints();
   vtkIdType np = (order[0] + 1) * (order[1] + 1) * (order[2] + 1);
-  pts->SetNumberOfPoints(existing + np);
+  pts->SetNumberOfPoints(np);
   // Insert corner points
-  vtkIdType sn = existing;
+  vtkIdType sn = 0;
   for (int ii = 0; ii < 8; ++ii)
   {
     pts->SetPoint(sn++, hexCorner[ii]);
@@ -1182,12 +1179,11 @@ void vtkHigherOrderInterpolation::AppendWedgeCollocationPoints(
     pts = vtkSmartPointer<vtkPoints>::New();
   }
 
-  vtkIdType existing = pts->GetNumberOfPoints();
   vtkIdType np =
     (order[0] + 1) * (order[1] + 2) * (order[2] + 1) / 2; // NB: assert(order[0] == order[1])
-  pts->SetNumberOfPoints(existing + np);
+  pts->SetNumberOfPoints(np);
   // Insert corner points
-  vtkIdType sn = existing;
+  vtkIdType sn = 0;
   for (int ii = 0; ii < 6; ++ii)
   {
     pts->SetPoint(sn++, wedgeCorner[ii]);

Common/DataModel/vtkHigherOrderQuadrilateral.cxx

@@ -591,6 +591,8 @@ void vtkHigherOrderQuadrilateral::SetUniformOrderFromNumPoints(const vtkIdType n
 
 void vtkHigherOrderQuadrilateral::SetOrder(const int s, const int t)
 {
+  if (this->PointParametricCoordinates && (Order[0] != s || Order[1] != t))
+    this->PointParametricCoordinates->Reset();
   Order[0] = s;
   Order[1] = t;
   Order[2] = (s + 1) * (t + 1);

Common/DataModel/vtkHigherOrderWedge.cxx

@@ -1239,6 +1239,9 @@ void vtkHigherOrderWedge::SetOrder(const int s, const int t, const int u, const
 {
   if (s != t)
     vtkErrorMacro("For wedges, the first two degrees should be equals.");
+
+  if (this->PointParametricCoordinates && (Order[0] != s || Order[2] != u))
+    this->PointParametricCoordinates->Reset();
   Order[0] = s;
   Order[1] = s;
   Order[2] = u;

IO/Core/vtkUpdateCellsV8toV9.h

@@ -21,6 +21,7 @@
 #define vtkUpdateCellsV8toV9_h
 
 #include "vtkCellArray.h"
+#include "vtkCellData.h"
 #include "vtkCellType.h"
 #include "vtkCellTypes.h"
 #include "vtkHigherOrderHexahedron.h"
@@ -41,11 +42,26 @@ inline void vtkUpdateCellsV8toV9(vtkUnstructuredGrid* output)
     {
       output->GetCells()->GetCellAtId(i, oldpts);
       newpts->DeepCopy(oldpts);
-      int order =
-        static_cast<int>(round(std::cbrt(static_cast<int>(oldpts->GetNumberOfIds())))) - 1;
+
+      int degs[3];
+      if (output->GetCellData()->SetActiveAttribute(
+            "HigherOrderDegrees", vtkDataSetAttributes::AttributeTypes::HIGHERORDERDEGREES) != -1)
+      {
+        vtkDataArray* v = output->GetCellData()->GetHigherOrderDegrees();
+        double degs_double[3];
+        v->GetTuple(i, degs_double);
+        for (int ii = 0; ii < 3; ii++)
+          degs[ii] = static_cast<int>(degs_double[ii]);
+      }
+      else
+      {
+        int order =
+          static_cast<int>(round(std::cbrt(static_cast<int>(oldpts->GetNumberOfIds())))) - 1;
+        degs[0] = degs[1] = degs[2] = order;
+      }
       for (int j = 0; j < oldpts->GetNumberOfIds(); j++)
       {
-        int newid = vtkHigherOrderHexahedron::NodeNumberingMappingFromVTK8To9(order, j);
+        int newid = vtkHigherOrderHexahedron::NodeNumberingMappingFromVTK8To9(degs, j);
         if (j != newid)
         {
           newpts->SetId(j, oldpts->GetId(newid));