fixed boundary plot for EG (is like Lagrange plot)

src/elasticity/assemble/EGElasticity.cpp

@@ -927,30 +927,54 @@ void EGElasticity::GetInvariant(const Vector &displacement, Vector &iota4) const
   iota4.Collect();
 }
 
-/*void EGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
-    if (vtkplot < 1) return;
-    Vector bnd(0.0, U.GetDisc());
-
-    auto &plot = mpp::plot("Boundary");
-    for (cell c = U.cells(); c != U.cells_end(); c++) {
-        MixedEGVectorFieldElement E(U, *c);
-        RowBndValues u_c(bnd, *c);
-        if (!u_c.onBnd()) continue;
-
-        for (int face = 0; face < c.Faces(); ++face) {
-            MixedEGVectorFieldFaceElement FE(U, *c, face);
-            int shape_size = FE.dimension()-1;
-            int bc = u_c.bc(face);
-
-            for (int i = 0; i < shape_size; ++i) {
-                bnd(c(), i) = max(bnd(c(), i), (double) bc);
-            }
-
+void EGElasticity::PlotBoundary(const Vector &U, const Vector &scal) const {
+  if (vtkplot < 1) return;
+  Vector bnd(scal);
+  bnd.Clear();
+  Vector count(bnd);
+  count.Clear();
+  auto &plot = mpp::plot("Boundary");
+  for (cell c = U.cells(); c != U.cells_end(); c++) {
+    RowBndValues u_c(bnd, *c);
+    RowBndValues u_count(count, *c);
+
+    if (!u_c.onBnd()) continue;
+
+    for (int i = 0; i < c.Faces(); ++i) {
+      for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
+        int id = U.IdOfNodalPointOnFace(*c, i, j);
+        if (u_c.bc(i) > u_c(id, 0)) {
+          u_c(id, 0) = (double) u_c.bc(i);
+          u_count(id, 0) = 1;
         }
+      }
     }
-    plot.AddData("Boundary", bnd);
-    plot.PlotFile("Boundary");
-}*/
+  }
+  bnd.Accumulate();
+  count.Accumulate();
+  for (size_t i = 0; i < count.size(); i++) {
+    if (count[i] > 0) {
+      bnd[i] /= count[i];
+    }
+  }
+  plot.AddData("Boundary", bnd, 0);
+
+  bnd = 0;
+  for (cell c = U.cells(); c != U.cells_end(); c++) {
+    RowBndValues u_c(bnd, *c);
+    if (!u_c.onBnd()) continue;
+
+    for (int i = 0; i < c.Faces(); ++i) {
+      for (int j = 0; j < U.NumberOfNodalPointsOnFace(*c, i); ++j) {
+        int id = U.IdOfNodalPointOnFace(*c, i, j);
+        u_c(id, 0) = max(u_c(id, 0), u_c.bc(i) == 199 ? 1.0 : 0.0);
+      }
+    }
+  }
+  plot.AddData("Fixation", bnd);
+
+  plot.PlotFile("Boundary");
+}
 
 void EGElasticity::SetDisplacement(Vector &U) const {
   for (row r = U.rows(); r != U.rows_end(); ++r) {

src/elasticity/assemble/EGElasticity.hpp

@@ -92,7 +92,7 @@ public:
 
   double NormPrestress();
 
-  /*void PlotBoundary(const Vector &U, const Vector &scal) const override;*/
+  void PlotBoundary(const Vector &U, const Vector &scal) const override;
 
   /*std::array<double, 4> CalculateVolume(const Vector &U) const override;*/