1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Author: Nicolas Kowalski
7 
8 #include "ThinLayer.h"
9 #include "GModel.h"
10 #include "robustPredicates.h"
11 #include "GRegion.h"
12 
13 std::map<MVertex *, std::vector<MTetrahedron *> > ThinLayer::VertexToTets;
14 std::map<MTetrahedron *, MTet4 *> ThinLayer::TetToTet4;
15 std::map<MVertex *, std::vector<CorrespVertices *> > ThinLayer::VertexToCorresp;
16 std::vector<std::vector<CorrespVertices *> > ThinLayer::vecOfThinSheets;
17 const double ThinLayer::epsilon = 0.00001;
18 const double ThinLayer::angleMax = 0.9;
19 const double ThinLayer::distP2PMax = 0.3;
20 
CorrespVertices()21 CorrespVertices::CorrespVertices() {}
~CorrespVertices()22 CorrespVertices::~CorrespVertices() {}
setStartPoint(MVertex * v)23 void CorrespVertices::setStartPoint(MVertex *v) { this->StartPoint = v; }
setEndPoint(const SPoint3 & p)24 void CorrespVertices::setEndPoint(const SPoint3 &p) { this->EndPoint = p; }
setStartNormal(const SVector3 & v)25 void CorrespVertices::setStartNormal(const SVector3 &v)
26 {
27   this->StartNormal = v;
28 }
setEndNormal(const SVector3 & v)29 void CorrespVertices::setEndNormal(const SVector3 &v) { this->EndNormal = v; }
setEndTriangle(faceXtet f)30 void CorrespVertices::setEndTriangle(faceXtet f) { this->EndTriangle = f; }
setdistP2P(double d)31 void CorrespVertices::setdistP2P(double d) { this->distP2P = d; }
setangleProd(double a)32 void CorrespVertices::setangleProd(double a) { this->angleProd = a; }
setActive(bool b)33 void CorrespVertices::setActive(bool b) { this->Active = b; }
setEndTriangleActive(bool b)34 void CorrespVertices::setEndTriangleActive(bool b)
35 {
36   this->EndTriangleActive = b;
37 }
setIsMaster(bool b)38 void CorrespVertices::setIsMaster(bool b) { this->IsMaster = b; }
setTagMaster(int i)39 void CorrespVertices::setTagMaster(int i) { this->tagMaster = i; }
getStartPoint()40 MVertex *CorrespVertices::getStartPoint() { return StartPoint; }
getEndPoint()41 SPoint3 CorrespVertices::getEndPoint() { return EndPoint; }
getStartNormal()42 SVector3 CorrespVertices::getStartNormal() { return StartNormal; }
getEndNormal()43 SVector3 CorrespVertices::getEndNormal() { return EndNormal; }
getEndTriangle()44 faceXtet CorrespVertices::getEndTriangle() { return EndTriangle; }
getdistP2P()45 double CorrespVertices::getdistP2P() { return distP2P; }
getangleProd()46 double CorrespVertices::getangleProd() { return angleProd; }
getActive()47 bool CorrespVertices::getActive() { return Active; }
getEndTriangleActive()48 bool CorrespVertices::getEndTriangleActive() { return EndTriangleActive; }
getIsMaster()49 bool CorrespVertices::getIsMaster() { return IsMaster; }
getTagMaster()50 int CorrespVertices::getTagMaster() { return tagMaster; }
51 
ThinLayer()52 ThinLayer::ThinLayer() {}
53 
~ThinLayer()54 ThinLayer::~ThinLayer() {}
55 
perform()56 void ThinLayer::perform()
57 {
58   ThinLayer::fillVertexToTets();
59   ThinLayer::fillTetToTet4();
60   std::map<MVertex *, double> AllDist = ThinLayer::computeAllDistToOppSide();
61   ThinLayer::checkOppositeTriangles();
62   ThinLayer::fillvecOfThinSheets();
63   std::set<MVertex *> constr_vertices;
64 }
65 
checkOppositeTriangles()66 void ThinLayer::checkOppositeTriangles()
67 {
68   // all endTriangle will be set to active or not
69   for(std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it1 =
70         VertexToCorresp.begin();
71       it1 != VertexToCorresp.end(); it1++) {
72     std::vector<CorrespVertices *> vecCorr = (*it1).second;
73     for(unsigned int i = 0; i < vecCorr.size(); i++) {
74       CorrespVertices *currentCorr = vecCorr[i];
75       faceXtet currentEndTri = currentCorr->getEndTriangle();
76       MVertex *endP0 = currentEndTri.v[0];
77       MVertex *endP1 = currentEndTri.v[1];
78       MVertex *endP2 = currentEndTri.v[2];
79       std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it2 =
80         VertexToCorresp.find(endP0);
81       std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it3 =
82         VertexToCorresp.find(endP1);
83       std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it4 =
84         VertexToCorresp.find(endP2);
85       (*it1).second[i]->setEndTriangleActive(false);
86       if(it2 != VertexToCorresp.end()) {
87         if(it3 != VertexToCorresp.end()) {
88           if(it4 != VertexToCorresp.end()) {
89             if((*it2).second[0]->getActive()) {
90               if((*it3).second[0]->getActive()) {
91                 if((*it4).second[0]->getActive()) {
92                   (*it1).second[i]->setEndTriangleActive(true);
93                 }
94               }
95             }
96           }
97         }
98       }
99     }
100   }
101 }
102 
fillvecOfThinSheets()103 void ThinLayer::fillvecOfThinSheets()
104 {
105   for(std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it1 =
106         VertexToCorresp.begin();
107       it1 != VertexToCorresp.end(); it1++) {
108     std::vector<CorrespVertices *> vecCorr = (*it1).second;
109     for(unsigned int i = 0; i < vecCorr.size(); i++) {
110       CorrespVertices *currentCorr = vecCorr[i];
111       if((currentCorr->getStartPoint()->onWhat()->dim() == 2) &&
112          (currentCorr->getActive()) && (currentCorr->getEndTriangleActive()) &&
113          (currentCorr->getTagMaster() == (-2))) {
114         // Found the first node of a new master sheet
115         std::vector<CorrespVertices *> MasterSheet;
116         MasterSheet.clear();
117         (*it1).second[i]->setTagMaster(-1);
118         faceXtet faceEndSlave = (*it1).second[i]->getEndTriangle();
119         for(unsigned int j = 0; j < 3; j++) {
120           std::map<MVertex *, std::vector<CorrespVertices *> >::iterator it2 =
121             VertexToCorresp.find(faceEndSlave.v[j]);
122           if(it2 != VertexToCorresp.end()) {
123             if(faceEndSlave.v[j]->onWhat()->dim() == 2) {
124               (*it2).second[0]->setTagMaster(
125                 (*it1).second[i]->getStartPoint()->onWhat()->tag());
126             }
127           }
128         }
129         MasterSheet.push_back((*it1).second[i]);
130         std::set<MVertex *> CurrentSheet;
131         CurrentSheet.clear();
132         CurrentSheet.insert((*it1).second[i]->getStartPoint());
133         while(CurrentSheet.size() != 0) {
134           MVertex *VToDo = (*CurrentSheet.begin());
135           std::vector<MTetrahedron *> surroundingTet = VertexToTets[VToDo];
136           for(unsigned int j = 0; j < surroundingTet.size(); j++) {
137             for(std::size_t k = 0; k < surroundingTet[j]->getNumVertices();
138                 k++) {
139               MVertex *ToInsertTmp = surroundingTet[j]->getVertex(k);
140               std::map<MVertex *, std::vector<CorrespVertices *> >::iterator
141                 it2 = VertexToCorresp.find(ToInsertTmp);
142               if(ToInsertTmp->onWhat()->tag() == VToDo->onWhat()->tag()) {
143                 // TODO: OR that onwhat -> dim <, for edges
144                 if(it2 != VertexToCorresp.end()) {
145                   CorrespVertices *correspToInsert = ((*it2).second)[0];
146                   if((correspToInsert->getStartPoint()->onWhat()->dim() == 2) &&
147                      (correspToInsert->getActive()) &&
148                      (correspToInsert->getEndTriangleActive()) &&
149                      (correspToInsert->getTagMaster() == (-2))) {
150                     MasterSheet.push_back((*it2).second[0]);
151                     (*it2).second[0]->setTagMaster(-1);
152                     faceXtet faceEndSlave2 = (*it2).second[0]->getEndTriangle();
153                     for(unsigned int j = 0; j < 3; j++) {
154                       std::map<MVertex *,
155                                std::vector<CorrespVertices *> >::iterator it3 =
156                         VertexToCorresp.find(faceEndSlave2.v[j]);
157                       if(it3 != VertexToCorresp.end()) {
158                         if(faceEndSlave2.v[j]->onWhat()->dim() == 2) {
159                           (*it3).second[0]->setTagMaster(
160                             (*it2).second[i]->getStartPoint()->onWhat()->tag());
161                         }
162                       }
163                     }
164                     CurrentSheet.insert(ToInsertTmp);
165                   }
166                 }
167               }
168             }
169           }
170           CurrentSheet.erase(VToDo);
171         }
172         vecOfThinSheets.push_back(MasterSheet);
173       }
174     }
175   }
176 }
177 
computeAllDistToOppSide()178 std::map<MVertex *, double> ThinLayer::computeAllDistToOppSide()
179 {
180   std::map<MVertex *, double> AllDistToOppSide;
181   GModel *m = GModel::current();
182   // std::vector<MElement*> crackElements;
183   std::set<MVertex *> BoundaryVertices;
184 
185   for(GModel::riter itr = m->firstRegion(); itr != m->lastRegion(); itr++) {
186     GRegion *rTmp = (*itr);
187     for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++) {
188       MTet4 *tet4Tmp = TetToTet4[rTmp->tetrahedra[i]];
189       for(unsigned int j = 0; j < 4; j++) {
190         if(tet4Tmp->getNeigh(j) == 0) {
191           // find the 4th point,and fill the two vector of the boundary triangle
192           faceXtet fxtTmp(tet4Tmp, j);
193           for(int k = 0; k < 3; k++) {
194             MVertex *toTest = fxtTmp.v[k];
195             if(toTest->onWhat()->dim() == 2) {
196               if(BoundaryVertices.find(toTest) == BoundaryVertices.end()) {
197                 BoundaryVertices.insert(toTest);
198               }
199             }
200           }
201           // crackElements.push_back(rTmp->getMeshElement(j));
202         }
203       }
204     }
205   }
206   for(std::set<MVertex *>::iterator it = BoundaryVertices.begin();
207       it != BoundaryVertices.end(); it++) {
208     MVertex *toCompute = (*it);
209     double resultTmp = computeDistToOppSide(toCompute);
210     AllDistToOppSide[toCompute] = resultTmp;
211   }
212 
213   return AllDistToOppSide;
214 }
215 
computeDistToOppSide(MVertex * v)216 double ThinLayer::computeDistToOppSide(MVertex *v)
217 {
218   double DistToOppSide = 0.;
219   // We assume v is on the boundary
220   // First we need to get the internal normal
221   SVector3 InteriorNormal = ThinLayer::computeInteriorNormal(v);
222   // Then we find the first triangle
223   MTet4 *FirstTet = ThinLayer::getTetFromPoint(v, InteriorNormal);
224   MTet4 *CurrentTet = FirstTet;
225   MTet4 *PastTet = FirstTet;
226   SPoint3 CurrentPos = SPoint3(v->x(), v->y(), v->z());
227   SPoint3 LastPos = CurrentPos;
228   int *CurrentTri = 0;
229   CorrespVertices CVTemp;
230   CVTemp.setStartPoint(v);
231   CVTemp.setStartNormal(InteriorNormal);
232   FindNewPoint((&CurrentPos), CurrentTri, CurrentTet, InteriorNormal);
233   faceXtet fxtCV(CurrentTet, (*CurrentTri));
234   //	while(CurrentTet->getNeigh((*CurrentTri)) != 0){
235   while(CurrentTet != 0) {
236     PastTet = CurrentTet;
237     faceXtet fxtCVtmp(PastTet, (*CurrentTri));
238     FindNewPoint((&CurrentPos), CurrentTri, CurrentTet, InteriorNormal);
239     CurrentTet = CurrentTet->getNeigh((*CurrentTri));
240     DistToOppSide += CurrentPos.distance(LastPos);
241     LastPos = CurrentPos;
242   }
243   CVTemp.setEndPoint(LastPos);
244   CVTemp.setEndTriangle(fxtCV);
245   SVector3 EndDir1(fxtCV.v[1]->x() - fxtCV.v[0]->x(),
246                    fxtCV.v[1]->y() - fxtCV.v[0]->y(),
247                    fxtCV.v[1]->z() - fxtCV.v[0]->z());
248   SVector3 EndDir2(fxtCV.v[2]->x() - fxtCV.v[0]->x(),
249                    fxtCV.v[2]->y() - fxtCV.v[0]->y(),
250                    fxtCV.v[2]->z() - fxtCV.v[0]->z());
251   SVector3 EndNormal(EndDir1.y() * EndDir2.z() - EndDir1.z() * EndDir2.y(),
252                      EndDir1.z() * EndDir2.x() - EndDir1.x() * EndDir2.z(),
253                      EndDir1.x() * EndDir2.y() - EndDir1.y() * EndDir2.x());
254   EndNormal.normalize();
255   CVTemp.setEndNormal(EndNormal);
256   CVTemp.setangleProd(
257     fabs(CVTemp.getStartNormal().x() * CVTemp.getEndNormal().x() +
258          CVTemp.getStartNormal().y() * CVTemp.getEndNormal().y() +
259          CVTemp.getStartNormal().z() * CVTemp.getEndNormal().z()));
260   CVTemp.setdistP2P(DistToOppSide);
261   if((CVTemp.getangleProd() > angleMax) && (CVTemp.getdistP2P() < distP2PMax)) {
262     CVTemp.setActive(true);
263   }
264   else {
265     CVTemp.setActive(false);
266   }
267   CVTemp.setTagMaster(-2);
268   VertexToCorresp[v].push_back(&CVTemp);
269   return DistToOppSide;
270 }
271 
computeInteriorNormal(MVertex * v)272 SVector3 ThinLayer::computeInteriorNormal(MVertex *v)
273 {
274   SPoint3 InteriorNormal(0.0, 0.0, 0.0);
275   std::vector<MTetrahedron *> currentVecTet = VertexToTets[v];
276   std::vector<SPoint3> vecInteriorNodes;
277   std::vector<SPoint3> vecFirstDir;
278   std::vector<SPoint3> vecSecondDir;
279   for(unsigned int i = 0; i < currentVecTet.size(); i++) {
280     MTet4 *tet4Tmp = TetToTet4[currentVecTet[i]];
281     for(int j = 0; j < 4; j++) {
282       if(tet4Tmp->getNeigh(j) == 0) {
283         // find the 4th point,and fill the two vector of the boundary triangle
284         faceXtet fxtTmp(tet4Tmp, j);
285         for(int k = 0; k < 4; k++) {
286           bool foundInteriorPoint = true;
287           for(int l = 0; l < 3; l++) {
288             if(tet4Tmp->tet()->getVertex(k) == fxtTmp.v[l]) {
289               foundInteriorPoint = false;
290             }
291           }
292           if(foundInteriorPoint) {
293             SPoint3 pointTmp(tet4Tmp->tet()->getVertex(k)->x(),
294                              tet4Tmp->tet()->getVertex(k)->y(),
295                              tet4Tmp->tet()->getVertex(k)->z());
296             vecInteriorNodes.push_back(pointTmp);
297           }
298         }
299         SPoint3 pointTmp1(fxtTmp.v[1]->x() - fxtTmp.v[0]->x(),
300                           fxtTmp.v[1]->y() - fxtTmp.v[0]->y(),
301                           fxtTmp.v[1]->z() - fxtTmp.v[0]->z());
302         SPoint3 pointTmp2(fxtTmp.v[2]->x() - fxtTmp.v[0]->x(),
303                           fxtTmp.v[2]->y() - fxtTmp.v[0]->y(),
304                           fxtTmp.v[2]->z() - fxtTmp.v[0]->z());
305         vecFirstDir.push_back(pointTmp1);
306         vecSecondDir.push_back(pointTmp2);
307       }
308     }
309   }
310   // at this point we have all the info necessary.
311   SPoint3 pointInteriorAverage(0.0, 0.0, 0.0);
312   for(unsigned int i = 0; i < vecInteriorNodes.size(); i++) {
313     pointInteriorAverage +=
314       SPoint3(vecInteriorNodes[i].x(), vecInteriorNodes[i].y(),
315               vecInteriorNodes[i].z());
316     // pointInteriorAverage.x() += vecInteriorNodes[i].x();
317     // pointInteriorAverage.y() += vecInteriorNodes[i].y();
318     // pointInteriorAverage.z() += vecInteriorNodes[i].z();
319   }
320   pointInteriorAverage =
321     SPoint3(pointInteriorAverage.x() / (double(vecInteriorNodes.size())),
322             pointInteriorAverage.y() / (double(vecInteriorNodes.size())),
323             pointInteriorAverage.z() / (double(vecInteriorNodes.size())));
324   // pointInteriorAverage.x() = pointInteriorAverage.x() /
325   // (double(vecInteriorNodes.size())); pointInteriorAverage.y() =
326   // pointInteriorAverage.y() / (double(vecInteriorNodes.size()));
327   // pointInteriorAverage.z() = pointInteriorAverage.z() /
328   // (double(vecInteriorNodes.size()));
329   SPoint3 dirInteriorAverage(pointInteriorAverage.x() - v->x(),
330                              pointInteriorAverage.y() - v->y(),
331                              pointInteriorAverage.z() - v->z());
332   double norme = sqrt(dirInteriorAverage.x() * dirInteriorAverage.x() +
333                       dirInteriorAverage.y() * dirInteriorAverage.y() +
334                       dirInteriorAverage.z() * dirInteriorAverage.z());
335   dirInteriorAverage =
336     SPoint3(dirInteriorAverage.x() / norme, dirInteriorAverage.y() / norme,
337             dirInteriorAverage.z() / norme);
338   // dirInteriorAverage.x() = dirInteriorAverage.x() / norme;
339   // dirInteriorAverage.y() = dirInteriorAverage.y() / norme;
340   // dirInteriorAverage.z() = dirInteriorAverage.z() / norme;
341   std::vector<SPoint3> vecOrthogDir;
342   for(unsigned int i = 0; i < vecFirstDir.size(); i++) {
343     SPoint3 pointTmp(vecFirstDir[i].y() * vecSecondDir[i].z() -
344                        vecFirstDir[i].z() * vecSecondDir[i].y(),
345                      vecFirstDir[i].z() * vecSecondDir[i].x() -
346                        vecFirstDir[i].x() * vecSecondDir[i].z(),
347                      vecFirstDir[i].x() * vecSecondDir[i].y() -
348                        vecFirstDir[i].y() * vecSecondDir[i].x());
349     vecOrthogDir.push_back(pointTmp);
350   }
351   for(unsigned int i = 0; i < vecOrthogDir.size(); i++) {
352     double prodScalTmp = vecOrthogDir[i].x() * dirInteriorAverage.x() +
353                          vecOrthogDir[i].y() * dirInteriorAverage.y() +
354                          vecOrthogDir[i].z() * dirInteriorAverage.z();
355     if(prodScalTmp < 0.0) {
356       vecOrthogDir[i] = SPoint3(-vecOrthogDir[i].x(), -vecOrthogDir[i].y(),
357                                 -vecOrthogDir[i].z());
358       // vecOrthogDir[i].x() = - vecOrthogDir[i].x();
359       // vecOrthogDir[i].y() = - vecOrthogDir[i].y();
360       // vecOrthogDir[i].z() = - vecOrthogDir[i].z();
361     }
362     double normeTmp = sqrt(vecOrthogDir[i].x() * vecOrthogDir[i].x() +
363                            vecOrthogDir[i].y() * vecOrthogDir[i].y() +
364                            vecOrthogDir[i].z() * vecOrthogDir[i].z());
365     vecOrthogDir[i] =
366       SPoint3(vecOrthogDir[i].x() / normeTmp, vecOrthogDir[i].y() / normeTmp,
367               vecOrthogDir[i].z() / normeTmp);
368     // vecOrthogDir[i].x() = vecOrthogDir[i].x() / normeTmp;
369     // vecOrthogDir[i].y() = vecOrthogDir[i].y() / normeTmp;
370     // vecOrthogDir[i].z() = vecOrthogDir[i].z() / normeTmp;
371     InteriorNormal +=
372       SPoint3(vecOrthogDir[i].x(), vecOrthogDir[i].y(), vecOrthogDir[i].z());
373     // InteriorNormal.x() += vecOrthogDir[i].x();
374     // InteriorNormal.y() += vecOrthogDir[i].y();
375     // InteriorNormal.z() += vecOrthogDir[i].z();
376   }
377   norme = sqrt(InteriorNormal.x() * InteriorNormal.x() +
378                InteriorNormal.y() * InteriorNormal.y() +
379                InteriorNormal.z() * InteriorNormal.z());
380   InteriorNormal =
381     SPoint3(InteriorNormal.x() / norme, InteriorNormal.y() / norme,
382             InteriorNormal.z() / norme);
383   // InteriorNormal.x() = InteriorNormal.x() / norme;
384   // InteriorNormal.y() = InteriorNormal.y() / norme;
385   // InteriorNormal.z() = InteriorNormal.z() / norme;
386   return InteriorNormal;
387 }
388 
getTetFromPoint(MVertex * v,const SVector3 & InteriorNormal)389 MTet4 *ThinLayer::getTetFromPoint(MVertex *v, const SVector3 &InteriorNormal)
390 {
391   MTet4 *TetToGet = 0;
392   std::vector<MTetrahedron *> currentVecTet = VertexToTets[v];
393   for(unsigned int i = 0; i < currentVecTet.size(); i++) {
394     std::vector<SVector3> vecDir;
395     for(int j = 0; j < 4; j++) {
396       if(currentVecTet[i]->getVertex(j) != v) {
397         SVector3 DirTmp(currentVecTet[i]->getVertex(j)->x() - v->x(),
398                         currentVecTet[i]->getVertex(j)->y() - v->y(),
399                         currentVecTet[i]->getVertex(j)->z() - v->z());
400         vecDir.push_back(DirTmp);
401       }
402     }
403     bool IsPositiv =
404       ThinLayer::IsPositivOrientation(vecDir[0], vecDir[1], vecDir[2]);
405     if(!IsPositiv) {
406       SVector3 DirTmp1 = vecDir[1];
407       SVector3 DirTmp2 = vecDir[0];
408       SVector3 DirTmp3 = vecDir[2];
409       vecDir.clear();
410       vecDir.push_back(DirTmp1);
411       vecDir.push_back(DirTmp2);
412       vecDir.push_back(DirTmp3);
413     }
414     bool isPositiv1 =
415       ThinLayer::IsPositivOrientation(vecDir[0], vecDir[1], InteriorNormal);
416     bool isPositiv2 =
417       ThinLayer::IsPositivOrientation(vecDir[1], vecDir[2], InteriorNormal);
418     bool isPositiv3 =
419       ThinLayer::IsPositivOrientation(vecDir[2], vecDir[0], InteriorNormal);
420     if(isPositiv1) {
421       if(isPositiv2) {
422         if(isPositiv3) {
423           TetToGet = TetToTet4[currentVecTet[i]];
424         }
425       }
426     }
427   }
428   return TetToGet;
429 }
430 
IsPositivOrientation(const SVector3 & a,const SVector3 & b,const SVector3 & c)431 bool ThinLayer::IsPositivOrientation(const SVector3 &a, const SVector3 &b,
432                                      const SVector3 &c)
433 {
434   bool result = false;
435   SPoint3 ProdVec(a.y() * b.z() - a.z() * b.y(), a.z() * b.x() - a.x() * b.z(),
436                   a.x() * b.y() - a.y() * b.x());
437   double ProdScal =
438     ProdVec.x() * c.x() + ProdVec.y() * c.y() + ProdVec.z() * c.z();
439   if(ProdScal >= 0.0) {
440     result = true;
441   }
442   return result;
443 }
444 
FindNewPoint(SPoint3 * CurrentPoint,int * CurrentTri,MTet4 * CurrentTet,const SVector3 & InteriorNormal)445 void ThinLayer::FindNewPoint(SPoint3 *CurrentPoint, int *CurrentTri,
446                              MTet4 *CurrentTet, const SVector3 &InteriorNormal)
447 {
448   double distanceP2P = 0.0;
449   double alphaMax = 0.0;
450   double betaMax = 0.0;
451   SPoint3 ResultPoint;
452   int triToGet = 0;
453   for(unsigned int n = 0; n < 4; n++) {
454     // calculer matrice a inverser
455     faceXtet fxt(CurrentTet, n);
456     double a = fxt.v[1]->x() - fxt.v[0]->x();
457     double b = fxt.v[2]->x() - fxt.v[0]->x();
458     double c = InteriorNormal.x();
459     double d = fxt.v[1]->y() - fxt.v[0]->y();
460     double e = fxt.v[2]->y() - fxt.v[0]->y();
461     double f = InteriorNormal.y();
462     double g = fxt.v[1]->z() - fxt.v[0]->z();
463     double h = fxt.v[2]->z() - fxt.v[0]->z();
464     double i = InteriorNormal.z();
465     // produit matrice inverse par vecteur donne poids
466     double detMat =
467       a * e * i + b * f * g + c * d * h - c * e * g - f * h * a - i * b * d;
468     double ai = e * i - f * h;
469     double bi = c * h - b * i;
470     double ci = b * f - c * e;
471     double di = f * g - d * i;
472     double ei = a * i - c * g;
473     double fi = c * d - a * f;
474     // double gi = d * h - e * g;
475     // double hi = b * g - a * h;
476     // double ii = a * e - b * d;
477     double oppx = (*CurrentPoint).x() - fxt.v[0]->x();
478     double oppy = (*CurrentPoint).y() - fxt.v[0]->y();
479     double oppz = (*CurrentPoint).z() - fxt.v[0]->z();
480     double alpha = ai / detMat * oppx + bi / detMat * oppy + ci / detMat * oppz;
481     double beta = di / detMat * oppx + ei / detMat * oppy + fi / detMat * oppz;
482     // double gamma = gi / detMat * oppx + hi / detMat * oppy + ii / detMat *
483     // oppz;
484     // Test si poids entre 0 et 1 et si length maximale
485     if((alpha >= (0.0 - epsilon)) && (alpha <= (1.0 + epsilon))) {
486       if((beta >= (0.0 - epsilon)) && (beta <= (1.0 + epsilon))) {
487         if(((1.0 - alpha - beta) >= (0.0 - epsilon)) &&
488            ((1.0 - alpha - beta) <= (1.0 + epsilon))) {
489           SPoint3 PointTmp(
490             fxt.v[0]->x() + alpha * (fxt.v[1]->x() - fxt.v[0]->x()) +
491               beta * (fxt.v[2]->x() - fxt.v[0]->x()),
492             fxt.v[0]->y() + alpha * (fxt.v[1]->y() - fxt.v[0]->y()) +
493               beta * (fxt.v[2]->y() - fxt.v[0]->y()),
494             fxt.v[0]->z() + alpha * (fxt.v[1]->z() - fxt.v[0]->z()) +
495               beta * (fxt.v[2]->z() - fxt.v[0]->z()));
496           double distanceTmp = PointTmp.distance((*CurrentPoint));
497           if(distanceTmp > distanceP2P) {
498             distanceP2P = distanceTmp;
499             ResultPoint = PointTmp;
500             triToGet = n;
501             alphaMax = alpha;
502             betaMax = beta;
503           }
504         }
505       }
506     }
507   }
508   // test si trop proche d'un point / une arete
509   if(((alphaMax < epsilon) && (betaMax < epsilon)) ||
510      ((alphaMax < epsilon) && ((1.0 - alphaMax - betaMax) < epsilon)) ||
511      (((1.0 - alphaMax - betaMax) < epsilon) && (betaMax < epsilon))) {
512     // proche d'un point
513     double DistMinTmp = 10000000.0;
514     int indexMinTmp = 0;
515     for(unsigned int i = 0; i < 4; i++) {
516       double distanceTmp =
517         sqrt((CurrentTet->tet()->getVertex(i)->x() - ResultPoint.x()) *
518                (CurrentTet->tet()->getVertex(i)->x() - ResultPoint.x()) +
519              (CurrentTet->tet()->getVertex(i)->y() - ResultPoint.y()) *
520                (CurrentTet->tet()->getVertex(i)->y() - ResultPoint.y()) +
521              (CurrentTet->tet()->getVertex(i)->z() - ResultPoint.z()) *
522                (CurrentTet->tet()->getVertex(i)->z() - ResultPoint.z()));
523       if(distanceTmp < DistMinTmp) {
524         DistMinTmp = distanceTmp;
525         indexMinTmp = i;
526       }
527     }
528     MTet4 *NewTet = ThinLayer::getTetFromPoint(
529       CurrentTet->tet()->getVertex(indexMinTmp), InteriorNormal);
530     SPoint3 PointTmp(CurrentTet->tet()->getVertex(indexMinTmp)->x(),
531                      CurrentTet->tet()->getVertex(indexMinTmp)->y(),
532                      CurrentTet->tet()->getVertex(indexMinTmp)->z());
533     (*CurrentPoint) = PointTmp;
534     CurrentTet = NewTet;
535   }
536   else if((alphaMax < epsilon) || (betaMax < epsilon) ||
537           ((1.0 - alphaMax - betaMax) < epsilon)) {
538     // trop proche d'une arete
539   }
540   else {
541     (*CurrentPoint) = ResultPoint;
542     (*CurrentTri) = triToGet;
543     CurrentTet = CurrentTet->getNeigh(triToGet);
544   }
545 }
546 
fillVertexToTets()547 void ThinLayer::fillVertexToTets()
548 {
549   GModel *m = GModel::current();
550   for(GModel::riter itr = m->firstRegion(); itr != m->lastRegion(); itr++) {
551     GRegion *rTmp = (*itr);
552     for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++) {
553       MTetrahedron *elem = rTmp->tetrahedra[i];
554       for(unsigned int j = 0; j < 4; j++) {
555         std::vector<MTetrahedron *> emptyTetVec;
556         emptyTetVec.clear();
557         VertexToTets[elem->getVertex(j)] = emptyTetVec;
558         std::vector<CorrespVertices *> emptyCVVec;
559         emptyCVVec.clear();
560         VertexToCorresp[elem->getVertex(j)] = emptyCVVec;
561       }
562     }
563   }
564   for(GModel::riter itr = m->firstRegion(); itr != m->lastRegion(); itr++) {
565     GRegion *rTmp = (*itr);
566     for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++) {
567       MTetrahedron *elem = rTmp->tetrahedra[i];
568       for(unsigned int j = 0; j < 4; j++) {
569         VertexToTets[elem->getVertex(j)].push_back(elem);
570       }
571     }
572   }
573 }
574 
fillTetToTet4()575 void ThinLayer::fillTetToTet4()
576 {
577   GModel *m = GModel::current();
578   std::vector<MTet4 *> vecAllTet4;
579   vecAllTet4.clear();
580   for(GModel::riter itr = m->firstRegion(); itr != m->lastRegion(); itr++) {
581     GRegion *rTmp = (*itr);
582     for(unsigned int i = 0; i < rTmp->tetrahedra.size(); i++) {
583       MTetrahedron *elem = rTmp->tetrahedra[i];
584       MTet4 tet4tmp = MTet4(elem, 0.0);
585       MTet4 *currentTet4 = &tet4tmp;
586       TetToTet4[elem] = currentTet4;
587       vecAllTet4.clear();
588     }
589   }
590   connectTets(vecAllTet4);
591 }
592