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
34namespace Dune {
35
36 namespace GridGlue {
37
44template<int dim, int dimworld, typename T = double>
46 : public StandardMerge<T,dim,dim,dimworld>
47{
48
49public:
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
62private:
63
64 /* M E M B E R V A R I A B L E S */
65
67 T tolerance_;
68
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
85public:
86
87 ConformingMerge(T tolerance = 1E-4) :
88 tolerance_(tolerance)
89 {}
90
91private:
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
132template<int dim, int dimworld, typename T>
133void 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
239template<int dim, int dimworld, typename T>
240inline 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
246template<int dim, int dimworld, typename T>
247inline 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
254template<int dim, int dimworld, typename T>
255typename 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
261template<int dim, int dimworld, typename T>
262typename 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 gridglue.hh:30
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
Implementation of the Merger concept for conforming interfaces.
Definition conformingmerge.hh:47
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
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.
Definition standardmerge.hh:55