3 #ifndef DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
4 #define DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
9 #include <dune/common/version.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
13 #include <dune/geometry/referenceelements.hh>
14 #include <dune/geometry/type.hh>
22 template <
class T,
int dim>
23 static Dune::FieldVector<T,dim>
crossProduct(
const Dune::FieldVector<T,dim>& a,
24 const Dune::FieldVector<T,dim>& b)
27 DUNE_THROW(Dune::NotImplemented,
"crossProduct does not work for dimension " << dim);
29 Dune::FieldVector<T,dim> c;
30 c[0] = a[1]*b[2] - a[2]*b[1];
31 c[1] = a[2]*b[0] - a[0]*b[2];
32 c[2] = a[0]*b[1] - a[1]*b[0];
38 template <
int dim,
int dimworld,
class T=
double>
41 typedef Dune::FieldVector<T,dimworld> WorldCoords;
42 typedef Dune::FieldVector<T,dim> LocalCoords;
57 const std::vector<WorldCoords>& directions,
58 const WorldCoords& target, LocalCoords& preImage,
const T overlap=1e-1) {
60 DUNE_THROW(Dune::NotImplemented,
"inverseProjection is not implemented for dimworld=="<<dimworld
75 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
76 const std::vector<WorldCoords>& targetCorners,
77 LocalCoords& image,
const T overlap=1e-1) {
79 DUNE_THROW(Dune::NotImplemented,
"projection is not implemented for dimworld=="<<dimworld
89 typedef Dune::FieldVector<T,3> WorldCoords;
90 typedef Dune::FieldVector<T,2> LocalCoords;
93 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
94 const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
98 if (targetCorners.size() != 3 && targetCorners.size() != 4)
99 DUNE_THROW(Dune::Exception,
"projection is not implemented for polygone with "<<targetCorners.size()<<
" vertices");
102 if (targetCorners.size() == 4) {
105 std::vector<WorldCoords> triCorners(3);
106 for (
int i=0; i<3; i++)
107 triCorners[i] = targetCorners[i];
109 if (
projection(corner,direction,triCorners,image,overlap))
113 triCorners[0] = targetCorners[3];
115 if (
projection(corner,direction,triCorners,image,overlap)) {
118 image[0] = 1-image[1];
141 Dune::FieldMatrix<T,3, 3> mat(0);
142 std::array<WorldCoords, 2> directions;
143 std::array<T, 2> scales;
144 for (
unsigned i = 0; i < 2; ++i) {
145 directions[i] = targetCorners[i+1] - targetCorners[0];
146 scales[i] = directions[i].infinity_norm();
147 directions[i] /= scales[i];
149 for (
int i=0; i<3; i++) {
150 mat[i][0] = directions[0][i];
151 mat[i][1] = directions[1][i];
152 mat[i][2] = -direction[i];
156 WorldCoords rhs = corner - targetCorners[0];
160 for (
unsigned i = 0; i < 2; ++i)
167 if (x[0]<-eps || x[1]<-eps || (x[0] + x[1]>1+eps) )
178 const std::vector<WorldCoords>& directions,
179 const WorldCoords& target, LocalCoords& preImage,
180 const T overlap = 1e-1)
184 if (corners.size() != 3 && corners.size() != 4)
185 DUNE_THROW(Dune::NotImplemented,
"inverseProjection is not implemented for elements with "<<corners.size()<<
" corners!");
188 if (corners.size() == 4) {
191 std::vector<WorldCoords> triCorners(3),triDirections(3);
192 for (
int i=0; i<3; i++) {
193 triCorners[i] = corners[i];
194 triDirections[i] = directions[i];
201 triCorners[0] = corners[3];
202 triDirections[0] = directions[3];
207 preImage[0] = 1-preImage[1];
222 WorldCoords n02 = directions[0] - directions[2];
223 WorldCoords n12 = directions[1] - directions[2];
226 T cubic = (directions[2]*n02n12);
230 WorldCoords p02 = corners[0] - corners[2];
231 WorldCoords p12 = corners[1] - corners[2];
232 WorldCoords p2q = corners[2] -target;
236 T quadratic = (directions[2]*p02n12)+ (directions[2]*n02p12)+ (p2q*n02n12);
240 T constant = p2q*p02p12;
243 T linear = (directions[2]*p02p12) + (p2q*n02p12) + (p2q*p02n12);
246 std::vector<T> zeros;
248 if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
250 }
else if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10) {
253 zeros.push_back(-constant/linear);
255 }
else if(std::fabs(cubic)<1e-10) {
258 T p = linear/quadratic;
259 T q = constant/quadratic;
267 zeros.push_back(-0.5*p + std::sqrt(sqt));
268 zeros.push_back(-0.5*p -std::sqrt(sqt));
278 T p= linear - quadratic*quadratic/3;
279 T q=quadratic*(2*quadratic*quadratic/27 - linear/3) + constant;
282 T D = 0.25*q*q + std::pow(p,3)/27;
288 T nu = -q/2+std::sqrt(D);
289 T zer = std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
291 nu = -q/2-std::sqrt(D);
292 zer += std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
294 zeros.push_back(zer-quadratic/3);
296 }
else if (D<-1e-10) {
299 T a = std::sqrt(-4*p/3);
300 T b = std::acos(-0.5*q*std::sqrt(-27/(std::pow(p,3))));
302 for (
int i=0;i<3; i++)
303 zeros.push_back(std::pow(-1,i+1)*a*std::cos((b+(1-i)*M_PI)/3) -quadratic/3);
309 if (std::fabs(q)<1e-10) {
310 zeros.push_back(-quadratic/3);
313 zeros.push_back(std::sqrt(-p)-quadratic/3);
315 }
else if (std::fabs(p)<1e-10) {
317 T nu = std::pow(std::fabs(q),1.0/3.0) * ((q<-eps) ? -1 : 1);
318 zeros.push_back(nu-quadratic/3);
321 zeros.push_back(3*q/p - quadratic/3);
322 zeros.push_back(-1.5*q/p - quadratic/3);
330 const auto residual = [&](
const WorldCoords& x) -> WorldCoords {
331 WorldCoords residual = p2q;
332 residual.axpy(x[0],p02);
333 residual.axpy(x[1],p12);
334 residual.axpy(x[2]*x[0],n02);
335 residual.axpy(x[2]*x[1],n12);
336 residual.axpy(x[2],directions[2]);
340 for (
size_t i=0;i<zeros.size();i++) {
348 if (nu > zeros[index])
355 WorldCoords e = p2q; e.axpy(nu,directions[2]);
356 WorldCoords f = p12; f.axpy(nu,n12);
357 WorldCoords g = p02; g.axpy(nu,n02);
362 for (
int j=0;j<3; j++) {
366 if (isnan(r[1]) || isinf(r[1]))
369 r[0] = -(e[(j+1)%3]+r[1]*f[(j+1)%3])/g[(j+1)%3];
372 if (!(isnan(r[0]) || isinf(r[0])) && (residual(r).two_norm()<1e-5))
375 r[0] = -(e[(j+2)%3]+r[1]*f[(j+2)%3])/g[(j+2)%3];
378 if (!(isnan(r[0]) || isinf(r[0])) && (residual(r).two_norm()<1e-5))
383 if (r[0] > -eps && r[1]> -eps && (r[0]+r[1] < 1+eps)) {
390 if (residual(x).two_norm()<1e-6) {
393 preImage[1] = 1-x[0]-x[1];
406 return inexactInverseProjection(corners,directions, target, preImage,overlap);
415 const std::vector<WorldCoords>& directions,
416 const WorldCoords& target, LocalCoords& preImage,
const T overlap=1e-1)
418 assert(corners.size() == 3);
419 assert(directions.size() == 3);
422 const int nCorners = 3;
423 Dune::FieldVector<T,nCorners> x(1.0/((T) nCorners));
426 for (std::size_t i = 0; i < corners.size(); ++i)
433 for (
int i=0; i<30; i++) {
436 WorldCoords Fxk = target -corners[2];
438 for (
size_t i=0; i<corners.size()-1; i++) {
439 Fxk.axpy(x[i],corners[2]-corners[i]);
440 Fxk.axpy(x[2]*x[i],directions[2]-directions[i]);
442 Fxk.axpy(-x[2],directions[2]);
444 Dune::FieldMatrix<T,nCorners,nCorners> FPrimexk(0);
446 FPrimexk[0] = corners[0] - corners[2];
447 FPrimexk[0].axpy(x[2],directions[0]-directions[2]);
448 FPrimexk[1] = corners[1] - corners[2];
449 FPrimexk[1].axpy(x[2],directions[1]-directions[2]);
450 FPrimexk[2] = directions[2];
451 FPrimexk[2].axpy(x[0],directions[0]-directions[2]);
452 FPrimexk[2].axpy(x[1],directions[1]-directions[2]);
457 catch (
const Dune::FMatrixError&) {
461 WorldCoords newtonCorrection;
463 FPrimexk.mtv(Fxk, newtonCorrection);
465 x += newtonCorrection;
469 if (x[0]>-1e-6 && x[1]>-1e-6 && (x[0]+x[1] <1+1e-6)) {
471 WorldCoords residual = corners[2]-target;
472 residual.axpy(x[0],corners[0]-corners[2]);
473 residual.axpy(x[1],corners[1]-corners[2]);
474 residual.axpy(x[2]*x[0],directions[0]-directions[2]);
475 residual.axpy(x[2]*x[1],directions[1]-directions[2]);
476 residual.axpy(x[2],directions[2]);
479 if (residual.two_norm()>1e-6 || x[2]<-overlap)
483 preImage[1] = 1-x[0]-x[1];
496 typedef Dune::FieldVector<T,2> WorldCoords;
497 typedef Dune::FieldVector<T,1> LocalCoords;
500 static bool projection(
const WorldCoords& corner,
const WorldCoords& direction,
501 const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
const T overlap=1e-1)
507 Dune::FieldMatrix<T,2,2> mat;
508 mat[0][0] = direction[0];
509 mat[1][0] = direction[1];
510 mat[0][1] = targetCorners[0][0]-targetCorners[1][0];
511 mat[1][1] = targetCorners[0][1]-targetCorners[1][1];
513 WorldCoords rhs = targetCorners[0] - corner;
519 T detinv = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
520 if (std::abs(detinv)<1e-14)
525 x[0] = detinv*(mat[1][1]*rhs[0]-mat[0][1]*rhs[1]);
526 x[1] = detinv*(mat[0][0]*rhs[1]-mat[1][0]*rhs[0]);
530 if (x[1]<-eps || x[1] > 1+eps)
543 const std::vector<WorldCoords>& directions,
544 const WorldCoords& target, LocalCoords& preImage,
545 const T overlap=1e-1)
550 WorldCoords p10 = corners[1]-corners[0];
551 WorldCoords v10 = directions[1]-directions[0];
552 WorldCoords qp0 = target-corners[0];
555 T a = p10[1]*v10[0] - p10[0]*v10[1];
557 T b = -qp0[1]*v10[0] +p10[1]*directions[0][0]
558 + qp0[0]*v10[1] -p10[0]*directions[0][1];
560 T c = -qp0[1]*directions[0][0] +qp0[0]*directions[0][1];
563 if (std::abs(a) < 1e-10) {
566 if (x >= -eps && x <= 1+eps) {
569 T dist = (qp0[0]-x*p10[0])/(directions[0][0]+x*v10[0]);
580 T mu_0 = (-b + std::sqrt(b*b - 4*a*c))/(2*a);
581 T mu_1 = (-b - std::sqrt(b*b - 4*a*c))/(2*a);
583 if (mu_0 >= -eps && mu_0 <= 1+eps) {
585 T dist = (qp0[0]-mu_0*p10[0])/(directions[0][0]+mu_0*v10[0]);
592 }
else if (mu_1 >= -eps && mu_1 <= 1+eps) {
594 T dist = (qp0[0]-mu_1*p10[0])/(directions[0][0]+mu_1*v10[0]);
617 template <
int dim,
int dimworld,
class T>
618 static bool projection(
const Dune::FieldVector<T,dimworld>& corner,
const Dune::FieldVector<T,dimworld>& direction,
619 const std::vector<Dune::FieldVector<T,dimworld> >& targetCorners, Dune::FieldVector<T,dim>& image,
620 const T overlap=1e-1)
623 targetCorners,image,overlap);
637 template <
int dim,
int dimworld,
class T>
639 const std::vector<Dune::FieldVector<T,dimworld> >& directions,
640 const Dune::FieldVector<T,dimworld>& target, Dune::FieldVector<T,dim>& preImage,
641 const T overlap=1e-1)
644 target, preImage, overlap);
647 template <
int dim,
int dimworld,
class T>
651 typedef Dune::FieldVector<T,dimworld> WorldCoords;
652 typedef Dune::FieldVector<T,dim> LocalCoords;
670 const std::vector<WorldCoords>& corners2,
671 const std::vector<WorldCoords>& directions1,
672 const Dune::GeometryType& gt1, Dune::GeometryType& gt2,
673 std::vector<std::array<LocalCoords,2> >& polygonCorners,
674 const std::vector<int>& hitCorners, std::bitset<(1<<dim)>& neighborIntersects1,
675 std::bitset<(1<<dim)>& neighborIntersects2,
const T overlap = 1e-1)
677 DUNE_THROW(Dune::NotImplemented,
"addEdgeIntersections is not implemented for dimworld=="<<dimworld
678 <<
" and dim=="<<dim);
682 template <
typename T>
686 typedef Dune::FieldVector<T,2> WorldCoords;
687 typedef Dune::FieldVector<T,1> LocalCoords;
692 const std::vector<WorldCoords>& corners2,
693 const std::vector<WorldCoords>& directions1,
694 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
695 std::vector<std::array<LocalCoords,2> >& polygonCorners,
696 const std::vector<int>& hitCorners, std::bitset<2>& neighborIntersects1,
697 std::bitset<2>& neighborIntersects2,
const T overlap=1e-1)
702 template <
typename T>
706 typedef Dune::FieldVector<T,3> WorldCoords;
707 typedef Dune::FieldVector<T,2> LocalCoords;
711 const std::vector<WorldCoords>& corners2,
712 const std::vector<WorldCoords>& directions1,
713 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
714 std::vector<std::array<LocalCoords,2> >& polygonCorners,
715 const std::vector<int>& hitCorners, std::bitset<4>& neighborIntersects1,
716 std::bitset<4>& neighborIntersects2,
const T overlap=1e-1)
718 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
719 const Dune::ReferenceElement<T,2>& ref1 = Dune::ReferenceElements<T,2>::general(gt1);
720 const Dune::ReferenceElement<T,2>& ref2 = Dune::ReferenceElements<T,2>::general(gt2);
722 const Dune::GenericReferenceElement<T,2>& ref1 = Dune::GenericReferenceElements<T,2>::general(gt1);
723 const Dune::GenericReferenceElement<T,2>& ref2 = Dune::GenericReferenceElements<T,2>::general(gt2);
727 std::array<Dune::FieldVector<T,1>,2> intersection;
728 for (
int i=0; i<ref1.size(1); i++) {
730 std::array<int, 2> edgeCorners1;
731 edgeCorners1[0] = ref1.subEntity(i,1,0,2);
732 edgeCorners1[1] = ref1.subEntity(i,1,1,2);
735 for (
int j=0; j<ref2.size(1); j++) {
737 std::array<int, 2> edgeCorners2;
738 edgeCorners2[0] = ref2.subEntity(j,1,0,2);
739 edgeCorners2[1] = ref2.subEntity(j,1,1,2);
743 if (hitCorners[edgeCorners2[0]]==edgeCorners1[0] || hitCorners[edgeCorners2[0]]==edgeCorners1[1]
744 || hitCorners[edgeCorners2[1]]==edgeCorners1[0] || hitCorners[edgeCorners2[1]]==edgeCorners1[1])
747 if (edgeIntersection(corners1[edgeCorners1[0]],corners1[edgeCorners1[1]],
748 corners2[edgeCorners2[0]],corners2[edgeCorners2[1]],
749 directions1[edgeCorners1[0]],directions1[edgeCorners1[1]],
750 intersection, overlap) )
755 std::array<LocalCoords,2> corner;
756 corner[0] = ref1.template geometry<1>(i).global(intersection[0]);
757 corner[1] = ref2.template geometry<1>(j).global(intersection[1]);
758 polygonCorners.push_back(corner);
761 neighborIntersects2[j] =
true;
767 neighborIntersects1[i] =
true;
786 static bool edgeIntersection(
const WorldCoords& corner1,
const WorldCoords& corner2,
787 const WorldCoords& target1,
const WorldCoords& target2,
788 const WorldCoords& direction1,
const WorldCoords& direction2,
789 std::array<Dune::FieldVector<T,1>,2>& intersection,
const T overlap = 1e-1)
795 WorldCoords n21 = direction2 - direction1;
796 WorldCoords p21 = corner2 - corner1;
797 WorldCoords q21 = target2 - target1;
800 WorldCoords p1q1 = corner1 -target1;
803 T quadratic = direction1*q21n21;
806 T linear = (direction1*q21p21) + (p1q1*q21n21);
809 T constant = p1q1*q21p21;
812 std::vector<T> zeros;
814 if (std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
816 }
else if (std::fabs(quadratic)<1e-10) {
819 zeros.push_back(-constant/linear);
824 T p = linear/quadratic;
825 T q = constant/quadratic;
833 zeros.push_back(-0.5*p + std::sqrt(sqt));
834 zeros.push_back(-0.5*p -std::sqrt(sqt));
841 for (
size_t i=0;i<zeros.size();i++) {
853 WorldCoords dummy = p1q1; dummy.axpy(eta,direction1);
855 dummy = p21; dummy.axpy(eta,n21);
858 for (
int j=0;j<3; j++) {
861 if (isnan(r[0]) || isinf(r[0]))
864 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];
867 WorldCoords residual = p1q1;
868 residual.axpy(r[0],p21); residual.axpy(r[2],direction1);
869 residual.axpy(r[2]*r[0],n21); residual.axpy(-r[1],q21);
871 if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
874 residual.axpy(r[1],q21);
875 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];
877 residual.axpy(r[1],q21);
879 if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
883 if (r[0] >= -eps && r[1]>= -eps && (r[0]<=1+eps) && (r[1] <= 1+eps)) {
890 WorldCoords residual = p1q1;
891 residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
892 residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
894 if (residual.two_norm()<eps)
896 if (index >= 0 && x[0] >= -eps && x[1]>= -eps && (x[0]<=1+eps) && (x[1] <= 1+eps)) {
898 intersection[0] = x[0];
899 intersection[1] = x[1];
911 WorldCoords newtonCorrection(0);
913 for (
int i=0; i<30; i++) {
917 WorldCoords Fxk = target1 -corner1;
918 Fxk.axpy(-x[0], p21);
919 Fxk.axpy(-x[2],direction1);
920 Fxk.axpy(-x[2]*x[0],n21);
923 Dune::FieldMatrix<T,3,3> FPrimexk;
925 FPrimexk[0].axpy(x[2],n21);
926 FPrimexk[1] = target1-target2;
927 FPrimexk[2] = direction1;
928 FPrimexk[2].axpy(x[0],n21);
929 if (FPrimexk.determinant()<1e-10)
934 FPrimexk.mtv(Fxk, newtonCorrection);
936 x += newtonCorrection;
940 residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
941 residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
943 if (residual.two_norm()<=eps) {
945 if (x[0]>=-eps && x[0]<=(1+eps) && x[1]>=-eps && x[1]<=(1+eps)) {
949 intersection[0] = x[0];
950 intersection[1] = x[1];
978 template <
int dim,
int dimworld,
class T=
double>
980 const std::vector<Dune::FieldVector<T,dimworld> >& corners2,
981 const std::vector<Dune::FieldVector<T,dimworld> >& directions1,
982 const Dune::GeometryType& gt1,
const Dune::GeometryType& gt2,
983 std::vector<std::array<Dune::FieldVector<T,dim>,2> >& polygonCorners,
984 const std::vector<int>& hitCorners, std::bitset<(1<<dim)>& neighborIntersects1,
985 std::bitset<(1<<dim)>& neighborIntersects2,
const T overlap = 1e-1)
988 directions1, gt1, gt2, polygonCorners,
989 hitCorners, neighborIntersects1,
990 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:75
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:414
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< std::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:691
Helper class that provides static methods to compute the projection and inverse projection of a point...
Definition: projectionhelper.hh:39
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< std::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:710
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< std::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:669
Helper class that provides static methods for the computation of the intersection of surface element ...
Definition: projectionhelper.hh:648
This namespace contains helper functions for the projection of two triangular surface on each other...
Definition: projectionhelper.hh:19
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:177
static bool projection(const WorldCoords &corner, const WorldCoords &direction, const std::vector< WorldCoords > &targetCorners, LocalCoords &image, const T overlap=1e-1)
Definition: projectionhelper.hh:93
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:638
static bool projection(const WorldCoords &corner, const WorldCoords &direction, const std::vector< WorldCoords > &targetCorners, LocalCoords &image, const T overlap=1e-1)
Definition: projectionhelper.hh:500
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:542
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:618
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< std::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:979
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:56
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:23