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