diff --git a/CMakeLists.txt b/CMakeLists.txt
index 16393c6634f342eb2a737c93738a76cb5e16da48..6cb12f9dd3ef0976449b883c84cfbb356246930c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -133,7 +133,7 @@ set(GMSH_API
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
Mesh/meshGFaceElliptic.h Mesh/meshPartition.h Mesh/meshGFaceDelaunayInsertion.h
Mesh/simple3D.h Mesh/meshPartitionOptions.h Mesh/directions3D.h Mesh/yamakawa.h
- Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h Mesh/meshMetric.h
+ Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h Mesh/meshMetric.h Mesh/ThinLayer.h
Numeric/mathEvaluator.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/crossConfTerm.h Solver/orthogonalTerm.h
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index d68ec33f65abc395e944a4cd0780eaa6e24b05c2..6a3e985ea7e2c84d34d081718bec9afe45a40315 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -237,6 +237,7 @@
#define ALGO_2D_BAMG 7
#define ALGO_2D_FRONTAL_QUAD 8
#define ALGO_2D_PACK_PRLGRMS 9
+#define ALGO_2D_PACK_PRLGRMS_CSTR 10
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 93b55f9a2d27cae84b9f3e60d3f8fdbc8e95d129..97e21ab0c84daae2135d9746d83c9fa968d479a5 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -74,6 +74,10 @@ class GFace : public GEntity{
GFace(GModel *model, int tag);
virtual ~GFace();
+
+ std::vector<MVertex*> additionalVertices;
+ std::set<MVertex*> constr_vertices;
+
// delete mesh data
virtual void deleteMesh();
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index c1cef52f21b542f243cbe9f3e6d49f77e4430dee..53c2aaca3b37f43c256d3756d416df3a7ba88a1c 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -35,7 +35,7 @@ set(SRC
multiscalePartition.cpp
QuadTriUtils.cpp
QuadTriExtruded2D.cpp QuadTriExtruded3D.cpp QuadTriTransfinite3D.cpp
- simple3D.cpp
+ simple3D.cpp ThinLayer.cpp
DivideAndConquer.cpp
Voronoi3D.cpp
Levy3D.cpp
diff --git a/Mesh/ThinLayer.cpp b/Mesh/ThinLayer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b6749a692db01cb684f23ff064ec53c4339c395
--- /dev/null
+++ b/Mesh/ThinLayer.cpp
@@ -0,0 +1,529 @@
+/*
+ * ThinLayer.cpp
+ *
+ * Created on: Oct 13, 2014
+ * Author: nicolas
+ */
+
+#include "ThinLayer.h"
+#include "GModel.h"
+#include "robustPredicates.h"
+#include "GRegion.h"
+
+
+CorrespVertices::CorrespVertices(){}
+CorrespVertices::~CorrespVertices(){}
+void CorrespVertices::setStartPoint(MVertex* v){
+ this->StartPoint = v;
+}
+void CorrespVertices::setEndPoint(SPoint3 p){
+ this->EndPoint = p;
+}
+void CorrespVertices::setStartNormal(SVector3 v){
+ this->StartNormal = v;
+}
+void CorrespVertices::setEndNormal(SVector3 v){
+ this->EndNormal = v;
+}
+void CorrespVertices::setEndTriangle(faceXtet f){
+ this->EndTriangle = f;
+}
+void CorrespVertices::setdistP2P(double d){
+ this->distP2P = d;
+}
+void CorrespVertices::setangleProd(double a){
+ this->angleProd = a;
+}
+void CorrespVertices::setActive(bool b){
+ this->Active = b;
+}
+void CorrespVertices::setEndTriangleActive(bool b){
+ this->EndTriangleActive = b;
+}
+void CorrespVertices::setIsMaster(bool b){
+ this->IsMaster = b;
+}
+void CorrespVertices::setTagMaster(int i){
+ this->tagMaster = i;
+}
+MVertex* CorrespVertices::getStartPoint(){
+ return StartPoint;
+}
+SPoint3 CorrespVertices::getEndPoint(){
+ return EndPoint;
+}
+SVector3 CorrespVertices::getStartNormal(){
+ return StartNormal;
+}
+SVector3 CorrespVertices::getEndNormal(){
+ return EndNormal;
+}
+faceXtet CorrespVertices::getEndTriangle(){
+ return EndTriangle;
+}
+double CorrespVertices::getdistP2P(){
+ return distP2P;
+}
+double CorrespVertices::getangleProd(){
+ return angleProd;
+}
+bool CorrespVertices::getActive(){
+ return Active;
+}
+bool CorrespVertices::getEndTriangleActive(){
+ return EndTriangleActive;
+}
+bool CorrespVertices::getIsMaster(){
+ return IsMaster;
+}
+int CorrespVertices::getTagMaster(){
+ return tagMaster;
+}
+
+
+
+ThinLayer::ThinLayer(){}
+
+ThinLayer::~ThinLayer(){}
+
+void ThinLayer::perform(){
+ ThinLayer::fillVertexToTets();
+ ThinLayer::fillTetToTet4();
+ std::map<MVertex*,double> AllDist = ThinLayer::computeAllDistToOppSide();
+ ThinLayer::checkOppositeTriangles();
+ ThinLayer::fillvecOfThinSheets();
+ std::set<MVertex*> constr_vertices;
+}
+
+void ThinLayer::checkOppositeTriangles(){
+ //all endTriangle will be set to active or not
+ for (std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it1 = VertexToCorresp.begin();it1 != VertexToCorresp.end();it1++){
+ MVertex* vertTmp = (*it1).first;
+ std::vector<CorrespVertices*> vecCorr = (*it1).second;
+ for (unsigned int i = 0;i < vecCorr.size();i++){
+ CorrespVertices* currentCorr = vecCorr[i];
+ faceXtet currentEndTri = currentCorr->getEndTriangle();
+ MVertex* endP0 = currentEndTri.v[0];
+ MVertex* endP1 = currentEndTri.v[1];
+ MVertex* endP2 = currentEndTri.v[2];
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it2 = VertexToCorresp.find(endP0);
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it3 = VertexToCorresp.find(endP1);
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it4 = VertexToCorresp.find(endP2);
+ (*it1).second[i]->setEndTriangleActive(false);
+ if (it2 != VertexToCorresp.end()){
+ if (it3 != VertexToCorresp.end()){
+ if (it4 != VertexToCorresp.end()){
+ if ((*it2).second[0]->getActive()){
+ if ((*it3).second[0]->getActive()){
+ if ((*it4).second[0]->getActive()){
+ (*it1).second[i]->setEndTriangleActive(true);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+void ThinLayer::fillvecOfThinSheets(){
+ for (std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it1 = VertexToCorresp.begin();it1 != VertexToCorresp.end();it1++){
+ MVertex* vertTmp = (*it1).first;
+ std::vector<CorrespVertices*> vecCorr = (*it1).second;
+ for (unsigned int i = 0;i < vecCorr.size();i++){
+ CorrespVertices* currentCorr = vecCorr[i];
+ if ((currentCorr->getStartPoint()->onWhat()->dim() == 2) && (currentCorr->getActive()) && (currentCorr->getEndTriangleActive()) && (currentCorr->getTagMaster() == (-2))){
+ //Found the first node of a new master sheet
+ std::vector<CorrespVertices*> MasterSheet;
+ MasterSheet.clear();
+ (*it1).second[i]->setTagMaster(-1);
+ faceXtet faceEndSlave = (*it1).second[i]->getEndTriangle();
+ for (unsigned int j = 0;j < 3;j++){
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it2 = VertexToCorresp.find(faceEndSlave.v[j]);
+ if (it2 != VertexToCorresp.end()){
+ if (faceEndSlave.v[j]->onWhat()->dim() == 2){
+ (*it2).second[0]->setTagMaster((*it1).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ MasterSheet.push_back((*it1).second[i]);
+ std::set<MVertex*> CurrentSheet;
+ CurrentSheet.clear();
+ CurrentSheet.insert((*it1).second[i]->getStartPoint());
+ while (CurrentSheet.size() != 0){
+ MVertex* VToDo = (*CurrentSheet.begin());
+ std::vector<MTetrahedron*> surroundingTet = VertexToTets[VToDo];
+ for (unsigned int j = 0;j < surroundingTet.size();j++){
+ for (unsigned int k = 0;k < surroundingTet[j]->getNumVertices();k++){
+ MVertex* ToInsertTmp = surroundingTet[j]->getVertex(k);
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it2 = VertexToCorresp.find(ToInsertTmp);
+ if (ToInsertTmp->onWhat()->tag() == VToDo->onWhat()->tag()){//TODO: OR that onwhat -> dim <, for edges
+ if (it2 != VertexToCorresp.end()){
+ CorrespVertices* correspToInsert = ((*it2).second)[0];
+ if ((correspToInsert->getStartPoint()->onWhat()->dim() == 2) && (correspToInsert->getActive()) && (correspToInsert->getEndTriangleActive()) && (correspToInsert->getTagMaster() == (-2))){
+ MasterSheet.push_back((*it2).second[0]);
+ (*it2).second[0]->setTagMaster(-1);
+ faceXtet faceEndSlave2 = (*it2).second[0]->getEndTriangle();
+ for (unsigned int j = 0;j < 3;j++){
+ std::map<MVertex*,std::vector<CorrespVertices*> >::iterator it3 = VertexToCorresp.find(faceEndSlave2.v[j]);
+ if (it3 != VertexToCorresp.end()){
+ if (faceEndSlave2.v[j]->onWhat()->dim() == 2){
+ (*it3).second[0]->setTagMaster((*it2).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ CurrentSheet.insert(ToInsertTmp);
+ }
+ }
+ }
+ }
+ }
+ CurrentSheet.erase(VToDo);
+ }
+ vecOfThinSheets.push_back(MasterSheet);
+ }
+ }
+ }
+}
+std::map<MVertex*,double> ThinLayer::computeAllDistToOppSide(){
+ std::map<MVertex*,double> AllDistToOppSide;
+ GModel *m = GModel::current();
+// std::vector<MElement*> crackElements;
+ std::set<MVertex*> BoundaryVertices;
+
+
+
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++){
+ MTet4* tet4Tmp = TetToTet4[rTmp->tetrahedra[i]];
+ for (unsigned int j = 0;j < 4;j++){
+ if (tet4Tmp->getNeigh(j) == 0){
+ //find the 4th point,and fill the two vector of the boundary triangle
+ faceXtet fxtTmp(tet4Tmp,j);
+ for (int k = 0;k < 3;k++){
+ MVertex* toTest = fxtTmp.v[k];
+ if (toTest->onWhat()->dim() == 2){
+ if(BoundaryVertices.find(toTest) == BoundaryVertices.end()){
+ BoundaryVertices.insert(toTest);
+ }
+ }
+ }
+ // crackElements.push_back(rTmp->getMeshElement(j));
+ }
+ }
+ }
+ }
+ for(std::set<MVertex*>::iterator it = BoundaryVertices.begin(); it != BoundaryVertices.end(); it++){
+ MVertex* toCompute = (*it);
+ double resultTmp = computeDistToOppSide(toCompute);
+ AllDistToOppSide[toCompute] = resultTmp;
+ }
+
+
+ return AllDistToOppSide;
+}
+
+double ThinLayer::computeDistToOppSide(MVertex* v){
+ double DistToOppSide;
+ //We assume v is on the boundary
+ //First we need to get the internal normal
+ SVector3 InteriorNormal = ThinLayer::computeInteriorNormal(v);
+ //Then we find the first triangle
+ MTet4* FirstTet = ThinLayer::getTetFromPoint(v,InteriorNormal);
+ MTet4* CurrentTet = FirstTet;
+ MTet4* PastTet = FirstTet;
+ SPoint3 CurrentPos = SPoint3(v->x(),v->y(),v->z());
+ SPoint3 LastPos = CurrentPos;
+ int* CurrentTri = 0;
+ CorrespVertices CVTemp;
+ CVTemp.setStartPoint(v);
+ CVTemp.setStartNormal(InteriorNormal);
+ FindNewPoint((&CurrentPos),CurrentTri,CurrentTet,InteriorNormal);
+ faceXtet fxtCV(CurrentTet,(*CurrentTri));
+// while(CurrentTet->getNeigh((*CurrentTri)) != 0){
+ while(CurrentTet != 0){
+ PastTet = CurrentTet;
+ faceXtet fxtCVtmp(PastTet,(*CurrentTri));
+ FindNewPoint((&CurrentPos),CurrentTri,CurrentTet,InteriorNormal);
+ CurrentTet = CurrentTet->getNeigh((*CurrentTri));
+ DistToOppSide += CurrentPos.distance(LastPos);
+ LastPos = CurrentPos;
+ }
+ CVTemp.setEndPoint(LastPos);
+ CVTemp.setEndTriangle(fxtCV);
+ SVector3 EndDir1(fxtCV.v[1]->x() - fxtCV.v[0]->x(),fxtCV.v[1]->y() - fxtCV.v[0]->y(),fxtCV.v[1]->z() - fxtCV.v[0]->z());
+ SVector3 EndDir2(fxtCV.v[2]->x() - fxtCV.v[0]->x(),fxtCV.v[2]->y() - fxtCV.v[0]->y(),fxtCV.v[2]->z() - fxtCV.v[0]->z());
+ SVector3 EndNormal(EndDir1.y() * EndDir2.z() - EndDir1.z() * EndDir2.y(),EndDir1.z() * EndDir2.x() - EndDir1.x() * EndDir2.z(),EndDir1.x() * EndDir2.y() - EndDir1.y() * EndDir2.x());
+ EndNormal.normalize();
+ CVTemp.setEndNormal(EndNormal);
+ CVTemp.setangleProd(fabs(CVTemp.getStartNormal().x() * CVTemp.getEndNormal().x() + CVTemp.getStartNormal().y() * CVTemp.getEndNormal().y() + CVTemp.getStartNormal().z() * CVTemp.getEndNormal().z()));
+ CVTemp.setdistP2P(DistToOppSide);
+ if ((CVTemp.getangleProd() > angleMax) &&(CVTemp.getdistP2P() < distP2PMax)){
+ CVTemp.setActive(true);
+ }
+ else{
+ CVTemp.setActive(false);
+ }
+ CVTemp.setTagMaster(-2);
+ VertexToCorresp[v].push_back(&CVTemp);
+ return DistToOppSide;
+}
+
+SVector3 ThinLayer::computeInteriorNormal(MVertex* v){
+ SPoint3 InteriorNormal(0.0,0.0,0.0);
+ std::vector<MTetrahedron*> currentVecTet = VertexToTets[v];
+ std::vector<SPoint3> vecInteriorNodes;
+ std::vector<SPoint3> vecFirstDir;
+ std::vector<SPoint3> vecSecondDir;
+ for (unsigned int i = 0;i < currentVecTet.size();i++){
+ MTet4* tet4Tmp = TetToTet4[currentVecTet[i]];
+ for (int j = 0;j < 4 ; j++){
+ if (tet4Tmp->getNeigh(j) == 0){
+ //find the 4th point,and fill the two vector of the boundary triangle
+ faceXtet fxtTmp(tet4Tmp,j);
+ for (int k = 0;k < 4;k++){
+ bool foundInteriorPoint = true;
+ for (int l = 0;l < 3;l++){
+ if (tet4Tmp->tet()->getVertex(k) == fxtTmp.v[l]){
+ foundInteriorPoint = false;
+ }
+ }
+ if (foundInteriorPoint){
+ SPoint3 pointTmp(tet4Tmp->tet()->getVertex(k)->x(),tet4Tmp->tet()->getVertex(k)->y(),tet4Tmp->tet()->getVertex(k)->z());
+ vecInteriorNodes.push_back(pointTmp);
+ }
+ }
+ SPoint3 pointTmp1(fxtTmp.v[1]->x() - fxtTmp.v[0]->x(),fxtTmp.v[1]->y() - fxtTmp.v[0]->y(),fxtTmp.v[1]->z() - fxtTmp.v[0]->z());
+ SPoint3 pointTmp2(fxtTmp.v[2]->x() - fxtTmp.v[0]->x(),fxtTmp.v[2]->y() - fxtTmp.v[0]->y(),fxtTmp.v[2]->z() - fxtTmp.v[0]->z());
+ vecFirstDir.push_back(pointTmp1);
+ vecSecondDir.push_back(pointTmp2);
+ }
+ }
+ }
+ //at this point we have all the info necessary.
+ SPoint3 pointInteriorAverage(0.0,0.0,0.0);
+ for (unsigned int i = 0;i < vecInteriorNodes.size();i++){
+ pointInteriorAverage += SPoint3(vecInteriorNodes[i].x(),vecInteriorNodes[i].y(),vecInteriorNodes[i].z());
+// pointInteriorAverage.x() += vecInteriorNodes[i].x();
+// pointInteriorAverage.y() += vecInteriorNodes[i].y();
+// pointInteriorAverage.z() += vecInteriorNodes[i].z();
+ }
+ pointInteriorAverage = SPoint3(pointInteriorAverage.x() / (double(vecInteriorNodes.size())),pointInteriorAverage.y() / (double(vecInteriorNodes.size())),pointInteriorAverage.z() / (double(vecInteriorNodes.size())));
+// pointInteriorAverage.x() = pointInteriorAverage.x() / (double(vecInteriorNodes.size()));
+// pointInteriorAverage.y() = pointInteriorAverage.y() / (double(vecInteriorNodes.size()));
+// pointInteriorAverage.z() = pointInteriorAverage.z() / (double(vecInteriorNodes.size()));
+ SPoint3 dirInteriorAverage(pointInteriorAverage.x() - v->x(),pointInteriorAverage.y() - v->y(),pointInteriorAverage.z() - v->z());
+ double norme = sqrt(dirInteriorAverage.x() * dirInteriorAverage.x() + dirInteriorAverage.y() * dirInteriorAverage.y() + dirInteriorAverage.z() * dirInteriorAverage.z());
+ dirInteriorAverage = SPoint3(dirInteriorAverage.x() / norme,dirInteriorAverage.y() / norme,dirInteriorAverage.z() / norme);
+// dirInteriorAverage.x() = dirInteriorAverage.x() / norme;
+// dirInteriorAverage.y() = dirInteriorAverage.y() / norme;
+// dirInteriorAverage.z() = dirInteriorAverage.z() / norme;
+ std::vector<SPoint3> vecOrthogDir;
+ for(unsigned int i = 0;i < vecFirstDir.size();i++){
+ SPoint3 pointTmp(vecFirstDir[i].y() * vecSecondDir[i].z() - vecFirstDir[i].z() * vecSecondDir[i].y(),vecFirstDir[i].z() * vecSecondDir[i].x() - vecFirstDir[i].x() * vecSecondDir[i].z(),vecFirstDir[i].x() * vecSecondDir[i].y() - vecFirstDir[i].y() * vecSecondDir[i].x());
+ vecOrthogDir.push_back(pointTmp);
+ }
+ for(unsigned int i = 0;i < vecOrthogDir.size();i++){
+ double prodScalTmp = vecOrthogDir[i].x() * dirInteriorAverage.x() + vecOrthogDir[i].y() * dirInteriorAverage.y() + vecOrthogDir[i].z() * dirInteriorAverage.z();
+ if (prodScalTmp < 0.0){
+ vecOrthogDir[i] = SPoint3(- vecOrthogDir[i].x(),- vecOrthogDir[i].y(),- vecOrthogDir[i].z());
+// vecOrthogDir[i].x() = - vecOrthogDir[i].x();
+// vecOrthogDir[i].y() = - vecOrthogDir[i].y();
+// vecOrthogDir[i].z() = - vecOrthogDir[i].z();
+ }
+ double normeTmp = sqrt(vecOrthogDir[i].x() * vecOrthogDir[i].x() + vecOrthogDir[i].y() * vecOrthogDir[i].y() + vecOrthogDir[i].z() * vecOrthogDir[i].z());
+ vecOrthogDir[i] = SPoint3(vecOrthogDir[i].x() / normeTmp,vecOrthogDir[i].y() / normeTmp,vecOrthogDir[i].z() / normeTmp);
+// vecOrthogDir[i].x() = vecOrthogDir[i].x() / normeTmp;
+// vecOrthogDir[i].y() = vecOrthogDir[i].y() / normeTmp;
+// vecOrthogDir[i].z() = vecOrthogDir[i].z() / normeTmp;
+ InteriorNormal += SPoint3(vecOrthogDir[i].x(),vecOrthogDir[i].y(),vecOrthogDir[i].z());
+// InteriorNormal.x() += vecOrthogDir[i].x();
+// InteriorNormal.y() += vecOrthogDir[i].y();
+// InteriorNormal.z() += vecOrthogDir[i].z();
+ }
+ norme = sqrt(InteriorNormal.x() * InteriorNormal.x() + InteriorNormal.y() * InteriorNormal.y() + InteriorNormal.z() * InteriorNormal.z());
+ InteriorNormal = SPoint3(InteriorNormal.x() / norme,InteriorNormal.y() / norme,InteriorNormal.z() / norme);
+// InteriorNormal.x() = InteriorNormal.x() / norme;
+// InteriorNormal.y() = InteriorNormal.y() / norme;
+// InteriorNormal.z() = InteriorNormal.z() / norme;
+ return InteriorNormal;
+}
+
+MTet4* ThinLayer::getTetFromPoint(MVertex* v, SVector3 InteriorNormal){
+ MTet4* TetToGet = 0;
+ std::vector<MTetrahedron*> currentVecTet = VertexToTets[v];
+ for (unsigned int i = 0;i < currentVecTet.size();i++){
+ std::vector<SVector3> vecDir;
+ for (int j = 0;j < 4 ; j++){
+ if (currentVecTet[i]->getVertex(j) != v){
+ SVector3 DirTmp(currentVecTet[i]->getVertex(j)->x() - v->x(),currentVecTet[i]->getVertex(j)->y() - v->y(),currentVecTet[i]->getVertex(j)->z() - v->z());
+ vecDir.push_back(DirTmp);
+ }
+ }
+ bool IsPositiv = ThinLayer::IsPositivOrientation(vecDir[0],vecDir[1],vecDir[2]);
+ if (!IsPositiv){
+ SVector3 DirTmp1 = vecDir[1];
+ SVector3 DirTmp2 = vecDir[0];
+ SVector3 DirTmp3 = vecDir[2];
+ vecDir.clear();
+ vecDir.push_back(DirTmp1);
+ vecDir.push_back(DirTmp2);
+ vecDir.push_back(DirTmp3);
+ }
+ bool isPositiv1 = ThinLayer::IsPositivOrientation(vecDir[0],vecDir[1],InteriorNormal);
+ bool isPositiv2 = ThinLayer::IsPositivOrientation(vecDir[1],vecDir[2],InteriorNormal);
+ bool isPositiv3 = ThinLayer::IsPositivOrientation(vecDir[2],vecDir[0],InteriorNormal);
+ if (isPositiv1){
+ if (isPositiv2){
+ if (isPositiv3){
+ TetToGet = TetToTet4[currentVecTet[i]];
+ }
+ }
+ }
+ }
+ return TetToGet;
+}
+
+bool ThinLayer::IsPositivOrientation(SVector3 a, SVector3 b, SVector3 c){
+ bool result = false;
+ SPoint3 ProdVec(a.y() * b.z() - a.z() * b.y(),a.z() * b.x() - a.x() * b.z(),a.x() * b.y() - a.y() * b.x());
+ double ProdScal = ProdVec.x() * c.x() + ProdVec.y() * c.y() + ProdVec.z() * c.z();
+ if (ProdScal >= 0.0){
+ result = true;
+ }
+ return result;
+}
+
+
+void ThinLayer::FindNewPoint(SPoint3* CurrentPoint, int* CurrentTri, MTet4* CurrentTet, SVector3 InteriorNormal){
+ double distanceP2P = 0.0;
+ double alphaMax;
+ double betaMax;
+ SPoint3 ResultPoint;
+ int triToGet = 0;
+ for (unsigned int n = 0;n < 4 ; n++){
+ //calculer matrice a inverser
+ faceXtet fxt(CurrentTet,n);
+ double a = fxt.v[1]->x() - fxt.v[0]->x();
+ double b = fxt.v[2]->x() - fxt.v[0]->x();
+ double c = InteriorNormal.x();
+ double d = fxt.v[1]->y() - fxt.v[0]->y();
+ double e = fxt.v[2]->y() - fxt.v[0]->y();
+ double f = InteriorNormal.y();
+ double g = fxt.v[1]->z() - fxt.v[0]->z();
+ double h = fxt.v[2]->z() - fxt.v[0]->z();
+ double i = InteriorNormal.z();
+ //produit matrice inverse par vecteur donne poids
+ double detMat = a * e * i + b * f * g + c * d * h - c * e * g - f * h * a - i * b * d;
+ double ai = e * i - f * h;
+ double bi = c * h - b * i;
+ double ci = b * f - c * e;
+ double di = f * g - d * i;
+ double ei = a * i - c * g;
+ double fi = c * d - a * f;
+// double gi = d * h - e * g;
+// double hi = b * g - a * h;
+// double ii = a * e - b * d;
+ double oppx = (*CurrentPoint).x() - fxt.v[0]->x();
+ double oppy = (*CurrentPoint).y() - fxt.v[0]->y();
+ double oppz = (*CurrentPoint).z() - fxt.v[0]->z();
+ double alpha = ai / detMat * oppx + bi / detMat * oppy + ci / detMat * oppz;
+ double beta = di / detMat * oppx + ei / detMat * oppy + fi / detMat * oppz;
+// double gamma = gi / detMat * oppx + hi / detMat * oppy + ii / detMat * oppz;
+ //Test si poids entre 0 et 1 et si length maximale
+ if ((alpha >= (0.0-epsilon)) && (alpha <= (1.0 + epsilon))){
+ if ((beta >= (0.0-epsilon)) && (beta <= (1.0 + epsilon))){
+ if (((1.0 - alpha - beta) >= (0.0-epsilon)) && ((1.0 - alpha - beta) <= (1.0 + epsilon))){
+ SPoint3 PointTmp(fxt.v[0]->x() + alpha * (fxt.v[1]->x() - fxt.v[0]->x()) + beta * (fxt.v[2]->x() - fxt.v[0]->x()),fxt.v[0]->y() + alpha * (fxt.v[1]->y() - fxt.v[0]->y()) + beta * (fxt.v[2]->y() - fxt.v[0]->y()),fxt.v[0]->z() + alpha * (fxt.v[1]->z() - fxt.v[0]->z()) + beta * (fxt.v[2]->z() - fxt.v[0]->z()));
+ double distanceTmp = PointTmp.distance((*CurrentPoint));
+ if (distanceTmp > distanceP2P){
+ distanceP2P = distanceTmp;
+ ResultPoint = PointTmp;
+ triToGet = n;
+ alphaMax = alpha;
+ betaMax = beta;
+ }
+ }
+ }
+ }
+ }
+ //test si trop proche d'un point / une arete
+ if (((alphaMax < epsilon) && (betaMax < epsilon)) || ((alphaMax < epsilon) && ((1.0 - alphaMax - betaMax) < epsilon)) || (((1.0 - alphaMax - betaMax) < epsilon) && (betaMax < epsilon))){
+ //proche d'un point
+ double DistMinTmp = 10000000.0;
+ int indexMinTmp;
+ for (unsigned int i = 0;i < 4;i++){
+ double distanceTmp = sqrt((CurrentTet->tet()->getVertex(i)->x() - ResultPoint.x()) * (CurrentTet->tet()->getVertex(i)->x() - ResultPoint.x()) + (CurrentTet->tet()->getVertex(i)->y() - ResultPoint.y()) * (CurrentTet->tet()->getVertex(i)->y() - ResultPoint.y()) + (CurrentTet->tet()->getVertex(i)->z() - ResultPoint.z()) * (CurrentTet->tet()->getVertex(i)->z() - ResultPoint.z()));
+ if (distanceTmp < DistMinTmp){
+ DistMinTmp = distanceTmp;
+ indexMinTmp = i;
+ }
+ }
+ MTet4* NewTet = ThinLayer::getTetFromPoint(CurrentTet->tet()->getVertex(indexMinTmp),InteriorNormal);
+ SPoint3 PointTmp(CurrentTet->tet()->getVertex(indexMinTmp)->x(),CurrentTet->tet()->getVertex(indexMinTmp)->y(),CurrentTet->tet()->getVertex(indexMinTmp)->z());
+ (*CurrentPoint) = PointTmp;
+ CurrentTet = NewTet;
+ }
+ else if ((alphaMax < epsilon) || (betaMax < epsilon) || ((1.0 - alphaMax - betaMax) < epsilon)){
+ //trop proche d'une arete
+ }
+ else{
+ (*CurrentPoint) = ResultPoint;
+ (*CurrentTri) = triToGet;
+ CurrentTet = CurrentTet->getNeigh(triToGet);
+ }
+}
+
+void ThinLayer::fillVertexToTets(){
+ GModel *m = GModel::current();
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+ MTetrahedron* elem = rTmp->tetrahedra[i];
+ for (unsigned int j = 0; j < 4;j++){
+ std::vector<MTetrahedron*> emptyTetVec;
+ emptyTetVec.clear();
+ VertexToTets[elem->getVertex(j)] = emptyTetVec;
+ std::vector<CorrespVertices*> emptyCVVec;
+ emptyCVVec.clear();
+ VertexToCorresp[elem->getVertex(j)] = emptyCVVec;
+ }
+ }
+ }
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+ MTetrahedron* elem = rTmp->tetrahedra[i];
+ for (unsigned int j = 0; j < 4;j++){
+ VertexToTets[elem->getVertex(j)].push_back(elem);
+ }
+ }
+ }
+}
+
+void ThinLayer::fillTetToTet4(){
+ GModel *m = GModel::current();
+ std::vector<MTet4*> vecAllTet4;
+ vecAllTet4.clear();
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+ MTetrahedron* elem = rTmp->tetrahedra[i];
+ MTet4 tet4tmp = MTet4(elem,0.0);
+ MTet4* currentTet4 = &tet4tmp;
+ TetToTet4[elem] = currentTet4;
+ vecAllTet4.clear();
+ }
+ }
+ connectTets(vecAllTet4);
+}
+
+/****************static declarations****************/
+
+std::map<MVertex*,std::vector<MTetrahedron*> > ThinLayer::VertexToTets;
+std::map<MTetrahedron*,MTet4*> ThinLayer::TetToTet4;
+std::map<MVertex*,std::vector<CorrespVertices*> > ThinLayer::VertexToCorresp;
+std::vector<std::vector<CorrespVertices*> > ThinLayer::vecOfThinSheets;
+
diff --git a/Mesh/ThinLayer.h b/Mesh/ThinLayer.h
new file mode 100644
index 0000000000000000000000000000000000000000..de9e37c819f6429bc122b36ac262ee3ce127e30c
--- /dev/null
+++ b/Mesh/ThinLayer.h
@@ -0,0 +1,137 @@
+/*
+ * ThinLayer.h
+ *
+ * Created on: Oct 13, 2014
+ * Author: nicolas
+ */
+
+#ifndef THINLAYER_H_
+#define THINLAYER_H_
+
+#include "MVertex.h"
+#include "MTriangle.h"
+#include "meshGRegionDelaunayInsertion.h"
+
+
+static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}};
+
+struct faceXtet{
+ MVertex *v[3],*unsorted[3];
+ MTet4 *t1;
+ int i1;
+ faceXtet(MTet4 *_t=0, int iFac=0) : t1(_t), i1(iFac)
+ {
+ MVertex *v0 = t1->tet()->getVertex(faces[iFac][0]);
+ MVertex *v1 = t1->tet()->getVertex(faces[iFac][1]);
+ MVertex *v2 = t1->tet()->getVertex(faces[iFac][2]);
+
+ unsorted[0] = v0;
+ unsorted[1] = v1;
+ unsorted[2] = v2;
+
+ v[0] = std::min(std::min(v0,v1),v2);
+ v[2] = std::max(std::max(v0,v1),v2);
+ v[1] = (v0 != v[0] && v0 != v[2]) ? v0 : (v1 != v[0] && v1 != v[2]) ? v1 : v2;
+ //
+ // std::sort(v, v + 3);
+ }
+
+ inline MVertex * getVertex (int i) const { return unsorted[i];}
+
+ inline bool operator < (const faceXtet & other) const
+ {
+ if (v[0] < other.v[0]) return true;
+ if (v[0] > other.v[0]) return false;
+ if (v[1] < other.v[1]) return true;
+ if (v[1] > other.v[1]) return false;
+ if (v[2] < other.v[2]) return true;
+ return false;
+ }
+ inline bool operator == (const faceXtet & other) const
+ {
+ return (v[0] == other.v[0] &&
+ v[1] == other.v[1] &&
+ v[2] == other.v[2] );
+ }
+ bool visible (MVertex *v){
+ MVertex* v0 = t1->tet()->getVertex(faces[i1][0]);
+ MVertex* v1 = t1->tet()->getVertex(faces[i1][1]);
+ MVertex* v2 = t1->tet()->getVertex(faces[i1][2]);
+ double a[3] = {v0->x(),v0->y(),v0->z()};
+ double b[3] = {v1->x(),v1->y(),v1->z()};
+ double c[3] = {v2->x(),v2->y(),v2->z()};
+ double d[3] = {v->x(),v->y(),v->z()};
+ double o = robustPredicates :: orient3d(a,b,c,d);
+ return o < 0;
+ }
+};
+
+class CorrespVertices{
+private:
+ MVertex* StartPoint;
+ SPoint3 EndPoint;
+ SVector3 StartNormal;
+ SVector3 EndNormal;
+ faceXtet EndTriangle;
+ double distP2P;
+ double angleProd;
+ bool Active;
+ bool EndTriangleActive;
+ bool IsMaster;
+ int tagMaster;
+public:
+ CorrespVertices();
+ ~CorrespVertices();
+ void setStartPoint(MVertex* v);
+ void setEndPoint(SPoint3 p);
+ void setStartNormal(SVector3 v);
+ void setEndNormal(SVector3 v);
+ void setEndTriangle(faceXtet f);
+ void setdistP2P(double d);
+ void setangleProd(double a);
+ void setActive(bool b);
+ void setEndTriangleActive(bool b);
+ void setIsMaster(bool b);
+ void setTagMaster(int i);
+ MVertex* getStartPoint();
+ SPoint3 getEndPoint();
+ SVector3 getStartNormal();
+ SVector3 getEndNormal();
+ faceXtet getEndTriangle();
+ double getdistP2P();
+ double getangleProd();
+ bool getActive();
+ bool getEndTriangleActive();
+ bool getIsMaster();
+ int getTagMaster();
+};
+
+class ThinLayer{
+private:
+public:
+ ThinLayer();
+ ~ThinLayer();
+ static void perform();
+ static void checkOppositeTriangles();
+ static void fillvecOfThinSheets();
+ static std::map<MVertex*,double> computeAllDistToOppSide();
+ static double computeDistToOppSide(MVertex* v);
+ static SVector3 computeInteriorNormal(MVertex* v);
+ static MTet4* getTetFromPoint(MVertex* v, SVector3 InteriorNormal);
+ static bool IsPositivOrientation(SVector3 a, SVector3 b, SVector3 c);
+ static void FindNewPoint(SPoint3* CurrentPoint, int* CurrentTri, MTet4* CurrentTet, SVector3 InteriorNormal);
+ static std::map<MVertex*,std::vector<MTetrahedron*> > VertexToTets;
+ static std::map<MTetrahedron*,MTet4*> TetToTet4;
+ static std::map<MVertex*,std::vector<CorrespVertices*> > VertexToCorresp;
+ static std::vector<std::vector<CorrespVertices*> > vecOfThinSheets;
+ static const double epsilon = 0.00001;
+ static const double angleMax = 0.9;
+ static const double distP2PMax = 0.3;
+ static void fillVertexToTets();
+ static void fillTetToTet4();
+};
+
+
+
+
+#endif /* THINLAYER_H_ */
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 197e197ed7bcc641a836456201f984f6bd291bbb..fdb3717b6d9fbfa4300e7f2b16dd90b3ce013046 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -496,6 +496,7 @@ static bool algoDelaunay2D(GFace *gf)
gf->getMeshingAlgo() == ALGO_2D_FRONTAL ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR ||
gf->getMeshingAlgo() == ALGO_2D_BAMG)
return true;
@@ -1498,7 +1499,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
bool infty = false;
if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
- gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
infty = true;
if (!onlyInitialMesh) {
if (infty)
@@ -1519,6 +1521,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS){
bowyerWatsonParallelograms(gf);
}
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR){
+ bowyerWatsonParallelogramsConstrained(gf,gf->constr_vertices);
+ }
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf);
@@ -2313,7 +2318,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
bool infty = false;
if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
- gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
+ gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
infty = true;
if (infty)
buildBackGroundMesh (gf, &equivalence, ¶metricCoordinates);
@@ -2328,6 +2334,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
bowyerWatsonFrontalLayers(gf,true, &equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
bowyerWatsonParallelograms(gf,&equivalence, ¶metricCoordinates);
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
+ bowyerWatsonParallelogramsConstrained(gf,gf->constr_vertices,&equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf,1000000000, &equivalence, ¶metricCoordinates);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 14cad8fad2278fd4bce86292988023faba3f9cdb..106365a53a1fa78bdc44ab5261b0907d9158bd8a 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1812,6 +1812,113 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
#endif
}
+void bowyerWatsonParallelogramsConstrained(GFace *gf,
+ std::set<MVertex*> constr_vertices,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates)
+{
+ std::cout<<" entered bowyerWatsonParallelogramsConstrained"<<std::endl;
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ bidimMeshData DATA(equivalence, parametricCoordinates);
+ std::vector<MVertex*> packed;
+ std::vector<SMetric3> metrics;
+
+ // printf("creating the points\n");
+ std::cout<<" entering packingOfParallelogramsConstrained"<<std::endl;
+ packingOfParallelogramsConstrained(gf, constr_vertices, packed, metrics);
+ // printf("points created\n");
+ std::cout<<"out of packingOfParallelogramsConstrained"<<std::endl;
+
+ buildMeshGenerationDataStructures (gf, AllTris, DATA);
+
+ std::cout<<"out of buildMeshGenerationDataStructures"<<std::endl;
+ // delaunise the initial mesh
+ int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, DATA);
+ Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
+
+ std::sort(packed.begin(), packed.end(), MVertexLessThanLexicographic());
+ std::cout<<"out of sort"<<std::endl;
+
+ // printf("staring to insert points\n");
+ N_GLOBAL_SEARCH = 0;
+ N_SEARCH = 0;
+ DT_INSERT_VERTEX = 0;
+ // double t1 = Cpu();
+ MTri3 *oneNewTriangle = 0;
+ std::cout<<"entering for packed"<<std::endl;
+ for (unsigned int i=0;i<packed.size();){
+ std::cout<<" First stop for"<<std::endl;
+ MTri3 *worst = *AllTris.begin();
+ std::cout<<" got worst"<<std::endl;
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
+ else{
+ double newPoint[2] ;
+ packed[i]->getParameter(0,newPoint[0]);
+ packed[i]->getParameter(1,newPoint[1]);
+ delete packed[i];
+ double metric[3];
+ // buildMetric(gf, newPoint, metrics[i], metric);
+ buildMetric(gf, newPoint, metric);
+
+ bool success = insertAPoint(gf, AllTris.begin(), newPoint, metric, DATA , AllTris, 0, oneNewTriangle, &oneNewTriangle);
+ if (!success) oneNewTriangle = 0;
+ // if (!success)printf("success %d %d\n",success,AllTris.size());
+ i++;
+ }
+ std::cout<<" out of first if"<<std::endl;
+
+ if(1.0* AllTris.size() > 2.5 * DATA.vSizes.size()){
+ // int n1 = AllTris.size();
+ std::set<MTri3*,compareTri3Ptr>::iterator itd = AllTris.begin();
+ while(itd != AllTris.end()){
+ if((*itd)->isDeleted()){
+ delete *itd;
+ AllTris.erase(itd++);
+ }
+ else
+ itd++;
+ }
+ // Msg::Info("cleaning up the memory %d -> %d", n1, AllTris.size());
+ }
+ std::cout<<" out of second if"<<std::endl;
+
+
+ }
+ std::cout<<"out of for packed"<<std::endl;
+ // printf("%d vertices \n",(int)packed.size());
+ //clock_t t2 = clock();
+ //double DT = (double)(t2-t1)/CLOCKS_PER_SEC;
+ //if (packed.size())printf("points inserted DT %12.5E points per minut : %12.5E %d global searchs %d seachs per insertion\n",DT,60.*packed.size()/DT,N_GLOBAL_SEARCH,N_SEARCH / packed.size());
+ transferDataStructure(gf, AllTris, DATA);
+ std::cout<<"out of transferDataStructure"<<std::endl;
+ std::cout<<"testing all vertices of gf"<<std::endl;
+ for (unsigned int i = 0; i < gf->getNumMeshVertices();i++){
+ MVertex* vtest = gf->getMeshVertex(i);
+ std::cout<<"going to test out parameterisation of the point after pacjing and everything"<<std::endl;
+ double para0, para1;
+ vtest->getParameter(0,para0);
+ vtest->getParameter(1,para1);
+ std::cout<<" point tested: para 1 "<<para0<<" and para 2 "<<para1<<std::endl;
+ }
+ backgroundMesh::unset();
+#if defined(HAVE_ANN)
+ {
+ FieldManager *fields = gf->model()->getFields();
+ BoundaryLayerField *blf = 0;
+ if(fields->getBoundaryLayerField() > 0){
+ Field *bl_field = fields->get(fields->getBoundaryLayerField());
+ blf = dynamic_cast<BoundaryLayerField*> (bl_field);
+ if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
+ }
+ }
+#endif
+ std::cout<<"out of Everything"<<std::endl;
+}
+
static void initialSquare(std::vector<MVertex*> &v,
MVertex *box[4],
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index bd627df025c7c6c6d749fbca89be71b64836e3ee..79a8c389af307c7b0225ec359d9249a0baf9f494 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -157,6 +157,10 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
void bowyerWatsonParallelograms(GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
+void bowyerWatsonParallelogramsConstrained(GFace *gf,
+ std::set<MVertex*> constr_vertices,
+ std::map<MVertex* , MVertex*>* equivalence= 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
void buildBackGroundMesh (GFace *gf,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index efda80006986c55f5a62ffff3aeeb81617778c35..95ffaf3025345d821b515de6d1d2eef89d76cee0 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -17,6 +17,7 @@
#include "CenterlineField.h"
#include <algorithm>
#include "directions3D.h"
+#include "ThinLayer.h"
#include "Context.h"
#include <iostream>
#include <string>
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index e205176e72be10e2360e5807538cbfb0452d9059..d352c5bf98a332bd30183001370d364da9037e7b 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -1036,3 +1036,171 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
// delete rtree;
}
+
+// fills a surface with points in order to build a nice
+// quad mesh ------------
+void packingOfParallelogramsConstrained(GFace* gf, std::set<MVertex*> constr_vertices, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics){
+ //PE MODIF
+// packingOfParallelogramsSmoothness(gf,packed,metrics);
+// return;
+ // END PE MODIF
+#if defined(HAVE_RTREE)
+
+ std::cout<<" inside packingOfParallelogramsConstrained"<<std::endl;
+ const bool goNonLinear = true;
+
+ // FILE *f = Fopen ("parallelograms.pos","w");
+
+ // get all the boundary vertices
+ std::set<MVertex*> bnd_vertices;
+ for(unsigned int i=0;i<gf->getNumMeshElements();i++){
+ MElement* element = gf->getMeshElement(i);
+ for(int j=0;j<element->getNumVertices();j++){
+ MVertex *vertex = element->getVertex(j);
+ if (vertex->onWhat()->dim() < 2)bnd_vertices.insert(vertex);
+ }
+ }
+ std::cout<<" got all the boundary vertices"<<std::endl;
+
+ // put boundary vertices in a fifo queue
+ // std::queue<surfacePointWithExclusionRegion*> fifo;
+ std::set<surfacePointWithExclusionRegion*, compareSurfacePointWithExclusionRegionPtr> fifo;
+ std::vector<surfacePointWithExclusionRegion*> vertices;
+ // put the RTREE
+ RTree<surfacePointWithExclusionRegion*,double,2,double> rtree;
+ SMetric3 metricField(1.0);
+ SPoint2 newp[4][NUMDIR];
+ std::set<MVertex*>::iterator it = bnd_vertices.begin() ;
+
+ char NAME[345]; sprintf(NAME,"crossReal%d.pos",gf->tag());
+ FILE *crossf = Fopen (NAME,"w");
+ if (crossf)fprintf(crossf,"View \"\"{\n");
+ std::cout<<" entering first for"<<std::endl;
+ for (; it != bnd_vertices.end() ; ++it){
+ SPoint2 midpoint;
+ compute4neighbors (gf, *it, midpoint, goNonLinear, newp, metricField,crossf);
+ surfacePointWithExclusionRegion *sp =
+ new surfacePointWithExclusionRegion (*it, newp, midpoint,metricField);
+ // fifo.push(sp);
+ fifo.insert(sp);
+ vertices.push_back(sp);
+ double _min[2],_max[2];
+ sp->minmax(_min,_max);
+ // printf("%g %g .. %g %g\n",_min[0],_min[1],_max[0],_max[1]);
+ rtree.Insert(_min,_max,sp);
+ // sp->print(f);
+ }
+
+ std::set<MVertex*>::iterator it_constr = constr_vertices.begin() ;
+ std::cout<<" entering second for"<<std::endl;
+ for (; it_constr != constr_vertices.end() ; ++it_constr){
+ SPoint2 midpoint;
+ std::cout<<"going to test out parameterisation of the new point"<<std::endl;
+ double para0, para1;
+ (*it_constr)->getParameter(0,para0);
+ (*it_constr)->getParameter(1,para1);
+ std::cout<<" point tested: para 1 "<<para0<<" and para 2 "<<para1<<std::endl;
+ std::cout<<" going to compute4neighbors"<<std::endl;
+ compute4neighbors (gf, *it_constr, midpoint, goNonLinear, newp, metricField,crossf);
+ std::cout<<" going to surfacePointWithExclusionRegion"<<std::endl;
+ surfacePointWithExclusionRegion *sp =
+ new surfacePointWithExclusionRegion (*it_constr, newp, midpoint,metricField);
+ // fifo.push(sp);
+ std::cout<<" done surfacePointWithExclusionRegion"<<std::endl;
+ fifo.insert(sp);
+ vertices.push_back(sp);
+ double _min[2],_max[2];
+ sp->minmax(_min,_max);
+ // printf("%g %g .. %g %g\n",_min[0],_min[1],_max[0],_max[1]);
+ rtree.Insert(_min,_max,sp);
+ // sp->print(f);
+ }
+
+ // printf("initially : %d vertices in the domain\n",vertices.size());
+
+
+ std::cout<<" entering while"<<std::endl;
+ while(!fifo.empty()){
+ //surfacePointWithExclusionRegion & parent = fifo.top();
+ // surfacePointWithExclusionRegion * parent = fifo.front();
+ surfacePointWithExclusionRegion * parent = *fifo.begin();
+ // fifo.pop();
+ fifo.erase(fifo.begin());
+ for (int dir=0;dir<NUMDIR;dir++){
+ // printf("dir = %d\n",dir);
+ int countOK = 0;
+ for (int i=0;i<4;i++){
+ // printf("i = %d %12.5E %12.5E \n",i,parent._p[i][dir].x(),parent._p[i][dir].y());
+
+ // if (!w._tooclose){
+ if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){
+ countOK++;
+ GPoint gp = gf->point(parent->_p[i][dir]);
+ MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
+ std::cout<<"going to test out parameterisation of the new point"<<std::endl;
+ double para0, para1;
+ v->getParameter(0,para0);
+ v->getParameter(1,para1);
+// if (v->getNum() == 341){
+// std::cout<<" 341 dans surface filler !!!!"<<std::endl;
+// std::system("pause");
+// getchar();
+// }
+ std::cout<<" point tested: para 1 "<<para0<<" and para 2 "<<para1<<std::endl;
+ // printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
+ SPoint2 midpoint;
+ compute4neighbors (gf, v, midpoint, goNonLinear, newp, metricField,crossf);
+ surfacePointWithExclusionRegion *sp =
+ new surfacePointWithExclusionRegion (v, newp, midpoint, metricField, parent);
+ // fifo.push(sp);
+ fifo.insert(sp);
+ vertices.push_back(sp);
+ double _min[2],_max[2];
+ sp->minmax(_min,_max);
+ rtree.Insert(_min,_max,sp);
+ }
+ }
+ if (countOK)break;
+ }
+ // printf("%d\n",vertices.size());
+ }
+ std::cout<<" entering if"<<std::endl;
+ if (crossf){
+ fprintf(crossf,"};\n");
+ fclose (crossf);
+ }
+ // printf("done\n");
+
+ // add the vertices as additional vertices in the
+ // surface mesh
+ char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
+ FILE *f = Fopen(ccc,"w");
+ fprintf(f,"View \"\"{\n");
+ std::cout<<" entering another for"<<std::endl;
+ for (unsigned int i=0;i<vertices.size();i++){
+ // if(vertices[i]->_v->onWhat() != gf)
+ vertices[i]->print(f,i);
+ if(vertices[i]->_v->onWhat() == gf) {
+ packed.push_back(vertices[i]->_v);
+ metrics.push_back(vertices[i]->_meshMetric);
+ SPoint2 midpoint;
+ reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
+ std::cout<<"going to test out parameterisation of the REPARAM point"<<std::endl;
+ double para0, para1;
+ vertices[i]->_v->getParameter(0,para0);
+ vertices[i]->_v->getParameter(1,para1);
+ std::cout<<" point tested: para 1 "<<para0<<" and para 2 "<<para1<<std::endl;
+ // fprintf(f,"TP(%22.15E,%22.15E,%g){%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E};\n",vertices[i]->_v->x(),vertices[i]->_v->y(),vertices[i]->_v->z(),
+ // vertices[i]->_meshMetric(0,0),vertices[i]->_meshMetric(0,1),vertices[i]->_meshMetric(0,2),
+ // vertices[i]->_meshMetric(1,0),vertices[i]->_meshMetric(1,1),vertices[i]->_meshMetric(1,2),
+ // vertices[i]->_meshMetric(2,0),vertices[i]->_meshMetric(2,1),vertices[i]->_meshMetric(2,2));
+ //fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",midpoint.x(),midpoint.y(),0.0);
+ }
+ delete vertices[i];
+ }
+ fprintf(f,"};");
+ fclose(f);
+ // printf("packed.size = %d\n",packed.size());
+ // delete rtree;
+#endif
+}
diff --git a/Mesh/surfaceFiller.h b/Mesh/surfaceFiller.h
index dabb87040b6c5eed4b33f161a6ed3dbcd4d06a4a..d317cb01800f40bbee91a2c57394f55c9309ef7d 100644
--- a/Mesh/surfaceFiller.h
+++ b/Mesh/surfaceFiller.h
@@ -8,14 +8,16 @@
#define _SURFACEFILLER_H_
#include <vector>
+#include <set>
#include "STensor3.h"
class GFace;
class MVertex;
-void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed,
- std::vector<SMetric3> &metrics);
-void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed,
- std::vector<SMetric3> &metrics);
+void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
+void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
+void packingOfParallelogramsConstrained(GFace* gf, std::set<MVertex*> constr_vertices, std::vector<MVertex*> &packed, std::vector<SMetric3> &metrics );
+
+
#endif
diff --git a/Plugin/CMakeLists.txt b/Plugin/CMakeLists.txt
index e542921345ffaf8ea0e2c3f5ceadf9c8587a8903..1853783794d5e205f22b7f3e2d00c8a135151733 100644
--- a/Plugin/CMakeLists.txt
+++ b/Plugin/CMakeLists.txt
@@ -32,7 +32,7 @@ set(SRC
Scal2Tens.cpp Scal2Vec.cpp
CutMesh.cpp
NewView.cpp
- SimplePartition.cpp Crack.cpp DuplicateBoundaries.cpp
+ SimplePartition.cpp Crack.cpp DuplicateBoundaries.cpp ThinLayerFixMesh.cpp
FaultZone.cpp
MeshSubEntities.cpp
CVTRemesh.cpp
diff --git a/Plugin/DuplicateBoundaries.cpp b/Plugin/DuplicateBoundaries.cpp
index 04fdb1694376e414cb712839dd5f60345b6538d4..33e31b7f3e1dff4f0ecba697a31620515fb0e37f 100644
--- a/Plugin/DuplicateBoundaries.cpp
+++ b/Plugin/DuplicateBoundaries.cpp
@@ -177,26 +177,8 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
for (unsigned int i = 0; i < fTmp->getNumMeshElements();i++){
MElement* elem = fTmp->getMeshElement(i);
-// std::cout<<"TEST INITIAL"<<std::endl;
- if (elem->getVertex(0)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(0)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (elem->getVertex(1)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(1)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (elem->getVertex(2)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(2)->getNum()<<std::endl;
- while (1){
-
- }
- }
+
}
@@ -205,10 +187,8 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
}
else{
facesBound.push_back(fTmp);
- //ToDuplicateListBoundary.insert(fTmp);
}
}
- std::cout<<"FIN INIT"<<std::endl;
for(unsigned int i=0;i<facesBound.size();i++){
double x = 0.0;
@@ -250,63 +230,34 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
- std::cout<<"case 1"<<std::endl;
-// std::cout<<"on a p1.x() = "<<p1.x()<<" et p2.x() = "<<p2.x()<<" et abs((p2.x()-p1.x()-1.0)) = "<<abs((p2.x()-p1.x()-1.0))<<std::endl;
-// std::cout<<"on a abs(0.5) = "<<abs(0.5)<<" et abs(-0.5) = "<<abs(-0.5)<<" et abs(1.5) = "<<abs(1.5)<<" et abs(-1.5) = "<<abs(-1.5)<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if(abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 2"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 3"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p2.y()-p1.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 4"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 5"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 6"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 7"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 8"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 9"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}
}else if (abs((p1.x()-p2.x()-1.0))<0.0001){
@@ -340,21 +291,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 10"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 11"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 12"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
@@ -368,9 +310,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 13"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
//NOTHING
}
@@ -378,12 +317,10 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
}
}
}
- std::cout<<"on a une taille de facesBound = "<<facesBound.size()<<" et ToDuplicateListBoundary = "<<ToDuplicateListBoundary.size()<<std::endl;
int counterNbDone = 10000;
//int tagEdges = 10000;
for (std::set<GFace*>::iterator itr = ToDuplicateList.begin();itr != ToDuplicateList.end();itr++){
- std::cout<<"on est entre pour la fois numero "<<counterNbDone<<std::endl;
counterNbDone++;
GFace* fTmp = (*itr);
//on doit dupliquer la face, les noeuds, les aretes, puis creer une region entre les deux
@@ -395,10 +332,8 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
GRegion *rTmp = fTmp->getRegion(1);
std::list<GFace*> facesReg;
facesReg = rTmp->faces();
- std::cout<<"rmtp a un nombre INITIAL de faces de "<<rTmp->faces().size()<<std::endl;
for (std::list<GFace*>::iterator it2=facesReg.begin();it2!=facesReg.end();it2++){
GFace* fTemp = (*it2);
- std::cout<<" ID "<<fTemp->tag();
}
double xCenterReg = 0.0;
double yCenterReg = 0.0;
@@ -407,9 +342,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
- std::cout<<"nb de faces de la region 2: "<<rTmp->faces().size()<<std::endl;
- std::cout<<"nb de points de la face 1 de region 2: "<<rTmp->faces().front()->vertices().size()<<std::endl;
- std::cout<<"nb de points de la region 2: "<<rTmp->vertices().size()<<std::endl;
int counterPts = 0;
for(std::list<GFace*>::iterator it2=facesReg.begin();it2!=facesReg.end();it2++){
@@ -427,12 +359,9 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
yCenterReg = yCenterReg/counterPts;
zCenterReg = zCenterReg/counterPts;
- std::cout<<"debut boucle for "<<counterNbDone<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
- std::cout<<"point 1 "<<counterNbDone<<std::endl;
GVertex* vTmp = (*itv);
- std::cout<<"point 2 "<<counterNbDone<<std::endl;
double ponderatedX = 0.99999999999 * vTmp->x() + 0.00000000001 * xCenterReg;
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
@@ -441,42 +370,19 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
// double ponderatedZ = 1.0 * vTmp->z() + 0.0 * zCenterReg;
//GVertex* newv = new gmshVertex(m,vTmp->prescribedMeshSizeAtVertex());
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
- std::cout<<"point 3 "<<counterNbDone<<std::endl;
GVertexAssociation[vTmp] = newv;
- std::cout<<"test du ancien point: "<<vTmp->x()<<" et "<<vTmp->y()<<" et "<<vTmp->z()<<std::endl;
//creation des Gedge correspondantes
- std::cout<<"point 4 "<<counterNbDone<<std::endl;
- std::cout<<"test du nouveau point: "<<newv->x()<<" et "<<newv->y()<<" et "<<newv->z()<<std::endl;
//GEdge* newE = new GEdge(m,tagEdges,vTmp,newv);
//tagEdges++;
//m->add(newE);
- std::cout<<"point gloubigoulba "<<counterNbDone<<std::endl;
GEdge* newE = m->addLine(vTmp,newv);
- std::cout<<"point 5 "<<counterNbDone<<std::endl;
SurroudingEdges[vTmp] = newE;
//maintenant traitement mesh
- std::cout<<"point 6 "<<counterNbDone<<std::endl;
MVertex *vMesh = vTmp->mesh_vertices[0];
- std::cout<<"point 7 "<<counterNbDone<<std::endl;
MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newv);
- std::cout<<"point 8 "<<counterNbDone<<std::endl;
VertexAssociation[vMesh] = newMv;
- std::cout<<"point 9 "<<counterNbDone<<std::endl;
newv->addMeshVertex(newMv);
- if (newMv->getNum() > 100000){
- std::cout<<"on a un getnum de "<<newMv->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (vMesh->getNum() > 100000){
- std::cout<<"What ? on a un getnum vtmp de "<<newMv->getNum()<<std::endl;
- while (1){
-
- }
- }
}
- std::cout<<"apres points "<<counterNbDone<<std::endl;
std::list<GEdge*> elist = fTmp->edges();
std::vector<GEdge*> newEdgesVector;
std::vector<GFace*> SurroundingsFaces;
@@ -513,7 +419,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
newE->addLine(newLine);
}
}
- std::cout<<"apres lignes "<<counterNbDone<<std::endl;
std::vector<std::vector<GEdge*> > VecOfVec;
VecOfVec.push_back(newEdgesVector);
//creation de la nouvelle face
@@ -533,46 +438,9 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
MTriangle *newTri = new MTriangle(firstE,secondE,thirdE);
GFaceAssociation->addTriangle(newTri);
- if (elem->getVertex(0)->getNum() > 100000){
- std::cout<<"ANCIENT TRI on a un getnum de "<<elem->getVertex(0)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (elem->getVertex(1)->getNum() > 100000){
- std::cout<<"ANCIENT TRI on a un getnum de "<<elem->getVertex(1)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (elem->getVertex(2)->getNum() > 100000){
- std::cout<<"ANCIENT TRI on a un getnum de "<<elem->getVertex(2)->getNum()<<std::endl;
- while (1){
-
- }
- }
-
- if (firstE->getNum() > 100000){
- std::cout<<"CREAT TRI on a un getnum de "<<firstE->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (secondE->getNum() > 100000){
- std::cout<<"CREAT TRI on a un getnum de "<<secondE->getNum()<<std::endl;
- while (1){
- }
- }
- if (thirdE->getNum() > 100000){
- std::cout<<"CREAT TRI on a un getnum de "<<thirdE->getNum()<<std::endl;
- while (1){
-
- }
- }
}
- std::cout<<"apres faces "<<counterNbDone<<std::endl;
std::vector<GFace*> VectorFaces;
std::list<GFace*> listFaces;
VectorFaces.push_back(fTmp);
@@ -585,150 +453,66 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
listFaces.push_back(GFaceAssociation);
std::vector<std::vector<GFace*> > VecOfVecGFace;
VecOfVecGFace.push_back(VectorFaces);
- std::cout<<"creation de la nouvelle region "<<counterNbDone<<std::endl;
//creation de la nouvelle region
GRegion* createdRegion = new GRegion(m,counterNbDone);
createdRegion->addPhysicalEntity(PhysicalInterface);
createdRegion->set(listFaces);
m->add(createdRegion);
-// std::cout<<"plop 1"<<std::endl;
std::list<GFace*> RegFaces = rTmp->faces();
-// std::cout<<"plop 2"<<std::endl;
std::list<GFace*>::iterator it = std::find(RegFaces.begin(), RegFaces.end(), fTmp);
-// std::cout<<"plop 3"<<std::endl;
-// std::cout<<"RegFaces a un nombre de faces de "<<RegFaces.size()<<std::endl;
std::cout<<(*it)->tag()<<std::endl;
if(it != RegFaces.end()) RegFaces.erase(it);
-// std::cout<<"RegFaces a un nombre de faces de "<<RegFaces.size()<<std::endl;
-// if(it != RegFaces.end()) RegFaces.remove((*it));
-// std::cout<<"plop 4"<<std::endl;
RegFaces.push_back(GFaceAssociation);
-// std::cout<<GFaceAssociation->tag()<<std::endl;
-// std::cout<<"plop 5"<<std::endl;
rTmp->set(RegFaces);
-// std::cout<<"rmtp a un nombre de faces de "<<rTmp->faces().size()<<std::endl;
std::list<GFace*> NewFacesReg;
NewFacesReg = rTmp->faces();
for (std::list<GFace*>::iterator it2=NewFacesReg.begin();it2!=NewFacesReg.end();it2++){
GFace* fTemp = (*it2);
-// std::cout<<" ID "<<fTemp->tag();
for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
GEdge* eTmp = (*ite);
GEdge* newE = GEdgeAssociation[eTmp];
fTemp->replaceEdge(eTmp,newE);
}
}
- std::cout<<"plop 6"<<std::endl;
for (unsigned int i = 0; i < fTmp->getNumMeshElements();i++){
MElement* elem = fTmp->getMeshElement(i);
-// std::cout<<"on an element with a number of vertices of: "<<elem->getNumVertices()<<std::endl;
- if (elem->getVertex(0)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(0)->getNum()<<std::endl;
- while (1){
+
- }
- }
- if (elem->getVertex(1)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(1)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (elem->getVertex(2)->getNum() > 100000){
- std::cout<<"on a un getnum de "<<elem->getVertex(2)->getNum()<<std::endl;
- while (1){
-
- }
- }
-
-// std::cout<<"on a un element avec un nb de vertex de "<<elem->getNumVertices()<<std::endl;
MVertex *firstE = VertexAssociation.find(elem->getVertex(0))->second;
MVertex *secondE = VertexAssociation.find(elem->getVertex(1))->second;
MVertex *thirdE = VertexAssociation.find(elem->getVertex(2))->second;
- if (firstE->getNum() > 100000){
- std::cout<<"on a un getnum de "<<firstE->getNum()<<std::endl;
- while (1){
+
- }
- }
- if (secondE->getNum() > 100000){
- std::cout<<"on a un getnum de "<<secondE->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (thirdE->getNum() > 100000){
- std::cout<<"on a un getnum de "<<thirdE->getNum()<<std::endl;
- while (1){
-
- }
- }
-// std::cout<<"on a les 3 vertices: "<<elem->getVertex(0)->getNum()<<" "<<elem->getVertex(1)->getNum()<<" "<<elem->getVertex(2)->getNum()<<std::endl;
-// std::cout<<"et les 3 vertices: "<<firstE->getNum()<<" "<<secondE->getNum()<<" "<<thirdE->getNum()<<std::endl;
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(elem->getVertex(0),elem->getVertex(1),elem->getVertex(2),firstE,secondE,thirdE);
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
-// std::cout<<"fin de for"<<std::endl;
}
- std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < rTmp->getNumMeshElements();i++){
MElement* elem = rTmp->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<MVertex*,MVertex* >::iterator itMap = VertexAssociation.find(vert);
if (itMap != VertexAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
- std::cout<<"plop 7"<<std::endl;
-// std::cout<<"rmtp a un nombre de faces de "<<rTmp->faces().size()<<std::endl;
int countTmp = 0;
//maitenant refonte points dans faces
for (std::list<GFace*>::iterator itTmp = NewFacesReg.begin();itTmp != NewFacesReg.end();itTmp++){
GFace* GFaceTmp = (*itTmp);
-// std::cout<<GFaceTmp->tag()<<std::endl;
countTmp++;
-// std::cout<<" iteration "<<countTmp<<std::endl;
for (unsigned int i = 0; i < GFaceTmp->getNumMeshElements();i++){
-// std::cout<<" plop 7.1"<<std::endl;
MElement* elem = GFaceTmp->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
-// std::cout<<" plop 7.2"<<std::endl;
MVertex* vert = elem->getVertex(j);
-// std::cout<<" plop 7.3"<<std::endl;
std::map<MVertex*,MVertex* >::iterator itMap = VertexAssociation.find(vert);
-// std::cout<<" plop 7.4"<<std::endl;
if (itMap != VertexAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
-// if (elem->getNumVertices() != 3){
-// std::cout<<"on a un pb ! nb de vertices de "<<elem->getNumVertices()<<std::endl;
-// while (1){
-//
-// }
-// }
elem->setVertex(j,itMap->second);
- if (elem->getVertex(j)->getNum() > 100000){
- std::cout<<"PENDANT CHANGEMENT on a un getnum de "<<elem->getVertex(j)->getNum()<<std::endl;
- while (1){
-
- }
- }
- if (itMap->second->getNum() > 100000){
- std::cout<<"PENDANT CHANGEMENT on a un ITMAP getnum de "<<itMap->second->getNum()<<std::endl;
- while (1){
-
- }
- }
+
std::vector<MVertex*> FaceVerti = GFaceTmp->mesh_vertices;
@@ -740,20 +524,14 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
}
-// std::cout<<"rmtp a un nombre de faces de "<<rTmp->faces().size()<<std::endl;
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
-// std::cout<<"plop 8"<<std::endl;
}
}
//maintenant on attaque les faces au bord
for (std::set<GFace*>::iterator itr = ToDuplicateListBoundary.begin();itr != ToDuplicateListBoundary.end();itr++){
- std::cout<<"AU BORD on est entre pour la fois numero "<<counterNbDone<<std::endl;
counterNbDone++;
GFace* fTmp = (*itr);
//on doit dupliquer la face, les noeuds, les aretes, puis creer une region entre les deux
@@ -787,48 +565,23 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
yCenterReg = yCenterReg/counterPts;
zCenterReg = zCenterReg/counterPts;
-// std::cout<<"debut boucle for "<<counterNbDone<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
-// std::cout<<"point 1 "<<counterNbDone<<std::endl;
GVertex* vTmp = (*itv);
-// std::cout<<"point 2 "<<counterNbDone<<std::endl;
double ponderatedX = 0.99999999999 * vTmp->x() + 0.00000000001 * xCenterReg;
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
-// double ponderatedX = 1.0 * vTmp->x() + 0.0 * xCenterReg;
-// double ponderatedY = 1.0 * vTmp->y() + 0.0 * yCenterReg;
-// double ponderatedZ = 1.0 * vTmp->z() + 0.0 * zCenterReg;
- //GVertex* newv = new gmshVertex(m,vTmp->prescribedMeshSizeAtVertex());
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
-// std::cout<<"point 3 "<<counterNbDone<<std::endl;
GVertexAssociation[vTmp] = newv;
//creation des Gedge correspondantes
-// std::cout<<"point 4 "<<counterNbDone<<std::endl;
-// std::cout<<"test du nouveau point: "<<newv->x()<<" et "<<newv->y()<<" et "<<newv->z()<<std::endl;
- //GEdge* newE = new GEdge(m,tagEdges,vTmp,newv);
- //tagEdges++;
- //m->add(newE);
GEdge* newE = m->addLine(vTmp,newv);
-// std::cout<<"point 5 "<<counterNbDone<<std::endl;
SurroudingEdges[vTmp] = newE;
//maintenant traitement mesh
-// std::cout<<"point 6 "<<counterNbDone<<std::endl;
MVertex *vMesh = vTmp->mesh_vertices[0];
-// std::cout<<"point 7 "<<counterNbDone<<std::endl;
MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newv);
-// std::cout<<"point 8 "<<counterNbDone<<std::endl;
VertexAssociation[vMesh] = newMv;
-// std::cout<<"point 9 "<<counterNbDone<<std::endl;
newv->addMeshVertex(newMv);
- if (newMv->getNum() > 100000){
- std::cout<<"on a un getnum de "<<newMv->getNum()<<std::endl;
- while (1){
-
- }
- }
}
-// std::cout<<"apres points "<<counterNbDone<<std::endl;
std::list<GEdge*> elist = fTmp->edges();
std::vector<GEdge*> newEdgesVector;
std::vector<GFace*> SurroundingsFaces;
@@ -865,7 +618,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
newE->addLine(newLine);
}
}
-// std::cout<<"apres lignes "<<counterNbDone<<std::endl;
std::vector<std::vector<GEdge*> > VecOfVec;
VecOfVec.push_back(newEdgesVector);
//creation de la nouvelle face
@@ -885,7 +637,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
MTriangle *newTri = new MTriangle(firstE,secondE,thirdE);
GFaceAssociation->addTriangle(newTri);
}
-// std::cout<<"apres faces "<<counterNbDone<<std::endl;
std::vector<GFace*> VectorFaces;
std::list<GFace*> listFaces;
VectorFaces.push_back(fTmp);
@@ -898,7 +649,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
listFaces.push_back(GFaceAssociation);
std::vector<std::vector<GFace*> > VecOfVecGFace;
VecOfVecGFace.push_back(VectorFaces);
-// std::cout<<"creation de la nouvelle region "<<counterNbDone<<std::endl;
//creation de la nouvelle region
GRegion* createdRegion = new GRegion(m,counterNbDone);
@@ -927,16 +677,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
for (unsigned int i = 0; i < fTmp->getNumMeshElements();i++){
MElement* elem = fTmp->getMeshElement(i);
-// std::cout<<"on a un element avec un nb de vertex de "<<elem->getNumVertices()<<std::endl;
MVertex *firstE = VertexAssociation.find(elem->getVertex(0))->second;
MVertex *secondE = VertexAssociation.find(elem->getVertex(1))->second;
MVertex *thirdE = VertexAssociation.find(elem->getVertex(2))->second;
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(elem->getVertex(0),elem->getVertex(1),elem->getVertex(2),firstE,secondE,thirdE);
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
-// std::cout<<"fin de for"<<std::endl;
}
for (unsigned int i = 0; i < rTmp->getNumMeshElements();i++){
MElement* elem = rTmp->getMeshElement(i);
@@ -1004,84 +750,11 @@ PView *GMSH_DuplicateBoundariesPlugin::executeBis(PView *view)
-/*
- std::map<std::pair<MVertex*,GRegion*>,MVertex* > VertexAssociation;
- for (GModel::viter itv = m->firstVertex();itv != m->lastVertex();itv++){
-// std::cout<<"entering first for"<<std::endl;
- GVertex* vTmp = (*itv);
- std::list<GRegion*> r = vTmp->regions();
- std::cout<<"size of mesh_vertices: "<<vTmp->mesh_vertices.size()<<std::endl;
- MVertex *vMesh = vTmp->mesh_vertices[0];
- std::cout<<"size of r: "<<r.size()<<std::endl;
- // for all regions, create association
-// for (int i = 0; i < vTmp->getNumMeshVertices();i++){
-// std::cout<<"entering second for"<<std::endl;
-// MVertex* vMesh = vTmp->getMeshVertex(i);
-// for (GModel::riter itr = vTmp->regions().begin();itr != vTmp->regions().end();itr++){
- for (std::list<GRegion*>::iterator itr = r.begin();itr != r.end();itr++){
- std::cout<<"entering third for"<<std::endl;
- GRegion* vReg = (*itr);
- // std::cout<<"we have the region "<<vReg->
- MVertex *newv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)vReg);
- vReg->mesh_vertices.push_back(newv);
- VertexAssociation[std::make_pair(vMesh,vReg)] = newv;
- std::cout<<"end third for"<<std::endl;
- }
-// std::cout<<"end second for"<<std::endl;
-// }
-// std::cout<<"end first for"<<std::endl;
- }
- std::cout<<"Phase 1 done for DuplicateBoundaries"<<std::endl;
- for (GModel::eiter ite = m->firstEdge();ite != m->lastEdge();ite++){
- GEdge* eTmp = (*ite);
- for (int i = 0; i < eTmp->mesh_vertices.size();i++){
- MVertex* vMesh = eTmp->mesh_vertices[i];
- std::list<GRegion*> r = eTmp->regions();
- for (std::list<GRegion*>::iterator itr = r.begin();itr != r.end();itr++){
- GRegion* vReg = (*itr);
- MVertex *newv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)vReg);
- vReg->mesh_vertices.push_back(newv);
- VertexAssociation[std::make_pair(vMesh,vReg)] = newv;
- }
- }
- }
- std::cout<<"Phase 2 done for DuplicateBoundaries"<<std::endl;
- for (GModel::fiter itf = m->firstFace();itf != m->lastFace();itf++){
- GFace* fTmp = (*itf);
- for (int i = 0; i < fTmp->mesh_vertices.size();i++){
- MVertex* vMesh = fTmp->mesh_vertices[i];
- std::list<GRegion*> r = fTmp->regions();
- for (std::list<GRegion*>::iterator itr = r.begin();itr != r.end();itr++){
- GRegion* vReg = (*itr);
- MVertex *newv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)vReg);
- vReg->mesh_vertices.push_back(newv);
- VertexAssociation[std::make_pair(vMesh,vReg)] = newv;
- }
- }
- }
- std::cout<<"Phase 3 done for DuplicateBoundaries"<<std::endl;
- for (GModel::riter itr = m->firstRegion();itr != m->lastRegion();itr++){
- GRegion* rTmp = (*itr);
- for (int i = 0; i < rTmp->getNumMeshElements();i++){
- MElement* elem = rTmp->getMeshElement(i);
- for (int j = 0;j < elem->getNumVertices();j++){
- MVertex* vert = elem->getVertex(j);
- std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexAssociation.find(std::make_pair(vert,rTmp));
- if (itMap != VertexAssociation.end()){
- elem->setVertex(j,itMap->second);
- }
- }
- }
- }
- */
-// std::cout<<"number for seeds: "<<PhysicalSeeds<<" and for interface: "<<PhysicalInterface<<std::endl;
- std::cout<<"End of DuplicateBoundaries"<<std::endl;
return view;
}
PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
{
- std::cout<<"starting DuplicateBoundaries"<<std::endl;
GModel *m = GModel::current();
m->setFactory("geo");
std::set<GFace*> ToDuplicateList;
@@ -1152,63 +825,34 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 1"<<std::endl;
-// std::cout<<"on a p1.x() = "<<p1.x()<<" et p2.x() = "<<p2.x()<<" et abs((p2.x()-p1.x()-1.0)) = "<<abs((p2.x()-p1.x()-1.0))<<std::endl;
-// std::cout<<"on a abs(0.5) = "<<abs(0.5)<<" et abs(-0.5) = "<<abs(-0.5)<<" et abs(1.5) = "<<abs(1.5)<<" et abs(-1.5) = "<<abs(-1.5)<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if(abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 2"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 3"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p2.y()-p1.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 4"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 5"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 6"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 7"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 8"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 9"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}
}else if (abs((p1.x()-p2.x()-1.0))<0.0001){
@@ -1242,21 +886,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 10"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 11"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 12"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
@@ -1270,9 +905,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 13"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
//NOTHING
}
@@ -1280,7 +912,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
}
}
}
- std::cout<<"on a une taille de facesBound = "<<facesBound.size()<<" et ToDuplicateListBoundary = "<<ToDuplicateListBoundary.size()<<std::endl;
@@ -1295,49 +926,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
std::map<std::pair<GVertex*,GRegion*>,GVertex* > GVertexGlobalAssociation;
std::map<std::pair<GEdge*,GRegion*>,GEdge* > GEdgeGlobalAssociation;
std::map<std::pair<GFace*,GRegion*>,GFace* > GFaceGlobalAssociation;
-// for (GModel::viter itv= m->firstVertex();itv != m->lastVertex();itv++){
-// GVertex* vTmp = (*itv);
-// std::vector<GVertex* > vecTmp;
-// vecTmp.clear();
-// GVertexGlobalAssociation[vTmp] = vecTmp;
-// for (int i = 0; i < vTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = vTmp->mesh_vertices[i];
-// std::vector<MVertex* > vecTmpBis;
-// vecTmpBis.clear();
-// VertexGlobalAssociation[vMesh] = vecTmpBis;
-// }
-// }
-// for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
-// GEdge* eTmp = (*ite);
-// std::vector<GEdge* > vecTmp;
-// vecTmp.clear();
-// GEdgeGlobalAssociation[eTmp] = vecTmp;
-// for (int i = 0; i < eTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = eTmp->mesh_vertices[i];
-// std::vector<MVertex* > vecTmpBis;
-// vecTmpBis.clear();
-// VertexGlobalAssociation[vMesh] = vecTmpBis;
-// }
-// }
-// for (GModel::fiter itf= m->firstFace();itf != m->lastFace();itf++){
-// GFace* fTmp = (*itf);
-// for (int i = 0; i < fTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = fTmp->mesh_vertices[i];
-// std::vector<MVertex* > vecTmpBis;
-// vecTmpBis.clear();
-// VertexGlobalAssociation[vMesh] = vecTmpBis;
-// }
-// }
-// for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
-// GRegion* rTmp = (*itr);
-// for (int i = 0; i < rTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = rTmp->mesh_vertices[i];
-// std::vector<MVertex* > vecTmpBis;
-// vecTmpBis.clear();
-// VertexGlobalAssociation[vMesh] = vecTmpBis;
-// }
-// }
- std::cout<<"entree dans regions first pass"<<std::endl;
+
for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
GRegion* rTmp = (*itr);
std::list<GFace*> RegFaces = rTmp->faces();
@@ -1366,27 +955,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
std::list<GVertex*> vlist;
- std::cout<<"wut 1"<<std::endl;
std::list<GFace*> listFacesTmp = rTmp->faces();
for (std::list<GFace*>::iterator itTp = listFacesTmp.begin();itTp != listFacesTmp.end();itTp++){
- std::cout<<"wut 2"<<std::endl;
std::list<GVertex*> vlist2;
- std::cout<<"wut 3"<<std::endl;
vlist2 = (*itTp)->vertices();
- std::cout<<"wut 4"<<std::endl;
for (std::list<GVertex*>::iterator itTp2 = vlist2.begin();itTp2 != vlist2.end();itTp2++){
- std::cout<<"wut 5"<<std::endl;
if(std::find(vlist.begin(), vlist.end(), *itTp2) == vlist.end())
vlist.push_back(*itTp2);
- std::cout<<"wut 6"<<std::endl;
}
- std::cout<<"wut 7"<<std::endl;
}
- std::cout<<"wut 8"<<std::endl;
- std::cout<<"on a une taille vlist de "<<vlist.size()<<std::endl;
- std::cout<<"init done entering vertices"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
GVertex* vTmp = (*itv);
@@ -1394,7 +973,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
- std::cout<<"inserted correspondance between : "<<vTmp->tag()<<" and : "<<newv->tag()<<std::endl;
GVertexAssociation[vTmp] = newv;
GVertexGlobalAssociation[std::make_pair(vTmp,rTmp)] = newv;
//creation des Gedge correspondantes
@@ -1409,18 +987,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
}
//maintenant on soccupe de duppliquer les edges de la region
std::list<GEdge*> elist = rTmp->edges();
- std::cout<<"on a une taille elist de "<<elist.size()<<std::endl;
- std::cout<<"vertices done entering edges"<<std::endl;
std::vector<GFace*> SurroundingsFaces;
for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
//duplication Gedge
GEdge* eTmp = (*ite);
- std::cout<<"before addline"<<std::endl;
- std::cout<<"begin ID : "<<eTmp->getBeginVertex()->tag()<<" and end ID : "<<eTmp->getEndVertex()->tag()<<std::endl;
- std::cout<<"corresponding begin ID : "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<" and end ID : "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
GEdge* newE = m->addLine(GVertexAssociation[eTmp->getEndVertex()],GVertexAssociation[eTmp->getBeginVertex()]);
- std::cout<<"after addline"<<std::endl;
GEdgeAssociation[eTmp] = newE;
GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
//creation des GFace correspondantes
@@ -1443,7 +1015,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
newE->addMeshVertex(newMv);
}
- std::cout<<"after newvertices"<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
@@ -1451,9 +1022,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
MLine *newLine = new MLine(firstETmp,lastETmp);
newE->addLine(newLine);
}
- std::cout<<"after newMlines"<<std::endl;
}
- std::cout<<"edges done entering faces"<<std::endl;
for (std::list<GFace*>::iterator itf = RegFaces.begin();itf != RegFaces.end();itf++){
GFace* fTmp = (*itf);
std::vector<GEdge*> newEdgesVector;
@@ -1485,9 +1054,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
}
GFaceGlobalAssociation[std::make_pair(fTmp,rTmp)] = GFaceAssociation;
}
- std::cout<<"faces done for regions"<<std::endl;
}
- std::cout<<"Regions done first pass"<<std::endl;
int counterNbDone = 10000;
//maintenant on va traiter les faces initiales
for (std::set<GFace*>::iterator itf = ToDuplicateList.begin();itf != ToDuplicateList.end();itf++){
@@ -1618,25 +1185,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(v1,v2,v3,v4,v5,v6);
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
-// std::cout<<"fin de for"<<std::endl;
}
- std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < reg1->getNumMeshElements();i++){
MElement* elem = reg1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -1646,16 +1205,11 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg2));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
}
- std::cout<<"interior faces done"<<std::endl;
//maintenant on va traiter les faces du bord
for (std::set<GFace*>::iterator itf = ToDuplicateListBoundary.begin();itf != ToDuplicateListBoundary.end();itf++){
counterNbDone++;
@@ -1754,25 +1308,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeFourth(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v3 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(v1,v2,v3,elem->getVertex(0),elem->getVertex(1),elem->getVertex(2));
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
- // std::cout<<"fin de for"<<std::endl;
}
- std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < reg1->getNumMeshElements();i++){
MElement* elem = reg1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -2027,63 +1573,34 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 1"<<std::endl;
-// std::cout<<"on a p1.x() = "<<p1.x()<<" et p2.x() = "<<p2.x()<<" et abs((p2.x()-p1.x()-1.0)) = "<<abs((p2.x()-p1.x()-1.0))<<std::endl;
-// std::cout<<"on a abs(0.5) = "<<abs(0.5)<<" et abs(-0.5) = "<<abs(-0.5)<<" et abs(1.5) = "<<abs(1.5)<<" et abs(-1.5) = "<<abs(-1.5)<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if(abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 2"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 3"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p2.y()-p1.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 4"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 5"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 6"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 7"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 8"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 9"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}
}else if (abs((p1.x()-p2.x()-1.0))<0.0001){
@@ -2117,21 +1634,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (abs((p2.z()-p1.z()))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 10"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 11"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 12"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}
}else if (abs((p1.y()-p2.y()-1.0))<0.0001){
if (abs((p2.z()-p1.z()))<0.0001){
@@ -2145,9 +1653,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (abs((p2.z()-p1.z()-1.0))<0.0001){
pairs.push_back(std::pair<GFace*,GFace*>(facesBound[i],facesBound[j]));
ToDuplicateListBoundary.insert(facesBound[i]);
-// std::cout<<"case 13"<<std::endl;
-// std::cout<<"on a mis dans la liste: x "<<p1.x()<<" y "<<p1.y()<<" z "<<p1.z()<<std::endl;
-// std::cout<<"on avait pour second: x "<<p2.x()<<" y "<<p2.y()<<" z "<<p2.z()<<std::endl;
}else if (abs((p1.z()-p2.z()-1.0))<0.0001){
//NOTHING
}
@@ -2155,7 +1660,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
}
}
}
-// std::cout<<"on a une taille de facesBound = "<<facesBound.size()<<" et ToDuplicateListBoundary = "<<ToDuplicateListBoundary.size()<<std::endl;
@@ -2212,7 +1716,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
// VertexGlobalAssociation[vMesh] = vecTmpBis;
// }
// }
-// std::cout<<"entree dans regions first pass"<<std::endl;
for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
GRegion* rTmp = (*itr);
std::list<GFace*> RegFaces = rTmp->faces();
@@ -2241,27 +1744,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
std::list<GVertex*> vlist;
-// std::cout<<"wut 1"<<std::endl;
std::list<GFace*> listFacesTmp = rTmp->faces();
for (std::list<GFace*>::iterator itTp = listFacesTmp.begin();itTp != listFacesTmp.end();itTp++){
-// std::cout<<"wut 2"<<std::endl;
std::list<GVertex*> vlist2;
-// std::cout<<"wut 3"<<std::endl;
vlist2 = (*itTp)->vertices();
-// std::cout<<"wut 4"<<std::endl;
for (std::list<GVertex*>::iterator itTp2 = vlist2.begin();itTp2 != vlist2.end();itTp2++){
-// std::cout<<"wut 5"<<std::endl;
if(std::find(vlist.begin(), vlist.end(), *itTp2) == vlist.end())
vlist.push_back(*itTp2);
-// std::cout<<"wut 6"<<std::endl;
}
-// std::cout<<"wut 7"<<std::endl;
}
-// std::cout<<"wut 8"<<std::endl;
-// std::cout<<"on a une taille vlist de "<<vlist.size()<<std::endl;
-// std::cout<<"init done entering vertices"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
GVertex* vTmp = (*itv);
@@ -2269,7 +1762,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
-// std::cout<<"inserted correspondance between : "<<vTmp->tag()<<" and : "<<newv->tag()<<std::endl;
GVertexAssociation[vTmp] = newv;
GVertexGlobalAssociation[std::make_pair(vTmp,rTmp)] = newv;
//creation des Gedge correspondantes
@@ -2284,18 +1776,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
}
//maintenant on soccupe de duppliquer les edges de la region
std::list<GEdge*> elist = rTmp->edges();
-// std::cout<<"on a une taille elist de "<<elist.size()<<std::endl;
-// std::cout<<"vertices done entering edges"<<std::endl;
std::vector<GFace*> SurroundingsFaces;
for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
//duplication Gedge
GEdge* eTmp = (*ite);
-// std::cout<<"before addline"<<std::endl;
-// std::cout<<"begin ID : "<<eTmp->getBeginVertex()->tag()<<" and end ID : "<<eTmp->getEndVertex()->tag()<<std::endl;
-// std::cout<<"corresponding begin ID : "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<" and end ID : "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
GEdge* newE = m->addLine(GVertexAssociation[eTmp->getEndVertex()],GVertexAssociation[eTmp->getBeginVertex()]);
-// std::cout<<"after addline"<<std::endl;
GEdgeAssociation[eTmp] = newE;
GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
//creation des GFace correspondantes
@@ -2318,7 +1804,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
newE->addMeshVertex(newMv);
}
-// std::cout<<"after newvertices"<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
@@ -2327,9 +1812,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
MLine3 *newLine = new MLine3(firstETmp,lastETmp,midETmp);
newE->addLine(newLine);
}
-// std::cout<<"after newMlines"<<std::endl;
}
-// std::cout<<"edges done entering faces"<<std::endl;
for (std::list<GFace*>::iterator itf = RegFaces.begin();itf != RegFaces.end();itf++){
GFace* fTmp = (*itf);
std::vector<GEdge*> newEdgesVector;
@@ -2352,7 +1835,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
// nameSurf += ss.str();
// int PhysicalSurfaceTmp = m->setPhysicalName(nameSurf,2);
// GFaceAssociation->addPhysicalEntity(PhysicalSurfaceTmp);
-// std::cout<<"on a associe a la face : "<<fTmp->tag()<<" la face : "<<GFaceAssociation->tag()<<std::endl;
//maintenant traitement mesh
for (unsigned int i = 0; i < fTmp->mesh_vertices.size();i++){
MVertex *vMesh = fTmp->mesh_vertices[i];
@@ -2374,9 +1856,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
}
GFaceGlobalAssociation[std::make_pair(fTmp,rTmp)] = GFaceAssociation;
}
-// std::cout<<"faces done for regions"<<std::endl;
}
-// std::cout<<"Regions done first pass"<<std::endl;
int counterNbDone = 10000;
//maintenant on va traiter les faces initiales
for (std::set<GFace*>::iterator itf = ToDuplicateList.begin();itf != ToDuplicateList.end();itf++){
@@ -2507,12 +1987,9 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(v1,v3,v2,v4,v6,v5);
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
-// std::cout<<"fin de for"<<std::endl;
//second
itMap = VertexGlobalAssociation.find(std::make_pair(elem->getVertex(1),reg1));
@@ -2539,9 +2016,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri2 = new MPrism(v1,v3,v2,v4,v6,v5);
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri2);
//third
@@ -2569,9 +2044,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri3 = new MPrism(v1,v3,v2,v4,v6,v5);
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri3);
//fourth
@@ -2599,24 +2072,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri4 = new MPrism(v1,v3,v2,v4,v6,v5);
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri4);
}
-// std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < reg1->getNumMeshElements();i++){
MElement* elem = reg1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -2626,16 +2092,11 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg2));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
}
-// std::cout<<"interior faces done"<<std::endl;
//maintenant on va traiter les faces du bord
for (std::set<GFace*>::iterator itf = ToDuplicateListBoundary.begin();itf != ToDuplicateListBoundary.end();itf++){
counterNbDone++;
@@ -2734,12 +2195,9 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v3 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(v1,v3,v2,elem->getVertex(0),elem->getVertex(5),elem->getVertex(3));
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
- // std::cout<<"fin de for"<<std::endl;
//second
itMap = VertexGlobalAssociation.find(std::make_pair(elem->getVertex(1),reg1));
if (itMap != VertexGlobalAssociation.end()){
@@ -2753,9 +2211,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v3 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri2 = new MPrism(v1,v3,v2,elem->getVertex(1),elem->getVertex(3),elem->getVertex(4));
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri2);
//third
@@ -2771,9 +2227,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v3 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri3 = new MPrism(v1,v3,v2,elem->getVertex(2),elem->getVertex(4),elem->getVertex(5));
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri3);
//fourth
@@ -2789,24 +2243,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeDuplicate(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v3 = itMap->second;
}
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri4 = new MPrism(v1,v3,v2,elem->getVertex(3),elem->getVertex(5),elem->getVertex(4));
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri4);
}
-// std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < reg1->getNumMeshElements();i++){
MElement* elem = reg1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -3004,21 +2451,8 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
std::vector<std::pair<GFace*,GFace*> > pairs;
std::vector<std::pair<GFace*,GFace*> > newPairs;
int PhysicalInterface = m->setPhysicalName("Interface",2);
-// for (GModel::fiter itf= m->firstFace();itf != m->lastFace();itf++){
-// GFace* ftmp= (*itf);
-// std::list<GEdge*> edgesFace = ftmp->edges();
-//// std::cout<<"size de edgesFace "<<edgesFace.size()<<std::endl;
-// for (std::list<GEdge*>::iterator it2=edgesFace.begin();it2!=edgesFace.end();it2++){
-// GEdge* etmp = (*it2);
-//// std::cout<<"on a la edge "<<etmp->tag()<<std::endl;
-//// std::cout<<"size de face de etmp before : "<<etmp->faces().size()<<std::endl;
-// etmp->addFace(ftmp);
-//// std::cout<<"size de face de etmp after : "<<etmp->faces().size()<<std::endl;
-// }
-// }
for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
GEdge* eTmp = (*ite);
-// std::cout<<"on a size de face de etmp de "<<eTmp->faces().size()<<std::endl;
if (eTmp->faces().size() == 2){
ToDuplicateList.insert(eTmp);
}
@@ -3037,7 +2471,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
std::map<std::pair<GVertex*,GFace*>,GVertex* > GVertexGlobalAssociation;
std::map<std::pair<GEdge*,GFace*>,GEdge* > GEdgeGlobalAssociation;
std::map<std::pair<GFace*,GRegion*>,GFace* > GFaceGlobalAssociation;
- std::cout<<"entree dans regions first pass"<<std::endl;
for (GModel::fiter itr= m->firstFace();itr != m->lastFace();itr++){
GFace* rTmp = (*itr);
std::list<GEdge*> RegEdges = rTmp->edges();
@@ -3066,33 +2499,17 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
std::list<GVertex*> vlist;
-// std::cout<<"wut 1"<<std::endl;
std::list<GEdge*> listedgesTmp = rTmp->edges();
std::cout<<"listeEdge size "<<listedgesTmp.size()<<std::endl;
for (std::list<GEdge*>::iterator itTp = listedgesTmp.begin();itTp != listedgesTmp.end();itTp++){
-// std::cout<<"wut 2"<<std::endl;
std::list<GVertex*> vlist2;
-// std::cout<<"wut 3"<<std::endl;
if(std::find(vlist.begin(), vlist.end(),(*itTp)->getBeginVertex()) == vlist.end())
vlist.push_back((*itTp)->getBeginVertex());
if(std::find(vlist.begin(), vlist.end(),(*itTp)->getEndVertex()) == vlist.end())
vlist.push_back((*itTp)->getEndVertex());
-// vlist2 = (*itTp)->vertices();
-// std::cout<<"vlist2 size "<<vlist2.size()<<std::endl;
-//// std::cout<<"wut 4"<<std::endl;
-// for (std::list<GVertex*>::iterator itTp2 = vlist2.begin();itTp2 != vlist2.end();itTp2++){
-//// std::cout<<"wut 5"<<std::endl;
-// if(std::find(vlist.begin(), vlist.end(), *itTp2) == vlist.end())
-// vlist.push_back(*itTp2);
-//// std::cout<<"wut 6"<<std::endl;
-// }
-// std::cout<<"wut 7"<<std::endl;
}
-// std::cout<<"wut 8"<<std::endl;
-// std::cout<<"on a une taille vlist de "<<vlist.size()<<std::endl;
-// std::cout<<"init done entering vertices"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
GVertex* vTmp = (*itv);
@@ -3100,7 +2517,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
-// std::cout<<"inserted correspondance between : "<<vTmp->tag()<<" and : "<<newv->tag()<<std::endl;
GVertexAssociation[vTmp] = newv;
GVertexGlobalAssociation[std::make_pair(vTmp,rTmp)] = newv;
//creation des Gedge correspondantes
@@ -3115,86 +2531,19 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
}
//maintenant on soccupe de duppliquer les edges de la region
std::list<GEdge*> elist = rTmp->edges();
- std::cout<<"on a une taille elist de "<<elist.size()<<std::endl;
- std::cout<<"vertices done entering edges"<<std::endl;
std::vector<GFace*> SurroundingsFaces;
-// for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
-// //duplication Gedge
-// GEdge* eTmp = (*ite);
-//// std::cout<<"before addline"<<std::endl;
-//// std::cout<<"begin ID : "<<eTmp->getBeginVertex()->tag()<<" and end ID : "<<eTmp->getEndVertex()->tag()<<std::endl;
-//// std::cout<<"corresponding begin ID : "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<" and end ID : "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
-// GEdge* newE = m->addLine(GVertexAssociation[eTmp->getEndVertex()],GVertexAssociation[eTmp->getBeginVertex()]);
-//// std::cout<<"after addline"<<std::endl;
-// GEdgeAssociation[eTmp] = newE;
-// GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
-// //creation des GFace correspondantes
-//// GEdge* firstE = SurroudingEdges.find(eTmp->getEndVertex())->second;
-//// GEdge* lastE = SurroudingEdges.find(eTmp->getBeginVertex())->second;
-//// std::vector<GEdge*> VecEdgesTmp;
-//// VecEdgesTmp.push_back(eTmp);
-//// VecEdgesTmp.push_back(firstE);
-//// VecEdgesTmp.push_back(newE);
-//// VecEdgesTmp.push_back(lastE);
-//// std::vector<std::vector<GEdge*> > VecOfVecTmp;
-//// VecOfVecTmp.push_back(VecEdgesTmp);
-//// GFace* newFaceTmp = m->addPlanarFace(VecOfVecTmp);
-//// SurroundingsFaces.push_back(newFaceTmp);
-// //maintenant traitement mesh
-// for (int i = 0; i < eTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = eTmp->mesh_vertices[i];
-// MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newE);
-// VertexAssociation[vMesh] = newMv;
-// VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// newE->addMeshVertex(newMv);
-// }
-//// std::cout<<"after newvertices"<<std::endl;
-// for (int i = 0; i < eTmp->getNumMeshElements();i++){
-// MElement* elem = eTmp->getMeshElement(i);
-// MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
-// MVertex *lastETmp = VertexAssociation.find(elem->getVertex(1))->second;
-// MVertex *midETmp = VertexAssociation.find(elem->getVertex(2))->second;
-// MLine3 *newLine = new MLine3(firstETmp,lastETmp,midETmp);
-// newE->addLine(newLine);
-// }
-//// std::cout<<"after newMlines"<<std::endl;
-// }
- std::cout<<"edges done entering faces"<<std::endl;
- std::cout<<"RegEdges size : "<<RegEdges.size()<<std::endl;
for (std::list<GEdge*>::iterator itf = RegEdges.begin();itf != RegEdges.end();itf++){
GEdge* eTmp = (*itf);
-// std::vector<GEdge*> newEdgesVector;
-// std::list<GEdge*> elistFace = fTmp->edges();
-// for (std::list<GEdge*>::iterator ite = elistFace.begin();ite != elistFace.end();ite++){
-// GEdge* eTmp = (*ite);
-// GEdge* eToFind = GEdgeAssociation[eTmp];
-// newEdgesVector.push_back(eToFind);
-// }
-// std::vector<std::vector<GEdge*> > VecOfVec;
-// VecOfVec.push_back(newEdgesVector);
//creation de la nouvelle face
- std::cout<<"before newvertices"<<std::endl;
- std::cout<<"beginvertex : "<<eTmp->getBeginVertex()->tag()<<std::endl;
- std::cout<<"endvertex : "<<eTmp->getEndVertex()->tag()<<std::endl;
- std::cout<<"associated begin "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<std::endl;
- std::cout<<"associated end "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
GEdge* newE = m->addLine(GVertexAssociation[eTmp->getBeginVertex()],GVertexAssociation[eTmp->getEndVertex()]);
-// std::cout<<"plop1"<<std::endl;
for (unsigned int i = 0; i < eTmp->mesh_vertices.size();i++){
-// std::cout<<"plop2"<<std::endl;
MVertex *vMesh = eTmp->mesh_vertices[i];
-// std::cout<<"plop3"<<std::endl;
MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newE);
-// std::cout<<"plop4"<<std::endl;
VertexAssociation[vMesh] = newMv;
-// std::cout<<"plop5"<<std::endl;
VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// std::cout<<"plop6"<<std::endl;
newE->addMeshVertex(newMv);
-// std::cout<<"plop7"<<std::endl;
}
- std::cout<<"after newvertices"<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
@@ -3203,51 +2552,15 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
MLine3 *newLine = new MLine3(firstETmp,lastETmp,midETmp);
newE->addLine(newLine);
}
- std::cout<<"after addline"<<std::endl;
GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
-// GEdge* GFaceAssociation = m->addPlanarFace(VecOfVec);
-// int tagTmp = GFaceAssociation->tag();
-// std::ostringstream ss;
-// ss << tagTmp;
-//// char *intStr = itoa(tagTmp);
-// std::string nameSurf = "SURFACE";
-//// nameSurf += std::string(intStr);
-// nameSurf += ss.str();
-// int PhysicalSurfaceTmp = m->setPhysicalName(nameSurf,2);
-// GFaceAssociation->addPhysicalEntity(PhysicalSurfaceTmp);
-// std::cout<<"on a associe a la face : "<<fTmp->tag()<<" la face : "<<GFaceAssociation->tag()<<std::endl;
-// //maintenant traitement mesh
-// for (int i = 0; i < fTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = fTmp->mesh_vertices[i];
-// MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)GFaceAssociation);
-// VertexAssociation[vMesh] = newMv;
-// VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// GFaceAssociation->addMeshVertex(newMv);
-// }
-// for (int i = 0; i < fTmp->getNumMeshElements();i++){
-// MElement* elem = fTmp->getMeshElement(i);
-// MVertex *firstE = VertexAssociation.find(elem->getVertex(0))->second;
-// MVertex *secondE = VertexAssociation.find(elem->getVertex(1))->second;
-// MVertex *thirdE = VertexAssociation.find(elem->getVertex(2))->second;
-// MVertex *fourthE = VertexAssociation.find(elem->getVertex(3))->second;
-// MVertex *fifthE = VertexAssociation.find(elem->getVertex(4))->second;
-// MVertex *sisthE = VertexAssociation.find(elem->getVertex(5))->second;
-// MTriangle6 *newTri = new MTriangle6(firstE,secondE,thirdE,fourthE,fifthE,sisthE);
-// GFaceAssociation->addTriangle(newTri);
-// }
-// GFaceGlobalAssociation[std::make_pair(fTmp,rTmp)] = GFaceAssociation;
}
-// std::cout<<"faces done for regions"<<std::endl;
}
- std::cout<<"Regions done first pass"<<std::endl;
int counterNbDone = 10000;
//maintenant on va traiter les faces initiales
for (std::set<GEdge*>::iterator itf = ToDuplicateList.begin();itf != ToDuplicateList.end();itf++){
counterNbDone++;
- std::cout<<"entree dans toduplicateList numero "<<counterNbDone<<std::endl;
GEdge* eTmp = (*itf);
std::list<GFace*> listFacesTmpT = eTmp->faces();
- std::cout<<"on a une taille de list faces de "<<listFacesTmpT.size()<<std::endl;
std::list<GFace*>::iterator itTmpFace = listFacesTmpT.begin();
GFace* fac1 = (*itTmpFace);
itTmpFace++;
@@ -3258,33 +2571,21 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
vlist.push_back(eTmp->getBeginVertex());
vlist.push_back(eTmp->getEndVertex());
std::map<GVertex*,GEdge* > SurroudingEdges;
- std::cout<<"entre dans vlist.begin"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
GVertex* vTmp = (*itv);
- std::cout<<"test1"<<std::endl;
- std::cout<<"on a fac1 "<<fac1->tag()<<" et fac2 "<<fac2->tag()<<" et vTmp "<<vTmp->tag()<<std::endl;
GVertex* v1=0;
GVertex* v2=0;
std::map<std::pair<GVertex*,GFace*>,GVertex* >::iterator itMap = GVertexGlobalAssociation.find(std::make_pair(vTmp,fac1));
- std::cout<<"test2"<<std::endl;
if (itMap != GVertexGlobalAssociation.end()){
v1 = itMap->second;
- std::cout<<"assigned v1"<<std::endl;
}
- std::cout<<"test3"<<std::endl;
itMap = GVertexGlobalAssociation.find(std::make_pair(vTmp,fac2));
- std::cout<<"test4"<<std::endl;
if (itMap != GVertexGlobalAssociation.end()){
v2 = itMap->second;
- std::cout<<"assigned v2"<<std::endl;
}
- std::cout<<"test5"<<std::endl;
GEdge* newE = m->addLine(v1,v2);
- std::cout<<"test6"<<std::endl;
SurroudingEdges[vTmp] = newE;
- std::cout<<"test7"<<std::endl;
}
- std::cout<<"vertex traites"<<std::endl;
//ici tous les vertex sont traites
//on va traiter les edges
// std::list<GEdge*> elist = fTmp->edges();
@@ -3302,7 +2603,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
if (itMap != GEdgeGlobalAssociation.end()){
e2 = itMap->second;
}
- std::cout<<"creation des GFace correspondantes"<<std::endl;
//creation des GFace correspondantes
GEdge* firstE = SurroudingEdges.find(eTmp->getBeginVertex())->second;
GEdge* lastE = SurroudingEdges.find(eTmp->getEndVertex())->second;
@@ -3336,7 +2636,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
newFaceTmp->addPhysicalEntity(PhysicalInterface);
// createdRegion->set(listFaces);
// m->add(createdRegion);
- std::cout<<"eTmp->getNumMeshElements() "<<eTmp->getNumMeshElements()<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex* v1=0;
@@ -3369,13 +2668,9 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
// if (itMap != VertexGlobalAssociation.end()){
// v6 = itMap->second;
// }
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MQuadrangle *newQua = new MQuadrangle(v1,v2,v3,v4);
- std::cout<<"addition du quad "<<newQua->getNum()<<" avec les vertices "<<v1->getNum()<<" ; "<<v2->getNum()<<" ; "<<v3->getNum()<<" ; "<<v4->getNum()<<std::endl;
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
newFaceTmp->addQuadrangle(newQua);
-// std::cout<<"fin de for"<<std::endl;
//second
itMap = VertexGlobalAssociation.find(std::make_pair(elem->getVertex(2),fac1));
@@ -3402,10 +2697,7 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
// if (itMap != VertexGlobalAssociation.end()){
// v6 = itMap->second;
// }
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MQuadrangle *newQua2 = new MQuadrangle(v1,v2,v3,v4);
- std::cout<<"addition du quad "<<newQua2->getNum()<<" avec les vertices "<<v1->getNum()<<" ; "<<v2->getNum()<<" ; "<<v3->getNum()<<" ; "<<v4->getNum()<<std::endl;
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
newFaceTmp->addQuadrangle(newQua2);
@@ -3413,18 +2705,13 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
}
-// std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < fac1->getNumMeshElements();i++){
MElement* elem = fac1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GFace*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,fac1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -3434,16 +2721,11 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2D(PView *view)
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GFace*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,fac2));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
}
- std::cout<<"interior faces done"<<std::endl;
//maintenant on va traiter les faces du bord
// std::ofstream file("MicrostructurePolycrystal3D.pos");
// file << "View \"test\" {\n";
@@ -3652,21 +2934,9 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
std::vector<std::pair<GFace*,GFace*> > pairs;
std::vector<std::pair<GFace*,GFace*> > newPairs;
int PhysicalInterface = m->setPhysicalName("Interface",2);
-// for (GModel::fiter itf= m->firstFace();itf != m->lastFace();itf++){
-// GFace* ftmp= (*itf);
-// std::list<GEdge*> edgesFace = ftmp->edges();
-//// std::cout<<"size de edgesFace "<<edgesFace.size()<<std::endl;
-// for (std::list<GEdge*>::iterator it2=edgesFace.begin();it2!=edgesFace.end();it2++){
-// GEdge* etmp = (*it2);
-//// std::cout<<"on a la edge "<<etmp->tag()<<std::endl;
-//// std::cout<<"size de face de etmp before : "<<etmp->faces().size()<<std::endl;
-// etmp->addFace(ftmp);
-//// std::cout<<"size de face de etmp after : "<<etmp->faces().size()<<std::endl;
-// }
-// }
+
for (GModel::eiter ite= m->firstEdge();ite != m->lastEdge();ite++){
GEdge* eTmp = (*ite);
-// std::cout<<"on a size de face de etmp de "<<eTmp->faces().size()<<std::endl;
if (eTmp->faces().size() == 2){
ToDuplicateList.insert(eTmp);
}
@@ -3688,7 +2958,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
std::map<std::pair<GVertex*,GFace*>,GVertex* > GVertexGlobalAssociation;
std::map<std::pair<GEdge*,GFace*>,GEdge* > GEdgeGlobalAssociation;
std::map<std::pair<GFace*,GRegion*>,GFace* > GFaceGlobalAssociation;
-// std::cout<<"entree dans regions first pass"<<std::endl;
for (GModel::fiter itr= m->firstFace();itr != m->lastFace();itr++){
GFace* rTmp = (*itr);
std::list<GEdge*> RegEdges = rTmp->edges();
@@ -3717,33 +2986,16 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
std::list<GVertex*> vlist;
-// std::cout<<"wut 1"<<std::endl;
std::list<GEdge*> listedgesTmp = rTmp->edges();
-// std::cout<<"listeEdge size "<<listedgesTmp.size()<<std::endl;
for (std::list<GEdge*>::iterator itTp = listedgesTmp.begin();itTp != listedgesTmp.end();itTp++){
-// std::cout<<"wut 2"<<std::endl;
std::list<GVertex*> vlist2;
-// std::cout<<"wut 3"<<std::endl;
if(std::find(vlist.begin(), vlist.end(),(*itTp)->getBeginVertex()) == vlist.end())
vlist.push_back((*itTp)->getBeginVertex());
if(std::find(vlist.begin(), vlist.end(),(*itTp)->getEndVertex()) == vlist.end())
vlist.push_back((*itTp)->getEndVertex());
-// vlist2 = (*itTp)->vertices();
-// std::cout<<"vlist2 size "<<vlist2.size()<<std::endl;
-//// std::cout<<"wut 4"<<std::endl;
-// for (std::list<GVertex*>::iterator itTp2 = vlist2.begin();itTp2 != vlist2.end();itTp2++){
-//// std::cout<<"wut 5"<<std::endl;
-// if(std::find(vlist.begin(), vlist.end(), *itTp2) == vlist.end())
-// vlist.push_back(*itTp2);
-//// std::cout<<"wut 6"<<std::endl;
-// }
-// std::cout<<"wut 7"<<std::endl;
}
-// std::cout<<"wut 8"<<std::endl;
-// std::cout<<"on a une taille vlist de "<<vlist.size()<<std::endl;
-// std::cout<<"init done entering vertices"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
GVertex* vTmp = (*itv);
@@ -3751,7 +3003,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
-// std::cout<<"inserted correspondance between : "<<vTmp->tag()<<" and : "<<newv->tag()<<std::endl;
GVertexAssociation[vTmp] = newv;
GVertexGlobalAssociation[std::make_pair(vTmp,rTmp)] = newv;
//creation des Gedge correspondantes
@@ -3766,86 +3017,18 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
}
//maintenant on soccupe de duppliquer les edges de la region
std::list<GEdge*> elist = rTmp->edges();
-// std::cout<<"on a une taille elist de "<<elist.size()<<std::endl;
-// std::cout<<"vertices done entering edges"<<std::endl;
std::vector<GFace*> SurroundingsFaces;
-// for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
-// //duplication Gedge
-// GEdge* eTmp = (*ite);
-//// std::cout<<"before addline"<<std::endl;
-//// std::cout<<"begin ID : "<<eTmp->getBeginVertex()->tag()<<" and end ID : "<<eTmp->getEndVertex()->tag()<<std::endl;
-//// std::cout<<"corresponding begin ID : "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<" and end ID : "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
-// GEdge* newE = m->addLine(GVertexAssociation[eTmp->getEndVertex()],GVertexAssociation[eTmp->getBeginVertex()]);
-//// std::cout<<"after addline"<<std::endl;
-// GEdgeAssociation[eTmp] = newE;
-// GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
-// //creation des GFace correspondantes
-//// GEdge* firstE = SurroudingEdges.find(eTmp->getEndVertex())->second;
-//// GEdge* lastE = SurroudingEdges.find(eTmp->getBeginVertex())->second;
-//// std::vector<GEdge*> VecEdgesTmp;
-//// VecEdgesTmp.push_back(eTmp);
-//// VecEdgesTmp.push_back(firstE);
-//// VecEdgesTmp.push_back(newE);
-//// VecEdgesTmp.push_back(lastE);
-//// std::vector<std::vector<GEdge*> > VecOfVecTmp;
-//// VecOfVecTmp.push_back(VecEdgesTmp);
-//// GFace* newFaceTmp = m->addPlanarFace(VecOfVecTmp);
-//// SurroundingsFaces.push_back(newFaceTmp);
-// //maintenant traitement mesh
-// for (int i = 0; i < eTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = eTmp->mesh_vertices[i];
-// MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newE);
-// VertexAssociation[vMesh] = newMv;
-// VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// newE->addMeshVertex(newMv);
-// }
-//// std::cout<<"after newvertices"<<std::endl;
-// for (int i = 0; i < eTmp->getNumMeshElements();i++){
-// MElement* elem = eTmp->getMeshElement(i);
-// MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
-// MVertex *lastETmp = VertexAssociation.find(elem->getVertex(1))->second;
-// MVertex *midETmp = VertexAssociation.find(elem->getVertex(2))->second;
-// MLine3 *newLine = new MLine3(firstETmp,lastETmp,midETmp);
-// newE->addLine(newLine);
-// }
-//// std::cout<<"after newMlines"<<std::endl;
-// }
-// std::cout<<"edges done entering faces"<<std::endl;
-// std::cout<<"RegEdges size : "<<RegEdges.size()<<std::endl;
for (std::list<GEdge*>::iterator itf = RegEdges.begin();itf != RegEdges.end();itf++){
GEdge* eTmp = (*itf);
-// std::vector<GEdge*> newEdgesVector;
-// std::list<GEdge*> elistFace = fTmp->edges();
-// for (std::list<GEdge*>::iterator ite = elistFace.begin();ite != elistFace.end();ite++){
-// GEdge* eTmp = (*ite);
-// GEdge* eToFind = GEdgeAssociation[eTmp];
-// newEdgesVector.push_back(eToFind);
-// }
-// std::vector<std::vector<GEdge*> > VecOfVec;
-// VecOfVec.push_back(newEdgesVector);
- //creation de la nouvelle face
-// std::cout<<"before newvertices"<<std::endl;
-// std::cout<<"beginvertex : "<<eTmp->getBeginVertex()->tag()<<std::endl;
-// std::cout<<"endvertex : "<<eTmp->getEndVertex()->tag()<<std::endl;
-// std::cout<<"associated begin "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<std::endl;
-// std::cout<<"associated end "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
GEdge* newE = m->addLine(GVertexAssociation[eTmp->getBeginVertex()],GVertexAssociation[eTmp->getEndVertex()]);
-// std::cout<<"plop1"<<std::endl;
for (unsigned int i = 0; i < eTmp->mesh_vertices.size();i++){
-// std::cout<<"plop2"<<std::endl;
MVertex *vMesh = eTmp->mesh_vertices[i];
-// std::cout<<"plop3"<<std::endl;
MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)newE);
-// std::cout<<"plop4"<<std::endl;
VertexAssociation[vMesh] = newMv;
-// std::cout<<"plop5"<<std::endl;
VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// std::cout<<"plop6"<<std::endl;
newE->addMeshVertex(newMv);
-// std::cout<<"plop7"<<std::endl;
}
-// std::cout<<"after newvertices"<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
@@ -3854,51 +3037,15 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
MLine3 *newLine = new MLine3(firstETmp,lastETmp,midETmp);
newE->addLine(newLine);
}
-// std::cout<<"after addline"<<std::endl;
GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
-// GEdge* GFaceAssociation = m->addPlanarFace(VecOfVec);
-// int tagTmp = GFaceAssociation->tag();
-// std::ostringstream ss;
-// ss << tagTmp;
-//// char *intStr = itoa(tagTmp);
-// std::string nameSurf = "SURFACE";
-//// nameSurf += std::string(intStr);
-// nameSurf += ss.str();
-// int PhysicalSurfaceTmp = m->setPhysicalName(nameSurf,2);
-// GFaceAssociation->addPhysicalEntity(PhysicalSurfaceTmp);
-// std::cout<<"on a associe a la face : "<<fTmp->tag()<<" la face : "<<GFaceAssociation->tag()<<std::endl;
-// //maintenant traitement mesh
-// for (int i = 0; i < fTmp->mesh_vertices.size();i++){
-// MVertex *vMesh = fTmp->mesh_vertices[i];
-// MVertex *newMv = new MVertex(vMesh->x(), vMesh->y(), vMesh->z(), (GEntity*)GFaceAssociation);
-// VertexAssociation[vMesh] = newMv;
-// VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
-// GFaceAssociation->addMeshVertex(newMv);
-// }
-// for (int i = 0; i < fTmp->getNumMeshElements();i++){
-// MElement* elem = fTmp->getMeshElement(i);
-// MVertex *firstE = VertexAssociation.find(elem->getVertex(0))->second;
-// MVertex *secondE = VertexAssociation.find(elem->getVertex(1))->second;
-// MVertex *thirdE = VertexAssociation.find(elem->getVertex(2))->second;
-// MVertex *fourthE = VertexAssociation.find(elem->getVertex(3))->second;
-// MVertex *fifthE = VertexAssociation.find(elem->getVertex(4))->second;
-// MVertex *sisthE = VertexAssociation.find(elem->getVertex(5))->second;
-// MTriangle6 *newTri = new MTriangle6(firstE,secondE,thirdE,fourthE,fifthE,sisthE);
-// GFaceAssociation->addTriangle(newTri);
-// }
-// GFaceGlobalAssociation[std::make_pair(fTmp,rTmp)] = GFaceAssociation;
}
-// std::cout<<"faces done for regions"<<std::endl;
}
-// std::cout<<"Regions done first pass"<<std::endl;
int counterNbDone = 10000;
//maintenant on va traiter les faces initiales
for (std::set<GEdge*>::iterator itf = ToDuplicateList.begin();itf != ToDuplicateList.end();itf++){
counterNbDone++;
-// std::cout<<"entree dans toduplicateList numero "<<counterNbDone<<std::endl;
GEdge* eTmp = (*itf);
std::list<GFace*> listFacesTmpT = eTmp->faces();
-// std::cout<<"on a une taille de list faces de "<<listFacesTmpT.size()<<std::endl;
std::list<GFace*>::iterator itTmpFace = listFacesTmpT.begin();
GFace* fac1 = (*itTmpFace);
itTmpFace++;
@@ -3909,36 +3056,23 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
vlist.push_back(eTmp->getBeginVertex());
vlist.push_back(eTmp->getEndVertex());
std::map<GVertex*,GEdge* > SurroudingEdges;
-// std::cout<<"entre dans vlist.begin"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
GVertex* vTmp = (*itv);
-// std::cout<<"test1"<<std::endl;
-// std::cout<<"on a fac1 "<<fac1->tag()<<" et fac2 "<<fac2->tag()<<" et vTmp "<<vTmp->tag()<<std::endl;
GVertex* v1=0;
GVertex* v2=0;
std::map<std::pair<GVertex*,GFace*>,GVertex* >::iterator itMap = GVertexGlobalAssociation.find(std::make_pair(vTmp,fac1));
-// std::cout<<"test2"<<std::endl;
if (itMap != GVertexGlobalAssociation.end()){
v1 = itMap->second;
-// std::cout<<"assigned v1"<<std::endl;
}
-// std::cout<<"test3"<<std::endl;
itMap = GVertexGlobalAssociation.find(std::make_pair(vTmp,fac2));
-// std::cout<<"test4"<<std::endl;
if (itMap != GVertexGlobalAssociation.end()){
v2 = itMap->second;
-// std::cout<<"assigned v2"<<std::endl;
}
-// std::cout<<"test5"<<std::endl;
GEdge* newE = m->addLine(v1,v2);
-// std::cout<<"test6"<<std::endl;
SurroudingEdges[vTmp] = newE;
-// std::cout<<"test7"<<std::endl;
}
-// std::cout<<"vertex traites"<<std::endl;
//ici tous les vertex sont traites
//on va traiter les edges
-// std::list<GEdge*> elist = fTmp->edges();
std::vector<GEdge*> newEdgesVector;
std::vector<GFace*> SurroundingsFaces;
//duplication Gedge
@@ -3953,7 +3087,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
if (itMap != GEdgeGlobalAssociation.end()){
e2 = itMap->second;
}
-// std::cout<<"creation des GFace correspondantes"<<std::endl;
//creation des GFace correspondantes
GEdge* firstE = SurroudingEdges.find(eTmp->getBeginVertex())->second;
GEdge* lastE = SurroudingEdges.find(eTmp->getEndVertex())->second;
@@ -3987,7 +3120,6 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
newFaceTmp->addPhysicalEntity(PhysicalInterface);
// createdRegion->set(listFaces);
// m->add(createdRegion);
-// std::cout<<"eTmp->getNumMeshElements() "<<eTmp->getNumMeshElements()<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex* v1=0;
@@ -4020,13 +3152,9 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
// if (itMap != VertexGlobalAssociation.end()){
// v6 = itMap->second;
// }
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MQuadrangle *newQua = new MQuadrangle(v1,v2,v3,v4);
-// std::cout<<"addition du quad "<<newQua->getNum()<<" avec les vertices "<<v1->getNum()<<" ; "<<v2->getNum()<<" ; "<<v3->getNum()<<" ; "<<v4->getNum()<<std::endl;
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
newFaceTmp->addQuadrangle(newQua);
-// std::cout<<"fin de for"<<std::endl;
//second
itMap = VertexGlobalAssociation.find(std::make_pair(elem->getVertex(2),fac1));
@@ -4053,10 +3181,7 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
// if (itMap != VertexGlobalAssociation.end()){
// v6 = itMap->second;
// }
- // std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MQuadrangle *newQua2 = new MQuadrangle(v1,v2,v3,v4);
-// std::cout<<"addition du quad "<<newQua2->getNum()<<" avec les vertices "<<v1->getNum()<<" ; "<<v2->getNum()<<" ; "<<v3->getNum()<<" ; "<<v4->getNum()<<std::endl;
- // std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
newFaceTmp->addQuadrangle(newQua2);
@@ -4064,18 +3189,13 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
}
-// std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < fac1->getNumMeshElements();i++){
MElement* elem = fac1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GFace*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,fac1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -4085,16 +3205,11 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GFace*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,fac2));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
}
-// std::cout<<"interior faces done"<<std::endl;
@@ -4216,18 +3331,13 @@ PView *GMSH_DuplicateBoundariesPlugin::execute2DWithBound(PView *view)
MQuadrangle *newQua2 = new MQuadrangle(v1,v2,elem->getVertex(1),elem->getVertex(2));
newFaceTmp->addQuadrangle(newQua2);
}
-// std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < fac1->getNumMeshElements();i++){
MElement* elem = fac1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GFace*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,fac1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -4497,7 +3607,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
std::map<std::pair<GEdge*,GRegion*>,GEdge* > GEdgeGlobalAssociation;
std::map<std::pair<GFace*,GRegion*>,GFace* > GFaceGlobalAssociation;
- std::cout<<"entree dans regions first pass"<<std::endl;
for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
GRegion* rTmp = (*itr);
std::list<GFace*> RegFaces = rTmp->faces();
@@ -4526,27 +3635,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
std::list<GVertex*> vlist;
- std::cout<<"wut 1"<<std::endl;
std::list<GFace*> listFacesTmp = rTmp->faces();
for (std::list<GFace*>::iterator itTp = listFacesTmp.begin();itTp != listFacesTmp.end();itTp++){
- std::cout<<"wut 2"<<std::endl;
std::list<GVertex*> vlist2;
- std::cout<<"wut 3"<<std::endl;
vlist2 = (*itTp)->vertices();
- std::cout<<"wut 4"<<std::endl;
for (std::list<GVertex*>::iterator itTp2 = vlist2.begin();itTp2 != vlist2.end();itTp2++){
- std::cout<<"wut 5"<<std::endl;
if(std::find(vlist.begin(), vlist.end(), *itTp2) == vlist.end())
vlist.push_back(*itTp2);
- std::cout<<"wut 6"<<std::endl;
}
- std::cout<<"wut 7"<<std::endl;
}
- std::cout<<"wut 8"<<std::endl;
- std::cout<<"on a une taille vlist de "<<vlist.size()<<std::endl;
- std::cout<<"init done entering vertices"<<std::endl;
for (std::list<GVertex*>::iterator itv = vlist.begin();itv != vlist.end();itv++){
//duplication Gvertex
GVertex* vTmp = (*itv);
@@ -4554,7 +3653,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
double ponderatedY = 0.99999999999 * vTmp->y() + 0.00000000001 * yCenterReg;
double ponderatedZ = 0.99999999999 * vTmp->z() + 0.00000000001 * zCenterReg;
GVertex* newv = m->addVertex(ponderatedX,ponderatedY,ponderatedZ,vTmp->prescribedMeshSizeAtVertex());
- std::cout<<"inserted correspondance between : "<<vTmp->tag()<<" and : "<<newv->tag()<<std::endl;
GVertexAssociation[vTmp] = newv;
GVertexGlobalAssociation[std::make_pair(vTmp,rTmp)] = newv;
//creation des Gedge correspondantes
@@ -4569,18 +3667,12 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
}
//maintenant on soccupe de duppliquer les edges de la region
std::list<GEdge*> elist = rTmp->edges();
- std::cout<<"on a une taille elist de "<<elist.size()<<std::endl;
- std::cout<<"vertices done entering edges"<<std::endl;
std::vector<GFace*> SurroundingsFaces;
for (std::list<GEdge*>::iterator ite = elist.begin();ite != elist.end();ite++){
//duplication Gedge
GEdge* eTmp = (*ite);
- std::cout<<"before addline"<<std::endl;
- std::cout<<"begin ID : "<<eTmp->getBeginVertex()->tag()<<" and end ID : "<<eTmp->getEndVertex()->tag()<<std::endl;
- std::cout<<"corresponding begin ID : "<<GVertexAssociation[eTmp->getBeginVertex()]->tag()<<" and end ID : "<<GVertexAssociation[eTmp->getEndVertex()]->tag()<<std::endl;
GEdge* newE = m->addLine(GVertexAssociation[eTmp->getEndVertex()],GVertexAssociation[eTmp->getBeginVertex()]);
- std::cout<<"after addline"<<std::endl;
GEdgeAssociation[eTmp] = newE;
GEdgeGlobalAssociation[std::make_pair(eTmp,rTmp)] = newE;
//creation des GFace correspondantes
@@ -4603,7 +3695,6 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
VertexGlobalAssociation[std::make_pair(vMesh,rTmp)] = newMv;
newE->addMeshVertex(newMv);
}
- std::cout<<"after newvertices"<<std::endl;
for (unsigned int i = 0; i < eTmp->getNumMeshElements();i++){
MElement* elem = eTmp->getMeshElement(i);
MVertex *firstETmp = VertexAssociation.find(elem->getVertex(0))->second;
@@ -4611,9 +3702,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
MLine *newLine = new MLine(firstETmp,lastETmp);
newE->addLine(newLine);
}
- std::cout<<"after newMlines"<<std::endl;
}
- std::cout<<"edges done entering faces"<<std::endl;
for (std::list<GFace*>::iterator itf = RegFaces.begin();itf != RegFaces.end();itf++){
GFace* fTmp = (*itf);
std::vector<GEdge*> newEdgesVector;
@@ -4645,9 +3734,7 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
}
GFaceGlobalAssociation[std::make_pair(fTmp,rTmp)] = GFaceAssociation;
}
- std::cout<<"faces done for regions"<<std::endl;
}
- std::cout<<"Regions done first pass"<<std::endl;
int counterNbDone = 10000;
//maintenant on va traiter les faces initiales
for (std::set<GFace*>::iterator itf = ToDuplicateList.begin();itf != ToDuplicateList.end();itf++){
@@ -4778,25 +3865,17 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
if (itMap != VertexGlobalAssociation.end()){
v6 = itMap->second;
}
-// std::cout<<"creation de prisme "<<counterNbDone<<std::endl;
MPrism *newPri = new MPrism(v1,v2,v3,v4,v5,v6);
-// std::cout<<"addition de prisme "<<counterNbDone<<std::endl;
createdRegion->addPrism(newPri);
-// std::cout<<"fin de for"<<std::endl;
}
- std::cout<<"apres region, refonte points "<<counterNbDone<<std::endl;
for (unsigned int i = 0; i < reg1->getNumMeshElements();i++){
MElement* elem = reg1->getMeshElement(i);
for (int j = 0;j < elem->getNumVertices();j++){
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg1));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
@@ -4806,16 +3885,11 @@ PView *GMSH_DuplicateBoundariesPlugin::executeTer(PView *view)
MVertex* vert = elem->getVertex(j);
std::map<std::pair<MVertex*,GRegion*>,MVertex* >::iterator itMap = VertexGlobalAssociation.find(std::make_pair(vert,reg2));
if (itMap != VertexGlobalAssociation.end()){
-// std::cout<<"on va changer un point"<<std::endl;
elem->setVertex(j,itMap->second);
-// std::cout<<"on a introduit "<<itMap->second->getNum()<<std::endl;
-// std::cout<<"maintenant elem j "<<elem->getVertex(j)->getNum()<<std::endl;
-// std::cout<<"youpi"<<std::endl;
}
}
}
}
- std::cout<<"interior faces done"<<std::endl;
for (std::set<GFace*>::iterator itf = ToDuplicateListBoundary.begin();itf != ToDuplicateListBoundary.end();itf++){
GFace* fTmp = (*itf);
diff --git a/Plugin/PluginManager.cpp b/Plugin/PluginManager.cpp
index 14c1de6efc6fda92034f6981c798f16468bfc22e..c460a099d608cd82747823163cea5e8e7cca8676 100644
--- a/Plugin/PluginManager.cpp
+++ b/Plugin/PluginManager.cpp
@@ -26,6 +26,7 @@
#include "SimplePartition.h"
#include "Crack.h"
#include "DuplicateBoundaries.h"
+#include "ThinLayerFixMesh.h"
#include "HarmonicToTime.h"
#include "ModulusPhase.h"
#include "Integrate.h"
@@ -260,6 +261,8 @@ void PluginManager::registerDefaultPlugins()
("FaultZone", GMSH_RegisterFaultZonePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("DuplicateBoundaries", GMSH_RegisterDuplicateBoundariesPlugin()));
+ allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
+ ("ThinLayerFixMesh", GMSH_RegisterThinLayerFixMeshPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("MeshSubEntities", GMSH_RegisterMeshSubEntitiesPlugin()));
#if defined(HAVE_REVOROPT)
diff --git a/Plugin/ThinLayerFixMesh.cpp b/Plugin/ThinLayerFixMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc13bc63e5ad3da24bb79491d98f9f75d11219a4
--- /dev/null
+++ b/Plugin/ThinLayerFixMesh.cpp
@@ -0,0 +1,1162 @@
+/*
+ * ThinLayerFixMesh.cpp
+ *
+ * Created on: Oct 13, 2014
+ * Author: nicolas
+ */
+
+#include "ThinLayerFixMesh.h"
+#include "GModel.h"
+#include "robustPredicates.h"
+#include "GRegion.h"
+#include "meshGFaceDelaunayInsertion.h"
+//#include "meshGFace.h"
+
+
+
+CorrespVerticesFixMesh::CorrespVerticesFixMesh(){
+// std::cout<<"started init CorrespVerticesFixMesh"<<std::endl;
+// this->EndTriangle = faceXtetFM();
+// this->StartPoint = 0;
+// this->EndPoint = SPoint3(0.0,0.0,0.0);
+// this->StartNormal = SVector3(0.0,0.0,0.0);
+// this->EndNormal = SVector3(0.0,0.0,0.0);
+// this->distP2P = 0.0;
+// this->angleProd = 0.0;
+// this->Active = false;
+// this->EndTriangleActive = false;
+// this->IsMaster = false;
+// this->tagMaster = 0;
+// std::cout<<"completed init CorrespVerticesFixMesh"<<std::endl;
+}
+CorrespVerticesFixMesh::~CorrespVerticesFixMesh(){}
+void CorrespVerticesFixMesh::setStartPoint(MVertex* v){
+ this->StartPoint = v;
+}
+void CorrespVerticesFixMesh::setEndPoint(SPoint3 p){
+ this->EndPoint = p;
+}
+void CorrespVerticesFixMesh::setStartNormal(SVector3 v){
+ this->StartNormal = v;
+}
+void CorrespVerticesFixMesh::setEndNormal(SVector3 v){
+ this->EndNormal = v;
+}
+//void CorrespVerticesFixMesh::setEndTriangle(faceXtetFM f){
+// this->EndTriangle(f.t1,f.i1);
+//}
+void CorrespVerticesFixMesh::setEndTrianglePoint1(MVertex* v){
+ this->EndTrianglePoint1 = v;
+}
+void CorrespVerticesFixMesh::setEndTrianglePoint2(MVertex* v){
+ this->EndTrianglePoint2 = v;
+}
+void CorrespVerticesFixMesh::setEndTrianglePoint3(MVertex* v){
+ this->EndTrianglePoint3 = v;
+}
+void CorrespVerticesFixMesh::setdistP2P(double d){
+ this->distP2P = d;
+}
+void CorrespVerticesFixMesh::setangleProd(double a){
+ this->angleProd = a;
+}
+void CorrespVerticesFixMesh::setActive(bool b){
+ this->Active = b;
+}
+void CorrespVerticesFixMesh::setEndTriangleActive(bool b){
+ this->EndTriangleActive = b;
+}
+void CorrespVerticesFixMesh::setIsMaster(bool b){
+ this->IsMaster = b;
+}
+void CorrespVerticesFixMesh::setTagMaster(int i){
+ this->tagMaster = i;
+}
+MVertex* CorrespVerticesFixMesh::getStartPoint(){
+ return StartPoint;
+}
+SPoint3 CorrespVerticesFixMesh::getEndPoint(){
+ return EndPoint;
+}
+SVector3 CorrespVerticesFixMesh::getStartNormal(){
+ return StartNormal;
+}
+SVector3 CorrespVerticesFixMesh::getEndNormal(){
+ return EndNormal;
+}
+//faceXtetFM CorrespVerticesFixMesh::getEndTriangle(){
+// return EndTriangle;
+//}
+MVertex* CorrespVerticesFixMesh::getEndTrianglePoint1(){
+ return EndTrianglePoint1;
+}
+MVertex* CorrespVerticesFixMesh::getEndTrianglePoint2(){
+ return EndTrianglePoint2;
+}
+MVertex* CorrespVerticesFixMesh::getEndTrianglePoint3(){
+ return EndTrianglePoint3;
+}
+double CorrespVerticesFixMesh::getdistP2P(){
+ return distP2P;
+}
+double CorrespVerticesFixMesh::getangleProd(){
+ return angleProd;
+}
+bool CorrespVerticesFixMesh::getActive(){
+ return Active;
+}
+bool CorrespVerticesFixMesh::getEndTriangleActive(){
+ return EndTriangleActive;
+}
+bool CorrespVerticesFixMesh::getIsMaster(){
+ return IsMaster;
+}
+int CorrespVerticesFixMesh::getTagMaster(){
+ return tagMaster;
+}
+
+StringXNumber ThingLayerFixMeshOptions_Number[] = {
+// {GMSH_FULLRC, "Dimension", NULL, 1.},
+// {GMSH_FULLRC, "PhysicalGroup", NULL, 1.},
+// {GMSH_FULLRC, "OpenBoundaryPhysicalGroup", NULL, 0.},
+};
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterThinLayerFixMeshPlugin()
+ {
+ return new GMSH_ThinLayerFixMeshPlugin();
+ }
+}
+
+//GMSH_ThinLayerFixMeshPlugin::GMSH_ThinLayerFixMeshPlugin(){}
+
+//GMSH_ThinLayerFixMeshPlugin::~GMSH_ThinLayerFixMeshPlugin(){}
+
+std::string GMSH_ThinLayerFixMeshPlugin::getHelp() const
+{
+ return "Fix the mesh in thin parts";
+}
+
+int GMSH_ThinLayerFixMeshPlugin::getNbOptions() const
+{
+ return sizeof(ThingLayerFixMeshOptions_Number) / sizeof(StringXNumber);
+}
+
+StringXNumber *GMSH_ThinLayerFixMeshPlugin::getOption(int iopt)
+{
+ return &ThingLayerFixMeshOptions_Number[iopt];
+}
+
+PView *GMSH_ThinLayerFixMeshPlugin::execute(PView *view)
+{
+ GModel *m = GModel::current();
+ GMSH_ThinLayerFixMeshPlugin::perform();
+// if (m->getDim() == 3){
+// view = GMSH_DuplicateBoundariesPlugin::executeDuplicate(view);
+// }
+// else if (m->getDim() == 2){
+// view = GMSH_DuplicateBoundariesPlugin::execute2DWithBound(view);
+// }
+ return view;
+}
+
+void GMSH_ThinLayerFixMeshPlugin::perform(){
+ VertexToTets.clear();
+ TetToTet4.clear();
+ VertexToCorresp.clear();
+ vecOfThinSheets.clear();
+ GMSH_ThinLayerFixMeshPlugin::fillVertexToTets();
+ GMSH_ThinLayerFixMeshPlugin::fillTetToTet4();
+ //std::cout<<"computeAllDistToOppSide"<<std::endl;
+ std::map<MVertex*,double> AllDist = GMSH_ThinLayerFixMeshPlugin::computeAllDistToOppSide();
+ for (std::map<MVertex*,double>::iterator allDistIt = AllDist.begin();allDistIt != AllDist.end();allDistIt++){
+ //std::cout<<"allDist of point "<<(*allDistIt).first->getNum()<<" with Pos "<<(*allDistIt).first->x()<<" ; "<<(*allDistIt).first->z()<<" ; "<<(*allDistIt).first->y()<<" is "<<(*allDistIt).second<<std::endl;
+ //std::cout<<" Size of vertexToCorresp "<<VertexToCorresp[(*allDistIt).first].size()<<std::endl;
+// //std::cout<<" Testing FaceXTet out of while fourth time "<<VertexToCorresp[(*allDistIt).first][VertexToCorresp[(*allDistIt).first].size() - 1]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ //std::cout<<" Testing StartPoint "<<VertexToCorresp[(*allDistIt).first][VertexToCorresp[(*allDistIt).first].size() - 1]->getStartPoint()->getNum()<<std::endl;
+ //std::cout<<" Testing StartNormal "<<VertexToCorresp[(*allDistIt).first][VertexToCorresp[(*allDistIt).first].size() - 1]->getStartNormal().norm()<<std::endl;
+ }
+ //std::cout<<"checkOppositeTriangles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
+ GMSH_ThinLayerFixMeshPlugin::checkOppositeTriangles();
+ //std::cout<<"fillvecOfThinSheets !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
+ GMSH_ThinLayerFixMeshPlugin::fillvecOfThinSheets();
+ //std::cout<<"Out of fillvecOfThinSheets"<<std::endl;
+// std::set<MVertex*> constr_vertices;
+ for (unsigned int i = 0;i < vecOfThinSheets.size();i++){
+ if (vecOfThinSheets[i].size() > 1){
+ GFace* OnSurf;
+ for (unsigned int j = 0;j < vecOfThinSheets[i].size();j++){
+ //find a point on the surface
+ MVertex* vertOnSurf = vecOfThinSheets[i][j]->getEndTrianglePoint1();
+ if (vertOnSurf->onWhat()->dim() < 2){
+ vertOnSurf = vecOfThinSheets[i][j]->getEndTrianglePoint2();
+ }
+ if (vertOnSurf->onWhat()->dim() < 2){
+ vertOnSurf = vecOfThinSheets[i][j]->getEndTrianglePoint3();
+ }
+ OnSurf = dynamic_cast<GFace*>(vertOnSurf->onWhat());
+ SPoint2 ParOnSurf = OnSurf->parFromPoint(vecOfThinSheets[i][j]->getEndPoint(),1);
+ MVertex* StartPo = vecOfThinSheets[i][j]->getStartPoint();
+ double param1 = 0.0;
+ double param2 = 0.0;
+ StartPo->getParameter(0,param1);
+ StartPo->getParameter(1,param2);
+ //std::cout<<" PointBegin is "<<StartPo->x()<<" ; "<<StartPo->y()<<" ; "<<StartPo->z()<<" with param "<<param1<<" ; "<<param2<<std::endl;
+ //std::cout<<"insertion of point "<<vecOfThinSheets[i][j]->getEndPoint().x()<<" ; "<<vecOfThinSheets[i][j]->getEndPoint().y()<<" ; "<<vecOfThinSheets[i][j]->getEndPoint().z()<<" with param "<<ParOnSurf.x()<<" ; "<<ParOnSurf.y()<<std::endl;
+ MFaceVertex *v = new MFaceVertex(vecOfThinSheets[i][j]->getEndPoint().x(),vecOfThinSheets[i][j]->getEndPoint().y(),vecOfThinSheets[i][j]->getEndPoint().z(),OnSurf,ParOnSurf.x(),ParOnSurf.y());
+ OnSurf->setMeshingAlgo(ALGO_2D_PACK_PRLGRMS_CSTR);
+ OnSurf->constr_vertices.insert(v);
+// OnSurf->addMeshVertex(v);
+// constr_vertices.insert(v);
+ //std::cout<<"inserted point with tag "<<v->getNum()<<" on surface "<<OnSurf->tag()<<std::endl;
+ }
+// OnSurf->setMeshingAlgo(ALGO_2D_PACK_PRLGRMS_CSTR);
+// OnSurf->GFace::deleteMesh();
+// buildBackGroundMesh (OnSurf);
+// meshGFace::modifyInitialMeshForTakingIntoAccountBoundaryLayers(OnSurf);
+// OnSurf->meshStatistics.status = GFace::PENDING;
+// OnSurf->meshStatistics.nbTriangle = OnSurf->meshStatistics.nbEdge = 0;
+// OnSurf->correspondingVertices.clear();
+// std::map<MVertex* , MVertex*>* equivalence;
+// std::map<MVertex*, SPoint2> * parametricCoordinates;
+// bowyerWatsonParallelogramsConstrained(OnSurf,constr_vertices);
+ }
+// constr_vertices.clear();
+ }
+ for (std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it1 = VertexToCorresp.begin();it1 != VertexToCorresp.end();it1++){
+ std::vector<CorrespVerticesFixMesh*> vecCorr = (*it1).second;
+ for (unsigned int i = 0;i < vecCorr.size();i++){
+ delete vecCorr[i];
+ }
+ }
+}
+
+void GMSH_ThinLayerFixMeshPlugin::checkOppositeTriangles(){
+ //all endTriangle will be set to active or not
+ for (std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it1 = VertexToCorresp.begin();it1 != VertexToCorresp.end();it1++){
+// std::cout<<" Entering For"<<std::endl;
+ MVertex* vertTmp = (*it1).first;
+ //std::cout<<" Vert Tested is "<<vertTmp->getNum()<<" and its vector size is "<<(*it1).second.size()<<" pos "<<vertTmp->x()<<" ; "<<vertTmp->y()<<" ; "<<vertTmp->z()<<std::endl;
+ std::vector<CorrespVerticesFixMesh*> vecCorr = (*it1).second;
+ for (unsigned int i = 0;i < vecCorr.size();i++){
+ //std::cout<<" Entering deeper For"<<std::endl;
+ CorrespVerticesFixMesh* currentCorr = vecCorr[i];
+// std::cout<<" Step 1"<<std::endl;
+// faceXtetFM currentEndTri = (*(currentCorr->getEndTriangle()));
+// std::cout<<" Step 2"<<std::endl;
+ MVertex* endP0 = currentCorr->getEndTrianglePoint1();
+// std::cout<<" Step 3"<<std::endl;
+ MVertex* endP1 = currentCorr->getEndTrianglePoint2();
+// std::cout<<" Step 4"<<std::endl;
+ MVertex* endP2 = currentCorr->getEndTrianglePoint3();
+// std::cout<<" Step 5"<<std::endl;
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it2 = VertexToCorresp.find(endP0);
+// std::cout<<" Step 6 ?"<<std::endl;
+// std::cout<<" Tet is "<<currentCorr->getEndTriangle().t1->tet()->getNum()<<std::endl;
+// std::cout<<" Face number is "<<currentCorr->getEndTriangle().i1<<std::endl;
+// std::cout<<" Tet is made of vertex 0 "<<currentCorr->getEndTriangle().t1->tet()->getVertex(0)->getNum()<<std::endl;
+// std::cout<<" Tet is made of vertex 1 "<<currentCorr->getEndTriangle().t1->tet()->getVertex(1)->getNum()<<std::endl;
+// std::cout<<" Tet is made of vertex 2 "<<currentCorr->getEndTriangle().t1->tet()->getVertex(2)->getNum()<<std::endl;
+// std::cout<<" Tet is made of vertex 3 "<<currentCorr->getEndTriangle().t1->tet()->getVertex(3)->getNum()<<std::endl;
+// std::cout<<" Adresses of endP0 "<<endP0<<" and endP1 "<<endP1<<" and endP2 "<<endP2<<std::endl;
+ //std::cout<<" endP0 is "<<endP0->getNum()<<" pos "<<endP0->x()<<" ; "<<endP0->y()<<" ; "<<endP0->z()<<std::endl;
+// std::cout<<" Step 6"<<std::endl;
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it3 = VertexToCorresp.find(endP1);
+ //std::cout<<" endP1 is "<<endP1->getNum()<<" pos "<<endP1->x()<<" ; "<<endP1->y()<<" ; "<<endP1->z()<<std::endl;
+// std::cout<<" Step 7"<<std::endl;
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it4 = VertexToCorresp.find(endP2);
+ //std::cout<<" endP2 is "<<endP2->getNum()<<" pos "<<endP2->x()<<" ; "<<endP2->y()<<" ; "<<endP2->z()<<std::endl;
+// std::cout<<" Step 8"<<std::endl;
+ (*it1).second[i]->setEndTriangleActive(false);
+ //std::cout<<" Starting tests"<<std::endl;
+ bool test0 = false;
+ bool test1 = false;
+ bool test2 = false;
+ if (endP0->onWhat()->dim() < 2){
+ test0 = true;
+ //std::cout<<" test0 true by dim"<<std::endl;
+ }
+ if (endP1->onWhat()->dim() < 2){
+ test1 = true;
+ //std::cout<<" test1 true by dim"<<std::endl;
+ }
+ if (endP2->onWhat()->dim() < 2){
+ test2 = true;
+ //std::cout<<" test2 true by dim"<<std::endl;
+ }
+ if (it2 != VertexToCorresp.end()){
+ if ((*it2).second.size() > 0){
+ if ((*it2).second[0]->getActive()){
+ test0 = true;
+ //std::cout<<" test0 true by active"<<std::endl;
+ }
+ }
+ }
+ if (it3 != VertexToCorresp.end()){
+ if ((*it3).second.size() > 0){
+ if ((*it3).second[0]->getActive()){
+ test1 = true;
+ //std::cout<<" test1 true by active"<<std::endl;
+ }
+ }
+ }
+ if (it4 != VertexToCorresp.end()){
+ if ((*it4).second.size() > 0){
+ if ((*it4).second[0]->getActive()){
+ test2 = true;
+ //std::cout<<" test2 true by active"<<std::endl;
+ }
+ }
+ }
+ if (test0){
+ if (test1){
+ if (test2){
+ (*it1).second[i]->setEndTriangleActive(true);
+ //std::cout<<" EndTriangle Activated"<<std::endl;
+ }
+ }
+ }
+// if (it2 != VertexToCorresp.end()){
+//// std::cout<<" Passed Number 1"<<std::endl;
+// if (it3 != VertexToCorresp.end()){
+//// std::cout<<" Passed Number 2"<<std::endl;
+// if (it4 != VertexToCorresp.end()){
+//// std::cout<<" Passed Number 3"<<std::endl;
+// if (((*it2).second.size() > 0) || (endP0->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 1 Bis"<<std::endl;
+// if (((*it3).second.size() > 0) || (endP1->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 2 Bis"<<std::endl;
+// if (((*it4).second.size() > 0) || (endP2->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 3 Bis"<<std::endl;
+// if (((*it2).second[0]->getActive()) || (endP0->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 4"<<std::endl;
+// if (((*it3).second[0]->getActive()) || (endP1->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 5"<<std::endl;
+// if (((*it4).second[0]->getActive()) || (endP2->onWhat()->dim() < 2)){
+// std::cout<<" Passed Number 6"<<std::endl;
+// (*it1).second[i]->setEndTriangleActive(true);
+// }
+// }
+// }
+// }
+// }
+// }
+// }
+// }
+// }
+ //std::cout<<" Exiting out of deeper For"<<std::endl;
+ }
+// std::cout<<" Getting Out Of For"<<std::endl;
+ }
+}
+
+void GMSH_ThinLayerFixMeshPlugin::fillvecOfThinSheets(){
+ for (std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it1 = VertexToCorresp.begin();it1 != VertexToCorresp.end();it1++){
+ MVertex* vertTmp = (*it1).first;
+ std::vector<CorrespVerticesFixMesh*> vecCorr = (*it1).second;
+ for (unsigned int i = 0;i < vecCorr.size();i++){
+ CorrespVerticesFixMesh* currentCorr = vecCorr[i];
+ //std::cout<<"going to test vecCorr["<<i<<"]"<<" vertex "<<currentCorr->getStartPoint()->getNum()<<" and pos "<<currentCorr->getStartPoint()->x()<<" ; "<<currentCorr->getStartPoint()->y()<<" ; "<<currentCorr->getStartPoint()->z()<<std::endl;
+ if (currentCorr->getStartPoint()->onWhat()->dim() == 2)
+ //std::cout<<"On a surface ; ";
+ if (currentCorr->getActive())
+ //std::cout<<"Is active ; ";
+ if (currentCorr->getEndTriangleActive())
+ //std::cout<<"End Triangle active ; ";
+ if (currentCorr->getTagMaster() == (-2))
+ //std::cout<<"Has yet to be used ; ";
+ //std::cout<<std::endl;
+ if ((currentCorr->getStartPoint()->onWhat()->dim() == 2) && (currentCorr->getActive()) && (currentCorr->getEndTriangleActive()) && (currentCorr->getTagMaster() == (-2))){
+ //Found the first node of a new master sheet
+ //std::cout<<"entering a new master sheet !"<<std::endl;
+ std::vector<CorrespVerticesFixMesh*> MasterSheet;
+ MasterSheet.clear();
+ (*it1).second[i]->setTagMaster(-1);
+// faceXtetFM faceEndSlave = (*((*it1).second[i]->getEndTriangle()));
+// faceXtetFM faceEndSlave = ((*it1).second[i]->getEndTriangle());
+ {
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it2 = VertexToCorresp.find((*it1).second[i]->getEndTrianglePoint1());
+ if (it2 != VertexToCorresp.end()){
+ if ((*it2).second.size() > 0){
+ if ((*it1).second[i]->getEndTrianglePoint1()->onWhat()->dim() == 2){
+ (*it2).second[0]->setTagMaster((*it1).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ it2 = VertexToCorresp.find((*it1).second[i]->getEndTrianglePoint2());
+ if (it2 != VertexToCorresp.end()){
+ if ((*it2).second.size() > 0){
+ if ((*it1).second[i]->getEndTrianglePoint2()->onWhat()->dim() == 2){
+ (*it2).second[0]->setTagMaster((*it1).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ it2 = VertexToCorresp.find((*it1).second[i]->getEndTrianglePoint3());
+ if (it2 != VertexToCorresp.end()){
+ if ((*it2).second.size() > 0){
+ if ((*it1).second[i]->getEndTrianglePoint3()->onWhat()->dim() == 2){
+ (*it2).second[0]->setTagMaster((*it1).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ }
+// for (unsigned int j = 0;j < 3;j++){
+// std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it2 = VertexToCorresp.find(faceEndSlave.v[j]);
+// if (it2 != VertexToCorresp.end()){
+// if (faceEndSlave.v[j]->onWhat()->dim() == 2){
+// (*it2).second[0]->setTagMaster((*it1).second[i]->getStartPoint()->onWhat()->tag());
+// }
+// }
+// }
+ MasterSheet.push_back((*it1).second[i]);
+ std::set<MVertex*> CurrentSheet;
+ CurrentSheet.clear();
+ CurrentSheet.insert((*it1).second[i]->getStartPoint());
+ while (CurrentSheet.size() != 0){
+ MVertex* VToDo = (*CurrentSheet.begin());
+ std::vector<MTetrahedron*> surroundingTet = VertexToTets[VToDo];
+ for (unsigned int j = 0;j < surroundingTet.size();j++){
+ for (unsigned int k = 0;k < surroundingTet[j]->getNumVertices();k++){
+ MVertex* ToInsertTmp = surroundingTet[j]->getVertex(k);
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it2 = VertexToCorresp.find(ToInsertTmp);
+ if (ToInsertTmp->onWhat()->tag() == VToDo->onWhat()->tag()){//TODO: OR that onwhat -> dim <, for edges
+ if (it2 != VertexToCorresp.end()){
+ if ((*it2).second.size() > 0){
+ CorrespVerticesFixMesh* correspToInsert = ((*it2).second)[0];
+ //std::cout<<" Testing "<<((*it2).second)[0]->getStartPoint()->getNum()<<" with ";
+ if (correspToInsert->getStartPoint()->onWhat()->dim() == 2)
+ //std::cout<<"On a surface ; ";
+ if (correspToInsert->getActive())
+ //std::cout<<"Is active ; ";
+ if (correspToInsert->getEndTriangleActive())
+ //std::cout<<"End Triangle active ; ";
+ if (correspToInsert->getTagMaster() == (-2))
+ //std::cout<<"Has yet to be used ; ";
+ //std::cout<<std::endl;
+ if ((correspToInsert->getStartPoint()->onWhat()->dim() == 2) && (correspToInsert->getActive()) && (correspToInsert->getEndTriangleActive()) && (correspToInsert->getTagMaster() == (-2))){
+ MasterSheet.push_back((*it2).second[0]);
+ (*it2).second[0]->setTagMaster(-1);
+ // faceXtetFM faceEndSlave2 = (*((*it2).second[0]->getEndTriangle()));
+ // faceXtetFM faceEndSlave2 = ((*it2).second[0]->getEndTriangle());
+ {
+ std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it3 = VertexToCorresp.find((*it2).second[0]->getEndTrianglePoint1());
+ if (it3 != VertexToCorresp.end()){
+ if ((*it3).second.size() > 0){
+ if ((*it2).second[0]->getEndTrianglePoint1()->onWhat()->dim() == 2){
+ (*it3).second[0]->setTagMaster((*it2).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ it3 = VertexToCorresp.find((*it2).second[0]->getEndTrianglePoint2());
+ if (it3 != VertexToCorresp.end()){
+ if ((*it3).second.size() > 0){
+ if ((*it2).second[0]->getEndTrianglePoint2()->onWhat()->dim() == 2){
+ (*it3).second[0]->setTagMaster((*it2).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ it3 = VertexToCorresp.find((*it2).second[0]->getEndTrianglePoint3());
+ if (it3 != VertexToCorresp.end()){
+ if ((*it3).second.size() > 0){
+ if ((*it2).second[0]->getEndTrianglePoint3()->onWhat()->dim() == 2){
+ (*it3).second[0]->setTagMaster((*it2).second[i]->getStartPoint()->onWhat()->tag());
+ }
+ }
+ }
+ }
+ // for (unsigned int j = 0;j < 3;j++){
+ // std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> >::iterator it3 = VertexToCorresp.find(faceEndSlave2.v[j]);
+ // if (it3 != VertexToCorresp.end()){
+ // if (faceEndSlave2.v[j]->onWhat()->dim() == 2){
+ // (*it3).second[0]->setTagMaster((*it2).second[i]->getStartPoint()->onWhat()->tag());
+ // }
+ // }
+ // }
+ CurrentSheet.insert(ToInsertTmp);
+ }
+ }
+ }
+ }
+ }
+ }
+ CurrentSheet.erase(VToDo);
+ }
+ vecOfThinSheets.push_back(MasterSheet);
+ //std::cout<<"describing new MasterSheet !"<<std::endl;
+ for (unsigned int j = 0;j < MasterSheet.size();j++){
+ //std::cout<<"Number "<<j<<" is "<<MasterSheet[j]->getStartPoint()->getNum()<<" with position "<<MasterSheet[j]->getStartPoint()->x()<<" ; "<<MasterSheet[j]->getStartPoint()->y()<<" ; "<<MasterSheet[j]->getStartPoint()->z()<<std::endl;
+ }
+ //std::cout<<"exiting the master sheet !"<<std::endl;
+ }
+ }
+ }
+}
+std::map<MVertex*,double> GMSH_ThinLayerFixMeshPlugin::computeAllDistToOppSide(){
+ std::map<MVertex*,double> AllDistToOppSide;
+ AllDistToOppSide.clear();
+ GModel *m = GModel::current();
+// std::vector<MElement*> crackElements;
+ std::set<MVertex*> BoundaryVertices;
+ int countdown = 0;
+
+
+
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ //std::cout<<" Entering region "<<rTmp->tag()<<std::endl;
+ for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++){
+ countdown++;
+ //std::cout<<" Entering tet "<<rTmp->tetrahedra[i]->getNum()<<" its the "<<countdown<<"st"<<std::endl;
+ MTet4* tet4Tmp = TetToTet4[rTmp->tetrahedra[i]];
+ //std::cout<<" Neighbours 0: "<<tet4Tmp->getNeigh(0)<<" Neighbours 1: "<<tet4Tmp->getNeigh(1)<<" Neighbours 2: "<<tet4Tmp->getNeigh(2)<<" Neighbours 3: "<<tet4Tmp->getNeigh(3)<<std::endl;
+ for (unsigned int j = 0;j < 4;j++){
+ //std::cout<<" Entering neighbour "<<j<<std::endl;
+ if (tet4Tmp->getNeigh(j) == 0){
+ //std::cout<<" Test Passed "<<std::endl;
+ //find the 4th point,and fill the two vector of the boundary triangle
+ faceXtetFM fxtTmp(tet4Tmp,j);
+ //std::cout<<" fxtTmp created "<<std::endl;
+ for (int k = 0;k < 3;k++){
+ //std::cout<<" Entering Point "<<k<<std::endl;
+ MVertex* toTest = fxtTmp.v[k];
+ if (toTest->onWhat()->dim() == 2){
+ //std::cout<<" Passed First test "<<std::endl;
+ if(BoundaryVertices.find(toTest) == BoundaryVertices.end()){
+ //std::cout<<" Passed Second test "<<std::endl;
+ BoundaryVertices.insert(toTest);
+ }
+ }
+ //std::cout<<" Exiting Point "<<k<<std::endl;
+ }
+ // crackElements.push_back(rTmp->getMeshElement(j));
+ }
+ //std::cout<<" Exiting neighbour "<<j<<std::endl;
+ }
+ //std::cout<<" Exiting tet "<<rTmp->tetrahedra[i]->getNum()<<std::endl;
+ }
+ //std::cout<<" Exiting region "<<rTmp->tag()<<std::endl;
+ }
+ //std::cout<<"Entering computeDistToOppSide"<<std::endl;
+ countdown = 0;
+ for(std::set<MVertex*>::iterator it = BoundaryVertices.begin(); it != BoundaryVertices.end(); it++){
+ countdown++;
+ //std::cout<<" entering Bound Vert "<<(*it)->getNum()<<" it's the "<<countdown<<"st"<<std::endl;
+ MVertex* toCompute = (*it);
+ double resultTmp = computeDistToOppSide(toCompute);
+ int lastTmp = VertexToCorresp[toCompute].size() - 1;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while third time "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ //std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+ //std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while popo "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ //std::cout<<" size of AllDistToOppSide is "<<AllDistToOppSide.size()<<std::endl;
+ //std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+ //std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+// std::cout<<" Testing FaceXTet out of while popo 2 "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ AllDistToOppSide.insert(std::make_pair(toCompute, resultTmp));
+ //std::cout<<" AllDistToOppSide[toCompute] is "<<AllDistToOppSide[toCompute]<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().x()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().y()<<" ; "<<VertexToCorresp[toCompute][lastTmp]->getEndNormal().z()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" AllDistToOppSide[toCompute] adress is "<<&AllDistToOppSide[toCompute]<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1 adress is "<<&(VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1)<<std::endl;
+ //std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+ //std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+ //std::cout<<" ToCompute adress is "<<toCompute<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while popo 3 "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+// AllDistToOppSide[toCompute] = 0.0;
+// std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+// std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+// std::cout<<" size of AllDistToOppSide is "<<AllDistToOppSide.size()<<std::endl;
+// std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+// std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+ AllDistToOppSide[toCompute] = resultTmp;
+ //std::cout<<" observation "<<countdown<<std::endl;
+ //std::cout<<" LastTmp is "<<lastTmp<<std::endl;
+ //std::cout<<" ToCompute is "<<toCompute->getNum()<<std::endl;
+ //std::cout<<" VertexToCorresp[toCompute][lastTmp] is "<<VertexToCorresp[toCompute][lastTmp]<<std::endl;
+// std::cout<<" getEndTriangle is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle()<<std::endl;
+// std::cout<<" t1 is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1<<std::endl;
+// std::cout<<" tet is "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while bis repetita "<<VertexToCorresp[toCompute][lastTmp]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ for (std::map<MVertex*,double>::iterator allDistIt = AllDistToOppSide.begin();allDistIt != AllDistToOppSide.end();allDistIt++){
+ //std::cout<<"allDist of point "<<(*allDistIt).first->getNum()<<" with Pos "<<(*allDistIt).first->x()<<" ; "<<(*allDistIt).first->z()<<" ; "<<(*allDistIt).first->y()<<" is "<<(*allDistIt).second<<std::endl;
+ //std::cout<<" Size of vertexToCorresp "<<VertexToCorresp[(*allDistIt).first].size()<<std::endl;
+ MVertex* toComputeTest = (*allDistIt).first;
+ //std::cout<<" ToComputeTest is "<<toComputeTest->getNum()<<std::endl;
+ int lastTmpTest = VertexToCorresp[toComputeTest].size() - 1;
+// std::cout<<" Testing FaceXTet out of while ter repetita "<<VertexToCorresp[toComputeTest][lastTmpTest]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while fourth time "<<VertexToCorresp[(*allDistIt).first][VertexToCorresp[(*allDistIt).first].size() - 1]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ }
+ }
+
+ //std::cout<<"first observation"<<std::endl;
+ for (std::map<MVertex*,double>::iterator allDistIt = AllDistToOppSide.begin();allDistIt != AllDistToOppSide.end();allDistIt++){
+ //std::cout<<"allDist of point "<<(*allDistIt).first->getNum()<<" with Pos "<<(*allDistIt).first->x()<<" ; "<<(*allDistIt).first->z()<<" ; "<<(*allDistIt).first->y()<<" is "<<(*allDistIt).second<<std::endl;
+ //std::cout<<" Size of vertexToCorresp "<<VertexToCorresp[(*allDistIt).first].size()<<std::endl;
+// std::cout<<" Testing FaceXTet out of while fourth time "<<VertexToCorresp[(*allDistIt).first][VertexToCorresp[(*allDistIt).first].size() - 1]->getEndTriangle().t1->tet()->getNum()<<std::endl;
+ }
+
+ return AllDistToOppSide;
+}
+
+double GMSH_ThinLayerFixMeshPlugin::computeDistToOppSide(MVertex* v){
+ double DistToOppSide = 0.0;
+ //We assume v is on the boundary
+ //First we need to get the internal normal
+ //std::cout<<" Entering ComputeDistToOppSide Function"<<std::endl;
+ SVector3 InteriorNormal = GMSH_ThinLayerFixMeshPlugin::computeInteriorNormal(v);
+ //std::cout<<" Normal computed"<<std::endl;
+ //Then we find the first triangle
+ MTet4* FirstTet = GMSH_ThinLayerFixMeshPlugin::getTetFromPoint(v,InteriorNormal);
+ //std::cout<<" FirstTet is tet "<<FirstTet->tet()->getNum()<<std::endl;
+ //std::cout<<" got FirstTet"<<std::endl;
+ MTet4* CurrentTet = FirstTet;
+ //std::cout<<" set currentTet"<<std::endl;
+ MTet4* PastTet = FirstTet;
+ //std::cout<<" set PastTet"<<std::endl;
+ SPoint3 CurrentPos = SPoint3(v->x(),v->y(),v->z());
+ //std::cout<<" set CurrentPos"<<std::endl;
+ SPoint3 LastPos = CurrentPos;
+ //std::cout<<" set LasPos"<<std::endl;
+ int CurrentTriValue = 0;
+ int* CurrentTri = &CurrentTriValue;
+ //std::cout<<" set CurrentTri"<<std::endl;
+ CorrespVerticesFixMesh* CVTemp = new CorrespVerticesFixMesh();
+ //std::cout<<" Created CVTemp"<<std::endl;
+ CorrespVerticesFixMesh CVTemp2;
+ //std::cout<<" Created CVTemp2"<<std::endl;
+ CVTemp->setStartPoint(v);
+ CVTemp2.setStartPoint(v);
+ //std::cout<<" setStartPoint"<<std::endl;
+ CVTemp->setStartNormal(InteriorNormal);
+ CVTemp2.setStartNormal(InteriorNormal);
+ //std::cout<<" setStartNormal"<<std::endl;
+ FindNewPoint((&CurrentPos),CurrentTri,&CurrentTet,InteriorNormal);
+ DistToOppSide += CurrentPos.distance(LastPos);
+ LastPos = CurrentPos;
+ //std::cout<<" Found New Point"<<std::endl;
+// faceXtetFM fxtCV(PastTet,(*CurrentTri));
+// CVTemp->setEndTriangle(&fxtCV);
+// while(CurrentTet->getNeigh((*CurrentTri)) != 0){
+ int countWhile = 0;
+ int prevprevtet = 1;
+ int prevtet = 2;
+ while((CurrentTet != 0) && (countWhile < 50)){
+ //std::cout<<" Entering While"<<std::endl;
+ //std::cout<<"CurrentTet is "<<CurrentTet->tet()->getNum()<< "adress "<<CurrentTet<< " and currentTri "<<(*CurrentTri)<<std::endl;
+ //std::cout<<"currentPos is "<<CurrentPos.x()<<" ; "<<CurrentPos.y()<<" ; "<<CurrentPos.z()<<std::endl;
+ countWhile++;
+ faceXtetFM fxtCVtmp(PastTet,(*CurrentTri));
+ PastTet = CurrentTet;
+ CVTemp->setEndTrianglePoint1(PastTet->tet()->getVertex(faces[(*CurrentTri)][0]));
+ CVTemp->setEndTrianglePoint2(PastTet->tet()->getVertex(faces[(*CurrentTri)][1]));
+ CVTemp->setEndTrianglePoint3(PastTet->tet()->getVertex(faces[(*CurrentTri)][2]));
+ CVTemp2.setEndTrianglePoint1(PastTet->tet()->getVertex(faces[(*CurrentTri)][0]));
+ CVTemp2.setEndTrianglePoint2(PastTet->tet()->getVertex(faces[(*CurrentTri)][1]));
+ CVTemp2.setEndTrianglePoint3(PastTet->tet()->getVertex(faces[(*CurrentTri)][2]));
+// CVTemp->setEndTriangle(fxtCVtmp);
+// CVTemp2.setEndTriangle(fxtCVtmp);
+ //std::cout<<" FaceXTetCreated"<<std::endl;
+ //std::cout<<" Testing FaceXTet in while "<<fxtCVtmp.t1->tet()->getNum()<<std::endl;
+ //std::cout<<" CurrentTet before FindNewPoint "<<CurrentTet<<std::endl;
+ FindNewPoint((&CurrentPos),CurrentTri,&CurrentTet,InteriorNormal);
+ //std::cout<<" CurrentTet after FindNewPoint "<<CurrentTet<<std::endl;
+ //std::cout<<" FoundNewPoint While"<<std::endl;
+// CurrentTet = CurrentTet->getNeigh((*CurrentTri));
+// //std::cout<<" GotNeigh While"<<std::endl;
+ DistToOppSide += CurrentPos.distance(LastPos);
+ LastPos = CurrentPos;
+
+ if (CurrentTet != 0){
+ if ((CurrentTet->tet()->getNum()) == prevprevtet){
+ //std::cout<<"reached standstill"<<std::endl;
+ while (1){
+
+ }
+ }
+ prevprevtet = prevtet;
+ prevtet=CurrentTet->tet()->getNum();
+ }
+ //std::cout<<" Exiting While"<<std::endl;
+ }
+ //std::cout<<" Out Of While"<<std::endl;
+ CVTemp->setEndPoint(LastPos);
+ CVTemp2.setEndPoint(LastPos);
+ faceXtetFM fxtCV(PastTet,(*CurrentTri));
+ //std::cout<<" Testing FaceXTet out of while "<<fxtCV.t1->tet()->getNum()<<std::endl;
+ CVTemp->setEndTrianglePoint1(PastTet->tet()->getVertex(faces[(*CurrentTri)][0]));
+ CVTemp->setEndTrianglePoint2(PastTet->tet()->getVertex(faces[(*CurrentTri)][1]));
+ CVTemp->setEndTrianglePoint3(PastTet->tet()->getVertex(faces[(*CurrentTri)][2]));
+ CVTemp2.setEndTrianglePoint1(PastTet->tet()->getVertex(faces[(*CurrentTri)][0]));
+ CVTemp2.setEndTrianglePoint2(PastTet->tet()->getVertex(faces[(*CurrentTri)][1]));
+ CVTemp2.setEndTrianglePoint3(PastTet->tet()->getVertex(faces[(*CurrentTri)][2]));
+// CVTemp->setEndTriangle(fxtCV);
+// CVTemp2.setEndTriangle(fxtCV);
+// //std::cout<<" Testing FaceXTet out of while second time "<<CVTemp->getEndTriangle().t1->tet()->getNum()<<std::endl;
+// //std::cout<<" Testing FaceXTet out of while second time CVTemp2 "<<CVTemp2.getEndTriangle().t1->tet()->getNum()<<std::endl;
+// CVTemp->setEndTriangle(&fxtCV);
+// SVector3 EndDir1(fxtCV.v[1]->x() - fxtCV.v[0]->x(),fxtCV.v[1]->y() - fxtCV.v[0]->y(),fxtCV.v[1]->z() - fxtCV.v[0]->z());
+// SVector3 EndDir2(fxtCV.v[2]->x() - fxtCV.v[0]->x(),fxtCV.v[2]->y() - fxtCV.v[0]->y(),fxtCV.v[2]->z() - fxtCV.v[0]->z());
+// SVector3 EndNormal(EndDir1.y() * EndDir2.z() - EndDir1.z() * EndDir2.y(),EndDir1.z() * EndDir2.x() - EndDir1.x() * EndDir2.z(),EndDir1.x() * EndDir2.y() - EndDir1.y() * EndDir2.x());
+ SVector3 EndNormal(-InteriorNormal.x(),-InteriorNormal.y(),-InteriorNormal.z());
+ EndNormal.normalize();
+ CVTemp->setEndNormal(EndNormal);
+ CVTemp2.setEndNormal(EndNormal);
+ CVTemp->setangleProd(fabs(CVTemp->getStartNormal().x() * CVTemp->getEndNormal().x() + CVTemp->getStartNormal().y() * CVTemp->getEndNormal().y() + CVTemp->getStartNormal().z() * CVTemp->getEndNormal().z()));
+ CVTemp2.setangleProd(fabs(CVTemp2.getStartNormal().x() * CVTemp2.getEndNormal().x() + CVTemp2.getStartNormal().y() * CVTemp2.getEndNormal().y() + CVTemp2.getStartNormal().z() * CVTemp2.getEndNormal().z()));
+ CVTemp->setdistP2P(DistToOppSide);
+ CVTemp2.setdistP2P(DistToOppSide);
+ if ((CVTemp->getangleProd() > angleMax) &&(CVTemp->getdistP2P() < distP2PMax)){
+ CVTemp->setActive(true);
+ }
+ else{
+ CVTemp->setActive(false);
+ }
+ if ((CVTemp2.getangleProd() > angleMax) &&(CVTemp2.getdistP2P() < distP2PMax)){
+ CVTemp2.setActive(true);
+ }
+ else{
+ CVTemp2.setActive(false);
+ }
+ CVTemp->setTagMaster(-2);
+ CVTemp2.setTagMaster(-2);
+ //std::cout<<" getEndNormal is "<<CVTemp2.getEndNormal().x()<<" ; "<<CVTemp2.getEndNormal().y()<<" ; "<<CVTemp2.getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<CVTemp2.getEndNormal().x()<<" ; "<<CVTemp2.getEndNormal().y()<<" ; "<<CVTemp2.getEndNormal().z()<<std::endl;
+
+ VertexToCorresp[v].push_back(CVTemp);
+// VertexToCorresp[v].push_back(&CVTemp2);
+ //std::cout<<" getEndNormal test is "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().x()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().y()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().x()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().y()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().x()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().y()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().z()<<std::endl;
+ //std::cout<<" getEndNormal is "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().x()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().y()<<" ; "<<VertexToCorresp[v][VertexToCorresp[v].size() - 1]->getEndNormal().z()<<std::endl;
+
+ return DistToOppSide;
+}
+
+SVector3 GMSH_ThinLayerFixMeshPlugin::computeInteriorNormal(MVertex* v){
+ SPoint3 InteriorNormal(0.0,0.0,0.0);
+ //std::cout<<" Entering computeInteriorNormal"<<std::endl;
+ std::vector<MTetrahedron*> currentVecTet = VertexToTets[v];
+ std::vector<SPoint3> vecInteriorNodes;
+ std::vector<SPoint3> vecFirstDir;
+ std::vector<SPoint3> vecSecondDir;
+ for (unsigned int i = 0;i < currentVecTet.size();i++){
+ //std::cout<<" Entering for "<<i<<std::endl;
+ MTet4* tet4Tmp = TetToTet4[currentVecTet[i]];
+ for (int j = 0;j < 4 ; j++){
+ //std::cout<<" Entering for "<<j<<std::endl;
+ if (tet4Tmp->getNeigh(j) == 0){
+ //std::cout<<" Test Passed neigh zero "<<std::endl;
+ //find the 4th point,and fill the two vector of the boundary triangle
+ faceXtetFM fxtTmp(tet4Tmp,j);
+ //std::cout<<" Face Created "<<j<<std::endl;
+ for (int k = 0;k < 4;k++){
+ //std::cout<<" Entering for "<<k<<std::endl;
+ bool foundInteriorPoint = true;
+ for (int l = 0;l < 3;l++){
+ if (tet4Tmp->tet()->getVertex(k) == fxtTmp.v[l]){
+ foundInteriorPoint = false;
+ }
+ }
+ if (foundInteriorPoint){
+ SPoint3 pointTmp(tet4Tmp->tet()->getVertex(k)->x(),tet4Tmp->tet()->getVertex(k)->y(),tet4Tmp->tet()->getVertex(k)->z());
+ vecInteriorNodes.push_back(pointTmp);
+ }
+ }
+ SPoint3 pointTmp1(fxtTmp.v[1]->x() - fxtTmp.v[0]->x(),fxtTmp.v[1]->y() - fxtTmp.v[0]->y(),fxtTmp.v[1]->z() - fxtTmp.v[0]->z());
+ SPoint3 pointTmp2(fxtTmp.v[2]->x() - fxtTmp.v[0]->x(),fxtTmp.v[2]->y() - fxtTmp.v[0]->y(),fxtTmp.v[2]->z() - fxtTmp.v[0]->z());
+ vecFirstDir.push_back(pointTmp1);
+ vecSecondDir.push_back(pointTmp2);
+ }
+ }
+ }
+ //std::cout<<" Out of loop for "<<std::endl;
+ //at this point we have all the info necessary.
+ SPoint3 pointInteriorAverage(0.0,0.0,0.0);
+ for (unsigned int i = 0;i < vecInteriorNodes.size();i++){
+ pointInteriorAverage += SPoint3(vecInteriorNodes[i].x(),vecInteriorNodes[i].y(),vecInteriorNodes[i].z());
+// pointInteriorAverage.x() += vecInteriorNodes[i].x();
+// pointInteriorAverage.y() += vecInteriorNodes[i].y();
+// pointInteriorAverage.z() += vecInteriorNodes[i].z();
+ }
+ pointInteriorAverage = SPoint3(pointInteriorAverage.x() / (double(vecInteriorNodes.size())),pointInteriorAverage.y() / (double(vecInteriorNodes.size())),pointInteriorAverage.z() / (double(vecInteriorNodes.size())));
+// pointInteriorAverage.x() = pointInteriorAverage.x() / (double(vecInteriorNodes.size()));
+// pointInteriorAverage.y() = pointInteriorAverage.y() / (double(vecInteriorNodes.size()));
+// pointInteriorAverage.z() = pointInteriorAverage.z() / (double(vecInteriorNodes.size()));
+ SPoint3 dirInteriorAverage(pointInteriorAverage.x() - v->x(),pointInteriorAverage.y() - v->y(),pointInteriorAverage.z() - v->z());
+ double norme = sqrt(dirInteriorAverage.x() * dirInteriorAverage.x() + dirInteriorAverage.y() * dirInteriorAverage.y() + dirInteriorAverage.z() * dirInteriorAverage.z());
+ dirInteriorAverage = SPoint3(dirInteriorAverage.x() / norme,dirInteriorAverage.y() / norme,dirInteriorAverage.z() / norme);
+// dirInteriorAverage.x() = dirInteriorAverage.x() / norme;
+// dirInteriorAverage.y() = dirInteriorAverage.y() / norme;
+// dirInteriorAverage.z() = dirInteriorAverage.z() / norme;
+ std::vector<SPoint3> vecOrthogDir;
+ for(unsigned int i = 0;i < vecFirstDir.size();i++){
+ SPoint3 pointTmp(vecFirstDir[i].y() * vecSecondDir[i].z() - vecFirstDir[i].z() * vecSecondDir[i].y(),vecFirstDir[i].z() * vecSecondDir[i].x() - vecFirstDir[i].x() * vecSecondDir[i].z(),vecFirstDir[i].x() * vecSecondDir[i].y() - vecFirstDir[i].y() * vecSecondDir[i].x());
+ vecOrthogDir.push_back(pointTmp);
+ }
+ for(unsigned int i = 0;i < vecOrthogDir.size();i++){
+ double prodScalTmp = vecOrthogDir[i].x() * dirInteriorAverage.x() + vecOrthogDir[i].y() * dirInteriorAverage.y() + vecOrthogDir[i].z() * dirInteriorAverage.z();
+ if (prodScalTmp < 0.0){
+ vecOrthogDir[i] = SPoint3(- vecOrthogDir[i].x(),- vecOrthogDir[i].y(),- vecOrthogDir[i].z());
+// vecOrthogDir[i].x() = - vecOrthogDir[i].x();
+// vecOrthogDir[i].y() = - vecOrthogDir[i].y();
+// vecOrthogDir[i].z() = - vecOrthogDir[i].z();
+ }
+ double normeTmp = sqrt(vecOrthogDir[i].x() * vecOrthogDir[i].x() + vecOrthogDir[i].y() * vecOrthogDir[i].y() + vecOrthogDir[i].z() * vecOrthogDir[i].z());
+ vecOrthogDir[i] = SPoint3(vecOrthogDir[i].x() / normeTmp,vecOrthogDir[i].y() / normeTmp,vecOrthogDir[i].z() / normeTmp);
+// vecOrthogDir[i].x() = vecOrthogDir[i].x() / normeTmp;
+// vecOrthogDir[i].y() = vecOrthogDir[i].y() / normeTmp;
+// vecOrthogDir[i].z() = vecOrthogDir[i].z() / normeTmp;
+ InteriorNormal += SPoint3(vecOrthogDir[i].x(),vecOrthogDir[i].y(),vecOrthogDir[i].z());
+// InteriorNormal.x() += vecOrthogDir[i].x();
+// InteriorNormal.y() += vecOrthogDir[i].y();
+// InteriorNormal.z() += vecOrthogDir[i].z();
+ }
+ norme = sqrt(InteriorNormal.x() * InteriorNormal.x() + InteriorNormal.y() * InteriorNormal.y() + InteriorNormal.z() * InteriorNormal.z());
+ InteriorNormal = SPoint3(InteriorNormal.x() / norme,InteriorNormal.y() / norme,InteriorNormal.z() / norme);
+// InteriorNormal.x() = InteriorNormal.x() / norme;
+// InteriorNormal.y() = InteriorNormal.y() / norme;
+// InteriorNormal.z() = InteriorNormal.z() / norme;
+ return InteriorNormal;
+}
+
+MTet4* GMSH_ThinLayerFixMeshPlugin::getTetFromPoint(MVertex* v, SVector3 InteriorNormal){
+ MTet4* TetToGet = 0;
+ std::vector<MTetrahedron*> currentVecTet = VertexToTets[v];
+ //std::cout<<" entering getTetFromPoint, vertex "<<v->x()<<" ; "<<v->y()<<" ; "<<v->z()<<" and dir "<<InteriorNormal.x()<<" ; "<<InteriorNormal.y()<<" ; "<<InteriorNormal.z()<<std::endl;
+ for (unsigned int i = 0;i < currentVecTet.size();i++){
+ std::vector<SVector3> vecDir;
+ for (int j = 0;j < 4 ; j++){
+ if (currentVecTet[i]->getVertex(j) != v){
+ SVector3 DirTmp(currentVecTet[i]->getVertex(j)->x() - v->x(),currentVecTet[i]->getVertex(j)->y() - v->y(),currentVecTet[i]->getVertex(j)->z() - v->z());
+ vecDir.push_back(DirTmp);
+ }
+ }
+ bool IsPositiv = GMSH_ThinLayerFixMeshPlugin::IsPositivOrientation(vecDir[0],vecDir[1],vecDir[2]);
+ if (!IsPositiv){
+ SVector3 DirTmp1 = vecDir[1];
+ SVector3 DirTmp2 = vecDir[0];
+ SVector3 DirTmp3 = vecDir[2];
+ vecDir.clear();
+ vecDir.push_back(DirTmp1);
+ vecDir.push_back(DirTmp2);
+ vecDir.push_back(DirTmp3);
+ }
+ bool isPositiv1 = GMSH_ThinLayerFixMeshPlugin::IsPositivOrientation(vecDir[0],vecDir[1],InteriorNormal);
+ bool isPositiv2 = GMSH_ThinLayerFixMeshPlugin::IsPositivOrientation(vecDir[1],vecDir[2],InteriorNormal);
+ bool isPositiv3 = GMSH_ThinLayerFixMeshPlugin::IsPositivOrientation(vecDir[2],vecDir[0],InteriorNormal);
+ if (isPositiv1){
+ if (isPositiv2){
+ if (isPositiv3){
+ TetToGet = TetToTet4[currentVecTet[i]];
+ //std::cout<<" Found one fitting ! "<<TetToGet->tet()->getNum()<<std::endl;
+ }
+ }
+ }
+ }
+ //std::cout<<" exiting getTetFromPoint with result "<<TetToGet<<std::endl;
+ //std::cout<<" getTetFromPoint is tet "<<TetToGet->tet()->getNum()<<std::endl;
+ return TetToGet;
+}
+
+bool GMSH_ThinLayerFixMeshPlugin::IsPositivOrientation(SVector3 a, SVector3 b, SVector3 c){
+ bool result = false;
+ SPoint3 ProdVec(a.y() * b.z() - a.z() * b.y(),a.z() * b.x() - a.x() * b.z(),a.x() * b.y() - a.y() * b.x());
+ double ProdScal = ProdVec.x() * c.x() + ProdVec.y() * c.y() + ProdVec.z() * c.z();
+ if (ProdScal >= 0.0){
+ result = true;
+ }
+ return result;
+}
+
+
+void GMSH_ThinLayerFixMeshPlugin::FindNewPoint(SPoint3* CurrentPoint, int* CurrentTri, MTet4** CurrentTet, SVector3 InteriorNormal){
+ double distanceP2P = 0.0;
+ double alphaMax;
+ double betaMax;
+ SPoint3 ResultPoint;
+ int triToGet = -1;
+ //std::cout<<" Entered FindNewPoint"<<std::endl;
+ for (unsigned int n = 0;n < 4 ; n++){
+ //calculer matrice a inverser
+ faceXtetFM fxt((*CurrentTet),n);
+ double a = fxt.v[1]->x() - fxt.v[0]->x();
+ double b = fxt.v[2]->x() - fxt.v[0]->x();
+ double c = InteriorNormal.x();
+ double d = fxt.v[1]->y() - fxt.v[0]->y();
+ double e = fxt.v[2]->y() - fxt.v[0]->y();
+ double f = InteriorNormal.y();
+ double g = fxt.v[1]->z() - fxt.v[0]->z();
+ double h = fxt.v[2]->z() - fxt.v[0]->z();
+ double i = InteriorNormal.z();
+ //produit matrice inverse par vecteur donne poids
+ double detMat = a * e * i + b * f * g + c * d * h - c * e * g - f * h * a - i * b * d;
+ double ai = e * i - f * h;
+ double bi = c * h - b * i;
+ double ci = b * f - c * e;
+ double di = f * g - d * i;
+ double ei = a * i - c * g;
+ double fi = c * d - a * f;
+// double gi = d * h - e * g;
+// double hi = b * g - a * h;
+// double ii = a * e - b * d;
+ double oppx = (*CurrentPoint).x() - fxt.v[0]->x();
+ double oppy = (*CurrentPoint).y() - fxt.v[0]->y();
+ double oppz = (*CurrentPoint).z() - fxt.v[0]->z();
+ double alpha = ai / detMat * oppx + bi / detMat * oppy + ci / detMat * oppz;
+ double beta = di / detMat * oppx + ei / detMat * oppy + fi / detMat * oppz;
+// double gamma = gi / detMat * oppx + hi / detMat * oppy + ii / detMat * oppz;
+ //Test si poids entre 0 et 1 et si length maximale
+ if ((alpha >= (0.0-epsilon)) && (alpha <= (1.0 + epsilon))){
+ if ((beta >= (0.0-epsilon)) && (beta <= (1.0 + epsilon))){
+ if (((1.0 - alpha - beta) >= (0.0-epsilon)) && ((1.0 - alpha - beta) <= (1.0 + epsilon))){
+ SPoint3 PointTmp(fxt.v[0]->x() + alpha * (fxt.v[1]->x() - fxt.v[0]->x()) + beta * (fxt.v[2]->x() - fxt.v[0]->x()),fxt.v[0]->y() + alpha * (fxt.v[1]->y() - fxt.v[0]->y()) + beta * (fxt.v[2]->y() - fxt.v[0]->y()),fxt.v[0]->z() + alpha * (fxt.v[1]->z() - fxt.v[0]->z()) + beta * (fxt.v[2]->z() - fxt.v[0]->z()));
+ double distanceTmp = PointTmp.distance((*CurrentPoint));
+ if (distanceTmp > distanceP2P){
+ distanceP2P = distanceTmp;
+ ResultPoint = PointTmp;
+ triToGet = n;
+ alphaMax = alpha;
+ betaMax = beta;
+ }
+ }
+ }
+ }
+ }
+ //std::cout<<" End of For loop"<<std::endl;
+ //test si trop proche d'un point / une arete
+ if (((alphaMax < epsilon) && (betaMax < epsilon)) || ((alphaMax < epsilon) && ((1.0 - alphaMax - betaMax) < epsilon)) || (((1.0 - alphaMax - betaMax) < epsilon) && (betaMax < epsilon))){
+ //proche d'un point
+ //std::cout<<" Close to point"<<std::endl;
+ double DistMinTmp = 10000000.0;
+ int indexMinTmp;
+ for (unsigned int i = 0;i < 4;i++){
+ double distanceTmp = sqrt(((*CurrentTet)->tet()->getVertex(i)->x() - ResultPoint.x()) * ((*CurrentTet)->tet()->getVertex(i)->x() - ResultPoint.x()) + ((*CurrentTet)->tet()->getVertex(i)->y() - ResultPoint.y()) * ((*CurrentTet)->tet()->getVertex(i)->y() - ResultPoint.y()) + ((*CurrentTet)->tet()->getVertex(i)->z() - ResultPoint.z()) * ((*CurrentTet)->tet()->getVertex(i)->z() - ResultPoint.z()));
+ if (distanceTmp < DistMinTmp){
+ DistMinTmp = distanceTmp;
+ indexMinTmp = i;
+ }
+ }
+ ////std::cout<<"NewTet before is "<<NewTet<<std::endl;
+ MTet4* NewTet = GMSH_ThinLayerFixMeshPlugin::getTetFromPoint((*CurrentTet)->tet()->getVertex(indexMinTmp),InteriorNormal);
+ //std::cout<<"NewTet after is "<<NewTet<<std::endl;
+ SPoint3 PointTmp((*CurrentTet)->tet()->getVertex(indexMinTmp)->x(),(*CurrentTet)->tet()->getVertex(indexMinTmp)->y(),(*CurrentTet)->tet()->getVertex(indexMinTmp)->z());
+ (*CurrentPoint) = PointTmp;
+ (*CurrentTet) = NewTet;
+ }
+ else if ((alphaMax < epsilon) || (betaMax < epsilon) || ((1.0 - alphaMax - betaMax) < epsilon)){
+ //trop proche d'une arete
+ //std::cout<<" Close to edge"<<std::endl;
+ }
+ else{
+ //std::cout<<" Close to nothing"<<std::endl;
+ //std::cout<<" ResultPoint is "<<ResultPoint.x()<<"; "<<ResultPoint.y()<<"; "<<ResultPoint.z()<<"; "<<" and tritoget is "<<triToGet<<std::endl;
+ //std::cout<<" CurrentPoint is "<<(*CurrentPoint).x()<<"; "<<(*CurrentPoint).y()<<"; "<<(*CurrentPoint).z()<<std::endl;
+ (*CurrentPoint) = ResultPoint;
+ //std::cout<<" test 1"<<std::endl;
+ //std::cout<<" CurrentTri is "<<(*CurrentTri)<<std::endl;
+ (*CurrentTri) = triToGet;
+ //std::cout<<" test 2"<<std::endl;
+ //std::cout<<" CurrentTet is "<<(*CurrentTet)<<" and has neighbours "<< (*CurrentTet)->getNeigh(0)<<" ; "<< (*CurrentTet)->getNeigh(1)<<" ; "<< (*CurrentTet)->getNeigh(2)<<" ; "<< (*CurrentTet)->getNeigh(3)<<std::endl;
+ (*CurrentTet) = (*CurrentTet)->getNeigh(triToGet);
+ //std::cout<<" CurrentTet has been changed to "<<(*CurrentTet)<<std::endl;
+ //std::cout<<" test 3"<<std::endl;
+ }
+ //std::cout<<" Exit FindNewPoint"<<std::endl;
+}
+
+void GMSH_ThinLayerFixMeshPlugin::fillVertexToTets(){
+ GModel *m = GModel::current();
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+ MTetrahedron* elem = rTmp->tetrahedra[i];
+ for (unsigned int j = 0; j < 4;j++){
+ std::vector<MTetrahedron*> emptyTetVec;
+ emptyTetVec.clear();
+ VertexToTets[elem->getVertex(j)] = emptyTetVec;
+ std::vector<CorrespVerticesFixMesh*> emptyCVVec;
+ emptyCVVec.clear();
+ VertexToCorresp[elem->getVertex(j)] = emptyCVVec;
+ }
+ }
+ }
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* rTmp = (*itr);
+ for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+ MTetrahedron* elem = rTmp->tetrahedra[i];
+ for (unsigned int j = 0; j < 4;j++){
+ VertexToTets[elem->getVertex(j)].push_back(elem);
+ }
+ }
+ }
+}
+
+static void setLcsFM(MTriangle *t, std::map<MVertex*, double> &vSizes,
+ std::set<MVertex*> &bndVertices)
+{
+ for(int i = 0; i < 3; i++){
+ bndVertices.insert(t->getVertex(i));
+ MEdge e = t->getEdge(i);
+ MVertex *vi = e.getVertex(0);
+ MVertex *vj = e.getVertex(1);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ // use largest edge length
+ if (iti == vSizes.end() || iti->second < l) vSizes[vi] = l;
+ if (itj == vSizes.end() || itj->second < l) vSizes[vj] = l;
+ }
+}
+
+static void setLcsFM(MTetrahedron *t, std::map<MVertex*, double> &vSizes,
+ std::set<MVertex*> &bndVertices)
+{
+ for (int i = 0; i < 4; i++){
+ for (int j = i + 1; j < 4; j++){
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ std::set<MVertex*>::iterator itvi = bndVertices.find(vi);
+ std::set<MVertex*>::iterator itvj = bndVertices.find(vj);
+ // smallest tet edge
+ if (itvi == bndVertices.end() &&
+ (iti == vSizes.end() || iti->second > l)) vSizes[vi] = l;
+ if (itvj == bndVertices.end() &&
+ (itj == vSizes.end() || itj->second > l)) vSizes[vj] = l;
+ }
+ }
+}
+
+void GMSH_ThinLayerFixMeshPlugin::fillTetToTet4(){
+ GModel *m = GModel::current();
+ std::vector<double> vSizes;
+ std::vector<double> vSizesBGM;
+ MTet4Factory myFactory(1600000);
+ std::set<MTet4*, compareTet4Ptr> &allTets = myFactory.getAllTets();
+ std::set<MTet4*, compareTet4Ptr> activeTets;
+ int NUM = 0;
+
+ for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+ GRegion* gr = (*itr);
+ { // leave this in a block so the map gets deallocated directly
+ std::map<MVertex*, double> vSizesMap;
+ std::set<MVertex*> bndVertices;
+ for(GModel::fiter it = gr->model()->firstFace(); it != gr->model()->lastFace(); ++it){
+ GFace *gf = *it;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ setLcsFM(gf->triangles[i], vSizesMap, bndVertices);
+ }
+ }
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
+ setLcsFM(gr->tetrahedra[i], vSizesMap, bndVertices);
+ for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
+ it != vSizesMap.end(); ++it){
+ it->first->setIndex(NUM++);
+ vSizes.push_back(it->second);
+ vSizesBGM.push_back(it->second);
+ }
+ }
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
+ {
+ MTet4* currentTet4 = myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM);
+ TetToTet4[gr->tetrahedra[i]] = currentTet4;
+ allTets.insert(currentTet4);
+ }
+ }
+
+
+ std::vector<MTet4*> vecAllTet4;
+ vecAllTet4.clear();
+ for (std::set<MTet4*, compareTet4Ptr>::iterator itTp = allTets.begin();itTp != allTets.end();itTp++){
+ vecAllTet4.push_back((*itTp));
+// //std::cout<<"inserted "<<(*itTp)->tet()->getNum()<<std::endl;
+ }
+ connectTets(vecAllTet4);
+
+
+
+// GModel *m = GModel::current();
+// std::vector<MTet4*> vecAllTet4;
+// vecAllTet4.clear();
+// for (GModel::riter itr= m->firstRegion();itr != m->lastRegion();itr++){
+// GRegion* rTmp = (*itr);
+// for (unsigned int i = 0; i < rTmp->tetrahedra.size();i++){
+// MTetrahedron* elem = rTmp->tetrahedra[i];
+// MTet4 tet4tmp = MTet4(elem,0.0);
+// MTet4* currentTet4 = &tet4tmp;
+// TetToTet4[elem] = currentTet4;
+// vecAllTet4.push_back(currentTet4);
+// }
+// }
+// connectTets(vecAllTet4);
+}
+
+/****************static declarations****************/
+
+std::map<MVertex*,std::vector<MTetrahedron*> > GMSH_ThinLayerFixMeshPlugin::VertexToTets;
+std::map<MTetrahedron*,MTet4*> GMSH_ThinLayerFixMeshPlugin::TetToTet4;
+std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> > GMSH_ThinLayerFixMeshPlugin::VertexToCorresp;
+std::vector<std::vector<CorrespVerticesFixMesh*> > GMSH_ThinLayerFixMeshPlugin::vecOfThinSheets;
+
diff --git a/Plugin/ThinLayerFixMesh.h b/Plugin/ThinLayerFixMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad412189247c0534e2d5bc7f62e3edcf17a35d74
--- /dev/null
+++ b/Plugin/ThinLayerFixMesh.h
@@ -0,0 +1,178 @@
+/*
+ * ThinLayerFixMesh.h
+ *
+ * Created on: Oct 13, 2014
+ * Author: nicolas
+ */
+
+#ifndef THINLAYERFIXMESH_H_
+#define THINLAYERFIXMESH_H_
+
+#include "Plugin.h"
+#include "MVertex.h"
+#include "MTriangle.h"
+#include "meshGRegionDelaunayInsertion.h"
+
+extern "C"
+{
+ GMSH_Plugin *GMSH_RegisterThinLayerFixMeshPlugin();
+}
+
+static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}};
+
+struct faceXtetFM{
+ MVertex *v[3],*unsorted[3];
+ MTet4 *t1;
+ int i1;
+ faceXtetFM(MTet4 *_t=0, int iFac=0) : t1(_t), i1(iFac)
+ {
+// std::cout<<"entering faceXTet"<<std::endl;
+// std::cout<<"tag of tet is "<<t1->tet()->getNum()<<std::endl;
+// std::cout<<"first vec of tet is "<<t1->tet()->getVertex(0)->getNum()<<std::endl;
+// std::cout<<"second vec of tet is "<<t1->tet()->getVertex(1)->getNum()<<std::endl;
+// std::cout<<"third vec of tet is "<<t1->tet()->getVertex(2)->getNum()<<std::endl;
+// std::cout<<"fourth vec of tet is "<<t1->tet()->getVertex(3)->getNum()<<std::endl;
+ MVertex *v0 = t1->tet()->getVertex(faces[iFac][0]);
+// std::cout<<"after first vert"<<std::endl;
+ MVertex *v1 = t1->tet()->getVertex(faces[iFac][1]);
+// std::cout<<"after second vert"<<std::endl;
+ MVertex *v2 = t1->tet()->getVertex(faces[iFac][2]);
+// std::cout<<"after third vert"<<std::endl;
+
+ unsorted[0] = v0;
+// std::cout<<"after first unsorted"<<std::endl;
+ unsorted[1] = v1;
+// std::cout<<"after second unsorted"<<std::endl;
+ unsorted[2] = v2;
+// std::cout<<"after third unsorted"<<std::endl;
+
+ v[0] = std::min(std::min(v0,v1),v2);
+// std::cout<<"after min"<<std::endl;
+ v[2] = std::max(std::max(v0,v1),v2);
+// std::cout<<"after max"<<std::endl;
+ v[1] = (v0 != v[0] && v0 != v[2]) ? v0 : (v1 != v[0] && v1 != v[2]) ? v1 : v2;
+// std::cout<<"after minimax"<<std::endl;
+ //
+ // std::sort(v, v + 3);
+ }
+
+ inline MVertex * getVertex (int i) const { return unsorted[i];}
+
+ inline bool operator < (const faceXtetFM & other) const
+ {
+ if (v[0] < other.v[0]) return true;
+ if (v[0] > other.v[0]) return false;
+ if (v[1] < other.v[1]) return true;
+ if (v[1] > other.v[1]) return false;
+ if (v[2] < other.v[2]) return true;
+ return false;
+ }
+ inline bool operator == (const faceXtetFM & other) const
+ {
+ return (v[0] == other.v[0] &&
+ v[1] == other.v[1] &&
+ v[2] == other.v[2] );
+ }
+ bool visible (MVertex *v){
+ MVertex* v0 = t1->tet()->getVertex(faces[i1][0]);
+ MVertex* v1 = t1->tet()->getVertex(faces[i1][1]);
+ MVertex* v2 = t1->tet()->getVertex(faces[i1][2]);
+ double a[3] = {v0->x(),v0->y(),v0->z()};
+ double b[3] = {v1->x(),v1->y(),v1->z()};
+ double c[3] = {v2->x(),v2->y(),v2->z()};
+ double d[3] = {v->x(),v->y(),v->z()};
+ double o = robustPredicates :: orient3d(a,b,c,d);
+ return o < 0;
+ }
+};
+
+class CorrespVerticesFixMesh{
+private:
+ MVertex* StartPoint;
+ SPoint3 EndPoint;
+ SVector3 StartNormal;
+ SVector3 EndNormal;
+ MVertex* EndTrianglePoint1;
+ MVertex* EndTrianglePoint2;
+ MVertex* EndTrianglePoint3;
+// faceXtetFM EndTriangle;
+ double distP2P;
+ double angleProd;
+ bool Active;
+ bool EndTriangleActive;
+ bool IsMaster;
+ int tagMaster;
+public:
+ CorrespVerticesFixMesh();
+ ~CorrespVerticesFixMesh();
+ void setStartPoint(MVertex* v);
+ void setEndPoint(SPoint3 p);
+ void setStartNormal(SVector3 v);
+ void setEndNormal(SVector3 v);
+// void setEndTriangle(faceXtetFM f);
+ void setEndTrianglePoint1(MVertex* v);
+ void setEndTrianglePoint2(MVertex* v);
+ void setEndTrianglePoint3(MVertex* v);
+ void setdistP2P(double d);
+ void setangleProd(double a);
+ void setActive(bool b);
+ void setEndTriangleActive(bool b);
+ void setIsMaster(bool b);
+ void setTagMaster(int i);
+ MVertex* getStartPoint();
+ SPoint3 getEndPoint();
+ SVector3 getStartNormal();
+ SVector3 getEndNormal();
+// faceXtetFM getEndTriangle();
+ MVertex* getEndTrianglePoint1();
+ MVertex* getEndTrianglePoint2();
+ MVertex* getEndTrianglePoint3();
+ double getdistP2P();
+ double getangleProd();
+ bool getActive();
+ bool getEndTriangleActive();
+ bool getIsMaster();
+ int getTagMaster();
+};
+
+class GMSH_ThinLayerFixMeshPlugin : public GMSH_PostPlugin
+{
+private:
+public:
+ GMSH_ThinLayerFixMeshPlugin(){}
+ //~GMSH_ThinLayerFixMeshPlugin();
+ std::string getName() const { return "ThinLayerFixMesh"; }
+ std::string getShortHelp() const
+ {
+ return "Fix the mesh in thin parts";
+ }
+ std::string getHelp() const;
+ int getNbOptions() const;
+ StringXNumber* getOption(int iopt);
+ PView *execute(PView *);
+ static void perform();
+ static void checkOppositeTriangles();
+ static void fillvecOfThinSheets();
+ static std::map<MVertex*,double> computeAllDistToOppSide();
+ static double computeDistToOppSide(MVertex* v);
+ static SVector3 computeInteriorNormal(MVertex* v);
+ static MTet4* getTetFromPoint(MVertex* v, SVector3 InteriorNormal);
+ static bool IsPositivOrientation(SVector3 a, SVector3 b, SVector3 c);
+ static void FindNewPoint(SPoint3* CurrentPoint, int* CurrentTri, MTet4** CurrentTet, SVector3 InteriorNormal);
+ static std::map<MVertex*,std::vector<MTetrahedron*> > VertexToTets;
+ static std::map<MTetrahedron*,MTet4*> TetToTet4;
+ static std::map<MVertex*,std::vector<CorrespVerticesFixMesh*> > VertexToCorresp;
+ static std::vector<std::vector<CorrespVerticesFixMesh*> > vecOfThinSheets;
+// static const double epsilon = 0.00001;
+ static const double epsilon = 0.00000000001;
+ static const double angleMax = 0.9;
+// static const double distP2PMax = 0.3;
+ static const double distP2PMax = 5.0;
+ static void fillVertexToTets();
+ static void fillTetToTet4();
+};
+
+
+
+
+#endif /* THINLAYERFIXMESH_H_ */