Forgot to add tests.

tests/test_basics.cpp

+#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

tests/test_differentiation.cpp

+#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

tests/test_differentiation_gpu.cpp

+#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

tests/test_lifting.cpp

+#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

tests/test_mass.cpp

+#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