Improve support for 2d meshes embedded in 3d space

examples/vector_fe/vector_fe_ex10/grad_div_exact_solution.h

@@ -29,34 +29,42 @@ class GradDivExactSolution
   GradDivExactSolution() = default;
   ~GradDivExactSolution() = default;
 
-  RealGradient operator() (Real x, Real y, Real z)
+  RealGradient operator() (Point p)
   {
-    libmesh_ignore(z);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
     const Real ux = cos(k*x)*sin(k*y);
     const Real uy = sin(k*x)*cos(k*y);
 
-    return RealGradient(ux, uy);
+    return R*RealGradient(ux, uy)/R.det();
   }
 
-  RealTensor grad(Real x, Real y, Real z)
+  RealTensor grad(Point p)
   {
-    libmesh_ignore(z);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
     const Real dux_dx = -k*sin(k*x)*sin(k*y);
     const Real dux_dy =  k*cos(k*x)*cos(k*y);
     const Real duy_dx = dux_dy;
     const Real duy_dy = dux_dx;
 
-    return RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy);
+    return R*RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy)*R.transpose()/R.det();
   }
 
-  RealGradient forcing(Real x, Real y, Real z)
+  RealGradient forcing(Point p)
   {
-    return (2*k*k + 1)*operator()(x, y, z);
+    return (2*k*k + 1)*operator()(p);
+  }
+
+  static void RM(RealTensor T)
+  {
+    R = T;
   }
 
 private:
+  static RealTensor R;
   const Real k = pi;
 };
 

examples/vector_fe/vector_fe_ex10/solution_function.h

@@ -45,10 +45,9 @@ class SolutionFunction : public FunctionBase<Number>
                            DenseVector<Number> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1), z=p(2);
-    output(0) = soln(x, y, z)(0);
-    output(1) = soln(x, y, z)(1);
-    output(2) = soln(x, y, z)(2);
+    output(0) = soln(p)(0);
+    output(1) = soln(p)(1);
+    output(2) = soln(p)(2);
   }
 
   virtual Number component(unsigned int component_in,
@@ -83,10 +82,9 @@ class SolutionGradient : public FunctionBase<Gradient>
                            DenseVector<Gradient> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1), z=p(2);
-    output(0) = soln.grad(x, y, z).row(0);
-    output(1) = soln.grad(x, y, z).row(1);
-    output(2) = soln.grad(x, y, z).row(2);
+    output(0) = soln.grad(p).row(0);
+    output(1) = soln.grad(p).row(1);
+    output(2) = soln.grad(p).row(2);
   }
 
   virtual Gradient component(unsigned int component_in,

examples/vector_fe/vector_fe_ex10/vector_fe_ex10.C

@@ -65,10 +65,12 @@
 #include "libmesh/getpot.h"
 #include "libmesh/exodusII_io.h"
 
-
 // Bring in everything from the libMesh namespace.
 using namespace libMesh;
 
+// Define static data member holding (optional) rotation matrix
+RealTensor GradDivExactSolution::R;
+
 // Function prototype. This is the function that will assemble
 // the linear system for our grad-div problem. Note that the
 // function will take the EquationSystems object and the
@@ -135,6 +137,12 @@ int main (int argc, char ** argv)
   // Make sure the code is robust against nodal reorderings.
   MeshTools::Modification::permute_elements(mesh);
 
+  // Make sure the code is robust against solves on 2d meshes rotated out of
+  // the xy plane. By default, all Euler angles are zero, the rotation matrix
+  // is the identity, and the mesh stays in place.
+  const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.);
+  GradDivExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi));
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
@@ -376,14 +384,9 @@ void assemble_graddiv(EquationSystems & es,
           // This involves a single loop in which we integrate the "forcing
           // function" in the PDE against the vector test functions (k).
           {
-            // The location of the current quadrature point.
-            const Real x = q_point[qp](0);
-            const Real y = q_point[qp](1);
-            const Real z = q_point[qp](2);
-
             // "f" is the forcing function, given by the well-known
             // "method of manufactured solutions".
-            RealGradient f = GradDivExactSolution().forcing(x, y, z);
+            RealGradient f = GradDivExactSolution().forcing(q_point[qp]);
 
             // Loop to integrate the vector test functions (k) against the
             // forcing function.
@@ -428,14 +431,8 @@ void assemble_graddiv(EquationSystems & es,
               // Loop over the face quadrature points for integration.
               for (unsigned int qp=0; qp<qface.n_points(); qp++)
                 {
-                  // The location on the boundary of the current
-                  // face quadrature point.
-                  const Real xf = qface_point[qp](0);
-                  const Real yf = qface_point[qp](1);
-                  const Real zf = qface_point[qp](2);
-
                   // The boundary value for the vector variable.
-                  RealGradient vector_value = GradDivExactSolution()(xf, yf, zf);
+                  RealGradient vector_value = GradDivExactSolution()(qface_point[qp]);
 
                   // We use the penalty method to set the flux of the vector
                   // variable at the boundary, i.e. the RT vector boundary dof.

examples/vector_fe/vector_fe_ex3/curl_curl_exact_solution.h

@@ -29,39 +29,43 @@ class CurlCurlExactSolution
   CurlCurlExactSolution() = default;
   ~CurlCurlExactSolution() = default;
 
-  RealGradient operator() (Real x, Real y)
+  RealGradient operator() (Point p)
   {
-    const Real ux =  cos(pi*x)*sin(pi*y);
-    const Real uy = -sin(pi*x)*cos(pi*y);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
 
-    return RealGradient(ux, uy);
+    const Real ux =  cos(k*x)*sin(k*y);
+    const Real uy = -sin(k*x)*cos(k*y);
+
+    return R*RealGradient(ux, uy);
   }
 
-  RealTensor grad(Real x, Real y)
+  RealTensor grad(Point p)
   {
-    const Real dux_dx = -pi*sin(pi*x)*sin(pi*y);
-    const Real dux_dy = pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dx = -pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dy = pi*sin(pi*x)*sin(pi*y);
+    Point pp = R.transpose()*p;
+    Real x = pp(0), y = pp(1);
+
+    const Real dux_dx = -k*sin(k*x)*sin(k*y);
+    const Real dux_dy =  k*cos(k*x)*cos(k*y);
+    const Real duy_dx = -k*cos(k*x)*cos(k*y);
+    const Real duy_dy =  k*sin(k*x)*sin(k*y);
 
-    return RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy);
+    return R*RealTensor(dux_dx, dux_dy, Real(0), duy_dx, duy_dy)*R.transpose();
   }
 
-  RealGradient curl(Real x, Real y)
+  RealGradient forcing(Point p)
   {
-    const Real dux_dy =  pi*cos(pi*x)*cos(pi*y);
-    const Real duy_dx = -pi*cos(pi*x)*cos(pi*y);
-
-    return RealGradient(Real(0), Real(0), duy_dx - dux_dy);
+    return (2*k*k + 1)*operator()(p);
   }
 
-  RealGradient forcing(Real x, Real y)
+  static void RM(RealTensor T)
   {
-    const Real fx =  (2*pi*pi + 1)*cos(pi*x)*sin(pi*y);
-    const Real fy = -(2*pi*pi + 1)*sin(pi*x)*cos(pi*y);
-
-    return RealGradient(fx, fy);
+    R = T;
   }
+
+private:
+  static RealTensor R;
+  const Real k = pi;
 };
 
 #endif // CURL_CURL_EXACT_SOLUTION_H

examples/vector_fe/vector_fe_ex3/curl_curl_system.C

@@ -151,7 +151,7 @@ bool CurlCurlSystem::element_time_derivative (bool request_jacobian,
       c.interior_curl(u_var, qp, curl_u);
 
       // Value of the forcing function at this quadrature point
-      RealGradient f = this->forcing(qpoint[qp]);
+      RealGradient f = soln.forcing(qpoint[qp]);
 
       // First, an i-loop over the degrees of freedom.
       for (unsigned int i=0; i != n_u_dofs; i++)
@@ -209,9 +209,9 @@ bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
       Gradient u;
       c.side_value(u_var, qp, u);
 
-      RealGradient N(normals[qp](0), normals[qp](1));
+      RealGradient N(normals[qp](0), normals[qp](1), normals[qp](2));
 
-      Gradient u_exact = this->exact_solution(qpoint[qp]);
+      Gradient u_exact = soln(qpoint[qp]);
 
       Gradient Ncu = (u - u_exact).cross(N);
 
@@ -227,17 +227,3 @@ bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
 
   return request_jacobian;
 }
-
-RealGradient CurlCurlSystem::forcing(const Point & p)
-{
-  Real x = p(0); Real y = p(1);
-
-  return soln.forcing(x, y);
-}
-
-RealGradient CurlCurlSystem::exact_solution(const Point & p)
-{
-  Real x = p(0); Real y = p(1);
-
-  return soln(x, y);
-}

examples/vector_fe/vector_fe_ex3/solution_function.h

@@ -47,19 +47,18 @@ class SolutionFunction : public FunctionBase<Number>
                            DenseVector<Number> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1);
     // libMesh assumes each component of the vector-valued variable is stored
     // contiguously.
-    output(_u_var)   = soln(x, y)(0);
-    output(_u_var+1) = soln(x, y)(1);
+    output(_u_var)   = soln(p)(0);
+    output(_u_var+1) = soln(p)(1);
+    output(_u_var+2) = soln(p)(2);
   }
 
   virtual Number component(unsigned int component_in,
                            const Point & p,
                            const Real)
   {
-    const Real x=p(0), y=p(1);
-    return soln(x, y)(component_in);
+    return soln(p)(component_in);
   }
 
   virtual std::unique_ptr<FunctionBase<Number>> clone() const
@@ -88,16 +87,15 @@ class SolutionGradient : public FunctionBase<Gradient>
                            DenseVector<Gradient> & output)
   {
     output.zero();
-    const Real x=p(0), y=p(1);
-    output(_u_var)   = soln.grad(x, y).row(0);
-    output(_u_var+1) = soln.grad(x, y).row(1);
+    output(_u_var)   = soln.grad(p).row(0);
+    output(_u_var+1) = soln.grad(p).row(1);
+    output(_u_var+2) = soln.grad(p).row(2);
   }
 
   virtual Gradient component(unsigned int component_in, const Point & p,
                              const Real)
   {
-    const Real x=p(0), y=p(1);
-    return soln.grad(x, y).row(component_in);
+    return soln.grad(p).row(component_in);
   }
 
   virtual std::unique_ptr<FunctionBase<Gradient>> clone() const

examples/vector_fe/vector_fe_ex3/vector_fe_ex3.C

@@ -42,10 +42,12 @@
 #include "libmesh/steady_solver.h"
 #include "solution_function.h"
 
-
 // Bring in everything from the libMesh namespace
 using namespace libMesh;
 
+// Define static data member holding (optional) rotation matrix
+RealTensor CurlCurlExactSolution::R;
+
 // The main program.
 int main (int argc, char ** argv)
 {
@@ -91,6 +93,12 @@ int main (int argc, char ** argv)
   // Make sure the code is robust against nodal reorderings.
   MeshTools::Modification::permute_elements(mesh);
 
+  // Make sure the code is robust against solves on 2d meshes rotated out of
+  // the xy plane. By default, all Euler angles are zero, the rotation matrix
+  // is the identity, and the mesh stays in place.
+  const Real phi = infile("phi", 0.), theta = infile("theta", 0.), psi = infile("psi", 0.);
+  CurlCurlExactSolution::RM(MeshTools::Modification::rotate(mesh, phi, theta, psi));
+
   // Print information about the mesh to the screen.
   mesh.print_info();
 
@@ -132,6 +140,7 @@ int main (int argc, char ** argv)
   // Print information about the system to the screen.
   equation_systems.print_info();
 
+  // Solve the system.
   system.solve();
 
   ExactSolution exact_sol(equation_systems);

include/fe/fe.h

@@ -389,7 +389,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
   static void nodal_soln(const Elem * elem, const Order o,
                          const std::vector<Number> & elem_soln,
                          std::vector<Number> & nodal_soln,
-                         bool add_p_level = true);
+                         bool add_p_level = true,
+                         const unsigned vdim = 1);
 
   /**
    * Build the nodal soln on one side from the (full) element soln.
@@ -401,7 +402,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
                               std::vector<Number> & nodal_soln_on_side,
-                              bool add_p_level = true);
+                              bool add_p_level = true,
+                              const unsigned vdim = 1);
 
   /**
    * \returns The number of shape functions associated with
@@ -791,7 +793,8 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
                                       const unsigned int side,
                                       const std::vector<Number> & elem_soln,
                                       std::vector<Number> & nodal_soln_on_side,
-                                      bool add_p_level = true);
+                                      bool add_p_level = true,
+                                      const unsigned vdim = 1);
 
   /**
    * An array of the node locations on the last

include/fe/fe_interface.h

@@ -289,7 +289,8 @@ class FEInterface
                          const Elem * elem,
                          const std::vector<Number> & elem_soln,
                          std::vector<Number> & nodal_soln,
-                         const bool add_p_level = true);
+                         const bool add_p_level = true,
+                         const unsigned int vdim = 1);
 
 
   /**
@@ -305,7 +306,8 @@ class FEInterface
                               const unsigned int side,
                               const std::vector<Number> & elem_soln,
                               std::vector<Number> & nodal_soln,
-                              const bool add_p_level = true);
+                              const bool add_p_level = true,
+                              const unsigned int vdim = 1);
 
 #ifdef LIBMESH_ENABLE_DEPRECATED
   /**