diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2cc0e1268f32904b68135746f6dab20b53ed06e6..bb86e28beb514b219c04b40a6678a36c48696ef7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -233,10 +233,12 @@ endif()
include_directories(contrib/sgr)
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+set(COMMON_SOURCES gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp connectivity.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 test_basics.cpp test_mass.cpp test_differentiation.cpp gmsh_io.cpp reference_element.cpp physical_element.cpp entity.cpp kernels_cpu.cpp)
+add_executable(test test.cpp test_basics.cpp test_mass.cpp test_differentiation.cpp test_lifting.cpp ${COMMON_SOURCES})
target_link_libraries(test ${LINK_LIBS})
add_executable(dump_orient dump_orient.cpp)
diff --git a/connectivity.cpp b/connectivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b42f125fb4b973026fd5c77cdbeea410e93f63af
--- /dev/null
+++ b/connectivity.cpp
@@ -0,0 +1,53 @@
+#include <iostream>
+
+#include "connectivity.h"
+
+/***************************************************************************/
+neighbour_descriptor::neighbour_descriptor()
+{}
+
+neighbour_descriptor::neighbour_descriptor(size_t p_iE, size_t p_iT, size_t p_via_iF)
+ : iE(p_iE), iT(p_iT), via_iF(p_via_iF)
+{}
+
+bool
+neighbour_descriptor::operator<(const neighbour_descriptor& other) const
+{
+ std::array<size_t, 3> arr_me{iE, iT, via_iF};
+ std::array<size_t, 3> arr_other{other.iE, other.iT, other.via_iF};
+ return arr_me < arr_other;
+}
+
+std::ostream& operator<<(std::ostream& os, const neighbour_descriptor& nd)
+{
+ os << "(" << nd.iE << ", " << nd.iT << ", " << nd.via_iF << ")";
+ return os;
+}
+
+/***************************************************************************/
+cell_descriptor::cell_descriptor()
+{}
+
+cell_descriptor::cell_descriptor(size_t p_iE, size_t p_iT)
+ : iE(p_iE), iT(p_iT)
+{}
+
+bool
+cell_descriptor::operator<(const cell_descriptor& other) const
+{
+ if (iE < other.iE)
+ return true;
+
+ if (iE == other.iE)
+ return (iT < other.iT);
+
+ return false;
+}
+
+std::ostream& operator<<(std::ostream& os, const cell_descriptor& cd)
+{
+ os << "(" << cd.iE << "," << cd.iT << ")";
+ return os;
+}
+
+/***************************************************************************/
diff --git a/connectivity.h b/connectivity.h
new file mode 100644
index 0000000000000000000000000000000000000000..6df1b85b43c9d5632a9635a80babc984e2edbad6
--- /dev/null
+++ b/connectivity.h
@@ -0,0 +1,227 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <array>
+#include <map>
+#include <set>
+
+#include "sgr.hpp"
+
+/***************************************************************************/
+struct neighbour_descriptor
+{
+ size_t iE; /* Entity offset in data structures */
+ size_t iT; /* Element offset in entity in data structures */
+ size_t via_iF; /* Element face index (i.e. 0,1,2 for triangles) */
+
+ neighbour_descriptor();
+ neighbour_descriptor(size_t, size_t, size_t);
+ bool operator<(const neighbour_descriptor&) const;
+};
+
+std::ostream& operator<<(std::ostream&, const neighbour_descriptor&);
+
+/***************************************************************************/
+struct cell_descriptor
+{
+ size_t iE;
+ size_t iT;
+
+ cell_descriptor();
+ cell_descriptor(size_t, size_t);
+ bool operator<(const cell_descriptor&) const;
+};
+
+std::ostream& operator<<(std::ostream&, const cell_descriptor&);
+
+/***************************************************************************/
+template<typename NT>
+class neighbour_pair
+{
+ NT n_first;
+ NT n_second;
+ size_t n_used;
+
+public:
+ neighbour_pair()
+ : n_used(0)
+ {}
+
+ void
+ add_neighbour(const NT& n)
+ {
+ assert(n_used == 0 or n_used == 1);
+ if (n_used == 0)
+ {
+ n_first = n;
+ n_used = 1;
+ return;
+ }
+
+ if (n_used == 1)
+ {
+ n_second = n;
+ n_used = 2;
+ return;
+ }
+
+ throw std::out_of_range("Pair has already two neighbours");
+ }
+
+ NT first() const
+ {
+ assert(n_used == 1 or n_used == 2);
+ if (n_used < 1)
+ throw std::out_of_range("Pair does not have n_first");
+
+ return n_first;
+ }
+
+ NT second() const
+ {
+ assert(n_used == 2);
+ if (n_used < 2)
+ throw std::out_of_range("Pair does not have n_second");
+
+ return n_second;
+ }
+
+ size_t num_neighbours() const
+ {
+ assert(n_used < 3);
+ return n_used;
+ }
+};
+
+/***************************************************************************/
+template<typename FaceKey>
+class element_connectivity
+{
+public:
+ using ConnData = std::pair<FaceKey, neighbour_descriptor>;
+ using NeighPair = neighbour_pair<neighbour_descriptor>;
+
+private:
+ std::map<FaceKey, NeighPair> face_cell_adj;
+ std::map<cell_descriptor, std::set<ConnData>> cell_cell_adj;
+
+public:
+ element_connectivity()
+ {}
+
+ void
+ connect(const FaceKey& fk, const neighbour_descriptor& nd)
+ {
+ /* Insert the information that one the two neighbours around
+ * face `e' is `nd'. */
+ auto& nr = face_cell_adj[fk];
+ nr.add_neighbour(nd);
+
+ /* If we've got both neighbours around `e`, add a connection
+ * between them in the cell-cell adjacency list */
+ if ( nr.num_neighbours() == 2 )
+ {
+ auto fst = nr.first();
+ auto snd = nr.second();
+ cell_cell_adj[{fst.iE, fst.iT}].insert({fk, snd});
+ cell_cell_adj[{snd.iE, snd.iT}].insert({fk, fst});
+ }
+ }
+
+ /* Find the neighbour of the element (iE,iT) via face fk */
+ std::pair<neighbour_descriptor, bool>
+ neighbour_via(size_t iE, size_t iT, const FaceKey& fk) const
+ {
+ auto nr = face_cell_adj.at(fk);
+
+ if ( nr.num_neighbours() == 1 )
+ {
+#ifndef _NDEBUG
+ auto fst = nr.first();
+ assert(fst.iE == iE);
+ assert(fst.iT == iT);
+#endif
+ return std::make_pair(neighbour_descriptor(), false);
+ }
+
+ if ( nr.num_neighbours() == 2 )
+ {
+ auto fst = nr.first();
+ auto snd = nr.second();
+
+ assert( (fst.iE == iE and fst.iT == iT) or
+ (snd.iE == iE and snd.iT == iT)
+ );
+
+ if ( fst.iE == iE and fst.iT == iT)
+ return std::make_pair(snd, true);
+
+ return std::make_pair(fst, true);
+ }
+
+ assert(true && "Shouldn't have arrived here");
+ return std::make_pair(neighbour_descriptor(), false); /* Silence the compiler */
+ }
+
+ std::set<ConnData>
+ neighbours(size_t iE, size_t iT) const
+ {
+ return cell_cell_adj.at({iE, iT});
+ }
+
+ std::set<ConnData>
+ neighbours(const cell_descriptor& cd) const
+ {
+ return cell_cell_adj.at(cd);
+ }
+
+ std::set<cell_descriptor>
+ neighbours_as_cell_descriptors(const cell_descriptor& cd) const
+ {
+ auto filter = [](const ConnData& connd) -> auto {
+ cell_descriptor ret;
+ ret.iE = connd.second.iE;
+ ret.iT = connd.second.iT;
+ return ret;
+ };
+
+ auto neighs = cell_cell_adj.at(cd);
+ std::set<cell_descriptor> ret;
+ std::transform(neighs.begin(), neighs.end(),
+ std::inserter(ret, ret.begin()), filter);
+ return ret;
+ }
+
+ bool
+ are_connected(const cell_descriptor& a, const cell_descriptor& b) const
+ {
+ auto neighs = neighbours(a);
+ for (const auto& [fk, nd] : neighs)
+ if (nd.iE == b.iE and nd.iT == b.iT)
+ return true;
+
+ return false;
+ }
+
+ std::set<cell_descriptor>
+ element_keys() const
+ {
+ using PT = std::pair<cell_descriptor, std::set<ConnData>>;
+ auto get_node = [](const PT& adjpair) -> auto {
+ return adjpair.first;
+ };
+
+ auto cc_begin = cell_cell_adj.begin();
+ auto cc_end = cell_cell_adj.end();
+ std::set<cell_descriptor> ret;
+ std::transform(cc_begin, cc_end, std::inserter(ret, ret.begin()), get_node);
+ return ret;
+ }
+
+ void clear(void)
+ {
+ face_cell_adj.clear();
+ }
+};
+
diff --git a/entity.cpp b/entity.cpp
index 32b509b847aea858b7ff0457db3eb070fd51ebfe..15154488c1dd44f6a67599ba96e4348c85ecf610 100644
--- a/entity.cpp
+++ b/entity.cpp
@@ -138,14 +138,14 @@ entity::num_cells_per_orientation(void) const
}
const physical_element&
-entity::cell(size_t iT)
+entity::cell(size_t iT) const
{
assert(iT < physical_cells.size());
return physical_cells[iT];
}
const physical_element&
-entity::face(size_t iF)
+entity::face(size_t iF) const
{
assert(iF < physical_faces.size());
return physical_faces[iF];
@@ -153,7 +153,7 @@ entity::face(size_t iF)
matxd
-entity::mass_matrix(size_t iT)
+entity::mass_matrix(size_t iT) const
{
const auto& pe = physical_cells[iT];
assert( pe.orientation() < reference_cells.size() );
@@ -177,6 +177,9 @@ entity::mass_matrix(size_t iT)
void
entity::sort_by_orientation(void)
{
+ if (cur_elem_ordering == entity_ordering::BY_ORIENTATION)
+ return;
+
/* 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();
@@ -191,6 +194,40 @@ entity::sort_by_orientation(void)
};
std::sort(physical_cells.begin(), physical_cells.end(), comp);
+
+ std::vector<physical_element> new_pf;
+ new_pf.resize( physical_faces.size() );
+ std::vector<size_t> new_ft;
+ new_ft.resize( faceTags.size() );
+ std::vector<size_t> new_fnt;
+ new_fnt.resize( faceNodesTags.size() );
+
+ size_t nodes_per_face = faceNodesTags.size()/faceTags.size();
+ assert(nodes_per_face * faceTags.size() == faceNodesTags.size());
+
+ for (size_t new_iT = 0; new_iT < physical_cells.size(); new_iT++)
+ {
+ auto old_iT = physical_cells[new_iT].original_position();
+ for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
+ {
+ auto old_ofs = NUM_TET_FACES*old_iT + iF;
+ auto new_ofs = NUM_TET_FACES*new_iT + iF;
+ new_pf[new_ofs] = physical_faces[old_ofs];
+ new_ft[new_ofs] = faceTags[old_ofs];
+
+ for (size_t iN = 0; iN < nodes_per_face; iN++)
+ {
+ auto old_ofs_nt = nodes_per_face*old_ofs + iN;
+ auto new_ofs_nt = nodes_per_face*new_ofs + iN;
+ new_fnt[new_ofs_nt] = faceNodesTags[old_ofs_nt];
+ }
+ }
+ }
+
+ std::swap(new_pf, physical_faces);
+ std::swap(new_ft, faceTags);
+ std::swap(new_fnt, faceNodesTags);
+
cur_elem_ordering = entity_ordering::BY_ORIENTATION;
}
@@ -210,6 +247,7 @@ entity::sort_by_gmsh(void)
std::sort(physical_cells.begin(), physical_cells.end(), comp);
cur_elem_ordering = entity_ordering::GMSH;
+ throw std::logic_error("NOT SORTED");
}
int
@@ -343,8 +381,13 @@ entity::populate_lifting_matrices_planar(entity_data_cpu& ed,
for (size_t iO = 0; iO < ed.num_orientations; iO++)
{
+ if (dofs_f2c.find(iO) == dofs_f2c.end())
+ continue;
+
auto& f2c = dofs_f2c.at(iO);
assert(f2c.size() == NUM_TET_FACES*ed.num_fluxes);
+ matxd lm = matxd::Zero(ed.num_bf, lm_cols);
+
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
auto fO = FACE_ORIENTATION(iO, iF);
@@ -358,22 +401,33 @@ entity::populate_lifting_matrices_planar(entity_data_cpu& ed,
for (size_t i = 0; i < mm.rows(); i++)
{
assert(iF*ed.num_fluxes + i < f2c.size());
- auto lm_row = iO*ed.num_bf + f2c[iF*ed.num_fluxes + i];
- assert(lm_row < lm_rows);
+ auto lm_row = f2c[iF*ed.num_fluxes + i];
+ assert(lm_row < ed.num_bf);
for (size_t j = 0; j < mm.cols(); j++)
{
auto lm_col = iF*ed.num_fluxes + j;
assert(lm_col < lm_cols);
- ed.lifting_matrices(lm_row, lm_col) = mm(i, j);
+ lm(lm_row, lm_col) = mm(i, j);
}
}
}
+
+ auto& re = reference_cells[iO];
+ 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.lifting_matrices.block(base_row, 0, ed.num_bf, lm_cols) =
+ mass_ldlt.solve(lm);
}
}
void
entity::populate_lifting_matrices_curved(entity_data_cpu&, const dofs_f2c_t&) const
-{}
+{
+}
void
entity::populate_lifting_matrices(entity_data_cpu& ed) const
@@ -387,26 +441,38 @@ entity::populate_lifting_matrices(entity_data_cpu& ed) const
{
auto& pe = physical_cells[iT];
auto orient = pe.orientation();
- auto cnt = pe.node_tags();
+ auto cnt = pe.bf_keys();
+#ifndef EXPENSIVE_ASSERTS
if ( dofs_f2c.find(orient) != dofs_f2c.end() )
continue;
+#endif
- std::vector<size_t> f2c( NUM_TET_FACES*num_face_bf );
+ std::vector<size_t> f2c;
+ f2c.reserve( NUM_TET_FACES*num_face_bf );
for (size_t iF = 0; iF < NUM_TET_FACES; iF++)
{
- auto origpos = pe.original_position();
- auto face_ofs = NUM_TET_FACES*origpos+iF;
+ auto face_ofs = NUM_TET_FACES*iT+iF;
auto& pf = physical_faces[face_ofs];
- auto fnt = pf.node_tags();
+ auto fnt = pf.bf_keys();
for (size_t i = 0; i < fnt.size(); i++)
for (size_t j = 0; j < cnt.size(); j++)
- if (cnt[i] == fnt[j])
- f2c[iF*num_face_bf + i] = j;
+ if (cnt[j] == fnt[i])
+ {
+ assert(f2c.size() == iF*num_face_bf + i);
+ f2c.push_back(j);
+ }
- dofs_f2c[orient] = f2c;
}
+ assert(f2c.size() == NUM_TET_FACES*num_face_bf);
+
+#ifdef EXPENSIVE_ASSERTS
+ if ( dofs_f2c.find(orient) != dofs_f2c.end() )
+ assert(dofs_f2c[orient] == f2c);
+#endif
+
+ dofs_f2c[orient] = f2c;
}
if (g_order == 1)
@@ -423,10 +489,20 @@ entity::populate_jacobians(entity_data_cpu& ed) const
if (g_order == 1)
{
ed.jacobians = matxd::Zero(3*num_elems, 3);
+ ed.cell_determinants = vecxd::Zero(num_elems);
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();
+ ed.cell_determinants(iT) = pe.determinant();
+ }
+
+ assert(num_faces() == 4*num_elems);
+ ed.face_determinants = vecxd(4*num_elems);
+ for (size_t iF = 0; iF < num_faces(); iF++)
+ {
+ auto& pf = physical_faces[iF];
+ ed.face_determinants(iF) = pf.determinant();
}
}
else
@@ -483,8 +559,10 @@ entity::populate_normals(entity_data_cpu& ed) const
void
entity::populate_entity_data(entity_data_cpu& ed) const
{
+ assert(cur_elem_ordering == entity_ordering::BY_ORIENTATION);
assert(reference_faces.size() > 0);
assert(reference_cells.size() > 0);
+ ed.ordering = cur_elem_ordering;
ed.g_order = g_order;
ed.a_order = a_order;
ed.num_elems = num_cells_per_orientation();
@@ -493,7 +571,8 @@ entity::populate_entity_data(entity_data_cpu& ed) const
ed.num_fluxes = reference_faces[0].num_basis_functions();
ed.num_qp = (g_order == 1) ? 1 : reference_cells[0].num_integration_points();
ed.num_face_qp = (g_order == 1) ? 1 : reference_faces[0].num_integration_points();
- ed.dofs_base = 0;
+ ed.dof_base = 0;
+ ed.flux_base = 0;
populate_differentiation_matrices(ed);
populate_jacobians(ed);
@@ -503,8 +582,9 @@ entity::populate_entity_data(entity_data_cpu& ed) const
ed.invmass_matrices = matxd::Zero(ed.num_bf*num_cells(), ed.num_bf);
for (size_t iT = 0; iT < num_cells(); iT++)
{
+ matxd mass = mass_matrix(iT);
ed.invmass_matrices.block(iT*ed.num_bf, 0, ed.num_bf, ed.num_bf) =
- mass_matrix(iT).inverse();
+ mass.inverse();
}
}
diff --git a/entity.h b/entity.h
index 7e23a26a7e6c453eb76c436d22cd577b4a11f936..fc9219b97b32f734339db41f9ddabe8e1f3ca5c3 100644
--- a/entity.h
+++ b/entity.h
@@ -4,16 +4,11 @@
#include "reference_element.h"
#include "physical_element.h"
-
+#include "connectivity.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;
@@ -29,10 +24,10 @@ class entity
std::vector<physical_element> physical_cells;
std::vector<reference_element> reference_faces;
- std::vector<physical_element> physical_faces; /* Always ordered by GMSH */
+ std::vector<physical_element> physical_faces;
- std::vector<size_t> faceTags; /* Always ordered by GMSH */
- std::vector<size_t> faceNodesTags; /* Always ordered by GMSH */
+ std::vector<size_t> faceTags;
+ std::vector<size_t> faceNodesTags;
using dofs_f2c_t = std::map<size_t, std::vector<size_t>>;
void populate_lifting_matrices_planar(entity_data_cpu&, const dofs_f2c_t&) const;
@@ -66,19 +61,19 @@ public:
size_t num_cell_orientations(void) const;
std::vector<size_t> num_cells_per_orientation(void) const;
+ const physical_element& cell(size_t) const;
+ const reference_element& cell_refelem(size_t) const;
+ const reference_element& cell_refelem(const physical_element&) const;
+
size_t num_faces(void) const;
std::vector<size_t> face_tags() const;
- const reference_element& cell_refelem(const physical_element&) const;
+ const physical_element& face(size_t) const;
+ const reference_element& face_refelem(size_t) const;
const reference_element& face_refelem(const physical_element&) const;
- const physical_element& cell(size_t);
- const physical_element& cell_bytag(size_t);
- const physical_element& face(size_t);
- const physical_element& face_bytag(size_t);
-
- matxd mass_matrix(size_t);
+ matxd mass_matrix(size_t) const;
void sort_by_orientation(void);
void sort_by_gmsh(void);
@@ -89,6 +84,8 @@ public:
void populate_entity_data(entity_data_cpu&) const;
+ constexpr size_t num_faces_per_elem(void) const { return 4; }
+
std::vector<physical_element>::const_iterator begin() const { return physical_cells.begin(); }
std::vector<physical_element>::const_iterator end() const { return physical_cells.end(); }
@@ -97,4 +94,3 @@ public:
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
index c447bdf0d216870835e76cfded32d35c9e8fb04c..d0cec3424dab1207bd0d3f1c7fcdcaf6cc9e24cd 100644
--- a/entity_data.h
+++ b/entity_data.h
@@ -2,6 +2,11 @@
#include "eigen.h"
+enum class entity_ordering {
+ GMSH,
+ BY_ORIENTATION
+};
+
struct entity_data_cpu
{
int g_order;
@@ -13,7 +18,12 @@ struct entity_data_cpu
size_t num_fluxes; /* No. of flux dofs per elem */
size_t num_face_qp; /* No. of integration pts per face */
size_t num_faces; /* No. of faces */
- size_t dofs_base;
+ entity_ordering ordering;
+
+ size_t dof_base; /* Offset in the global dof vector */
+ size_t flux_base; /* Offset in the global flux vector */
+ std::vector<size_t> fluxdofs_mine;
+ std::vector<size_t> fluxdofs_neighbour;
/* Inverse of mass matrices. Populated only if geometric order > 1.
* num_bf*num_elems rows, num_bf columns. Matrix at offset elem*num_bf
@@ -39,4 +49,5 @@ struct entity_data_cpu
matxd normals;
vecxd cell_determinants;
+ vecxd face_determinants;
};
diff --git a/gmsh_io.cpp b/gmsh_io.cpp
index dbde114eeedeb164b7463c44f9438e759d440d89..6b13afd9f291dbd5e229b8086311e7b564ae7a3d 100644
--- a/gmsh_io.cpp
+++ b/gmsh_io.cpp
@@ -117,12 +117,28 @@ model::entity_at(size_t pos)
return entities.at(pos);
}
+void
+model::update_connectivity(const entity& e, size_t entity_index)
+{
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < e.num_faces_per_elem(); iF++)
+ {
+ size_t findex = e.num_faces_per_elem() * iT + iF;
+ neighbour_descriptor cd(entity_index, iT, iF);
+ element_key fk( e.face(findex) );
+ conn.connect( fk, cd );
+ }
+ }
+}
+
void
model::import_gmsh_entities(void)
{
gvp_t gmsh_entities;
gm::getEntities(gmsh_entities, DIMENSION(3));
+ size_t entity_index = 0;
for (auto [dim, tag] : gmsh_entities)
{
std::vector<int> elemTypes;
@@ -136,7 +152,12 @@ model::import_gmsh_entities(void)
{
entity e({.dim = dim, .tag = tag, .etype = elemType,
.gorder = geometric_order, .aorder = approximation_order});
+
+ e.sort_by_orientation();
+ update_connectivity(e, entity_index);
+
entities.push_back( std::move(e) );
+ entity_index++;
}
}
}
diff --git a/gmsh_io.h b/gmsh_io.h
index 6bb94845742d9a4ffeb2e7ebbc7a788d28974cc0..190515764b364e82bae4e35a25d73e428d82af91 100644
--- a/gmsh_io.h
+++ b/gmsh_io.h
@@ -5,6 +5,7 @@
#include <gmsh.h>
#include "entity.h"
+#include "physical_element.h"
namespace g = gmsh;
namespace go = gmsh::option;
@@ -32,11 +33,15 @@ class model
std::vector<entity> entities;
- std::map<size_t, std::vector<element_key>> boundary_map;
+ std::map<size_t, std::vector<element_key>> boundary_map;
+
+ element_connectivity<element_key> conn;
void map_boundaries(void);
void import_gmsh_entities(void);
+ void update_connectivity(const entity&, size_t);
+
public:
model();
model(int, int);
diff --git a/kernels_cpu.cpp b/kernels_cpu.cpp
index 0dc5b6ee9e2157eb17073f260877f26d1f815ab8..71becb6da678ff70168af3f30e07f4b43a020bba 100644
--- a/kernels_cpu.cpp
+++ b/kernels_cpu.cpp
@@ -53,4 +53,57 @@ void compute_field_derivatives(const entity_data_cpu& ed,
default:
std::cout << "compute_field_derivatives: invalid order" << std::endl;
}
+}
+
+/* Field derivatives kernel, case for elements with planar faces */
+void compute_flux_lifting(const entity_data_cpu& ed,
+ const vecxd& flux, vecxd& field)
+{
+ switch (ed.a_order)
+ {
+ case 1:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<1>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<1>(ed, flux, field);
+ break;
+
+ case 2:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<2>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<2>(ed, flux, field);
+ break;
+
+ case 3:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<3>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<3>(ed, flux, field);
+ break;
+
+ case 4:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<4>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<4>(ed, flux, field);
+ break;
+
+ case 5:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<5>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<5>(ed, flux, field);
+ break;
+
+ case 6:
+ if (ed.g_order == 1)
+ compute_lifting_kernel_planar<6>(ed, flux, field);
+ else
+ compute_lifting_kernel_curved<6>(ed, flux, field);
+ break;
+
+ default:
+ std::cout << "compute_flux_lifting: invalid order" << std::endl;
+ }
}
\ No newline at end of file
diff --git a/kernels_cpu.h b/kernels_cpu.h
index b3d33a481c229097932dbc3bd5ebb94c5ed366f0..b71bd44a1aded5bf83254e9b84c73d89af6a9325 100644
--- a/kernels_cpu.h
+++ b/kernels_cpu.h
@@ -45,7 +45,7 @@ void compute_field_derivatives_kernel_planar(const entity_data_cpu& ed,
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;
+ size_t dofofs = ed.dof_base + orient_base + iT*ed.num_bf;
for (size_t odof = 0; odof < KS::num_bf; odof++)
{
@@ -111,7 +111,7 @@ void compute_field_derivatives_kernel_curved(const entity_data_cpu& ed,
{
//printf( "Thread %d works with idx %lu\n", omp_get_thread_num(), iT);
size_t ent_iT = orient_elem_base + iT;
- size_t dofofs = ed.dofs_base + orient_base + iT*ed.num_bf;
+ size_t dofofs = ed.dof_base + orient_base + iT*ed.num_bf;
df_dx.segment(dofofs, ed.num_bf) = vecxd::Zero(ed.num_bf);
df_dy.segment(dofofs, ed.num_bf) = vecxd::Zero(ed.num_bf);
@@ -173,4 +173,56 @@ void compute_field_derivatives_kernel_curved(const entity_data_cpu& ed,
void compute_field_derivatives(const entity_data_cpu& ed,
- const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz);
\ No newline at end of file
+ const vecxd& f, vecxd& df_dx, vecxd& df_dy, vecxd& df_dz);
+
+
+template<size_t AK>
+void compute_lifting_kernel_planar(const entity_data_cpu& ed,
+ const vecxd& flux, vecxd& field)
+{
+ /* Flop count: 21*ed.num_bf*ed.num_bf*num_total_elems */
+ /* The elements must be ordered by orientation -> entity::sort_by_orientation() */
+
+ assert(ed.g_order == 1);
+
+ using KS = kernel_cpu_sizes<AK>;
+ assert(ed.num_bf == KS::num_bf);
+ assert(ed.num_fluxes == KS::num_fluxes);
+
+ 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 dof_base = ed.dof_base + orient_elem_base * ed.num_bf;
+ size_t flux_base = ed.flux_base + orient_elem_base * 4*ed.num_fluxes;
+
+ for (size_t iT = 0; iT < num_elems; iT++)
+ {
+ size_t dof_ofs = dof_base + iT*ed.num_bf;
+ size_t flux_ofs = flux_base + 4*iT*KS::num_fluxes;
+ auto inv_det = 1./ed.cell_determinants[orient_elem_base+iT];
+
+ for (size_t i = 0; i < KS::num_bf; i++)
+ {
+ auto facc = field(dof_ofs+i);
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ auto face_det = ed.face_determinants[4*(orient_elem_base+iT)+iF];
+ for (size_t k = 0; k < KS::num_fluxes; k++)
+ facc += face_det * inv_det * ed.lifting_matrices(iO*KS::num_bf+i, KS::num_fluxes*iF + k) * flux(flux_ofs+KS::num_fluxes*iF + k);
+ }
+ field(dof_ofs+i) = facc;
+ }
+ }
+
+ orient_elem_base += num_elems;
+ }
+}
+
+template<size_t AK>
+void compute_lifting_kernel_curved(const entity_data_cpu& ed,
+ const vecxd& flux, vecxd& field)
+{}
+
+void compute_flux_lifting(const entity_data_cpu& ed,
+ const vecxd& flux, vecxd& field);
\ No newline at end of file
diff --git a/physical_element.cpp b/physical_element.cpp
index d806625292a9b953bc88421c3f39b030b89beb4f..651a94bd97bcb8b6a1e55586948a8de57e77026b 100644
--- a/physical_element.cpp
+++ b/physical_element.cpp
@@ -127,6 +127,13 @@ physical_elements_factory::physical_elements_factory(const entity_params& ep)
gmm::getBasisFunctionsOrientationForElements(elemType, basis_func_name(approx_order),
cellOrientations, tag);
assert(cellOrientations.size() == cellTags.size());
+
+ std::vector<double> coord;
+ gmm::getKeysForElements(elemType, basis_func_name(approx_order),
+ keypairs, coord, tag, false);
+
+ keys_per_elem = gmm::getNumberOfKeysForElements(elemType, basis_func_name(approx_order));
+ assert(keys_per_elem*cellTags.size() == keypairs.size());
}
std::vector<physical_element>
@@ -155,6 +162,13 @@ physical_elements_factory::get_elements()
new_pe.m_measure = 0.0;
new_pe.m_gmsh_elemtype = elemType;
+ new_pe.m_bf_keys.resize(keys_per_elem);
+ for (size_t i = 0; i < keys_per_elem; i++)
+ {
+ auto [vi, vu] = keypairs[keys_per_elem*elem + i];
+ new_pe.m_bf_keys[i] = bf_key(vi,vu);
+ }
+
new_pe.m_node_tags.resize(num_nodes);
for (size_t i = 0; i < num_nodes; i++)
new_pe.m_node_tags[i] = cellNodesTags[num_nodes*elem+i];
diff --git a/physical_element.h b/physical_element.h
index 819fa91d69c68845f6a3f102a63c044a6589f1b4..1ee451e02a53da5299248817eaa16331f08878be 100644
--- a/physical_element.h
+++ b/physical_element.h
@@ -1,8 +1,33 @@
#pragma once
+#include "gmsh.h"
#include "point.hpp"
#include "types.h"
+/* This is a key for a basis function. Did this
+ * class to protect the code from the GMSH api,
+ * which is apparently changing */
+class bf_key
+{
+ int wtf_int; // No idea of what this integer represents
+ size_t wtf_uint; // This is some kind of bf numbering
+
+public:
+ bf_key() = default;
+ bf_key(const bf_key&) = default;
+ bf_key(bf_key&&) = default;
+ bf_key(int pi, size_t pu)
+ : wtf_int(pi), wtf_uint(pu)
+ {}
+
+ bf_key& operator=(const bf_key&) = default;
+
+ bool operator==(const bf_key& other) const
+ {
+ return (wtf_int == other.wtf_int) and (wtf_uint == other.wtf_uint);
+ }
+};
+
class physical_element
{
size_t m_original_position; /* Position in GMSH ordering (relative to entity) */
@@ -19,6 +44,7 @@ class physical_element
double m_measure; /* Element measure */
int m_geometric_order;
int m_approximation_order;
+ std::vector<bf_key> m_bf_keys;
public:
physical_element();
@@ -31,6 +57,8 @@ public:
int element_type(void) const;
vec_size_t node_tags(void) const;
+ const std::vector<bf_key> bf_keys(void) const { return m_bf_keys; }
+
double determinant(void) const;
double determinant(size_t) const;
vec_double determinants(void) const;
@@ -64,6 +92,8 @@ class physical_elements_factory
std::vector<double> cellJacs;
std::vector<double> cellDets;
std::vector<int> cellOrientations;
+ gmsh::vectorpair keypairs;
+ int keys_per_elem;
public:
physical_elements_factory(const entity_params& ep);
diff --git a/test.cpp b/test.cpp
index a18a62bd50d49e6900c1816bb594c010d9de2376..ed5506a90bb036eb6d1bc4bdb0c6f3d495b99d88 100644
--- a/test.cpp
+++ b/test.cpp
@@ -28,7 +28,6 @@ int main(void)
int failed_tests = 0;
-
#if 0
std::cout << Bmagentafg << " *** TESTING: BASIC OPERATIONS ***" << reset << std::endl;
for (size_t go = 1; go < 5; go++)
@@ -39,6 +38,9 @@ int main(void)
}
}
+ std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION ***" << reset << std::endl;
+ for (size_t i = 1; i < 5; i++)
+ test_normals(i,i);
std::cout << Bmagentafg << " *** TESTING: PROJECTIONS ***" << reset << std::endl;
for (size_t go = 1; go < 5; go++)
@@ -46,19 +48,17 @@ int main(void)
for (size_t ao = go; ao < 5; ao++)
failed_tests += test_mass_convergence(go, ao);
}
-#endif
std::cout << Bmagentafg << " *** TESTING: DIFFERENTIATION ***" << reset << std::endl;
for (size_t go = 1; go < 5; go++)
for (size_t ao = go; ao < 5; ao++)
failed_tests += test_differentiation_convergence(go, ao);
+#endif
+
+ std::cout << Bmagentafg << " *** TESTING: LIFTING ***" << reset << std::endl;
+ for (size_t ao = 1; ao < 2; ao++)
+ test_lifting(1, ao);
- /*****************************************/
-
- std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION ***" << reset << std::endl;
- for (size_t i = 1; i < 5; i++)
- test_normals(i,i);
-
gmsh::finalize();
return failed_tests;
diff --git a/test.h b/test.h
index b7e7ee54895010f74f8f65d2394917e20d481d92..8e376469e1fa0bfbd76eff5b8bcd34f3d5ca6dd0 100644
--- a/test.h
+++ b/test.h
@@ -28,4 +28,4 @@ int test_basics(int, int);
int test_normals(int, int);
int test_mass_convergence(int, int);
int test_differentiation_convergence(int, int);
-
+int test_lifting(int, int);
diff --git a/test_basics.cpp b/test_basics.cpp
index 4f1d63c149509a35c8e56fd5875bb2c27c475e60..806303d3b28679d3332bed19389018a610fa86d0 100644
--- a/test_basics.cpp
+++ b/test_basics.cpp
@@ -1,13 +1,15 @@
#include <iostream>
-
#include <unistd.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
#include "test.h"
#include "gmsh_io.h"
-#include "entity_data.h"
-
#include "sgr.hpp"
+
using namespace sgr;
static void
@@ -189,6 +191,7 @@ test_normals(int geometric_order, int approximation_order)
std::sort(faceTags.begin(), faceTags.end());
auto& e0 = mod.entity_at(0);
+
entity_data_cpu ed;
e0.populate_entity_data(ed);
auto ft = e0.face_tags();
@@ -231,4 +234,4 @@ test_normals(int geometric_order, int approximation_order)
}
return 0;
-}
\ No newline at end of file
+}
diff --git a/test_lifting.cpp b/test_lifting.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..61158f1765b756159ece0369b50655374d754525
--- /dev/null
+++ b/test_lifting.cpp
@@ -0,0 +1,144 @@
+#include <iostream>
+#include <iomanip>
+
+#include <unistd.h>
+
+#include "test.h"
+#include "gmsh_io.h"
+#include "sgr.hpp"
+#include "entity_data.h"
+#include "kernels_cpu.h"
+#include "timecounter.h"
+
+using namespace sgr;
+
+static void
+make_geometry(int order, double mesh_h)
+{
+ gm::add("difftest");
+ gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0);
+ gmo::synchronize();
+
+ gvp_t vp;
+ gm::getEntities(vp);
+ gmm::setSize(vp, mesh_h);
+}
+
+
+int test_lifting(int geometric_order, int approximation_order)
+{
+ std::vector<double> sizes({ 0.32, 0.16, 0.08, 0.04 });
+ std::vector<double> errors;
+
+ std::cout << cyanfg << "Testing geometric order " << geometric_order;
+ std::cout << ", approximation order = " << approximation_order << nofg;
+ std::cout << std::endl;
+
+ auto F = [](const point_3d& pt) -> vecxd {
+ vec3d Fret;
+ Fret(0) = std::sin(M_PI*pt.y());
+ Fret(1) = std::sin(M_PI*pt.z());
+ Fret(2) = std::sin(M_PI*pt.x());
+ return Fret;
+ };
+
+ auto div_F = [](const point_3d& pt) -> double {
+ return 0;
+ auto dFx_dx = M_PI*std::cos(M_PI*pt.x());
+ auto dFy_dy = M_PI*std::cos(M_PI*pt.y());
+ auto dFz_dz = M_PI*std::cos(M_PI*pt.z());
+ return dFx_dx + dFy_dy + dFz_dz;
+ };
+
+ for (size_t 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();
+
+ entity_data_cpu ed;
+ e0.populate_entity_data(ed);
+
+ vecxd exp_field = e0.project(div_F);
+ vecxd lift_field = vecxd::Zero( e0.num_cells() * ed.num_bf );
+ vecxd flux = vecxd::Zero( e0.num_cells() * 4*ed.num_fluxes );
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ for (size_t iF = 0; iF < 4; iF++)
+ {
+ auto& pf = e0.face(4*iT+iF);
+ auto& rf = e0.face_refelem(pf);
+ auto num_fluxes = rf.num_basis_functions();
+ const auto pqps = pf.integration_points();
+
+ matxd mass = matxd::Zero(num_fluxes, num_fluxes);
+ vecxd rhs = vecxd::Zero(num_fluxes);
+ for (size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ vec3d n;
+ if (ed.g_order == 1)
+ n = ed.normals.row(4*iT+iF);
+ else
+ n = ed.normals.row( (4*iT+iF)*pqps.size() + iQp );
+
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = rf.basis_functions(iQp);
+ mass += pqp.weight() * phi * phi.transpose();
+ rhs += pqp.weight() * F(pqp.point()).dot(n) * phi;
+ }
+ auto fluxofs = (4*iT+iF)*num_fluxes;
+ flux.segment(fluxofs, num_fluxes) = mass.ldlt().solve(rhs);
+ }
+ }
+
+ timecounter tc;
+ tc.tic();
+ compute_flux_lifting(ed, flux, lift_field);
+ double time = tc.toc();
+
+ double err = 0.0;
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
+ {
+ auto& pe = e0.cell(iT);
+ auto& re = e0.cell_refelem(pe);
+ auto num_bf = re.num_basis_functions();
+ matxd mass = matxd::Zero(num_bf, num_bf);
+ vecxd vol = vecxd::Zero(num_bf);
+
+ const auto pqps = pe.integration_points();
+ for(size_t iQp = 0; iQp < pqps.size(); iQp++)
+ {
+ const auto& pqp = pqps[iQp];
+ const vecxd phi = re.basis_functions(iQp);
+ const matxd dphi = re.basis_gradients(iQp);
+ mass += pqp.weight() * phi * phi.transpose();
+ vol += pqp.weight() * dphi * F(pqp.point());
+ }
+
+ vecxd comp_field = lift_field.segment(iT*num_bf, num_bf) - mass.ldlt().solve(vol);
+ vecxd diff = comp_field - exp_field.segment(iT*num_bf, num_bf);
+ err += diff.dot(mass*diff);
+ }
+
+ std::cout << std::sqrt(err) << std::endl;
+
+ if (geometric_order == 1)
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = 3*(ed.num_bf)*4*ed.num_fluxes*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ else
+ {
+ std::cout << "Kernel runtime: " << time << " seconds. Estimated performance: ";
+ double flops = ((21*ed.num_bf+6)*ed.num_bf*ed.num_qp + 3*(2*ed.num_bf-1)*ed.num_bf)*e0.num_cells();
+ std::cout << flops/(1e9*time) << " GFlops/s" << std::endl;
+ }
+ }
+
+ return 0;
+}
\ No newline at end of file
diff --git a/types.h b/types.h
index 100670e37e55a9ebd0880dee4c684ab05dd1fc49..cb89d6cccb39cfb2cf3ae102278b8fc2be7faaf1 100644
--- a/types.h
+++ b/types.h
@@ -8,11 +8,13 @@ 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 vec3d = Eigen::Matrix<double, 3, 1>;
using mat3d = Eigen::Matrix<double, 3, 3, Eigen::RowMajor>;
using vecxd = Eigen::Matrix<double, Eigen::Dynamic, 1>;
using matxd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
using vec_mat3d = eigen_compatible_stdvector<mat3d>;
using vec_matxd = eigen_compatible_stdvector<matxd>;
+
struct entity_params
{
int dim;