Merge branch '407-fix-deformed-dg-plotting' into 'feature'

[407-fix-deformed-dg-plotting] Fixed deforming grid

Closes #407

See merge request mpp/mpp!396

src/lib4_fem/plot/VtuPlot.cpp

@@ -86,6 +86,8 @@ void VtuPlot::pointDataStream(std::ostream &out, int k) {
     out << aArray;
   }
   out << ">" << endl;
+
+
   for (const auto &pData: pointData[plotMesh().t(k)]) {
     if (pData.first == "ID") continue;
 
@@ -101,6 +103,22 @@ void VtuPlot::pointDataStream(std::ostream &out, int k) {
     }
     out << endl << indentString(4) << "</DataArray>" << endl;
   }
+
+  if (isDeformed) {
+    out << indentString(4)
+        << R"(<DataArray type="Float32" Name="U)"
+        << R"(" NumberOfComponents="3)"
+        << R"(" format="ascii">)";
+    for (int i = 0; i < points[plotMesh().t(k)].size(); i++) {
+      for (int j = 0; j < 3; j++) {
+        Point displacement = displacements[plotMesh().t(k)].at(i);
+        out << lineBreakString(i + j, lineBreakZ, 5) << to_string(displacement[j]);
+      }
+
+    }
+    out << endl << indentString(4) << "</DataArray>" << endl;
+  }
+
   out << indentString(3) << "</PointData>" << endl;
 }
 
@@ -124,8 +142,8 @@ void VtuPlot::cellDataStream(std::ostream &out, int k) {
   if (!activeArrays.second.empty()) {
     string aArray = " ";
     aArray +=
-      (cellDataSize[plotMesh().slice(k)].at(activeArrays.second) > 1 ? R"(Vectors=")"
-                                                                     : R"(Scalars=")");
+        (cellDataSize[plotMesh().slice(k)].at(activeArrays.second) > 1 ? R"(Vectors=")"
+                                                                       : R"(Scalars=")");
     aArray += activeArrays.second;
     aArray += R"(")";
     out << aArray;
@@ -250,12 +268,12 @@ void VtuPlot::initCells(const Mesh &init) {
           try {
             if (i == 0) {
               cellConnectivity[center.t()] +=
-                lineBreakString(cells[center.t()].size() - 1, lineBreakC, 5);
+                  lineBreakString(cells[center.t()].size() - 1, lineBreakC, 5);
               cellConnectivity[center.t()] +=
-                to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
+                  to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
             } else {
               cellConnectivity[center.t()] +=
-                " " + to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
+                  " " + to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
             }
           }
           catch (...) {
@@ -264,11 +282,11 @@ void VtuPlot::initCells(const Mesh &init) {
         }
       }
       cellTypes[center.t()] +=
-        lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
+          lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
       cellTypes[center.t()] += to_string(VtkCellType(cType));
       totalOffset[center.t()] += m;
       cellOffset[center.t()] +=
-        lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
+          lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
       if (init.IsSTMesh())
         cellOffset[center.t()] += to_string(totalOffset[center.t()] / 2);
       else
@@ -336,7 +354,7 @@ void VtuPlot::addData(const string &name, const Vector &data, std::vector<int> i
     }
   }
   if (typeid(discToPlot) == typeid(LagrangeDiscretization) ||
-    typeid(discToPlot) == typeid(MultiPartDiscretization)) {
+      typeid(discToPlot) == typeid(MultiPartDiscretization)) {
     int d = ((const LagrangeDiscretization &) discToPlot).Degree();
     if (d > 0) addPointData(name, data, indices);
     else addCellData(name, data, indices);
@@ -345,25 +363,25 @@ void VtuPlot::addData(const string &name, const Vector &data, std::vector<int> i
   } else if (typeid(discToPlot) == typeid(DoFDiscretization<VertexDoF>)) {
     addPointData(name, data, indices);
   } else if (typeid(discToPlot) == typeid(DGDiscretization) ||
-    typeid(discToPlot) == typeid(DoFDiscretization<DGDoF>)) {
+             typeid(discToPlot) == typeid(DoFDiscretization<DGDoF>)) {
     addCellData(name, data, indices);
   } else if (typeid(discToPlot) == typeid(ArgyrisDiscretization)) {
     addArgyrisData(name, data, indices);
   } else if (typeid(discToPlot) == typeid(DivergenceFreeDiscretization)) {
     addDivergenceFreeData(name, data);
 #if USE_SPACETIME
-  } else if (typeid(discToPlot) == typeid(STDiscretization_DGDG)) {
-    addSpaceTimeData(name, data, indices);
-  } else if (typeid(discToPlot) == typeid(STDiscretization_DGDG_GLGL) ||
-    typeid(discToPlot) == typeid(STDiscretization_PGDG)) {
-    Warning("Plotting not implemented for"
-            " STDiscretization_DGDG_GLGL or"
-            " STDiscretization_PGDG")
+    } else if (typeid(discToPlot) == typeid(STDiscretization_DGDG)) {
+      addSpaceTimeData(name, data, indices);
+    } else if (typeid(discToPlot) == typeid(STDiscretization_DGDG_GLGL) ||
+      typeid(discToPlot) == typeid(STDiscretization_PGDG)) {
+      Warning("Plotting not implemented for"
+              " STDiscretization_DGDG_GLGL or"
+              " STDiscretization_PGDG")
 #endif
   } else {
     Warning("Plot type for Discretization "
-              + std::string(typeid(discToPlot).name())
-              + " not defined.")
+            + std::string(typeid(discToPlot).name())
+            + " not defined.")
   }
 }
 
@@ -515,6 +533,7 @@ void VtuPlot::DeformGrid(const Vector &displacement) {
   } else {
     deformLagrangeGrid(displacement);
   }
+  isDeformed = true;
 }
 
 void VtuPlot::Erase() {
@@ -536,7 +555,7 @@ void VtuPlot::ClearData() {
   dataToDrop.clear();
   for (const auto &cData: cellData[0.0]) {
     if (std::string{"IDSubdomainProcLoad"}.find(cData.first)
-      != std::string::npos)
+        != std::string::npos)
       continue;
     dataToDrop.push_back(cData.first);
   }
@@ -614,7 +633,7 @@ void VtuPlot::deformDGGrid(const Vector &displacement) {
       DGVectorFieldElement E(displacement, c);
       for (int i = 0; i < c.Corners(); ++i) {
         exBuffer.Send(0) << c.Corner(i);
-        VectorField u = E.VectorValue(c.Corner(i), displacement);
+        VectorField u = E.VectorValue(c.LocalCorner(i), displacement);
         for (int j = 0; j < dim; ++j)
           exBuffer.Send(0) << u[j];
       }
@@ -643,9 +662,9 @@ void VtuPlot::deformDGGrid(const Vector &displacement) {
     }
 
 //  for (auto &p: points) {
-  for (int i = 0; i < points.size(); i++) {
+  for (int i = 0; i < points[0.0].size(); i++) {
     displacements[0.0][i] =
-      deformationData[points[0.0][i]] / deformationCount[points[0.0][i]];
+        deformationData[points[0.0][i]] / deformationCount[points[0.0][i]];
   }
 //  }
 
@@ -678,7 +697,7 @@ void VtuPlot::deformLagrangeGrid(const Vector &displacement) {
       deformationData->insert({z, d});
     }
 
-  for (int i = 0; i < points.size(); i++) {
+  for (int i = 0; i < points[0.0].size(); i++) {
     displacements[0.0][i] = deformationData->at(points[0.0].at(i));
   }
 

src/lib4_fem/plot/VtuPlot.hpp

@@ -36,6 +36,7 @@ protected:
 
   std::map<double, std::vector<Point>> points{};
 
+  bool isDeformed{false};
   std::map<double, std::vector<Point>> displacements{};
 
   std::map<double, std::unordered_map<std::string, int>> cellDataSize{};