dune-grid-glue  2.3.0
contactmerge.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:
8 #ifndef DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
9 #define DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
10 
11 
12 #include <iostream>
13 #include <fstream>
14 #include <iomanip>
15 #include <vector>
16 #include <algorithm>
17 #include <limits>
18 
19 #include <dune/common/fvector.hh>
20 #include <dune/common/function.hh>
21 #include <dune/common/exceptions.hh>
22 #include <dune/common/bitsetvector.hh>
23 
24 #include <dune/grid/common/grid.hh>
25 
29 
35 template<int dimworld, typename T = double>
37 : public StandardMerge<T,dimworld-1,dimworld-1,dimworld>
38 {
39  enum {dim = dimworld-1};
40 
41  static_assert( dim==1 || dim==2,
42  "ContactMerge yet only handles the cases dim==1 and dim==2!");
43 
45 public:
46 
47  /* 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 */
48 
50  typedef T ctype;
51 
53  typedef Dune::FieldVector<T, dimworld> WorldCoords;
54 
56  typedef Dune::FieldVector<T, dim> LocalCoords;
57 
58  ContactMerge(const T allowedOverlap=T(0),
59  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections = NULL,
60  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections = NULL)
61  : domainDirections_(domainDirections), targetDirections_(targetDirections),
62  overlap_(allowedOverlap)
63  {}
64 
73  inline
74  void setSurfaceDirections(const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections,
75  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections)
76  {
77  domainDirections_ = domainDirections;
78  targetDirections_ = targetDirections;
79  this->valid = false;
80  }
81 
83  void setOverlap(T overlap)
84  {
85  overlap_ = overlap;
86  }
87 
88 protected:
89  typedef typename StandardMerge<T,dimworld-1,dimworld-1,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
90 
91 private:
95  const Dune::VirtualFunction<WorldCoords,WorldCoords>* domainDirections_;
96  std::vector<WorldCoords> nodalDomainDirections_;
97 
106  const Dune::VirtualFunction<WorldCoords,WorldCoords>* targetDirections_;
107  std::vector<WorldCoords> nodalTargetDirections_;
108 
110  T overlap_;
111 
116  void computeIntersections(const Dune::GeometryType& grid1ElementType,
117  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
118  std::bitset<(1<<dim)>& neighborIntersects1,
119  unsigned int grid1Index,
120  const Dune::GeometryType& grid2ElementType,
121  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
122  std::bitset<(1<<dim)>& neighborIntersects2,
123  unsigned int grid2Index,
124  std::vector<RemoteSimplicialIntersection>& intersections);
125 
129  void build(const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
130  const std::vector<unsigned int>& grid1Elements,
131  const std::vector<Dune::GeometryType>& grid1ElementTypes,
132  const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
133  const std::vector<unsigned int>& grid2Elements,
134  const std::vector<Dune::GeometryType>& grid2ElementTypes)
135  {
136  std::cout<<"ContactMerge building grid!\n";
137  // setup the nodal direction vectors
138  setupNodalDirections(grid1Coords, grid1Elements, grid1ElementTypes,
139  grid2Coords, grid2Elements, grid2ElementTypes);
140 
141  Base::build(grid1Coords, grid1Elements, grid1ElementTypes,
142  grid2Coords, grid2Elements, grid2ElementTypes);
143 
144  }
145 
146 private:
147 
149  static LocalCoords localCornerCoords(int i, const Dune::GeometryType& gt)
150  {
151  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
152  return ref.position(i,dim);
153  }
154 
156  static int isCorner(const Dune::GeometryType& gt, const LocalCoords& local)
157  {
158  const Dune::ReferenceElement<T,dim>& ref = Dune::ReferenceElements<T,dim>::general(gt);
159 
160  for (int i=0; i<ref.size(dim); i++)
161  if ((ref.position(i,dim)-local).two_norm()<1e-8)
162  return i;
163  return -1;
164  }
165 
167  static std::vector<T> localToBarycentric(const Dune::GeometryType& gt, const LocalCoords& local)
168  {
169  std::vector<T> bar(Dune::ReferenceElements<T,dim>::general(gt).size(dim));
170  if (bar.size()<4) {
171  bar[0] = 1.0;
172  for (int i=0; i<dim; i++) {
173  bar[i+1] = local[i];
174  bar[0] -= local[i];
175  }
176  } else {
177  // quadrilateral case
178  for (int i=0; i<4; i++) {
179  bar[i] = 1;
180 
181  for (int j=0; j<dim; j++)
182  // if j-th bit of i is set multiply with in[j], else with 1-in[j]
183  bar[i] *= (i & (1<<j)) ? local[j] : 1-local[j];
184  }
185  }
186  return bar;
187  }
188 
190  static WorldCoords interpolate(const std::vector<WorldCoords>& p,
191  const Dune::GeometryType& gt, const LocalCoords& local)
192  {
193  // Compute barycentric coordinates
194  std::vector<T> bar = localToBarycentric(gt, local);
195 
196  WorldCoords global(0);
197  for (size_t i=0; i<p.size(); i++)
198  global.axpy(bar[i],p[i]);
199 
200  return global;
201  }
202 
204  void computeCyclicOrder(const std::vector<Dune::array<LocalCoords,2> >& polytopeCorners,
205  const LocalCoords& center, std::vector<int>& ordering) const;
206 
208  void setupNodalDirections(const std::vector<WorldCoords>& coords1,
209  const std::vector<unsigned int>& elements1,
210  const std::vector<Dune::GeometryType>& elementTypes1,
211  const std::vector<WorldCoords>& coords2,
212  const std::vector<unsigned int>& elements2,
213  const std::vector<Dune::GeometryType>& elementTypes2);
214 
216  void computeOuterNormalField(const std::vector<WorldCoords>& coords,
217  const std::vector<unsigned int>& elements,
218  const std::vector<Dune::GeometryType>& elementTypes,
219  std::vector<WorldCoords>& normals);
220 
222  void removeDoubles(std::vector<Dune::array<LocalCoords,2> >& polytopeCorners);
223 };
224 
225 #include "contactmerge.cc"
226 
227 #endif // DUNE_GRIDGLUE_MERGING_CONTACTMERGE_HH
Central component of the module implementing the coupling of two grids.
Merge two codimension-1 surfaces that may be a positive distance apart.
Definition: contactmerge.hh:36
virtual void computeIntersections(const Dune::GeometryType &grid1ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid1ElementCorners, std::bitset<(1<< grid1Dim)> &neighborIntersects1, unsigned int grid1Index, const Dune::GeometryType &grid2ElementType, const std::vector< Dune::FieldVector< T, dimworld > > &grid2ElementCorners, std::bitset<(1<< grid2Dim)> &neighborIntersects2, unsigned int grid2Index, std::vector< RemoteSimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:59
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:50
virtual void build(const std::vector< Dune::FieldVector< T, dimworld > > &grid1_Coords, const std::vector< unsigned int > &grid1_elements, const std::vector< Dune::GeometryType > &grid1_element_types, const std::vector< Dune::FieldVector< T, dimworld > > &grid2_coords, const std::vector< unsigned int > &grid2_elements, const std::vector< Dune::GeometryType > &grid2_element_types)
builds the merged grid
Common base class for many merger implementations: produce pairs of entities that may intersect...
Base class for predicates selecting the part of a grid to be extracted.
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:68
bool valid
Definition: standardmerge.hh:72