dune-grid-glue  2.5-git
conformingmerge.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 /*
4  * Filename: conformingmerge.hh
5  * Version: 1.0
6  * Created on: Sep 14, 2009
7  * Author: Oliver Sander
8  * ---------------------------------
9  * Project: dune-grid-glue
10  * Description: implementation of the Merger concept for conforming interfaces
11  *
12  */
19 #ifndef DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
20 #define DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
21 
22 #include <iomanip>
23 #include <vector>
24 #include <algorithm>
25 #include <bitset>
26 
27 #include <dune/common/fmatrix.hh>
28 #include <dune/common/fvector.hh>
29 
30 #include <dune/geometry/referenceelements.hh>
31 
33 
34 namespace Dune {
35 
36  namespace GridGlue {
37 
44 template<int dim, int dimworld, typename T = double>
46  : public StandardMerge<T,dim,dim,dimworld>
47 {
48 
49 public:
50 
51  /* 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 */
52 
54  typedef T ctype;
55 
57  typedef Dune::FieldVector<T, dimworld> WorldCoords;
58 
60  typedef Dune::FieldVector<T, dim> LocalCoords;
61 
62 private:
63 
64  /* M E M B E R V A R I A B L E S */
65 
67  T tolerance_;
68 
69  typedef typename StandardMerge<T,dim,dim,dimworld>::RemoteSimplicialIntersection RemoteSimplicialIntersection;
70 
75  void computeIntersections(const Dune::GeometryType& grid1ElementType,
76  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
77  std::bitset<(1<<dim)>& neighborIntersects1,
78  unsigned int grid1Index,
79  const Dune::GeometryType& grid2ElementType,
80  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
81  std::bitset<(1<<dim)>& neighborIntersects2,
82  unsigned int grid2Index,
83  std::vector<RemoteSimplicialIntersection>& intersections);
84 
85 public:
86 
87  ConformingMerge(T tolerance = 1E-4) :
88  tolerance_(tolerance)
89  {}
90 
91 private:
92 
93  /* M A P P I N G O N I N D E X B A S I S */
94 
100  unsigned int grid1Parent(unsigned int idx, unsigned int parId = 0) const;
101 
107  unsigned int grid2Parent(unsigned int idx, unsigned int parId = 0) const;
108 
109  /* G E O M E T R I C A L I N F O R M A T I O N */
110 
118  LocalCoords grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
119 
127  LocalCoords grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId = 0) const;
128 
129 };
130 
131 
132 template<int dim, int dimworld, typename T>
133 void ConformingMerge<dim, dimworld, T>::computeIntersections(const Dune::GeometryType& grid1ElementType,
134  const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
135  std::bitset<(1<<dim)>& neighborIntersects1,
136  unsigned int grid1Index,
137  const Dune::GeometryType& grid2ElementType,
138  const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
139  std::bitset<(1<<dim)>& neighborIntersects2,
140  unsigned int grid2Index,
141  std::vector<RemoteSimplicialIntersection>& intersections)
142 {
143  this->counter++;
144 
145  // A few consistency checks
146  assert((unsigned int)(Dune::ReferenceElements<T,dim>::general(grid1ElementType).size(dim)) == grid1ElementCorners.size());
147  assert((unsigned int)(Dune::ReferenceElements<T,dim>::general(grid2ElementType).size(dim)) == grid2ElementCorners.size());
148  // any intersection we may find will be the entire elements.
149  neighborIntersects1.reset();
150  neighborIntersects2.reset();
151 
152  // the intersection is either conforming or empty, hence the GeometryTypes have to match
153  if (grid1ElementType != grid2ElementType)
154  return;
155 
156  // ////////////////////////////////////////////////////////////
157  // Find correspondences between the different corners
158  // ////////////////////////////////////////////////////////////
159  std::vector<int> other(grid1ElementCorners.size(), -1);
160 
161  for (unsigned int i=0; i<grid1ElementCorners.size(); i++) {
162 
163  for (unsigned int j=0; j<grid2ElementCorners.size(); j++) {
164 
165  if ( (grid1ElementCorners[i]-grid2ElementCorners[j]).two_norm() < tolerance_ ) {
166 
167  other[i] = j;
168  break;
169 
170  }
171 
172  }
173 
174  // No corresponding grid2 vertex found for this grid1 vertex
175  if (other[i] == -1)
176  return;
177 
178  }
179 
180  // ////////////////////////////////////////////////////////////
181  // Set up the new remote intersection
182  // ////////////////////////////////////////////////////////////
183 
184  const Dune::ReferenceElement<T,dim>& refElement = Dune::ReferenceElements<T,dim>::general(grid1ElementType);
185 
187  if (grid1ElementType.isSimplex()) {
188 
189  intersections.push_back(RemoteSimplicialIntersection(grid1Index, grid2Index));
190 
191  for (int i=0; i<refElement.size(dim); i++) {
192  intersections.back().grid1Local_[0][i] = refElement.position(i,dim);
193  intersections.back().grid2Local_[0][i] = refElement.position(other[i],dim);
194  }
195 
196  } else if (grid1ElementType.isQuadrilateral()) {
197 
198  // split the quadrilateral into two triangles
199  const unsigned int subVertices[2][3] = {{0,1,3}, {0,3,2}};
200 
201  for (int i=0; i<2; i++) {
202 
203  RemoteSimplicialIntersection newSimplicialIntersection(grid1Index, grid2Index);
204 
205  for (int j=0; j<dim+1; j++) {
206  newSimplicialIntersection.grid1Local_[0][j] = refElement.position(subVertices[i][j],dim);
207  newSimplicialIntersection.grid2Local_[0][j] = refElement.position(subVertices[i][other[j]],dim);
208  }
209 
210  intersections.push_back(newSimplicialIntersection);
211 
212  }
213 
214  } else if (grid1ElementType.isHexahedron()) {
215 
216  // split the hexahedron into five tetrahedra
217  // This can be removed if ever we allow RemoteIntersections that are not simplices
218  const unsigned int subVertices[5][4] = {{0,1,3,5}, {0,3,2,6}, {4,5,0,6}, {6,7,6,3}, {6,0,5,3}};
219 
220  for (int i=0; i<5; i++) {
221 
222  RemoteSimplicialIntersection newSimplicialIntersection(grid1Index, grid2Index);
223 
224  for (int j=0; j<dim+1; j++) {
225  newSimplicialIntersection.grid1Local_[0][j] = refElement.position(subVertices[i][j],dim);
226  newSimplicialIntersection.grid2Local_[0][j] = refElement.position(subVertices[i][other[j]],dim);
227  }
228 
229  intersections.push_back(newSimplicialIntersection);
230 
231  }
232 
233  } else
234  DUNE_THROW(Dune::GridError, "Unsupported element type");
235 
236 }
237 
238 
239 template<int dim, int dimworld, typename T>
240 inline unsigned int ConformingMerge<dim, dimworld, T>::grid1Parent(unsigned int idx, unsigned int parId) const
241 {
242  return this->intersections_[idx].grid1Entities_[parId];
243 }
244 
245 
246 template<int dim, int dimworld, typename T>
247 inline unsigned int ConformingMerge<dim, dimworld, T>::grid2Parent(unsigned int idx, unsigned int parId) const
248 {
249  // Warning: Be careful to use the ACTUAL indexing here defined in the array sorted after grid1 parent indices!!
250  return this->intersections_[idx].grid2Entities_[parId];
251 }
252 
253 
254 template<int dim, int dimworld, typename T>
255 typename ConformingMerge<dim, dimworld, T>::LocalCoords ConformingMerge<dim, dimworld, T>::grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
256 {
257  return this->intersections_[idx].grid1Local_[parId][corner];
258 }
259 
260 
261 template<int dim, int dimworld, typename T>
262 typename ConformingMerge<dim, dimworld, T>::LocalCoords ConformingMerge<dim, dimworld, T>::grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
263 {
264  return this->intersections_[idx].grid2Local_[parId][corner];
265 }
266 
267 } // namespace GridGlue
268 
269 } // namespace Dune
270 
271 #endif // DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:53
Definition: gridglue.hh:31
std::vector< RemoteSimplicialIntersection > intersections_
The computed intersections.
Definition: standardmerge.hh:150
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Coordinate corner(unsigned c)
Definition: projection_impl.hh:22
T ctype
the numeric type used in this interface
Definition: conformingmerge.hh:54
Common base class for many merger implementations: produce pairs of entities that may intersect...
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: conformingmerge.hh:60
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: conformingmerge.hh:57
ConformingMerge(T tolerance=1E-4)
Definition: conformingmerge.hh:87
Implementation of the Merger concept for conforming interfaces.
Definition: conformingmerge.hh:45
unsigned int counter
Counts the number of times the computeIntersection method has been called.
Definition: merger.hh:195