366 surface mesh reorient (#372)

* Init surface reorientation

* Surface reorientation implemented in 3D

* update doc for surfaceReorient

---------

Co-authored-by: franck.ledoux <[email protected]>

docs/mkdocs/user-guide/igalgo.md

@@ -5,4 +5,7 @@ This module provides some basic algorithms for the incidence graph (*ig*) data s
 - Most important algorithms are relative to retrieving and building some boundary informations: boundaryOperator, boundaryExtractor.
 - The *GridBuilder* class provides algoriths to build 2D and 3D structured grids for unit testing purposes mainly.
 - The *THexBuilder* create a hexahedral mesh by splitting each tet of an input tetrahedral mesh.
-- The *MeshQualityCompute* go through all the mesh cells and assign variables depending on a set of selected quality criteria
+- The *MeshQualityCompute* go through all the mesh cells and assign variables depending on a set of selected quality criteria
+- The *SurfaceReorient* class reorient surface mesh (in 2D and 3D) in order to have all faces oriented in the same direction. If
+the mesh is split in different connex components, the orientation between different components can be different.
+In practice, the orientation is consistent in 2D.

igalgo/CMakeLists.txt

@@ -12,13 +12,15 @@ set(GMDS_INC
         inc/gmds/igalgo/BoundaryOperator.h
         inc/gmds/igalgo/BoundaryOperator2D.h
         inc/gmds/igalgo/GridBuilder.h
+        inc/gmds/igalgo/SurfaceReorient.h
         inc/gmds/igalgo/THexBuilder.h)
 set(GMDS_SRC
         src/BoundaryExtractor3D.cpp
         src/BoundaryExtractor2D.cpp
         src/BoundaryOperator.cpp
         src/BoundaryOperator2D.cpp
         src/GridBuilder.cpp
+        src/SurfaceReorient.cpp
         src/THexBuilder.cpp)
 
 #==============================================================================

igalgo/inc/gmds/igalgo/SurfaceReorient.h

@@ -0,0 +1,62 @@
+/*----------------------------------------------------------------------------*/
+#ifndef GMDS_SURFACE_REORIENT_H
+#define GMDS_SURFACE_REORIENT_H
+/*----------------------------------------------------------------------------*/
+#include <gmds/ig/Mesh.h>
+#include "GMDSIgAlgo_export.h"
+/*----------------------------------------------------------------------------*/
+#include <map>
+#include <vector>
+/*----------------------------------------------------------------------------*/
+namespace gmds{
+
+    /*----------------------------------------------------------------------------*/
+    /** @class  SurfaceReorient
+     *  @brief  Class that provides an algorithm to reorient surfaces in a consistent
+     *  manner
+     */
+    class GMDSIgAlgo_API SurfaceReorient
+    {
+    public:
+
+        /*------------------------------------------------------------------------*/
+        /** @brief Constructor.
+         *
+         * @param Amesh the mesh to work on
+         * @param ADim the grid dimension: 2 or 3 (default). Dimension of 2 means
+         * that we know we work in the 2D plane. In this case the orientation is
+         * consistent with FEM requirements.
+         */
+	   SurfaceReorient(Mesh* AMesh, const TInt ADim=3);
+
+        /*------------------------------------------------------------------------*/
+        /** @brief  Destructor.	*/
+        virtual ~SurfaceReorient();
+
+        /*------------------------------------------------------------------------*/
+        /** @brief Check if the mesh fits algorithm requirements, which are:
+         *         - 3D mesh with R, N and R2N
+         *         - 2D mesh with F, N and F2N
+         */
+        bool isValid() const;
+        /*------------------------------------------------------------------------*/
+        /** @brief  Performs the reorientation algorithm.
+         * @return the number of faces that have been reoriented
+         */
+        int execute();
+
+	   private:
+	     int orient2d();
+	     int orient3d();
+	     bool orient2d(gmds::Face& AF);
+	     static TCoord isLeft(Node& AN1, Node& AN2, Node& AN3);
+    private:
+        /** a mesh */
+        Mesh* m_mesh;
+        /** problem dimension. 2 means we have a pure (X,Y,0) plan*/
+        TInt  m_dim;
+    };
+    /*----------------------------------------------------------------------------*/
+}
+/*----------------------------------------------------------------------------*/
+#endif //GMDS_SURFACE_REORIENT_H

igalgo/src/SurfaceReorient.cpp

@@ -0,0 +1,165 @@
+/*----------------------------------------------------------------------------*/
+#include <gmds/igalgo/SurfaceReorient.h>
+/*----------------------------------------------------------------------------*/
+#include <sstream>
+/*----------------------------------------------------------------------------*/
+using namespace gmds;
+/*----------------------------------------------------------------------------*/
+SurfaceReorient::SurfaceReorient(Mesh* AMesh, const TInt ADim)
+        :m_mesh(AMesh), m_dim(ADim)
+{}
+/*----------------------------------------------------------------------------*/
+SurfaceReorient::~SurfaceReorient()
+{}
+/*----------------------------------------------------------------------------*/
+bool SurfaceReorient::isValid() const
+{
+    if(m_dim==3)
+        return (m_mesh->getModel()==(DIM3|F|N|F2N|N2F));
+    else if(m_dim==2)
+        return (m_mesh->getModel()==(DIM3|F|N|F2N) ||
+                m_mesh->getModel()==(DIM2|F|N|F2N));
+
+    //dimension error
+    return false;
+}
+/*----------------------------------------------------------------------------*/
+int SurfaceReorient::execute() {
+	 return (m_dim==2)?orient2d():orient3d();
+}
+/*----------------------------------------------------------------------------*/
+int SurfaceReorient::orient2d()
+{
+	 auto nb_reorientation =0;
+	 for(auto f_id:m_mesh->faces()){
+		  Face f = m_mesh->get<Face>(f_id);
+		  if (orient2d(f))
+			  nb_reorientation++;
+	 }
+	 return nb_reorientation;
+}
+/*------------------------------------------------------------------------*/
+bool SurfaceReorient::orient2d(Face& AF)
+{
+	 bool isReoriented = false;
+	 std::vector<Node> nodes = AF.get<Node>();
+	 TCoord orientation=0;
+	 if(AF.type()==GMDS_TRIANGLE) {
+		  orientation = isLeft(nodes[0],nodes[1],nodes[2]);
+	 }
+	 else {
+		  //find the rightmost lowest vertex of the polygon
+		  unsigned int index_min=0;
+		  TCoord x_min = nodes[0].X();
+		  TCoord y_min = nodes[0].Y();
+		  for(unsigned int
+		          i=0;i<nodes.size();i++)
+		  {
+			  if(nodes[i].Y()>y_min)
+				  continue;
+			  if(nodes[i].Y()==y_min) {	// just as low
+				  if(nodes[i].X()<x_min)  // and to left
+					  continue;
+			  }
+
+			  index_min =i;
+			  x_min = nodes[i].X();
+			  y_min = nodes[i].Y();
+		  }
+
+		  if(index_min==0)
+			  orientation = isLeft(nodes[nodes.size()-1],nodes[0],nodes[1]);
+		  else if (index_min==nodes.size()-1)
+			  orientation = isLeft(nodes[index_min-1],nodes[index_min],nodes[0]);
+		  else
+			  orientation = isLeft(nodes[index_min-1],nodes[index_min],nodes[index_min+1]);
+	 }
+	 if(orientation>0.0) // clockwise or degenerated (=0)
+	 {
+		  isReoriented= true;
+		  std::vector<Node> nodes_inv;
+		  nodes_inv.resize(nodes.size());
+		  auto node_size = nodes.size();
+		  for(unsigned int i=0;i<node_size;i++)
+			  nodes_inv[i] = nodes[node_size-1-i];
+		  AF.set<Node>(nodes_inv);
+	 }
+	 return isReoriented;
+}
+/*----------------------------------------------------------------------------*/
+TCoord SurfaceReorient::isLeft(Node& AN1, Node& AN2, Node& AN3)
+{
+	 return ( (AN2.X()-AN1.X()) * (AN3.Y()-AN1.Y()) -
+	         (AN3.X()-AN1.X()) * (AN2.Y()-AN1.Y()) );
+}
+/*----------------------------------------------------------------------------*/
+int SurfaceReorient::orient3d()
+{
+	 auto nb_reorientation =0;
+	 auto mark_done = m_mesh->newMark<Face>();
+
+	 bool keep_working =true;
+	 while (keep_working) {
+		 //We look for the first unmarked face
+
+		 auto it_face = m_mesh->faces_begin();
+		 while (m_mesh->isMarked<Face>(*it_face,mark_done) &&
+		    it_face!=m_mesh->faces_end()){
+			 ++it_face;
+		 }
+		 if(it_face==m_mesh->faces_end()){
+			 keep_working=false;
+		 }
+		 if(keep_working) {
+			 // We take the first face as the orientation reference seed
+			 auto seed_id = *it_face;
+			 Face seed = m_mesh->get<Face>(seed_id);
+			    // And now, we go through all the faces in an advancing-front manner to
+			    // orient faces according to the seed
+			    std::vector<TCellID> front;
+			 front.push_back(seed.id());
+			 m_mesh->mark(seed, mark_done);
+			 while (!front.empty()) {
+				 auto current_id = front.back();
+				 front.pop_back();
+				 auto current_face = m_mesh->get<Face>(current_id);
+				 // the current face is considered as well oriented and we orient the face sharing
+				 //  an edge with it in a valid manner
+				 std::vector<Node> node_faces = current_face.get<Node>();
+				 for (auto i = 0; i < node_faces.size(); i++) {
+					 Node ni = node_faces[i];
+					 Node nj = node_faces[(i + 1) % node_faces.size()];
+					 auto fij_ids = m_mesh->getCommonFaces(ni, nj);
+					 if (fij_ids.size() == 2) {
+						 // means there is another face
+						 auto other_face_id = (fij_ids[0] == current_id) ? fij_ids[1] : fij_ids[0];
+						 auto other_face = m_mesh->get<Face>(other_face_id);
+						 if (!m_mesh->isMarked(other_face, mark_done)) {
+							 std::vector<TCellID> other_face_node_ids = other_face.getIDs<Node>();
+							 // we check if the node ni and nj are traversed in the same way.
+							 bool same_direction = false;
+							 for (auto k = 0; k < other_face_node_ids.size(); k++) {
+
+								 auto id_k = other_face_node_ids[k];
+								 auto id_l = other_face_node_ids[(k + 1) % other_face_node_ids.size()];
+								 if (id_k == ni.id() && id_l == nj.id()) same_direction = true;
+							 }
+							 if (same_direction) {
+								 std::reverse(other_face_node_ids.begin(), other_face_node_ids.end());
+								 other_face.set<Node>(other_face_node_ids);
+								 nb_reorientation++;
+							 }
+							 m_mesh->mark(other_face, mark_done);
+							 front.push_back(other_face_id);
+						 }
+					 }
+				 }
+			 }
+		 }
+	 }
+	 m_mesh->negateMaskMark<Node>(mark_done);
+	 m_mesh->freeMark<Node>(mark_done);
+	 //warning, we must ensure that the mesh is reoriented even if it has several
+	 //connex part
+	 return nb_reorientation;
+}

igalgo/tst/SurfaceReorientTestSuite.h

@@ -0,0 +1,39 @@
+/*----------------------------------------------------------------------------*/
+#include <gtest/gtest.h>
+#include <gmds/ig/Mesh.h>
+#include <gmds/ig/MeshDoctor.h>
+#include <gmds/io/VTKWriter.h>
+#include <gmds/io/VTKReader.h>
+#include <gmds/igalgo/SurfaceReorient.h>
+#include <iostream>
+#include <gmds/io/IGMeshIOService.h>
+#include <unit_test_config.h>
+using namespace gmds;
+/*----------------------------------------------------------------------------*/
+TEST(SurfaceReorientTestClass, testReorient)
+{
+    // WE WRITE
+    Mesh m(MeshModel(DIM3|F|N|F2N|N2F));
+
+    std::string dir(TEST_SAMPLES_DIR);
+    std::string vtk_file = dir+"/reorient3D_test.vtk";
+    gmds::IGMeshIOService ioService(&m);
+    gmds::VTKReader vtkReader(&ioService);
+    vtkReader.setCellOptions(gmds::N|gmds::F);
+    vtkReader.read(vtk_file);
+
+    gmds::MeshDoctor doc(&m);
+    doc.updateUpwardConnectivity();
+    gmds::SurfaceReorient reorient(&m);
+
+    ASSERT_TRUE(reorient.isValid());
+
+	 auto nb_reorient = reorient.execute();
+    ASSERT_TRUE(nb_reorient!=0);
+
+	 gmds::VTKWriter w(&ioService);
+	 w.setCellOptions(gmds::N|gmds::F);
+	 w.write("toto.vtk");
+
+
+}

igalgo/tst/main_test.cpp

@@ -6,6 +6,7 @@
 #include "BoundaryOpTestSuite.h"
 #include "BoundaryExtractorTestSuite.h"
 #include "GridBuilderTestSuite.h"
+#include "SurfaceReorientTestSuite.h"
 /*----------------------------------------------------------------------------*/
 int main(int argc, char ** argv) {
   ::testing::InitGoogleTest(&argc, argv);