Fixed lifting test.

CMakeLists.txt

@@ -121,6 +121,7 @@ endif()
 if (OPT_ENABLE_GPU_SOLVER)
     option(OPT_DEBUG_GPU "Enable GPU debugging" OFF)
     if (OPT_DEBUG_GPU)
+        add_definitions(-DCUDA_DEBUG)
         add_definitions(-DGPU_DEBUG)
     endif()
 endif()
@@ -212,19 +213,17 @@ if (COMPILER_IS_CLANG OR COMPILER_IS_GNU)
     endif()
 endif()
 
-option(OPT_ENABLE_CUDA_DEBUG_GUARDS "Enable CUDA debug guards")
-if (OPT_ENABLE_CUDA_DEBUG_GUARDS AND HAVE_CUDA)
-    add_definitions(-DCUDA_DEBUG)
+option(OPT_WRITE_TEST_OUTPUTS "Produce SILO files in tests" OFF)
+if (OPT_WRITE_TEST_OUTPUTS)
+    add_definitions(-DWRITE_TEST_OUTPUTS)
 endif()
 
+option(OPT_EXPENSIVE_ASSERTS "Enable expensive assert()s" OFF)
+if (OPT_EXPENSIVE_ASSERTS)
+    add_definitions(-DEXPENSIVE_ASSERTS)
+endif()
 ######################################################################
 ## Enable/disable solver modules
-
-option(OPT_ENABLE_CPU_SOLVER "Enable CPU DG solver")
-if (OPT_ENABLE_CPU_SOLVER)
-    add_definitions(-DENABLE_CPU_SOLVER)
-endif()
-
 option(OPT_DEBUG_PRINT "Print debug information")
 if (OPT_DEBUG_PRINT)
     add_definitions(-DENABLE_DEBUG_PRINT)
@@ -261,6 +260,9 @@ target_link_libraries(test_mass gmshdg)
 add_executable(test_differentiation tests/test_differentiation.cpp)
 target_link_libraries(test_differentiation gmshdg)
 
+add_executable(test_lifting tests/test_lifting.cpp)
+target_link_libraries(test_lifting gmshdg)
+
 add_executable(test_differentiation_gpu tests/test_differentiation_gpu.cpp)
 target_link_libraries(test_differentiation_gpu gmshdg)
 

src/entity.cpp

@@ -532,6 +532,9 @@ entity::populate_lifting_matrices_planar(entity_data_cpu& ed,
         size_t base_row = iO * ed.num_bf;
         ed.lifting_matrices.block(base_row, 0, ed.num_bf, lm_cols) =
             mass_ldlt.solve(lm);
+
+        //std::cout << "Orientation : " << iO << std::endl;
+        //std::cout << mass_ldlt.solve(lm) << std::endl;
     }
 }
 
@@ -703,6 +706,7 @@ entity::populate_entity_data(entity_data_cpu& ed) const
     {
         auto& pe = physical_cells[iT];
         auto cell_keys = pe.bf_keys();
+        auto cell_base = iT * ed.num_bf;
 
         for (size_t iF = 0; iF < 4; iF++)
         {
@@ -710,8 +714,7 @@ entity::populate_entity_data(entity_data_cpu& ed) const
             auto& pf = physical_faces[face_num];
             auto face_keys = pf.bf_keys();
 
-            auto face_base = (4*iT + iF) * ed.num_fluxes;
-            auto cell_base = iT * ed.num_bf;
+            auto face_base = face_num * ed.num_fluxes;
 
             for (size_t i = 0; i < face_keys.size(); i++)
                 for (size_t j = 0; j < cell_keys.size(); j++)
@@ -719,7 +722,6 @@ entity::populate_entity_data(entity_data_cpu& ed) const
                         ed.fluxdofs_mine.at(face_base + i) = cell_base + j;
         }
     }
-
 }
 
 

tests/test_basics.cpp

@@ -151,7 +151,7 @@ test_basics(int geometric_order, int approximation_order)
 
     return 0;
 }
-#include <fstream>
+
 int
 test_normals(int geometric_order, int approximation_order)
 {
@@ -165,6 +165,8 @@ test_normals(int geometric_order, int approximation_order)
     std::cout << cyanfg << "Testing normals (order " << geometric_order;
     std::cout << ")" << reset << std::endl;
 
+    int errors = 0;
+
     for (auto [dim, pgtag] : pgs)
     {
         std::vector<int> tags;
@@ -200,6 +202,7 @@ test_normals(int geometric_order, int approximation_order)
                     auto pf = e0.face(4*iT+iF);
                     auto rf = e0.face_refelem(pf);
                     auto qps = pf.integration_points();
+                    auto fbar = pf.barycenter();
 
                     auto ek = element_key(pf);
                     if ( !std::binary_search(bm.begin(), bm.end(), ek) )
@@ -210,7 +213,7 @@ test_normals(int geometric_order, int approximation_order)
                     for (size_t iQp = 0; iQp < num_qp; iQp++)
                     {
                         Eigen::Vector3d n = ed.normals.row((4*iT+iF)*num_qp + iQp);
-                        point_3d p = qps[iQp].point();
+                        point_3d p = (geometric_order == 1) ? fbar : qps[iQp].point();
                         point_3d v = p - point_3d(0.5, 0.5, 0.5);
                         Eigen::Vector3d nref;
                         nref(0) = v.x(); nref(1) = v.y(); nref(2) = v.z();
@@ -218,16 +221,53 @@ test_normals(int geometric_order, int approximation_order)
 
                         if ( n.cross(nref).norm() > 1e-2 )
                         {
-                            std::cout << "Wrong normal. Expected: " << nref.transpose() << ", ";
-                            std::cout << "got: " << n.transpose() << std::endl;
-                            return 1;
+                            std::cout << "Wrong normal. Expected: " << greenfg << nref.transpose();
+                            std::cout << reset << ", " << "got: " << redfg << n.transpose();
+                            std::cout << reset << " diff: " << (nref-n).transpose() << std::endl;
+                            errors++;
+                        }
                     }
                 }
             }
         }
     }
+
+    return errors;
+}
+
+int
+test_normals_2()
+{
+    make_geometry(0, 0.1);
+    model mod(1, 1);
+
+    vec3d totflux = vec3d::Zero();
+    double tf = 0.0;
+
+    for (const auto& e : mod)
+    {
+        entity_data_cpu ed;
+        e.populate_entity_data(ed);
+
+        for (size_t iT = 0; iT < e.num_cells(); iT++)
+        {
+            vec3d flux = vec3d::Zero();
+
+            for (size_t iF = 0; iF < 4; iF++)
+            {
+                auto face_num = 4*iT+iF;
+                auto& pf = e.face(face_num);
+                vec3d n = ed.normals.row(face_num);
+                flux += n * pf.measure();
+            }
+            tf += flux.norm();
+            totflux += flux;
+        }
     }
 
+    COMPARE_VALUES_ABSOLUTE("Volume flux", tf, 0, 1e-14);
+    COMPARE_VALUES_ABSOLUTE("Volume flux", totflux.norm(), 0, 1e-14);
+
     return 0;
 }
 
@@ -240,7 +280,7 @@ int main(void)
 
 
     int failed_tests = 0;
-
+/*
     std::cout << Bmagentafg << " *** TESTING: BASIC OPERATIONS ***" << reset << std::endl; 
     for (size_t go = 1; go < 5; go++)
     {
@@ -249,10 +289,13 @@ int main(void)
             failed_tests += test_basics(go, ao);
         }
     }
+*/
+    //std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION ***" << reset << std::endl;
+    //for (size_t i = 1; i < 5; i++)
+    //    test_normals(i,i);
 
-    std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION ***" << reset << std::endl;
-    for (size_t i = 1; i < 5; i++)
-        test_normals(i,i);
+    std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION (2) ***" << reset << std::endl;
+    test_normals_2();
 
     return failed_tests;
 }
\ No newline at end of file

tests/test_differentiation_gpu.cpp

@@ -41,8 +41,6 @@ make_geometry(int order, double mesh_h)
     gmm::setSize(vp, mesh_h);
 }
 
-#define WRITE_TEST_OUTPUTS
-
 int test_differentiation_convergence(int geometric_order, int approximation_order)
 {
     std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });

tests/test_lifting.cpp

@@ -70,15 +70,15 @@ int test_lifting(int geometric_order, int approximation_order)
 
     /* Test field */
     auto Fx = [](const point_3d& pt) -> double {
-        return 1; //std::sin(M_PI*pt.x()); 
+        return std::sin(M_PI*pt.x()); 
     };
 
     auto Fy = [](const point_3d& pt) -> double {
-        return 1; //std::sin(M_PI*pt.y());
+        return std::sin(M_PI*pt.y());
     };
 
     auto Fz = [](const point_3d& pt) -> double {
-        return 1; //std::sin(M_PI*pt.z());
+        return std::sin(M_PI*pt.z());
     };
 
     auto F = [&](const point_3d& pt) -> vecxd {
@@ -91,7 +91,7 @@ int test_lifting(int geometric_order, int approximation_order)
 
     /* Expected divergence */
     auto div_F = [](const point_3d& pt) -> double {
-        return 0;
+        //return 0;
         auto dFx_dx = M_PI*std::cos(M_PI*pt.x());
         auto dFy_dy = M_PI*std::cos(M_PI*pt.y());
         auto dFz_dz = M_PI*std::cos(M_PI*pt.z());
@@ -209,14 +209,20 @@ int test_lifting(int geometric_order, int approximation_order)
                 const auto& pqp = pqps[iQp];
                 const vecxd phi = re.basis_functions(iQp);
                 const matxd dphi = re.basis_gradients(iQp);
+                const mat3d J = pe.jacobian(iQp);
                 mass += pqp.weight() * phi * phi.transpose();
-                vol += pqp.weight() * dphi * F(pqp.point());
+                vol += pqp.weight() * (dphi*J.inverse()) * F(pqp.point());
             }
 
             vecxd expected = exp_field.segment(iT*num_bf, num_bf);
             vecxd lf = lift_field.segment(iT*num_bf, num_bf);
+            vecxd lfp = mass*lf;
             vecxd volp = mass.ldlt().solve(vol);
 
+            //std::cout << "***" << std::endl;
+            //std::cout << lfp.transpose() << std::endl;
+            //std::cout << vol.transpose() << std::endl;
+
             vecxd comp_field = lf - volp;
             vecxd diff = comp_field - expected;
             double loc_err = diff.dot(mass*diff);
@@ -258,3 +264,24 @@ int test_lifting(int geometric_order, int approximation_order)
 
     return 0;
 }
+
+int main(void)
+{
+    gmsh::initialize();
+    gmsh::option::setNumber("General.Terminal", 0);
+    gmsh::option::setNumber("Mesh.Algorithm", 1);
+    gmsh::option::setNumber("Mesh.Algorithm3D", 1);
+
+
+    int failed_tests = 0;
+
+    std::cout << Bmagentafg << " *** TESTING: LIFTING ***" << reset << std::endl;
+    for (size_t ao = 1; ao < 5; ao++)
+            test_lifting(1, ao);
+
+    gmsh::finalize();
+
+    return failed_tests;
+}
+
+