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);
29template<
int dimWorld,
int dim1,
int dim2,
typename T>
31 static_assert(dim1 > dim2,
"Specialization missing");
34 typedef FieldVector<T, dimWorld>
Vector;
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)
49 std::vector<std::vector<unsigned int> >& subElements,
50 std::vector<std::vector<int> >& faceIds)
52 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
53 elementCorners,subElements, faceIds);
57 std::vector<std::vector<unsigned int> >& subElements,
58 std::vector<std::vector<int> >& faceIds)
60 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
61 elementCorners, subElements, faceIds);
68template<
int dimWorld,
typename T>
73 typedef FieldVector<T, dimWorld>
Vector;
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)
85 assert(X.size() == 1 && Y.size() == 1);
87 P.clear(); SX.clear(); SY.clear();
91 T a = X[0].infinity_norm();
92 T b = Y[0].infinity_norm();
93 T c = (X[0] - Y[0]).infinity_norm();
95 if (c <= eps*a || c <= eps*b ||
96 (a<eps && b< eps && c < 0.5*eps) ) {
105 std::vector<std::vector<unsigned int> >& subElements,
106 std::vector<std::vector<int> >& faceIds)
108 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
109 elementCorners,subElements, faceIds);
113 std::vector<std::vector<unsigned int> >& subElements,
114 std::vector<std::vector<int> >& faceIds)
116 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
117 elementCorners, subElements, faceIds);
124template<
int dimWorld,
typename T>
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)
140 assert(X.size() == 1 && Y.size() == 2);
142 P.clear(); SX.clear(); SY.clear();
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]));
149 if (lowerBound <= upperBound) {
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];
163 T t = v0.two_norm()/v2.two_norm();
168 if (v2.infinity_norm() < eps && s<=eps && t<=1+eps) {
171 if (s < eps && t < eps)
173 else if (s < eps && t>1-eps)
182 std::vector<std::vector<unsigned int> >& subElements,
183 std::vector<std::vector<int> >& faceIds)
185 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
186 elementCorners,subElements, faceIds);
190 std::vector<std::vector<unsigned int> >& subElements,
191 std::vector<std::vector<int> >& faceIds)
193 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
194 elementCorners, subElements, faceIds);
201template<
int dimWorld,
typename T>
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)
217 assert(X.size() == 1 && Y.size() == 3 && dimWorld > 1);
219 P.clear(); SX.clear(); SY.clear();
226 FieldVector<T,dimWorld> v0,v1,v2,r;
234 d = ((v0.dot(v0))*(v1.dot(v1)) - (v0.dot(v1))*(v0.dot(v1)));
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;
243 if (s > -eps && t > -eps && (s+t)< 1+eps && (r-X[0]).infinity_norm() < eps) {
263 std::vector<std::vector<unsigned int> >& subElements,
264 std::vector<std::vector<int> >& faceIds)
266 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
267 elementCorners,subElements, faceIds);
271 std::vector<std::vector<unsigned int> >& subElements,
272 std::vector<std::vector<int> >& faceIds)
274 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
275 elementCorners, subElements, faceIds);
281template<
int dimWorld,
typename T>
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)
297 assert(X.size() == 1 && Y.size() == 4 && dimWorld == 3);
299 P.clear(); SX.clear(); SY.clear();
304 FieldMatrix<T,dimWorld+1,dimWorld+1> D,DD ;
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 ;
311 std::array<T, 5> detD;
312 detD[0] = D.determinant();
314 for(
unsigned i = 1; i < detD.size(); ++i) {
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]))
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);
333 std::vector<std::vector<unsigned int> >& subElements,
334 std::vector<std::vector<int> >& faceIds)
336 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
337 elementCorners,subElements, faceIds);
341 std::vector<std::vector<unsigned int> >& subElements,
342 std::vector<std::vector<int> >& faceIds)
344 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
345 elementCorners, subElements, faceIds);
350template<
int dimWorld,
typename T>
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)
366 assert(X.size() == 2 && Y.size() == 2);
368 P.clear(); SX.clear(); SY.clear();
372 int k,idX_min=-1,idX_max=-1,idY_min=-1,idY_max=-1;
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])));
381 if (lowerbound[0] < upperbound[0]) {
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));
385 idY_min = ((Y[0][0]<Y[1][0])?(0):(1));
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));
389 idY_max = ((Y[0][0]>Y[1][0])?(0):(1));
393 SX[idX_min].push_back(k);
395 SY[idY_min].push_back(k);
399 SX[idX_max].push_back(k);
401 SY[idY_max].push_back(k);
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];
414 if (std::abs(A.determinant())>eps) {
416 FieldVector<T,dimWorld> p,r,b = Y[0] - X[0];
419 if ((r[0]>-eps)&&(r[0]<=1+eps)&&(r[1]>-eps)&&(r[1]<1+eps)) {
428 else if(r[0] > 1-eps) {
436 else if(r[1] > 1-eps) {
442 }
else if ((X[1]-X[0]).infinity_norm() > eps && (Y[1]-Y[0]).infinity_norm() > eps) {
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) {
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) {
466 { FieldVector<T,dimWorld> dX, dY, dZ, cXY, cYZ;
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];
476 if (fabs(dZ.dot(cXY)) < eps*1e+3 && cXY.infinity_norm()>eps) {
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];
482 T s = -cYZ.dot(cXY) / cXY.two_norm2();
484 if (s > -eps && s < 1+eps) {
487 T o = (dX - Y[0]).two_norm() + (dX- Y[1]).two_norm();
489 if (std::abs(o-dY.two_norm()) < eps) {
495 }
else if(s > 1-eps) {
498 }
else if ((dX - Y[0]).two_norm() < eps) {
501 }
else if((dX - Y[1]).two_norm() < eps) {
509 }
else if (cXY.infinity_norm() <= eps) {
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())
517 - std::abs((X[1]-X[0]).two_norm())) < eps) &&
518 (std::abs((Y[i]-X[0]).dot((Y[i]-X[1]))) > eps
519 || (Y[i]-X[0]).infinity_norm() < eps || (Y[i]-X[1]).infinity_norm() < eps)) {
524 if (std::abs((X[i]-Y[0]).two_norm() + std::abs((X[i]-Y[1]).two_norm())
525 - std::abs((Y[1]-Y[0]).two_norm())) < eps &&
526 (std::abs((X[i]-Y[0]).dot((X[i]-Y[1]))) > eps
527 || (X[i]-Y[0]).infinity_norm() < eps || (X[i]-Y[1]).infinity_norm() < eps)){
543 std::vector<std::vector<unsigned int> >& subElements,
544 std::vector<std::vector<int> >& faceIds)
546 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
547 elementCorners,subElements, faceIds);
551 std::vector<std::vector<unsigned int> >& subElements,
552 std::vector<std::vector<int> >& faceIds)
554 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
555 elementCorners, subElements, faceIds);
560template<
int dimWorld,
typename T>
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)
577 assert(X.size() == 2 && Y.size() == 3 && dimWorld > 1);
578 P.clear(); SX.clear(); SY.clear();
583 std::vector<FieldVector<T,dimWorld> > surfPts, edge(2), pni(1);
584 std::vector<std::vector<int> > hSX, hSY;
587 for (
unsigned ni = 0; ni < 2; ++ni) {
593 for (
unsigned e=0; e < 3; ++e)
594 if (hSY[e].size() > 0)
598 surfPts.clear(); hSX.clear(); hSY.clear();
604 unsigned int faces[3] = {0,2,1};
606 for (
unsigned ni = 0; ni < 3; ++ni) {
608 edge[1] = Y[(ni+1)%3];
611 for (
unsigned ne = 0; ne < surfPts.size(); ++ne) {
613 SY[faces[ni]].push_back(k);
614 if (hSX[0].size() > 0)
616 if (hSX[1].size() > 0)
623 surfPts.clear(); hSX.clear(); hSY.clear();
634 Dune::FieldVector<T,dimWorld> B,r,p ;
635 Dune::FieldMatrix<T,dimWorld,dimWorld> A ;
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]);
645 if (std::abs(A.determinant())>eps) {
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)) {
658 if (std::abs(r[0]) < eps)
660 else if (std::abs(r[0]) > 1-eps)
662 else if (std::abs(r[1]) < eps && std::abs(r[2]) < eps)
664 else if (std::abs(r[1]) < eps && std::abs(r[2]) > 1-eps)
666 else if (std::abs(r[1]) > 1-eps && std::abs(r[2]) < eps)
669 if (std::abs(r[1]) < eps)
671 if (std::fabs(r[dimWorld-1]) < eps)
673 if (std::fabs(r[1]+r[dimWorld-1] - 1) < eps)
685 std::vector<std::vector<unsigned int> >& subElements,
686 std::vector<std::vector<int> >& faceIds)
688 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
689 elementCorners,subElements, faceIds);
693 std::vector<std::vector<unsigned int> >& subElements,
694 std::vector<std::vector<int> >& faceIds)
696 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
697 elementCorners, subElements, faceIds);
702template<
int dimWorld,
typename T>
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)
718 assert(X.size() == 2 && Y.size() == 4 && dimWorld > 2);
720 std::vector<int> indices;
721 std::vector<FieldVector<T,dimWorld> > surfPts;
722 std::vector<std::vector<int> > hSX, hSY;
724 P.clear(); SX.clear(); SY.clear();
728 std::vector<FieldVector<T,dimWorld> > pni(1);
733 for (
unsigned int ci = 0; ci < 2; ++ ci) {
742 hSX.clear(); hSY.clear();
748 unsigned int faces[4] = {0,3,2,1};
750 std::vector<FieldVector<T,dimWorld> > triangle(3);
751 for (
unsigned int ci = 0; ci < 4; ++ci) {
753 triangle[1] = Y[(ci+1)%4];
754 triangle[2] = Y[(ci+2)%4];
757 std::vector<int> indices(surfPts.size());
758 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
761 SY[faces[ci]].push_back(k);
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]]);
773 hSX.clear(); hSY.clear();
780 std::vector<std::vector<unsigned int> >& subElements,
781 std::vector<std::vector<int> >& faceIds)
783 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
784 elementCorners,subElements, faceIds);
788 std::vector<std::vector<unsigned int> >& subElements,
789 std::vector<std::vector<int> >& faceIds)
791 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
792 elementCorners, subElements, faceIds);
797template<
int dimWorld,
typename T>
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)
813 assert(X.size() == 3 && Y.size() == 3 && dimWorld > 1);
815 P.clear(); SX.clear(); SY.clear();
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;
827 unsigned int faces[3] = {0,2,1};
829 for (
unsigned int ni = 0; ni < 3; ++ni) {
832 edge[1] = Y[(ni+1)%3];
836 indices.resize(surfPts.size());
838 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
841 SY[faces[ni]].push_back(k);
849 surfPts.clear(); hSX.clear(); hSY.clear();
852 edge[1] = X[(ni+1)%3];
856 indices.resize(surfPts.size());
858 for (
unsigned int np = 0; np < surfPts.size(); ++np) {
861 SX[faces[ni]].push_back(k);
869 surfPts.clear(); hSX.clear(); hSY.clear();
876 std::vector<std::vector<unsigned int> >& subElements,
877 std::vector<std::vector<int> >& faceIds)
879 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
880 elementCorners,subElements, faceIds);
884 std::vector<std::vector<unsigned int> >& subElements,
885 std::vector<std::vector<int> >& faceIds)
887 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
888 elementCorners, subElements, faceIds);
893template<
int dimWorld,
typename T>
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)
909 assert(X.size() == 3 && Y.size() == 4 && dimWorld > 2);
910 P.clear(); SX.clear(); SY.clear();
916 std::vector<FieldVector<T,dimWorld> > surfPts, xni(1);
917 std::vector<std::vector<int> > hSX, hSY;
919 unsigned int fiX[3][2];
921 fiX[0][0] = 0; fiX[1][0] = 0; fiX[2][0] = 1;
922 fiX[0][1] = 1; fiX[1][1] = 2; fiX[2][1] = 2;
924 for (ni = 0; ni < 3; ++ni) {
928 std::vector<int> indices(surfPts.size());
929 for (np = 0; np < surfPts.size(); ++np) {
932 SX[fiX[ni][0]].push_back(k);
933 SX[fiX[ni][1]].push_back(k);
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]]);
943 hSX.clear(); hSY.clear(); surfPts.clear();
952 unsigned int facesY[4] = {0,3,2,1};
954 std::vector<FieldVector<T,dimWorld> > triangle(3);
955 for (ni = 0; ni < 4; ++ni) {
958 triangle[1] = Y[(ni+1)%4];
959 triangle[2] = Y[(ni+2)%4];
962 std::vector<int> indices(surfPts.size());
963 for (np = 0; np < surfPts.size(); ++np) {
966 SY[facesY[ni]].push_back(k);
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]]);
978 hSX.clear(); hSY.clear(); surfPts.clear();
984 std::vector<std::vector<unsigned int> >& subElements,
985 std::vector<std::vector<int> >& faceIds)
987 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
988 elementCorners,subElements, faceIds);
992 std::vector<std::vector<unsigned int> >& subElements,
993 std::vector<std::vector<int> >& faceIds)
995 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
996 elementCorners, subElements, faceIds);
1000template<
int dimWorld,
typename T>
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)
1016 assert(X.size() == 4 && Y.size() == 4 && dimWorld > 2);
1017 P.clear(); SX.clear(); SY.clear();
1023 int ci,np,ne,k,ni[4][3];
1025 ni[0][0]= 0 ; ni[0][1]= 1 ; ni[0][2]= 2 ;
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 ;
1031 std::vector<FieldVector<T,dimWorld> > surfPts, pci(1);
1032 std::vector<std::vector<int> > hSX, hSY;
1035 for (ci = 0; ci < 3; ++ci) {
1039 std::vector<int> indices(surfPts.size());
1041 for (np = 0; np < surfPts.size(); ++np) {
1045 SX[ni[ci][0]].push_back(k);
1046 SX[ni[ci][1]].push_back(k);
1047 SX[ni[ci][2]].push_back(k);
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]]);
1059 hSX.clear(); hSY.clear(); surfPts.clear();
1065 std::vector<int> indices(surfPts.size());
1067 for (np = 0; np < surfPts.size(); ++np) {
1070 SY[ni[ci][0]].push_back(k);
1071 SY[ni[ci][1]].push_back(k);
1072 SY[ni[ci][2]].push_back(k);
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]]);
1083 hSX.clear(); hSY.clear(); surfPts.clear();
1088 unsigned int faces[4] = {0,3,2,1};
1090 std::vector<FieldVector<T,dimWorld> > triangle(3);
1091 for (ci = 0; ci < 4; ++ci) {
1093 triangle[0] = Y[ci];
1094 triangle[1] = Y[(ci+1)%4];
1095 triangle[2] = Y[(ci+2)%4];
1100 for (np = 0; np < surfPts.size(); ++np) {
1102 SY[faces[ci]].push_back(k);
1106 hSX.clear(); hSY.clear(); surfPts.clear();
1108 triangle[0] = X[ci];
1109 triangle[1] = X[(ci+1)%4];
1110 triangle[2] = X[(ci+2)%4];
1115 for (np = 0; np < surfPts.size(); ++np) {
1117 SX[faces[ci]].push_back(k);
1121 hSX.clear(); hSY.clear(); surfPts.clear();
1128 std::vector<std::vector<unsigned int> >& subElements,
1129 std::vector<std::vector<int> >& faceIds)
1131 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid1Dimension>(),
1132 elementCorners,subElements, faceIds);
1136 std::vector<std::vector<unsigned int> >& subElements,
1137 std::vector<std::vector<int> >& faceIds)
1139 simplexSubdivision<dimWorld,T>(std::integral_constant<int,grid2Dimension>(),
1140 elementCorners, subElements, faceIds);
1144template <
int dimworld,
typename T>
1146 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1147 std::vector<std::vector<unsigned int> >& subElements,
1148 std::vector<std::vector<int> >& faceIds)
1150 subElements.resize(1);
1153 subElements[0].push_back(0);
1156template <
int dimworld,
typename T>
1158 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1159 std::vector<std::vector<unsigned int> >& subElements,
1160 std::vector<std::vector<int> >& faceIds)
1162 subElements.resize(1);
1165 subElements[0].push_back(0);
1166 subElements[0].push_back(1);
1168 faceIds[0].push_back(0);
1169 faceIds[0].push_back(1);
1172template <
int dimworld,
typename T>
1174 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1175 std::vector<std::vector<unsigned int> >& subElements,
1176 std::vector<std::vector<int> >& faceIds)
1178 subElements.clear();
1181 if (elementCorners.size() == 3) {
1182 subElements.resize(1);
1185 subElements[0].push_back(0);
1186 subElements[0].push_back(1);
1187 subElements[0].push_back(2);
1189 faceIds[0].push_back(0);
1190 faceIds[0].push_back(1);
1191 faceIds[0].push_back(2);
1192 }
else if (elementCorners.size() == 4) {
1193 subElements.resize(2);
1196 subElements[0].push_back(0);
1197 subElements[0].push_back(1);
1198 subElements[0].push_back(2);
1200 subElements[1].push_back(1);
1201 subElements[1].push_back(2);
1202 subElements[1].push_back(3);
1204 faceIds[0].push_back(2);
1205 faceIds[0].push_back(0);
1206 faceIds[0].push_back(-1);
1208 faceIds[1].push_back(-1);
1209 faceIds[1].push_back(1);
1210 faceIds[1].push_back(3);
1214template <
int dimworld,
typename T>
1216 const std::vector<Dune::FieldVector<T, dimworld> >& elementCorners,
1217 std::vector<std::vector<unsigned int> >& subElements,
1218 std::vector<std::vector<int> >& faceIds)
1220 subElements.clear();
1223 if (elementCorners.size() == 4) {
1224 subElements.resize(1);
1227 subElements[0].push_back(0);
1228 subElements[0].push_back(1);
1229 subElements[0].push_back(2);
1230 subElements[0].push_back(3);
1232 faceIds[0].push_back(0);
1233 faceIds[0].push_back(1);
1234 faceIds[0].push_back(2);
1235 faceIds[0].push_back(3);
1237 }
else if (elementCorners.size() == 8) {
1238 subElements.resize(5);
1241 subElements[0].push_back(0);
1242 subElements[0].push_back(2);
1243 subElements[0].push_back(3);
1244 subElements[0].push_back(6);
1246 subElements[1].push_back(0);
1247 subElements[1].push_back(1);
1248 subElements[1].push_back(3);
1249 subElements[1].push_back(5);
1251 subElements[2].push_back(0);
1252 subElements[2].push_back(3);
1253 subElements[2].push_back(5);
1254 subElements[2].push_back(6);
1256 subElements[3].push_back(0);
1257 subElements[3].push_back(4);
1258 subElements[3].push_back(5);
1259 subElements[3].push_back(6);
1261 subElements[4].push_back(3);
1262 subElements[4].push_back(5);
1263 subElements[4].push_back(6);
1264 subElements[4].push_back(7);
1266 faceIds[0].push_back(4);
1267 faceIds[0].push_back(0);
1268 faceIds[0].push_back(-1);
1269 faceIds[0].push_back(3);
1271 faceIds[1].push_back(4);
1272 faceIds[1].push_back(2);
1273 faceIds[1].push_back(-1);
1274 faceIds[1].push_back(1);
1276 faceIds[2].push_back(-1);
1277 faceIds[2].push_back(-1);
1278 faceIds[2].push_back(-1);
1279 faceIds[2].push_back(-1);
1281 faceIds[3].push_back(2);
1282 faceIds[3].push_back(0);
1283 faceIds[3].push_back(-1);
1284 faceIds[3].push_back(5);
1286 faceIds[4].push_back(-1);
1287 faceIds[4].push_back(1);
1288 faceIds[4].push_back(3);
1289 faceIds[4].push_back(5);
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