diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6c600cd624e1da76d5e4113e0a6456265c125fd9..1abd4e1cd6b8ebe181c8bed224f602fe76a945d4 100644
--- a/CMakeLists.txt
+++ b/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})
diff --git a/dump_orient.cpp b/dump_orient.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c941697d58f9e5177245f2c4b805986279d2fe79
--- /dev/null
+++ b/dump_orient.cpp
@@ -0,0 +1,123 @@
+#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;
+}
diff --git a/entity.cpp b/entity.cpp
index 86b2aae57065a2d8f540c87ffc89918ee519f71a..e8cad719a1ae7b246350f72e90a8c764957a6178 100644
--- a/entity.cpp
+++ b/entity.cpp
@@ -1,49 +1,362 @@
#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);
+}
+
diff --git a/entity.h b/entity.h
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..71d3ab964bf786846ca851f0d1a27c6fc07560ec 100644
--- a/entity.h
+++ b/entity.h
@@ -0,0 +1,75 @@
+#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);
+
diff --git a/entity_data.h b/entity_data.h
new file mode 100644
index 0000000000000000000000000000000000000000..130ce394efb3290982793e9f9191e3883f725aa0
--- /dev/null
+++ b/entity_data.h
@@ -0,0 +1,23 @@
+#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;
+};
+
diff --git a/gmsh_io.cpp b/gmsh_io.cpp
index ef5af96157df81584b0b3304b0adf9c280398980..3b4defe325d301b297bdfed5e7f9217522170884 100644
--- a/gmsh_io.cpp
+++ b/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) );
}
}
}
diff --git a/gmsh_io.h b/gmsh_io.h
index 4e542960e8e64d894df152d5d8fdbfca242796e0..5c18fd0cbd2e49ec4b50b350929bfe1852071a04 100644
--- a/gmsh_io.h
+++ b/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;
+
};
diff --git a/kernels_cpu.cpp b/kernels_cpu.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..11dc87d0b96a526521213c567177dfc9199f504f
--- /dev/null
+++ b/kernels_cpu.cpp
@@ -0,0 +1,46 @@
+#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
diff --git a/kernels_cpu.h b/kernels_cpu.h
new file mode 100644
index 0000000000000000000000000000000000000000..08b5abd085b5ec60a12e318892543f1581902d9e
--- /dev/null
+++ b/kernels_cpu.h
@@ -0,0 +1,5 @@
+#pragma once
+
+void compute_field_derivatives(const entity_data_cpu<1>& ed,
+ const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz);
+
diff --git a/main.cpp b/main.cpp
index ac49ce9f3ed2c3b0b3ff7f1078003de64ba67f1c..ed7832b1dbcbcf633eb09ed96013632993bdaaec 100644
--- a/main.cpp
+++ b/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;
diff --git a/physical_element.cpp b/physical_element.cpp
index c4ad9c0ef181c8559c6139ab0f63415ef7eab6d3..a8b34ea87482438d440eac4db615962ebc9fe795 100644
--- a/physical_element.cpp
+++ b/physical_element.cpp
@@ -1,23 +1,113 @@
+#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;
}
-
diff --git a/physical_element.h b/physical_element.h
index 8a977a2aedd9c662a3b904ede63eadf9b94fd318..867bec1045bd7e420c8789feb85df0bb244344ab 100644
--- a/physical_element.h
+++ b/physical_element.h
@@ -1,29 +1,46 @@
#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();
};
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/point.hpp b/point.hpp
index 5a53b129536e80fe7713c46f2c52d8678a97f774..03662cbd4d6dff47db40cb5298e4a52727c189ec 100644
--- a/point.hpp
+++ b/point.hpp
@@ -217,6 +217,7 @@ class quadrature_point
point<T, DIM> p;
T w;
+public:
quadrature_point()
{}
diff --git a/reference_element.cpp b/reference_element.cpp
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..b91ae6aeed23b6cd6833420544a1b247fae211d1 100644
--- a/reference_element.cpp
+++ b/reference_element.cpp
@@ -0,0 +1,198 @@
+#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;
+}
+
diff --git a/reference_element.h b/reference_element.h
index 990a29db63914cd9ce5358733cf5a874563feca0..446b94c499efde0e52f945d8c9842491756e30c6 100644
--- a/reference_element.h
+++ b/reference_element.h
@@ -1,8 +1,70 @@
#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();
+};
diff --git a/silo_output.hpp b/silo_output.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..daac379a41959cb29f7a51c59c5ed85d081153c9
--- /dev/null
+++ b/silo_output.hpp
@@ -0,0 +1,182 @@
+#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;
+ }
+};
diff --git a/test.cpp b/test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b162a87fc8beeead3061fb39d5df004363a74f3
--- /dev/null
+++ b/test.cpp
@@ -0,0 +1,345 @@
+#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;
+}
+
+
diff --git a/test.h b/test.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a9dd89706de43d661379dba383df2b931b6154f
--- /dev/null
+++ b/test.h
@@ -0,0 +1,24 @@
+#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; }
+
diff --git a/types.h b/types.h
new file mode 100644
index 0000000000000000000000000000000000000000..64116d921fca3de6269b269b3d125f4cffc05f21
--- /dev/null
+++ b/types.h
@@ -0,0 +1,24 @@
+#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;
+};
+