dune-grid-glue  2.3.0
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 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
147  assert((unsigned int)(Dune::ReferenceElements<T,dim>::general(grid1ElementType).size(dim)) == grid1ElementCorners.size());
148  assert((unsigned int)(Dune::ReferenceElements<T,dim>::general(grid2ElementType).size(dim)) == grid2ElementCorners.size());
149 #else
150  assert((unsigned int)(Dune::GenericReferenceElements<T,dim>::general(grid1ElementType).size(dim)) == grid1ElementCorners.size());
151  assert((unsigned int)(Dune::GenericReferenceElements<T,dim>::general(grid2ElementType).size(dim)) == grid2ElementCorners.size());
152 #endif
153  // any intersection we may find will be the entire elements.
154  neighborIntersects1.reset();
155  neighborIntersects2.reset();
156 
157  // the intersection is either conforming or empty, hence the GeometryTypes have to match
158  if (grid1ElementType != grid2ElementType)
159  return;
160 
161  // ////////////////////////////////////////////////////////////
162  // Find correspondences between the different corners
163  // ////////////////////////////////////////////////////////////
164  std::vector<int> other(grid1ElementCorners.size(), -1);
165 
166  for (unsigned int i=0; i<grid1ElementCorners.size(); i++) {
167 
168  for (unsigned int j=0; j<grid2ElementCorners.size(); j++) {
169 
170  if ( (grid1ElementCorners[i]-grid2ElementCorners[j]).two_norm() < tolerance_ ) {
171 
172  other[i] = j;
173  break;
174 
175  }
176 
177  }
178 
179  // No corresponding grid2 vertex found for this grid1 vertex
180  if (other[i] == -1)
181  return;
182 
183  }
184 
185  // ////////////////////////////////////////////////////////////
186  // Set up the new remote intersection
187  // ////////////////////////////////////////////////////////////
188 
189 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
190  const Dune::ReferenceElement<T,dim>& refElement = Dune::ReferenceElements<T,dim>::general(grid1ElementType);
191 #else
192  const Dune::GenericReferenceElement<T,dim>& refElement = Dune::GenericReferenceElements<T,dim>::general(grid1ElementType);
193 #endif
194 
196  if (grid1ElementType.isSimplex()) {
197 
198  intersections.push_back(RemoteSimplicialIntersection(grid1Index, grid2Index));
199 
200  for (int i=0; i<refElement.size(dim); i++) {
201  intersections.back().grid1Local_[0][i] = refElement.position(i,dim);
202  intersections.back().grid2Local_[0][i] = refElement.position(other[i],dim);
203  }
204 
205  } else if (grid1ElementType.isQuadrilateral()) {
206 
207  // split the quadrilateral into two triangles
208  const unsigned int subVertices[2][3] = {{0,1,3}, {0,3,2}};
209 
210  for (int i=0; i<2; i++) {
211 
212  RemoteSimplicialIntersection newSimplicialIntersection(grid1Index, grid2Index);
213 
214  for (int j=0; j<dim+1; j++) {
215  newSimplicialIntersection.grid1Local_[0][j] = refElement.position(subVertices[i][j],dim);
216  newSimplicialIntersection.grid2Local_[0][j] = refElement.position(subVertices[i][other[j]],dim);
217  }
218 
219  intersections.push_back(newSimplicialIntersection);
220 
221  }
222 
223  } else if (grid1ElementType.isHexahedron()) {
224 
225  // split the hexahedron into five tetrahedra
226  // This can be removed if ever we allow RemoteIntersections that are not simplices
227  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}};
228 
229  for (int i=0; i<5; i++) {
230 
231  RemoteSimplicialIntersection newSimplicialIntersection(grid1Index, grid2Index);
232 
233  for (int j=0; j<dim+1; j++) {
234  newSimplicialIntersection.grid1Local_[0][j] = refElement.position(subVertices[i][j],dim);
235  newSimplicialIntersection.grid2Local_[0][j] = refElement.position(subVertices[i][other[j]],dim);
236  }
237 
238  intersections.push_back(newSimplicialIntersection);
239 
240  }
241 
242  } else
243  DUNE_THROW(Dune::GridError, "Unsupported element type");
244 
245 }
246 
247 
248 template<int dim, int dimworld, typename T>
249 inline unsigned int ConformingMerge<dim, dimworld, T>::grid1Parent(unsigned int idx, unsigned int parId) const
250 {
251  return this->intersections_[idx].grid1Entities_[parId];
252 }
253 
254 
255 template<int dim, int dimworld, typename T>
256 inline unsigned int ConformingMerge<dim, dimworld, T>::grid2Parent(unsigned int idx, unsigned int parId) const
257 {
258  // Warning: Be careful to use the ACTUAL indexing here defined in the array sorted after grid1 parent indices!!
259  return this->intersections_[idx].grid2Entities_[parId];
260 }
261 
262 
263 template<int dim, int dimworld, typename T>
264 typename ConformingMerge<dim, dimworld, T>::LocalCoords ConformingMerge<dim, dimworld, T>::grid1ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
265 {
266  return this->intersections_[idx].grid1Local_[parId][corner];
267 }
268 
269 
270 template<int dim, int dimworld, typename T>
271 typename ConformingMerge<dim, dimworld, T>::LocalCoords ConformingMerge<dim, dimworld, T>::grid2ParentLocal(unsigned int idx, unsigned int corner, unsigned int parId) const
272 {
273  return this->intersections_[idx].grid2Local_[parId][corner];
274 }
275 
276 } // namespace GridGlue
277 
278 } // namespace Dune
279 
280 #endif // DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
T ctype
the numeric type used in this interface
Definition: conformingmerge.hh:54
Definition: gridglue.hh:34
ConformingMerge(T tolerance=1E-4)
Definition: conformingmerge.hh:87
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: conformingmerge.hh:57
Implementation of the Merger concept for conforming interfaces.
Definition: conformingmerge.hh:45
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:50
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition: conformingmerge.hh:60
Common base class for many merger implementations: produce pairs of entities that may intersect...