dune-grid-glue  2.4-git
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 #include <iostream>
18 #include <fstream>
19 #include <iomanip>
20 #include <memory>
21 #include <vector>
22 #include <algorithm>
23 #include <limits>
24 
25 #include <dune/common/fvector.hh>
26 #include <dune/common/exceptions.hh>
27 #include <dune/common/bitsetvector.hh>
28 #include <dune/common/shared_ptr.hh>
29 #include <dune/common/version.hh>
30 
31 #include <dune/grid/common/grid.hh>
32 
34 
35 #include <tr1/array>
36 #if HAVE_PSURFACE
37 #include <psurface/DirectionFunction.h>
38 #include <psurface/IntersectionPrimitive.h>
39 #else
40 // forward declaration of PSurface classes
41 template <int dim, typename ctype> class DirectionFunction;
42 // switch off the macro that contains (in certain versions) the psurface namespace prefix
43 #define PSURFACE_NAMESPACE
44 #endif
45 
46 
53 template<int dim, int dimworld, typename T = double>
55  : public Merger<T,dim,dim,dimworld>
56 {
57  static_assert( dim==1 || dim==2,
58  "PSurface can only handle the cases dim==1 and dim==2!");
59 
60  static_assert( dim==dimworld || dim+1==dimworld,
61  "PSurface can only handle the cases dim==dimworld and dim+1==dimworld!");
62 
63  // The psurface library itself always expects dimworld to be dim+1
64  // To be able to handle the case dim==dimworld we keep an artificial world
65  // around with the additional dimension
66  enum {psurfaceDimworld = dimworld + (dim==dimworld)};
67 
68 public:
69 
70  /* 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 */
71 
73  typedef T ctype;
74 
76  typedef Dune::FieldVector<T, dimworld> WorldCoords;
77 
79  typedef Dune::FieldVector<T, dim> LocalCoords;
80 
81 
82 private:
83 
84 #if HAVE_PSURFACE
85 
86  /* P R I V A T E H E L P E R S */
87 
91  static Dune::FieldVector<ctype, dim> barycentricToReference(const Dune::FieldVector<ctype, dim+1>& bar)
92  {
93  Dune::FieldVector<ctype, dim> result;
94  for (int i=0; i<dim; i++)
95  result[i] = bar[i+1];
96 
97  return result;
98  }
99 
100 
104  static Dune::FieldVector<ctype, dim+1> referenceToBarycentric(const Dune::FieldVector<ctype, dim>& ref)
105  {
106  Dune::FieldVector<ctype, dim+1> result;
107  result[0] = 1.0;
108  for (int i=0; i<dim; i++) {
109  result[i+1] = ref[i];
110  result[0] -= ref[i];
111  }
112 
113  return result;
114  }
115 
116 
123  static bool domainParentSmaller(const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& a, const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& b)
124  {
125  return a.tris[0] < b.tris[0];
126  }
127 
134  static bool targetParentSmaller(const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* a, const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* b)
135  {
136  return a->tris[1] < b->tris[1];
137  }
138 
139 
140  /* P R I V A T E S U B C L A S S E S */
141 
146  struct OverlapManager
147  {
148 
149  OverlapManager()
150  {}
151 
152 
158  void setOverlaps(const std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype> >& unordered);
159 
163  unsigned int nOverlaps() const
164  {
165  return this->domOrder.size();
166  }
167 
174  unsigned int firstDomainParent(unsigned int parent) const;
175 
182  unsigned int firstTargetParent(unsigned int parent) const;
183 
189  const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& domain(unsigned int idx) const
190  {
191  return this->domOrder[idx];
192  }
193 
199  const PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>& target(unsigned int idx) const
200  {
201  return *this->tarOrder[idx];
202  }
203 
209  unsigned int domainIndex(unsigned int idx) const
210  {
211  return (this->tarOrder[idx] - this->baseptr);
212  }
213 
216  void clear()
217  {
218  purge(domOrder);
219  purge(tarOrder);
220  }
221 
222  private:
223 
224  template<typename V>
225  static void purge(V & v)
226  {
227  v.clear();
228  V v2(v);
229  v.swap(v2);
230  }
231 
234  std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype> > domOrder;
235 
238  std::vector<PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>*> tarOrder;
239 
241  PSURFACE_NAMESPACE IntersectionPrimitive<dim,ctype>* baseptr;
242  };
243 
249  template <int dir>
250  struct ConstantDirection : public PSURFACE_NAMESPACE AnalyticDirectionFunction<psurfaceDimworld,ctype>
251  {
252  PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld> operator()(const PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld>& position) const
253  {
254  PSURFACE_NAMESPACE StaticVector<ctype,psurfaceDimworld> result;
255  for (size_t i=0; i<psurfaceDimworld-1; i++)
256  result[i] = 0;
257  result[psurfaceDimworld-1] = dir;
258  return result;
259  }
260  };
261 
262 private:
263 
264  /* M E M B E R V A R I A B L E S */
265 
266  /* geometric data for both domain and target */
267 
269  std::vector<std::tr1::array<double,psurfaceDimworld> > domc_;
270 
272  std::vector<std::tr1::array<double,psurfaceDimworld> > tarc_;
273 
274 
275  /* topologic information for domain and target */
276 
278  std::vector<std::tr1::array<int,dim+1> > domi_;
279 
281  std::vector<std::tr1::array<int,dim+1> > tari_;
282 
283 
284  /* members associated with the merged grid */
285 
288  OverlapManager olm_;
289 
293  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections_;
294 
303  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections_;
304 
305  bool valid;
306 
307 #endif // if HAVE_PSURFACE
308 
309 public:
310 
311  PSurfaceMerge(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
312  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections);
313 
314  PSurfaceMerge(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections = nullptr,
315  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections = nullptr);
316 
317  /* 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 */
318 
328  inline
329  void setSurfaceDirections(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
330  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections);
331 
332  inline
333  void setSurfaceDirections(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
334  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections);
335 
336  /* 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 */
337 
341  void build(const std::vector<Dune::FieldVector<ctype,dimworld> >& grid1_coords,
342  const std::vector<unsigned int>& grid1_elements,
343  const std::vector<Dune::GeometryType>& grid1_element_types,
344  const std::vector<Dune::FieldVector<ctype,dimworld> >& grid2_coords,
345  const std::vector<unsigned int>& grid2_elements,
346  const std::vector<Dune::GeometryType>& grid2_element_types);
347 
348 
349  /* Q U E S T I O N I N G T H E M E R G E D G R I D */
350 
353  unsigned int nSimplices() const;
354 
356  void clear()
357  {
358 #ifdef HAVE_PSURFACE
359  // Delete old internal data, from a possible previous run
360  purge(domi_);
361  purge(tari_);
362  purge(domc_);
363  purge(domi_);
364 
365  olm_.clear();
366 
367  valid = false;
368 #endif
369  }
370 
371 private:
372 
373  template<typename V>
374  static void purge(V & v)
375  {
376  v.clear();
377  V v2(v);
378  v.swap(v2);
379  }
380 
381 
382  /* M A P P I N G O N I N D E X B A S I S */
383  unsigned int grid1Parents(unsigned int idx) const;
384  unsigned int grid2Parents(unsigned int idx) const;
390  unsigned int grid1Parent(unsigned int idx, unsigned int parId = 0) const;
391 
397  unsigned int grid2Parent(unsigned int idx, unsigned int parId = 0) const;
398 
399 
400  /* G E O M E T R I C A L I N F O R M A T I O N */
401 
409  LocalCoords grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
410 
418  LocalCoords grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
419 
420 };
421 
422 #if HAVE_PSURFACE
423 
424 /* 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 */
425 
426 template<int dim, int dimworld, typename T>
427 PSurfaceMerge<dim, dimworld, T>::PSurfaceMerge(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
428  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections)
429  : PSurfaceMerge(Dune::stackobject_to_shared_ptr(*domainDirections), Dune::stackobject_to_shared_ptr(*targetDirections))
430 {}
431 
432 template<int dim, int dimworld, typename T>
433 PSurfaceMerge<dim, dimworld, T>::PSurfaceMerge(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
434  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections)
435  : domainDirections_(domainDirections)
436  , targetDirections_(targetDirections)
437  , valid(false)
438 {}
439 
440 template<int dim, int dimworld, typename T>
441 inline void PSurfaceMerge<dim, dimworld, T>::setSurfaceDirections(const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* domainDirections,
442  const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype>* targetDirections)
443 {
444  setSurfaceDirections(Dune::stackobject_to_shared_ptr(*domainDirections), Dune::stackobject_to_shared_ptr(*targetDirections));
445 }
446 
447 template<int dim, int dimworld, typename T>
448 inline void PSurfaceMerge<dim, dimworld, T>::setSurfaceDirections(std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > domainDirections,
449  std::shared_ptr<const PSURFACE_NAMESPACE DirectionFunction<psurfaceDimworld,ctype> > targetDirections)
450 {
451  domainDirections_ = domainDirections;
452  targetDirections_ = targetDirections;
453  valid = false;
454 }
455 
456 
457 template<int dim, int dimworld, typename T>
459 {
460  assert(valid);
461  return this->olm_.nOverlaps();
462 }
463 
464 template<int dim, int dimworld, typename T>
465 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid1Parents(unsigned int idx) const
466 {
467  assert(valid);
468  return 1;
469 }
470 
471 
472 template<int dim, int dimworld, typename T>
473 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid2Parents(unsigned int idx) const
474 {
475  assert(valid);
476  return 1;
477 }
478 
479 template<int dim, int dimworld, typename T>
480 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid1Parent(unsigned int idx, unsigned int parId) const
481 {
482  assert(valid);
483  return this->olm_.domain(idx).tris[0];
484 }
485 
486 
487 template<int dim, int dimworld, typename T>
488 inline unsigned int PSurfaceMerge<dim, dimworld, T>::grid2Parent(unsigned int idx, unsigned int parId) const
489 {
490  assert(valid);
491  // Warning: Be careful to use the ACTUAL indexing here defined in the array sorted after domain parent indices!!
492  return this->olm_.domain(idx).tris[1];
493 }
494 
495 #endif // HAVE_PSURFACE
496 
497 #ifdef PSURFACE_EXTRA_TYPES
498 #define PSURFACE_EXTERN
499 #include "psurfacemerge.cc"
500 #undef PSURFACE_EXTERN
501 #endif
502 
503 #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:441
T ctype
the numeric type used in this interface
Definition: psurfacemerge.hh:73
unsigned int nSimplices() const
get the number of simplices in the merged grid The indices are then in 0..nSimplices()-1 ...
Definition: psurfacemerge.hh:458
void clear()
clear the internal state, so that we can run a new merging process
Definition: psurfacemerge.hh:356
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:427
Definition: gridglue.hh:33
Standard implementation of the SurfaceMerge concept using the psurface library.
Definition: psurfacemerge.hh:54
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: psurfacemerge.hh:79
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: psurfacemerge.hh:76