refine one quad to 4 triangles

NonLinearSolver/nlsolver/meshModification.cpp

@@ -6,6 +6,8 @@
 //
 #include "meshModification.h"
 #include "Context.h"
+#include "MElement.h"
+#include "MTriangle.h"
 
 TrivialOperation::TrivialOperation():modificationOperation(){};
 
@@ -117,7 +119,7 @@ void meshModification::applyModification(const std::string fname, double version
 	Msg::Info("End modifying mesh");
 }
 
-void makePhysicalByBox::write_MSH2(const std::map<int, elementGroup>& gMap,  const std::string filename)
+void write_MSH2(const std::map<int, elementGroup>& gMap,  const std::string filename)
 {
   std::set<MVertex*> allVertices;
   int totalNumberOfElements = 0;
@@ -149,7 +151,14 @@ void makePhysicalByBox::write_MSH2(const std::map<int, elementGroup>& gMap,  con
     Msg::Info("num of eelement in gr %d = %ld",itg->first,g->size());
     for (elementGroup::elementContainer::const_iterator it = g->begin(); it!= g->end(); it++){
       MElement* ele = it->second;
-      ele->writeMSH2(fp,2.2,0,0,itg->first+1,itg->first+1);
+      fprintf(fp, "%ld %d 2 %ld %d",ele->getNum(),ele->getTypeForMSH(),itg->first,itg->first);
+      std::vector<MVertex*> vers;
+      ele->getVertices(vers);
+      for (int i=0; i< vers.size(); i++)
+      {
+        fprintf(fp, " %ld",vers[i]->getNum());
+      }
+      fprintf(fp, "\n");
     };
   }
 	fprintf(fp, "$EndElements\n");
@@ -285,3 +294,111 @@ void makePhysicalByBox::run(const std::string inputMeshFile, const std::string o
   
   write_MSH2(allG,outputMeshFile);
 };
+
+
+QuadToTriangle::QuadToTriangle(): _pModel(NULL){};
+QuadToTriangle::~QuadToTriangle()
+{
+  if (_pModel) delete _pModel;
+};
+
+void QuadToTriangle::loadModel(const std::string meshFileName)
+{
+  if (_pModel) delete _pModel;
+  _pModel = new GModel();
+  _pModel->readMSH(meshFileName.c_str());
+};
+
+void QuadToTriangle::remesh(MElement* ele, std::vector<MElement*>& newEle) const
+{
+  if (ele->getTypeForMSH() == MSH_QUA_4)
+  {
+    std::vector<MVertex*> vv;
+    ele->getVertices(vv);
+    
+    std::vector<double > x(3,0.);
+    for (int i=0; i< 4; i++)
+    {
+      x[0] += 0.25*vv[i]->x();
+      x[1] += 0.25*vv[i]->y();
+      x[2] += 0.25*vv[i]->z();
+    }
+    MVertex* newVer = new MVertex(x[0],x[1],x[2]);
+    
+    Msg::Info("add vertex %d",newVer->getNum());
+    
+    MElement* e1 = new MTriangle(vv[0],vv[1],newVer);
+    MElement* e2 = new MTriangle(vv[1],vv[2],newVer);
+    MElement* e3 = new MTriangle(vv[2],vv[3],newVer);
+    MElement* e4 = new MTriangle(vv[3],vv[0],newVer);
+    
+    newEle.push_back(e1);
+    newEle.push_back(e2);
+    newEle.push_back(e3);
+    newEle.push_back(e4);
+  }
+  else
+  {
+    Msg::Error("MSH_QUA_4 must be used");
+  }
+};
+
+void QuadToTriangle::apply(const std::string fname)
+{
+  Msg::Info("start modifying mesh");
+  
+  int dim = _pModel->getNumRegions() ? 3 : 2;
+  std::map<int, std::vector<GEntity*> > groups[4];
+  _pModel->getPhysicalGroups(groups);
+  std::map<int, std::vector<GEntity*> > &entmap = groups[dim];
+  
+  std::map<int, elementGroup> allG;  
+
+  for (std::map<int, std::vector<GEntity*> >::iterator it = entmap.begin(); it!= entmap.end(); it++){
+    std::vector<GEntity*> &ent = it->second;
+    for (unsigned int i = 0; i < ent.size(); i++){
+      for (unsigned int j = 0; j < ent[i]->getNumMeshElements(); j++){
+        MElement *e = ent[i]->getMeshElement(j);
+        std::vector<MElement*> enew;
+        remesh(e, enew);
+        for (int j=0; j< enew.size(); j++)
+        {
+          elementGroup& gr = allG[it->first];
+          gr.insert(enew[j]);
+        }
+        
+
+      }
+    }
+  }
+  
+  std::map<int, std::vector<GEntity*> > &entmap1D = groups[1];
+  Msg::Info("size of entmap1D = %d",entmap1D.size());
+  for (std::map<int, std::vector<GEntity*> >::iterator it = entmap1D.begin(); it!= entmap1D.end(); it++){
+    std::vector<GEntity*> &ent = it->second;
+    for (unsigned int i = 0; i < ent.size(); i++){
+      for (unsigned int j = 0; j < ent[i]->getNumMeshElements(); j++){
+        MElement *e = ent[i]->getMeshElement(j);
+        elementGroup& gr = allG[it->first];
+        gr.insert(e);          
+      }
+    }
+  }
+  
+  std::map<int, std::vector<GEntity*> > &entmap0D = groups[0];
+  Msg::Info("size of entmap0D = %d",entmap0D.size());
+  for (std::map<int, std::vector<GEntity*> >::iterator it = entmap0D.begin(); it!= entmap0D.end(); it++){
+    std::vector<GEntity*> &ent = it->second;
+    for (unsigned int i = 0; i < ent.size(); i++){
+      for (unsigned int j = 0; j < ent[i]->getNumMeshElements(); j++){
+        MElement *e = ent[i]->getMeshElement(j);
+        elementGroup& gr = allG[it->first];
+        gr.insert(e);  
+      }
+    }
+  }
+  
+  write_MSH2(allG,fname);
+  
+  Msg::Info("done modifying mesh");
+};
\ No newline at end of file

NonLinearSolver/nlsolver/meshModification.h

@@ -85,6 +85,7 @@ class meshModification
     void applyModification(const std::string fname = "",double ver=4.1); 
 };
 
+void write_MSH2(const std::map<int, elementGroup>& gMap,  const std::string filename);
 
 class makePhysicalByBox
 {
@@ -103,10 +104,9 @@ class makePhysicalByBox
           printf("box: xmin = %e xmax = %e y min = %e, ymax = %e\n",xmin, xmax, ymin, ymax);
         }
     };
-    
-    void write_MSH2(const std::map<int, elementGroup>& gMap,  const std::string filename);
     #endif //SWIG
 
+    
   public:
     makePhysicalByBox(){};
     ~makePhysicalByBox(){};
@@ -117,4 +117,22 @@ class makePhysicalByBox
              double ymin, double ymax, int ny);
 };
 
+class QuadToTriangle 
+{
+  #ifndef SWIG
+  protected:
+    GModel* _pModel;
+  
+  protected:
+    void remesh(MElement* ele, std::vector<MElement*>& newEle) const;
+  #endif // SWIG
+    
+  public:
+    QuadToTriangle();
+    ~QuadToTriangle();
+    void loadModel(const std::string meshFileName);
+    void apply(const std::string fname); 
+    
+};
+
 #endif //MESH_MODIFICATION_H_
\ No newline at end of file