Rewrite of DG solver, initial commit.

.gitmodules

+[submodule "contrib/sgr"]
+	path = contrib/sgr
+	url = https://github.com/datafl4sh/sgr
+[submodule "contrib/sol2"]
+	path = contrib/sol2
+	url = https://github.com/ThePHD/sol2

CMakeLists.txt

+cmake_minimum_required(VERSION 3.18)
+project(gmsh_gpu_dg)
+include(CheckLanguage)
+
+
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CUDA_ARCHITECTURES 35)
+set(CMAKE_CUDA_HOST_COMPILER "/usr/bin/g++-6")
+
+# Additional modules path for cmake
+set (CMAKE_MODULE_PATH
+     ${CMAKE_MODULE_PATH}
+     ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+
+
+if(NOT DEFINED HIP_PATH)
+    if(NOT DEFINED ENV{HIP_PATH})
+        set(HIP_PATH "/opt/rocm/hip" CACHE PATH "Path to which HIP has been installed")
+    else()
+        set(HIP_PATH $ENV{HIP_PATH} CACHE PATH "Path to which HIP has been installed")
+    endif()
+endif()
+set(CMAKE_MODULE_PATH "${HIP_PATH}/cmake" ${CMAKE_MODULE_PATH})
+
+######################################################################
+## Find required libraries
+
+find_package(Eigen3 REQUIRED)
+set(LINK_LIBS ${LINK_LIBS} Eigen3::Eigen)
+
+set(GMSH_DIR "/opt/uliege/gmsh" CACHE PATH "GMSH install path")
+if (DEFINED ENV{GMSH_DIR})
+    set (GMSH_DIR $ENV{GMSH_DIR})
+endif()
+
+find_package(SILO)
+if (SILO_FOUND)
+    add_definitions(-DHAVE_SILO)
+    include_directories("${SILO_INCLUDE_DIRS}")
+    set(LINK_LIBS ${LINK_LIBS} ${SILO_LIBRARIES})
+endif()
+
+find_package(Lua REQUIRED)
+include_directories("${LUA_INCLUDE_DIR}")
+set(LINK_LIBS ${LINK_LIBS} ${LUA_LIBRARIES})
+set(SOL2_BUILD_LUA OFF CACHE BOOL "Override sol2 build lua option")
+add_subdirectory(contrib/sol2)
+include_directories(contrib/sol2/include)
+
+######################################################################
+## Set GMSH paths
+
+if (GMSH_DIR STREQUAL "" OR NOT EXISTS ${GMSH_DIR})
+    message(FATAL_ERROR "GMSH not found")
+else()
+    set(GMSH_INCLUDE ${GMSH_DIR}/include)
+    set(GMSH_LIB ${GMSH_DIR}/lib)
+    set(LINK_LIBS ${LINK_LIBS} gmsh)
+endif()
+
+include_directories(${GMSH_INCLUDE})
+link_directories(${GMSH_LIB})
+
+
+
+
+######################################################################
+## Use or not use GPU solvers
+
+option(OPT_ENABLE_GPU_SOLVER "Enable GPU solvers" OFF)
+if (OPT_ENABLE_GPU_SOLVER)
+    add_definitions(-DENABLE_GPU_SOLVER)
+
+    ######################################################################
+    ## CUDA stuff, find toolkit, compiler & NVML
+    find_package(CUDAToolkit)
+    if (CUDAToolkit_FOUND)
+        set(CMAKE_CUDA_COMPILER ${CUDAToolkit_NVCC_EXECUTABLE})
+        enable_language(CUDA)
+        set(CMAKE_CUDA_STANDARD 14)
+        set(CMAKE_CUDA_STANDARD_REQUIRED TRUE)
+        #set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -lineinfo --resource-usage -Xptxas -dlcm=ca")
+        set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -lineinfo --resource-usage -prec-div=true")
+
+        if (TARGET CUDA::cudart)
+            set(HAVE_CUDA TRUE)
+            set(LINK_LIBS ${LINK_LIBS} CUDA::cudart)
+            add_definitions(-DHAVE_CUDA)
+        endif()
+
+        if (TARGET CUDA::nvml)
+            set(HAVE_CUDA_NVML TRUE)
+            set(LINK_LIBS ${LINK_LIBS} CUDA::nvml)
+            add_definitions(-DHAVE_CUDA_NVML)
+        endif()
+	link_directories("/usr/local/cuda/lib64")
+    endif()
+
+    ######################################################################
+    ## Find HIP stuff
+    find_package(HIP QUIET MODULE)
+    if (HIP_FOUND)
+        add_definitions(-DHAVE_HIP)
+        include_directories("/opt/rocm/include")
+        set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+        if (DEFINED ENV{HIP_PLATFORM}) 
+            if ($ENV{HIP_PLATFORM} STREQUAL "nvcc")
+                #set(HIP_NVCC_FLAGS "-ccbin g++-7" ${HIP_NVCC_FLAGS})
+                set(CMAKE_HIP_LINK_EXECUTABLE "${CUDAToolkit_NVCC_EXECUTABLE} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+            endif()
+        endif()
+
+        set(HIP_NVCC_FLAGS "-ccbin g++-8" ${HIP_NVCC_FLAGS})
+        option(OPT_USE_HIP "Use HIP GPU code path" OFF)
+        if (OPT_USE_HIP)
+            add_definitions(-DGPU_USE_HIP)
+            add_definitions(-D__HIP_PLATFORM_HCC__)
+        endif()
+    endif()
+
+endif()
+
+
+######################################################################
+## Compiler identification stuff
+if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+    set(COMPILER_IS_CLANG TRUE)
+elseif (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
+    set(COMPILER_IS_GNU TRUE)
+elseif (CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
+    set(COMPILER_IS_INTEL TRUE)
+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_DEBUG "-O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -g -DNDEBUG")
+set(CMAKE_CXX_FLAGS_RELEASEASSERT "-O3 -g -fpermissive")
+
+option(OPT_ENABLE_OPENMP "Enable OpenMP parallelization" OFF)
+if (OPT_ENABLE_OPENMP)
+    find_package(OpenMP)
+    if (OpenMP_FOUND)
+        set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    endif()
+endif()
+
+if (COMPILER_IS_CLANG OR COMPILER_IS_GNU OR COMPILER_IS_INTEL)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
+endif()
+
+######################################################################
+## Optimization stuff
+
+option(OPT_PREFER_512bit "Prefer 512 bit vectors with AVX512 (only Clang > 10 & GCC)" OFF)
+if (OPT_PREFER_512bit)
+    # https://reviews.llvm.org/D67259
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -mprefer-vector-width=512")
+endif()
+
+option(OPT_AGGRESSIVE_FP "Enable DAZ, FTZ and -ffast-math" ON)
+if (OPT_AGGRESSIVE_FP)
+    add_definitions(-DDISALLOW_DENORMALS)
+    if (COMPILER_IS_CLANG OR COMPILER_IS_GNU)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffast-math")
+    endif()
+endif()
+
+option(OPT_VECTORIZER_REMARKS "Enable vectorizer remarks" OFF)
+if (OPT_VECTORIZER_REMARKS)
+    if (COMPILER_IS_CLANG)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Rpass=loop-vectorize -Rpass-missed=loop-vectorize -Rpass-analysis=loop-vectorize")
+    endif()
+
+    if (COMPILER_IS_INTEL)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -qopt-report-phase=vec -qopt-report=2")
+    endif()
+
+    if (COMPILER_IS_GNU)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopt-info-vec-optimized")
+    endif()
+endif()
+
+option(OPT_MEASURE_PERFORMANCE "Enable performance measurement (produces a ton of output)" OFF)
+if (OPT_MEASURE_PERFORMANCE)
+    add_definitions(-DMEASURE_PERFORMANCE)
+endif()
+
+option(OPT_CRASH_ON_NAN "Crash code if NaNs are generated" OFF)
+if (OPT_CRASH_ON_NAN)
+    add_definitions(-DCRASH_ON_NAN)
+endif()
+
+######################################################################
+## Development & debug stuff
+
+if (COMPILER_IS_CLANG OR COMPILER_IS_GNU)
+    option(OPT_ASAN "Enable Address Sanitizer" OFF)
+    if (OPT_ASAN)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fsanitize=address")
+    endif()
+endif()
+
+option(OPT_ENABLE_CUDA_DEBUG_GUARDS "Enable CUDA debug guards")
+if (OPT_ENABLE_CUDA_DEBUG_GUARDS AND HAVE_CUDA)
+    add_definitions(-DCUDA_DEBUG)
+endif()
+
+######################################################################
+## Enable/disable solver modules
+
+option(OPT_ENABLE_CPU_SOLVER "Enable CPU DG solver")
+if (OPT_ENABLE_CPU_SOLVER)
+    add_definitions(-DENABLE_CPU_SOLVER)
+endif()
+
+option(OPT_DEBUG_PRINT "Print debug information")
+if (OPT_DEBUG_PRINT)
+    add_definitions(-DENABLE_DEBUG_PRINT)
+endif()
+
+######################################################################
+## Configure targets 
+
+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)
+target_link_libraries(gmsh_gpu_dg ${LINK_LIBS})

cmake/FindSILO.cmake

+include(FindPackageHandleStandardArgs)
+
+
+find_path(SILO_INCLUDE_DIRS
+          NAMES silo.h
+          PATHS /usr/include /usr/local/include)
+
+find_library(SILO_LIBRARIES
+             NAMES silo siloh5
+             PATHS /usr/lib /usr/local/include)
+
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(SILO DEFAULT_MSG SILO_LIBRARIES SILO_INCLUDE_DIRS)

sgr @ af618523

+Subproject commit af618523a0c76d5ae4cebde9f5f07c702833d924

sol2 @ f56b3c69

+Subproject commit f56b3c698c2116f41aad3bf731ba8400e68c498b

eigen.h

+#pragma once
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/StdVector>
+
+template<typename T>
+using eigen_compatible_stdvector = std::vector<T, Eigen::aligned_allocator<T>>;
+

entity.cpp

+#include <iostream>
+
+#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)
+{
+    /* Get all the elements belonging to this group */
+    gmm::getElementsByType(g_elemType, g_elemTags, g_elemNodesTags, g_tag);
+    auto num_elems = g_elemTags.size();
+
+#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());
+
+    gmm::getJacobians(g_elemType, ips, jacs, dets, pts, g_tag);
+
+    assert(jacs.size() == 9*dets.size());
+    assert(pts.size() == 3*dets.size());
+    assert(dets.size() == num_elems*iws.size());
+
+    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);
+
+    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);
+
+    assert(g_num_bfun_orient == g_num_bgrad_orient);
+
+    gmm::getBarycenters(g_elemType, g_tag, false, false, barycenters);
+    assert(barycenters.size() == 3*num_elems);
+
+#ifdef ENABLE_DEBUG_PRINT
+    std::cout << yellowfg << "Orientations: " << g_num_bfun_orient << " " << g_num_bgrad_orient << reset << std::endl;
+#endif
+}
+

gmsh_io.cpp

+#include <iostream>
+#include <sstream>
+#include <cassert>
+
+#include "gmsh_io.h"
+
+std::string
+quadrature_name(int order)
+{
+    std::stringstream ss;
+    ss << "Gauss" << order;
+    return ss.str();
+}
+
+std::string
+basis_func_name(int order)
+{
+    std::stringstream ss;
+    ss << "H1Legendre" << order;
+    return ss.str();
+}
+
+std::string
+basis_grad_name(int order)
+{
+    std::stringstream ss;
+    ss << "GradH1Legendre" << order;
+    return ss.str();
+}
+
+model::model()
+    : geometric_order(1), approximation_order(1)
+{
+    read_gmsh_entities();
+}
+
+model::model(int pgo, int pao)
+    : geometric_order(pgo), approximation_order(pao)
+{
+    assert(geometric_order >= 1 and geometric_order <= 3);
+    assert(approximation_order > 0);
+
+    read_gmsh_entities();
+}
+
+model::model(const char *filename)
+    : geometric_order(1), approximation_order(1)
+{
+    gmsh::open( std::string(filename) );
+    read_gmsh_entities();
+}
+
+model::model(const char *filename, int pgo, int pao)
+    : geometric_order(pgo), approximation_order(pao)
+{
+    assert(geometric_order >= 1 and geometric_order <= 3);
+    assert(approximation_order > 0);
+
+    gmsh::open( std::string(filename) );
+    read_gmsh_entities();
+}
+
+void
+model::read_gmsh_entities(void)
+{
+    gvp_t entities;
+    gm::getEntities(entities, DIMENSION(3));
+
+    for (auto [dim, tag] : 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)
+        {
+            
+        }
+    }
+}

gmsh_io.h

+#pragma once
+
+#include <vector>
+
+#include <gmsh.h>
+
+#include "reference_element.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;
+
+/* Macros used to annotate GMSH integer inputs */
+#define DIMENSION(x) x
+
+std::string quadrature_name(int);
+std::string basis_func_name(int);
+std::string basis_grad_name(int);
+
+class model
+{
+    int                             geometric_order;
+    int                             approximation_order;
+
+    void    read_gmsh_entities(void);
+
+public:
+    model();
+    model(int, int);
+    model(const char *);
+    model(const char *, int, int);
+};
+

main.cpp

+#include <iostream>
+
+#include <unistd.h>
+
+#include "gmsh_io.h"
+
+static void usage(const char *progname)
+{
+    std::cout << "Usage: " << progname << "[params]" << std::endl;
+    std::cout << "<TODO>" << std::endl;
+}
+
+int main(int argc, char *argv[])
+{
+    const char *model_filename = nullptr;
+
+    int ch;
+    while ( (ch = getopt(argc, argv, "m:")) != -1 )
+    {
+        switch(ch)
+        {
+        case 'm':
+            model_filename = optarg;
+            break;
+
+        default:
+            usage(argv[0]);
+            return 1;
+        }
+    }
+
+    if (!model_filename)
+    {
+        std::cout << "Model filename (-m) not specified" << std::endl;
+        return 1;
+    }
+
+    gmsh::initialize(argc, argv);
+
+    
+
+    model mod;
+
+    gmsh::finalize();
+    return 0;
+}
+
+

physical_element.cpp

+
+#include "physical_element.h"
+#include "gmsh_io.h"
+
+physical_element::physical_element()
+{}
+
+physical_elements_factory::physical_elements_factory(int pdim, int ptag, int pElemType)
+    : dim(pdim), tag(ptag), elemType(pElemType)
+{
+    
+}
+
+std::vector<physical_element>
+physical_elements_factory::get_elements()
+{
+    auto order = 2;
+    std::vector<double>     ips;
+    std::vector<double>     iws;
+    gmm::getIntegrationPoints(elemType, quadrature_name(2*order), ips, iws);
+
+    auto num_gf = iws.size();
+
+    std::vector<double>     ppts;
+    std::vector<double>     jacs;
+    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<physical_element> ret;
+    ret.reserve(num_elems);
+
+    for (size_t elem = 0; elem < num_elems; elem++)
+    {
+        physical_element new_pe;
+        new_pe.m_original_position = elem;
+        new_pe.m_dimension = dim;
+        new_pe.m_parent_entity_tag = tag;
+
+        for (size_t gf = 0; gf < num_gf; gf++)
+        {
+            const auto JSIZE = 3*3; /* Jacobian matrix size */
+            auto jacs_base = elem*num_gf*JSIZE + gf*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];
+            new_pe.m_jacobians.push_back(jac);
+
+            auto dets_ofs = elem*num_gf + gf;
+            assert( dets_ofs < dets.size() );
+            new_pe.m_determinants.push_back(dets[dets_ofs]);
+
+            const auto PSIZE = 3;
+            auto pts_base = elem*num_gf*PSIZE + gf*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);
+        }
+    }
+
+    return std::vector<physical_element>();
+}
+

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>;
+
+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;
+
+public:
+    physical_element();
+
+    friend class physical_elements_factory;
+};
+
+class physical_elements_factory
+{
+    int         dim;
+    int         tag;
+    int         elemType;
+
+public:
+    physical_elements_factory(int, int, int);
+
+    std::vector<physical_element> get_elements();
+};
+

point.hpp

+/*
+ * Yaourt-FEM - Yet AnOther Useful Resource for Teaching FEM and DG.
+ *
+ * Matteo Cicuttin (C) 2019
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <iostream>
+#include <array>
+#include <cmath>
+
+template<typename T, size_t DIM>
+class point
+{
+    std::array<T, DIM>     m_coords;
+
+public:
+    typedef T               value_type;
+    const static size_t     dimension = DIM;
+
+    point()
+    {
+        for (size_t i = 0; i < DIM; i++)
+            m_coords[i] = T(0);
+    }
+
+    point(const point& other)
+        : m_coords(other.m_coords)
+    {}
+    
+    point operator=(const point& other)
+    {
+        m_coords = other.m_coords;
+        return *this;
+    }
+    
+    template<typename U = T>
+    point(const typename std::enable_if<DIM == 1, U>::type& x)
+    {
+        m_coords[0] = x;
+    }
+    
+    template<typename U = T>
+    point(const typename std::enable_if<DIM == 2, U>::type& x, const U& y)
+    {
+        m_coords[0] = x;
+        m_coords[1] = y;
+    }
+    
+    template<typename U = T>
+    point(const typename std::enable_if<DIM == 3, U>::type& x, const U& y, const U& z)
+    {
+        m_coords[0] = x;
+        m_coords[1] = y;
+        m_coords[2] = z;
+    }
+
+    T   at(size_t pos) const { return m_coords.at(pos); }
+    T&  at(size_t pos)       { return m_coords.at(pos); }
+
+    T   operator[](size_t pos) const { return m_coords[pos]; }
+    T&  operator[](size_t pos)       { return m_coords[pos]; }
+
+    point operator-() const {
+        auto ret = -1.0 * (*this);
+        return ret;
+    }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 1 || DIM == 2 || DIM == 3, U>::type
+    x() const { return m_coords[0]; }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 1 || DIM == 2 || DIM == 3, U>::type&
+    x() { return m_coords[0]; }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 2 || DIM == 3, U>::type
+    y() const { return m_coords[1]; }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 2 || DIM == 3, U>::type&
+    y() { return m_coords[1]; }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 3, U>::type
+    z() const { return m_coords[2]; }
+
+    template<typename U = T>
+    typename std::enable_if<DIM == 3, U>::type&
+    z() { return m_coords[2]; }
+
+    point&
+    operator+=(const point& other)
+    {
+        for (size_t i = 0; i < DIM; i++)
+            m_coords[i] += other.m_coords[i];
+
+        return *this;
+    }
+
+    point
+    operator+(const point& other) const
+    {
+        point ret = *this;
+        ret += other;
+        return ret;
+    }
+
+    point&
+    operator-=(const point& other)
+    {
+        for (size_t i = 0; i < DIM; i++)
+            m_coords[i] -= other.m_coords[i];
+
+        return *this;
+    }
+
+    point
+    operator-(const point& other) const
+    {
+        point ret = *this;
+        ret -= other;
+        return ret;
+    }
+
+    point&
+    operator*=(const T& scalefactor)
+    {
+        for (size_t i = 0; i < DIM; i++)
+            m_coords[i] *= scalefactor;
+
+        return *this;
+    }
+
+    point
+    operator*(const T& scalefactor) const
+    {
+        point ret = *this;
+        ret *= scalefactor;
+        return ret;
+    }
+
+    friend point
+    operator*(T scalefactor, const point& p)
+    {
+        return p * scalefactor;
+    }
+
+    point
+    operator/=(const T& scalefactor)
+    {
+        point ret;
+        for (size_t i = 0; i < DIM; i++)
+            m_coords[i] /= scalefactor;
+        
+        return ret;
+    }
+
+    point
+    operator/(const T& scalefactor) const
+    {
+        point ret = *this;
+        ret /= scalefactor;
+        return ret;
+    }
+};
+
+template<typename T, size_t DIM>
+T
+distance(const point<T,DIM>& p1, const point<T,DIM>& p2)
+{
+    auto acc = 0.0;
+
+    for (size_t i = 0; i < DIM; i++)
+        acc += (p1[i] - p2[i])*(p1[i] - p2[i]);
+
+    return std::sqrt(acc);
+}
+
+template<typename T, size_t DIM>
+std::ostream&
+operator<<(std::ostream& os, const point<T, DIM>& pt)
+{
+    os << "( ";
+    for (size_t i = 0; i < DIM; i++)
+    {
+        os << pt[i];
+
+        if (i < DIM-1)
+            os << ", ";
+    }
+    os << " )";
+    return os;
+}
+
+using point_3f = point<float,3>;
+using point_3d = point<double,3>;
+
+template<typename T, size_t DIM>
+class quadrature_point
+{
+    point<T, DIM>   p;
+    T               w;
+
+    quadrature_point()
+    {}
+
+    quadrature_point(const point<T,DIM>& pp, const T& pw)
+        : p(pp), w(pw)
+    {}
+
+    point<T,DIM> point() const
+    {
+        return p;
+    }
+
+    T weight() const
+    {
+        return w;
+    }
+};
+
+using quadrature_point_3f = quadrature_point<float, 3>;
+using quadrature_point_3d = quadrature_point<double, 3>;
+

reference_element.h

+#pragma once
+
+class reference_cell
+{};
+
+class reference_face
+{};
+