Merge pull request #3397 from roystgnr/added_sides_fixes

Fix for added_sides, more robust HEX+SIDE_HIERARCHIC, additional unit tests

Author: roystgnr

src/fe/fe_hierarchic_shape_3D.C

@@ -874,6 +874,9 @@ Real FE<3,SIDE_HIERARCHIC>::shape(const Elem * elem,
         libmesh_assert_less(i, 6*dofs_per_side);
 
         const unsigned int sidenum = cube_side(p);
+        if (sidenum > 5)
+          return std::numeric_limits<Real>::quiet_NaN();
+
         const unsigned int dof_offset = sidenum * dofs_per_side;
 
         if (i < dof_offset) // i is on a previous side
@@ -1095,6 +1098,9 @@ Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
         libmesh_assert_less(i, 6*dofs_per_side);
 
         const unsigned int sidenum = cube_side(p);
+        if (sidenum > 5)
+          return std::numeric_limits<Real>::quiet_NaN();
+
         const unsigned int dof_offset = sidenum * dofs_per_side;
 
         if (i < dof_offset) // i is on a previous side
@@ -1346,6 +1352,9 @@ Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
         libmesh_assert_less(i, 6*dofs_per_side);
 
         const unsigned int sidenum = cube_side(p);
+        if (sidenum > 5)
+          return std::numeric_limits<Real>::quiet_NaN();
+
         const unsigned int dof_offset = sidenum * dofs_per_side;
 
         if (i < dof_offset) // i is on a previous side
@@ -1578,8 +1587,10 @@ unsigned int cube_side (const Point & p)
     return 4;
   else if (zeta > maxabs_xi_eta)
     return 5;
-  else
-    libmesh_error_msg("Cannot determine side to evaluate");
+
+  // We need to be able to return invalid values for cases where
+  // mixed FE are being evaluated together on edges and vertices
+  return 65535;
 }
 
 

src/systems/equation_systems.C

@@ -900,13 +900,13 @@ EquationSystems::build_parallel_solution_vector(const std::set<std::string> * sy
                                 {
                                   // Retrieve index into global solution vector.
                                   std::size_t node_index =
-                                    libmesh_map_find(discontinuous_node_indices,
-                                                     std::make_tuple(elem->id(), s, n));
+                                    nv * libmesh_map_find(discontinuous_node_indices,
+                                                          std::make_tuple(elem->id(), s, n));
 
                                   for (unsigned int d=0; d < n_vec_dim; d++)
                                     {
-                                      parallel_soln.add(node_index*n_vec_dim+d, nodal_soln[n_vec_dim*n+d]);
-                                      repeat_count.add(node_index*n_vec_dim+d, 1);
+                                      parallel_soln.add(node_index + (var_inc+d + var_num), nodal_soln[n_vec_dim*n+d]);
+                                      repeat_count.add(node_index + (var_inc+d + var_num), 1);
                                     }
                                 }
                             }

tests/mesh/mesh_input.C

@@ -114,14 +114,18 @@ public:
 #endif // !LIBMESH_USE_COMPLEX_NUMBERS
 
   CPPUNIT_TEST( testExodusWriteAddedSidesEdgeC0 );
+  CPPUNIT_TEST( testExodusWriteAddedSidesMixedEdgeC0 );
   // CPPUNIT_TEST( testExodusWriteAddedSidesEdgeDisc ); // need is_on_face fixes
   CPPUNIT_TEST( testExodusWriteAddedSidesTriC0 );
+  CPPUNIT_TEST( testExodusWriteAddedSidesMixedTriC0 );
   // CPPUNIT_TEST( testExodusWriteAddedSidesTriDisc ); // Need aligned faces
   CPPUNIT_TEST( testExodusWriteAddedSidesQuadC0 );
+  CPPUNIT_TEST( testExodusWriteAddedSidesMixedQuadC0 );
   // CPPUNIT_TEST( testExodusWriteAddedSidesQuadDisc ); // need is_on_face fixes
   // CPPUNIT_TEST( testExodusWriteAddedSidesTetC0 ); // BROKEN!?!  WHY!?!
   // CPPUNIT_TEST( testExodusWriteAddedSidesTetDisc );
   CPPUNIT_TEST( testExodusWriteAddedSidesHexC0 );
+  CPPUNIT_TEST( testExodusWriteAddedSidesMixedHexC0 );
   CPPUNIT_TEST( testExodusWriteAddedSidesHexDisc );
 
   CPPUNIT_TEST( testExodusFileMappingsPlateWithHole);
@@ -748,7 +752,9 @@ public:
     (Number (*exact_sol)(const Point &, const Parameters &, const
                          std::string &, const std::string &),
      const ElemType elem_type,
-     const Order order)
+     const Order order,
+     const std::vector<FEType> earlier_vars = {},
+     const std::vector<FEType> later_vars = {})
   {
     constexpr unsigned int nx = added_sides_nxyz[0],
                            ny = added_sides_nxyz[1],
@@ -770,8 +776,16 @@ public:
 
       EquationSystems es(mesh);
       System & sys = es.add_system<System> ("SimpleSystem");
+      int varnum = 1;
+      for (auto vartype : earlier_vars)
+        sys.add_variable("earlier_"+std::to_string(varnum++), vartype);
+
       sys.add_variable("u", order, SIDE_HIERARCHIC);
 
+      varnum = 1;
+      for (auto vartype : later_vars)
+        sys.add_variable("later_"+std::to_string(varnum++), vartype);
+
       if (dim == 3)
         MeshTools::Generation::build_cube
           (mesh, nx, ny, nz, 0., 1., 0., 1., 0., 1., elem_type);
@@ -965,6 +979,12 @@ public:
     testExodusWriteAddedSides(six_x_plus_sixty_y, EDGE3, SECOND);
   }
 
+  void testExodusWriteAddedSidesMixedEdgeC0()
+  {
+    testExodusWriteAddedSides(six_x_plus_sixty_y, EDGE3, FIRST, {{FIRST, LAGRANGE}});
+    testExodusWriteAddedSides(six_x_plus_sixty_y, EDGE3, SECOND, {}, {{FIRST, LAGRANGE}});
+  }
+
   void testExodusWriteAddedSidesEdgeDisc()
   {
     testExodusWriteAddedSides(designed_for_side_elems, EDGE3, SECOND);
@@ -976,6 +996,12 @@ public:
     testExodusWriteAddedSides(six_x_plus_sixty_y, TRI6, SECOND);
   }
 
+  void testExodusWriteAddedSidesMixedTriC0()
+  {
+    testExodusWriteAddedSides(six_x_plus_sixty_y, TRI6, FIRST, {{SECOND, HIERARCHIC}});
+    testExodusWriteAddedSides(six_x_plus_sixty_y, TRI6, SECOND, {}, {{SECOND, SZABAB}});
+  }
+
   void testExodusWriteAddedSidesTriDisc()
   {
     testExodusWriteAddedSides(designed_for_side_elems, TRI6, SECOND);
@@ -987,6 +1013,12 @@ public:
     testExodusWriteAddedSides(six_x_plus_sixty_y, QUAD9, SECOND);
   }
 
+  void testExodusWriteAddedSidesMixedQuadC0()
+  {
+    testExodusWriteAddedSides(six_x_plus_sixty_y, QUAD9, FIRST, {{SECOND, LAGRANGE}});
+    testExodusWriteAddedSides(six_x_plus_sixty_y, QUAD9, SECOND, {}, {{FIRST, LAGRANGE}});
+  }
+
   void testExodusWriteAddedSidesQuadDisc()
   {
     testExodusWriteAddedSides(designed_for_side_elems, QUAD9, SECOND);
@@ -1009,6 +1041,12 @@ public:
     testExodusWriteAddedSides(six_x_plus_sixty_y, HEX27, SECOND);
   }
 
+  void testExodusWriteAddedSidesMixedHexC0()
+  {
+    testExodusWriteAddedSides(six_x_plus_sixty_y, HEX27, FIRST, {{FIRST, LAGRANGE}});
+    testExodusWriteAddedSides(six_x_plus_sixty_y, HEX27, SECOND, {}, {{SECOND, HIERARCHIC}});
+  }
+
   void testExodusWriteAddedSidesHexDisc()
   {
     testExodusWriteAddedSides(designed_for_side_elems, HEX27, SECOND);