diff --git a/entity.cpp b/entity.cpp
index c7d9b88dcf68844784d487e15c10bcde53b3c74f..fb246a2b5f4e5a7ba5d7d5f5b5217537d317b7e6 100644
--- a/entity.cpp
+++ b/entity.cpp
@@ -386,7 +386,7 @@ entity::populate_normals(entity_data_cpu& ed) const
{
auto fnum = 4*iT+iF;
assert(fnum < physical_faces.size());
- Eigen::Vector3d n = physical_faces[fnum].jacobian(iQp).inverse()*ref_normal;
+ Eigen::Vector3d n = physical_faces[fnum].jacobian(iQp)*ref_normal;
auto norm_offset = (4*iT+iF)*num_face_qps + iQp;
assert(norm_offset < ed.normals.rows());
ed.normals.block(norm_offset, 0, 1, 3) = n.transpose()/n.norm();
diff --git a/gmsh_io.cpp b/gmsh_io.cpp
index 59663b0cf1187ea398bcd27d36377fcc4d4cd55a..dbde114eeedeb164b7463c44f9438e759d440d89 100644
--- a/gmsh_io.cpp
+++ b/gmsh_io.cpp
@@ -77,6 +77,7 @@ model::remesh()
gmm::setOrder( geometric_order );
entities.clear();
import_gmsh_entities();
+ map_boundaries();
}
void
@@ -155,18 +156,8 @@ model::map_boundaries(void)
assert(elemTypes.size() == 1);
for (auto& elemType : elemTypes)
{
- std::vector<size_t> nTags;
- gmm::getElementFaceNodes(elemType, 3, nTags, tag, true);
-
- std::vector<face_key> fkeys;
- for (size_t i = 0; i < nTags.size(); i+= 3)
- {
- auto fk = face_key({nTags[i+0], nTags[i+1], nTags[i+2]});
- fkeys.push_back( fk );
- }
- std::sort(fkeys.begin(), fkeys.end());
-
- boundaries_elemTags[tag] = std::move(fkeys);
+ element_key_factory ekf(DIMENSION(2), tag, elemType);
+ boundary_map[tag] = ekf.get_keys();
}
}
}
\ No newline at end of file
diff --git a/gmsh_io.h b/gmsh_io.h
index ee27d25528c5d510bb9d9234eeac6a06b31086d1..6bb94845742d9a4ffeb2e7ebbc7a788d28974cc0 100644
--- a/gmsh_io.h
+++ b/gmsh_io.h
@@ -31,7 +31,8 @@ class model
int approximation_order;
- std::vector<entity> entities;
+ std::vector<entity> entities;
+ std::map<size_t, std::vector<element_key>> boundary_map;
void map_boundaries(void);
void import_gmsh_entities(void);
@@ -46,6 +47,11 @@ public:
void remesh(void);
void remesh(int);
+ std::vector<element_key> get_bnd(size_t which)
+ {
+ return boundary_map.at(which);
+ }
+
std::vector<entity>::const_iterator begin() const;
std::vector<entity>::const_iterator end() const;
std::vector<entity>::iterator begin();
diff --git a/physical_element.cpp b/physical_element.cpp
index 6a866c171ed3aee97fcf4c54ebd5a4cb3369d132..d806625292a9b953bc88421c3f39b030b89beb4f 100644
--- a/physical_element.cpp
+++ b/physical_element.cpp
@@ -44,6 +44,18 @@ physical_element::element_tag() const
return m_element_tag;
}
+int
+physical_element::element_type() const
+{
+ return m_gmsh_elemtype;
+}
+
+vec_size_t
+physical_element::node_tags(void) const
+{
+ return m_node_tags;
+}
+
/* Return only one determinant: this is used for geometric order 1
* where the determinants and the jacobians are constant */
double
@@ -105,13 +117,6 @@ physical_element::jacobians(void) const
-
-
-
-
-
-
-
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)
@@ -148,6 +153,7 @@ physical_elements_factory::get_elements()
new_pe.m_orientation = cellOrientations[elem];
new_pe.m_element_tag = cellTags[elem];
new_pe.m_measure = 0.0;
+ new_pe.m_gmsh_elemtype = elemType;
new_pe.m_node_tags.resize(num_nodes);
for (size_t i = 0; i < num_nodes; i++)
@@ -194,3 +200,103 @@ physical_elements_factory::get_elements()
return ret;
}
+
+element_key::element_key()
+{}
+
+element_key::element_key(const physical_element& pe)
+{
+ m_dim = pe.dimension();
+ m_elemType = pe.element_type();
+
+ auto nodes = pe.node_tags();
+
+ if (m_dim == 2)
+ {
+ m_key_data[0] = 3;
+ m_key_data[1] = nodes[0];
+ m_key_data[2] = nodes[1];
+ m_key_data[3] = nodes[2];
+ std::sort(std::next(m_key_data.begin()), m_key_data.end());
+ return;
+ }
+
+ std::stringstream ss;
+ ss << "Key not implemented for element type " << m_elemType;
+ throw std::invalid_argument(ss.str());
+}
+
+bool
+element_key::operator<(const element_key& other) const
+{
+ if (m_dim != other.m_dim or m_elemType != other.m_elemType)
+ return false;
+
+ if (m_key_data[0] != other.m_key_data[0])
+ return false;
+
+ auto mb = std::next(m_key_data.begin());
+ auto me = m_key_data.end();
+ auto ob = std::next(other.m_key_data.begin());
+ auto oe = other.m_key_data.end();
+
+ return std::lexicographical_compare(mb, me, ob, oe);
+}
+
+bool
+element_key::operator==(const element_key& other) const
+{
+ if (m_dim != other.m_dim or m_elemType != other.m_elemType)
+ return false;
+
+ if (m_key_data[0] != other.m_key_data[0])
+ return false;
+
+ auto mb = std::next(m_key_data.begin());
+ auto me = m_key_data.end();
+ auto ob = std::next(other.m_key_data.begin());
+
+ return std::equal(mb, me, ob);
+}
+
+
+std::ostream& operator<<(std::ostream& os, const element_key& ek)
+{
+ os << "ek: " << ek.m_dim << " " << ek.m_elemType << " (";
+ for (size_t i = 0; i < 8; i++)
+ os << ek.m_key_data[i] << " ";
+ os << ")";
+ return os;
+}
+
+
+
+element_key_factory::element_key_factory(int dim, int tag, int etype)
+{
+ assert(dim == 2);
+
+ std::vector<size_t> nTags;
+ gmm::getElementFaceNodes(etype, 3, nTags, tag, true);
+
+ for (size_t i = 0; i < nTags.size(); i+= 3)
+ {
+ element_key ek;
+ ek.m_dim = dim;
+ ek.m_elemType = etype;
+ ek.m_key_data[0] = 3;
+ ek.m_key_data[1] = nTags[i+0];
+ ek.m_key_data[2] = nTags[i+1];
+ ek.m_key_data[3] = nTags[i+2];
+
+ std::sort(std::next(ek.m_key_data.begin()), ek.m_key_data.end());
+
+ ekeys.push_back( ek );
+ }
+ std::sort(ekeys.begin(), ekeys.end());
+}
+
+std::vector<element_key>
+element_key_factory::get_keys(void) const
+{
+ return ekeys;
+}
diff --git a/physical_element.h b/physical_element.h
index c86f54611c9cc7e6123a4d3790173a9a849d1784..819fa91d69c68845f6a3f102a63c044a6589f1b4 100644
--- a/physical_element.h
+++ b/physical_element.h
@@ -9,6 +9,7 @@ class physical_element
int m_dimension; /* Geometric dimension */
int m_orientation; /* GMSH orientation */
int m_parent_entity_tag; /* Tag of the parent entity */
+ int m_gmsh_elemtype; /* GMSH element type */
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 */
@@ -27,6 +28,7 @@ public:
int orientation(void) const;
int parent_entity_tag(void) const;
int element_tag(void) const;
+ int element_type(void) const;
vec_size_t node_tags(void) const;
double determinant(void) const;
@@ -70,8 +72,38 @@ public:
};
+class element_key
+{
+ int m_dim;
+ int m_elemType;
+ std::array<size_t, 8> m_key_data{};
+public:
+ element_key();
+ element_key(const element_key&) = default;
+ element_key(element_key&&) = default;
+ element_key(const physical_element&);
+
+ element_key& operator=(const element_key&) = default;
+
+ bool operator<(const element_key&) const;
+ bool operator==(const element_key&) const;
+
+ friend class element_key_factory;
+ friend std::ostream& operator<<(std::ostream& os, const element_key& ek);
+};
+std::ostream& operator<<(std::ostream& os, const element_key& ek);
+
+class element_key_factory
+{
+ std::vector<element_key> ekeys;
+
+public:
+ element_key_factory(int dim, int tag, int etype);
+
+ std::vector<element_key> get_keys(void) const;
+};
diff --git a/test.cpp b/test.cpp
index 8d0922766150cfc445d95595468817378d5f5676..1aaf28a71fc7fc7eb0816f6499ac7db49b3cc937 100644
--- a/test.cpp
+++ b/test.cpp
@@ -29,7 +29,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++)
{
@@ -52,10 +51,11 @@ int main(void)
for (size_t ao = go; ao < 5; ao++)
failed_tests += test_differentiation_convergence(go, ao);
-#endif
/*****************************************/
- test_normals(2,2);
+ std::cout << Bmagentafg << " *** TESTING: NORMAL COMPUTATION ***" << reset << std::endl;
+ for (size_t i = 1; i < 5; i++)
+ test_normals(i,i);
gmsh::finalize();
diff --git a/test_basics.cpp b/test_basics.cpp
index fe78195fb571bfe271fcc7e210f3776788f27744..4f1d63c149509a35c8e56fd5875bb2c27c475e60 100644
--- a/test_basics.cpp
+++ b/test_basics.cpp
@@ -165,6 +165,9 @@ test_normals(int geometric_order, int approximation_order)
gm::getPhysicalGroups(pgs);
assert(pgs.size() == 1);
+ std::cout << cyanfg << "Testing normals (order " << geometric_order;
+ std::cout << ")" << reset << std::endl;
+
for (auto [dim, pgtag] : pgs)
{
std::vector<int> tags;
@@ -180,7 +183,6 @@ test_normals(int geometric_order, int approximation_order)
for (auto& elemType : elemTypes)
{
- std::ofstream ofs_vecs("ofs_vecs.txt");
std::vector<size_t> faceTags, faceNodesTags;
gmm::getElementsByType(elemType, faceTags, faceNodesTags, tags[0]);
@@ -191,7 +193,9 @@ test_normals(int geometric_order, int approximation_order)
e0.populate_entity_data(ed);
auto ft = e0.face_tags();
- for (size_t iT = 0; iT < 10; iT++)
+ auto bm = mod.get_bnd(7);
+
+ for (size_t iT = 0; iT < e0.num_cells(); iT++)
{
for (size_t iF = 0; iF < 4; iF++)
{
@@ -199,23 +203,32 @@ test_normals(int geometric_order, int approximation_order)
auto rf = e0.face_refelem(pf);
auto qps = pf.integration_points();
- std::cout << pf.element_tag() << std::endl;
+ auto ek = element_key(pf);
+ if ( !std::binary_search(bm.begin(), bm.end(), ek) )
+ continue;
+
+ size_t num_qp = (geometric_order == 1) ? 1 : qps.size();
- for (size_t iQp = 0; iQp < qps.size(); iQp++)
+ for (size_t iQp = 0; iQp < num_qp; iQp++)
{
- Eigen::Vector3d n = ed.normals.row((4*iT+iF)*qps.size() + iQp);
+ Eigen::Vector3d n = ed.normals.row((4*iT+iF)*num_qp + iQp);
point_3d p = qps[iQp].point();
point_3d v = p - point_3d(0.5, 0.5, 0.5);
Eigen::Vector3d nref;
nref(0) = v.x(); nref(1) = v.y(); nref(2) = v.z();
nref /= nref.norm();
- std::cout << n.cross(nref).transpose() << std::endl;
- ofs_vecs << p.x() << " " << p.y() << " " << p.z() << " ";
- ofs_vecs << v.x() << " " << v.y() << " " << v.z() << std::endl;
+ if ( n.cross(nref).norm() > 1e-2 )
+ {
+ std::cout << "Wrong normal. Expected: " << nref.transpose() << ", ";
+ std::cout << "got: " << n.transpose() << std::endl;
+ return 1;
+ }
}
}
}
}
}
+
+ return 0;
}
\ No newline at end of file