diff --git a/tests/test_basics.cpp b/tests/test_basics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5eee4db9208ef620be4e340bb815f6d02b74bc19
--- /dev/null
+++ b/tests/test_basics.cpp
@@ -0,0 +1,258 @@
+#include <iostream>
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("test");
+
+ std::vector<std::pair<int,int>> objects;
+ objects.push_back(
+ std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0) )
+ );
+ objects.push_back(
+ std::make_pair(3, gmo::addSphere(0.5, 0.5, 0.5, SPHERE_RADIUS) )
+ );
+
+ std::vector<std::pair<int, int>> tools;
+ gmsh::vectorpair odt;
+ std::vector<gmsh::vectorpair> odtm;
+ gmo::fragment(objects, tools, odt, odtm);
+
+ gmo::synchronize();
+
+ int pgtag = gm::addPhysicalGroup(2, {7});
+ gm::setPhysicalName(2, pgtag, "SphereSurface");
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+int
+test_basics(int geometric_order, int approximation_order)
+{
+ make_geometry(0, 0.1);
+
+ model mod(geometric_order, approximation_order);
+
+ auto idx = geometric_order - 1;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ double vol_e0 = (4.0/3.0)*M_PI*SPHERE_RADIUS*SPHERE_RADIUS*SPHERE_RADIUS;
+ double vol_e1 = 1.0 - vol_e0;
+
+ /*****************************************/
+ std::cout << yellowfg << "Volumes by reference quadrature points and determinants";
+ std::cout << nofg << std::endl;
+ std::vector<double> t1_vals { 0.07, 0.00145, 0.00024, 1.83e-5 };
+ COMPARE_VALUES_RELATIVE("Entity 0 vol r+d", mod.entity_at(0).measure(), vol_e0, t1_vals[idx]);
+ std::vector<double> t2_vals { 0.0024, 5.02e-5, 8.21e-06, 6.35e-7 };
+ COMPARE_VALUES_RELATIVE("Entity 1 vol r+d", mod.entity_at(1).measure(), vol_e1, t2_vals[idx]);
+
+ /*****************************************/
+ std::cout << yellowfg << "Volumes by physical quadrature points" << nofg << std::endl;
+ double vol_e0_pqp = 0.0;
+ auto& e0 = mod.entity_at(0);
+ for (const auto& e : e0)
+ {
+ const auto& pqps = e.integration_points();
+ for (const auto& qp : pqps)
+ vol_e0_pqp += qp.weight();
+ }
+ COMPARE_VALUES_RELATIVE("Entity 0 vol pqp", vol_e0_pqp, vol_e0, t1_vals[idx]);
+
+ double vol_e1_pqp = 0.0;
+ auto& e1 = mod.entity_at(1);
+ for (const auto& e : e1)
+ {
+ const auto& pqps = e.integration_points();
+ for (const auto& qp : pqps)
+ vol_e1_pqp += qp.weight();
+ }
+ COMPARE_VALUES_RELATIVE("Entity 1 vol pqp", vol_e1_pqp, vol_e1, t2_vals[idx]);
+
+ /*****************************************/
+ std::cout << yellowfg << "Volumes by jacobians" << nofg << std::endl;
+ double vol_e0_j = 0.0;
+ for (const auto& e : e0)
+ {
+ const auto& js = e.jacobians();
+ auto& re = e0.cell_refelem(e);
+ auto rqps = re.integration_points();
+ auto dets = e.determinants();
+ for (size_t i = 0; i < js.size(); i++)
+ {
+ double dj = js[i].determinant();
+ double dg = dets[i];
+ auto relerr = dj/dg - 1;
+ assert( std::abs(relerr) < 3e-14 );
+ vol_e0_j += dj*rqps[i].weight();
+ }
+ }
+ COMPARE_VALUES_RELATIVE("Entity 0 vol j", vol_e0_j, vol_e0, t1_vals[idx]);
+
+ double vol_e1_j = 0.0;
+ for (const auto& e : e1)
+ {
+ const auto& js = e.jacobians();
+ auto& re = e1.cell_refelem(e);
+ auto rqps = re.integration_points();
+ auto dets = e.determinants();
+ for (size_t i = 0; i < js.size(); i++)
+ {
+ double dj = js[i].determinant();
+ double dg = dets[i];
+ auto relerr = dj/dg - 1;
+ assert( std::abs(relerr) < 3e-14 );
+ vol_e1_j += dj*rqps[i].weight();
+ }
+ }
+ COMPARE_VALUES_RELATIVE("Entity 1 vol j", vol_e1_j, vol_e1, t2_vals[idx]);
+
+ /*****************************************/
+ std::cout << yellowfg << "Volumes by integrating f(x) = 1 over the domain";
+ std::cout << nofg << std::endl;
+ double vol_e0_mass = 0.0;
+ auto f = [](const point_3d& pt) -> double { return 1; };
+ vecxd f_e0 = e0.project(f);
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ auto num_bf = re.num_basis_functions();
+ matxd mass = e0.mass_matrix(iT);
+ vecxd f_pe = f_e0.segment(num_bf*iT, num_bf);
+ vol_e0_mass += f_pe.transpose() * mass * f_pe;
+ }
+ COMPARE_VALUES_RELATIVE("Entity 0 vol mass", vol_e0_mass, vol_e0, t1_vals[idx]);
+
+ double vol_e1_mass = 0.0;
+ vecxd f_e1 = e1.project(f);
+ for (size_t iT = 0; iT < e1.num_cells(); iT++)
+ {
+ auto& pe = e1.cell(iT);
+ auto& re = e1.cell_refelem(pe);
+ auto num_bf = re.num_basis_functions();
+ matxd mass = e1.mass_matrix(iT);
+ vecxd f_pe = f_e1.segment(num_bf*iT, num_bf);
+ vol_e1_mass += f_pe.transpose() * mass * f_pe;
+ }
+ COMPARE_VALUES_RELATIVE("Entity 1 vol mass", vol_e1_mass, vol_e1, t2_vals[idx]);
+
+ return 0;
+}
+#include <fstream>
+int
+test_normals(int geometric_order, int approximation_order)
+{
+ make_geometry(0, 0.1);
+ model mod(geometric_order, approximation_order);
+
+ gmsh::vectorpair pgs;
+ gm::getPhysicalGroups(pgs);
+ assert(pgs.size() == 1);
+
+ std::cout << cyanfg << "Testing normals (order " << geometric_order;
+ std::cout << ")" << reset << std::endl;
+
+ for (auto [dim, pgtag] : pgs)
+ {
+ std::vector<int> tags;
+ gm::getEntitiesForPhysicalGroup(dim, pgtag, tags);
+ assert(tags.size() == 1);
+
+ std::vector<int> elemTypes;
+ gmm::getElementTypes(elemTypes, dim, tags[0]);
+
+ assert(elemTypes.size() == 1);
+ if (elemTypes.size() != 1)
+ throw std::invalid_argument("Only one element type per entity is allowed");
+
+ for (auto& elemType : elemTypes)
+ {
+ std::vector<size_t> faceTags, faceNodesTags;
+ gmm::getElementsByType(elemType, faceTags, faceNodesTags, tags[0]);
+
+ std::sort(faceTags.begin(), faceTags.end());
+
+ auto& e0 = mod.entity_at(0);
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+ auto ft = e0.face_tags();
+
+ auto bm = mod.get_bnd(7);
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ auto pf = e0.face(4*iT+iF);
+ auto rf = e0.face_refelem(pf);
+ auto qps = pf.integration_points();
+
+ auto ek = element_key(pf);
+ if ( !std::binary_search(bm.begin(), bm.end(), ek) )
+ continue;
+
+ size_t num_qp = (geometric_order == 1) ? 1 : qps.size();
+
+ 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 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();
+ nref /= nref.norm();
+
+ if ( n.cross(nref).norm() > 1e-2 )
+ {
+ std::cout << "Wrong normal. Expected: " << nref.transpose() << ", ";
+ std::cout << "got: " << n.transpose() << std::endl;
+ return 1;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ 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: BASIC OPERATIONS ***" << reset << std::endl;
+ for (size_t go = 1; go < 5; go++)
+ {
+ for (size_t ao = go; ao < 5; ao++)
+ {
+ 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);
+
+ return failed_tests;
+}
\ No newline at end of file
diff --git a/tests/test_differentiation.cpp b/tests/test_differentiation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ffc813b127b0feef2576bdcec5d26925f5e3010
--- /dev/null
+++ b/tests/test_differentiation.cpp
@@ -0,0 +1,152 @@
+#include <iostream>
+#include <iomanip>
+
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+#include "entity_data.h"
+#include "kernels_cpu.h"
+#include "timecounter.h"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("difftest");
+ gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
+ gmo::synchronize();
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+int test_differentiation_convergence(int geometric_order, int approximation_order)
+{
+ std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
+ std::vector<double> errors;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ auto f = [](const point_3d& pt) -> double {
+ return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dx = [](const point_3d& pt) -> double {
+ return M_PI*std::cos(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dy = [](const point_3d& pt) -> double {
+ return M_PI*std::sin(M_PI*pt.x())*std::cos(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dz = [](const point_3d& pt) -> double {
+ return M_PI*std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::cos(M_PI*pt.z());
+ };
+
+ std::vector<double> errs_x;
+ std::vector<double> errs_y;
+ std::vector<double> errs_z;
+
+ for (size_t i = 0; i < sizes.size(); i++)
+ {
+ double h = sizes[i];
+ make_geometry(0,h);
+ model mod(geometric_order, approximation_order);
+
+ auto& e0 = mod.entity_at(0);
+ e0.sort_by_orientation();
+ vecxd Pf_e0 = e0.project(f);
+ vecxd df_dx_e0 = vecxd::Zero( Pf_e0.rows() );
+ vecxd df_dy_e0 = vecxd::Zero( Pf_e0.rows() );
+ vecxd df_dz_e0 = vecxd::Zero( Pf_e0.rows() );
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ timecounter tc;
+ tc.tic();
+ compute_field_derivatives(ed, Pf_e0, df_dx_e0, df_dy_e0, df_dz_e0);
+ double time = tc.toc();
+
+ if (geometric_order == 1)
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = 21*ed.num_bf*ed.num_bf*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ else
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = ((21*ed.num_bf+6)*ed.num_bf*ed.num_qp + 3*(2*ed.num_bf-1)*ed.num_bf)*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+
+ vecxd Pdf_dx_e0 = e0.project(df_dx);
+ vecxd Pdf_dy_e0 = e0.project(df_dy);
+ vecxd Pdf_dz_e0 = e0.project(df_dz);
+
+ double err_x = 0.0;
+ double err_y = 0.0;
+ double err_z = 0.0;
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ matxd mass = e0.mass_matrix(iT);
+ auto num_bf = re.num_basis_functions();
+ vecxd diff_x = Pdf_dx_e0.segment(iT*num_bf, num_bf) - df_dx_e0.segment(iT*num_bf, num_bf);
+ vecxd diff_y = Pdf_dy_e0.segment(iT*num_bf, num_bf) - df_dy_e0.segment(iT*num_bf, num_bf);
+ vecxd diff_z = Pdf_dz_e0.segment(iT*num_bf, num_bf) - df_dz_e0.segment(iT*num_bf, num_bf);
+
+ err_x += diff_x.dot(mass*diff_x);
+ err_y += diff_y.dot(mass*diff_y);
+ err_z += diff_z.dot(mass*diff_z);
+ }
+ std::cout << "Errors: " << std::sqrt(err_x) << " " << std::sqrt(err_y);
+ std::cout << " " << std::sqrt(err_z) << std::endl;
+
+ errs_x.push_back( std::sqrt(err_x) );
+ errs_y.push_back( std::sqrt(err_y) );
+ errs_z.push_back( std::sqrt(err_z) );
+ }
+
+ double rate_x = 0.0;
+ double rate_y = 0.0;
+ double rate_z = 0.0;
+
+ std::cout << Byellowfg << "rate df/dx rate df/dy rate df/dz" << reset << std::endl;
+ for (size_t i = 1; i < sizes.size(); i++)
+ {
+ std::cout << (rate_x = std::log2(errs_x[i-1]/errs_x[i]) ) << " ";
+ std::cout << (rate_y = std::log2(errs_y[i-1]/errs_y[i]) ) << " ";
+ std::cout << (rate_z = std::log2(errs_z[i-1]/errs_z[i]) ) << std::endl;
+ }
+
+ COMPARE_VALUES_ABSOLUTE("df/dx", rate_x, double(approximation_order), 0.2);
+ COMPARE_VALUES_ABSOLUTE("df/dy", rate_y, double(approximation_order), 0.2);
+ COMPARE_VALUES_ABSOLUTE("df/dz", rate_z, double(approximation_order), 0.2);
+
+ 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: DIFFERENTIATION ***" << reset << std::endl;
+ for (size_t go = 1; go < 5; go++)
+ for (size_t ao = go; ao < 5; ao++)
+ failed_tests += test_differentiation_convergence(go, ao);
+
+ return failed_tests;
+}
\ No newline at end of file
diff --git a/tests/test_differentiation_gpu.cpp b/tests/test_differentiation_gpu.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c4e65650307234c019ec9412d4822b50fa57aa03
--- /dev/null
+++ b/tests/test_differentiation_gpu.cpp
@@ -0,0 +1,158 @@
+#include <iostream>
+#include <iomanip>
+
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+#include "entity_data.h"
+#include "kernels_cpu.h"
+#include "timecounter.h"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("difftest");
+ gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
+ gmo::synchronize();
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+int test_differentiation_convergence(int geometric_order, int approximation_order)
+{
+ std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
+ std::vector<double> errors;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ auto f = [](const point_3d& pt) -> double {
+ return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dx = [](const point_3d& pt) -> double {
+ return M_PI*std::cos(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dy = [](const point_3d& pt) -> double {
+ return M_PI*std::sin(M_PI*pt.x())*std::cos(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+ auto df_dz = [](const point_3d& pt) -> double {
+ return M_PI*std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::cos(M_PI*pt.z());
+ };
+
+ std::vector<double> errs_x;
+ std::vector<double> errs_y;
+ std::vector<double> errs_z;
+
+ for (size_t i = 0; i < sizes.size(); i++)
+ {
+ double h = sizes[i];
+ make_geometry(0,h);
+ model mod(geometric_order, approximation_order);
+
+ auto& e0 = mod.entity_at(0);
+ e0.sort_by_orientation();
+ vecxd Pf_e0 = e0.project(f);
+ vecxd df_dx_e0 = vecxd::Zero( Pf_e0.rows() );
+ vecxd df_dy_e0 = vecxd::Zero( Pf_e0.rows() );
+ vecxd df_dz_e0 = vecxd::Zero( Pf_e0.rows() );
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ entity_data_gpu edg(ed);
+
+#if 0
+ timecounter tc;
+ tc.tic();
+ compute_field_derivatives(ed, Pf_e0, df_dx_e0, df_dy_e0, df_dz_e0);
+ double time = tc.toc();
+
+ if (geometric_order == 1)
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = 21*ed.num_bf*ed.num_bf*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ else
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = ((21*ed.num_bf+6)*ed.num_bf*ed.num_qp + 3*(2*ed.num_bf-1)*ed.num_bf)*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+
+ vecxd Pdf_dx_e0 = e0.project(df_dx);
+ vecxd Pdf_dy_e0 = e0.project(df_dy);
+ vecxd Pdf_dz_e0 = e0.project(df_dz);
+
+ double err_x = 0.0;
+ double err_y = 0.0;
+ double err_z = 0.0;
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ matxd mass = e0.mass_matrix(iT);
+ auto num_bf = re.num_basis_functions();
+ vecxd diff_x = Pdf_dx_e0.segment(iT*num_bf, num_bf) - df_dx_e0.segment(iT*num_bf, num_bf);
+ vecxd diff_y = Pdf_dy_e0.segment(iT*num_bf, num_bf) - df_dy_e0.segment(iT*num_bf, num_bf);
+ vecxd diff_z = Pdf_dz_e0.segment(iT*num_bf, num_bf) - df_dz_e0.segment(iT*num_bf, num_bf);
+
+ err_x += diff_x.dot(mass*diff_x);
+ err_y += diff_y.dot(mass*diff_y);
+ err_z += diff_z.dot(mass*diff_z);
+ }
+ std::cout << "Errors: " << std::sqrt(err_x) << " " << std::sqrt(err_y);
+ std::cout << " " << std::sqrt(err_z) << std::endl;
+
+ errs_x.push_back( std::sqrt(err_x) );
+ errs_y.push_back( std::sqrt(err_y) );
+ errs_z.push_back( std::sqrt(err_z) );
+#endif
+ }
+
+ double rate_x = 0.0;
+ double rate_y = 0.0;
+ double rate_z = 0.0;
+
+#if 0
+ std::cout << Byellowfg << "rate df/dx rate df/dy rate df/dz" << reset << std::endl;
+ for (size_t i = 1; i < sizes.size(); i++)
+ {
+ std::cout << (rate_x = std::log2(errs_x[i-1]/errs_x[i]) ) << " ";
+ std::cout << (rate_y = std::log2(errs_y[i-1]/errs_y[i]) ) << " ";
+ std::cout << (rate_z = std::log2(errs_z[i-1]/errs_z[i]) ) << std::endl;
+ }
+
+ COMPARE_VALUES_ABSOLUTE("df/dx", rate_x, double(approximation_order), 0.2);
+ COMPARE_VALUES_ABSOLUTE("df/dy", rate_y, double(approximation_order), 0.2);
+ COMPARE_VALUES_ABSOLUTE("df/dz", rate_z, double(approximation_order), 0.2);
+#endif
+
+ 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: DIFFERENTIATION ***" << reset << std::endl;
+ for (size_t go = 1; go < 5; go++)
+ for (size_t ao = go; ao < 5; ao++)
+ failed_tests += test_differentiation_convergence(go, ao);
+
+ return failed_tests;
+}
\ No newline at end of file
diff --git a/tests/test_lifting.cpp b/tests/test_lifting.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..154cccee435232ee017419bcf649732786d424f6
--- /dev/null
+++ b/tests/test_lifting.cpp
@@ -0,0 +1,260 @@
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+#include "entity_data.h"
+#include "kernels_cpu.h"
+#include "timecounter.h"
+#include "silo_output.hpp"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("difftest");
+ gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
+ gmo::synchronize();
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+
+int test_lifting_xx(int geometric_order, int approximation_order)
+{
+ make_geometry(0, 0.25);
+ model mod(geometric_order, approximation_order);
+
+ auto& e0 = mod.entity_at(0);
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ for (size_t iO = 0; iO < e0.num_cell_orientations(); iO++)
+ {
+ auto& re = e0.cell_refelem(iO);
+ auto num_bf = re.num_basis_functions();
+
+ auto lift_cols = ed.lifting_matrices.cols();
+
+ matxd mass = re.mass_matrix();
+ matxd curlift = mass * ed.lifting_matrices.block(iO*num_bf, 0, num_bf, lift_cols);
+
+ vecxd ones_cell = vecxd::Ones(num_bf);
+ vecxd ones_face = vecxd::Ones(lift_cols/4);
+
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ matxd liftblock = curlift.block(0,iF*lift_cols/4,num_bf,lift_cols/4);
+ std::cout << ones_cell.dot(liftblock*ones_face) << std::endl;
+ }
+ }
+
+ return 0;
+}
+
+int test_lifting(int geometric_order, int approximation_order)
+{
+ std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
+ std::vector<double> errors;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ /* Test field */
+ auto Fx = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.x());
+ };
+
+ auto Fy = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.y());
+ };
+
+ auto Fz = [](const point_3d& pt) -> double {
+ return 1; //std::sin(M_PI*pt.z());
+ };
+
+ auto F = [&](const point_3d& pt) -> vecxd {
+ vec3d Fret;
+ Fret(0) = Fx(pt);
+ Fret(1) = Fy(pt);
+ Fret(2) = Fz(pt);
+ return Fret;
+ };
+
+ /* Expected divergence */
+ auto div_F = [](const point_3d& pt) -> double {
+ 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());
+ return dFx_dx + dFy_dy + dFz_dz;
+ };
+
+ for (size_t sz = 0; sz < sizes.size(); sz++)
+ {
+ double h = sizes[sz];
+ make_geometry(0,h);
+ model mod(geometric_order, approximation_order);
+
+ std::stringstream ss;
+ ss << "lifting_test_" << sz << ".silo";
+
+ silo s;
+ s.create_db( ss.str() );
+ s.import_mesh_from_gmsh();
+ s.write_mesh();
+
+ auto& e0 = mod.entity_at(0);
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ vecxd exp_field = e0.project(div_F);
+ vecxd Fx_proj = e0.project(Fx);
+ vecxd Fy_proj = e0.project(Fy);
+ vecxd Fz_proj = e0.project(Fz);
+
+ auto flx_size = e0.num_cells() * 4*ed.num_fluxes;
+ vecxd Fx_flux = vecxd::Zero( flx_size );
+ vecxd Fy_flux = vecxd::Zero( flx_size );
+ vecxd Fz_flux = vecxd::Zero( flx_size );
+
+ for (size_t i = 0; i < flx_size; i++)
+ {
+ Fx_flux(i) = Fx_proj( ed.fluxdofs_mine[i] );
+ Fy_flux(i) = Fy_proj( ed.fluxdofs_mine[i] );
+ Fz_flux(i) = Fz_proj( ed.fluxdofs_mine[i] );
+ }
+
+ vecxd flux = vecxd::Zero( flx_size );
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ vec3d n = ed.normals.row(4*iT+iF);
+ for (size_t k = 0; k < ed.num_fluxes; k++)
+ {
+ auto dof_ofs = (4*iT+iF)*ed.num_fluxes + k;
+ flux(dof_ofs) = Fx_flux(dof_ofs)*n(0) +
+ Fy_flux(dof_ofs)*n(1) +
+ Fz_flux(dof_ofs)*n(2);
+ }
+ }
+ }
+
+ vecxd lift_field = vecxd::Zero( e0.num_cells() * ed.num_bf );
+#if 0
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ auto& pf = e0.face(4*iT+iF);
+ auto& rf = e0.face_refelem(pf);
+ auto num_fluxes = rf.num_basis_functions();
+ const auto pqps = pf.integration_points();
+
+ matxd mass = matxd::Zero(num_fluxes, num_fluxes);
+ vecxd rhs = vecxd::Zero(num_fluxes);
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ vec3d n;
+ if (ed.g_order == 1)
+ n = ed.normals.row(4*iT+iF);
+ else
+ n = ed.normals.row( (4*iT+iF)*pqps.size() + iQp );
+
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = rf.basis_functions(iQp);
+ mass += pqp.weight() * phi * phi.transpose();
+ rhs += pqp.weight() * F(pqp.point()).dot(n) * phi;
+ }
+ auto fluxofs = (4*iT+iF)*num_fluxes;
+ flux.segment(fluxofs, num_fluxes) = mass.ldlt().solve(rhs);
+ }
+ }
+#endif
+ timecounter tc;
+ tc.tic();
+ compute_flux_lifting(ed, flux, lift_field);
+ double time = tc.toc();
+
+ std::vector<double> var_div_F( e0.num_cells() );
+ std::vector<double> var_lift( e0.num_cells() );
+ std::vector<double> var_volp( e0.num_cells() );
+ std::vector<double> var_elift( e0.num_cells() );
+ std::vector<double> var_err( e0.num_cells() );
+ std::vector<double> var_orient( e0.num_cells() );
+
+ double err = 0.0;
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ auto num_bf = re.num_basis_functions();
+ matxd mass = matxd::Zero(num_bf, num_bf);
+ vecxd vol = vecxd::Zero(num_bf);
+
+ const auto pqps = pe.integration_points();
+ for(size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = re.basis_functions(iQp);
+ const matxd dphi = re.basis_gradients(iQp);
+ mass += pqp.weight() * phi * phi.transpose();
+ vol += pqp.weight() * dphi * F(pqp.point());
+ }
+
+ vecxd expected = exp_field.segment(iT*num_bf, num_bf);
+ vecxd lf = lift_field.segment(iT*num_bf, num_bf);
+ vecxd volp = mass.ldlt().solve(vol);
+
+ vecxd comp_field = lf - volp;
+ vecxd diff = comp_field - expected;
+ double loc_err = diff.dot(mass*diff);
+ err += loc_err;
+
+ auto bar = pe.barycenter();
+ vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
+
+ var_div_F.at( pe.original_position() ) = expected.dot(phi_bar);
+ var_lift.at( pe.original_position() ) = lf.dot(phi_bar);
+ var_volp.at( pe.original_position() ) = volp.dot(phi_bar);
+ var_elift.at( pe.original_position() ) = (volp + expected).dot(phi_bar);
+ var_err.at( pe.original_position() ) = loc_err;
+ var_orient.at( pe.original_position() ) = pe.orientation();
+ }
+
+ s.write_zonal_variable("div_F", var_div_F);
+ s.write_zonal_variable("lift", var_lift);
+ s.write_zonal_variable("vol", var_volp);
+ s.write_zonal_variable("exp_lift", var_elift);
+ s.write_zonal_variable("err", var_err);
+ s.write_zonal_variable("orient", var_orient);
+
+ std::cout << std::sqrt(err) << std::endl;
+
+ if (geometric_order == 1)
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = 3*(ed.num_bf)*4*ed.num_fluxes*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ else
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = ((21*ed.num_bf+6)*ed.num_bf*ed.num_qp + 3*(2*ed.num_bf-1)*ed.num_bf)*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ }
+
+ return 0;
+}
\ No newline at end of file
diff --git a/tests/test_mass.cpp b/tests/test_mass.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5690cefb3844db6b61bfa757d00f465210b5f60d
--- /dev/null
+++ b/tests/test_mass.cpp
@@ -0,0 +1,115 @@
+#include <iostream>
+#include <iomanip>
+
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("difftest");
+ gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
+ gmo::synchronize();
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+#define OVERINTEGRATE
+int test_mass_convergence(int geometric_order, int approximation_order)
+{
+ std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
+ std::vector<double> errors;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ auto test_f = [](const point_3d& pt) -> double {
+ return std::sin(M_PI*pt.x())*std::sin(M_PI*pt.y())*std::sin(M_PI*pt.z());
+ };
+
+ for (size_t i = 0; i < sizes.size(); i++)
+ {
+ double h = sizes[i];
+ make_geometry(0,h);
+ model mod(geometric_order, approximation_order);
+
+ auto& e0 = mod.entity_at(0);
+ e0.sort_by_orientation();
+ vecxd test_f_e0 = e0.project(test_f);
+ double proj_error_e0 = 0.0;
+#ifdef OVERINTEGRATE
+ auto [rps, pqps] = integrate(e0, 2*approximation_order+1);
+ auto num_qp = rps.size();
+#endif
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ auto num_bf = re.num_basis_functions();
+
+ vecxd pvals = test_f_e0.segment(iT*num_bf, num_bf);
+#ifdef OVERINTEGRATE
+ for (size_t iQp = 0; iQp < num_qp; iQp++)
+ {
+ vecxd phi = re.basis_functions( rps[iQp] );
+ auto pqpofs = num_qp*iT + iQp;
+ double diff = (pvals.dot(phi) - test_f(pqps[pqpofs].point()));
+ proj_error_e0 += std::abs(pqps[pqpofs].weight())*diff*diff;
+ }
+#else
+ auto pqps = pe.integration_points();
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ vecxd phi = re.basis_functions( iQp );
+ double diff = (pvals.dot(phi) - test_f(pqps[iQp].point()));
+ proj_error_e0 += std::abs(pqps[iQp].weight())*diff*diff;
+ }
+#endif
+ }
+
+ errors.push_back( std::sqrt(proj_error_e0) );
+ }
+
+ std::cout << Byellowfg << "Errors: " << reset;
+ for (auto& error : errors)
+ std::cout << error << " ";
+ std::cout << std::endl;
+
+ std::cout << Byellowfg << "Rates: " << reset;
+ double rate = 0.0;
+ for (size_t i = 1; i < errors.size(); i++)
+ std::cout << (rate = std::log2(errors[i-1]/errors[i])) << " ";
+ std::cout << std::endl;
+
+ COMPARE_VALUES_ABSOLUTE("Projection", rate, double(approximation_order+1), 0.2);
+
+ 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: PROJECTIONS ***" << reset << std::endl;
+ for (size_t go = 1; go < 5; go++)
+ {
+ for (size_t ao = go; ao < 5; ao++)
+ failed_tests += test_mass_convergence(go, ao);
+ }
+
+ return failed_tests;
+}
\ No newline at end of file