diff --git a/CMakeLists.txt b/CMakeLists.txt
index 81f61794abc979ec3abd1f5a3aba88c3679750c4..909b1e5e4d4e274d0a631b1f9bb3421123709bc4 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/src/entity.cpp b/src/entity.cpp
index 71f63439a6ffbfb5fa0fee9537211e6aea05d037..d6801d4006661d498773a8fcf4a763611f6eb895 100644
--- a/src/entity.cpp
+++ b/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;
}
}
-
}
diff --git a/tests/test_basics.cpp b/tests/test_basics.cpp
index 5eee4db9208ef620be4e340bb815f6d02b74bc19..1235eb6d2c6630a6a06e90447f9160d594bf2da4 100644
--- a/tests/test_basics.cpp
+++ b/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,9 +221,10 @@ 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++;
}
}
}
@@ -228,6 +232,42 @@ test_normals(int geometric_order, int approximation_order)
}
}
+ 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
diff --git a/tests/test_differentiation_gpu.cpp b/tests/test_differentiation_gpu.cpp
index 5b368d136c42f106f0c0a8b63029781131fa77ac..f074f0d7e51e4fa34bbe86e7c4505de995ec1c82 100644
--- a/tests/test_differentiation_gpu.cpp
+++ b/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 });
diff --git a/tests/test_lifting.cpp b/tests/test_lifting.cpp
index 154cccee435232ee017419bcf649732786d424f6..d6bd47cc96eb9a68813e049d8ea36d5789c9996b 100644
--- a/tests/test_lifting.cpp
+++ b/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);
@@ -257,4 +263,25 @@ int test_lifting(int geometric_order, int approximation_order)
}
return 0;
-}
\ No newline at end of file
+}
+
+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;
+}
+
+