3 #ifndef DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
4 #define DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
8 #include <dune/common/array.hh>
9 #include <dune/common/fvector.hh>
10 #include <dune/common/fmatrix.hh>
18 template <
class T,
int dim>
19 static Dune::FieldVector<T,dim>
crossProduct(
const Dune::FieldVector<T,dim>& a,
20 const Dune::FieldVector<T,dim>& b)
23 DUNE_THROW(Dune::NotImplemented,
"crossProduct does not work for dimension " << dim);
25 Dune::FieldVector<T,dim> c;
26 c[0] = a[1]*b[2] - a[2]*b[1];
27 c[1] = a[2]*b[0] - a[0]*b[2];
28 c[2] = a[0]*b[1] - a[1]*b[0];
34 template <
int dim,
int dimworld,
class T=
double>
37 typedef Dune::FieldVector<T,dimworld> WorldCoords;
38 typedef Dune::FieldVector<T,dim> LocalCoords;
53 const std::vector<WorldCoords>& directions,
54 const WorldCoords& target, LocalCoords& preImage,
const T overlap=1e-1) {
56 DUNE_THROW(Dune::NotImplemented,
"inverseProjection is not implemented for dimworld=="<<dimworld
71 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
72 const std::vector<WorldCoords>& targetCorners,
73 LocalCoords& image,
const T overlap=1e-1) {
75 DUNE_THROW(Dune::NotImplemented,
"projection is not implemented for dimworld=="<<dimworld
85 typedef Dune::FieldVector<T,3> WorldCoords;
86 typedef Dune::FieldVector<T,2> LocalCoords;
89 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
90 const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
94 if (targetCorners.size() != 3 && targetCorners.size() != 4)
95 DUNE_THROW(Dune::Exception,
"projection is not implemented for polygone with "<<targetCorners.size()<<
" vertices");
98 if (targetCorners.size() == 4) {
101 std::vector<WorldCoords> triCorners(3);
102 for (
int i=0; i<3; i++)
103 triCorners[i] = targetCorners[i];
105 if (
projection(corner,direction,triCorners,image,overlap))
109 triCorners[0] = targetCorners[3];
111 if (
projection(corner,direction,triCorners,image,overlap)) {
114 image[0] = 1-image[1];
136 Dune::FieldMatrix<T,3, 3> mat(0);
137 for (
int i=0; i<3; i++) {
138 mat[i][0] = targetCorners[1][i]-targetCorners[0][i];
139 mat[i][1] = targetCorners[2][i]-targetCorners[0][i];
140 mat[i][2] = -direction[i];
142 if (std::fabs(mat.determinant())<eps)
146 WorldCoords rhs = corner - targetCorners[0];
154 if (x[0]<-eps || x[1]<-eps || (x[0] + x[1]>1+eps) )
165 const std::vector<WorldCoords>& directions,
166 const WorldCoords& target, LocalCoords& preImage,
167 const T overlap = 1e-1)
171 if (corners.size() != 3 && corners.size() != 4)
172 DUNE_THROW(Dune::NotImplemented,
"inverseProjection is not implemented for elements with "<<corners.size()<<
" corners!");
175 if (corners.size() == 4) {
178 std::vector<WorldCoords> triCorners(3),triDirections(3);
179 for (
int i=0; i<3; i++) {
180 triCorners[i] = corners[i];
181 triDirections[i] = directions[i];
188 triCorners[0] = corners[3];
189 triDirections[0] = directions[3];
194 preImage[0] = 1-preImage[1];
207 WorldCoords n02 = directions[0] - directions[2];
208 WorldCoords n12 = directions[1] - directions[2];
211 T cubic = (directions[2]*n02n12);
215 WorldCoords p02 = corners[0] - corners[2];
216 WorldCoords p12 = corners[1] - corners[2];
217 WorldCoords p2q = corners[2] -target;
221 T quadratic = (directions[2]*p02n12)+ (directions[2]*n02p12)+ (p2q*n02n12);
225 T constant = p2q*p02p12;
228 T linear = (directions[2]*p02p12) + (p2q*n02p12) + (p2q*p02n12);
231 std::vector<T> zeros;
233 if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
235 }
else if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10) {
238 zeros.push_back(-constant/linear);
240 }
else if(std::fabs(cubic)<1e-10) {
243 T p = linear/quadratic;
244 T q = constant/quadratic;
252 zeros.push_back(-0.5*p + std::sqrt(sqt));
253 zeros.push_back(-0.5*p -std::sqrt(sqt));
263 T p= linear - quadratic*quadratic/3;
264 T q=quadratic*(2*quadratic*quadratic/27 - linear/3) + constant;
267 T D = 0.25*q*q + std::pow(p,3)/27;
273 T nu = -q/2+std::sqrt(D);
274 T zer = std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
276 nu = -q/2-std::sqrt(D);
277 zer += std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
279 zeros.push_back(zer-quadratic/3);
281 }
else if (D<-1e-10) {
284 T a = std::sqrt(-4*p/3);
285 T b = std::acos(-0.5*q*std::sqrt(-27/(std::pow(p,3))));
287 for (
int i=0;i<3; i++)
288 zeros.push_back(std::pow(-1,i+1)*a*std::cos((b+(1-i)*M_PI)/3) -quadratic/3);
294 if (std::fabs(q)<1e-10) {
295 zeros.push_back(-quadratic/3);
298 zeros.push_back(std::sqrt(-p)-quadratic/3);
300 }
else if (std::fabs(p)<1e-10) {
302 T nu = std::pow(std::fabs(q),1.0/3.0) * ((q<-eps) ? -1 : 1);
303 zeros.push_back(nu-quadratic/3);
306 zeros.push_back(3*q/p - quadratic/3);
307 zeros.push_back(-1.5*q/p - quadratic/3);
315 for (
size_t i=0;i<zeros.size();i++) {
323 if (nu > zeros[index])
330 WorldCoords e = p2q; e.axpy(nu,directions[2]);
331 WorldCoords f = p12; f.axpy(nu,n12);
332 WorldCoords g = p02; g.axpy(nu,n02);
337 for (
int j=0;j<3; j++) {
341 if (isnan(r[1]) || isinf(r[1]))
344 r[0] = -(e[(j+1)%3]+r[1]*f[(j+1)%3])/g[(j+1)%3];
347 WorldCoords residual = p2q;
348 residual.axpy(r[0],p02);
349 residual.axpy(r[1],p12);
350 residual.axpy(r[2]*r[0],n02);
351 residual.axpy(r[2]*r[1],n12);
352 residual.axpy(r[2],directions[2]);
355 if (!(isnan(r[0]) || isinf(r[0])) && (residual.two_norm()<1e-5))
358 r[0] = -(e[(j+2)%3]+r[1]*f[(j+2)%3])/g[(j+2)%3];
362 residual.axpy(r[0],p02);
363 residual.axpy(r[1],p12);
364 residual.axpy(r[2]*r[0],n02);
365 residual.axpy(r[2]*r[1],n12);
366 residual.axpy(r[2],directions[2]);
368 if (!(isnan(r[0]) || isinf(r[0])) && (residual.two_norm()<1e-5))
373 if (r[0] > -eps && r[1]> -eps && (r[0]+r[1] < 1+eps)) {
380 WorldCoords residual = p2q;
381 residual.axpy(x[0],p02);
382 residual.axpy(x[1],p12);
383 residual.axpy(x[2]*x[0],n02);
384 residual.axpy(x[2]*x[1],n12);
385 residual.axpy(x[2],directions[2]);
387 if (residual.two_norm()<1e-6) {
390 preImage[1] = 1-x[0]-x[1];
403 return inexactInverseProjection(corners,directions, target, preImage,overlap);
412 const std::vector<WorldCoords>& directions,
413 const WorldCoords& target, LocalCoords& preImage,
const T overlap=1e-1)
417 const int nCorners = 3;
418 Dune::FieldVector<T,nCorners> x(1.0/((T) nCorners));
420 for (
int i=0; i<30; i++) {
423 WorldCoords Fxk = target -corners[2];
425 for (
size_t i=0; i<corners.size()-1; i++) {
426 Fxk.axpy(x[i],corners[2]-corners[i]);
427 Fxk.axpy(x[2]*x[i],directions[2]-directions[i]);
429 Fxk.axpy(-x[2],directions[2]);
431 Dune::FieldMatrix<T,nCorners,nCorners> FPrimexk(0);
433 FPrimexk[0] = corners[0] - corners[2];
434 FPrimexk[0].axpy(x[2],directions[0]-directions[2]);
435 FPrimexk[1] = corners[1] - corners[2];
436 FPrimexk[1].axpy(x[2],directions[1]-directions[2]);
437 FPrimexk[2] = directions[2];
438 FPrimexk[2].axpy(x[0],directions[0]-directions[2]);
439 FPrimexk[2].axpy(x[1],directions[1]-directions[2]);
443 WorldCoords newtonCorrection;
445 FPrimexk.mtv(Fxk, newtonCorrection);
447 x += newtonCorrection;
451 if (x[0]>-1e-6 && x[1]>-1e-6 && (x[0]+x[1] <1+1e-6)) {
453 WorldCoords residual = corners[2]-target;
454 residual.axpy(x[0],corners[0]-corners[2]);
455 residual.axpy(x[1],corners[1]-corners[2]);
456 residual.axpy(x[2]*x[0],directions[0]-directions[2]);
457 residual.axpy(x[2]*x[1],directions[1]-directions[2]);
458 residual.axpy(x[2],directions[2]);
461 if (residual.two_norm()>1e-6 || x[2]<-overlap)
465 preImage[1] = 1-x[0]-x[1];
478 typedef Dune::FieldVector<T,2> WorldCoords;
479 typedef Dune::FieldVector<T,1> LocalCoords;
482 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
483 const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
const T overlap=1e-1)
489 Dune::FieldMatrix<T,2,2> mat;
490 mat[0][0] = direction[0];
491 mat[1][0] = direction[1];
492 mat[0][1] = targetCorners[0][0]-targetCorners[1][0];
493 mat[1][1] = targetCorners[0][1]-targetCorners[1][1];
495 WorldCoords rhs = targetCorners[0] - corner;
501 T detinv = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
502 if (std::abs(detinv)<1e-14)
507 x[0] = detinv*(mat[1][1]*rhs[0]-mat[0][1]*rhs[1]);
508 x[1] = detinv*(mat[0][0]*rhs[1]-mat[1][0]*rhs[0]);
512 if (x[1]<-eps || x[1] > 1+eps)
525 const std::vector<WorldCoords>& directions,
526 const WorldCoords& target, LocalCoords& preImage,
527 const T overlap=1e-1)
532 WorldCoords p10 = corners[1]-corners[0];
533 WorldCoords v10 = directions[1]-directions[0];
534 WorldCoords qp0 = target-corners[0];
537 T a = p10[1]*v10[0] - p10[0]*v10[1];
539 T b = -qp0[1]*v10[0] +p10[1]*directions[0][0]
540 + qp0[0]*v10[1] -p10[0]*directions[0][1];
542 T c = -qp0[1]*directions[0][0] +qp0[0]*directions[0][1];
545 if (std::abs(a) < 1e-10) {
548 if (x >= -eps && x <= 1+eps) {
551 T dist = (qp0[0]-x*p10[0])/(directions[0][0]+x*v10[0]);
562 T mu_0 = (-b + std::sqrt(b*b - 4*a*c))/(2*a);
563 T mu_1 = (-b - std::sqrt(b*b - 4*a*c))/(2*a);
565 if (mu_0 >= -eps && mu_0 <= 1+eps) {
567 T dist = (qp0[0]-mu_0*p10[0])/(directions[0][0]+mu_0*v10[0]);
574 }
else if (mu_1 >= -eps && mu_1 <= 1+eps) {
576 T dist = (qp0[0]-mu_1*p10[0])/(directions[0][0]+mu_1*v10[0]);
599 template <
int dim,
int dimworld,
class T>
600 static bool projection(
const Dune::FieldVector<T,dimworld>& corner,
const Dune::FieldVector<T,dimworld>& direction,
601 const std::vector<Dune::FieldVector<T,dimworld> >& targetCorners, Dune::FieldVector<T,dim>& image,
602 const T overlap=1e-1)
605 targetCorners,image,overlap);
619 template <
int dim,
int dimworld,
class T>
621 const std::vector<Dune::FieldVector<T,dimworld> >& directions,
622 const Dune::FieldVector<T,dimworld>& target, Dune::FieldVector<T,dim>& preImage,
623 const T overlap=1e-1)
626 target, preImage, overlap);
629 template <
int dim,
int dimworld,
class T>
633 typedef Dune::FieldVector<T,dimworld> WorldCoords;
634 typedef Dune::FieldVector<T,dim> LocalCoords;
652 const std::vector<WorldCoords>& corners2,
653 const std::vector<WorldCoords>& directions1,
654 const Dune::GeometryType& gt1, Dune::GeometryType& gt2,
655 std::vector<Dune::array<LocalCoords,2> >& polygonCorners,
656 const std::vector<int>& hitCorners, std::bitset<(1<<dim)>& neighborIntersects1,
657 std::bitset<(1<<dim)>& neighborIntersects2,
const T overlap = 1e-1)
659 DUNE_THROW(Dune::NotImplemented,
"addEdgeIntersections is not implemented for dimworld=="<<dimworld
660 <<
" and dim=="<<dim);
664 template <
typename T>
668 typedef Dune::FieldVector<T,2> WorldCoords;
669 typedef Dune::FieldVector<T,1> LocalCoords;
674 const std::vector<WorldCoords>& corners2,
675 const std::vector<WorldCoords>& directions1,
676 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
677 std::vector<Dune::array<LocalCoords,2> >& polygonCorners,
678 const std::vector<int>& hitCorners, std::bitset<2>& neighborIntersects1,
679 std::bitset<2>& neighborIntersects2,
const T overlap=1e-1)
684 template <
typename T>
688 typedef Dune::FieldVector<T,3> WorldCoords;
689 typedef Dune::FieldVector<T,2> LocalCoords;
693 const std::vector<WorldCoords>& corners2,
694 const std::vector<WorldCoords>& directions1,
695 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
696 std::vector<Dune::array<LocalCoords,2> >& polygonCorners,
697 const std::vector<int>& hitCorners, std::bitset<4>& neighborIntersects1,
698 std::bitset<4>& neighborIntersects2,
const T overlap=1e-1)
700 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
701 const Dune::ReferenceElement<T,2>& ref1 = Dune::ReferenceElements<T,2>::general(gt1);
702 const Dune::ReferenceElement<T,2>& ref2 = Dune::ReferenceElements<T,2>::general(gt2);
704 const Dune::GenericReferenceElement<T,2>& ref1 = Dune::GenericReferenceElements<T,2>::general(gt1);
705 const Dune::GenericReferenceElement<T,2>& ref2 = Dune::GenericReferenceElements<T,2>::general(gt2);
709 Dune::array<Dune::FieldVector<T,1>,2> intersection;
710 for (
int i=0; i<ref1.size(1); i++) {
712 std::vector<int> edgeCorners1(2);
713 edgeCorners1[0] = ref1.subEntity(i,1,0,2);
714 edgeCorners1[1] = ref1.subEntity(i,1,1,2);
717 for (
int j=0; j<ref2.size(1); j++) {
719 std::vector<int> edgeCorners2(2);
720 edgeCorners2[0] = ref2.subEntity(j,1,0,2);
721 edgeCorners2[1] = ref2.subEntity(j,1,1,2);
725 if (hitCorners[edgeCorners2[0]]==edgeCorners1[0] || hitCorners[edgeCorners2[0]]==edgeCorners1[1]
726 || hitCorners[edgeCorners2[1]]==edgeCorners1[0] || hitCorners[edgeCorners2[1]]==edgeCorners1[1])
729 if (edgeIntersection(corners1[edgeCorners1[0]],corners1[edgeCorners1[1]],
730 corners2[edgeCorners2[0]],corners2[edgeCorners2[1]],
731 directions1[edgeCorners1[0]],directions1[edgeCorners1[1]],
732 intersection, overlap) )
737 Dune::array<LocalCoords,2> corner;
738 corner[0] = ref1.template geometry<1>(i).global(intersection[0]);
739 corner[1] = ref2.template geometry<1>(j).global(intersection[1]);
740 polygonCorners.push_back(corner);
743 neighborIntersects2[j] =
true;
749 neighborIntersects1[i] =
true;
768 static bool edgeIntersection(
const WorldCoords& corner1,
const WorldCoords& corner2,
769 const WorldCoords& target1,
const WorldCoords& target2,
770 const WorldCoords& direction1,
const WorldCoords& direction2,
771 Dune::array<Dune::FieldVector<T,1>,2>& intersection,
const T overlap = 1e-1)
777 WorldCoords n21 = direction2 - direction1;
778 WorldCoords p21 = corner2 - corner1;
779 WorldCoords q21 = target2 - target1;
782 WorldCoords p1q1 = corner1 -target1;
785 T quadratic = direction1*q21n21;
788 T linear = (direction1*q21p21) + (p1q1*q21n21);
791 T constant = p1q1*q21p21;
794 std::vector<T> zeros;
796 if (std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
798 }
else if (std::fabs(quadratic)<1e-10) {
801 zeros.push_back(-constant/linear);
806 T p = linear/quadratic;
807 T q = constant/quadratic;
815 zeros.push_back(-0.5*p + std::sqrt(sqt));
816 zeros.push_back(-0.5*p -std::sqrt(sqt));
823 for (
size_t i=0;i<zeros.size();i++) {
835 WorldCoords dummy = p1q1; dummy.axpy(eta,direction1);
837 dummy = p21; dummy.axpy(eta,n21);
840 for (
int j=0;j<3; j++) {
843 if (isnan(r[0]) || isinf(r[0]))
846 r[1] = (p1q1[(j+1)%3]+eta*direction1[(j+1)%3] + r[0]*(p21[(j+1)%3]+eta*n21[(j+1)%3]))/q21[(j+1)%3];
849 WorldCoords residual = p1q1;
850 residual.axpy(r[0],p21); residual.axpy(r[2],direction1);
851 residual.axpy(r[2]*r[0],n21); residual.axpy(-r[1],q21);
853 if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
856 residual.axpy(r[1],q21);
857 r[1] = (p1q1[(j+2)%3]+eta*direction1[(j+2)%3] + r[0]*(p21[(j+2)%3]+eta*n21[(j+2)%3]))/q21[(j+2)%3];
859 residual.axpy(r[1],q21);
861 if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
865 if (r[0] >= -eps && r[1]>= -eps && (r[0]<=1+eps) && (r[1] <= 1+eps)) {
872 WorldCoords residual = p1q1;
873 residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
874 residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
876 if (residual.two_norm()<eps)
878 if (index >= 0 && x[0] >= -eps && x[1]>= -eps && (x[0]<=1+eps) && (x[1] <= 1+eps)) {
880 intersection[0] = x[0];
881 intersection[1] = x[1];
893 WorldCoords newtonCorrection(0);
895 for (
int i=0; i<30; i++) {
899 WorldCoords Fxk = target1 -corner1;
900 Fxk.axpy(-x[0], p21);
901 Fxk.axpy(-x[2],direction1);
902 Fxk.axpy(-x[2]*x[0],n21);
905 Dune::FieldMatrix<T,3,3> FPrimexk;
907 FPrimexk[0].axpy(x[2],n21);
908 FPrimexk[1] = target1-target2;
909 FPrimexk[2] = direction1;
910 FPrimexk[2].axpy(x[0],n21);
911 if (FPrimexk.determinant()<1e-10)
916 FPrimexk.mtv(Fxk, newtonCorrection);
918 x += newtonCorrection;
922 residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
923 residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
925 if (residual.two_norm()<=eps) {
927 if (x[0]>=-eps && x[0]<=(1+eps) && x[1]>=-eps && x[1]<=(1+eps)) {
931 intersection[0] = x[0];
932 intersection[1] = x[1];
960 template <
int dim,
int dimworld,
class T=
double>
962 const std::vector<Dune::FieldVector<T,dimworld> >& corners2,
963 const std::vector<Dune::FieldVector<T,dimworld> >& directions1,
964 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
965 std::vector<Dune::array<Dune::FieldVector<T,dim>,2> >& polygonCorners,
966 const std::vector<int>& hitCorners, std::bitset<(1<<dim)>& neighborIntersects1,
967 std::bitset<(1<<dim)>& neighborIntersects2,
const T overlap = 1e-1)
970 directions1, gt1, gt2, polygonCorners,
971 hitCorners, neighborIntersects1,
972 neighborIntersects2, overlap);
static bool projection(const WorldCoords &corner, const WorldCoords &direction, const std::vector< WorldCoords > &targetCorners, LocalCoords &image, const T overlap=1e-1)
Compute the projection of a point along a given direction into the convex hull of some target points...
Definition: projectionhelper.hh:71
static bool inexactInverseProjection(const std::vector< WorldCoords > &corners, const std::vector< WorldCoords > &directions, const WorldCoords &target, LocalCoords &preImage, const T overlap=1e-1)
Compute the inverse projection using Newton's method. This is more stable but as the problem to solve...
Definition: projectionhelper.hh:411
Helper class that provides static methods to compute the projection and inverse projection of a point...
Definition: projectionhelper.hh:35
static bool projection(const Dune::FieldVector< T, dimworld > &corner, const Dune::FieldVector< T, dimworld > &direction, const std::vector< Dune::FieldVector< T, dimworld > > &targetCorners, Dune::FieldVector< T, dim > &image, const T overlap=1e-1)
Compute the projection of a point along a given direction into the convex hull of some target points...
Definition: projectionhelper.hh:600
Helper class that provides static methods for the computation of the intersection of surface element ...
Definition: projectionhelper.hh:630
static bool inverseProjection(const std::vector< Dune::FieldVector< T, dimworld > > &corners, const std::vector< Dune::FieldVector< T, dimworld > > &directions, const Dune::FieldVector< T, dimworld > &target, Dune::FieldVector< T, dim > &preImage, const T overlap=1e-1)
Compute the inverse projection of a point onto some surface element where the projection is done acro...
Definition: projectionhelper.hh:620
static bool inverseProjection(const std::vector< WorldCoords > &corners, const std::vector< WorldCoords > &directions, const WorldCoords &target, LocalCoords &preImage, const T overlap=1e-1)
Definition: projectionhelper.hh:164
static bool projection(const WorldCoords &corner, const WorldCoords &direction, const std::vector< WorldCoords > &targetCorners, LocalCoords &image, const T overlap=1e-1)
Definition: projectionhelper.hh:89
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
Compute the cross product of two vectors.
Definition: projectionhelper.hh:19
static void addEdgeIntersections(const std::vector< Dune::FieldVector< T, dimworld > > &corners1, const std::vector< Dune::FieldVector< T, dimworld > > &corners2, const std::vector< Dune::FieldVector< T, dimworld > > &directions1, const Dune::GeometryType >1, const Dune::GeometryType >2, std::vector< Dune::array< Dune::FieldVector< T, dim >, 2 > > &polygonCorners, const std::vector< int > &hitCorners, std::bitset<(1<< dim)> &neighborIntersects1, std::bitset<(1<< dim)> &neighborIntersects2, const T overlap=1e-1)
Compute the projection along given directions of surface element edges onto target edges...
Definition: projectionhelper.hh:961
This namespace contains helper functions for the projection of two triangular surface on each other...
Definition: projectionhelper.hh:15
static bool projection(const WorldCoords &corner, const WorldCoords &direction, const std::vector< WorldCoords > &targetCorners, LocalCoords &image, const T overlap=1e-1)
Definition: projectionhelper.hh:482
static void addEdgeIntersections(const std::vector< WorldCoords > &corners1, const std::vector< WorldCoords > &corners2, const std::vector< WorldCoords > &directions1, const Dune::GeometryType >1, const Dune::GeometryType >2, std::vector< Dune::array< LocalCoords, 2 > > &polygonCorners, const std::vector< int > &hitCorners, std::bitset< 2 > &neighborIntersects1, std::bitset< 2 > &neighborIntersects2, const T overlap=1e-1)
For 1D surfaces there are no edges.
Definition: projectionhelper.hh:673
static void addEdgeIntersections(const std::vector< WorldCoords > &corners1, const std::vector< WorldCoords > &corners2, const std::vector< WorldCoords > &directions1, const Dune::GeometryType >1, Dune::GeometryType >2, std::vector< Dune::array< LocalCoords, 2 > > &polygonCorners, const std::vector< int > &hitCorners, std::bitset<(1<< dim)> &neighborIntersects1, std::bitset<(1<< dim)> &neighborIntersects2, const T overlap=1e-1)
Compute the projection along given directions of surface element edges onto target edges...
Definition: projectionhelper.hh:651
static bool inverseProjection(const std::vector< WorldCoords > &corners, const std::vector< WorldCoords > &directions, const WorldCoords &target, LocalCoords &preImage, const T overlap=1e-1)
Definition: projectionhelper.hh:524
static void addEdgeIntersections(const std::vector< WorldCoords > &corners1, const std::vector< WorldCoords > &corners2, const std::vector< WorldCoords > &directions1, const Dune::GeometryType >1, const Dune::GeometryType >2, std::vector< Dune::array< LocalCoords, 2 > > &polygonCorners, const std::vector< int > &hitCorners, std::bitset< 4 > &neighborIntersects1, std::bitset< 4 > &neighborIntersects2, const T overlap=1e-1)
Definition: projectionhelper.hh:692
static bool inverseProjection(const std::vector< WorldCoords > &corners, const std::vector< WorldCoords > &directions, const WorldCoords &target, LocalCoords &preImage, const T overlap=1e-1)
Compute the inverse projection of a point onto some surface element where the projection is done acro...
Definition: projectionhelper.hh:52