dune-grid-glue  2.3.0
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 
8 #include <dune/common/array.hh>
9 #include <dune/common/fvector.hh>
10 #include <dune/common/fmatrix.hh>
11 
15 namespace Projection
16 {
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)
21  {
22  if (dim!=3)
23  DUNE_THROW(Dune::NotImplemented, "crossProduct does not work for dimension " << dim);
24 
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];
29 
30  return c;
31  }
32 
34  template <int dim, int dimworld, class T=double>
36  {
37  typedef Dune::FieldVector<T,dimworld> WorldCoords;
38  typedef Dune::FieldVector<T,dim> LocalCoords;
39 
40  public:
41 
52  static bool inverseProjection(const std::vector<WorldCoords>& corners,
53  const std::vector<WorldCoords>& directions,
54  const WorldCoords& target, LocalCoords& preImage, const T overlap=1e-1) {
55 
56  DUNE_THROW(Dune::NotImplemented, "inverseProjection is not implemented for dimworld=="<<dimworld
57  <<" and dim=="<<dim);
58  return false;
59  }
60 
71  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
72  const std::vector<WorldCoords>& targetCorners,
73  LocalCoords& image, const T overlap=1e-1) {
74 
75  DUNE_THROW(Dune::NotImplemented, "projection is not implemented for dimworld=="<<dimworld
76  <<" and dim=="<<dim);
77  return false;
78 
79  }
80  };
81 
82  template <class T>
83  class ProjectionHelper<2,3,T>
84  {
85  typedef Dune::FieldVector<T,3> WorldCoords;
86  typedef Dune::FieldVector<T,2> LocalCoords;
87 
88  public:
89  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
90  const std::vector<WorldCoords>& targetCorners, LocalCoords& image,
91  const T overlap=1e-1)
92  {
93 
94  if (targetCorners.size() != 3 && targetCorners.size() != 4)
95  DUNE_THROW(Dune::Exception, "projection is not implemented for polygone with "<<targetCorners.size()<<" vertices");
96 
97  // split up quadrilateral into triangles
98  if (targetCorners.size() == 4) {
99 
100  // lower triangle
101  std::vector<WorldCoords> triCorners(3);
102  for (int i=0; i<3; i++)
103  triCorners[i] = targetCorners[i];
104 
105  if (projection(corner,direction,triCorners,image,overlap))
106  return true;
107 
108  // upper triangle
109  triCorners[0] = targetCorners[3];
110 
111  if (projection(corner,direction,triCorners,image,overlap)) {
112  // Transform local coordinates of the upper triangle to quadrilateral ones
113  T a = image[0];
114  image[0] = 1-image[1];
115  image[1] = 1-a;
116  return true;
117  }
118  }
119 
120 
122  //
123  // Solve the linear system:
124  //
125  // (1-x-y)*t_0 + x*t_1 + y*t_2 + -z*direction = corner
126  //
127  // with (x,y) local coordinates in the target element
128  // z the 'direction' scaling
130 
131 
132  T eps = 1e-6;
133 
134  image = 0;
135 
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];
141  }
142  if (std::fabs(mat.determinant())<eps)
143  return false;
144 
145  // Solve the linear system
146  WorldCoords rhs = corner - targetCorners[0];
147  WorldCoords x(0);
148  mat.solve(x,rhs);
149 
150  // only allow a certain overlap (we solved for '-z')
151  if (x[2]<-overlap)
152  return false;
153 
154  if (x[0]<-eps || x[1]<-eps || (x[0] + x[1]>1+eps) )
155  return false;
156 
157  image[0] = x[0];
158  image[1] = x[1];
159 
160  return true;
161  }
162 
163 
164  static bool inverseProjection(const std::vector<WorldCoords>& corners,
165  const std::vector<WorldCoords>& directions,
166  const WorldCoords& target, LocalCoords& preImage,
167  const T overlap = 1e-1)
168  {
169  const T eps = 1e-6;
170 
171  if (corners.size() != 3 && corners.size() != 4)
172  DUNE_THROW(Dune::NotImplemented, "inverseProjection is not implemented for elements with "<<corners.size()<<" corners!");
173 
174  // Split quadrilateral into triangles
175  if (corners.size() == 4) {
176 
177  // lower triangle
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];
182  }
183 
184  if (inverseProjection(triCorners,triDirections,target,preImage,overlap))
185  return true;
186 
187  // upper triangle
188  triCorners[0] = corners[3];
189  triDirections[0] = directions[3];
190 
191  if (inverseProjection(triCorners,triDirections,target,preImage,overlap)) {
192  // Transform local coordinates to quadrilateral ones
193  T a=preImage[0];
194  preImage[0] = 1-preImage[1];
195  preImage[1] = 1-a;
196 
197  return true;
198  }
199  }
200 
201  // try to solve a cubic equation for the distance parameter, then compute the barycentric coordinates from it
202 
203  // the barycentric coordinates and the distance of the projected point
204  WorldCoords x(3);
205 
206  // cubic coefficient
207  WorldCoords n02 = directions[0] - directions[2];
208  WorldCoords n12 = directions[1] - directions[2];
209  WorldCoords n02n12 = crossProduct(n02,n12);
210 
211  T cubic = (directions[2]*n02n12);
212 
213  // quadratic coefficient
214 
215  WorldCoords p02 = corners[0] - corners[2];
216  WorldCoords p12 = corners[1] - corners[2];
217  WorldCoords p2q = corners[2] -target;
218  WorldCoords p02n12 = crossProduct(p02,n12);
219  WorldCoords n02p12 = crossProduct(n02,p12);
220 
221  T quadratic = (directions[2]*p02n12)+ (directions[2]*n02p12)+ (p2q*n02n12);
222 
223  // constant coefficient
224  WorldCoords p02p12 = crossProduct(p02,p12);
225  T constant = p2q*p02p12;
226 
227  // linear coefficient
228  T linear = (directions[2]*p02p12) + (p2q*n02p12) + (p2q*p02n12);
229 
230  // save all zeros we find
231  std::vector<T> zeros;
232 
233  if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
234  return false;
235  } else if (std::fabs(cubic) <1e-10 && std::fabs(quadratic)<1e-10) {
236 
237  // problem is linear
238  zeros.push_back(-constant/linear);
239 
240  } else if(std::fabs(cubic)<1e-10) {
241 
242  // problem is quadratic
243  T p = linear/quadratic;
244  T q = constant/quadratic;
245 
246  T sqt = 0.25*p*p -q;
247 
248  // no real solution
249  if (sqt<-1e-10)
250  return false;
251 
252  zeros.push_back(-0.5*p + std::sqrt(sqt));
253  zeros.push_back(-0.5*p -std::sqrt(sqt));
254 
255  } else {
256 
257  // problem is cubic
258  quadratic /= cubic;
259  linear /= cubic;
260  constant /= cubic;
261 
262  // Transform to reduced form z^3 + p*z + q = 0 where x = z-quadratic/3
263  T p= linear - quadratic*quadratic/3;
264  T q=quadratic*(2*quadratic*quadratic/27 - linear/3) + constant;
265 
266  // use Cardano's method to solve the problem
267  T D = 0.25*q*q + std::pow(p,3)/27;
268 
269  if (D>1e-10) {
270  // one real zero
271 
272  // be careful when computing the cubic roots
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);
275 
276  nu = -q/2-std::sqrt(D);
277  zer += std::pow(std::fabs(nu),1.0/3.0) * ((nu<-1e-10) ? -1 : 1);
278 
279  zeros.push_back(zer-quadratic/3);
280 
281  } else if (D<-1e-10) {
282 
283  // three real zeros, using trigonometric functions to compute them
284  T a = std::sqrt(-4*p/3);
285  T b = std::acos(-0.5*q*std::sqrt(-27/(std::pow(p,3))));
286 
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);
289 
290 
291  } else {
292  // one single and one double zero
293 
294  if (std::fabs(q)<1e-10) {
295  zeros.push_back(-quadratic/3);
296 
297  if (p<-1e-10)
298  zeros.push_back(std::sqrt(-p)-quadratic/3);
299 
300  } else if (std::fabs(p)<1e-10) { // is this case correct?
301 
302  T nu = std::pow(std::fabs(q),1.0/3.0) * ((q<-eps) ? -1 : 1);
303  zeros.push_back(nu-quadratic/3);
304 
305  } else {
306  zeros.push_back(3*q/p - quadratic/3);
307  zeros.push_back(-1.5*q/p - quadratic/3);
308  }
309  }
310  }
311 
312  int index = -1;
313  WorldCoords r;
314 
315  for (size_t i=0;i<zeros.size();i++) {
316 
317  T nu=zeros[i];
318  // only look in the direction of the outer normals
319  if (nu<-overlap) // allowed overlap
320  continue;
321 
322  if (index != -1)
323  if (nu > zeros[index]) // is this one really closer ?
324  continue;
325 
326  r[2] = nu;
327 
328  // the computation of the other components might lead to nan or inf
329  // if this happens use a different equation to compute them
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);
333  WorldCoords c = crossProduct(e,g);
334  WorldCoords d = crossProduct(g,f);
335 
336  // computation of the other components is unstable
337  for (int j=0;j<3; j++) {
338 
339  r[1] = c[j]/d[j];
340 
341  if (isnan(r[1]) || isinf(r[1]))
342  continue;
343 
344  r[0] = -(e[(j+1)%3]+r[1]*f[(j+1)%3])/g[(j+1)%3];
345 
346  // Compute the residual
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]);
353 
354 
355  if (!(isnan(r[0]) || isinf(r[0])) && (residual.two_norm()<1e-5))
356  break;
357 
358  r[0] = -(e[(j+2)%3]+r[1]*f[(j+2)%3])/g[(j+2)%3];
359 
360  // Compute the residual
361  residual = p2q;
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]);
367 
368  if (!(isnan(r[0]) || isinf(r[0])) && (residual.two_norm()<1e-5))
369  break;
370 
371  }
372 
373  if (r[0] > -eps && r[1]> -eps && (r[0]+r[1] < 1+eps)) {
374  index = i;
375  x = r;
376  }
377  }
378 
379  // Check if we found a feasible zero
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]);
386 
387  if (residual.two_norm()<1e-6) {
388  if (index >= 0) {
389  preImage[0] = x[1];
390  preImage[1] = 1-x[0]-x[1];
391 
392  return true;
393  }
394  return false;
395  }
396 
397 
398  //std::cout<<"Direct solution failed, use Newton method\n";
399  // In some cases the direct solution of the cubic equation is unstable, in this case we use
400  // Newton's method to compute at least one solution
401 
402  // Some problems have two solutions and the Newton converges to the wrong one
403  return inexactInverseProjection(corners,directions, target, preImage,overlap);
404  }
405 
411  static bool inexactInverseProjection(const std::vector<WorldCoords>& corners,
412  const std::vector<WorldCoords>& directions,
413  const WorldCoords& target, LocalCoords& preImage, const T overlap=1e-1)
414  {
415 
416  // feasible initial Newton iterate
417  const int nCorners = 3;
418  Dune::FieldVector<T,nCorners> x(1.0/((T) nCorners));
419 
420  for (int i=0; i<30; i++) {
421 
422  // compute Newton correction
423  WorldCoords Fxk = target -corners[2];
424 
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]);
428  }
429  Fxk.axpy(-x[2],directions[2]);
430 
431  Dune::FieldMatrix<T,nCorners,nCorners> FPrimexk(0);
432 
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]);
440 
441  FPrimexk.invert();
442 
443  WorldCoords newtonCorrection; // = (-1) * FPrimexk.inverse() * Fxk;
444 
445  FPrimexk.mtv(Fxk, newtonCorrection);
446 
447  x += newtonCorrection;
448  }
449 
450 
451  if (x[0]>-1e-6 && x[1]>-1e-6 && (x[0]+x[1] <1+1e-6)) {
452 
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]);
459 
460  // Newton did not converge
461  if (residual.two_norm()>1e-6 || x[2]<-overlap)
462  return false;
463 
464  preImage[0] = x[1];
465  preImage[1] = 1-x[0]-x[1];
466 
467  return true;
468  }
469 
470  return false;
471  }
472  };
473 
474  template <class T>
475  class ProjectionHelper<1,2,T>
476  {
477 
478  typedef Dune::FieldVector<T,2> WorldCoords;
479  typedef Dune::FieldVector<T,1> LocalCoords;
480 
481  public:
482  static bool projection(const WorldCoords& corner, const WorldCoords& direction,
483  const std::vector<WorldCoords>& targetCorners, LocalCoords& image, const T overlap=1e-1)
484  {
485  T eps = 1e-8;
486  // we solve the equation basePoint + x_0 * normal = a + x_1 * (b-a)
487  image = 0;
488 
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];
494 
495  WorldCoords rhs = targetCorners[0] - corner;
496  WorldCoords x(0);
497 
498  // Solve the system. If it is singular the normal and the segment
499  // are parallel and there is no intersection
500 
501  T detinv = mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
502  if (std::abs(detinv)<1e-14)
503  return false;
504 
505  detinv = 1/detinv;
506 
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]);
509 
510  // x[0] is the distance, x[1] is the intersection point
511  // in local coordinates on the segment
512  if (x[1]<-eps || x[1] > 1+eps)
513  return false;
514 
515  // allow some overlap
516  if (x[0]<-overlap)
517  return false;
518 
519  image[0] = x[1];
520 
521  return true;
522  }
523 
524  static bool inverseProjection(const std::vector<WorldCoords>& corners,
525  const std::vector<WorldCoords>& directions,
526  const WorldCoords& target, LocalCoords& preImage,
527  const T overlap=1e-1)
528  {
529  T eps = 1e-8;
530 
531  preImage = 0;
532  WorldCoords p10 = corners[1]-corners[0];
533  WorldCoords v10 = directions[1]-directions[0];
534  WorldCoords qp0 = target-corners[0];
535 
536  // Compute coefficients of the polynomial by eliminating the distance variable
537  T a = p10[1]*v10[0] - p10[0]*v10[1];
538 
539  T b = -qp0[1]*v10[0] +p10[1]*directions[0][0]
540  + qp0[0]*v10[1] -p10[0]*directions[0][1];
541 
542  T c = -qp0[1]*directions[0][0] +qp0[0]*directions[0][1];
543 
544  // Is the quadratic formula degenerated to a linear one?
545  if (std::abs(a) < 1e-10) {
546  T x = -c/b;
547 
548  if (x >= -eps && x <= 1+eps) {
549 
550  //check if projection is along the positive directions
551  T dist = (qp0[0]-x*p10[0])/(directions[0][0]+x*v10[0]);
552  if (dist <-overlap)
553  return false;
554 
555  preImage[0] = x;
556  return true;
557  }
558  return false;
559  }
560 
561  // The abc-formula
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);
564 
565  if (mu_0 >= -eps && mu_0 <= 1+eps) {
566 
567  T dist = (qp0[0]-mu_0*p10[0])/(directions[0][0]+mu_0*v10[0]);
568  if (dist <-overlap)
569  return false;
570 
571  preImage[0] = mu_0;
572  return true;
573 
574  } else if (mu_1 >= -eps && mu_1 <= 1+eps) {
575 
576  T dist = (qp0[0]-mu_1*p10[0])/(directions[0][0]+mu_1*v10[0]);
577  if (dist <-overlap)
578  return false;
579 
580  preImage[0] = mu_1;
581  return true;
582  }
583  return false;
584 
585  }
586  };
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)
603  {
604  return ProjectionHelper<dim,dimworld,T>::projection(corner,direction,
605  targetCorners,image,overlap);
606  }
607 
619  template <int dim, int dimworld, class T>
620  static bool inverseProjection(const std::vector<Dune::FieldVector<T,dimworld> >& corners,
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)
624  {
625  return ProjectionHelper<dim,dimworld,T>::inverseProjection(corners,directions,
626  target, preImage, overlap);
627  }
629  template <int dim, int dimworld, class T>
631  {
632 
633  typedef Dune::FieldVector<T,dimworld> WorldCoords;
634  typedef Dune::FieldVector<T,dim> LocalCoords;
635 
636  public:
651  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
658  {
659  DUNE_THROW(Dune::NotImplemented, "addEdgeIntersections is not implemented for dimworld=="<<dimworld
660  <<" and dim=="<<dim);
661  }
662  };
663 
664  template <typename T>
665  class EdgeIntersectionHelper<1, 2, T>
666  {
667 
668  typedef Dune::FieldVector<T,2> WorldCoords;
669  typedef Dune::FieldVector<T,1> LocalCoords;
670 
671  public:
673  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
680  {}
681 
682  };
683 
684  template <typename T>
685  class EdgeIntersectionHelper<2, 3, T>
686  {
687 
688  typedef Dune::FieldVector<T,3> WorldCoords;
689  typedef Dune::FieldVector<T,2> LocalCoords;
690 
691  public:
692  static void addEdgeIntersections(const std::vector<WorldCoords>& corners1,
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)
699  {
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);
703 #else
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);
706 #endif
707 
708  //loop over edges
709  Dune::array<Dune::FieldVector<T,1>,2> intersection;
710  for (int i=0; i<ref1.size(1); i++) {
711 
712  std::vector<int> edgeCorners1(2);
713  edgeCorners1[0] = ref1.subEntity(i,1,0,2);
714  edgeCorners1[1] = ref1.subEntity(i,1,1,2);
715 
716  int nIntersects(0);
717  for (int j=0; j<ref2.size(1); j++) {
718 
719  std::vector<int> edgeCorners2(2);
720  edgeCorners2[0] = ref2.subEntity(j,1,0,2);
721  edgeCorners2[1] = ref2.subEntity(j,1,1,2);
722 
723  // if any edge endpoints hit each other we don't have to compute the intersections
724  // either there is none or we already have computed it before
725  if (hitCorners[edgeCorners2[0]]==edgeCorners1[0] || hitCorners[edgeCorners2[0]]==edgeCorners1[1]
726  || hitCorners[edgeCorners2[1]]==edgeCorners1[0] || hitCorners[edgeCorners2[1]]==edgeCorners1[1])
727  continue;
728 
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) )
733  {
734  nIntersects++;
735 
736  // compute the local coordinates
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);
741 
742  // if the grid2 edge intersected then the neighbor will also intersect
743  neighborIntersects2[j] = true;
744  }
745  }
746 
747  // if grid1 edge intersected an other edge then the neighbor will also intersect
748  if (nIntersects>0)
749  neighborIntersects1[i] = true;
750  }
751  }
752 
753  private:
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)
772  {
773 
774  T eps = 1e-6;
775  // solve a quadratic scalar equation for the distance parameter eta, then compute the barycentric coordinates from it
776 
777  WorldCoords n21 = direction2 - direction1;
778  WorldCoords p21 = corner2 - corner1;
779  WorldCoords q21 = target2 - target1;
780  WorldCoords q21n21 = crossProduct(q21,n21);
781  WorldCoords q21p21 = crossProduct(q21,p21);
782  WorldCoords p1q1 = corner1 -target1;
783 
784  // quadratic coefficient
785  T quadratic = direction1*q21n21;
786 
787  // linear coefficient
788  T linear = (direction1*q21p21) + (p1q1*q21n21);
789 
790  // constant coefficient
791  T constant = p1q1*q21p21;
792 
793  // save all zeros we find
794  std::vector<T> zeros;
795 
796  if (std::fabs(quadratic)<1e-10 && std::fabs(linear)<1e-10) {
797  return false;
798  } else if (std::fabs(quadratic)<1e-10) {
799 
800  // problem is linear
801  zeros.push_back(-constant/linear);
802 
803  } else {
804 
805  // problem is quadratic
806  T p = linear/quadratic;
807  T q = constant/quadratic;
808 
809  T sqt = 0.25*p*p -q;
810 
811  // no real solution
812  if (sqt<-1e-10)
813  return false;
814 
815  zeros.push_back(-0.5*p + std::sqrt(sqt));
816  zeros.push_back(-0.5*p -std::sqrt(sqt));
817 
818  }
819 
820  int index = -1;
821  WorldCoords r,x(-1);
822 
823  for (size_t i=0;i<zeros.size();i++) {
824 
825  T eta=zeros[i];
826 
827  // only look in the direction of the outer normals
828  if (eta<-overlap)
829  continue;
830 
831  r[2] = eta;
832 
833  // the computation of the other components might lead to nan or inf
834  // if this happens use a different equation to compute them
835  WorldCoords dummy = p1q1; dummy.axpy(eta,direction1);
836  WorldCoords c =crossProduct(dummy,q21);
837  dummy = p21; dummy.axpy(eta,n21);
838  WorldCoords d =crossProduct(q21,dummy);
839 
840  for (int j=0;j<3; j++) {
841 
842  r[0] = c[j]/d[j];
843  if (isnan(r[0]) || isinf(r[0]))
844  continue;
845 
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];
847 
848  // computation of the other components can be instable
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);
852 
853  if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
854  break;
855 
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];
858 
859  residual.axpy(r[1],q21);
860  // computation of the other components can be instable
861  if (!(isnan(r[1]) || isinf(r[1])) && residual.two_norm()<1e-5)
862  break;
863 
864  }
865  if (r[0] >= -eps && r[1]>= -eps && (r[0]<=1+eps) && (r[1] <= 1+eps)) {
866  index = i;
867  x=r;
868  }
869  }
870 
871  // TODO CLEAN THIS UP
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);
875 
876  if (residual.two_norm()<eps)
877  {
878  if (index >= 0 && x[0] >= -eps && x[1]>= -eps && (x[0]<=1+eps) && (x[1] <= 1+eps)) {
879 
880  intersection[0] = x[0];
881  intersection[1] = x[1];
882  return true;
883  }
884  return false;
885  }
886 
887  // if the direct compuation failed, use a Newton method to compute at least one zero
888 
889  // Fix some initial value
890  // sometimes it only works when the initial value is an intersection...
891  x[0] = x[1] = 0.5;
892  x[2] = 0.5;
893  WorldCoords newtonCorrection(0);
894 
895  for (int i=0; i<30; i++) {
896 
897  // compute Newton correction
898 
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);
903  Fxk.axpy(x[1],q21);
904 
905  Dune::FieldMatrix<T,3,3> FPrimexk;
906  FPrimexk[0] = p21;
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)
912  return false;
913 
914  FPrimexk.invert();
915 
916  FPrimexk.mtv(Fxk, newtonCorrection);
917 
918  x += newtonCorrection;
919  }
920 
921  residual = p1q1;
922  residual.axpy(x[0],p21); residual.axpy(x[2],direction1);
923  residual.axpy(x[2]*x[0],n21); residual.axpy(-x[1],q21);
924 
925  if (residual.two_norm()<=eps) {
926 
927  if (x[0]>=-eps && x[0]<=(1+eps) && x[1]>=-eps && x[1]<=(1+eps)) {
928  if (x[2]<-overlap)
929  return false;
930 
931  intersection[0] = x[0];
932  intersection[1] = x[1];
933 
934  return true;
935  }
936  return false;
937 
938  }
939 
940  //std::cout<<"Newton did not converge either!\n";
941 
942  return false;
943  }
944  };
945 
960  template <int dim, int dimworld, class T=double>
961  static void addEdgeIntersections(const std::vector<Dune::FieldVector<T,dimworld> >& corners1,
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)
968  {
970  directions1, gt1, gt2, polygonCorners,
971  hitCorners, neighborIntersects1,
972  neighborIntersects2, overlap);
973  }
974 
975 } // end namespace
976 #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: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 &gt1, const Dune::GeometryType &gt2, 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 &gt1, const Dune::GeometryType &gt2, 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 &gt1, Dune::GeometryType &gt2, 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 &gt1, const Dune::GeometryType &gt2, 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