starting refactor of ref_elem to support 1D, 2D, and 3D ref spaces

.gitmodules

@@ -1,4 +1,4 @@
 [submodule "MATAR"]
 	path = matar
 	url = https://github.com/lanl/MATAR.git
-	branch = main
+	branch = main

examples/CMakeLists.txt

@@ -4,4 +4,5 @@ add_subdirectory(average)
 add_subdirectory(decomp_example)
 add_subdirectory(point_connectivity)
 add_subdirectory(mesh_form)
+add_subdirectory(reference_element)
 # add_subdirectory(other_example) # Future examples go here

examples/average/include/mesh_io.h

@@ -142,7 +142,7 @@ inline int PointIndexFromIJK(int i, int j, int k, const int* order)
 ///
 /////////////////////////////////////////////////////////////////////////////
 void build_3d_box(
-    swage::Mesh& mesh,
+    swage::Mesh_t& mesh,
     node_t&   node,
     double origin[3],
     double length[3],
@@ -212,7 +212,7 @@ void build_3d_box(
     node.coords.update_host();
 
     // initialize elem variables
-    mesh.initialize_elems(num_elems, num_dim);
+    mesh.initialize_elems(num_elems);
 
     // populate the point data structures
     FOR_ALL(k, 0, num_elems_k,
@@ -268,7 +268,7 @@ void build_3d_box(
 /// \param rank rank
 ///
 /////////////////////////////////////////////////////////////////////////////
-    void write_vtk(swage::Mesh& mesh,
+    void write_vtk(swage::Mesh_t& mesh,
         node_t& node,
         int rank)
     {
@@ -526,7 +526,7 @@ void build_3d_box(
 /// \param comm MPI communicator
 ///
 /////////////////////////////////////////////////////////////////////////////
-void write_vtu(swage::Mesh& mesh,
+void write_vtu(swage::Mesh_t& mesh,
                node_t& node,
                GaussPoint_t& gauss_point,
                int rank,
@@ -843,7 +843,7 @@ void write_vtu(swage::Mesh& mesh,
     /// \param Number of dimensions
     ///
     /////////////////////////////////////////////////////////////////////////////
-    void read_vtk_mesh(swage::Mesh& mesh,
+    void read_vtk_mesh(swage::Mesh_t& mesh,
         node_t&   node,
         int num_dims,
         std::string mesh_file_)
@@ -942,7 +942,7 @@ void write_vtu(swage::Mesh& mesh,
         printf("Number of elements read in %zu\n", num_elem);
 
         // initialize elem variables
-        mesh.initialize_elems(num_elem, num_dims);
+        mesh.initialize_elems(num_elem);
 
         found=true;
     } // end if
@@ -1001,9 +1001,9 @@ void write_vtu(swage::Mesh& mesh,
     mesh.nodes_in_elem.update_device();
 
 
-    // initialize corner variables
-    size_t num_corners = num_elem * num_nodes_in_elem;
-    mesh.initialize_corners(num_corners);
+    // initialize corner variables (it is set in the elems initialization)
+    //size_t num_corners = num_elem * num_nodes_in_elem;
+    //mesh.initialize_corners(num_corners);
 
 
     // Build connectivity

examples/average/src/average.cpp

@@ -44,14 +44,18 @@ ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "state.h"
 
 
+
 int main(int argc, char** argv) {
 
 MATAR_INITIALIZE(argc, argv);
 { // MATAR scope
     std::cout<<"Hello, Average Example!"<<std::endl;
 
+
+    printf("\nAverage test finished.\n");
+
 } // end MATAR scope
 MATAR_FINALIZE();
 
 return 0;
-}
+}

examples/decomp_example/include/mesh_io.h

@@ -164,7 +164,7 @@ int PointIndexFromIJK(int i, int j, int k, const int* order)
 ///
 /////////////////////////////////////////////////////////////////////////////
 void build_3d_box(
-    swage::Mesh& mesh,
+    swage::Mesh_t& mesh,
     MPICArrayKokkos<double>& node_coords,
     double origin[3],
     double length[3],
@@ -243,11 +243,11 @@ void build_3d_box(
     node_coords.update_host();
 
     // initialize elem variables
-    if (Pn_order == 0){
-        mesh.initialize_elems(num_elems, num_dim); 
+    if (Pn_order == 0){ // BUG: Pn_order=1 is needed for coding above to work correctly
+        mesh.initialize_elems(num_elems); 
         Pn_order = 1;
     } else {
-        mesh.initialize_elems_Pn(num_elems, num_dim, Pn_order);
+        mesh.initialize_elems_Pn(num_elems, Pn_order, 2*Pn_order);
     }
 
     // populate the point data structures
@@ -321,7 +321,7 @@ void build_3d_box(
 ///
 /////////////////////////////////////////////////////////////////////////////
 void build_2d_polar(
-    swage::Mesh& mesh,
+    swage::Mesh_t& mesh,
     MPICArrayKokkos<double>& node_coords,
     double& inner_radius,
     double& outer_radius,
@@ -398,7 +398,7 @@ void build_2d_polar(
     node_coords.update_device();
 
     // initialize elem variables
-    mesh.initialize_elems(num_elems, num_dim);
+    mesh.initialize_elems(num_elems);
 
     // populate the elem center data structures
     for (int j = 0; j < num_elems_j; j++) {
@@ -432,9 +432,9 @@ void build_2d_polar(
     // update device side
     mesh.nodes_in_elem.update_device();
 
-    // intialize corner variables
-    int num_corners = num_elems * mesh.num_nodes_in_elem;
-    mesh.initialize_corners(num_corners);
+    // intialize corner variables (this is initialized in initialize_elems)
+    //int num_corners = num_elems * mesh.num_nodes_in_elem;
+    //mesh.initialize_corners(num_corners);
     // corner.initialize(num_corners, num_dim);
 
     // Build connectivity
@@ -455,7 +455,7 @@ void build_2d_polar(
 /// \param comm MPI communicator
 ///
 /////////////////////////////////////////////////////////////////////////////
-void write_vtu(swage::Mesh& mesh,
+void write_vtu(swage::Mesh_t& mesh,
                node_t& node,
                GaussPoint_t& gauss_point,
                int rank,

examples/decomp_example/src/mesh_decomp_example.cpp

@@ -59,11 +59,13 @@ int main(int argc, char** argv) {
     int num_elems_j = 40;
 
     // Initial mesh built on rank zero
-    swage::Mesh initial_mesh;
+    swage::Mesh_t initial_mesh;
+    initial_mesh.num_dims = num_dims;
     MPICArrayKokkos<double> initial_node_coords;
 
     // Mesh partitioned by pt-scotch, including ghost
-    swage::Mesh final_mesh;
+    swage::Mesh_t final_mesh;
+    final_mesh.num_dims = num_dims;
     node_t final_node;
     MPICArrayKokkos<double> final_node_coords;
 

examples/mesh_form/include/mesh_io.h

@@ -164,7 +164,7 @@ int PointIndexFromIJK(int i, int j, int k, const int* order)
 ///
 /////////////////////////////////////////////////////////////////////////////
 void build_3d_box(
-    swage::Mesh& mesh,
+    swage::Mesh_t& mesh,
     DCArrayKokkos<double>& node_coords,
     double origin[3],
     double length[3],
@@ -251,9 +251,9 @@ void build_3d_box(
 
     // initialize elem variables
     if (is_linear){
-        mesh.initialize_elems(num_elems, num_dim); 
+        mesh.initialize_elems(num_elems); 
     } else {
-        mesh.initialize_elems_Pn(num_elems, num_dim, Pn_order);
+        mesh.initialize_elems_Pn(num_elems, Pn_order, 2*Pn_order);
     }
 
     // populate the point data structures
@@ -310,8 +310,9 @@ void build_3d_box(
 
     // initialize corner variables (corner = element-node pair)
     // used for per-corner data like corner_delta in mesh_mold.cpp
-    int num_corners = num_elems * mesh.num_nodes_in_elem;
-    mesh.initialize_corners(num_corners);
+    // the corners are initialized in initialize_elems
+    //int num_corners = num_elems * mesh.num_nodes_in_elem;
+    //mesh.initialize_corners(num_corners);
 
     // Build connectivity
     mesh.build_connectivity();
@@ -332,7 +333,7 @@ void build_3d_box(
 ///
 /////////////////////////////////////////////////////////////////////////////
 void build_2d_polar(
-    swage::Mesh& mesh,
+    swage::Mesh_t& mesh,
     MPICArrayKokkos<double>& node_coords,
     double& inner_radius,
     double& outer_radius,
@@ -409,7 +410,7 @@ void build_2d_polar(
     node_coords.update_device();
 
     // initialize elem variables
-    mesh.initialize_elems(num_elems, num_dim);
+    mesh.initialize_elems(num_elems);
 
     // populate the elem center data structures
     for (int j = 0; j < num_elems_j; j++) {
@@ -444,8 +445,9 @@ void build_2d_polar(
     mesh.nodes_in_elem.update_device();
 
     // intialize corner variables
-    int num_corners = num_elems * mesh.num_nodes_in_elem;
-    mesh.initialize_corners(num_corners);
+    // this is initialized in initialize_elems
+    //int num_corners = num_elems * mesh.num_nodes_in_elem;
+    //mesh.initialize_corners(num_corners);
     // corner.initialize(num_corners, num_dim);
 
     // Build connectivity
@@ -466,7 +468,7 @@ void build_2d_polar(
 /// \param comm MPI communicator
 ///
 /////////////////////////////////////////////////////////////////////////////
-void write_vtu(swage::Mesh& mesh,
+void write_vtu(swage::Mesh_t& mesh,
                node_t& node,
                GaussPoint_t& gauss_point,
                int rank,
@@ -864,7 +866,7 @@ void write_vtu(swage::Mesh& mesh,
 /// \param Number of dimensions
 ///
 /////////////////////////////////////////////////////////////////////////////
-void read_vtk_mesh(swage::Mesh& mesh,
+void read_vtk_mesh(swage::Mesh_t& mesh,
     DCArrayKokkos<double>& node_coords,
     int num_dims,
     std::string mesh_file_)
@@ -964,7 +966,7 @@ found=false;
                     throw std::invalid_argument("Failed to parse CELLS header line: \"" + line + "\"");
                 }
                 printf("Number of elements read in %zu\n", num_elem);
-                mesh.initialize_elems(num_elem, num_dims);
+                mesh.initialize_elems(num_elem);
                 found = true;
                 break;
             }
@@ -1097,8 +1099,9 @@ mesh.nodes_in_elem.update_device();
 
 
 // initialize corner variables
-size_t num_corners = num_elem * num_nodes_in_elem;
-mesh.initialize_corners(num_corners);
+// this is initialized in initialize_elems
+//size_t num_corners = num_elem * num_nodes_in_elem;
+//mesh.initialize_corners(num_corners);
 
 
 // Build connectivity

examples/mesh_form/src/mesh_form.cpp

@@ -68,7 +68,7 @@ int main(int argc, char** argv) {
 
 
     // Initial mesh built on rank zero
-    swage::Mesh mesh;
+    swage::Mesh_t mesh;
     node_t node;
     DCArrayKokkos<double> node_coords;
 

examples/reference_element/CMakeLists.txt

@@ -0,0 +1,37 @@
+# examples/average/CMakeLists.txt
+
+# 1. Define the executable
+#    Point directly to the source file inside 'src/'
+add_executable(remap_fv_test src/remap_fv_test.cpp)
+add_executable(remap_dg_test src/remap_dg_test.cpp)
+add_executable(ref_flux_test src/ref_flux_test.cpp)
+add_executable(ref_plus_mesh_test src/ref_plus_mesh_test.cpp)
+add_executable(integration_test src/integration_test.cpp)
+add_executable(gradient_test src/gradient_test.cpp)
+add_executable(interpolation_test src/interpolation_test.cpp)
+add_executable(kronecker_delta_test src/kronecker_delta_test.cpp)
+add_executable(partition_unity_test src/partition_unity_test.cpp)
+
+# 2. Add this example's specific include path
+#    This allows main.cpp to find headers in examples/point_connectivity/include/
+target_include_directories(remap_fv_test PRIVATE include)
+target_include_directories(remap_dg_test PRIVATE include)
+target_include_directories(ref_flux_test PRIVATE include)
+target_include_directories(ref_plus_mesh_test PRIVATE include)
+target_include_directories(integration_test PRIVATE include)
+target_include_directories(gradient_test PRIVATE include)
+target_include_directories(interpolation_test PRIVATE include)
+target_include_directories(kronecker_delta_test PRIVATE include)
+target_include_directories(partition_unity_test PRIVATE include)
+
+# 3. Link against the main library (ELEMENTS)
+#    This pulls in Kokkos, MATAR, and the main library headers automatically.
+target_link_libraries(remap_fv_test PRIVATE ELEMENTS)
+target_link_libraries(remap_dg_test PRIVATE ELEMENTS)
+target_link_libraries(ref_flux_test PRIVATE ELEMENTS)
+target_link_libraries(ref_plus_mesh_test PRIVATE ELEMENTS)
+target_link_libraries(integration_test PRIVATE ELEMENTS)
+target_link_libraries(gradient_test PRIVATE ELEMENTS)
+target_link_libraries(interpolation_test PRIVATE ELEMENTS)
+target_link_libraries(kronecker_delta_test PRIVATE ELEMENTS)
+target_link_libraries(partition_unity_test PRIVATE ELEMENTS)