WIP: Refactor Core FEM

data/meshes/one_tri.mesh

@@ -0,0 +1,29 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+1
+1 2 0 1 2
+
+boundary
+0
+
+vertices
+3
+2
+0 0
+1 0
+0 1

data/meshes/two_hexes.mesh

@@ -0,0 +1,39 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+2
+1 5 0 1 3 2 4 5 7 6
+1 5 4 5 7 6 8 9 11 10
+
+boundary
+0
+
+vertices
+12
+3
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1
+0 0 2
+1 0 2
+0 1 2
+1 1 2

data/meshes/two_tets.mesh

@@ -0,0 +1,32 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+3
+
+elements
+2
+1 4 0 1 2 3
+1 4 1 2 3 4
+
+boundary
+0
+
+vertices
+5
+3
+0 0 0
+1 0 0
+0 1 0
+0 0 1
+1 1 1

data/meshes/two_tris.mesh

@@ -0,0 +1,35 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+2
+1 2 0 1 2
+1 2 1 3 2
+
+boundary
+4
+1 1 0 1
+1 1 1 3
+1 1 3 2
+1 1 2 0
+
+vertices
+4
+2
+0 0
+1 0
+0 1
+1 1

src/serac/numerics/CMakeLists.txt

@@ -6,6 +6,8 @@
 
 add_subdirectory(functional)
 
+add_subdirectory(refactor)
+
 set(numerics_headers
     equation_solver.hpp
     odes.hpp

src/serac/numerics/functional/detail/hexahedron_H1.inl

@@ -209,7 +209,7 @@ struct finite_element<mfem::Geometry::CUBE, H1<p, c> > {
     union {
       tensor<qf_input_type, q * q * q>                                   one_dimensional;
       tensor<tuple<tensor<double, c>, tensor<double, c, dim> >, q, q, q> three_dimensional;
-    } output;
+    } output{};
 
     for (int qz = 0; qz < q; qz++) {
       for (int qy = 0; qy < q; qy++) {

src/serac/numerics/functional/detail/metaprogramming.hpp

@@ -98,3 +98,16 @@ SERAC_HOST_DEVICE constexpr void for_constexpr(const lambda& f)
 {
   detail::for_constexpr(f, std::make_integer_sequence<int, n>{}...);
 }
+
+
+namespace impl {
+  template < auto x >
+  struct value{
+    constexpr operator decltype(x)() { return x; }
+  };
+}
+
+template < auto ... args, typename T >
+void foreach_constexpr(T && function) {
+  (function(impl::value<args>{}), ...);
+}

src/serac/numerics/functional/detail/quadrilateral_H1.inl

@@ -226,7 +226,7 @@ struct finite_element<mfem::Geometry::SQUARE, H1<p, c> > {
     union {
       tensor<qf_input_type, q * q>                                    one_dimensional;
       tensor<tuple<tensor<double, c>, tensor<double, c, dim> >, q, q> two_dimensional;
-    } output;
+    } output{};
 
     for (int qy = 0; qy < q; qy++) {
       for (int qx = 0; qx < q; qx++) {

src/serac/numerics/functional/domain.hpp

@@ -182,6 +182,7 @@ struct Domain {
     exit(1);
   }
 
+
   /**
    * @brief returns how many elements of any type belong to this domain
    */
@@ -246,6 +247,28 @@ struct Domain {
                    mfem::Geometry::Type element_geometry);
 };
 
+inline uint32_t spatial_dimension(const Domain & domain) {
+  return static_cast<uint32_t>(domain.mesh_.SpaceDimension());
+}
+
+inline uint32_t geometry_dimension(const Domain & domain) {
+  return static_cast<uint32_t>(domain.dim_);
+}
+
+enum class DomainType { ISOPARAMETRIC, SPATIAL };
+
+struct DomainWithType {
+  const Domain & domain;
+  const DomainType type;
+
+  DomainWithType(const Domain & d) : domain(d), type(DomainType::SPATIAL) {};
+  DomainWithType(DomainType t, const Domain & d) : domain(d), type(t) {};
+};
+
+inline DomainWithType isoparametric(const Domain & domain) {
+  return DomainWithType(DomainType::ISOPARAMETRIC, domain);
+}
+
 /// @brief constructs a domain from all the elements in a mesh
 Domain EntireDomain(const mesh_t& mesh);
 

src/serac/numerics/functional/element_restriction.hpp

@@ -61,7 +61,7 @@ struct DoF {
   DoF() : bits{} {}
 
   /// copy ctor
-  DoF(const DoF& other) : bits{other.bits} {}
+  //DoF(const DoF& other) : bits{other.bits} {}
 
   /// create a `DoF` from the given index, sign and orientation values
   DoF(uint64_t index, uint64_t sign = 0, uint64_t orientation = 0)
@@ -70,7 +70,7 @@ struct DoF {
   }
 
   /// copy assignment operator
-  void operator=(const DoF& other) { bits = other.bits; }
+  //void operator=(const DoF& other) { bits = other.bits; }
 
   /// get the sign field of this `DoF`
   int sign() const { return (bits & sign_mask) ? -1 : 1; }

src/serac/numerics/functional/family.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+namespace serac {
+
+/**
+ * @brief Element conformity
+ *
+ * QOI   denotes a "quantity of interest", implying integration with the test function "1"
+ * H1    denotes a function space where values are continuous across element boundaries
+ * HCURL denotes a vector-valued function space where only the tangential component is continuous across element
+ * boundaries HDIV  denotes a vector-valued function space where only the normal component is continuous across element
+ * boundaries L2    denotes a function space where values are discontinuous across element boundaries
+ */
+enum class Family
+{
+  QOI,
+  H1,
+  HCURL,
+  HDIV,
+  L2
+};
+
+}

src/serac/numerics/functional/finite_element.hpp

@@ -12,6 +12,7 @@
  */
 #pragma once
 
+#include "family.hpp"
 #include "tuple.hpp"
 #include "tensor.hpp"
 #include "geometry.hpp"
@@ -169,23 +170,7 @@ SERAC_HOST_DEVICE constexpr int elements_per_block(int q)
   }
 }
 
-/**
- * @brief Element conformity
- *
- * QOI   denotes a "quantity of interest", implying integration with the test function "1"
- * H1    denotes a function space where values are continuous across element boundaries
- * HCURL denotes a vector-valued function space where only the tangential component is continuous across element
- * boundaries HDIV  denotes a vector-valued function space where only the normal component is continuous across element
- * boundaries L2    denotes a function space where values are discontinuous across element boundaries
- */
-enum class Family
-{
-  QOI,
-  H1,
-  HCURL,
-  HDIV,
-  L2
-};
+
 
 /**
  * @brief H1 elements of order @p p