From 7bcee9e14c900556dc79862ef5550285af773c29 Mon Sep 17 00:00:00 2001 From: Nicolas Kowalski <kowalski.nico@gmail.com> Date: Thu, 26 Nov 2015 22:58:29 +0000 Subject: [PATCH] Addition of the Thin Layer plugin. Requires the volume to be meshed with tetra, and remesh the surface, constraining vertices. After, the mesh can be cleaned and rebuilt, and the new constrained vertices will be used instead of the original ones. Constrained vertices appears in thin parts of the domain to ensure each vertex has a corresponding one on the other surface. --- CMakeLists.txt | 2 +- Common/GmshDefines.h | 1 + Geo/GFace.h | 4 + Mesh/CMakeLists.txt | 2 +- Mesh/ThinLayer.cpp | 529 ++++++++++++ Mesh/ThinLayer.h | 137 ++++ Mesh/meshGFace.cpp | 12 +- Mesh/meshGFaceDelaunayInsertion.cpp | 107 +++ Mesh/meshGFaceDelaunayInsertion.h | 4 + Mesh/simple3D.cpp | 1 + Mesh/surfaceFiller.cpp | 168 ++++ Mesh/surfaceFiller.h | 10 +- Plugin/CMakeLists.txt | 2 +- Plugin/DuplicateBoundaries.cpp | 938 +-------------------- Plugin/PluginManager.cpp | 3 + Plugin/ThinLayerFixMesh.cpp | 1162 +++++++++++++++++++++++++++ Plugin/ThinLayerFixMesh.h | 178 ++++ 17 files changed, 2319 insertions(+), 941 deletions(-) create mode 100644 Mesh/ThinLayer.cpp create mode 100644 Mesh/ThinLayer.h create mode 100644 Plugin/ThinLayerFixMesh.cpp create mode 100644 Plugin/ThinLayerFixMesh.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 16393c6634..6cb12f9dd3 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 d68ec33f65..6a3e985ea7 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 93b55f9a2d..97e21ab0c8 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 c1cef52f21..53c2aaca3b 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 0000000000..8b6749a692 --- /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 0000000000..de9e37c819 --- /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 197e197ed7..fdb3717b6d 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 14cad8fad2..106365a53a 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 bd627df025..79a8c389af 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 efda800069..95ffaf3025 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 e205176e72..d352c5bf98 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 dabb87040b..d317cb0180 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 e542921345..1853783794 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 04fdb16943..33e31b7f3e 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 14c1de6efc..c460a099d6 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 0000000000..fc13bc63e5 --- /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 0000000000..ad41218924 --- /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_ */ -- GitLab