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