dune-grid-glue 2.9
Loading...
Searching...
No Matches
simplexintersection.cc
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
5
6namespace Dune {
7namespace GridGlue {
8
9template <int dimworld, typename T>
10inline void simplexSubdivision(std::integral_constant<int,0>,const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
11 std::vector<std::vector<unsigned int> >& subElements,
12 std::vector<std::vector<int> >& faceIds);
13template <int dimworld, typename T>
14inline void simplexSubdivision(std::integral_constant<int,1>,const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
15 std::vector<std::vector<unsigned int> >& subElements,
16 std::vector<std::vector<int> >& faceIds);
17template <int dimworld, typename T>
18inline void simplexSubdivision(std::integral_constant<int,2>,const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
19 std::vector<std::vector<unsigned int> >& subElements,
20 std::vector<std::vector<int> >& faceIds);
21template <int dimworld, typename T>
22inline void simplexSubdivision(std::integral_constant<int,3>,const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
23 std::vector<std::vector<unsigned int> >& subElements,
24 std::vector<std::vector<int> >& faceIds);
25
26
27
28// *****************SIMPLEX INTERSECTION COMPUTATION METHODS ***************************
29template<int dimWorld,int dim1, int dim2, typename T>
30class SimplexMethod : public ComputationMethod<dimWorld,dim1,dim2,T>{
31 static_assert(dim1 > dim2, "Specialization missing");
32 friend class ComputationMethod<dimWorld,dim1,dim2,T>;
33public:
34 typedef FieldVector<T, dimWorld> Vector;
35 static const int grid1Dimension = dim1;
36 static const int grid2Dimension = dim2;
37 static const int intersectionDimension = dim2;
38
39 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
40 const std::vector<FieldVector<T,dimWorld> >& Y,
41 std::vector<std::vector<int> > & SX,
42 std::vector<std::vector<int> > & SY,
43 std::vector<FieldVector<T,dimWorld> > & P)
44 {
46 }
47
48 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
49 std::vector<std::vector<unsigned int> >& subElements,
50 std::vector<std::vector<int> >& faceIds)
51 {
52 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
53 elementCorners,subElements, faceIds);
54 }
55
56 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
57 std::vector<std::vector<unsigned int> >& subElements,
58 std::vector<std::vector<int> >& faceIds)
59 {
60 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
61 elementCorners, subElements, faceIds);
62 }
63};
64
65
66
67// POINTS ARE EQUAL
68template<int dimWorld,typename T>
69class SimplexMethod<dimWorld,0,0,T> : public ComputationMethod<dimWorld,0,0,T>{
70 friend class ComputationMethod<dimWorld,0,0,T>;
71
72public:
73 typedef FieldVector<T, dimWorld> Vector;
74 static const int grid1Dimension = 0;
75 static const int grid2Dimension = 0;
76 static const int intersectionDimension = 0;
77
79 const std::vector<FieldVector<T,dimWorld> >& X,
80 const std::vector<FieldVector<T,dimWorld> >& Y,
81 std::vector<std::vector<int> > & SX,
82 std::vector<std::vector<int> > & SY,
83 std::vector<FieldVector<T,dimWorld> > & P)
84 {
85 assert(X.size() == 1 && Y.size() == 1);
86
87 P.clear(); SX.clear(); SY.clear();
88 int k;
89
90 T eps = 1e-8;
91 T a = X[0].infinity_norm();
92 T b = Y[0].infinity_norm();
93 T c = (X[0] - Y[0]).infinity_norm();
94
95 if (c <= eps*a || c <= eps*b ||
96 (a<eps && b< eps && c < 0.5*eps) ) {
97 k = insertPoint(X[0],P);
98
99 return true;
100 }
101 return false;
102 }
103
104 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
105 std::vector<std::vector<unsigned int> >& subElements,
106 std::vector<std::vector<int> >& faceIds)
107 {
108 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
109 elementCorners,subElements, faceIds);
110 }
111
112 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
113 std::vector<std::vector<unsigned int> >& subElements,
114 std::vector<std::vector<int> >& faceIds)
115 {
116 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
117 elementCorners, subElements, faceIds);
118 }
119
120};
121
122
123// POINT ON LINE SEGMENT - :)
124template<int dimWorld,typename T>
125class SimplexMethod<dimWorld,0,1,T> : public ComputationMethod<dimWorld,0,1,T>{
126 friend class ComputationMethod<dimWorld,0,1,T>;
127
128public:
129 typedef FieldVector<T, dimWorld> Vector;
130 static const int grid1Dimension = 0;
131 static const int grid2Dimension = 1;
132 static const int intersectionDimension = 0;
133
134 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
135 const std::vector<FieldVector<T,dimWorld> >& Y,
136 std::vector<std::vector<int> > & SX,
137 std::vector<std::vector<int> > & SY,
138 std::vector<FieldVector<T,dimWorld> > & P)
139 {
140 assert(X.size() == 1 && Y.size() == 2);
141
142 P.clear(); SX.clear(); SY.clear();
143 SY.resize(2);
144
145 if (dimWorld == 1) {
146 T lowerBound = std::max(X[0][0], std::min(Y[0][0],Y[1][0]));
147 T upperBound = std::min(X[0][0], std::max(Y[0][0],Y[1][0]));
148
149 if (lowerBound <= upperBound) { // Intersection is non-empty
150 insertPoint(X[0],P);
151 return true;
152 }
153 } else {
154
155 T eps = 1e-8;
156
157 // check whether the point is on the segment
158 FieldVector<T,dimWorld> v0 = X[0] - Y[0];
159 FieldVector<T,dimWorld> v1 = X[0] - Y[1];
160 FieldVector<T,dimWorld> v2 = Y[1] - Y[0];
161
162 T s = v0.dot(v1);
163 T t = v0.two_norm()/v2.two_norm();
164 v2*=t;
165 v2+=Y[0];
166 v2-=X[0];
167
168 if (v2.infinity_norm() < eps && s<=eps && t<=1+eps) {
169
170 int k = insertPoint(X[0],P);
171 if (s < eps && t < eps)
172 SY[0].push_back(k);
173 else if (s < eps && t>1-eps)
174 SY[1].push_back(k);
175 return true;
176 }
177 }
178 return false;
179 }
180
181 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
182 std::vector<std::vector<unsigned int> >& subElements,
183 std::vector<std::vector<int> >& faceIds)
184 {
185 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
186 elementCorners,subElements, faceIds);
187 }
188
189 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
190 std::vector<std::vector<unsigned int> >& subElements,
191 std::vector<std::vector<int> >& faceIds)
192 {
193 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
194 elementCorners, subElements, faceIds);
195 }
196
197};
198
199
200// POINT IN TRIANGLE - :)
201template<int dimWorld,typename T>
202class SimplexMethod<dimWorld,0,2,T> : public ComputationMethod<dimWorld,0,2,T>{
203 friend class ComputationMethod<dimWorld,0,2,T>;
204
205public:
206 typedef FieldVector<T, dimWorld> Vector;
207 static const int grid1Dimension = 0;
208 static const int grid2Dimension = 2;
209 static const int intersectionDimension = 0;
210
211 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
212 const std::vector<FieldVector<T,dimWorld> >& Y,
213 std::vector<std::vector<int> > & SX,
214 std::vector<std::vector<int> > & SY,
215 std::vector<FieldVector<T,dimWorld> > & P)
216 {
217 assert(X.size() == 1 && Y.size() == 3 && dimWorld > 1);
218
219 P.clear(); SX.clear(); SY.clear();
220 SY.resize(3);
221 int k;
222
223 // If not, check whether it is inside the triangle
224 double eps= 1e-8 ; // tolerance for relative error
225
226 FieldVector<T,dimWorld> v0,v1,v2,r;
227
228 v0 = Y[1] - Y[0];
229 v1 = Y[2] - Y[0];
230 v2 = X[0] - Y[0];
231
232 T s,t,d;
233
234 d = ((v0.dot(v0))*(v1.dot(v1)) - (v0.dot(v1))*(v0.dot(v1)));
235
236 s = ((v1.dot(v1))*(v0.dot(v2)) - (v0.dot(v1))*(v1.dot(v2))) / d;
237 t = ((v0.dot(v0))*(v1.dot(v2)) - (v0.dot(v1))*(v0.dot(v2))) / d;
238
239 v0*=s;
240 v1*=t;
241 r = Y[0] + v0 + v1;
242
243 if (s > -eps && t > -eps && (s+t)< 1+eps && (r-X[0]).infinity_norm() < eps) {
244 k = insertPoint(X[0],P);
245
246 if (t < eps) { // t ~ 0, s varies -> edge 0
247 SY[0].push_back(k);
248 }
249 if (s < eps) { // s ~ 0, t varies -> edge 1
250 SY[1].push_back(k);
251 }
252 if (s+t > 1-eps) { // s+t ~ 1 -> edge 2
253 SY[2].push_back(k);
254 }
255
256 return true;
257 }
258
259 return false;
260 }
261
262 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
263 std::vector<std::vector<unsigned int> >& subElements,
264 std::vector<std::vector<int> >& faceIds)
265 {
266 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
267 elementCorners,subElements, faceIds);
268 }
269
270 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
271 std::vector<std::vector<unsigned int> >& subElements,
272 std::vector<std::vector<int> >& faceIds)
273 {
274 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
275 elementCorners, subElements, faceIds);
276 }
277};
278
279
280// POINT IN TETRAHEDRON - : )
281template<int dimWorld,typename T>
282class SimplexMethod<dimWorld,0,3,T> : public ComputationMethod<dimWorld,0,3,T>{
283 friend class ComputationMethod<dimWorld,0,3,T>;
284
285public:
286 typedef FieldVector<T, dimWorld> Vector;
287 static const int grid1Dimension = 0;
288 static const int grid2Dimension = 3;
289 static const int intersectionDimension = 0;
290
291 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
292 const std::vector<FieldVector<T,dimWorld> >& Y,
293 std::vector<std::vector<int> > & SX,
294 std::vector<std::vector<int> > & SY,
295 std::vector<FieldVector<T,dimWorld> > & P)
296 {
297 assert(X.size() == 1 && Y.size() == 4 && dimWorld == 3);
298
299 P.clear(); SX.clear(); SY.clear();
300 SY.resize(4);
301
302 T eps = 1e-8;
303 // if not, check whether its inside the tetrahedron
304 FieldMatrix<T,dimWorld+1,dimWorld+1> D,DD ;
305
306 D[0][0] = Y[0][0] ; D[0][1] = Y[1][0] ; D[0][2] = Y[2][0] ; D[0][3] = Y[3][0] ;
307 D[1][0] = Y[0][1] ; D[1][1] = Y[1][1] ; D[1][2] = Y[2][1] ; D[1][3] = Y[3][1] ;
308 D[2][0] = Y[0][2] ; D[2][1] = Y[1][2] ; D[2][2] = Y[2][2] ; D[2][3] = Y[3][2] ;
309 D[3][0] = 1 ; D[3][1] = 1 ; D[3][2] = 1 ; D[3][3] = 1 ;
310
311 std::array<T, 5> detD;
312 detD[0] = D.determinant();
313
314 for(unsigned i = 1; i < detD.size(); ++i) {
315 DD = D;
316 for (unsigned d = 0; d < dimWorld; ++d)
317 DD[d][i-1] = X[0][d];
318 detD[i] = DD.determinant();
319 if (std::abs(detD[i]) > eps && std::signbit(detD[0]) != std::signbit(detD[i]))
320 return false; // We are outside.
321 }
322
323 int k = insertPoint(X[0],P);
324 unsigned int faces[4] = {3,2,1,0};
325 for (unsigned i = 1; i < detD.size(); ++i)
326 if(std::abs(detD[i]) < eps)
327 SY[faces[i-1]].push_back(k); // on triangle not containing node i-1
328
329 return true;
330 }
331
332 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
333 std::vector<std::vector<unsigned int> >& subElements,
334 std::vector<std::vector<int> >& faceIds)
335 {
336 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
337 elementCorners,subElements, faceIds);
338 }
339
340 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
341 std::vector<std::vector<unsigned int> >& subElements,
342 std::vector<std::vector<int> >& faceIds)
343 {
344 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
345 elementCorners, subElements, faceIds);
346 }
347};
348
349// SEGEMENT-SEGMENT INTERSECTION - :)
350template<int dimWorld,typename T>
351class SimplexMethod<dimWorld,1,1,T> : public ComputationMethod<dimWorld,1,1,T>{
352 friend class ComputationMethod<dimWorld,1,1,T>;
353
354public:
355 typedef FieldVector<T, dimWorld> Vector;
356 static const int grid1Dimension = 1;
357 static const int grid2Dimension = 1;
358 static const int intersectionDimension = 1;
359
360 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
361 const std::vector<FieldVector<T,dimWorld> >& Y,
362 std::vector<std::vector<int> > & SX,
363 std::vector<std::vector<int> > & SY,
364 std::vector<FieldVector<T,dimWorld> > & P)
365 {
366 assert(X.size() == 2 && Y.size() == 2);
367
368 P.clear(); SX.clear(); SY.clear();
369 SX.resize(2);
370 SY.resize(2);
371 T eps = 1e-8;
372 int k,idX_min=-1,idX_max=-1,idY_min=-1,idY_max=-1;
373
374 // compute intersections
375 switch (dimWorld) {
376 case 1: // d
377 {
378 FieldVector<T,dimWorld> lowerbound(std::max(std::min(X[0][0],X[1][0]), std::min(Y[0][0],X[1][0])));
379 FieldVector<T,dimWorld> upperbound(std::min(std::max(X[0][0],X[1][0]), std::max(Y[0][0],Y[1][0])));
380
381 if (lowerbound[0] < upperbound[0]) { // Intersection is non-empty
382
383 idX_min = (std::min(X[0][0],X[1][0]) < std::min(Y[0][0],Y[1][0]))?(-1):((X[0][0]<X[1][0])?(0):(1));
384 if (idX_min < 0)
385 idY_min = ((Y[0][0]<Y[1][0])?(0):(1));
386
387 idX_max = (std::max(X[0][0],X[1][0]) > std::max(Y[0][0],Y[1][0]))?(-1):((X[0][0]>X[1][0])?(0):(1));
388 if (idX_max < 0)
389 idY_max = ((Y[0][0]>Y[1][0])?(0):(1));
390
391 k = insertPoint(lowerbound,P);
392 if (idX_min >= 0)
393 SX[idX_min].push_back(k);
394 else
395 SY[idY_min].push_back(k);
396
397 k = insertPoint(upperbound,P);
398 if (idX_max >= 0)
399 SX[idX_max].push_back(k);
400 else
401 SY[idY_max].push_back(k);
402
403 return true;
404 }
405 return false;
406 }
407 case 2: // solve X0 + r_0 * (X1 - X0) = Y0 + r_1 * (Y1 - Y0)
408 { // get size_type for all the vectors we are using
409
410 FieldMatrix<T,dimWorld, dimWorld> A;
411 A[0][0] = X[1][0] - X[0][0]; A[0][1] = Y[0][0] - Y[1][0];
412 A[1][0] = X[1][1] - X[0][1]; A[1][1] = Y[0][1] - Y[1][1];
413
414 if (std::abs(A.determinant())>eps) {
415 // lines are non parallel and not degenerated
416 FieldVector<T,dimWorld> p,r,b = Y[0] - X[0];
417 A.solve(r,b) ;
418
419 if ((r[0]>-eps)&&(r[0]<=1+eps)&&(r[1]>-eps)&&(r[1]<1+eps)) {
420 p = X[1] - X[0];
421 p *= r[0] ;
422 p += X[0] ;
423 k = insertPoint(p,P);
424 if(r[0] < eps) { // X = X_0 + r_0 (X_1 - X_0) = X_0
425 SX[0].push_back(k);
426 P[k] = X[0];
427 }
428 else if(r[0] > 1-eps) { // X = X_0 + r_0 (X_1 - X_0) = X_1
429 SX[1].push_back(k);
430 P[k] = X[1];
431 }
432 if(r[1] < eps){ // Y = Y_0 + r_1 (Y_1 - Y_0) = Y_0
433 SY[0].push_back(k);
434 P[k] = Y[0];
435 }
436 else if(r[1] > 1-eps) { // Y = Y_0 + r_1 (Y_1 - Y_0) = Y_1
437 SY[1].push_back(k);
438 P[k] = Y[1];
439 }
440 return true;
441 }
442 } else if ((X[1]-X[0]).infinity_norm() > eps && (Y[1]-Y[0]).infinity_norm() > eps) {
443 // lines are paralles, but non degenerated
444 bool found = false;
445
446 // use triangle equality ||a - b||_2 = || a -c ||_2 + || c - b ||_2 for non perpendicular lines
447 for (unsigned i = 0; i < 2; ++i) {
448 if (std::abs((Y[i]-X[0]).two_norm() + std::abs((Y[i]-X[1]).two_norm())
449 - std::abs((X[1]-X[0]).two_norm())) < eps) {
450 k = insertPoint(Y[i],P);
451 SY[i].push_back(k);
452 found = true;
453 }
454 if (std::abs((X[i]-Y[0]).two_norm() + std::abs((X[i]-Y[1]).two_norm())
455 - std::abs((Y[1]-Y[0]).two_norm())) < eps) {
456 k = insertPoint(X[i],P);
457 SX[i].push_back(k);
458 found = true;
459 }
460 }
461 return found;
462 }
463 return false;
464 }
465 case 3: // solve X0 + r_0 * (X1 - X0) = Y0 + r_1 * (Y1 - Y0)
466 { FieldVector<T,dimWorld> dX, dY, dZ, cXY, cYZ;
467
468 dX = X[1]-X[0];
469 dY = Y[1]-Y[0];
470 dZ = Y[0]-X[0];
471
472 cXY[0] = dX[1]* dY[2] - dX[2]* dY[1];
473 cXY[1] = dX[2]* dY[0] - dX[0]* dY[2];
474 cXY[2] = dX[0]* dY[1] - dX[1]* dY[0];
475
476 if (fabs(dZ.dot(cXY)) < eps*1e+3 && cXY.infinity_norm()>eps) { // coplanar, but not aligned
477
478 cYZ[0] = dY[1]* dZ[2] - dY[2]* dZ[1];
479 cYZ[1] = dY[2]* dZ[0] - dY[0]* dZ[2];
480 cYZ[2] = dY[0]* dZ[1] - dY[1]* dZ[0];
481
482 T s = -cYZ.dot(cXY) / cXY.two_norm2();
483
484 if (s > -eps && s < 1+eps) {
485 dX*= s;
486 dX+= X[0];
487 T o = (dX - Y[0]).two_norm() + (dX- Y[1]).two_norm();
488
489 if (std::abs(o-dY.two_norm()) < eps) {
490 k = insertPoint(dX,P);
491 \
492 if (s<eps) {
493 P[k] = X[0];
494 SX[0].push_back(k);
495 } else if(s > 1-eps) {
496 P[k] = X[1];
497 SX[1].push_back(k);
498 } else if ((dX - Y[0]).two_norm() < eps) {
499 P[k] = Y[0];
500 SY[0].push_back(k);
501 } else if((dX - Y[1]).two_norm() < eps) {
502 P[k] = Y[1];
503 SY[1].push_back(k);
504 }
505
506 return true;
507 }
508 }
509 } else if (cXY.infinity_norm() <= eps) {// lines are parallel
510
511 bool found = false;
512
513 // use triangle equality ||a - b||_2 = || a -c ||_2 + || c - b ||_2,
514 // under the assumption (a-c)*(c-b) > 0 or (a-c) = 0 or (c-b) = 0
515 for (unsigned i = 0; i < 2; ++i) {
516 if ((std::abs((Y[i]-X[0]).two_norm() + std::abs((Y[i]-X[1]).two_norm()) // triangle equality
517 - std::abs((X[1]-X[0]).two_norm())) < eps) &&
518 (std::abs((Y[i]-X[0]).dot((Y[i]-X[1]))) > eps // assumption
519 || (Y[i]-X[0]).infinity_norm() < eps || (Y[i]-X[1]).infinity_norm() < eps)) {
520 k = insertPoint(Y[i],P);
521 SY[i].push_back(k);
522 found = true;
523 }
524 if (std::abs((X[i]-Y[0]).two_norm() + std::abs((X[i]-Y[1]).two_norm()) // triangle equality
525 - std::abs((Y[1]-Y[0]).two_norm())) < eps &&
526 (std::abs((X[i]-Y[0]).dot((X[i]-Y[1]))) > eps // assumption
527 || (X[i]-Y[0]).infinity_norm() < eps || (X[i]-Y[1]).infinity_norm() < eps)){
528 k = insertPoint(X[i],P);
529 SX[i].push_back(k);
530 found = true;
531 }
532 }
533 return found;
534 }
535
536 return false;
537 }
538 }
539 return false;
540 }
541
542 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
543 std::vector<std::vector<unsigned int> >& subElements,
544 std::vector<std::vector<int> >& faceIds)
545 {
546 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
547 elementCorners,subElements, faceIds);
548 }
549
550 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
551 std::vector<std::vector<unsigned int> >& subElements,
552 std::vector<std::vector<int> >& faceIds)
553 {
554 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
555 elementCorners, subElements, faceIds);
556 }
557};
558
559// SEGEMENT-TRIANGLE INTERSECTION - : )
560template<int dimWorld,typename T>
561class SimplexMethod<dimWorld,1,2,T> : public ComputationMethod<dimWorld,1,2,T>{
562 friend class ComputationMethod<dimWorld,1,2,T>;
563
564public:
565 typedef FieldVector<T, dimWorld> Vector;
566 static const int grid1Dimension = 1;
567 static const int grid2Dimension = 2;
568 static const int intersectionDimension = 1;
569
570 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
571 const std::vector<FieldVector<T,dimWorld> >& Y,
572 std::vector<std::vector<int> > & SX,
573 std::vector<std::vector<int> > & SY,
574 std::vector<FieldVector<T,dimWorld> > & P)
575 {
576
577 assert(X.size() == 2 && Y.size() == 3 && dimWorld > 1);
578 P.clear(); SX.clear(); SY.clear();
579 SX.resize(2);
580 SY.resize(3);
581
582 int k;
583 std::vector<FieldVector<T,dimWorld> > surfPts, edge(2), pni(1);
584 std::vector<std::vector<int> > hSX, hSY;
585
586 // is any segment point inside the triangle?
587 for (unsigned ni = 0; ni < 2; ++ni) {
588 pni[0] = X[ni];
589
591 k = insertPoint(X[ni],P);
592 SX[ni].push_back(k);
593 for (unsigned e=0; e < 3; ++e)
594 if (hSY[e].size() > 0)
595 SY[e].push_back(k);
596
597 }
598 surfPts.clear(); hSX.clear(); hSY.clear();
599 }
600
601 if (P.size() >= 2) // we cannot have more than two intersection points
602 return true;
603
604 unsigned int faces[3] = {0,2,1};
605 // do triangle faces intersect with the segment?
606 for (unsigned ni = 0; ni < 3; ++ni) {
607 edge[0] = Y[ni];
608 edge[1] = Y[(ni+1)%3];
609
611 for (unsigned ne = 0; ne < surfPts.size(); ++ne) {
612 k = insertPoint(surfPts[ne],P);
613 SY[faces[ni]].push_back(k);
614 if (hSX[0].size() > 0)
615 SX[0].push_back(k);
616 if (hSX[1].size() > 0)
617 SX[1].push_back(k);
618 }
619
620 if (P.size() >= 2) // we cannot have more than two intersection points
621 return true;
622
623 surfPts.clear(); hSX.clear(); hSY.clear();
624 }
625 }
626
627 if (P.size() >= 2) // we cannot have more than two intersection points
628 return true;
629
630 // if line and triangle are not coplanar in 3d and do not intersect at boundaries
631 if (dimWorld == 3) {
632 T eps = 1e-8;
633
634 Dune::FieldVector<T,dimWorld> B,r,p ;
635 Dune::FieldMatrix<T,dimWorld,dimWorld> A ;
636
637 B = Y[0] - X[0] ;
638
639 for (unsigned i = 0; i < dimWorld; ++i) {
640 A[i][0] = (X[1][i] - X[0][i]);
641 A[i][1] = (Y[0][i] - Y[1][i]);
642 A[i][dimWorld-1] = (Y[0][i] - Y[dimWorld-1][i]);
643 }
644
645 if (std::abs(A.determinant())>eps) {
646
647 A.solve(r,B) ;
648
649 if ((r[0]>=-eps)&&(r[0]<=1+eps)
650 &&(r[1]>=-eps)&&(r[1]<=1+eps)
651 &&(r[2]>=-eps)&&(r[2]<=1+eps)
652 &&(r[1]+r[2]>=-eps) &&(r[1]+r[2]<=1+eps)) {
653 p = X[1] - X[0] ;
654 p *= r[0] ;
655 p += X[0] ;
656 k = insertPoint(p,P);
657
658 if (std::abs(r[0]) < eps) // we prefer exact locations
659 P[k] = X[0];
660 else if (std::abs(r[0]) > 1-eps)
661 P[k] = X[1];
662 else if (std::abs(r[1]) < eps && std::abs(r[2]) < eps)
663 P[k] = Y[0];
664 else if (std::abs(r[1]) < eps && std::abs(r[2]) > 1-eps)
665 P[k] = Y[2];
666 else if (std::abs(r[1]) > 1-eps && std::abs(r[2]) < eps)
667 P[k] = Y[1];
668
669 if (std::abs(r[1]) < eps)
670 SY[1].push_back(k);
671 if (std::fabs(r[dimWorld-1]) < eps)
672 SY[0].push_back(k);
673 if (std::fabs(r[1]+r[dimWorld-1] - 1) < eps)
674 SY[2].push_back(k);
675
676
677 return true ;
678 }
679 }
680 }
681 return false ;
682 }
683
684 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
685 std::vector<std::vector<unsigned int> >& subElements,
686 std::vector<std::vector<int> >& faceIds)
687 {
688 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
689 elementCorners,subElements, faceIds);
690 }
691
692 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
693 std::vector<std::vector<unsigned int> >& subElements,
694 std::vector<std::vector<int> >& faceIds)
695 {
696 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
697 elementCorners, subElements, faceIds);
698 }
699};
700
701// SEGEMENT-TETRAHEDRON INTERSECTION - : )
702template<int dimWorld,typename T>
703class SimplexMethod<dimWorld,1,3,T> : public ComputationMethod<dimWorld,1,3,T>{
704 friend class ComputationMethod<dimWorld,1,3,T>;
705
706public:
707 typedef FieldVector<T, dimWorld> Vector;
708 static const int grid1Dimension = 1;
709 static const int grid2Dimension = 3;
710 static const int intersectionDimension = 1;
711
712 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
713 const std::vector<FieldVector<T,dimWorld> >& Y,
714 std::vector<std::vector<int> > & SX,
715 std::vector<std::vector<int> > & SY,
716 std::vector<FieldVector<T,dimWorld> > & P)
717 {
718 assert(X.size() == 2 && Y.size() == 4 && dimWorld > 2);
719
720 std::vector<int> indices;
721 std::vector<FieldVector<T,dimWorld> > surfPts;
722 std::vector<std::vector<int> > hSX, hSY;
723
724 P.clear(); SX.clear(); SY.clear();
725 SX.resize(2);
726 SY.resize(4);
727
728 std::vector<FieldVector<T,dimWorld> > pni(1);
729
730 bool b = false;
731
732 // check whether the corners of the segment are contained in the tetrahedron
733 for (unsigned int ci = 0; ci < 2; ++ ci) {
734
735 pni[0] = X[ci];
737 int k = insertPoint(X[ci],P);
738 SX[ci].push_back(k);
739 b=true;
740 }
741 surfPts.clear();
742 hSX.clear(); hSY.clear();
743 }
744
745 if (P.size() == 2)
746 return true;
747
748 unsigned int faces[4] = {0,3,2,1};
749 // check whether tetrahedron faces intersect with segment
750 std::vector<FieldVector<T,dimWorld> > triangle(3);
751 for (unsigned int ci = 0; ci < 4; ++ci) { // iterate over all faces
752 triangle[0] = Y[ci];
753 triangle[1] = Y[(ci+1)%4];
754 triangle[2] = Y[(ci+2)%4];
755
756 if (SimplexMethod<dimWorld,1,2,T>::computeIntersectionPoints(X,triangle,hSX,hSY,surfPts)) { // seg - triangle intersection
757 std::vector<int> indices(surfPts.size());
758 for (unsigned int np = 0; np < surfPts.size(); ++np) {
759 int k = insertPoint(surfPts[np],P);
760 indices[np]=k;
761 SY[faces[ci]].push_back(k);
762 }
763
764 // hSX[*] is nonempty if the intersection point is on an edge of the current face of Y
765 for (unsigned int np = 0; np < hSX[0].size(); ++np)
766 SX[0].push_back(indices[hSX[0][np]]);
767 for (unsigned int np = 0; np < hSX[1].size(); ++np)
768 SX[0].push_back(indices[hSX[1][np]]);
769
770 b = true;
771 }
772 surfPts.clear();
773 hSX.clear(); hSY.clear();
774 }
775
776 return b;
777 }
778
779 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
780 std::vector<std::vector<unsigned int> >& subElements,
781 std::vector<std::vector<int> >& faceIds)
782 {
783 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
784 elementCorners,subElements, faceIds);
785 }
786
787 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
788 std::vector<std::vector<unsigned int> >& subElements,
789 std::vector<std::vector<int> >& faceIds)
790 {
791 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
792 elementCorners, subElements, faceIds);
793 }
794};
795
796// TRIANGLE -TRIANGLE INTERSECTION
797template<int dimWorld,typename T>
798class SimplexMethod<dimWorld,2,2,T> : public ComputationMethod<dimWorld,2,2,T>{
799 friend class ComputationMethod<dimWorld,2,2,T>;
800
801public:
802 typedef FieldVector<T, dimWorld> Vector;
803 static const int grid1Dimension = 2;
804 static const int grid2Dimension = 2;
805 static const int intersectionDimension = 2;
806
807 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
808 const std::vector<FieldVector<T,dimWorld> >& Y,
809 std::vector<std::vector<int> > & SX,
810 std::vector<std::vector<int> > & SY,
811 std::vector<FieldVector<T,dimWorld> > & P)
812 {
813 assert(X.size() == 3 && Y.size() == 3 && dimWorld > 1);
814
815 P.clear(); SX.clear(); SY.clear();
816
817 SX.resize(3);
818 SY.resize(3);
819
820 bool b = false;
821
822 std::vector<FieldVector<T,dimWorld> > edge(2);
823 std::vector<FieldVector<T,dimWorld> > surfPts;
824 std::vector<std::vector<int> > hSX, hSY;
825 std::vector<int> indices;
826
827 unsigned int faces[3] = {0,2,1};
828
829 for (unsigned int ni = 0; ni < 3; ++ni) {
830 // check whether the faces of triangle Y intersect the triangle X
831 edge[0] = Y[ni];
832 edge[1] = Y[(ni+1)%3];
833
835
836 indices.resize(surfPts.size());
837 // add intersections of edges of Y with triangle X
838 for (unsigned int np = 0; np < surfPts.size(); ++np) {
839 int k = insertPoint(surfPts[np],P);
840 indices[np] = k;
841 SY[faces[ni]].push_back(k); // add edge data
842 }
843
844 b=true;
845 }
846 if (P.size() >= 6)
847 return true;
848
849 surfPts.clear(); hSX.clear(); hSY.clear();
850 // check whether the faces of triangle X intersect the triangle Y
851 edge[0] = X[ni];
852 edge[1] = X[(ni+1)%3];
853
855
856 indices.resize(surfPts.size());
857 // add intersections of edges of X with triangle Y
858 for (unsigned int np = 0; np < surfPts.size(); ++np) {
859 int k = insertPoint(surfPts[np],P);
860 indices[np] = k;
861 SX[faces[ni]].push_back(k); // add edge data
862 }
863
864 b=true;
865 }
866 if (P.size() >= 6)
867 return true;
868
869 surfPts.clear(); hSX.clear(); hSY.clear();
870 }
871
872 return b;
873 }
874
875 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
876 std::vector<std::vector<unsigned int> >& subElements,
877 std::vector<std::vector<int> >& faceIds)
878 {
879 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
880 elementCorners,subElements, faceIds);
881 }
882
883 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
884 std::vector<std::vector<unsigned int> >& subElements,
885 std::vector<std::vector<int> >& faceIds)
886 {
887 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
888 elementCorners, subElements, faceIds);
889 }
890};
891
892// TRIANGLE -TETRAHEDRON INTERSECTION -: )
893template<int dimWorld,typename T>
894class SimplexMethod<dimWorld,2,3,T> : public ComputationMethod<dimWorld,2,3,T>{
895 friend class ComputationMethod<dimWorld,2,3,T>;
896
897public:
898 typedef FieldVector<T, dimWorld> Vector;
899 static const int grid1Dimension = 2;
900 static const int grid2Dimension = 3;
901 static const int intersectionDimension = 2;
902
903 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
904 const std::vector<FieldVector<T,dimWorld> >& Y,
905 std::vector<std::vector<int> > & SX,
906 std::vector<std::vector<int> > & SY,
907 std::vector<FieldVector<T,dimWorld> > & P)
908 {
909 assert(X.size() == 3 && Y.size() == 4 && dimWorld > 2);
910 P.clear(); SX.clear(); SY.clear();
911 SX.resize(3);
912 SY.resize(4);
913
914 bool b = false;
915 int k,ni,np, ep;
916 std::vector<FieldVector<T,dimWorld> > surfPts, xni(1);
917 std::vector<std::vector<int> > hSX, hSY;
918
919 unsigned int fiX[3][2];
920
921 fiX[0][0] = 0; fiX[1][0] = 0; fiX[2][0] = 1; // faces to node
922 fiX[0][1] = 1; fiX[1][1] = 2; fiX[2][1] = 2;
923 // 1st step: check whether the points of the triangle are contained in the tetrahedron
924 for (ni = 0; ni < 3; ++ni) {
925
926 xni[0] = X[ni];
928 std::vector<int> indices(surfPts.size());
929 for (np = 0; np < surfPts.size(); ++np) {
930 k = insertPoint(X[ni],P);
931 indices[np] = k;
932 SX[fiX[ni][0]].push_back(k); // the corresponding edges to the node are ni and (ni+2)%3
933 SX[fiX[ni][1]].push_back(k);
934 }
935
936 for (ep = 0; ep < 4; ++ep) {
937 for (np = 0; np < hSY[ep].size();++np) {
938 SY[ep].push_back(indices[hSY[ep][np]]);
939 }
940 }
941 b = true;
942 }
943 hSX.clear(); hSY.clear(); surfPts.clear();
944 }
945
946 if (P.size() == 3) // intersection is given by all three corners of the triangle
947 return true;
948
949 // note: points of Y in X is called indirectly via triangle-triangle intesection
950
951 // check whether the triangles of the one tetrahedron intersect the triangle
952 unsigned int facesY[4] = {0,3,2,1}; // face ordering
953
954 std::vector<FieldVector<T,dimWorld> > triangle(3);
955 for (ni = 0; ni < 4; ++ni) {
956
957 triangle[0] = Y[ni];
958 triangle[1] = Y[(ni+1)%4];
959 triangle[2] = Y[(ni+2)%4];
960
961 if (SimplexMethod<dimWorld,2,2,T>::computeIntersectionPoints(X,triangle,hSX,hSY,surfPts)) {
962 std::vector<int> indices(surfPts.size());
963 for (np = 0; np < surfPts.size(); ++np) {
964 k = insertPoint(surfPts[np],P);
965 indices[np]=k;
966 SY[facesY[ni]].push_back(k);
967 }
968
969 // SX[*] is nonempty if the face * of X is intersected
970 for (np = 0; np < 3; ++np) {
971 for (ep = 0; ep < hSX[np].size(); ++ep) {
972 SX[np].push_back(indices[hSX[np][ep]]);
973 }
974 }
975
976 b = true;
977 }
978 hSX.clear(); hSY.clear(); surfPts.clear();
979 }
980 return b;
981 }
982
983 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
984 std::vector<std::vector<unsigned int> >& subElements,
985 std::vector<std::vector<int> >& faceIds)
986 {
987 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
988 elementCorners,subElements, faceIds);
989 }
990
991 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
992 std::vector<std::vector<unsigned int> >& subElements,
993 std::vector<std::vector<int> >& faceIds)
994 {
995 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
996 elementCorners, subElements, faceIds);
997 }
998};
999
1000template<int dimWorld,typename T>
1001class SimplexMethod<dimWorld,3,3,T> : public ComputationMethod<dimWorld,3,3,T>{
1002 friend class ComputationMethod<dimWorld,3,3,T>;
1003
1004public:
1005 typedef FieldVector<T, dimWorld> Vector;
1006 static const int grid1Dimension = 3;
1007 static const int grid2Dimension = 3;
1008 static const int intersectionDimension = 3;
1009
1010 static bool computeIntersectionPoints(const std::vector<FieldVector<T,dimWorld> >& X,
1011 const std::vector<FieldVector<T,dimWorld> >& Y,
1012 std::vector<std::vector<int> > & SX,
1013 std::vector<std::vector<int> > & SY,
1014 std::vector<FieldVector<T,dimWorld> > & P)
1015 {
1016 assert(X.size() == 4 && Y.size() == 4 && dimWorld > 2);
1017 P.clear(); SX.clear(); SY.clear();
1018
1019 SX.resize(4);
1020 SY.resize(4);
1021
1022 bool b = false;
1023 int ci,np,ne,k,ni[4][3];
1024
1025 ni[0][0]= 0 ; ni[0][1]= 1 ; ni[0][2]= 2 ; // faces touching each node
1026 ni[1][0]= 0 ; ni[1][1]= 1 ; ni[1][2]= 3 ;
1027 ni[2][0]= 0 ; ni[2][1]= 2 ; ni[2][2]= 3 ;
1028 ni[3][0]= 1 ; ni[3][1]= 2 ; ni[3][2]= 3 ;
1029
1030 // temporal data container
1031 std::vector<FieldVector<T,dimWorld> > surfPts, pci(1);
1032 std::vector<std::vector<int> > hSX, hSY;
1033
1034 // check whether the points of the one tetrahedron are contained in the other tetrahedron
1035 for (ci = 0; ci < 3; ++ci) {
1036 pci[0] = X[ci];
1037 // point of X is inside Y
1039 std::vector<int> indices(surfPts.size());
1040
1041 for (np = 0; np < surfPts.size(); ++np) {
1042
1043 k = insertPoint(X[ci],P);
1044 indices[np]=k;
1045 SX[ni[ci][0]].push_back(k); // add face data
1046 SX[ni[ci][1]].push_back(k);
1047 SX[ni[ci][2]].push_back(k);
1048 }
1049
1050 // hSY[*] is nonempty if X[ci] is on the face * of Y
1051 for (ne = 0; ne < 4; ++ne) {
1052 for (np = 0; np < hSY[ne].size(); ++np) {
1053 SY[ne].push_back(indices[hSY[ne][np]]);
1054 }
1055 }
1056
1057 b = true;
1058 }
1059 hSX.clear(); hSY.clear(); surfPts.clear();
1060
1061 // probably one point of Y is inside X
1062 surfPts.resize(0);
1063 pci[0]=Y[ci];
1065 std::vector<int> indices(surfPts.size());
1066
1067 for (np = 0; np < surfPts.size(); ++np) {
1068 k = insertPoint(Y[ci],P);
1069 indices[np]=k;
1070 SY[ni[ci][0]].push_back(k); // add face data
1071 SY[ni[ci][1]].push_back(k);
1072 SY[ni[ci][2]].push_back(k);
1073 }
1074
1075 // hSX[*] is nonempty if the point Y[ci] is on the face * of X
1076 for (ne = 0; ne < 4; ++ne) {
1077 for (np = 0; np < hSX[ne].size(); ++np) {
1078 SX[ne].push_back(indices[hSX[ne][np]]);
1079 }
1080 }
1081 b = true;
1082 }
1083 hSX.clear(); hSY.clear(); surfPts.clear();
1084 }
1085
1086 // check whether the triangles of the one tetrahedron intersect the triangles
1087 // of the other tetrahedron
1088 unsigned int faces[4] = {0,3,2,1}; // face ordering
1089
1090 std::vector<FieldVector<T,dimWorld> > triangle(3);
1091 for (ci = 0; ci < 4; ++ci) { // iterate over faces of Y
1092
1093 triangle[0] = Y[ci];
1094 triangle[1] = Y[(ci+1)%4];
1095 triangle[2] = Y[(ci+2)%4];
1096
1097 if(SimplexMethod<dimWorld,3,2,T>::computeIntersectionPoints(X, triangle,hSX,hSY,surfPts)) {
1098
1099 // add Triangle of Y intersects tetrahedron Y data
1100 for (np = 0; np < surfPts.size(); ++np) {
1101 k = insertPoint(surfPts[np],P);
1102 SY[faces[ci]].push_back(k); // add face data
1103 }
1104 b = true;
1105 }
1106 hSX.clear(); hSY.clear(); surfPts.clear();
1107
1108 triangle[0] = X[ci];
1109 triangle[1] = X[(ci+1)%4];
1110 triangle[2] = X[(ci+2)%4];
1111
1112 if(SimplexMethod<dimWorld,3,2,T>::computeIntersectionPoints(Y, triangle,hSY,hSX,surfPts)) {
1113
1114 // add Triangle of Y intersects tetrahedron Y data
1115 for (np = 0; np < surfPts.size(); ++np) {
1116 k = insertPoint(surfPts[np],P);
1117 SX[faces[ci]].push_back(k); // add face data
1118 }
1119 b = true;
1120 }
1121 hSX.clear(); hSY.clear(); surfPts.clear();
1122 }
1123
1124 return b;
1125 }
1126
1127 static void grid1_subdivisions(const std::vector<Vector>& elementCorners,
1128 std::vector<std::vector<unsigned int> >& subElements,
1129 std::vector<std::vector<int> >& faceIds)
1130 {
1131 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
1132 elementCorners,subElements, faceIds);
1133 }
1134
1135 static void grid2_subdivisions(const std::vector<Vector>& elementCorners,
1136 std::vector<std::vector<unsigned int> >& subElements,
1137 std::vector<std::vector<int> >& faceIds)
1138 {
1139 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
1140 elementCorners, subElements, faceIds);
1141 }
1142};
1143
1144template <int dimworld, typename T>
1145inline void simplexSubdivision(std::integral_constant<int,0>,
1146 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1147 std::vector<std::vector<unsigned int> >& subElements,
1148 std::vector<std::vector<int> >& faceIds)
1149{
1150 subElements.resize(1);
1151 faceIds.resize(0);
1152
1153 subElements[0].push_back(0);
1154}
1155
1156template <int dimworld, typename T>
1157inline void simplexSubdivision(std::integral_constant<int,1>,
1158 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1159 std::vector<std::vector<unsigned int> >& subElements,
1160 std::vector<std::vector<int> >& faceIds)
1161{
1162 subElements.resize(1);
1163 faceIds.resize(1);
1164
1165 subElements[0].push_back(0);
1166 subElements[0].push_back(1);
1167
1168 faceIds[0].push_back(0);
1169 faceIds[0].push_back(1);
1170}
1171
1172template <int dimworld, typename T>
1173inline void simplexSubdivision(std::integral_constant<int,2>,
1174 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1175 std::vector<std::vector<unsigned int> >& subElements,
1176 std::vector<std::vector<int> >& faceIds)
1177{
1178 subElements.clear();
1179 faceIds.clear();
1180
1181 if (elementCorners.size() == 3) { // triangle
1182 subElements.resize(1);
1183 faceIds.resize(1);
1184
1185 subElements[0].push_back(0);
1186 subElements[0].push_back(1);
1187 subElements[0].push_back(2);
1188
1189 faceIds[0].push_back(0);
1190 faceIds[0].push_back(1);
1191 faceIds[0].push_back(2);
1192 } else if (elementCorners.size() == 4) { // quadrilateral => 2 triangles
1193 subElements.resize(2);
1194 faceIds.resize(2);
1195
1196 subElements[0].push_back(0);
1197 subElements[0].push_back(1);
1198 subElements[0].push_back(2);
1199
1200 subElements[1].push_back(1);
1201 subElements[1].push_back(2);
1202 subElements[1].push_back(3);
1203
1204 faceIds[0].push_back(2);
1205 faceIds[0].push_back(0);
1206 faceIds[0].push_back(-1);
1207
1208 faceIds[1].push_back(-1);
1209 faceIds[1].push_back(1);
1210 faceIds[1].push_back(3);
1211 }
1212}
1213
1214template <int dimworld, typename T>
1215inline void simplexSubdivision(std::integral_constant<int,3>,
1216 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1217 std::vector<std::vector<unsigned int> >& subElements,
1218 std::vector<std::vector<int> >& faceIds)
1219{
1220 subElements.clear();
1221 faceIds.clear();
1222
1223 if (elementCorners.size() == 4) { // tetrahedron
1224 subElements.resize(1);
1225 faceIds.resize(1);
1226
1227 subElements[0].push_back(0);
1228 subElements[0].push_back(1);
1229 subElements[0].push_back(2);
1230 subElements[0].push_back(3);
1231
1232 faceIds[0].push_back(0);
1233 faceIds[0].push_back(1);
1234 faceIds[0].push_back(2);
1235 faceIds[0].push_back(3);
1236
1237 } else if (elementCorners.size() == 8) { // cube => 5 tetrahedra
1238 subElements.resize(5);
1239 faceIds.resize(5);
1240
1241 subElements[0].push_back(0);
1242 subElements[0].push_back(2);
1243 subElements[0].push_back(3);
1244 subElements[0].push_back(6);
1245
1246 subElements[1].push_back(0);
1247 subElements[1].push_back(1);
1248 subElements[1].push_back(3);
1249 subElements[1].push_back(5);
1250
1251 subElements[2].push_back(0);
1252 subElements[2].push_back(3);
1253 subElements[2].push_back(5);
1254 subElements[2].push_back(6);
1255
1256 subElements[3].push_back(0);
1257 subElements[3].push_back(4);
1258 subElements[3].push_back(5);
1259 subElements[3].push_back(6);
1260
1261 subElements[4].push_back(3);
1262 subElements[4].push_back(5);
1263 subElements[4].push_back(6);
1264 subElements[4].push_back(7);
1265
1266 faceIds[0].push_back(4);
1267 faceIds[0].push_back(0);
1268 faceIds[0].push_back(-1);
1269 faceIds[0].push_back(3);
1270
1271 faceIds[1].push_back(4);
1272 faceIds[1].push_back(2);
1273 faceIds[1].push_back(-1);
1274 faceIds[1].push_back(1);
1275
1276 faceIds[2].push_back(-1);
1277 faceIds[2].push_back(-1);
1278 faceIds[2].push_back(-1);
1279 faceIds[2].push_back(-1);
1280
1281 faceIds[3].push_back(2);
1282 faceIds[3].push_back(0);
1283 faceIds[3].push_back(-1);
1284 faceIds[3].push_back(5);
1285
1286 faceIds[4].push_back(-1);
1287 faceIds[4].push_back(1);
1288 faceIds[4].push_back(3);
1289 faceIds[4].push_back(5);
1290 }
1291}
1292
1293} /* namespace Dune::GridGlue */
1294} /* namespace Dune */
Definition gridglue.hh:37
int insertPoint(const V p, std::vector< V > &P)
Definition computeintersection.hh:164
void simplexSubdivision(std::integral_constant< int, 0 >, const std::vector< Dune::FieldVector< T, dimworld > > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:1145
Definition computeintersection.hh:13
Definition simplexintersection.cc:30
static const int grid1Dimension
Definition simplexintersection.cc:35
static const int grid2Dimension
Definition simplexintersection.cc:36
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:34
static const int intersectionDimension
Definition simplexintersection.cc:37
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:39
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:48
static void grid2_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:56
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:104
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:78
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:73
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:181
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:129
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:134
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:206
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:211
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:262
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:291
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:286
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:332
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:360
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:542
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:355
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:570
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:684
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:565
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:779
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:712
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:707
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:802
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:875
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:807
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:903
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:983
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:898
static void grid1_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:1127
static bool computeIntersectionPoints(const std::vector< FieldVector< T, dimWorld > > &X, const std::vector< FieldVector< T, dimWorld > > &Y, std::vector< std::vector< int > > &SX, std::vector< std::vector< int > > &SY, std::vector< FieldVector< T, dimWorld > > &P)
Definition simplexintersection.cc:1010
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:1005
static void grid2_subdivisions(const std::vector< Vector > &elementCorners, std::vector< std::vector< unsigned int > > &subElements, std::vector< std::vector< int > > &faceIds)
Definition simplexintersection.cc:1135