Tests, differentiation kernel.

CMakeLists.txt

@@ -3,7 +3,7 @@ project(gmsh_gpu_dg)
 include(CheckLanguage)
 
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CUDA_ARCHITECTURES 35)
 set(CMAKE_CUDA_HOST_COMPILER "/usr/bin/g++-6")
@@ -135,7 +135,7 @@ endif ()
 ######################################################################
 ## Compiler optimization options
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Weverything -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-zero-as-null-pointer-constant")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Weverything -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-zero-as-null-pointer-constant -Wno-global-constructors -Wno-padded -Wno-shorten-64-to-32 -Wno-sign-conversion -Wno-old-style-cast -Wno-sign-compare -Wno-c99-extensions -Wno-extra-semi-stmt")
 
 set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g -DNDEBUG")
@@ -230,5 +230,11 @@ endif()
 include_directories(contrib/sgr)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 
-add_executable(gmsh_gpu_dg main.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp)
+add_executable(gmsh_gpu_dg main.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
 target_link_libraries(gmsh_gpu_dg ${LINK_LIBS})
+
+add_executable(test test.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
+target_link_libraries(test ${LINK_LIBS})
+
+add_executable(dump_orient dump_orient.cpp)
+target_link_libraries(dump_orient ${LINK_LIBS})

dump_orient.cpp

+#include <iostream>
+#include <cassert>
+
+#include "gmsh.h"
+
+namespace g     = gmsh;
+namespace go    = gmsh::option;
+namespace gm    = gmsh::model;
+namespace gmm   = gmsh::model::mesh;
+namespace gmo   = gmsh::model::occ;
+namespace gmg   = gmsh::model::geo;
+namespace gv    = gmsh::view;
+
+using gvp_t = gmsh::vectorpair;
+
+#define SPHERE_RADIUS 0.2
+#define DIMENSION(x) x
+
+static void
+make_geometry(int order)
+{
+    gm::add("difftest");
+
+    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();
+
+    gvp_t vp;
+    gm::getEntities(vp);
+    gmm::setSize(vp, 0.5);
+
+    gmm::generate( DIMENSION(3) );
+    gmm::setOrder( order );
+}
+
+int main(void)
+{
+    gmsh::initialize();
+
+    int gorder = 2;
+
+    make_geometry(gorder);
+
+    gvp_t gmsh_entities;
+    gm::getEntities(gmsh_entities, DIMENSION(3));
+
+    for (auto [dim, tag] : gmsh_entities)
+    {
+        std::vector<int> elemTypes;
+        gmm::getElementTypes(elemTypes, dim, tag);
+        
+        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> elemTags;
+            std::vector<size_t> elemNodesTags;
+            gmm::getElementsByType(elemType, elemTags, elemNodesTags, tag);
+
+            std::vector<int> orientations;
+            gmm::getBasisFunctionsOrientationForElements(elemType, "H1Legendre3",
+                orientations, tag);
+
+            std::vector<std::size_t> face_nodes;
+            gmm::getElementFaceNodes(elemType, 3, face_nodes, tag); // 3 is for triangles
+            int fe_tag = gm::addDiscreteEntity(dim-1);
+            int eleType2D = gmsh::model::mesh::getElementType("triangle", gorder);
+            gmsh::model::mesh::addElementsByType(fe_tag, eleType2D, {}, face_nodes);
+
+            std::vector<int> orientations2D;
+            gmm::getBasisFunctionsOrientationForElements(eleType2D, "H1Legendre3",
+                orientations2D, fe_tag);
+
+            std::cout << "3D: " << orientations.size() << std::endl;
+            std::cout << "2D: " << orientations2D.size() << std::endl;
+            assert(orientations2D.size() == 4*orientations.size());
+
+            using op = std::pair<int, std::array<int, 4>>;
+            std::vector<op> ordermap;
+
+            for (size_t i = 0; i < orientations.size(); i++)
+            {
+                op p;
+                p.first = orientations[i];
+                std::array<int, 4> or2d;
+
+                for (size_t j = 0; j < 4; j++)
+                    or2d[j] = orientations2D[4*i+j];
+                
+                p.second = or2d;
+                ordermap.push_back(p);
+            }
+
+            auto comp = [](const op& a, const op& b) -> bool { return a.first < b.first; };
+
+            std::sort(ordermap.begin(), ordermap.end(), comp);
+
+            for (auto& om : ordermap)
+            {
+                std::cout << om.first << ": ";
+                for (auto& o : om.second)
+                    std::cout << o << " ";
+                std::cout << std::endl;
+            }
+        }
+    }
+
+    gmsh::finalize();
+    return 0;
+}

entity.cpp

 #include <iostream>
-
+#include <cassert>
 #include "sgr.hpp"
 
 #include "gmsh_io.h"
 #include "entity.h"
 
-entity::entity(int pdim, int ptag, int pelemType, int porder)
-    : g_dim(pdim), g_tag(ptag), g_elemType(pelemType), order(porder)
+
+entity::entity(const entity_params& ep)
+    : dim(ep.dim), tag(ep.tag), elemType(ep.etype), g_order(ep.gorder),
+      a_order(ep.aorder), cur_elem_ordering(entity_ordering::GMSH)
 {
-    /* Get all the elements belonging to this group */
-    gmm::getElementsByType(g_elemType, g_elemTags, g_elemNodesTags, g_tag);
-    auto num_elems = g_elemTags.size();
+    /* Prepare reference elements */
+    reference_elements_factory ref(ep);
+    reference_cells = ref.get_elements();
 
-#ifdef ENABLE_DEBUG_PRINT
-    std::cout << yellowfg << "entity::entity(): elemType " << g_elemType
-              << ", elems: " << num_elems << std::endl;
-    std::cout << quadrature_name(2*order) << " " << basis_func_name(order) << " "
-              << basis_grad_name(order) << reset << std::endl;
-#endif
-    /* Get integration points and weights on refelem */
-    gmm::getIntegrationPoints(g_elemType, quadrature_name(2*order), ips, iws);
-    assert(ips.size() == 3*iws.size());
+    physical_elements_factory pef(ep);
+    physical_cells = pef.get_elements();
+
+    std::vector<std::size_t> face_nodes;
+    gmm::getElementFaceNodes(elemType, 3, face_nodes, tag); // 3 is for triangles
+    int tag_2D = gm::addDiscreteEntity(dim-1);
+    elemType_2D = gmsh::model::mesh::getElementType("triangle", g_order);
+    gmsh::model::mesh::addElementsByType(tag_2D, elemType_2D, {}, face_nodes);
+
+    entity_params ep_2D;
+    ep_2D.dim       = 2;
+    ep_2D.tag       = tag_2D;
+    ep_2D.etype     = elemType_2D;
+    ep_2D.gorder    = g_order; 
+    ep_2D.aorder    = a_order; 
+
+    reference_elements_factory ref_2D(ep_2D);
+    reference_faces = ref_2D.get_elements();
+
+    physical_elements_factory pef_2D(ep_2D);
+    physical_faces = pef_2D.get_elements();
+}
+
+double
+entity::measure(void) const
+{
+    auto meas = 0.0;
+    for (auto& pe : physical_cells)
+        meas += pe.measure();
+    
+    return meas;
+}
+
+vecxd
+entity::project(const scalar_function& function) const
+{
+    auto num_bf = reference_cells[0].num_basis_functions();
+
+    vecxd ret = vecxd::Zero( physical_cells.size() * num_bf );
+    for (size_t iT = 0; iT < physical_cells.size(); iT++)
+    {
+        const auto& pe = physical_cells[iT];
+        //assert( pe.orientation() == 0 );
+        assert( pe.orientation() < reference_cells.size() );
+        const auto& re = reference_cells[pe.orientation()];
+        const auto pqps = pe.integration_points();
+        const auto dets = pe.determinants();
+
+        assert(pqps.size() == dets.size());
+
+        matxd mass = matxd::Zero(num_bf, num_bf);
+        vecxd rhs = vecxd::Zero(num_bf);
+        for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+        {
+            const auto& pqp = pqps[iQp];
+            const vecxd phi = re.basis_functions(iQp);
+            mass += pqp.weight() * phi * phi.transpose();
+            rhs += pqp.weight() * function(pqp.point()) * phi; 
+        }
+
+        ret.segment(iT*num_bf, num_bf) = mass.ldlt().solve(rhs);
+    }
+
+    return ret;
+}
+
+const reference_element&
+entity::refelem(const physical_element& pe) const
+{
+    assert(pe.orientation() < reference_cells.size());
+    return reference_cells[pe.orientation()];
+}
+
+size_t
+entity::num_orientations(void) const
+{
+    return reference_cells.size();
+}
+
+size_t
+entity::num_elements(void) const
+{
+    return physical_cells.size();
+}
+
+std::vector<size_t>
+entity::num_elements_per_orientation(void) const
+{
+    std::vector<size_t> ret( num_orientations(), 0 );
+    for (auto& pe : physical_cells)
+    {
+        assert( pe.orientation() < ret.size() );
+        ret[ pe.orientation() ]++;
+    }
+
+    return ret;
+}
+
+const physical_element&
+entity::elem(size_t iT)
+{
+    assert(iT < physical_cells.size());
+    return physical_cells[iT];
+}
 
-    gmm::getJacobians(g_elemType, ips, jacs, dets, pts, g_tag);
+matxd
+entity::mass_matrix(size_t iT)
+{
+    const auto& pe = physical_cells[iT];
+    assert( pe.orientation() < reference_cells.size() );
+    const auto& re = reference_cells[pe.orientation()];
+    const auto rqps = re.integration_points();
+    const auto dets = pe.determinants();
+    auto num_bf = re.num_basis_functions();
 
-    assert(jacs.size() == 9*dets.size());
-    assert(pts.size() == 3*dets.size());
-    assert(dets.size() == num_elems*iws.size());
+        matxd mass = matxd::Zero(num_bf, num_bf);
 
-    int bf_nc;
-    gmm::getBasisFunctions(g_elemType, ips, basis_func_name(order), bf_nc, basisFunctions, g_num_bfun_orient);
-    assert(bf_nc == 1);
-    gmm::getBasisFunctionsOrientationForElements(g_elemType, basis_func_name(order), basisFunctionsOrientation);
+        for (size_t iQp = 0; iQp < rqps.size(); iQp++)
+        {
+            const auto& qp = rqps[iQp];
+            const vecxd phi = re.basis_functions(iQp);
+            mass += dets[iQp] * qp.weight() * phi * phi.transpose();
+        }
+    
+    return mass;
+}
 
-    gmm::getBasisFunctions(g_elemType, ips, basis_grad_name(order), bf_nc, basisGradients, g_num_bgrad_orient);
-    assert(bf_nc == 3);
-    gmm::getBasisFunctionsOrientationForElements(g_elemType, basis_grad_name(order), basisGradientsOrientation);
+void
+entity::sort_by_orientation(void)
+{
+    /* Sort by orientation and keep GMSH ordering within the same orientation */
+    auto comp = [](const physical_element& e1, const physical_element& e2) -> bool {
+        auto oe1 = e1.orientation();
+        auto oe2 = e2.orientation();
+        auto op1 = e1.original_position();
+        auto op2 = e2.original_position();
 
-    assert(g_num_bfun_orient == g_num_bgrad_orient);
+        if ( oe1 < oe2 and op1 < op2 )
+            return true;
 
-    gmm::getBarycenters(g_elemType, g_tag, false, false, barycenters);
-    assert(barycenters.size() == 3*num_elems);
+        return false;
+    };
 
-#ifdef ENABLE_DEBUG_PRINT
-    std::cout << yellowfg << "Orientations: " << g_num_bfun_orient << " " << g_num_bgrad_orient << reset << std::endl;
+    std::sort(physical_cells.begin(), physical_cells.end(), comp);
+    cur_elem_ordering = entity_ordering::BY_ORIENTATION;
+}
+
+void
+entity::sort_by_gmsh(void)
+{
+    /* Sort by GMSH original ordering */
+    auto comp = [](const physical_element& e1, const physical_element& e2) -> bool {
+        auto op1 = e1.original_position();
+        auto op2 = e2.original_position();
+
+        if ( op1 < op2 )
+            return true;
+
+        return false;
+    };
+
+    std::sort(physical_cells.begin(), physical_cells.end(), comp);
+    cur_elem_ordering = entity_ordering::GMSH;
+}
+
+int
+entity::gmsh_tag(void) const
+{
+    return tag;
+}
+
+int
+entity::gmsh_elem_type(void) const
+{
+    return elemType;
+}
+
+entity_ordering
+entity::current_elem_ordering(void) const
+{
+    return cur_elem_ordering;
+}
+
+
+void
+entity::populate_differentiation_matrices(entity_data_cpu<1>& ed) const
+{
+    /* In the case of geometric order = 1, the inverse of the
+     * mass matrix is embedded in the differentiation matrix */
+
+    auto dm_rows = ed.num_orientations * ed.num_bf;
+    auto dm_cols = 3 * ed.num_bf;
+    ed.differentiation_matrix = matxd::Zero(dm_rows, dm_cols);
+
+    for (size_t iO = 0; iO < ed.num_orientations; iO++)
+    {
+        matxd dm = matxd::Zero(ed.num_bf, 3*ed.num_bf);
+
+        auto& re = reference_cells[iO];
+        auto qps = re.integration_points();
+        auto num_qp = qps.size();
+        for (size_t iQp = 0; iQp < num_qp; iQp++)
+        {
+            auto w = qps[iQp].weight();
+            vecxd phi = re.basis_functions(iQp);
+            matxd dphi = re.basis_gradients(iQp);
+            vecxd du_phi = dphi.col(0);
+            vecxd dv_phi = dphi.col(1);
+            vecxd dw_phi = dphi.col(2);
+
+            dm.block(0,           0, ed.num_bf, ed.num_bf) += w * phi * du_phi.transpose();
+            dm.block(0,   ed.num_bf, ed.num_bf, ed.num_bf) += w * phi * dv_phi.transpose();
+            dm.block(0, 2*ed.num_bf, ed.num_bf, ed.num_bf) += w * phi * dw_phi.transpose();
+        }
+
+        matxd mass = re.mass_matrix();
+        Eigen::LDLT<matxd> mass_ldlt(mass);
+        if (mass_ldlt.info() != Eigen::Success )
+            throw std::invalid_argument("LDLT failed");
+
+        size_t base_row = iO * ed.num_bf;
+        ed.differentiation_matrix.block(base_row, 0, ed.num_bf, 3*ed.num_bf) =
+            mass_ldlt.solve(dm);
+    }
+}
+
+void
+entity::populate_jacobians(entity_data_cpu<1>& ed) const
+{
+    auto num_elems = num_elements();
+
+    ed.jacobians = matxd::Zero(3*num_elems, 3);
+    for (size_t iT = 0; iT < num_elems; iT++)
+    {
+        auto& pe = physical_cells[iT];
+        ed.jacobians.block(3*iT, 0, 3, 3) = pe.jacobian().inverse().transpose();
+    }
+}
+
+void
+entity::populate_entity_data(entity_data_cpu<1>& ed) const
+{
+    assert(g_order == 1);
+    assert(reference_cells.size() > 0);
+    ed.num_elems = num_elements_per_orientation();
+    ed.num_orientations = num_orientations();
+    ed.num_bf = reference_cells[0].num_basis_functions();
+    ed.dofs_base = 0;
+
+    populate_differentiation_matrices(ed);
+    populate_jacobians(ed);
+
+#if 0
+    if (new_re.m_geometric_order == 1)
+    {
+        new_re.m_
+
+    }
+    else
+    {
+        new_re.m_mass_matrices = matxd::Zero(num_bf*num_qp, 3*num_bf);
+        for (size_t iQp = 0; iQp < num_qp; iQp++)
+        {
+            auto w = new_re.m_quadpoints[iQp].weight();
+            vecxd phi = new_re.m_basis_functions.segment(num_bf*iQp, num_bf);
+            new_re.m_mass_matrices.block(num_bf*iQp, 0, num_bf, num_bf) =
+                w * phi * phi.transpose();
+        }
+    }
 #endif
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+std::pair<std::vector<point_3d>, std::vector<quadrature_point_3d>>
+integrate(const entity& e, int order)
+{
+    auto tag = e.gmsh_tag();
+    auto elemType = e.gmsh_elem_type();
+
+    std::vector<double>     ips;
+    std::vector<double>     iws;
+    gmm::getIntegrationPoints(elemType, quadrature_name(order), ips, iws);
+
+    auto num_qp = iws.size();
+
+    std::vector<double>     ppts;
+    std::vector<double>     jacs;
+    std::vector<double>     dets;
+    gmm::getJacobians(elemType, ips, jacs, dets, ppts, tag);
+
+    std::vector<point_3d> retr;
+
+    for (size_t iQp = 0; iQp < num_qp; iQp++)
+    {
+        double u = ips[3*iQp + 0];
+        double v = ips[3*iQp + 1];
+        double w = ips[3*iQp + 2];
+        point_3d p(u,v,w);
+        retr.push_back(p);
+    }
+
+    std::vector<quadrature_point_3d> retp;
+    retp.reserve(dets.size());
+
+    for (auto& pe : e)
+    {
+        auto g_pos = pe.original_position();
+        auto dets_base = g_pos*num_qp;
+        auto ppts_base = 3*g_pos*num_qp;
+
+        for (size_t iQp = 0; iQp < num_qp; iQp++)
+        {
+            double w = dets[dets_base+iQp]*iws[iQp];
+            double x = ppts[ppts_base + 3*iQp + 0];
+            double y = ppts[ppts_base + 3*iQp + 1];
+            double z = ppts[ppts_base + 3*iQp + 2];
+            point_3d p(x,y,z);
+            quadrature_point_3d q(p,w);
+            retp.push_back(q);
+        }
+    }
+    return std::make_pair(retr, retp);
+}
+

entity.h

+#pragma once
+
+#include <iostream>
+
+#include "reference_element.h"
+#include "physical_element.h"
+
+#include "entity_data.h"
+
+using scalar_function = std::function<double(const point_3d&)>;
+
+enum class entity_ordering {
+    GMSH,
+    BY_ORIENTATION
+};
+
+class entity
+{
+    int     dim;
+    int     tag;
+    int     elemType;
+    int     elemType_2D;
+    int     g_order;
+    int     a_order;
+
+    entity_ordering                     cur_elem_ordering;
+
+    std::vector<reference_element>      reference_cells;
+    std::vector<physical_element>       physical_cells;
+
+    std::vector<reference_element>      reference_faces;
+    std::vector<physical_element>       physical_faces;
+
+    void populate_differentiation_matrices(entity_data_cpu<1>&) const;
+    void populate_jacobians(entity_data_cpu<1>&) const;
+
+public:
+    entity() = delete;
+    entity(const entity&) = delete;
+    entity(entity&&) = default;
+    entity& operator=(const entity&) = delete;
+    entity& operator=(entity&&) = delete;
+
+    entity(const entity_params& ep);
+
+    double                      measure(void) const;
+    vecxd                       project(const scalar_function&) const;
+
+    size_t                      num_elements(void) const;
+    size_t                      num_orientations(void) const;
+    std::vector<size_t>         num_elements_per_orientation(void) const;
+
+    const reference_element&    refelem(const physical_element&) const;
+    const physical_element&     elem(size_t);
+
+    matxd                       mass_matrix(size_t);
+    
+    void                        sort_by_orientation(void);
+    void                        sort_by_gmsh(void);
+
+    int                         gmsh_tag(void) const;
+    int                         gmsh_elem_type(void) const;
+    entity_ordering             current_elem_ordering(void) const;
+
+    void                        populate_entity_data(entity_data_cpu<1>&) const;
+
+
+    std::vector<physical_element>::const_iterator  begin() const { return physical_cells.begin(); }
+    std::vector<physical_element>::const_iterator    end() const { return physical_cells.end(); }
+
+};
+
+std::pair<std::vector<point_3d>, std::vector<quadrature_point_3d>>
+integrate(const entity& e, int order);
+

entity_data.h

+#pragma once
+
+#include "eigen.h"
+
+template<size_t GK>
+struct entity_data_cpu;
+
+template<>
+struct entity_data_cpu<1>
+{
+    std::vector<size_t> num_elems;          /* No. of elements in the entity */
+    size_t              num_orientations;   /* No. of bf. orientations */
+    size_t              num_bf;             /* No. of dofs per elem */
+    size_t              num_fluxes;         /* No. of flux dofs per elem */
+    size_t              num_faces;          /* No. of faces */
+    size_t              dofs_base;
+
+    matxd               mass_matrix;
+    matxd               differentiation_matrix;
+    matxd               lifting_matrix;
+    matxd               jacobians;
+};
+

gmsh_io.cpp

@@ -3,6 +3,7 @@
 #include <cassert>
 
 #include "gmsh_io.h"
+#include "entity.h"
 
 std::string
 quadrature_name(int order)
@@ -15,7 +16,9 @@ quadrature_name(int order)
 std::string
 basis_func_name(int order)
 {
+    (void) order;
     std::stringstream ss;
+    //ss << "Lagrange";
     ss << "H1Legendre" << order;
     return ss.str();
 }
@@ -23,7 +26,9 @@ basis_func_name(int order)
 std::string
 basis_grad_name(int order)
 {
+    (void) order;
     std::stringstream ss;
+    //ss << "GradLagrange";
     ss << "GradH1Legendre" << order;
     return ss.str();
 }
@@ -31,42 +36,93 @@ basis_grad_name(int order)
 model::model()
     : geometric_order(1), approximation_order(1)
 {
-    read_gmsh_entities();
+    remesh(geometric_order);
 }
 
 model::model(int pgo, int pao)
     : geometric_order(pgo), approximation_order(pao)
 {
-    assert(geometric_order >= 1 and geometric_order <= 3);
+    assert(geometric_order > 0);
     assert(approximation_order > 0);
 
-    read_gmsh_entities();
+    remesh(geometric_order);
 }
 
 model::model(const char *filename)
     : geometric_order(1), approximation_order(1)
 {
     gmsh::open( std::string(filename) );
-    read_gmsh_entities();
+    remesh(geometric_order);
 }
 
 model::model(const char *filename, int pgo, int pao)
     : geometric_order(pgo), approximation_order(pao)
 {
-    assert(geometric_order >= 1 and geometric_order <= 3);
+    assert(geometric_order > 0);
     assert(approximation_order > 0);
 
     gmsh::open( std::string(filename) );
-    read_gmsh_entities();
+    remesh(geometric_order);
+}
+
+model::~model()
+{
+    gmsh::clear();
+}
+
+void
+model::remesh()
+{
+    gmm::generate( DIMENSION(3) );
+    gmm::setOrder( geometric_order );
+    entities.clear();
+    import_gmsh_entities();
+}
+
+void
+model::remesh(int pgo)
+{
+    geometric_order = pgo;
+    remesh();
+}
+
+std::vector<entity>::const_iterator
+model::begin() const
+{
+    return entities.begin();
+}
+
+std::vector<entity>::const_iterator
+model::end() const
+{
+    return entities.begin();
+}
+
+std::vector<entity>::iterator
+model::begin()
+{
+    return entities.begin();
+}
+
+std::vector<entity>::iterator
+model::end()
+{
+    return entities.begin();
+}
+
+entity&
+model::entity_at(size_t pos)
+{
+    return entities.at(pos);
 }
 
 void
-model::read_gmsh_entities(void)
+model::import_gmsh_entities(void)
 {
-    gvp_t entities;
-    gm::getEntities(entities, DIMENSION(3));
+    gvp_t gmsh_entities;
+    gm::getEntities(gmsh_entities, DIMENSION(3));
 
-    for (auto [dim, tag] : entities)
+    for (auto [dim, tag] : gmsh_entities)
     {
         std::vector<int> elemTypes;
         gmm::getElementTypes(elemTypes, dim, tag);
@@ -77,7 +133,9 @@ model::read_gmsh_entities(void)
 
         for (auto& elemType : elemTypes)
         {
-            
+            entity e({.dim = dim, .tag = tag, .etype = elemType,
+                      .aorder = approximation_order, .gorder = geometric_order});
+            entities.push_back( std::move(e) );
         }
     }
 }

gmsh_io.h

@@ -4,7 +4,7 @@
 
 #include <gmsh.h>
 
-#include "reference_element.h"
+#include "entity.h"
 
 namespace g     = gmsh;
 namespace go    = gmsh::option;
@@ -28,12 +28,29 @@ class model
     int                             geometric_order;
     int                             approximation_order;
 
-    void    read_gmsh_entities(void);
+
+    std::vector<entity>             entities;
+
+
+    void    import_gmsh_entities(void);
 
 public:
     model();
     model(int, int);
     model(const char *);
     model(const char *, int, int);
+    ~model();
+
+    void    remesh(void);
+    void    remesh(int);
+
+    std::vector<entity>::const_iterator begin() const;
+    std::vector<entity>::const_iterator end() const;
+    std::vector<entity>::iterator begin();
+    std::vector<entity>::iterator end();
+
+    entity&     entity_at(size_t);
+    entity      entity_at(size_t) const;
+
 };
 

kernels_cpu.cpp

+#include "types.h"
+#include "entity_data.h"
+
+void compute_field_derivatives(const entity_data_cpu<1>& ed,
+    const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz)
+{
+    /* The elements must be ordered by orientation -> entity::sort_by_orientation() */
+    size_t orient_elem_base = 0;
+    for (size_t iO = 0; iO < ed.num_orientations; iO++)
+    {
+        size_t num_elems = ed.num_elems[iO];
+        size_t orient_base = orient_elem_base * ed.num_bf;
+
+        for (size_t iT = 0; iT < num_elems; iT++)
+        {
+            size_t ent_iT = orient_elem_base + iT;
+            size_t dofofs = ed.dofs_base + orient_base + iT*ed.num_bf;
+
+            for (size_t odof = 0; odof < ed.num_bf; odof++)
+            {
+                double acc_df_dx = 0.0;
+                double acc_df_dy = 0.0;
+                double acc_df_dz = 0.0;
+
+                for (size_t refsys_dir = 0; refsys_dir < 3; refsys_dir++)
+                {
+                    for (size_t idof = 0; idof < ed.num_bf; idof++)
+                    {
+                        auto dm_row = iO*ed.num_bf + odof;
+                        auto dm_col = refsys_dir*ed.num_bf + idof;
+                        double D = ed.differentiation_matrix(dm_row, dm_col);
+                        acc_df_dx += ed.jacobians(3*ent_iT + 0, refsys_dir) * D * f(dofofs+idof);
+                        acc_df_dy += ed.jacobians(3*ent_iT + 1, refsys_dir) * D * f(dofofs+idof);
+                        acc_df_dz += ed.jacobians(3*ent_iT + 2, refsys_dir) * D * f(dofofs+idof);
+                    }
+                }
+
+                df_dx(dofofs+odof) = acc_df_dx;
+                df_dy(dofofs+odof) = acc_df_dy;
+                df_dz(dofofs+odof) = acc_df_dz;
+            }
+        }
+
+        orient_elem_base += num_elems;
+    }
+}
\ No newline at end of file

kernels_cpu.h

+#pragma once
+
+void compute_field_derivatives(const entity_data_cpu<1>& ed,
+    const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz);
+

main.cpp

@@ -3,6 +3,7 @@
 #include <unistd.h>
 
 #include "gmsh_io.h"
+#include "silo_output.hpp"
 
 static void usage(const char *progname)
 {
@@ -13,9 +14,11 @@ static void usage(const char *progname)
 int main(int argc, char *argv[])
 {
     const char *model_filename = nullptr;
+    int geometric_order = 1;
+    int approximation_order = 1;
 
     int ch;
-    while ( (ch = getopt(argc, argv, "m:")) != -1 )
+    while ( (ch = getopt(argc, argv, "m:o:O:")) != -1 )
     {
         switch(ch)
         {
@@ -23,12 +26,25 @@ int main(int argc, char *argv[])
             model_filename = optarg;
             break;
 
+        case 'o':
+            approximation_order = atoi(optarg);
+            /* TODO: check approximation order */
+            break;
+
+        case 'O':
+            geometric_order = atoi(optarg);
+            /* TODO: check geometric order */
+            break;
+
         default:
             usage(argv[0]);
             return 1;
         }
     }
 
+    argc -= optind;
+    argv += optind;
+
     if (!model_filename)
     {
         std::cout << "Model filename (-m) not specified" << std::endl;
@@ -37,9 +53,13 @@ int main(int argc, char *argv[])
 
     gmsh::initialize(argc, argv);
 
-    
+    geometric_order = approximation_order;
+    model mod(model_filename, geometric_order, approximation_order);
 
-    model mod;
+    silo silo;
+    silo.create_db("test.silo");
+    silo.import_mesh_from_gmsh();
+    silo.write_mesh();
 
     gmsh::finalize();
     return 0;

physical_element.cpp

+#include <iostream>
+#include <cassert>
 
 #include "physical_element.h"
 #include "gmsh_io.h"
+#include "sgr.hpp"
 
 physical_element::physical_element()
 {}
 
-physical_elements_factory::physical_elements_factory(int pdim, int ptag, int pElemType)
-    : dim(pdim), tag(ptag), elemType(pElemType)
+size_t
+physical_element::original_position(void) const
 {
-    
+    return m_original_position;
 }
 
+int
+physical_element::dimension(void) const
+{
+    return m_dimension;
+}
+
+int
+physical_element::orientation(void) const
+{
+    return m_orientation;
+}
+
+double
+physical_element::measure(void) const
+{
+    return m_measure;
+}
+
+vec_quadpt_3d
+physical_element::integration_points(void) const
+{
+    return m_phys_quadpoints;
+}
+
+/* Return only one determinant: this is used for geometric order 1
+ * where the determinants and the jacobians are constant */
+double
+physical_element::determinant(void) const
+{
+    assert(m_geometric_order == 1);
+    return m_determinants[0];
+}
+
+/* Return all the determinants in all the integration points: this
+ * is for elements with geometric order > 1 */
+vec_double
+physical_element::determinants(void) const
+{
+    assert(m_geometric_order >= 1);
+    return m_determinants;
+}
+
+/* Return only one jacobian: this is used for geometric order 1
+ * where the determinants and the jacobians are constant */
+mat3d
+physical_element::jacobian(void) const
+{
+    assert(m_geometric_order == 1);
+    return m_jacobians[0];
+}
+
+/* Return all the jacobians in all the integration points: this
+ * is for elements with geometric order > 1 */
+vec_mat3d
+physical_element::jacobians(void) const
+{
+    assert(m_geometric_order >= 1);
+    return m_jacobians;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+physical_elements_factory::physical_elements_factory(const entity_params& ep)
+    : dim(ep.dim), tag(ep.tag), elemType(ep.etype), geom_order(ep.gorder),
+      approx_order(ep.aorder)
+{}
+
 std::vector<physical_element>
 physical_elements_factory::get_elements()
 {
-    auto order = 2;
+    std::vector<size_t> elemTags;
+    std::vector<size_t> elemNodesTags;
+    gmm::getElementsByType(elemType, elemTags, elemNodesTags, tag);
+    auto num_elems = elemTags.size();
+    auto num_nodes = elemNodesTags.size()/elemTags.size();
+    assert( elemTags.size()*num_nodes == elemNodesTags.size() );
+
     std::vector<double>     ips;
     std::vector<double>     iws;
-    gmm::getIntegrationPoints(elemType, quadrature_name(2*order), ips, iws);
+    gmm::getIntegrationPoints(elemType, quadrature_name(2*approx_order), ips, iws);
 
     auto num_gf = iws.size();
 
@@ -26,50 +116,70 @@ physical_elements_factory::get_elements()
     std::vector<double>     dets;
     gmm::getJacobians(elemType, ips, jacs, dets, ppts, tag);
 
-    auto num_elems = dets.size()/num_gf;
     assert( num_gf * num_elems == dets.size() );
 
+    std::vector<int> orientations;
+    gmm::getBasisFunctionsOrientationForElements(elemType, basis_func_name(approx_order),
+        orientations, tag);
+    assert(orientations.size() == num_elems);
+
     std::vector<physical_element> ret;
     ret.reserve(num_elems);
 
     for (size_t elem = 0; elem < num_elems; elem++)
     {
         physical_element new_pe;
+        new_pe.m_geometric_order = geom_order;
+        new_pe.m_approximation_order = approx_order;
         new_pe.m_original_position = elem;
         new_pe.m_dimension = dim;
         new_pe.m_parent_entity_tag = tag;
+        new_pe.m_orientation = orientations[elem];
+        new_pe.m_element_tag = elemTags[elem];
+        new_pe.m_measure = 0.0;
+
+        new_pe.m_node_tags.resize(num_nodes);
+        for (size_t i = 0; i < num_nodes; i++)
+            new_pe.m_node_tags[i] = elemNodesTags[num_nodes*elem+i];
+
+        auto elem_base = elem*num_gf;
 
         for (size_t gf = 0; gf < num_gf; gf++)
-        {
+        {   
+            auto gf_offset = elem_base + gf;
+ 
             const auto JSIZE = 3*3; /* Jacobian matrix size */
-            auto jacs_base = elem*num_gf*JSIZE + gf*JSIZE;
+            auto jacs_base = gf_offset*JSIZE;
             assert(jacs_base+8 < jacs.size());
 
             /* GMSH returns jacobians by column */
             Eigen::Matrix3d jac;
-            jac(0,0) = jacs[jacs_base+0];
-            jac(0,1) = jacs[jacs_base+1];
-            jac(0,2) = jacs[jacs_base+2];
-            jac(1,0) = jacs[jacs_base+3];
-            jac(1,1) = jacs[jacs_base+4];
-            jac(1,2) = jacs[jacs_base+5];
-            jac(2,0) = jacs[jacs_base+6];
-            jac(2,1) = jacs[jacs_base+7];
-            jac(2,2) = jacs[jacs_base+8];
+            jac(0,0) = jacs[jacs_base+0];   /* dx/du */
+            jac(1,0) = jacs[jacs_base+1];   /* dy/du */
+            jac(2,0) = jacs[jacs_base+2];   /* dz/du */
+            jac(0,1) = jacs[jacs_base+3];   /* dx/dv */
+            jac(1,1) = jacs[jacs_base+4];   /* dy/dv */
+            jac(2,1) = jacs[jacs_base+5];   /* dz/dv */
+            jac(0,2) = jacs[jacs_base+6];   /* dx/dw */
+            jac(1,2) = jacs[jacs_base+7];   /* dy/dw */
+            jac(2,2) = jacs[jacs_base+8];   /* dz/dw */
             new_pe.m_jacobians.push_back(jac);
 
-            auto dets_ofs = elem*num_gf + gf;
+            auto dets_ofs = gf_offset;
             assert( dets_ofs < dets.size() );
             new_pe.m_determinants.push_back(dets[dets_ofs]);
+            new_pe.m_measure += dets[dets_ofs]*iws[gf];
 
             const auto PSIZE = 3;
-            auto pts_base = elem*num_gf*PSIZE + gf*PSIZE;
+            auto pts_base = gf_offset*PSIZE;
             assert(pts_base + 2 < ppts.size());
             point_3d pt(ppts[pts_base+0], ppts[pts_base+1], ppts[pts_base+2]);
-            new_pe.m_phys_quadpoints.push_back(pt);
+            quadrature_point_3d qpt(pt, iws[gf]*dets[dets_ofs]);
+            new_pe.m_phys_quadpoints.push_back(qpt);
         }
+
+        ret.push_back( std::move(new_pe) );
     }
 
-    return std::vector<physical_element>();
+    return ret;
 }
-

physical_element.h

 #pragma once
 
 #include "point.hpp"
-#include "eigen.h"
-
-using vec_int = std::vector<int>;
-using vec_double = std::vector<double>;
-using vec_point_3d = std::vector<point_3d>;
-using vec_mat3d = eigen_compatible_stdvector<Eigen::Matrix3d>;
+#include "types.h"
 
 class physical_element
 {
-    size_t          m_original_position;
-    int             m_dimension;
-    int             m_orientation;
-    int             m_parent_entity_tag;
-    int             m_element_tag;
-    vec_int         m_node_tags;
-    vec_double      m_determinants;
-    vec_mat3d       m_jacobians;
-    vec_point_3d    m_phys_quadpoints;
-    point_3d        m_barycenter;
+    size_t          m_original_position;    /* Position in GMSH ordering (relative to entity) */
+    int             m_dimension;            /* Geometric dimension */
+    int             m_orientation;          /* GMSH orientation */
+    int             m_parent_entity_tag;    /* Tag of the parent entity */
+    size_t          m_element_tag;          /* Tag of the element */
+    vec_size_t      m_node_tags;            /* Tags of the element nodes */
+    vec_double      m_determinants;         /* SIGNED determinants in the quadrature points */
+    vec_mat3d       m_jacobians;            /* Jacobians in the quadrature points */
+    vec_quadpt_3d   m_phys_quadpoints;      /* Quadrature points in the physical element */
+    point_3d        m_barycenter;           /* Element barycenter */
+    double          m_measure;              /* Element measure */
+    int             m_geometric_order;
+    int             m_approximation_order;
 
 public:
     physical_element();
 
+    size_t          original_position(void) const;
+    int             dimension(void) const;
+    int             orientation(void) const;
+    int             parent_entity_tag(void) const;
+    int             element_tag(void) const;
+    vec_int         node_tags(void) const;
+
+    double          determinant(void) const;
+    vec_double      determinants(void) const;
+
+    mat3d           jacobian(void) const;
+    vec_mat3d       jacobians(void) const;
+
+    vec_quadpt_3d   integration_points(void) const;
+    point_3d        barycenter(void) const;
+    double          measure(void) const;
+
+    /* This class is not intended to be constructed
+     * by the user but by the appropriate factory */
     friend class physical_elements_factory;
 };
 
@@ -32,10 +49,24 @@ class physical_elements_factory
     int         dim;
     int         tag;
     int         elemType;
+    int         geom_order;
+    int         approx_order;
 
 public:
-    physical_elements_factory(int, int, int);
+    physical_elements_factory(const entity_params& ep);
 
     std::vector<physical_element> get_elements();
 };
 
+
+
+
+
+
+
+
+
+
+
+
+

point.hpp

@@ -217,6 +217,7 @@ class quadrature_point
     point<T, DIM>   p;
     T               w;
 
+public:
     quadrature_point()
     {}
 

reference_element.cpp

+#include "gmsh_io.h"
+#include "reference_element.h"
+
+reference_element::reference_element()
+{}
+
+size_t
+reference_element::num_basis_functions(void) const
+{
+    return m_num_bf;
+}
+
+size_t
+reference_element::num_integration_points(void) const
+{
+    return m_quadpoints.size();
+}
+
+vecxd
+reference_element::basis_functions(size_t iQp) const
+{
+    assert(iQp < num_integration_points());
+    return m_basis_functions.segment(m_num_bf*iQp, m_num_bf);
+}
+
+matxd
+reference_element::basis_gradients(size_t iQp) const
+{
+    assert(iQp < num_integration_points());
+    return m_basis_gradients.block(m_num_bf * iQp, 0, m_num_bf, 3);
+}
+
+vecxd
+reference_element::basis_functions(const point_3d& p) const
+{
+    std::vector<double> ip;
+    ip.push_back(p.x());
+    ip.push_back(p.y());
+    ip.push_back(p.z());
+
+    int bf_nc, bf_no;
+    std::vector<double> basisFunctions;
+    gmm::getBasisFunctions(m_gmsh_elem_type, ip, basis_func_name(m_approx_order), bf_nc,
+        basisFunctions, bf_no);
+
+    assert( basisFunctions.size() == m_num_bf*bf_no );
+    auto bf_ofs = m_num_bf * m_orientation;
+
+    vecxd ret = vecxd::Zero(m_num_bf);
+    for (size_t i = 0; i < m_num_bf; i++)
+        ret(i) = basisFunctions[bf_ofs + i];
+
+    return ret;
+}
+
+matxd
+reference_element::mass_matrix(void) const
+{
+    assert(m_geometric_order == 1);
+    return m_mass_matrix;
+}
+
+matxd
+reference_element::mass_matrix(size_t iQp) const
+{
+    assert(m_geometric_order > 1);
+    return m_mass_matrices.block(m_num_bf*iQp, 0, m_num_bf, m_num_bf);
+}
+
+matxd
+reference_element::mass_matrices(void) const
+{
+    assert(m_geometric_order > 1);
+    return m_mass_matrices;
+}
+
+vec_quadpt_3d
+reference_element::integration_points(void) const
+{
+    return m_quadpoints;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+reference_elements_factory::reference_elements_factory(const entity_params& ep)
+    : dim(ep.dim), tag(ep.tag), elemType(ep.etype), geom_order(ep.gorder),
+      approx_order(ep.aorder)
+{}
+
+std::vector<reference_element>
+reference_elements_factory::get_elements()
+{
+    /* Get integration points and weights on refelem */
+    std::vector<double>     ips;
+    std::vector<double>     iws;
+    gmm::getIntegrationPoints(elemType, quadrature_name(2*approx_order), ips, iws);
+    assert(ips.size() == 3*iws.size());
+
+    auto num_qp = iws.size();
+
+    int bf_nc, num_bfun_orient;
+    std::vector<double> basisFunctions;
+    gmm::getBasisFunctions(elemType, ips, basis_func_name(approx_order), bf_nc,
+        basisFunctions, num_bfun_orient);
+    assert(bf_nc == 1);
+    
+    auto num_bf = basisFunctions.size() / (num_qp * num_bfun_orient);
+    assert( num_bf * num_qp * num_bfun_orient == basisFunctions.size() );
+
+    int num_bgrad_orient;
+    std::vector<double> basisGradients;
+    gmm::getBasisFunctions(elemType, ips, basis_grad_name(approx_order), bf_nc,
+        basisGradients, num_bgrad_orient);
+    assert(bf_nc == 3);
+
+    assert(num_bfun_orient == num_bgrad_orient);
+
+    std::vector<reference_element> ret;
+
+    for (size_t iO = 0; iO < num_bfun_orient; iO++)
+    {
+        /* Create one reference element per orientation */
+        reference_element new_re;
+        new_re.m_approx_order = approx_order;
+        new_re.m_geometric_order = geom_order;
+        new_re.m_dimension = dim;
+        new_re.m_orientation = iO;
+        new_re.m_parent_entity_tag = tag;
+        new_re.m_basis_functions = vecxd::Zero(num_qp*num_bf);
+        new_re.m_basis_gradients = matxd::Zero(num_bf*num_qp, 3);
+        new_re.m_num_bf = num_bf;
+        new_re.m_gmsh_elem_type = elemType;
+
+        auto bf_base = num_bf * num_qp * iO;
+        auto bg_base = 3 * bf_base;
+        for (size_t iQp = 0; iQp < num_qp; iQp++)
+        {
+            for (size_t ibf = 0; ibf < num_bf; ibf++)
+            {
+                new_re.m_basis_functions(num_bf*iQp + ibf) =
+                    basisFunctions[bf_base + num_bf*iQp + ibf];
+                new_re.m_basis_gradients(num_bf*iQp + ibf, 0) =
+                    basisGradients[bg_base + 3*num_bf*iQp + 3*ibf + 0];
+                new_re.m_basis_gradients(num_bf*iQp + ibf, 1) =
+                    basisGradients[bg_base + 3*num_bf*iQp + 3*ibf + 1];
+                new_re.m_basis_gradients(num_bf*iQp + ibf, 2) =
+                    basisGradients[bg_base + 3*num_bf*iQp + 3*ibf + 2];
+            }
+        }
+
+        for (size_t iQp = 0; iQp < num_qp; iQp++)
+        {
+            auto weight = iws[iQp];
+            auto qpt_u = ips[3*iQp+0];
+            auto qpt_v = ips[3*iQp+1];
+            auto qpt_w = ips[3*iQp+2];
+            auto pt = point_3d(qpt_u, qpt_v, qpt_w);
+            new_re.m_quadpoints.push_back( quadrature_point_3d(pt, weight) );       
+        }
+
+        if (new_re.m_geometric_order == 1)
+        {
+            new_re.m_mass_matrix = matxd::Zero(num_bf, num_bf);
+            for (size_t iQp = 0; iQp < num_qp; iQp++)
+            {
+                auto w = new_re.m_quadpoints[iQp].weight();
+                vecxd phi = new_re.m_basis_functions.segment(num_bf*iQp, num_bf);
+                new_re.m_mass_matrix += w * phi * phi.transpose();
+            }
+        }
+        else
+        {
+            new_re.m_mass_matrices = matxd::Zero(num_bf*num_qp, num_bf);
+            for (size_t iQp = 0; iQp < num_qp; iQp++)
+            {
+                auto w = new_re.m_quadpoints[iQp].weight();
+                vecxd phi = new_re.m_basis_functions.segment(num_bf*iQp, num_bf);
+                new_re.m_mass_matrices.block(num_bf*iQp, 0, num_bf, num_bf) =
+                    w * phi * phi.transpose();
+            }
+        }
+
+        ret.push_back(std::move(new_re));
+    }
+
+    return ret;
+}
+

reference_element.h

 #pragma once
 
-class reference_cell
-{};
+#include "types.h"
 
-class reference_face
-{};
+/*
+static size_t orientations[] = {
+    1, 0, 1, 4,     1, 0, 0, 2,     0, 0, 1, 5,
+    1, 1, 0, 0,     0, 1, 1, 3,     0, 1, 0, 1,
+    3, 2, 1, 4,     3, 2, 0, 2,     2, 0, 3, 5,
+    1, 3, 2, 0,     2, 1, 3, 3,     0, 3, 2, 1,
+    5, 2, 3, 4,     3, 4, 2, 2,     4, 2, 3, 5,
+    3, 5, 2, 0,     2, 3, 5, 3,     2, 3, 4, 1,
+    5, 4, 5, 4,     5, 4, 4, 2,     4, 4, 5, 5,
+    5, 5, 4, 0,     4, 5, 5, 3,     4, 5, 4, 1
+};
+
+#define FACE_ORIENTATION(co, fn) ( orientations[4*co + fn] )
+*/
+
+class reference_element
+{
+    int                 m_approx_order;
+    int                 m_geometric_order;
+    int                 m_dimension;
+    int                 m_orientation;
+    int                 m_gmsh_elem_type;
+    int                 m_parent_entity_tag;
+    size_t              m_num_bf;
+    vecxd               m_basis_functions;
+    matxd               m_basis_gradients;
+    matxd               m_mass_matrix;
+    matxd               m_mass_matrices;
+    vec_quadpt_3d       m_quadpoints;
+
+public:
+    reference_element();
+
+    size_t          num_basis_functions(void) const;
+    size_t          num_integration_points(void) const;
+
+    vec_quadpt_3d   integration_points(void) const;
+
+    vecxd           basis_functions(size_t iQp) const;
+    matxd           basis_gradients(size_t iQp) const;
+
+    vecxd           basis_functions(const point_3d&) const;
+
+    matxd           mass_matrix(void) const;
+    matxd           mass_matrix(size_t iQp) const;
+    matxd           mass_matrices(void) const;
+
+    /* This class is not intended to be constructed
+     * by the user but by the appropriate factory */
+    friend class reference_elements_factory;
+};
+
+class reference_elements_factory
+{
+    int         dim;
+    int         tag;
+    int         elemType;
+    int         geom_order;
+    int         approx_order;
+
+public:
+    reference_elements_factory(const entity_params& ep);
+
+    std::vector<reference_element> get_elements();
+};
 

silo_output.hpp

+#pragma once
+
+#include <cassert>
+#include <silo.h>
+#include "gmsh.h"
+
+#define INVALID_OFS ((size_t) (~0))
+
+class silo
+{
+    std::vector<size_t>     node_tag2ofs;
+    std::vector<double>     node_coords_x;
+    std::vector<double>     node_coords_y;
+    std::vector<double>     node_coords_z;
+
+    std::vector<int>        nodelist;
+    int                     num_cells;
+
+    DBfile                  *m_siloDb;
+
+    void gmsh_get_nodes()
+    {
+        std::vector<size_t>     nodeTags;
+        std::vector<double>     coords;
+        std::vector<double>     paraCoords;
+
+        gmsh::model::mesh::getNodes(nodeTags, coords, paraCoords, -1, -1, true, false);
+
+        auto maxtag_pos = std::max_element(nodeTags.begin(), nodeTags.end());
+        
+        size_t maxtag = 0;
+        if (maxtag_pos != nodeTags.end())
+            maxtag = (*maxtag_pos)+1;
+
+        node_coords_x.resize( nodeTags.size() );
+        node_coords_y.resize( nodeTags.size() );
+        node_coords_z.resize( nodeTags.size() );
+
+        for (size_t i = 0; i < nodeTags.size(); i++)
+        {
+            node_coords_x[i] = coords[3*i+0];
+            node_coords_y[i] = coords[3*i+1];
+            node_coords_z[i] = coords[3*i+2];
+        }
+
+        node_tag2ofs.resize(maxtag, INVALID_OFS);
+        for (size_t i = 0; i < nodeTags.size(); i++)
+            node_tag2ofs.at( nodeTags[i] ) = i;
+    }
+
+    void gmsh_get_elements()
+    {
+        gmsh::vectorpair entities;
+        num_cells = 0;
+        gmsh::model::getEntities(entities, 3);
+        for (auto [dim, tag] : entities)
+        {
+            std::vector<int> elemTypes;
+            gmsh::model::mesh::getElementTypes(elemTypes, dim, tag);
+            for (auto& elemType : elemTypes)
+            {
+                std::vector<size_t> elemTags;
+                std::vector<size_t> elemNodeTags;
+                gmsh::model::mesh::getElementsByType(elemType, elemTags, elemNodeTags, tag);
+                auto nodesPerElem = elemNodeTags.size()/elemTags.size();
+                assert( elemTags.size() * nodesPerElem == elemNodeTags.size() );
+            
+                for (size_t i = 0; i < elemTags.size(); i++)
+                {
+                    auto base = nodesPerElem * i;
+
+                    auto node0_tag = elemNodeTags[base + 0];
+                    assert(node0_tag < node_tag2ofs.size());
+                    auto node0_ofs = node_tag2ofs[node0_tag];
+                    assert(node0_ofs != INVALID_OFS);
+                    nodelist.push_back(node0_ofs + 1);
+
+                    auto node1_tag = elemNodeTags[base + 1];
+                    assert(node1_tag < node_tag2ofs.size());
+                    auto node1_ofs = node_tag2ofs[node1_tag];
+                    assert(node1_ofs != INVALID_OFS);
+                    nodelist.push_back(node1_ofs + 1);
+
+                    auto node2_tag = elemNodeTags[base + 2];
+                    assert(node2_tag < node_tag2ofs.size());
+                    auto node2_ofs = node_tag2ofs[node2_tag];
+                    assert(node2_ofs != INVALID_OFS);
+                    nodelist.push_back(node2_ofs + 1);
+
+                    auto node3_tag = elemNodeTags[base + 3];
+                    assert(node3_tag < node_tag2ofs.size());
+                    auto node3_ofs = node_tag2ofs[node3_tag];
+                    assert(node3_ofs != INVALID_OFS);
+                    nodelist.push_back(node3_ofs + 1);
+                }
+
+                num_cells += elemTags.size();
+            }
+        }
+        std::cout << nodelist.size() << " " << num_cells << std::endl;
+
+        assert(nodelist.size() == 4*num_cells);
+    }
+
+public:
+    silo()
+    {}
+
+    void import_mesh_from_gmsh()
+    {
+        gmsh_get_nodes();
+        gmsh_get_elements();        
+    }
+
+    bool create_db(const std::string& dbname)
+    {
+        m_siloDb = DBCreate(dbname.c_str(), DB_CLOBBER, DB_LOCAL, NULL, DB_HDF5);
+        if (m_siloDb)
+            return true;
+
+        std::cout << "SILO: Error creating database" << std::endl;
+        return false;
+    }
+
+    bool write_mesh(double time = -1.0, int cycle = 0)
+    {
+        DBoptlist   *optlist = nullptr;
+
+        if (time >= 0)
+        {
+            optlist = DBMakeOptlist(2);
+            DBAddOption(optlist, DBOPT_CYCLE, &cycle);
+            DBAddOption(optlist, DBOPT_DTIME, &time);
+        }
+
+        double *coords[] = {node_coords_x.data(), node_coords_y.data(), node_coords_z.data()};
+
+        int lnodelist = nodelist.size();
+
+        int shapetype[] = { DB_ZONETYPE_TET };
+        int shapesize[] = {4};
+        int shapecounts[] = { num_cells };
+        int nshapetypes = 1;
+        int nnodes = node_coords_x.size();
+        int nzones = num_cells;
+        int ndims = 3;
+
+        DBPutZonelist2(m_siloDb, "zonelist", nzones, ndims,
+            nodelist.data(), lnodelist, 1, 0, 0, shapetype, shapesize,
+            shapecounts, nshapetypes, NULL);
+
+        DBPutUcdmesh(m_siloDb, "mesh", ndims, NULL, coords, nnodes, nzones,
+            "zonelist", NULL, DB_DOUBLE, optlist);
+
+        if (optlist)
+            DBFreeOptlist(optlist);
+
+        return true;
+    }
+
+    bool write_zonal_variable(const std::string& name, const std::vector<double>& data)
+    {
+        if (data.size() != num_cells)
+        {
+            std::cout << "Invalid number of cells" << std::endl;
+            return false;
+        }
+
+        DBPutUcdvar1(m_siloDb, name.c_str(), "mesh", data.data(), data.size(), NULL,
+                     0, DB_DOUBLE, DB_ZONECENT, NULL);
+        
+        return true;
+    }
+
+    void close_db()
+    {
+        if (m_siloDb)
+            DBClose(m_siloDb);
+
+        m_siloDb = nullptr;
+    }
+};

test.cpp

+#include <iostream>
+
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "silo_output.hpp"
+#include "sgr.hpp"
+#include "entity_data.h"
+#include "kernels_cpu.h"
+
+using namespace sgr;
+
+#define SPHERE_RADIUS 0.2
+
+static void
+resize_mesh(double mesh_h)
+{
+    gvp_t vp;
+    gm::getEntities(vp);
+    gmm::setSize(vp, mesh_h);
+}
+
+static void
+make_geometry(int order, double mesh_h)
+{
+    gm::add("difftest");
+
+    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();
+
+    resize_mesh(mesh_h);
+}
+
+
+int test_basics(int geometric_order, int approximation_order)
+{
+    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 << "Volumes by reference quadrature points and determinants" << std::endl;
+    std::vector<double> t1_vals { 0.07, 0.00145, 0.00021, 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.024, 5.02e-5, 7.003e-06, 6.35e-7 };
+    COMPARE_VALUES_RELATIVE("Entity 1 vol r+d", mod.entity_at(1).measure(), vol_e1, t2_vals[idx]);
+
+    /*****************************************/
+    std::cout << "Volumes by physical quadrature points" << 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 << "Volumes by jacobians" << std::endl;
+    double vol_e0_j = 0.0;
+    for (const auto& e : e0)
+    {
+        const auto& js = e.jacobians();
+        auto& re = e0.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.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 << "Volumes by integrating f(x) = 1 over the domain" << 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_elements(); iT++)
+    {
+        auto& pe = e0.elem(iT);
+        auto& re = e0.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_elements(); iT++)
+    {
+        auto& pe = e1.elem(iT);
+        auto& re = e1.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;
+}
+
+#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_elements(); iT++)
+        {
+            auto& pe = e0.elem(iT);
+            auto& re = e0.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) );
+    }
+
+    for (auto& error : errors)
+        std::cout << error << " ";
+    std::cout << std::endl;
+
+    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 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()); 
+    };
+
+    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<1> ed;
+        e0.populate_entity_data(ed);
+
+        compute_field_derivatives(ed, Pf_e0, df_dx_e0, df_dy_e0, df_dz_e0);
+
+        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_elements(); iT++)
+        {
+            auto& pe = e0.elem(iT);
+            auto& re = e0.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 << std::sqrt(err_x) << " " << std::sqrt(err_y) << " " << std::sqrt(err_z) << std::endl;
+    }
+
+    return 0;
+}
+
+int main(int argc, char **argv)
+{
+    gmsh::initialize();
+    gmsh::option::setNumber("General.Terminal", 0);
+    gmsh::option::setNumber("Mesh.Algorithm", 1);
+    gmsh::option::setNumber("Mesh.Algorithm3D", 1);
+
+    int failed_tests = 0;
+
+    /*
+    for (size_t go = 1; go < 5; go++)
+    {
+        for (size_t ao = go; ao < 5; ao++)
+        {
+            make_geometry(0, 0.1);
+            failed_tests += test_basics(go, ao);
+        }
+    }
+    */
+   
+    /*
+    for (size_t go = 1; go < 5; go++)
+    {
+        for (size_t ao = go; ao < 5; ao++)
+            failed_tests += test_mass_convergence(go, ao);
+    }
+    */
+
+    for (size_t ao = 1; ao < 5; ao++)
+        failed_tests += test_differentiation_convergence(1, ao);
+
+    /*****************************************/
+    
+    //COMPARE_VALUES_RELATIVE("Entity 0 vol mass", vol_e0_mass, vol_e0, t1_vals[order-1]);
+
+    #if 0
+
+
+
+    #endif
+
+    gmsh::finalize();
+    return failed_tests;
+}
+
+

test.h

+#pragma once
+
+#define COMPARE_VALUES_ABSOLUTE(name, val, ref, tol)                    \
+    std::cout << "Test: " << name << "...";                             \
+    if ( double real_tol = std::abs(val - ref); real_tol < tol )        \
+    { std::cout << greenfg << "PASS" << nofg << " (val = " << val;      \
+      std::cout << ", ref = " << ref << ", tol = " << tol;              \
+      std::cout << ", real_tol = " << real_tol << std::endl; }          \
+    else                                                                \
+    { std::cout << redfg << "FAIL" << nofg << " (val = " << val;        \
+      std::cout << ", ref = " << ref << ", tol = " << tol;              \
+      std::cout << ", real_tol = " << real_tol << std::endl; return 1; }
+
+#define COMPARE_VALUES_RELATIVE(name, val, ref, tol)                                \
+    std::cout << "Test: " << name << "...";                                         \
+    if ( double real_tol = (std::abs(val - ref)/std::abs(ref)); real_tol < tol )    \
+    { std::cout << greenfg << "PASS" << nofg << " (val = " << val;                  \
+      std::cout << ", ref = " << ref << ", tol = " << tol;                          \
+      std::cout << ", real_tol = " << real_tol << std::endl; }                      \
+    else                                                                            \
+    { std::cout << redfg << "FAIL" << nofg << " (val = " << val;        \
+      std::cout << ", ref = " << ref << ", tol = " << tol;              \
+      std::cout << ", real_tol = " << real_tol << std::endl; return 1; }
+

types.h

+#pragma once
+
+#include "point.hpp"
+#include "eigen.h"
+
+using vec_int = std::vector<int>;
+using vec_size_t = std::vector<size_t>;
+using vec_double = std::vector<double>;
+using vec_point_3d = std::vector<point_3d>;
+using vec_quadpt_3d = std::vector<quadrature_point_3d>;
+using mat3d = Eigen::Matrix<double, 3, 3, Eigen::RowMajor>;
+using vec_mat3d = eigen_compatible_stdvector<mat3d>;
+using vecxd = Eigen::Matrix<double, Eigen::Dynamic, 1>;
+using matxd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+struct entity_params
+{
+    int     dim;
+    int     tag;
+    int     etype;
+    int     gorder;
+    int     aorder;
+};
+