dune-grid-glue  2.4-git
projectionhelper.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
4 #define DUNE_GRIDGLUE_COMMON_PROJECTIONHELPER_HH
5 
6 #include <bitset>
7 #include <memory>
8 
9 #include <dune/common/version.hh>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 
13 #include <dune/geometry/referenceelements.hh>
14 #include <dune/geometry/type.hh>
15 
19 namespace Projection
20 {
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)
25  {
26  if (dim!=3)
27  DUNE_THROW(Dune::NotImplemented, "crossProduct does not work for dimension " << dim);
28 
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];
33 
34  return c;
35  }
36 
38  template <int dim, int dimworld, class T=double>
40  {
41  typedef Dune::FieldVector<T,dimworld> WorldCoords;
42  typedef Dune::FieldVector<T,dim> LocalCoords;
43 
44  public:
45 
56  static bool inverseProjection(const std::vector<WorldCoords>& corners,
57  const std::vector<WorldCoords>& directions,
58  const WorldCoords& target, LocalCoords& preImage, const T overlap=1e-1) {
59 
60  DUNE_THROW(Dune::NotImplemented, "inverseProjection is not implemented for dimworld=="<<dimworld
61  <<" and dim=="<<dim);
62  return false;
63  }
64 
75  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
76  const std::vector<WorldCoords>& targetCorners,
77  LocalCoords& image, const T overlap=1e-1) {
78 
79  DUNE_THROW(Dune::NotImplemented, "projection is not implemented for dimworld=="<<dimworld
80  <<" and dim=="<<dim);
81  return false;
82 
83  }
84  };
85 
86  template <class T>
87  class ProjectionHelper<2,3,T>
88  {
89  typedef Dune::FieldVector<T,3> WorldCoords;
90  typedef Dune::FieldVector<T,2> LocalCoords;
91 
92  public:
93  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
94  const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
95  const T overlap=1e-1)
96  {
97 
98  if (targetCorners.size() != 3 && targetCorners.size() != 4)
99  DUNE_THROW(Dune::Exception, "projection is not implemented for polygone with "<<targetCorners.size()<<" vertices");
100 
101  // split up quadrilateral into triangles
102  if (targetCorners.size() == 4) {
103 
104  // lower triangle
105  std::vector<WorldCoords> triCorners(3);
106  for (int i=0; i<3; i++)
107  triCorners[i] = targetCorners[i];
108 
109  if (projection(corner,direction,triCorners,image,overlap))
110  return true;
111 
112  // upper triangle
113  triCorners[0] = targetCorners[3];
114 
115  if (projection(corner,direction,triCorners,image,overlap)) {
116  // Transform local coordinates of the upper triangle to quadrilateral ones
117  T a = image[0];
118  image[0] = 1-image[1];
119  image[1] = 1-a;
120  return true;
121  }
122 
123  return false;
124  }
125 
127  //
128  // Solve the linear system:
129  //
130  // (1-x-y)*t_0 + x*t_1 + y*t_2 + -z*direction = corner
131  //
132  // with (x,y) local coordinates in the target element
133  // z the 'direction' scaling
135 
136 
137  T eps = 1e-6;
138 
139  image = 0;
140 
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];
148  }
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];
153  }
154 
155  // Solve the linear system
156  WorldCoords rhs = corner - targetCorners[0];
157  WorldCoords x(0);
158  mat.solve(x,rhs);
159 
160  for (unsigned i = 0; i < 2; ++i)
161  x[i] /= scales[i];
162 
163  // only allow a certain overlap (we solved for '-z')
164  if (x[2]<-overlap)
165  return false;
166 
167  if (x[0]<-eps || x[1]<-eps || (x[0] + x[1]>1+eps) )
168  return false;
169 
170  image[0] = x[0];
171  image[1] = x[1];
172 
173  return true;
174  }
175 
176 
177  static bool inverseProjection(const std::vector<WorldCoords>& corners,
178  const std::vector<WorldCoords>& directions,
179  const WorldCoords& target, LocalCoords& preImage,
180  const T overlap = 1e-1)
181  {
182  const T eps = 1e-6;
183 
184  if (corners.size() != 3 && corners.size() != 4)
185  DUNE_THROW(Dune::NotImplemented, "inverseProjection is not implemented for elements with "<<corners.size()<<" corners!");
186 
187  // Split quadrilateral into triangles
188  if (corners.size() == 4) {
189 
190  // lower triangle
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];
195  }
196 
197  if (inverseProjection(triCorners,triDirections,target,preImage,overlap))
198  return true;
199 
200  // upper triangle
201  triCorners[0] = corners[3];
202  triDirections[0] = directions[3];
203 
204  if (inverseProjection(triCorners,triDirections,target,preImage,overlap)) {
205  // Transform local coordinates to quadrilateral ones
206  T a=preImage[0];
207  preImage[0] = 1-preImage[1];
208  preImage[1] = 1-a;
209 
210  return true;
211  }
212 
213  return false;
214  }
215 
216  // try to solve a cubic equation for the distance parameter, then compute the barycentric coordinates from it
217 
218  // the barycentric coordinates and the distance of the projected point
219  WorldCoords x(3);
220 
221  // cubic coefficient
222  WorldCoords n02 = directions[0] - directions[2];
223  WorldCoords n12 = directions[1] - directions[2];
224  WorldCoords n02n12 = crossProduct(n02,n12);
225 
226  T cubic = (directions[2]*n02n12);
227 
228  // quadratic coefficient
229 
230  WorldCoords p02 = corners[0] - corners[2];
231  WorldCoords p12 = corners[1] - corners[2];
232  WorldCoords p2q = corners[2] -target;
233  WorldCoords p02n12 = crossProduct(p02,n12);
234  WorldCoords n02p12 = crossProduct(n02,p12);
235 
236  T quadratic = (directions[2]*p02n12)+ (directions[2]*n02p12)+ (p2q*n02n12);
237 
238  // constant coefficient
239  WorldCoords p02p12 = crossProduct(p02,p12);
240  T constant = p2q*p02p12;
241 
242  // linear coefficient
243  T linear = (directions[2]*p02p12) + (p2q*n02p12) + (p2q*p02n12);
244 
245  // save all zeros we find
246  std::vector<T> zeros;
247 
248  if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
249  return false;
250  } else if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10) {
251 
252  // problem is linear
253  zeros.push_back(-constant/linear);
254 
255  } else if(std::fabs(cubic)<1e-10) {
256 
257  // problem is quadratic
258  T p = linear/quadratic;
259  T q = constant/quadratic;
260 
261  T sqt = 0.25*p*p -q;
262 
263  // no real solution
264  if (sqt<-1e-10)
265  return false;
266 
267  zeros.push_back(-0.5*p + std::sqrt(sqt));
268  zeros.push_back(-0.5*p -std::sqrt(sqt));
269 
270  } else {
271 
272  // problem is cubic
273  quadratic /= cubic;
274  linear /= cubic;
275  constant /= cubic;
276 
277  // Transform to reduced form z^3 + p*z + q = 0 where x = z-quadratic/3
278  T p= linear - quadratic*quadratic/3;
279  T q=quadratic*(2*quadratic*quadratic/27 - linear/3) + constant;
280 
281  // use Cardano's method to solve the problem
282  T D = 0.25*q*q + std::pow(p,3)/27;
283 
284  if (D>1e-10) {
285  // one real zero
286 
287  // be careful when computing the cubic roots
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);
290 
291  nu = -q/2-std::sqrt(D);
292  zer += std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
293 
294  zeros.push_back(zer-quadratic/3);
295 
296  } else if (D<-1e-10) {
297 
298  // three real zeros, using trigonometric functions to compute them
299  T a = std::sqrt(-4*p/3);
300  T b = std::acos(-0.5*q*std::sqrt(-27/(std::pow(p,3))));
301 
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);
304 
305 
306  } else {
307  // one single and one double zero
308 
309  if (std::fabs(q)<1e-10) {
310  zeros.push_back(-quadratic/3);
311 
312  if (p<-1e-10)
313  zeros.push_back(std::sqrt(-p)-quadratic/3);
314 
315  } else if (std::fabs(p)<1e-10) { // is this case correct?
316 
317  T nu = std::pow(std::fabs(q),1.0/3.0) * ((q<-eps) ? -1 : 1);
318  zeros.push_back(nu-quadratic/3);
319 
320  } else {
321  zeros.push_back(3*q/p - quadratic/3);
322  zeros.push_back(-1.5*q/p - quadratic/3);
323  }
324  }
325  }
326 
327  int index = -1;
328  WorldCoords r;
329 
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]);
337  return residual;
338  };
339 
340  for (size_t i=0;i<zeros.size();i++) {
341 
342  T nu=zeros[i];
343  // only look in the direction of the outer normals
344  if (nu<-overlap) // allowed overlap
345  continue;
346 
347  if (index != -1)
348  if (nu > zeros[index]) // is this one really closer ?
349  continue;
350 
351  r[2] = nu;
352 
353  // the computation of the other components might lead to nan or inf
354  // if this happens use a different equation to compute them
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);
358  WorldCoords c = crossProduct(e,g);
359  WorldCoords d = crossProduct(g,f);
360 
361  // computation of the other components is unstable
362  for (int j=0;j<3; j++) {
363 
364  r[1] = c[j]/d[j];
365 
366  if (isnan(r[1]) || isinf(r[1]))
367  continue;
368 
369  r[0] = -(e[(j+1)%3]+r[1]*f[(j+1)%3])/g[(j+1)%3];
370 
371  // Check residual
372  if (!(isnan(r[0]) || isinf(r[0])) && (residual(r).two_norm()<1e-5))
373  break;
374 
375  r[0] = -(e[(j+2)%3]+r[1]*f[(j+2)%3])/g[(j+2)%3];
376 
377  // Check residual
378  if (!(isnan(r[0]) || isinf(r[0])) && (residual(r).two_norm()<1e-5))
379  break;
380 
381  }
382 
383  if (r[0] > -eps && r[1]> -eps && (r[0]+r[1] < 1+eps)) {
384  index = i;
385  x = r;
386  }
387  }
388 
389  // Check if we found a feasible zero
390  if (residual(x).two_norm()<1e-6) {
391  if (index >= 0) {
392  preImage[0] = x[1];
393  preImage[1] = 1-x[0]-x[1];
394 
395  return true;
396  }
397  return false;
398  }
399 
400 
401  //std::cout<<"Direct solution failed, use Newton method\n";
402  // In some cases the direct solution of the cubic equation is unstable, in this case we use
403  // Newton's method to compute at least one solution
404 
405  // Some problems have two solutions and the Newton converges to the wrong one
406  return inexactInverseProjection(corners,directions, target, preImage,overlap);
407  }
408 
414  static bool inexactInverseProjection(const std::vector<WorldCoords>& corners,
415  const std::vector<WorldCoords>& directions,
416  const WorldCoords& target, LocalCoords& preImage, const T overlap=1e-1)
417  {
418  assert(corners.size() == 3);
419  assert(directions.size() == 3);
420 
421  // feasible initial Newton iterate
422  const int nCorners = 3;
423  Dune::FieldVector<T,nCorners> x(1.0/((T) nCorners));
424  {
425  WorldCoords d(0);
426  for (std::size_t i = 0; i < corners.size(); ++i)
427  d += corners[i];
428  d *= 1./3;
429  d -= target;
430  x[2] = d.two_norm();
431  }
432 
433  for (int i=0; i<30; i++) {
434 
435  // compute Newton correction
436  WorldCoords Fxk = target -corners[2];
437 
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]);
441  }
442  Fxk.axpy(-x[2],directions[2]);
443 
444  Dune::FieldMatrix<T,nCorners,nCorners> FPrimexk(0);
445 
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]);
453 
454  try {
455  FPrimexk.invert();
456  }
457  catch (const Dune::FMatrixError&) {
458  return false;
459  }
460 
461  WorldCoords newtonCorrection; // = (-1) * FPrimexk.inverse() * Fxk;
462 
463  FPrimexk.mtv(Fxk, newtonCorrection);
464 
465  x += newtonCorrection;
466  }
467 
468 
469  if (x[0]>-1e-6 && x[1]>-1e-6 && (x[0]+x[1] <1+1e-6)) {
470 
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]);
477 
478  // Newton did not converge
479  if (residual.two_norm()>1e-6 || x[2]<-overlap)
480  return false;
481 
482  preImage[0] = x[1];
483  preImage[1] = 1-x[0]-x[1];
484 
485  return true;
486  }
487 
488  return false;
489  }
490  };
491 
492  template <class T>
493  class ProjectionHelper<1,2,T>
494  {
495 
496  typedef Dune::FieldVector<T,2> WorldCoords;
497  typedef Dune::FieldVector<T,1> LocalCoords;
498 
499  public:
500  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
501  const std::vector<WorldCoords>& targetCorners, LocalCoords& image, const T overlap=1e-1)
502  {
503  T eps = 1e-8;
504  // we solve the equation basePoint + x_0 * normal = a + x_1 * (b-a)
505  image = 0;
506 
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];
512 
513  WorldCoords rhs = targetCorners[0] - corner;
514  WorldCoords x(0);
515 
516  // Solve the system. If it is singular the normal and the segment
517  // are parallel and there is no intersection
518 
519  T detinv = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
520  if (std::abs(detinv)<1e-14)
521  return false;
522 
523  detinv = 1/detinv;
524 
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]);
527 
528  // x[0] is the distance, x[1] is the intersection point
529  // in local coordinates on the segment
530  if (x[1]<-eps || x[1] > 1+eps)
531  return false;
532 
533  // allow some overlap
534  if (x[0]<-overlap)
535  return false;
536 
537  image[0] = x[1];
538 
539  return true;
540  }
541 
542  static bool inverseProjection(const std::vector<WorldCoords>& corners,
543  const std::vector<WorldCoords>& directions,
544  const WorldCoords& target, LocalCoords& preImage,
545  const T overlap=1e-1)
546  {
547  T eps = 1e-8;
548 
549  preImage = 0;
550  WorldCoords p10 = corners[1]-corners[0];
551  WorldCoords v10 = directions[1]-directions[0];
552  WorldCoords qp0 = target-corners[0];
553 
554  // Compute coefficients of the polynomial by eliminating the distance variable
555  T a = p10[1]*v10[0] - p10[0]*v10[1];
556 
557  T b = -qp0[1]*v10[0] +p10[1]*directions[0][0]
558  + qp0[0]*v10[1] -p10[0]*directions[0][1];
559 
560  T c = -qp0[1]*directions[0][0] +qp0[0]*directions[0][1];
561 
562  // Is the quadratic formula degenerated to a linear one?
563  if (std::abs(a) < 1e-10) {
564  T x = -c/b;
565 
566  if (x >= -eps && x <= 1+eps) {
567 
568  //check if projection is along the positive directions
569  T dist = (qp0[0]-x*p10[0])/(directions[0][0]+x*v10[0]);
570  if (dist <-overlap)
571  return false;
572 
573  preImage[0] = x;
574  return true;
575  }
576  return false;
577  }
578 
579  // The abc-formula
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);
582 
583  if (mu_0 >= -eps && mu_0 <= 1+eps) {
584 
585  T dist = (qp0[0]-mu_0*p10[0])/(directions[0][0]+mu_0*v10[0]);
586  if (dist <-overlap)
587  return false;
588 
589  preImage[0] = mu_0;
590  return true;
591 
592  } else if (mu_1 >= -eps && mu_1 <= 1+eps) {
593 
594  T dist = (qp0[0]-mu_1*p10[0])/(directions[0][0]+mu_1*v10[0]);
595  if (dist <-overlap)
596  return false;
597 
598  preImage[0] = mu_1;
599  return true;
600  }
601  return false;
602 
603  }
604  };
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)
621  {
622  return ProjectionHelper<dim,dimworld,T>::projection(corner,direction,
623  targetCorners,image,overlap);
624  }
625 
637  template <int dim, int dimworld, class T>
638  static bool inverseProjection(const std::vector<Dune::FieldVector<T,dimworld> >& corners,
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)
642  {
643  return ProjectionHelper<dim,dimworld,T>::inverseProjection(corners,directions,
644  target, preImage, overlap);
645  }
647  template <int dim, int dimworld, class T>
649  {
650 
651  typedef Dune::FieldVector<T,dimworld> WorldCoords;
652  typedef Dune::FieldVector<T,dim> LocalCoords;
653 
654  public:
669  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
676  {
677  DUNE_THROW(Dune::NotImplemented, "addEdgeIntersections is not implemented for dimworld=="<<dimworld
678  <<" and dim=="<<dim);
679  }
680  };
681 
682  template <typename T>
683  class EdgeIntersectionHelper<1, 2, T>
684  {
685 
686  typedef Dune::FieldVector<T,2> WorldCoords;
687  typedef Dune::FieldVector<T,1> LocalCoords;
688 
689  public:
691  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
698  {}
699 
700  };
701 
702  template <typename T>
703  class EdgeIntersectionHelper<2, 3, T>
704  {
705 
706  typedef Dune::FieldVector<T,3> WorldCoords;
707  typedef Dune::FieldVector<T,2> LocalCoords;
708 
709  public:
710  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
717  {
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);
721 #else
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);
724 #endif
725 
726  //loop over edges
727  std::array<Dune::FieldVector<T,1>,2> intersection;
728  for (int i=0; i<ref1.size(1); i++) {
729 
730  std::array<int, 2> edgeCorners1;
731  edgeCorners1[0] = ref1.subEntity(i,1,0,2);
732  edgeCorners1[1] = ref1.subEntity(i,1,1,2);
733 
734  int nIntersects(0);
735  for (int j=0; j<ref2.size(1); j++) {
736 
737  std::array<int, 2> edgeCorners2;
738  edgeCorners2[0] = ref2.subEntity(j,1,0,2);
739  edgeCorners2[1] = ref2.subEntity(j,1,1,2);
740 
741  // if any edge endpoints hit each other we don't have to compute the intersections
742  // either there is none or we already have computed it before
743  if (hitCorners[edgeCorners2[0]]==edgeCorners1[0] || hitCorners[edgeCorners2[0]]==edgeCorners1[1]
744  || hitCorners[edgeCorners2[1]]==edgeCorners1[0] || hitCorners[edgeCorners2[1]]==edgeCorners1[1])
745  continue;
746 
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) )
751  {
752  nIntersects++;
753 
754  // compute the local coordinates
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);
759 
760  // if the grid2 edge intersected then the neighbor will also intersect
761  neighborIntersects2[j] = true;
762  }
763  }
764 
765  // if grid1 edge intersected an other edge then the neighbor will also intersect
766  if (nIntersects>0)
767  neighborIntersects1[i] = true;
768  }
769  }
770 
771  private:
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)
790  {
791 
792  T eps = 1e-6;
793  // solve a quadratic scalar equation for the distance parameter eta, then compute the barycentric coordinates from it
794 
795  WorldCoords n21 = direction2 - direction1;
796  WorldCoords p21 = corner2 - corner1;
797  WorldCoords q21 = target2 - target1;
798  WorldCoords q21n21 = crossProduct(q21,n21);
799  WorldCoords q21p21 = crossProduct(q21,p21);
800  WorldCoords p1q1 = corner1 -target1;
801 
802  // quadratic coefficient
803  T quadratic = direction1*q21n21;
804 
805  // linear coefficient
806  T linear = (direction1*q21p21) + (p1q1*q21n21);
807 
808  // constant coefficient
809  T constant = p1q1*q21p21;
810 
811  // save all zeros we find
812  std::vector<T> zeros;
813 
814  if (std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
815  return false;
816  } else if (std::fabs(quadratic)<1e-10) {
817 
818  // problem is linear
819  zeros.push_back(-constant/linear);
820 
821  } else {
822 
823  // problem is quadratic
824  T p = linear/quadratic;
825  T q = constant/quadratic;
826 
827  T sqt = 0.25*p*p -q;
828 
829  // no real solution
830  if (sqt<-1e-10)
831  return false;
832 
833  zeros.push_back(-0.5*p + std::sqrt(sqt));
834  zeros.push_back(-0.5*p -std::sqrt(sqt));
835 
836  }
837 
838  int index = -1;
839  WorldCoords r,x(-1);
840 
841  for (size_t i=0;i<zeros.size();i++) {
842 
843  T eta=zeros[i];
844 
845  // only look in the direction of the outer normals
846  if (eta<-overlap)
847  continue;
848 
849  r[2] = eta;
850 
851  // the computation of the other components might lead to nan or inf
852  // if this happens use a different equation to compute them
853  WorldCoords dummy = p1q1; dummy.axpy(eta,direction1);
854  WorldCoords c =crossProduct(dummy,q21);
855  dummy = p21; dummy.axpy(eta,n21);
856  WorldCoords d =crossProduct(q21,dummy);
857 
858  for (int j=0;j<3; j++) {
859 
860  r[0] = c[j]/d[j];
861  if (isnan(r[0]) || isinf(r[0]))
862  continue;
863 
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];
865 
866  // computation of the other components can be instable
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);
870 
871  if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
872  break;
873 
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];
876 
877  residual.axpy(r[1],q21);
878  // computation of the other components can be instable
879  if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
880  break;
881 
882  }
883  if (r[0] >= -eps && r[1]>= -eps && (r[0]<=1+eps) && (r[1] <= 1+eps)) {
884  index = i;
885  x=r;
886  }
887  }
888 
889  // TODO CLEAN THIS UP
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);
893 
894  if (residual.two_norm()<eps)
895  {
896  if (index >= 0 && x[0] >= -eps && x[1]>= -eps && (x[0]<=1+eps) && (x[1] <= 1+eps)) {
897 
898  intersection[0] = x[0];
899  intersection[1] = x[1];
900  return true;
901  }
902  return false;
903  }
904 
905  // if the direct compuation failed, use a Newton method to compute at least one zero
906 
907  // Fix some initial value
908  // sometimes it only works when the initial value is an intersection...
909  x[0] = x[1] = 0.5;
910  x[2] = 0.5;
911  WorldCoords newtonCorrection(0);
912 
913  for (int i=0; i<30; i++) {
914 
915  // compute Newton correction
916 
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);
921  Fxk.axpy(x[1],q21);
922 
923  Dune::FieldMatrix<T,3,3> FPrimexk;
924  FPrimexk[0] = p21;
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)
930  return false;
931 
932  FPrimexk.invert();
933 
934  FPrimexk.mtv(Fxk, newtonCorrection);
935 
936  x += newtonCorrection;
937  }
938 
939  residual = p1q1;
940  residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
941  residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
942 
943  if (residual.two_norm()<=eps) {
944 
945  if (x[0]>=-eps && x[0]<=(1+eps) && x[1]>=-eps && x[1]<=(1+eps)) {
946  if (x[2]<-overlap)
947  return false;
948 
949  intersection[0] = x[0];
950  intersection[1] = x[1];
951 
952  return true;
953  }
954  return false;
955 
956  }
957 
958  //std::cout<<"Newton did not converge either!\n";
959 
960  return false;
961  }
962  };
963 
978  template <int dim, int dimworld, class T=double>
979  static void addEdgeIntersections(const std::vector<Dune::FieldVector<T,dimworld> >& corners1,
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)
986  {
988  directions1, gt1, gt2, polygonCorners,
989  hitCorners, neighborIntersects1,
990  neighborIntersects2, overlap);
991  }
992 
993 } // end namespace
994 #endif
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 &gt1, const Dune::GeometryType &gt2, 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 &gt1, const Dune::GeometryType &gt2, 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 &gt1, Dune::GeometryType &gt2, 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 &gt1, const Dune::GeometryType &gt2, 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