MeshTools::all_rbb() to generate IGA rational quadratic meshes

include/mesh/mesh_generation.h

@@ -114,8 +114,8 @@ void build_square (UnstructuredMesh & mesh,
  * Meshes a spherical or mapped-spherical domain.
  */
 void build_sphere (UnstructuredMesh & mesh,
-                   const Real rad=1,
-                   const unsigned int nr=2,
+                   const Real radius=1,
+                   const unsigned int n_refinements=2,
                    const ElemType type=INVALID_ELEM,
                    const unsigned int n_smooth=2,
                    const bool flat=true);

include/mesh/mesh_modification.h

@@ -141,6 +141,14 @@ void scale (MeshBase & mesh,
  */
 void all_tri (MeshBase & mesh);
 
+/**
+ * Converts all element geometric mappings from the default Lagrange
+ * to the more flexible Rational-Bezier-Bernstein.  When elements have
+ * curved edges and/or faces, node weights are chosen so that the new
+ * edges interpolate the old edge node locations with a circular arc.
+ */
+void all_rbb (MeshBase & mesh);
+
 /**
  * Smooth the mesh with a simple Laplace smoothing algorithm.  The mesh is
  * smoothed \p n_iterations times.  If the parameter \p power is 0, each

src/geom/edge_edge3.C

@@ -223,6 +223,10 @@ Edge3::second_order_child_vertex (const unsigned int) const
 
 Real Edge3::volume () const
 {
+  // This specialization is good for Lagrange mappings only
+  if (this->mapping_type() != LAGRANGE_MAP)
+    return this->Elem::volume();
+
   // Finding the (exact) length of a general quadratic element
   // is a surprisingly complicated formula.
   Point A = this->point(0) + this->point(1) - 2*this->point(2);

src/geom/edge_edge4.C

@@ -298,6 +298,10 @@ dof_id_type Edge4::key () const
 
 Real Edge4::volume () const
 {
+  // This specialization is good for Lagrange mappings only
+  if (this->mapping_type() != LAGRANGE_MAP)
+    return this->Elem::volume();
+
   // Make copies of our points.  It makes the subsequent calculations a bit
   // shorter and avoids dereferencing the same pointer multiple times.
   Point

src/mesh/mesh_modification.C

@@ -1452,6 +1452,336 @@ void MeshTools::Modification::all_tri (MeshBase & mesh)
 }
 
 
+
+void MeshTools::Modification::all_rbb (MeshBase & mesh)
+{
+  // By default, use 1.0 as the weight on every RATIONAL_BERNSTEIN
+  // mapped node
+  const Real default_weight = 1.0;
+
+  const auto weight_index =
+    (mesh.add_node_datum<Real>("rational_weight", true,
+                               &default_weight));
+
+  mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
+  mesh.set_default_mapping_data(weight_index);
+
+  // Out of loop to reduce heap allocations
+  std::unique_ptr<Elem> edge_ptr, face_ptr;
+
+  // Utility for checking extrusion directions later
+  auto almost_equal = [](Real a, Real b) {
+    return (std::abs(a-b) < TOLERANCE*TOLERANCE);
+  };
+
+  for (auto & elem : mesh.element_ptr_range())
+    {
+      if (elem->level())
+        libmesh_not_implemented_msg
+          ("all_rbb() currently only supports flat meshes with no refinement levels");
+
+#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
+      if (elem->infinite())
+        libmesh_not_implemented_msg
+          ("all_rbb() currently only supports finite geometric elements");
+#endif
+
+      elem->set_mapping_type(RATIONAL_BERNSTEIN_MAP);
+      elem->set_mapping_data(weight_index);
+
+      // Nothing to do unless we have curves to correct
+      if (elem->default_order() == FIRST)
+        continue;
+
+      // Modify the center node of an "edge" - possibly an actual
+      // edge, possibly a center node between edge points - to fit a
+      // circular curve.
+      auto make_edge_rbb = [default_weight, weight_index, almost_equal]
+        (const Node & n0, const Node & n1, Node & n_center)
+      {
+        // Skip edges we've already modified; the center node for
+        // these is no longer at the curve point we wish to
+        // interpolate, it should already be at the control point that
+        // accomplishes the interpolation.
+        const Real old_weight = n_center.get_extra_datum<Real>(weight_index);
+        if (old_weight != default_weight)
+          return;
+
+        const Point & p0 = n0;
+        const Point & p1 = n1;
+        Point & p2 = n_center;
+        const Point midchord = (p0+p1)/2;
+        const Real edge_chord_len_sq = (p1-p0).norm_sq();
+
+        const Real w0 = n0.get_extra_datum<Real>(weight_index);
+        const Real w1 = n1.get_extra_datum<Real>(weight_index);
+
+        // If we see edges that are unevenly parameterized, not just
+        // curved, I'm not sure what we want to do with those.  We
+        // can't isogeometrically represent a circular arc in this
+        // case unless we change the weights on the endpoints, but
+        // then that cascades to requiring changes in every other
+        // element sharing those endpoints.
+        //
+        // Presumably we want to maintain a somewhat similar uneven
+        // parameterization, for whatever boundary layer grading the
+        // mesh user wanted?  For now just scream and die.
+        const Point e02 = p2-p0,
+                    e21 = p1-p2;
+        const Real chord_02_len_sq = e02.norm_sq(),
+                   chord_21_len_sq = e21.norm_sq();
+        if (std::abs(chord_02_len_sq-chord_21_len_sq) > edge_chord_len_sq * TOLERANCE*TOLERANCE)
+          libmesh_not_implemented_msg
+            ("all_rbb() currently doesn't support unevenly parameterized edges");
+
+        // For straight edges we'll just want the middle node weight
+        // to match the endpoint nodes
+        const Point displacement_vec = p2-midchord;
+        if (displacement_vec.norm_sq() <
+            edge_chord_len_sq*TOLERANCE*TOLERANCE*TOLERANCE)
+        {
+          if (!almost_equal(w0, w1))
+            libmesh_not_implemented();
+
+          n_center.set_extra_datum<Real>(weight_index, w0);
+
+          return;
+        }
+
+        // Circularize the curve on curved edges
+        //
+        // First find the cosine of phi, the angle between our two
+        // subchords (turning from the direction of one to the
+        // direction of the other; this is the supplementary angle to
+        // the angle at the midpoint).  This is the same as half of
+        // the angle of our circular arc, which nicely enough is also
+        // the angle we take cos and sec of in NURBS formulae
+        const Real cos_phi = (e02*e21)/std::sqrt(chord_02_len_sq*chord_21_len_sq);
+        n_center.set_extra_datum<Real>(weight_index, cos_phi);
+
+        // There's a way to do really large arcs using negative
+        // weights, but we're going to get lousy approximation quality
+        // from isoparametric elements if we go too low, as well as
+        // bad numerics here, so let's just disallow it.
+        if (cos_phi < 0.5)
+          libmesh_not_implemented_msg
+            ("all_rbb() is not recommended for extremely sharp curves on one edge");
+
+        // Next find the radius of the circle our arc is from
+        const Real r = std::sqrt(edge_chord_len_sq)/2/std::sqrt(1-cos_phi*cos_phi);
+
+        // And perturb our center node so that it becomes a control
+        // point for that arc rather than a midpoint of that arc.
+        const Real R = r/cos_phi;
+
+        p2 += displacement_vec.unit() * (R-r);
+      };
+
+      auto make_face_rbb = [weight_index] (Elem & face)
+      {
+        // Prisms and pyramids may need to skip some faces while
+        // adjusting others
+        if (face.type() == TRI6)
+          return;
+
+        if (face.type() != QUAD9)
+          libmesh_not_implemented_msg
+            ("all_rbb() currently only supports mid-face nodes on Quad9 faces");
+
+        // We only use [4,8) but matching indices is nice and stack is
+        // cheap.
+        // that we're only using what we set, though.
+        Real w[9];
+
+        for (unsigned int i : make_range(4u, 8u))
+          w[i] = face.node_ref(i).get_extra_datum<Real>(weight_index);
+
+        // We can't currently handle arbitrary vertex weights
+#ifndef NDEBUG
+        for (unsigned int i : make_range(4u))
+          libmesh_assert_equal_to
+            (face.node_ref(i).get_extra_datum<Real>(weight_index), 1);
+#endif
+
+        // For the mid-face point, if we want to exactly match
+        // any cylinders and cones and spheres, we're actually already
+        // entirely constrained by the other points.
+        //
+        // This formula gives the minimum-energy Steiner surface based
+        // on the outer 8 points.
+        //
+        // That's an isogeometric representation of a cylinder aligned
+        // to either axis, or of a sphere where the quad edges are on
+        // latitude/longitude lines, or of a cone where two edges are
+        // segments of cone generating lines and the other two are
+        // arcs perpendicular to the axis.
+        //
+        // It's not perfectly isogeometric for the spheres we generate
+        // (where the quad edges are all great circles), but it should
+        // still converge asymptotically faster than non-rational
+        // quadratic Lagrange.
+        const Point xi_avg = (face.point(7) + face.point(5))/2;
+        const Point eta_avg = (face.point(4) + face.point(6))/2;
+        const Point vertex_avg = (face.point(0) + face.point(1) +
+                                  face.point(2) + face.point(3))/4;
+
+        const Real w_xi  = (w[7] + w[5])/2;
+        const Real w_eta = (w[4] + w[6])/2;
+        const Real w_mid = w_xi * w_eta;
+
+        Node & midnode = face.node_ref(8);
+        midnode.set_extra_datum<Real>(weight_index, w_mid);
+        midnode = ((1+w_mid)/(w_xi+w_eta) * (w_xi*xi_avg + w_eta*eta_avg) - vertex_avg)/w_mid;
+      };
+
+      // Check each edge for a curve, and adjust it if needed.
+      for (auto e : elem->edge_index_range())
+        {
+          elem->build_edge_ptr(edge_ptr, e);
+
+          // We should add EDGE4 once we have QUAD16/TRI10/HEX64 to
+          // use it
+          if (edge_ptr->type() != EDGE3)
+            libmesh_not_implemented_msg
+              ("all_rbb() currently only supports meshes with 2- and/or 3-node edges");
+
+          make_edge_rbb(edge_ptr->node_ref(0), edge_ptr->node_ref(1),
+                        edge_ptr->node_ref(2));
+
+        }
+
+      // If we're in 3D, we may have face nodes that also need to be
+      // adjusted to replace an interpolated curve with a spline
+      // curve.  We know what to do with a quad face, but we'll have
+      // to scream and die if we see a Tri7 face node.
+      bool check_face_points = (elem->dim() > 2) &&
+        (elem->n_nodes() > elem->n_edges() + elem->n_vertices());
+
+      if (check_face_points)
+        for (auto f : elem->side_index_range())
+          {
+            // Prisms and pyramids may need to skip some faces while
+            // adjusting others
+            if (elem->side_type(f) == TRI6)
+              continue;
+
+            elem->build_side_ptr(face_ptr, f);
+
+            make_face_rbb(*face_ptr);
+          }
+
+      bool check_interior_points =
+        elem->n_nodes() > elem->n_edges() + elem->n_vertices() + elem->n_faces();
+
+      if (check_interior_points)
+        {
+          if (elem->type() == EDGE3)
+            {
+              make_edge_rbb(elem->node_ref(0), elem->node_ref(1), elem->node_ref(2));
+            }
+          else if (elem->dim() == 2)
+            {
+              make_face_rbb(*elem);
+            }
+          else if (elem->type() == HEX27)
+            {
+              // For now we just support 2.5D (extrusions of 2D)
+              // meshes.  Let's look for an extrusion direction.
+              int extrude_dir = -1;
+
+              Real w[27];
+              for (unsigned int i : make_range(27u))
+                w[i] = elem->node_ref(i).get_extra_datum<Real>(weight_index);
+
+              if (almost_equal(w[0], w[1]) &&
+                  almost_equal(w[0], w[8]) &&
+                  almost_equal(w[2], w[3]) &&
+                  almost_equal(w[2], w[10]) &&
+                  almost_equal(w[4], w[5]) &&
+                  almost_equal(w[4], w[16]) &&
+                  almost_equal(w[6], w[7]) &&
+                  almost_equal(w[6], w[18]) &&
+                  almost_equal(w[9], w[11]) &&
+                  almost_equal(w[9], w[20]) &&
+                  almost_equal(w[12], w[13]) &&
+                  almost_equal(w[12], w[21]) &&
+                  almost_equal(w[14], w[15]) &&
+                  almost_equal(w[14], w[23]) &&
+                  almost_equal(w[17], w[19]) &&
+                  almost_equal(w[17], w[25]) &&
+                  almost_equal(w[22], w[24]))
+                extrude_dir = 0;
+
+              if (almost_equal(w[0], w[3]) &&
+                  almost_equal(w[0], w[11]) &&
+                  almost_equal(w[1], w[2]) &&
+                  almost_equal(w[1], w[9]) &&
+                  almost_equal(w[4], w[7]) &&
+                  almost_equal(w[4], w[19]) &&
+                  almost_equal(w[5], w[6]) &&
+                  almost_equal(w[5], w[17]) &&
+                  almost_equal(w[8], w[10]) &&
+                  almost_equal(w[8], w[20]) &&
+                  almost_equal(w[12], w[15]) &&
+                  almost_equal(w[12], w[24]) &&
+                  almost_equal(w[13], w[14]) &&
+                  almost_equal(w[13], w[22]) &&
+                  almost_equal(w[16], w[18]) &&
+                  almost_equal(w[16], w[25]) &&
+                  almost_equal(w[21], w[23]))
+                extrude_dir = 1;
+
+              if (almost_equal(w[0], w[4]) &&
+                  almost_equal(w[0], w[12]) &&
+                  almost_equal(w[1], w[5]) &&
+                  almost_equal(w[1], w[13]) &&
+                  almost_equal(w[2], w[6]) &&
+                  almost_equal(w[2], w[14]) &&
+                  almost_equal(w[3], w[7]) &&
+                  almost_equal(w[3], w[15]) &&
+                  almost_equal(w[8], w[16]) &&
+                  almost_equal(w[8], w[21]) &&
+                  almost_equal(w[9], w[17]) &&
+                  almost_equal(w[9], w[22]) &&
+                  almost_equal(w[10], w[18]) &&
+                  almost_equal(w[10], w[23]) &&
+                  almost_equal(w[11], w[19]) &&
+                  almost_equal(w[11], w[24]) &&
+                  almost_equal(w[20], w[25]))
+                extrude_dir = 2;
+
+              // If we're extruding, then we can derive the center point
+              // values from the two perpendicular faces' centers.
+              Node & n_center = elem->node_ref(26);
+              Point & p_center = n_center;
+              if (extrude_dir == 0)
+                {
+                  n_center.set_extra_datum<Real>(weight_index, w[22]);
+                  p_center = (elem->point(22) + elem->point(24))/2;
+                }
+              else if (extrude_dir == 1)
+                {
+                  n_center.set_extra_datum<Real>(weight_index, w[21]);
+                  p_center = (elem->point(21) + elem->point(23))/2;
+                }
+              else if (extrude_dir == 2)
+                {
+                  n_center.set_extra_datum<Real>(weight_index, w[20]);
+                  p_center = (elem->point(20) + elem->point(25))/2;
+                }
+              else
+                libmesh_not_implemented_msg
+                  ("all_rbb() currently only supports 2.5D (extrusions) meshes in 3D");
+            }
+          else
+            libmesh_not_implemented_msg
+              ("all_rbb() doesn't yet support " << elem->type());
+        }
+    }
+}
+
+
+
 void MeshTools::Modification::smooth (MeshBase & mesh,
                                       const unsigned int n_iterations,
                                       const Real power)

tests/Makefile.am

@@ -56,6 +56,7 @@ unit_tests_sources = \
   geom/side_vertex_average_normal_test.C \
   geom/which_node_am_i_test.C \
   mesh/all_second_order.C \
+  mesh/all_rbb.C \
   mesh/all_tri.C \
   mesh/distort.C \
   mesh/boundary_mesh.C \