Face keys, normals test.

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();

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

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();

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;
+}

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;
+};
 
 
 

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();
 

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