│ │ │
│ │ │
│ │ │
│ │ │ -
│ │ │ -
10#ifndef DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
│ │ │ -
11#define DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
21#ifndef DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
│ │ │ +
22#define DUNE_GRIDGLUE_MERGING_CONFORMINGMERGE_HH
│ │ │
│ │ │ -
24#include <dune/common/fvector.hh>
│ │ │ -
25#include <dune/common/bitsetvector.hh>
│ │ │ -
26#include <dune/common/stdstreams.hh>
│ │ │ -
27#include <dune/common/timer.hh>
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │
│ │ │ -
29#include <dune/geometry/referenceelements.hh>
│ │ │ -
30#include <dune/grid/common/grid.hh>
│ │ │ +
29#include <dune/common/fmatrix.hh>
│ │ │ +
30#include <dune/common/fvector.hh>
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
32#include <dune/geometry/referenceelements.hh>
│ │ │ +
│ │ │ +
│ │ │
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
55template<
class T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
│ │ │ -
│ │ │ -
57 :
public Merger<T,grid1Dim,grid2Dim,dimworld>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
46template<
int dim,
int dimworld,
typename T =
double>
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
100 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
│ │ │ -
101 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
│ │ │ -
102 unsigned int grid1Index,
│ │ │ -
103 const Dune::GeometryType& grid2ElementType,
│ │ │ -
104 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
│ │ │ -
105 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
│ │ │ -
106 unsigned int grid2Index,
│ │ │ -
107 std::vector<SimplicialIntersection>& intersections) = 0;
│ │ │ -
│ │ │ -
│ │ │ -
113 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
114 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
115 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
│ │ │ -
116 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
117 const std::vector<Dune::GeometryType>& grid2_element_types,
│ │ │ -
118 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
77 void computeIntersections(
const Dune::GeometryType& grid1ElementType,
│ │ │ +
78 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
│ │ │ +
79 std::bitset<(1<<dim)>& neighborIntersects1,
│ │ │ +
80 unsigned int grid1Index,
│ │ │ +
81 const Dune::GeometryType& grid2ElementType,
│ │ │ +
82 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
│ │ │ +
83 std::bitset<(1<<dim)>& neighborIntersects2,
│ │ │ +
84 unsigned int grid2Index,
│ │ │ +
85 std::vector<SimplicialIntersection>& intersections);
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
96template<
int dim,
int dimworld,
typename T>
│ │ │ +
│ │ │ +
│ │ │ +
99template<
int dim,
int dimworld,
typename T>
│ │ │ +
100void ConformingMerge<dim, dimworld, T>::computeIntersections(
const Dune::GeometryType& grid1ElementType,
│ │ │ +
101 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
│ │ │ +
102 std::bitset<(1<<dim)>& neighborIntersects1,
│ │ │ +
103 unsigned int grid1Index,
│ │ │ +
104 const Dune::GeometryType& grid2ElementType,
│ │ │ +
105 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
│ │ │ +
106 std::bitset<(1<<dim)>& neighborIntersects2,
│ │ │ +
107 unsigned int grid2Index,
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
113 assert((
unsigned int)(Dune::ReferenceElements<T,dim>::general(grid1ElementType).size(dim)) == grid1ElementCorners.size());
│ │ │ +
114 assert((
unsigned int)(Dune::ReferenceElements<T,dim>::general(grid2ElementType).size(dim)) == grid2ElementCorners.size());
│ │ │ +
│ │ │ +
116 neighborIntersects1.reset();
│ │ │ +
117 neighborIntersects2.reset();
│ │ │ +
│ │ │ +
│ │ │ +
120 if (grid1ElementType != grid2ElementType)
│ │ │ +
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
126 std::vector<int> other(grid1ElementCorners.size(), -1);
│ │ │ +
│ │ │ +
128 for (
unsigned int i=0; i<grid1ElementCorners.size(); i++) {
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
140 void build(
const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
│ │ │ -
141 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
142 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
143 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
│ │ │ -
144 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
145 const std::vector<Dune::GeometryType>& grid2_element_types)
override;
│ │ │ +
130 for (
unsigned int j=0; j<grid2ElementCorners.size(); j++) {
│ │ │ +
│ │ │ +
132 if ( (grid1ElementCorners[i]-grid2ElementCorners[j]).two_norm() < tolerance_ ) {
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
168 m_enableFallback = fallback;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
173 m_enableBruteForce = bruteForce;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
180 bool m_enableFallback =
false;
│ │ │ -
│ │ │ -
185 bool m_enableBruteForce =
false;
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
151 const auto& refElement = Dune::ReferenceElements<T,dim>::general(grid1ElementType);
│ │ │ +
│ │ │ +
154 if (grid1ElementType.isSimplex()) {
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
158 for (
int i=0; i<refElement.size(dim); i++) {
│ │ │ +
159 intersections.back().corners0[0][i] = refElement.position(i,dim);
│ │ │ +
160 intersections.back().corners1[0][i] = refElement.position(other[i],dim);
│ │ │ +
│ │ │ +
│ │ │ +
163 }
else if (dim == 2 && grid1ElementType.isQuadrilateral()) {
│ │ │ +
│ │ │ +
│ │ │ +
166 const unsigned int subVertices[2][3] = {{0,1,3}, {0,3,2}};
│ │ │ +
│ │ │ +
168 for (
int i=0; i<2; i++) {
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
172 for (
int j=0; j<dim+1; j++) {
│ │ │ +
173 newSimplicialIntersection.corners0[0][j] = refElement.position(subVertices[i][j],dim);
│ │ │ +
174 newSimplicialIntersection.corners1[0][j] = refElement.position(subVertices[i][other[j]],dim);
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
181 }
else if (grid1ElementType.isHexahedron()) {
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
185 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}};
│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
192 static void purge(V & v)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
203 void generateSeed(std::vector<int>& seeds,
│ │ │ -
204 Dune::BitSetVector<1>& isHandled2,
│ │ │ -
205 std::stack<unsigned>& candidates2,
│ │ │ -
206 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
│ │ │ -
207 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
208 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
│ │ │ -
209 const std::vector<Dune::GeometryType>& grid2_element_types);
│ │ │ -
│ │ │ -
214 int insertIntersections(
unsigned int candidate1,
unsigned int candidate2,std::vector<SimplicialIntersection>& intersections);
│ │ │ -
│ │ │ -
219 int bruteForceSearch(
int candidate1,
│ │ │ -
220 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
221 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
222 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
223 const std::vector<Dune::GeometryType>& grid2_element_types);
│ │ │ -
│ │ │ -
228 std::pair<bool, unsigned int>
│ │ │ -
229 intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
│ │ │ -
│ │ │ -
│ │ │ -
235 template <
int gr
idDim>
│ │ │ -
236 void computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
│ │ │ -
237 const std::vector<std::vector<unsigned int> >& gridElementCorners,
│ │ │ -
238 std::vector<std::vector<int> >& elementNeighbors);
│ │ │ -
│ │ │ -
240 void buildAdvancingFront(
│ │ │ -
241 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
│ │ │ -
242 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
243 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
244 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
│ │ │ -
245 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
246 const std::vector<Dune::GeometryType>& grid2_element_types
│ │ │ -
│ │ │ -
│ │ │ -
249 void buildBruteForce(
│ │ │ -
250 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
│ │ │ -
251 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
252 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
253 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
│ │ │ -
254 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
255 const std::vector<Dune::GeometryType>& grid2_element_types
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
262template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
│ │ │ -
│ │ │ -
264 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
265 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
266 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
│ │ │ -
267 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
268 const std::vector<Dune::GeometryType>& grid2_element_types,
│ │ │ -
269 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
274 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
│ │ │ -
275 for (
int i=0; i<grid1NumVertices; i++)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
280 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
│ │ │ -
281 for (
int i=0; i<grid2NumVertices; i++)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
292 neighborIntersects1, candidate0,
│ │ │ -
293 grid2_element_types[candidate1], grid2ElementCorners,
│ │ │ -
294 neighborIntersects2, candidate1,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
298 if(insert && !intersections.empty())
│ │ │ -
299 insertIntersections(candidate0,candidate1,intersections);
│ │ │ -
│ │ │ -
│ │ │ -
302 return !intersections.empty() || neighborIntersects1.any() || neighborIntersects2.any();
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
306template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
307int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::bruteForceSearch(
int candidate1,
│ │ │ -
308 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
309 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
310 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
311 const std::vector<Dune::GeometryType>& grid2_element_types)
│ │ │ -
│ │ │ -
313 std::bitset<(1<<grid1Dim)> neighborIntersects1;
│ │ │ -
314 std::bitset<(1<<grid2Dim)> neighborIntersects2;
│ │ │ -
315 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
│ │ │ -
│ │ │ -
317 bool intersectionFound = computeIntersection(i, candidate1,
│ │ │ -
318 grid1Coords, grid1_element_types, neighborIntersects1,
│ │ │ -
319 grid2Coords, grid2_element_types, neighborIntersects2,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
323 if (intersectionFound)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
332template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
│ │ │ -
334void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
│ │ │ -
335computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
│ │ │ -
336 const std::vector<std::vector<unsigned int> >& gridElementCorners,
│ │ │ -
337 std::vector<std::vector<int> >& elementNeighbors)
│ │ │ -
│ │ │ -
339 typedef std::vector<unsigned int> FaceType;
│ │ │ -
340 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
346 elementNeighbors.resize(gridElementTypes.size());
│ │ │ -
│ │ │ -
348 for (
size_t i=0; i<gridElementTypes.size(); i++)
│ │ │ -
349 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
│ │ │ -
│ │ │ -
351 for (
size_t i=0; i<gridElementTypes.size(); i++) {
│ │ │ -
352 const auto& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
│ │ │ -
│ │ │ -
354 for (
size_t j=0; j<(size_t)refElement.size(1); j++) {
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
358 for (
size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
│ │ │ -
359 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
│ │ │ -
│ │ │ -
│ │ │ -
362 std::sort(face.begin(), face.end());
│ │ │ -
│ │ │ -
364 typename FaceSetType::iterator faceHandle = faces.find(face);
│ │ │ -
│ │ │ -
366 if (faceHandle == faces.end()) {
│ │ │ -
│ │ │ -
│ │ │ -
369 faces.insert(std::make_pair(face, std::make_pair(i,j)));
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
374 elementNeighbors[i][j] = faceHandle->second.first;
│ │ │ -
375 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
│ │ │ -
│ │ │ -
377 faces.erase(faceHandle);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
391template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
│ │ │ -
│ │ │ -
393 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
394 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
395 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
396 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
397 const std::vector<Dune::GeometryType>& grid2_element_types
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
401 std::cout <<
"StandardMerge building merged grid..." << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
418 unsigned int grid1CornerCounter = 0;
│ │ │ -
│ │ │ -
420 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
│ │ │ -
│ │ │ -
│ │ │ -
423 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
│ │ │ -
│ │ │ -
425 for (
int j=0; j<numVertices; j++)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
433 unsigned int grid2CornerCounter = 0;
│ │ │ -
│ │ │ -
435 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
│ │ │ -
│ │ │ -
│ │ │ -
438 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
│ │ │ -
│ │ │ -
440 for (
int j=0; j<numVertices; j++)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
452 std::cout <<
"setup took " << watch.elapsed() <<
" seconds." << std::endl;
│ │ │ -
│ │ │ -
454 if (m_enableBruteForce)
│ │ │ -
455 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
│ │ │ -
│ │ │ -
457 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
│ │ │ -
│ │ │ -
│ │ │ -
460 std::cout <<
"intersection construction took " << watch.elapsed() <<
" seconds." << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
463template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
464void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildAdvancingFront(
│ │ │ -
465 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
466 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
467 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
468 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
469 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
470 const std::vector<Dune::GeometryType>& grid2_element_types
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
477 std::stack<unsigned int> candidates1;
│ │ │ -
478 std::stack<unsigned int> candidates2;
│ │ │ -
│ │ │ -
480 std::vector<int> seeds(grid2_element_types.size(), -1);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
488 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
│ │ │ -
│ │ │ -
│ │ │ -
491 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
│ │ │ -
│ │ │ -
493 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
499 std::set<unsigned int> isHandled1;
│ │ │ -
│ │ │ -
501 std::set<unsigned int> isCandidate1;
│ │ │ -
│ │ │ -
503 while (!candidates2.empty()) {
│ │ │ -
│ │ │ -
│ │ │ -
506 unsigned int currentCandidate2 = candidates2.top();
│ │ │ -
507 int seed = seeds[currentCandidate2];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
511 isHandled2[currentCandidate2] =
true;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
515 candidates1.push(seed);
│ │ │ -
│ │ │ -
│ │ │ -
518 isCandidate1.clear();
│ │ │ -
│ │ │ -
520 while (!candidates1.empty()) {
│ │ │ -
│ │ │ -
522 unsigned int currentCandidate1 = candidates1.top();
│ │ │ -
│ │ │ -
524 isHandled1.insert(currentCandidate1);
│ │ │ -
│ │ │ -
│ │ │ -
527 std::bitset<(1<<grid1Dim)> neighborIntersects1;
│ │ │ -
528 std::bitset<(1<<grid2Dim)> neighborIntersects2;
│ │ │ -
529 bool intersectionFound = computeIntersection(currentCandidate1, currentCandidate2,
│ │ │ -
530 grid1Coords,grid1_element_types, neighborIntersects1,
│ │ │ -
531 grid2Coords,grid2_element_types, neighborIntersects2);
│ │ │ -
│ │ │ -
533 for (
size_t i=0; i<neighborIntersects2.size(); i++)
│ │ │ -
534 if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
│ │ │ -
535 seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
│ │ │ -
│ │ │ -
│ │ │ -
538 if (intersectionFound) {
│ │ │ -
│ │ │ -
540 for (
size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
│ │ │ -
│ │ │ -
542 int neighbor = elementNeighbors1_[currentCandidate1][i];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
547 if (isHandled1.find(neighbor) == isHandled1.end()
│ │ │ -
548 && isCandidate1.find(neighbor) == isCandidate1.end()) {
│ │ │ -
549 candidates1.push(neighbor);
│ │ │ -
550 isCandidate1.insert(neighbor);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
564 bool seedFound = !candidates2.empty();
│ │ │ -
565 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
│ │ │ -
│ │ │ -
567 int neighbor = elementNeighbors2_[currentCandidate2][i];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
573 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
│ │ │ -
│ │ │ -
575 isCandidate2[neighbor][0] =
true;
│ │ │ -
576 candidates2.push(neighbor);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
581 if (seedFound || !m_enableFallback)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
586 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
│ │ │ -
│ │ │ -
588 int neighbor = elementNeighbors2_[currentCandidate2][i];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
593 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
600 for (
typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
│ │ │ -
601 seedIt != isHandled1.end(); ++seedIt) {
│ │ │ -
│ │ │ -
603 std::bitset<(1<<grid1Dim)> neighborIntersects1;
│ │ │ -
604 std::bitset<(1<<grid2Dim)> neighborIntersects2;
│ │ │ -
605 bool intersectionFound = computeIntersection(*seedIt, neighbor,
│ │ │ -
606 grid1Coords, grid1_element_types, neighborIntersects1,
│ │ │ -
607 grid2Coords, grid2_element_types, neighborIntersects2,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
611 if (intersectionFound) {
│ │ │ -
│ │ │ -
613 Dune::dwarn <<
"Algorithm entered first fallback method and found a new seed in the build algorithm." <<
│ │ │ -
614 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
623 seed = bruteForceSearch(neighbor,
│ │ │ -
624 grid1Coords,grid1_element_types,
│ │ │ -
625 grid2Coords,grid2_element_types);
│ │ │ -
626 Dune::dwarn <<
"Algorithm entered second fallback method. This probably should not happen." << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
631 isCandidate2[neighbor] =
true;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
638 candidates2.push(neighbor);
│ │ │ -
639 seeds[neighbor] = seed;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
649 if (!seedFound && candidates2.empty()) {
│ │ │ -
650 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
655template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
656void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildBruteForce(
│ │ │ -
657 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
│ │ │ -
658 const std::vector<unsigned int>& grid1_elements,
│ │ │ -
659 const std::vector<Dune::GeometryType>& grid1_element_types,
│ │ │ -
660 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
│ │ │ -
661 const std::vector<unsigned int>& grid2_elements,
│ │ │ -
662 const std::vector<Dune::GeometryType>& grid2_element_types
│ │ │ -
│ │ │ -
│ │ │ -
665 std::bitset<(1<<grid1Dim)> neighborIntersects1;
│ │ │ -
666 std::bitset<(1<<grid2Dim)> neighborIntersects2;
│ │ │ -
│ │ │ -
668 for (
unsigned i = 0; i < grid1_element_types.size(); ++i) {
│ │ │ -
669 for (
unsigned j = 0; j < grid2_element_types.size(); ++j) {
│ │ │ -
670 (void) computeIntersection(i, j,
│ │ │ -
671 grid1Coords, grid1_element_types, neighborIntersects1,
│ │ │ -
672 grid2Coords, grid2_element_types, neighborIntersects2);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
677template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
678void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::generateSeed(std::vector<int>& seeds, Dune::BitSetVector<1>& isHandled2, std::stack<unsigned>& candidates2,
const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
const std::vector<Dune::GeometryType>& grid1_element_types,
const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
const std::vector<Dune::GeometryType>& grid2_element_types)
│ │ │ -
│ │ │ -
680 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
│ │ │ -
│ │ │ -
682 if (seeds[j] > 0 || isHandled2[j][0])
│ │ │ -
│ │ │ -
│ │ │ -
685 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
692 isHandled2[j] =
true;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
696template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
697int StandardMerge<T,grid1Dim,grid2Dim,dimworld>::insertIntersections(
unsigned int candidate1,
unsigned int candidate2,
│ │ │ -
│ │ │ -
│ │ │ -
700 typedef typename std::vector<SimplicialIntersection>::size_type size_t;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
707 std::tie(found, index) = intersectionIndex(candidate1,candidate2,
intersections[i]);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
717 for (
size_t j = 0; j <
intersections[i].parents0.size(); ++j) {
│ │ │ -
718 intersection.parents0.push_back(candidate1);
│ │ │ -
719 intersection.corners0.push_back(
intersections[i].corners0[j]);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
723 for (
size_t j = 0; j <
intersections[i].parents1.size(); ++j) {
│ │ │ -
724 intersection.parents1.push_back(candidate2);
│ │ │ -
725 intersection.corners1.push_back(
intersections[i].corners1[j]);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
730 Dune::dwarn <<
"Computed the same intersection twice!" << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
736template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
│ │ │ -
737std::pair<bool, unsigned int>
│ │ │ -
738StandardMerge<T,grid1Dim,grid2Dim,dimworld>::intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
746 if (grid1Dim == grid2Dim)
│ │ │ -
747 return {
true, n_intersections};
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
751 for (std::size_t i = 0; i < n_intersections; ++i) {
│ │ │ -
│ │ │ -
│ │ │ -
754 for (std::size_t ei = 0; ei < this->
intersections()[i].parents0.size(); ++ei)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
758 for (std::size_t er = 0; er < intersection.parents0.size(); ++er)
│ │ │ -
│ │ │ -
760 bool found_all =
true;
│ │ │ -
│ │ │ -
762 for (std::size_t ci = 0; ci < this->
intersections()[i].corners0[ei].size(); ++ci)
│ │ │ -
│ │ │ -
764 Dune::FieldVector<T,grid1Dim> ni = this->
intersections()[i].corners0[ei][ci];
│ │ │ -
765 bool found_ni =
false;
│ │ │ -
766 for (std::size_t cr = 0; cr < intersection.corners0[er].size(); ++cr)
│ │ │ -
│ │ │ -
768 Dune::FieldVector<T,grid1Dim> nr = intersection.corners0[er][cr];
│ │ │ -
│ │ │ -
770 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
774 found_all = found_all && found_ni;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
780 if (found_all && (this->
intersections()[i].parents1[ei] != grid2Index))
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
789 for (std::size_t ei = 0; ei < this->
intersections()[i].parents1.size(); ++ei)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
793 for (std::size_t er = 0; er < intersection.parents1.size(); ++er)
│ │ │ -
│ │ │ -
795 bool found_all =
true;
│ │ │ -
│ │ │ -
797 for (std::size_t ci = 0; ci < this->
intersections()[i].corners1[ei].size(); ++ci)
│ │ │ -
│ │ │ -
799 Dune::FieldVector<T,grid2Dim> ni = this->
intersections()[i].corners1[ei][ci];
│ │ │ -
800 bool found_ni =
false;
│ │ │ -
801 for (std::size_t cr = 0; cr < intersection.corners1[er].size(); ++cr)
│ │ │ -
│ │ │ -
803 Dune::FieldVector<T,grid2Dim> nr = intersection.corners1[er][cr];
│ │ │ -
804 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
809 found_all = found_all && found_ni;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
815 if (found_all && (this->
intersections()[i].parents0[ei] != grid1Index))
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
824 return {
true, n_intersections};
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
828#define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \
│ │ │ -
│ │ │ -
830 void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \
│ │ │ -
831 const std::vector<unsigned int>& grid1_elements, \
│ │ │ -
832 const std::vector<Dune::GeometryType>& grid1_element_types, \
│ │ │ -
833 const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \
│ │ │ -
834 const std::vector<unsigned int>& grid2_elements, \
│ │ │ -
835 const std::vector<Dune::GeometryType>& grid2_element_types \
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
841#undef STANDARD_MERGE_INSTANTIATE
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
#define STANDARD_MERGE_INSTANTIATE(T, A, B, C)
Definition standardmerge.cc:13
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ +
187 for (
int i=0; i<5; i++) {
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
191 for (
int j=0; j<dim+1; j++) {
│ │ │ +
192 newSimplicialIntersection.corners0[0][j] = refElement.position(subVertices[i][j],dim);
│ │ │ +
193 newSimplicialIntersection.corners1[0][j] = refElement.position(subVertices[i][other[j]],dim);
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
201 DUNE_THROW(Dune::GridError,
"Unsupported element type");
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
Common base class for many merger implementations: produce pairs of entities that may intersect.
│ │ │
Definition gridglue.hh:37
│ │ │
Definition gridglue.hh:38
│ │ │
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
│ │ │ -
Definition intersectionlist.hh:207
│ │ │ -
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition merger.hh:27
│ │ │ -
Dune::FieldVector< T, dimworld > WorldCoords
Definition merger.hh:37
│ │ │ -
Dune::GridGlue::IntersectionList< Grid1Coords, Grid2Coords > IntersectionList
Definition merger.hh:39
│ │ │ -
Dune::FieldVector< T, grid1Dim > Grid1Coords
Definition merger.hh:31
│ │ │ -
unsigned int counter
Definition merger.hh:114
│ │ │ -
Dune::FieldVector< T, grid2Dim > Grid2Coords
Definition merger.hh:34
│ │ │ +
Dune::FieldVector< T, dim > LocalCoords
the coordinate type used in this interface
Definition conformingmerge.hh:62
│ │ │ +
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition conformingmerge.hh:59
│ │ │ +
static constexpr T default_tolerance
Definition conformingmerge.hh:89
│ │ │ +
T ctype
the numeric type used in this interface
Definition conformingmerge.hh:56
│ │ │ +
ConformingMerge(T tolerance=default_tolerance)
Definition conformingmerge.hh:91
│ │ │
│ │ │ -
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< SimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
│ │ │ -
std::shared_ptr< IntersectionList > intersectionList_
Definition standardmerge.hh:124
│ │ │ -
typename Base::Grid1Coords Grid1Coords
Definition standardmerge.hh:69
│ │ │ -
void enableFallback(bool fallback)
Definition standardmerge.hh:166
│ │ │ -
std::vector< std::vector< int > > elementNeighbors2_
Definition standardmerge.hh:131
│ │ │ -
std::vector< std::vector< int > > elementNeighbors1_
Definition standardmerge.hh:130
│ │ │ -
T ctype
Definition standardmerge.hh:66
│ │ │ -
void clear() override
Definition standardmerge.hh:150
│ │ │ -
typename Base::IntersectionList IntersectionList
Definition standardmerge.hh:77
│ │ │ -
bool computeIntersection(unsigned int candidate0, unsigned int candidate1, const std::vector< Dune::FieldVector< T, dimworld > > &grid1Coords, const std::vector< Dune::GeometryType > &grid1_element_types, std::bitset<(1<< grid1Dim)> &neighborIntersects1, const std::vector< Dune::FieldVector< T, dimworld > > &grid2Coords, const std::vector< Dune::GeometryType > &grid2_element_types, std::bitset<(1<< grid2Dim)> &neighborIntersects2, bool insert=true)
Compute the intersection between two overlapping elements.
Definition standardmerge.hh:263
│ │ │ -
void enableBruteForce(bool bruteForce)
Definition standardmerge.hh:171
│ │ │ -
std::shared_ptr< IntersectionListProvider > intersectionListProvider_
Definition standardmerge.hh:123
│ │ │ -
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) override
builds the merged grid
Definition standardmerge.hh:392
│ │ │ -
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Definition standardmerge.hh:127
│ │ │ -
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition standardmerge.hh:128
│ │ │ -
typename Base::Grid2Coords Grid2Coords
Definition standardmerge.hh:72
│ │ │ -
std::shared_ptr< IntersectionList > intersectionList() const final
Definition standardmerge.hh:160
│ │ │ -
virtual ~StandardMerge()=default
│ │ │ -
SimplicialIntersection RemoteSimplicialIntersection
Definition standardmerge.hh:84
│ │ │ -
SimplicialIntersectionListProvider< grid1Dim, grid2Dim > IntersectionListProvider
Definition standardmerge.hh:82
│ │ │ -
bool valid
Definition standardmerge.hh:86
│ │ │ -
StandardMerge()
Definition standardmerge.hh:88
│ │ │ -
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition standardmerge.hh:83
│ │ │ -
typename Base::WorldCoords WorldCoords
Definition standardmerge.hh:75
│ │ │ +
StandardMerge()
Definition standardmerge.hh:88
│ │ │ +
typename IntersectionListProvider::SimplicialIntersection SimplicialIntersection
Definition standardmerge.hh:83
│ │ │