7template <
int dimworld,
typename T>
8inline 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);
11template <
int dimworld,
typename T>
12inline 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);
15template <
int dimworld,
typename T>
16inline 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);
19template <
int dimworld,
typename T>
20inline 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);
27template<
int dimWorld,
int dim1,
int dim2,
typename T>
29 static_assert(dim1 > dim2,
"Specialization missing");
32 typedef FieldVector<T, dimWorld>
Vector;
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)
47 std::vector<std::vector<unsigned int> >& subElements,
48 std::vector<std::vector<int> >& faceIds)
50 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
51 elementCorners,subElements, faceIds);
55 std::vector<std::vector<unsigned int> >& subElements,
56 std::vector<std::vector<int> >& faceIds)
58 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
59 elementCorners, subElements, faceIds);
66template<
int dimWorld,
typename T>
71 typedef FieldVector<T, dimWorld>
Vector;
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)
83 assert(X.size() == 1 && Y.size() == 1);
85 P.clear(); SX.clear(); SY.clear();
89 T a = X[0].infinity_norm();
90 T b = Y[0].infinity_norm();
91 T c = (X[0] - Y[0]).infinity_norm();
93 if (c <= eps*a || c <= eps*b ||
94 (a<eps && b< eps && c < 0.5*eps) ) {
103 std::vector<std::vector<unsigned int> >& subElements,
104 std::vector<std::vector<int> >& faceIds)
106 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
107 elementCorners,subElements, faceIds);
111 std::vector<std::vector<unsigned int> >& subElements,
112 std::vector<std::vector<int> >& faceIds)
114 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
115 elementCorners, subElements, faceIds);
122template<
int dimWorld,
typename T>
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)
138 assert(X.size() == 1 && Y.size() == 2);
140 P.clear(); SX.clear(); SY.clear();
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]));
147 if (lowerBound <= upperBound) {
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];
161 T t = v0.two_norm()/v2.two_norm();
166 if (v2.infinity_norm() < eps && s<=eps && t<=1+eps) {
169 if (s < eps && t < eps)
171 else if (s < eps && t>1-eps)
180 std::vector<std::vector<unsigned int> >& subElements,
181 std::vector<std::vector<int> >& faceIds)
183 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
184 elementCorners,subElements, faceIds);
188 std::vector<std::vector<unsigned int> >& subElements,
189 std::vector<std::vector<int> >& faceIds)
191 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
192 elementCorners, subElements, faceIds);
199template<
int dimWorld,
typename T>
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)
215 assert(X.size() == 1 && Y.size() == 3 && dimWorld > 1);
217 P.clear(); SX.clear(); SY.clear();
224 FieldVector<T,dimWorld> v0,v1,v2,r;
232 d = ((v0.dot(v0))*(v1.dot(v1)) - (v0.dot(v1))*(v0.dot(v1)));
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;
241 if (s > -eps && t > -eps && (s+t)< 1+eps && (r-X[0]).infinity_norm() < eps) {
261 std::vector<std::vector<unsigned int> >& subElements,
262 std::vector<std::vector<int> >& faceIds)
264 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
265 elementCorners,subElements, faceIds);
269 std::vector<std::vector<unsigned int> >& subElements,
270 std::vector<std::vector<int> >& faceIds)
272 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
273 elementCorners, subElements, faceIds);
279template<
int dimWorld,
typename T>
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)
295 assert(X.size() == 1 && Y.size() == 4 && dimWorld == 3);
297 P.clear(); SX.clear(); SY.clear();
302 FieldMatrix<T,dimWorld+1,dimWorld+1> D,DD ;
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 ;
309 std::array<T, 5> detD;
310 detD[0] = D.determinant();
312 for(
unsigned i = 1; i < detD.size(); ++i) {
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]))
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);
331 std::vector<std::vector<unsigned int> >& subElements,
332 std::vector<std::vector<int> >& faceIds)
334 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
335 elementCorners,subElements, faceIds);
339 std::vector<std::vector<unsigned int> >& subElements,
340 std::vector<std::vector<int> >& faceIds)
342 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
343 elementCorners, subElements, faceIds);
348template<
int dimWorld,
typename T>
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)
364 assert(X.size() == 2 && Y.size() == 2);
366 P.clear(); SX.clear(); SY.clear();
370 int k,idX_min=-1,idX_max=-1,idY_min=-1,idY_max=-1;
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])));
379 if (lowerbound[0] < upperbound[0]) {
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));
383 idY_min = ((Y[0][0]<Y[1][0])?(0):(1));
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));
387 idY_max = ((Y[0][0]>Y[1][0])?(0):(1));
391 SX[idX_min].push_back(k);
393 SY[idY_min].push_back(k);
397 SX[idX_max].push_back(k);
399 SY[idY_max].push_back(k);
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];
412 if (std::abs(A.determinant())>eps) {
414 FieldVector<T,dimWorld> p,r,b = Y[0] - X[0];
417 if ((r[0]>-eps)&&(r[0]<=1+eps)&&(r[1]>-eps)&&(r[1]<1+eps)) {
426 else if(r[0] > 1-eps) {
434 else if(r[1] > 1-eps) {
440 }
else if ((X[1]-X[0]).infinity_norm() > eps && (Y[1]-Y[0]).infinity_norm() > eps) {
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) {
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) {
464 { FieldVector<T,dimWorld> dX, dY, dZ, cXY, cYZ;
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];
474 if (fabs(dZ.dot(cXY)) < eps*1e+3 && cXY.infinity_norm()>eps) {
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];
480 T s = -cYZ.dot(cXY) / cXY.two_norm2();
482 if (s > -eps && s < 1+eps) {
485 T o = (dX - Y[0]).two_norm() + (dX- Y[1]).two_norm();
487 if (std::abs(o-dY.two_norm()) < eps) {
493 }
else if(s > 1-eps) {
496 }
else if ((dX - Y[0]).two_norm() < eps) {
499 }
else if((dX - Y[1]).two_norm() < eps) {
507 }
else if (cXY.infinity_norm() <= eps) {
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())
515 - std::abs((X[1]-X[0]).two_norm())) < eps) &&
516 (std::abs((Y[i]-X[0]).dot((Y[i]-X[1]))) > eps
517 || (Y[i]-X[0]).infinity_norm() < eps || (Y[i]-X[1]).infinity_norm() < eps)) {
522 if (std::abs((X[i]-Y[0]).two_norm() + std::abs((X[i]-Y[1]).two_norm())
523 - std::abs((Y[1]-Y[0]).two_norm())) < eps &&
524 (std::abs((X[i]-Y[0]).dot((X[i]-Y[1]))) > eps
525 || (X[i]-Y[0]).infinity_norm() < eps || (X[i]-Y[1]).infinity_norm() < eps)){
541 std::vector<std::vector<unsigned int> >& subElements,
542 std::vector<std::vector<int> >& faceIds)
544 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
545 elementCorners,subElements, faceIds);
549 std::vector<std::vector<unsigned int> >& subElements,
550 std::vector<std::vector<int> >& faceIds)
552 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
553 elementCorners, subElements, faceIds);
558template<
int dimWorld,
typename T>
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)
575 assert(X.size() == 2 && Y.size() == 3 && dimWorld > 1);
576 P.clear(); SX.clear(); SY.clear();
581 std::vector<FieldVector<T,dimWorld> > surfPts, edge(2), pni(1);
582 std::vector<std::vector<int> > hSX, hSY;
585 for (
unsigned ni = 0; ni < 2; ++ni) {
591 for (
unsigned e=0; e < 3; ++e)
592 if (hSY[e].size() > 0)
596 surfPts.clear(); hSX.clear(); hSY.clear();
602 unsigned int faces[3] = {0,2,1};
604 for (
unsigned ni = 0; ni < 3; ++ni) {
606 edge[1] = Y[(ni+1)%3];
609 for (
unsigned ne = 0; ne < surfPts.size(); ++ne) {
611 SY[faces[ni]].push_back(k);
612 if (hSX[0].size() > 0)
614 if (hSX[1].size() > 0)
621 surfPts.clear(); hSX.clear(); hSY.clear();
632 Dune::FieldVector<T,dimWorld> B,r,p ;
633 Dune::FieldMatrix<T,dimWorld,dimWorld> A ;
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]);
643 if (std::abs(A.determinant())>eps) {
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)) {
656 if (std::abs(r[0]) < eps)
658 else if (std::abs(r[0]) > 1-eps)
660 else if (std::abs(r[1]) < eps && std::abs(r[2]) < eps)
662 else if (std::abs(r[1]) < eps && std::abs(r[2]) > 1-eps)
664 else if (std::abs(r[1]) > 1-eps && std::abs(r[2]) < eps)
667 if (std::abs(r[1]) < eps)
669 if (std::fabs(r[dimWorld-1]) < eps)
671 if (std::fabs(r[1]+r[dimWorld-1] - 1) < eps)
683 std::vector<std::vector<unsigned int> >& subElements,
684 std::vector<std::vector<int> >& faceIds)
686 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
687 elementCorners,subElements, faceIds);
691 std::vector<std::vector<unsigned int> >& subElements,
692 std::vector<std::vector<int> >& faceIds)
694 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
695 elementCorners, subElements, faceIds);
700template<
int dimWorld,
typename T>
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)
716 assert(X.size() == 2 && Y.size() == 4 && dimWorld > 2);
718 std::vector<int> indices;
719 std::vector<FieldVector<T,dimWorld> > surfPts;
720 std::vector<std::vector<int> > hSX, hSY;
722 P.clear(); SX.clear(); SY.clear();
726 std::vector<FieldVector<T,dimWorld> > pni(1);
731 for (
unsigned int ci = 0; ci < 2; ++ ci) {
740 hSX.clear(); hSY.clear();
746 unsigned int faces[4] = {0,3,2,1};
748 std::vector<FieldVector<T,dimWorld> > triangle(3);
749 for (
unsigned int ci = 0; ci < 4; ++ci) {
751 triangle[1] = Y[(ci+1)%4];
752 triangle[2] = Y[(ci+2)%4];
755 std::vector<int> indices(surfPts.size());
756 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
759 SY[faces[ci]].push_back(k);
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]]);
771 hSX.clear(); hSY.clear();
778 std::vector<std::vector<unsigned int> >& subElements,
779 std::vector<std::vector<int> >& faceIds)
781 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
782 elementCorners,subElements, faceIds);
786 std::vector<std::vector<unsigned int> >& subElements,
787 std::vector<std::vector<int> >& faceIds)
789 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
790 elementCorners, subElements, faceIds);
795template<
int dimWorld,
typename T>
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)
811 assert(X.size() == 3 && Y.size() == 3 && dimWorld > 1);
813 P.clear(); SX.clear(); SY.clear();
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;
825 unsigned int faces[3] = {0,2,1};
827 for (
unsigned int ni = 0; ni < 3; ++ni) {
830 edge[1] = Y[(ni+1)%3];
834 indices.resize(surfPts.size());
836 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
839 SY[faces[ni]].push_back(k);
847 surfPts.clear(); hSX.clear(); hSY.clear();
850 edge[1] = X[(ni+1)%3];
854 indices.resize(surfPts.size());
856 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
859 SX[faces[ni]].push_back(k);
867 surfPts.clear(); hSX.clear(); hSY.clear();
874 std::vector<std::vector<unsigned int> >& subElements,
875 std::vector<std::vector<int> >& faceIds)
877 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
878 elementCorners,subElements, faceIds);
882 std::vector<std::vector<unsigned int> >& subElements,
883 std::vector<std::vector<int> >& faceIds)
885 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
886 elementCorners, subElements, faceIds);
891template<
int dimWorld,
typename T>
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)
907 assert(X.size() == 3 && Y.size() == 4 && dimWorld > 2);
908 P.clear(); SX.clear(); SY.clear();
914 std::vector<FieldVector<T,dimWorld> > surfPts, xni(1);
915 std::vector<std::vector<int> > hSX, hSY;
917 unsigned int fiX[3][2];
919 fiX[0][0] = 0; fiX[1][0] = 0; fiX[2][0] = 1;
920 fiX[0][1] = 1; fiX[1][1] = 2; fiX[2][1] = 2;
922 for (ni = 0; ni < 3; ++ni) {
926 std::vector<int> indices(surfPts.size());
927 for (np = 0; np < surfPts.size(); ++np) {
930 SX[fiX[ni][0]].push_back(k);
931 SX[fiX[ni][1]].push_back(k);
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]]);
941 hSX.clear(); hSY.clear(); surfPts.clear();
950 unsigned int facesY[4] = {0,3,2,1};
952 std::vector<FieldVector<T,dimWorld> > triangle(3);
953 for (ni = 0; ni < 4; ++ni) {
956 triangle[1] = Y[(ni+1)%4];
957 triangle[2] = Y[(ni+2)%4];
960 std::vector<int> indices(surfPts.size());
961 for (np = 0; np < surfPts.size(); ++np) {
964 SY[facesY[ni]].push_back(k);
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]]);
976 hSX.clear(); hSY.clear(); surfPts.clear();
982 std::vector<std::vector<unsigned int> >& subElements,
983 std::vector<std::vector<int> >& faceIds)
985 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
986 elementCorners,subElements, faceIds);
990 std::vector<std::vector<unsigned int> >& subElements,
991 std::vector<std::vector<int> >& faceIds)
993 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
994 elementCorners, subElements, faceIds);
998template<
int dimWorld,
typename T>
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)
1014 assert(X.size() == 4 && Y.size() == 4 && dimWorld > 2);
1015 P.clear(); SX.clear(); SY.clear();
1021 int ci,np,ne,k,ni[4][3];
1023 ni[0][0]= 0 ; ni[0][1]= 1 ; ni[0][2]= 2 ;
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 ;
1029 std::vector<FieldVector<T,dimWorld> > surfPts, pci(1);
1030 std::vector<std::vector<int> > hSX, hSY;
1033 for (ci = 0; ci < 3; ++ci) {
1037 std::vector<int> indices(surfPts.size());
1039 for (np = 0; np < surfPts.size(); ++np) {
1043 SX[ni[ci][0]].push_back(k);
1044 SX[ni[ci][1]].push_back(k);
1045 SX[ni[ci][2]].push_back(k);
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]]);
1057 hSX.clear(); hSY.clear(); surfPts.clear();
1063 std::vector<int> indices(surfPts.size());
1065 for (np = 0; np < surfPts.size(); ++np) {
1068 SY[ni[ci][0]].push_back(k);
1069 SY[ni[ci][1]].push_back(k);
1070 SY[ni[ci][2]].push_back(k);
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]]);
1081 hSX.clear(); hSY.clear(); surfPts.clear();
1086 unsigned int faces[4] = {0,3,2,1};
1088 std::vector<FieldVector<T,dimWorld> > triangle(3);
1089 for (ci = 0; ci < 4; ++ci) {
1091 triangle[0] = Y[ci];
1092 triangle[1] = Y[(ci+1)%4];
1093 triangle[2] = Y[(ci+2)%4];
1098 for (np = 0; np < surfPts.size(); ++np) {
1100 SY[faces[ci]].push_back(k);
1104 hSX.clear(); hSY.clear(); surfPts.clear();
1106 triangle[0] = X[ci];
1107 triangle[1] = X[(ci+1)%4];
1108 triangle[2] = X[(ci+2)%4];
1113 for (np = 0; np < surfPts.size(); ++np) {
1115 SX[faces[ci]].push_back(k);
1119 hSX.clear(); hSY.clear(); surfPts.clear();
1126 std::vector<std::vector<unsigned int> >& subElements,
1127 std::vector<std::vector<int> >& faceIds)
1129 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
1130 elementCorners,subElements, faceIds);
1134 std::vector<std::vector<unsigned int> >& subElements,
1135 std::vector<std::vector<int> >& faceIds)
1137 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
1138 elementCorners, subElements, faceIds);
1142template <
int dimworld,
typename T>
1144 const std::vector<Dune::FieldVector<T, dimworld> > elementCorners,
1145 std::vector<std::vector<unsigned int> >& subElements,
1146 std::vector<std::vector<int> >& faceIds)
1148 subElements.resize(1);
1151 subElements[0].push_back(0);
1154template <
int dimworld,
typename T>
1156 const std::vector<Dune::FieldVector<T, dimworld> > elementCorners,
1157 std::vector<std::vector<unsigned int> >& subElements,
1158 std::vector<std::vector<int> >& faceIds)
1160 subElements.resize(1);
1163 subElements[0].push_back(0);
1164 subElements[0].push_back(1);
1166 faceIds[0].push_back(0);
1167 faceIds[0].push_back(1);
1170template <
int dimworld,
typename T>
1172 const std::vector<Dune::FieldVector<T, dimworld> > elementCorners,
1173 std::vector<std::vector<unsigned int> >& subElements,
1174 std::vector<std::vector<int> >& faceIds)
1176 subElements.clear();
1179 if (elementCorners.size() == 3) {
1180 subElements.resize(1);
1183 subElements[0].push_back(0);
1184 subElements[0].push_back(1);
1185 subElements[0].push_back(2);
1187 faceIds[0].push_back(0);
1188 faceIds[0].push_back(1);
1189 faceIds[0].push_back(2);
1190 }
else if (elementCorners.size() == 4) {
1191 subElements.resize(2);
1194 subElements[0].push_back(0);
1195 subElements[0].push_back(1);
1196 subElements[0].push_back(2);
1198 subElements[1].push_back(1);
1199 subElements[1].push_back(2);
1200 subElements[1].push_back(3);
1202 faceIds[0].push_back(2);
1203 faceIds[0].push_back(0);
1204 faceIds[0].push_back(-1);
1206 faceIds[1].push_back(-1);
1207 faceIds[1].push_back(1);
1208 faceIds[1].push_back(3);
1212template <
int dimworld,
typename T>
1214 const std::vector<Dune::FieldVector<T, dimworld> > elementCorners,
1215 std::vector<std::vector<unsigned int> >& subElements,
1216 std::vector<std::vector<int> >& faceIds)
1218 subElements.clear();
1221 if (elementCorners.size() == 4) {
1222 subElements.resize(1);
1225 subElements[0].push_back(0);
1226 subElements[0].push_back(1);
1227 subElements[0].push_back(2);
1228 subElements[0].push_back(3);
1230 faceIds[0].push_back(0);
1231 faceIds[0].push_back(1);
1232 faceIds[0].push_back(2);
1233 faceIds[0].push_back(3);
1235 }
else if (elementCorners.size() == 8) {
1236 subElements.resize(5);
1239 subElements[0].push_back(0);
1240 subElements[0].push_back(2);
1241 subElements[0].push_back(3);
1242 subElements[0].push_back(6);
1244 subElements[1].push_back(0);
1245 subElements[1].push_back(1);
1246 subElements[1].push_back(3);
1247 subElements[1].push_back(5);
1249 subElements[2].push_back(0);
1250 subElements[2].push_back(3);
1251 subElements[2].push_back(5);
1252 subElements[2].push_back(6);
1254 subElements[3].push_back(0);
1255 subElements[3].push_back(4);
1256 subElements[3].push_back(5);
1257 subElements[3].push_back(6);
1259 subElements[4].push_back(3);
1260 subElements[4].push_back(5);
1261 subElements[4].push_back(6);
1262 subElements[4].push_back(7);
1264 faceIds[0].push_back(4);
1265 faceIds[0].push_back(0);
1266 faceIds[0].push_back(-1);
1267 faceIds[0].push_back(3);
1269 faceIds[1].push_back(4);
1270 faceIds[1].push_back(2);
1271 faceIds[1].push_back(-1);
1272 faceIds[1].push_back(1);
1274 faceIds[2].push_back(-1);
1275 faceIds[2].push_back(-1);
1276 faceIds[2].push_back(-1);
1277 faceIds[2].push_back(-1);
1279 faceIds[3].push_back(2);
1280 faceIds[3].push_back(0);
1281 faceIds[3].push_back(-1);
1282 faceIds[3].push_back(5);
1284 faceIds[4].push_back(-1);
1285 faceIds[4].push_back(1);
1286 faceIds[4].push_back(3);
1287 faceIds[4].push_back(5);
Definition gridglue.hh:30
int insertPoint(const V p, std::vector< V > &P)
Definition computeintersection.hh:162
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:1143
Definition computeintersection.hh:11
Definition simplexintersection.cc:28
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:54
static const int grid1Dimension
Definition simplexintersection.cc:33
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:37
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:46
static const int grid2Dimension
Definition simplexintersection.cc:34
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:32
static const int intersectionDimension
Definition simplexintersection.cc:35
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:102
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:71
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:76
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:132
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:179
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:127
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:260
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:209
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:204
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:330
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:289
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:284
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:358
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:540
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:353
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:563
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:568
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:682
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:710
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:705
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:777
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:800
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:805
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:873
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:901
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:896
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:981
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:1133
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:1125
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:1008
FieldVector< T, dimWorld > Vector
Definition simplexintersection.cc:1003