dune-grid-glue  2.3.0
psurfacemerge.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:
14 #ifndef DUNE_GRIDGLUE_MERGING_PSURFACEMERGE_HH
15 #define DUNE_GRIDGLUE_MERGING_PSURFACEMERGE_HH
16 
17 
18 #include <iostream>
19 #include <fstream>
20 #include <iomanip>
21 #include <memory>
22 #include <vector>
23 #include <algorithm>
24 #include <limits>
25 
26 #include <dune/common/fvector.hh>
27 #include <dune/common/exceptions.hh>
28 #include <dune/common/bitsetvector.hh>
29 #include <dune/common/shared_ptr.hh>
30 #include <dune/common/version.hh>
31 
32 #include <dune/grid/common/grid.hh>
33 
35 
36 #include <tr1/array>
37 #if HAVE_PSURFACE
38 #include <psurface/DirectionFunction.h>
39 #include <psurface/IntersectionPrimitive.h>
40 #else
41 // forward declaration of PSurface classes
42 template <int dim, typename ctype> class DirectionFunction;
43 // switch off the macro that contains (in certain versions) the psurface namespace prefix
44 #define PSURFACE_NAMESPACE
45 #endif
46 
47 
54 template<int dim, int dimworld, typename T = double>
56  : public Merger<T,dim,dim,dimworld>
57 {
58 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,3,0)
59  static_assert( dim==1 || dim==2,
60  "PSurface can only handle the cases dim==1 and dim==2!");
61 
62  static_assert( dim==dimworld || dim+1==dimworld,
63  "PSurface can only handle the cases dim==dimworld and dim+1==dimworld!");
64 #else
65  dune_static_assert( dim==1 || dim==2,
66  "PSurface can only handle the cases dim==1 and dim==2!");
67 
68  dune_static_assert( dim==dimworld || dim+1==dimworld,
69  "PSurface can only handle the cases dim==dimworld and dim+1==dimworld!");
70 #endif
71 
72  // The psurface library itself always expects dimworld to be dim+1
73  // To be able to handle the case dim==dimworld we keep an artificial world
74  // around with the additional dimension
75  enum {psurfaceDimworld = dimworld + (dim==dimworld)};
76 
77 public:
78 
79  /* E X P O R T E D T Y P E S A N D C O N S T A N T S */
80 
82  typedef T ctype;
83 
85  typedef Dune::FieldVector<T, dimworld> WorldCoords;
86 
88  typedef Dune::FieldVector<T, dim> LocalCoords;
89 
90 
91 private:
92 
93 #if HAVE_PSURFACE
94 
95  /* P R I V A T E H E L P E R S */
96 
100  static Dune::FieldVector<ctype, dim> barycentricToReference(const Dune::FieldVector<ctype, dim+1>& bar)
101  {
102  Dune::FieldVector<ctype, dim> result;
103  for (int i=0; i<dim; i++)
104  result[i] = bar[i+1];
105 
106  return result;
107  }
108 
109 
113  static Dune::FieldVector<ctype, dim+1> referenceToBarycentric(const Dune::FieldVector<ctype, dim>& ref)
114  {
115  Dune::FieldVector<ctype, dim+1> result;
116  result[0] = 1.0;
117  for (int i=0; i<dim; i++) {
118  result[i+1] = ref[i];
119  result[0] -= ref[i];
120  }
121 
122  return result;
123  }
124 
125 
132  static bool domainParentSmaller(const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& a, const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& b)
133  {
134  return a.tris[0] < b.tris[0];
135  }
136 
143  static bool targetParentSmaller(const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* a, const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* b)
144  {
145  return a->tris[1] < b->tris[1];
146  }
147 
148 
149  /* P R I V A T E S U B C L A S S E S */
150 
155  struct OverlapManager
156  {
157 
158  OverlapManager()
159  {}
160 
161 
167  void setOverlaps(const std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype> >& unordered);
168 
172  unsigned int nOverlaps() const
173  {
174  return this->domOrder.size();
175  }
176 
183  unsigned int firstDomainParent(unsigned int parent) const;
184 
191  unsigned int firstTargetParent(unsigned int parent) const;
192 
198  const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& domain(unsigned int idx) const
199  {
200  return this->domOrder[idx];
201  }
202 
208  const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& target(unsigned int idx) const
209  {
210  return *this->tarOrder[idx];
211  }
212 
218  unsigned int domainIndex(unsigned int idx) const
219  {
220  return (this->tarOrder[idx] - this->baseptr);
221  }
222 
225  void clear()
226  {
227  purge(domOrder);
228  purge(tarOrder);
229  }
230 
231  private:
232 
233  template<typename V>
234  static void purge(V & v)
235  {
236  v.clear();
237  V v2(v);
238  v.swap(v2);
239  }
240 
243  std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype> > domOrder;
244 
247  std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>*> tarOrder;
248 
250  PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* baseptr;
251  };
252 
258  template <int dir>
259  struct ConstantDirection : public PSURFACE_NAMESPACE AnalyticDirectionFunction<psurfaceDimworld,ctype>
260  {
261  PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld> operator()(const PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld>& position) const
262  {
263  PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld> result;
264  for (size_t i=0; i<psurfaceDimworld-1; i++)
265  result[i] = 0;
266  result[psurfaceDimworld-1] = dir;
267  return result;
268  }
269  };
270 
271 private:
272 
273  /* M E M B E R V A R I A B L E S */
274 
275  /* geometric data for both domain and target */
276 
278  std::vector<std::tr1::array<double,psurfaceDimworld> > domc_;
279 
281  std::vector<std::tr1::array<double,psurfaceDimworld> > tarc_;
282 
283 
284  /* topologic information for domain and target */
285 
287  std::vector<std::tr1::array<int,dim+1> > domi_;
288 
290  std::vector<std::tr1::array<int,dim+1> > tari_;
291 
292 
293  /* members associated with the merged grid */
294 
297  OverlapManager olm_;
298 
302  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections_;
303 
312  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections_;
313 
314  bool valid;
315 
316 #endif // if HAVE_PSURFACE
317 
318 public:
319 
320  PSurfaceMerge(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
321  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections);
322 
323  PSurfaceMerge(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections = nullptr,
324  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections = nullptr);
325 
326  /* M O D E L S P E C I F I C E X T E N D I N G F U N C T I O N A L I T Y */
327 
337  inline
338  void setSurfaceDirections(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
339  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections);
340 
341  inline
342  void setSurfaceDirections(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
343  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections);
344 
345  /* C O N C E P T I M P L E M E N T I N G I N T E R F A C E */
346 
350  void build(const std::vector<Dune::FieldVector<ctype,dimworld> >& grid1_coords,
351  const std::vector<unsigned int>& grid1_elements,
352  const std::vector<Dune::GeometryType>& grid1_element_types,
353  const std::vector<Dune::FieldVector<ctype,dimworld> >& grid2_coords,
354  const std::vector<unsigned int>& grid2_elements,
355  const std::vector<Dune::GeometryType>& grid2_element_types);
356 
357 
358  /* Q U E S T I O N I N G T H E M E R G E D G R I D */
359 
362  unsigned int nSimplices() const;
363 
365  void clear()
366  {
367 #ifdef HAVE_PSURFACE
368  // Delete old internal data, from a possible previous run
369  purge(domi_);
370  purge(tari_);
371  purge(domc_);
372  purge(domi_);
373 
374  olm_.clear();
375 
376  valid = false;
377 #endif
378  }
379 
380 private:
381 
382  template<typename V>
383  static void purge(V & v)
384  {
385  v.clear();
386  V v2(v);
387  v.swap(v2);
388  }
389 
390 
391  /* M A P P I N G O N I N D E X B A S I S */
392  unsigned int grid1Parents(unsigned int idx) const;
393  unsigned int grid2Parents(unsigned int idx) const;
399  unsigned int grid1Parent(unsigned int idx, unsigned int parId = 0) const;
400 
406  unsigned int grid2Parent(unsigned int idx, unsigned int parId = 0) const;
407 
408 
409  /* G E O M E T R I C A L I N F O R M A T I O N */
410 
418  LocalCoords grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
419 
427  LocalCoords grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
428 
429 };
430 
431 #if HAVE_PSURFACE
432 
433 /* IMPLEMENTATION OF CLASS C O N T A C T M A P P I N G S U R F A C E M E R G E */
434 
435 template<int dim, int dimworld, typename T>
436 PSurfaceMerge<dim, dimworld, T>::PSurfaceMerge(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
437  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections)
438  : PSurfaceMerge(Dune::stackobject_to_shared_ptr(*domainDirections), Dune::stackobject_to_shared_ptr(*targetDirections))
439 {}
440 
441 template<int dim, int dimworld, typename T>
442 PSurfaceMerge<dim, dimworld, T>::PSurfaceMerge(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
443  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections)
444  : domainDirections_(domainDirections)
445  , targetDirections_(targetDirections)
446  , valid(false)
447 {}
448 
449 template<int dim, int dimworld, typename T>
450 inline void PSurfaceMerge<dim, dimworld, T>::setSurfaceDirections(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
451  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections)
452 {
453  setSurfaceDirections(Dune::stackobject_to_shared_ptr(*domainDirections), Dune::stackobject_to_shared_ptr(*targetDirections));
454 }
455 
456 template<int dim, int dimworld, typename T>
457 inline void PSurfaceMerge<dim, dimworld, T>::setSurfaceDirections(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
458  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections)
459 {
460  domainDirections_ = domainDirections;
461  targetDirections_ = targetDirections;
462  valid = false;
463 }
464 
465 
466 template<int dim, int dimworld, typename T>
468 {
469  assert(valid);
470  return this->olm_.nOverlaps();
471 }
472 
473 template<int dim, int dimworld, typename T>
474 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid1Parents(unsigned int idx) const
475 {
476  assert(valid);
477  return 1;
478 }
479 
480 
481 template<int dim, int dimworld, typename T>
482 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid2Parents(unsigned int idx) const
483 {
484  assert(valid);
485  return 1;
486 }
487 
488 template<int dim, int dimworld, typename T>
489 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid1Parent(unsigned int idx, unsigned int parId) const
490 {
491  assert(valid);
492  return this->olm_.domain(idx).tris[0];
493 }
494 
495 
496 template<int dim, int dimworld, typename T>
497 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid2Parent(unsigned int idx, unsigned int parId) const
498 {
499  assert(valid);
500  // Warning: Be careful to use the ACTUAL indexing here defined in the array sorted after domain parent indices!!
501  return this->olm_.domain(idx).tris[1];
502 }
503 
504 #endif // HAVE_PSURFACE
505 
506 #ifdef PSURFACE_EXTRA_TYPES
507 #define PSURFACE_EXTERN
508 #include "psurfacemerge.cc"
509 #undef PSURFACE_EXTERN
510 #endif
511 
512 #endif // DUNE_GRIDGLUE_MERGING_PSURFACEMERGE_HH
unsigned int parent(unsigned int idx, unsigned int parId=0) const
get index of grid-n's parent simplex for given merged grid simplex
Definition: merger.hh:156
void build(const std::vector< Dune::FieldVector< ctype, dimworld > > &grid1_coords, const std::vector< unsigned int > &grid1_elements, const std::vector< Dune::GeometryType > &grid1_element_types, const std::vector< Dune::FieldVector< ctype, dimworld > > &grid2_coords, const std::vector< unsigned int > &grid2_elements, const std::vector< Dune::GeometryType > &grid2_element_types)
builds the merged grid
Definition: psurfacemerge.cc:17
void setSurfaceDirections(const PSURFACE_NAMESPACE DirectionFunction< psurfaceDimworld, ctype > *domainDirections, const PSURFACE_NAMESPACE DirectionFunction< psurfaceDimworld, ctype > *targetDirections)
Set surface direction functions.
Definition: psurfacemerge.hh:450
T ctype
the numeric type used in this interface
Definition: psurfacemerge.hh:82
Definition: gridglue.hh:34
unsigned int nSimplices() const
get the number of simplices in the merged grid The indices are then in 0..nSimplices()-1 ...
Definition: psurfacemerge.hh:467
void clear()
clear the internal state, so that we can run a new merging process
Definition: psurfacemerge.hh:365
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:13
PSurfaceMerge(const PSURFACE_NAMESPACE DirectionFunction< psurfaceDimworld, ctype > *domainDirections, const PSURFACE_NAMESPACE DirectionFunction< psurfaceDimworld, ctype > *targetDirections)
Definition: psurfacemerge.hh:436
Standard implementation of the SurfaceMerge concept using the psurface library.
Definition: psurfacemerge.hh:55
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: psurfacemerge.hh:88
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: psurfacemerge.hh:85