Second order writeMatplotlib

examples/refinement/second_order_demo.cc

@@ -43,6 +43,10 @@ int main() {
           TikzOutputCtrl::RenderCells | TikzOutputCtrl::CellNumbering |
               TikzOutputCtrl::VerticeNumbering | TikzOutputCtrl::NodeNumbering |
               TikzOutputCtrl::EdgeNumbering | TikzOutputCtrl::SecondOrder);
+      lf::io::writeMatplotlib(*mesh,
+                              mesh_name.substr(0, mesh_name.find_last_of('.')) +
+                                  "_" + std::to_string(step) + ".txt",
+                              true);
       multi_mesh.RefineRegular();
     }
   }

lib/lf/io/plot_mesh.py

@@ -4,42 +4,51 @@
 import numpy as np
 
 if len(argv) < 2:
-    print('usage: python plot_mesh.py mesh.csv')
+    print('usage: python plot_mesh.py writeMatplotlib_output.txt')
     exit(-1)
 
-vertices = []
-segments = []
-triangles = []
-quads = []
+points = dict()
+segments = dict()
+triangles = dict()
+quadrilaterals = dict()
 
 with open(argv[1]) as f:
     for line in f:
-        array = np.fromstring(line.split()[0], sep=',')
-
-        if array[0] == 0:
-            if array.size == 5:
-                triangles.append(array)
-            elif array.size == 6:
-                quads.append(array)
-            else:
-                raise ValueError('triangles and quadrilaterals only')
-        if array[0] == 1:
-            segments.append(array)
-        if array[0] == 2:
-            vertices.append(array)
-
-vertices = np.vstack(vertices)
-segments = np.vstack(segments)
-if triangles != []:
-    triangles = np.vstack(triangles)
-if quads != []:
-    quads = np.vstack(quads)
+        data = line.rstrip().split(' ')
+        geometry_type = data[0]
+        index = int(data[1])
+        coordinates = [
+            np.array([float(x), float(y)])
+            for x, y in zip(data[2::2], data[3::2])
+        ]
+
+        if geometry_type == 'Point':
+            geometry_list = points
+        elif geometry_type.startswith('Segment'):
+            geometry_list = segments
+        elif geometry_type.startswith('Tria'):
+            geometry_list = triangles
+        elif geometry_type.startswith('Quad'):
+            geometry_list = quadrilaterals
+        else:
+            raise RuntimeError('Unknown geometry')
+
+        geometry_list[index] = coordinates
+
+x_min, x_max = float('inf'), float('-inf')
+y_min, y_max = float('inf'), float('-inf')
 
 # plot vertices
-for i, num in enumerate(vertices[:, 1]):
+for idx, coords in points.items():
+    for coord in coords:
+        x_min = min(coord[0], x_min)
+        x_max = max(coord[0], x_max)
+        y_min = min(coord[1], y_min)
+        y_max = max(coord[1], y_max)
+
     plt.annotate(
-        int(num),
-        vertices[np.where(vertices[:, 1] == num)][0, -2:],
+        idx,
+        coords[0],
         fontsize='xx-large',
         weight='bold',
         bbox=dict(boxstyle='circle', facecolor='none', edgecolor='red'),
@@ -48,56 +57,96 @@
         va='center'
     )
 
-# plot segments
-for segment in segments:
-    x_coords, y_coords = np.column_stack((
-        vertices[np.where(vertices[:, 1] == int(segment[2]))][0, -2:],
-        vertices[np.where(vertices[:, 1] == int(segment[3]))][0, -2:]
-    ))
-    plt.plot(x_coords, y_coords, 'k-')
-
-    midpoint = .5 * (
-            vertices[np.where(vertices[:, 1] == segment[2])][0, -2:] +
-            vertices[np.where(vertices[:, 1] == segment[3])][0, -2:]
-    )
 
-    plt.annotate(
-        int(segment[1]),
-        midpoint,
-        fontsize='xx-large',
-        ha='center',
-        va='center'
-    )
+def poly_coeffs(x, coeffs):
+    order = len(coeffs)
+    y = 0
+
+    for i in range(order):
+        y += coeffs[i] * x ** (order - (i + 1))
+
+    return y
+
+
+def collinear(p0, p1, p2):
+    x1, y1 = p1[0] - p0[0], p1[1] - p0[1]
+    x2, y2 = p2[0] - p0[0], p2[1] - p0[1]
+    return abs(x1 * y2 - x2 * y1) < 1e-12
 
 
 # plot cells
-def plot_cells(cells):
-    for cell in cells:
-        cell_vertices = []
-        for i in range(2, cell.size):
-            segment = segments[np.where(segments[:, 1] == cell[i])][0]
-            cell_vertices.append(
-                vertices[np.where(vertices[:, 1] == segment[2])][0, -2:]
+for idx, coords in segments.items():
+    for coord in coords:
+        x_min = min(coord[0], x_min)
+        x_max = max(coord[0], x_max)
+        y_min = min(coord[1], y_min)
+        y_max = max(coord[1], y_max)
+
+    # SegmentO1
+    if len(coords) == 2:
+        coordinates = np.vstack(coords)
+        x_coords = coordinates[:, 0]
+        y_coords = coordinates[:, 1]
+        plt.plot(x_coords, y_coords, 'k-')
+
+        plt.annotate(
+            idx,
+            [x_coords.mean(), y_coords.mean()],
+            fontsize='xx-large',
+            ha='center',
+            va='center'
+        )
+
+    # SegmentO2
+    elif len(coords) == 3:
+        vertex_0 = coords[0]
+        vertex_1 = coords[1]
+        midpoint = coords[2]
+
+        # check if points are collinear
+        if collinear(vertex_0, midpoint, vertex_1):
+            plt.plot(
+                [vertex_0[0], vertex_1[0]], [vertex_0[1], vertex_1[1]], 'k-'
             )
-            cell_vertices.append(
-                vertices[np.where(vertices[:, 1] == segment[3])][0, -2:]
+        else:
+            coefficients = np.polyfit(
+                [vertex_0[0], midpoint[0], vertex_1[0]],
+                [vertex_0[1], midpoint[1], vertex_1[1]],
+                2
             )
-
-        cell_vertices = np.vstack(cell_vertices)
+            x = np.linspace(vertex_0[0], vertex_1[0], 1000)
+            plt.plot(x, poly_coeffs(x, coefficients), 'k-')
 
         plt.annotate(
-            int(cell[1]),
-            np.mean(cell_vertices, axis=0),
+            idx,
+            midpoint,
             fontsize='xx-large',
-            color='deeppink',
-            weight='bold',
             ha='center',
             va='center'
         )
+# plot cells
+cell_data = {'limegreen': triangles, 'm': quadrilaterals}
 
+for color, cells in cell_data.items():
+    for idx, coords in cells.items():
+        for coord in coords:
+            x_min = min(coord[0], x_min)
+            x_max = max(coord[0], x_max)
+            y_min = min(coord[1], y_min)
+            y_max = max(coord[1], y_max)
 
-plot_cells(triangles)
-plot_cells(quads)
+        plt.annotate(
+            idx,
+            np.vstack(coords).mean(0),
+            fontsize='xx-large',
+            color=color,
+            weight='bold',
+            ha='center',
+            va='center'
+        )
 
+offset = 1e-2
+plt.xlim(left=x_min - offset, right=x_max + offset)
+plt.ylim(bottom=y_min - offset, top=y_max + offset)
 plt.axis('off')
 plt.show()

lib/lf/io/write_matplotlib.cc

@@ -9,16 +9,16 @@
 #include "write_matplotlib.h"
 
 #include <fstream>
-#include <iostream>
 
 namespace lf::io {
 
-void writeMatplotlib(const lf::mesh::Mesh &mesh, std::string filename) {
+void writeMatplotlib(const lf::mesh::Mesh &mesh, std::string filename,
+                     bool second_order) {
   using dim_t = lf::base::RefEl::dim_t;
 
-  // append .csv to filename if necessary
-  if (filename.compare(filename.size() - 4, 4, ".csv") != 0) {
-    filename += ".csv";
+  // append .txt to filename if necessary
+  if (filename.compare(filename.size() - 4, 4, ".txt") != 0) {
+    filename += ".txt";
   }
 
   std::ofstream file(filename);
@@ -28,50 +28,110 @@ void writeMatplotlib(const lf::mesh::Mesh &mesh, std::string filename) {
     LF_VERIFY_MSG(dim_mesh == 2,
                   "write_matplotlib() only available for 2D meshes");
 
-    // loop through all elements of every codimension
-    for (int codim = 0; codim <= dim_mesh; ++codim) {
-      for (const lf::mesh::Entity &obj : mesh.Entities(codim)) {
-        const size_t obj_idx = mesh.Index(obj);
-        const lf::base::RefEl obj_ref_el = obj.RefEl();
-        const lf::geometry::Geometry *obj_geo_ptr = obj.Geometry();
-        const Eigen::MatrixXd vertices =
-            obj_geo_ptr->Global(obj_ref_el.NodeCoords());
-
-        switch (obj_ref_el) {
-          case lf::base::RefEl::kPoint(): {
-            file << codim << ',' << obj_idx << ',' << vertices(0, 0) << ','
-                 << vertices(1, 0) << std::endl;
-            break;
-          }
-          case lf::base::RefEl::kSegment(): {
-            file << codim << ',' << obj_idx << ',';
-            // to access points of segment use SubEntities(1)
-            for (const auto &sub : obj.SubEntities(codim)) {
-              file << mesh.Index(sub) << ',';
-            }
-            file << std::endl;
-
-            break;
-          }
-          case lf::base::RefEl::kTria():
-          case lf::base::RefEl::kQuad(): {
-            file << codim << ',' << obj_idx << ',';
-            // to access points of cell use SubEntities(1)
-            for (const auto &sub : obj.SubEntities(codim + 1)) {
-              file << mesh.Index(sub) << ',';
-            }
-            file << std::endl;
-
-            break;
-          }
-          default: {
-            LF_VERIFY_MSG(false, "Error for object " + std::to_string(obj_idx) +
-                                     " in codim " + std::to_string(codim) +
-                                     " of type " + obj_ref_el.ToString());
-            break;
-          }
+    // store points to file
+    {
+      std::vector<std::pair<size_t, Eigen::Vector2d>> points;
+
+      for (const lf::mesh::Entity &point : mesh.Entities(2)) {
+        points.emplace_back(std::make_pair(
+            mesh.Index(point),
+            point.Geometry()->Global(point.RefEl().NodeCoords())));
+      }
+
+      for (const auto &point : points) {
+        const Eigen::Vector2d coords = point.second;
+        file << "Point"
+             << " " << point.first << " " << coords(0) << " " << coords(1)
+             << std::endl;
+      }
+    }
+
+    // store segments to file
+    {
+      std::vector<std::pair<size_t, Eigen::Matrix<double, 2, Eigen::Dynamic>>>
+          segments;
+
+      for (const lf::mesh::Entity &segment : mesh.Entities(1)) {
+        if (second_order) {
+          segments.emplace_back(std::make_pair(
+              mesh.Index(segment),
+              segment.Geometry()->Global(
+                  (Eigen::MatrixXd(1, 3) << 0, 1, 0.5).finished())));
+        } else {
+          segments.emplace_back(
+              std::make_pair(mesh.Index(segment),
+                             segment.Geometry()->Global(
+                                 (Eigen::MatrixXd(1, 2) << 0, 1).finished())));
+        }
+      }
+
+      for (const auto &segment : segments) {
+        file << (second_order ? "SegmentO2" : "SegmentO1") << " "
+             << segment.first;
+        const Eigen::Matrix<double, 2, Eigen::Dynamic> coords = segment.second;
+
+        for (int col = 0; col < coords.cols(); ++col) {
+          file << " " << coords(0, col) << " " << coords(1, col);
+        }
+
+        file << std::endl;
+      }
+    }
+
+    {
+      std::vector<std::pair<size_t, Eigen::Matrix<double, 2, Eigen::Dynamic>>>
+          cells;
+
+      for (const lf::mesh::Entity &cell : mesh.Entities(0)) {
+        const lf::geometry::Geometry *geometry_ptr = cell.Geometry();
+        const Eigen::MatrixXd local_vertices = cell.RefEl().NodeCoords();
+        const Eigen::MatrixXd vertices = geometry_ptr->Global(local_vertices);
+
+        if (!second_order) {
+          cells.emplace_back(std::make_pair(mesh.Index(cell), vertices));
+
+        } else {
+          // compute midpoints
+          const int num_points = vertices.cols();
+          Eigen::MatrixXd lower_diagonal =
+              Eigen::MatrixXd::Zero(num_points, num_points);
+          lower_diagonal.diagonal<-1>() = Eigen::VectorXd::Ones(num_points - 1);
+          lower_diagonal(0, num_points - 1) = 1;
+
+          const Eigen::MatrixXd local_midpoints =
+              local_vertices *
+              (.5 * (lower_diagonal +
+                     Eigen::MatrixXd::Identity(num_points, num_points)));
+          const Eigen::MatrixXd midpoints =
+              geometry_ptr->Global(local_midpoints);
+
+          cells.emplace_back(std::make_pair(
+              mesh.Index(cell),
+              (Eigen::MatrixXd(2, 2 * num_points) << vertices, midpoints)
+                  .finished()));
         }
       }
+
+      for (const auto &cell : cells) {
+        const Eigen::Matrix<double, 2, Eigen::Dynamic> coords = cell.second;
+
+        if (coords.cols() == 3 || coords.cols() == 6) {
+          file << (second_order ? "TriaO2" : "TriaO1");
+        } else if (coords.cols() == 4 || coords.cols() == 8) {
+          file << (second_order ? "QuadO2" : "QuadO1");
+
+        } else {
+          LF_VERIFY_MSG(false, "Unknown cell geometry");
+        }
+
+        file << " " << cell.first;
+
+        for (int col = 0; col < coords.cols(); ++col) {
+          file << " " << coords(0, col) << " " << coords(1, col);
+        }
+
+        file << std::endl;
+      }
     }
   }
 }