193 m_enableFallback = fallback;
198 m_enableBruteForce = bruteForce;
205 bool m_enableFallback =
false;
210 bool m_enableBruteForce =
false;
214 static void purge(V & v)
228 unsigned int grid1Parents(
unsigned int idx)
const;
235 unsigned int grid2Parents(
unsigned int idx)
const;
242 unsigned int grid1Parent(
unsigned int idx,
unsigned int parId = 0)
const;
249 unsigned int grid2Parent(
unsigned int idx,
unsigned int parId = 0)
const;
261 Grid1Coords grid1ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
270 Grid2Coords grid2ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
276 void generateSeed(std::vector<int>& seeds,
277 Dune::BitSetVector<1>& isHandled2,
278 std::stack<unsigned>& candidates2,
279 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
280 const std::vector<Dune::GeometryType>& grid1_element_types,
281 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
282 const std::vector<Dune::GeometryType>& grid2_element_types);
287 int insertIntersections(
unsigned int candidate1,
unsigned int candidate2,std::vector<RemoteSimplicialIntersection>&
intersections);
292 int bruteForceSearch(
int candidate1,
293 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
294 const std::vector<Dune::GeometryType>& grid1_element_types,
295 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
296 const std::vector<Dune::GeometryType>& grid2_element_types);
301 std::pair<bool, unsigned int>
302 intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
303 RemoteSimplicialIntersection& intersection);
308 template <
int gr
idDim>
309 void computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
310 const std::vector<std::vector<unsigned int> >& gridElementCorners,
311 std::vector<std::vector<int> >& elementNeighbors);
313 void buildAdvancingFront(
314 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
315 const std::vector<unsigned int>& grid1_elements,
316 const std::vector<Dune::GeometryType>& grid1_element_types,
317 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
318 const std::vector<unsigned int>& grid2_elements,
319 const std::vector<Dune::GeometryType>& grid2_element_types
322 void buildBruteForce(
323 const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
324 const std::vector<unsigned int>& grid1_elements,
325 const std::vector<Dune::GeometryType>& grid1_element_types,
326 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
327 const std::vector<unsigned int>& grid2_elements,
328 const std::vector<Dune::GeometryType>& grid2_element_types
335template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
337 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
338 const std::vector<Dune::GeometryType>& grid1_element_types,
339 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
340 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
341 const std::vector<Dune::GeometryType>& grid2_element_types,
342 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
346 int grid1NumVertices = grid1ElementCorners_[candidate0].size();
347 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
348 for (
int i=0; i<grid1NumVertices; i++)
349 grid1ElementCorners[i] = grid1Coords[grid1ElementCorners_[candidate0][i]];
352 int grid2NumVertices = grid2ElementCorners_[candidate1].size();
353 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
354 for (
int i=0; i<grid2NumVertices; i++)
355 grid2ElementCorners[i] = grid2Coords[grid2ElementCorners_[candidate1][i]];
364 computeIntersections(grid1_element_types[candidate0], grid1ElementCorners,
365 neighborIntersects1, candidate0,
366 grid2_element_types[candidate1], grid2ElementCorners,
367 neighborIntersects2, candidate1,
375 return (
intersections.size() > 0 || neighborIntersects1.any() || neighborIntersects2.any());
379template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
381 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
382 const std::vector<Dune::GeometryType>& grid1_element_types,
383 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
384 const std::vector<Dune::GeometryType>& grid2_element_types)
386 std::bitset<(1<<grid1Dim)> neighborIntersects1;
387 std::bitset<(1<<grid2Dim)> neighborIntersects2;
388 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
390 bool intersectionFound = computeIntersection(i, candidate1,
391 grid1Coords, grid1_element_types, neighborIntersects1,
392 grid2Coords, grid2_element_types, neighborIntersects2,
396 if (intersectionFound)
405template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
407void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::
408computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
409 const std::vector<std::vector<unsigned int> >& gridElementCorners,
410 std::vector<std::vector<int> >& elementNeighbors)
412 typedef std::vector<unsigned int> FaceType;
413 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
419 elementNeighbors.resize(gridElementTypes.size());
421 for (
size_t i=0; i<gridElementTypes.size(); i++)
422 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
424 for (
size_t i=0; i<gridElementTypes.size(); i++) {
425 const Dune::ReferenceElement<T,gridDim>& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
427 for (
size_t j=0; j<(size_t)refElement.size(1); j++) {
431 for (
size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
432 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
435 std::sort(face.begin(), face.end());
437 typename FaceSetType::iterator faceHandle = faces.find(face);
439 if (faceHandle == faces.end()) {
442 faces.insert(std::make_pair(face, std::make_pair(i,j)));
447 elementNeighbors[i][j] = faceHandle->second.first;
448 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
450 faces.erase(faceHandle);
464template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
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
474 std::cout <<
"StandardMerge building merged grid..." << std::endl;
479 intersections_.clear();
489 grid1ElementCorners_.resize(grid1_element_types.size());
491 unsigned int grid1CornerCounter = 0;
493 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
496 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
497 grid1ElementCorners_[i].resize(numVertices);
498 for (
int j=0; j<numVertices; j++)
499 grid1ElementCorners_[i][j] = grid1_elements[grid1CornerCounter++];
504 grid2ElementCorners_.resize(grid2_element_types.size());
506 unsigned int grid2CornerCounter = 0;
508 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
511 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
512 grid2ElementCorners_[i].resize(numVertices);
513 for (
int j=0; j<numVertices; j++)
514 grid2ElementCorners_[i][j] = grid2_elements[grid2CornerCounter++];
522 computeNeighborsPerElement<grid1Dim>(grid1_element_types, grid1ElementCorners_, elementNeighbors1_);
523 computeNeighborsPerElement<grid2Dim>(grid2_element_types, grid2ElementCorners_, elementNeighbors2_);
525 std::cout <<
"setup took " << watch.elapsed() <<
" seconds." << std::endl;
527 if (m_enableBruteForce)
528 buildBruteForce(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
530 buildAdvancingFront(grid1Coords, grid1_elements, grid1_element_types, grid2Coords, grid2_elements, grid2_element_types);
533 std::cout <<
"intersection construction took " << watch.elapsed() <<
" seconds." << std::endl;
536template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
538 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
539 const std::vector<unsigned int>& grid1_elements,
540 const std::vector<Dune::GeometryType>& grid1_element_types,
541 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
542 const std::vector<unsigned int>& grid2_elements,
543 const std::vector<Dune::GeometryType>& grid2_element_types
550 std::stack<unsigned int> candidates1;
551 std::stack<unsigned int> candidates2;
553 std::vector<int> seeds(grid2_element_types.size(), -1);
561 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
564 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
566 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
572 std::set<unsigned int> isHandled1;
574 std::set<unsigned int> isCandidate1;
576 while (!candidates2.empty()) {
579 unsigned int currentCandidate2 = candidates2.top();
580 int seed = seeds[currentCandidate2];
584 isHandled2[currentCandidate2] =
true;
588 candidates1.push(seed);
591 isCandidate1.clear();
593 while (!candidates1.empty()) {
595 unsigned int currentCandidate1 = candidates1.top();
597 isHandled1.insert(currentCandidate1);
600 std::bitset<(1<<grid1Dim)> neighborIntersects1;
601 std::bitset<(1<<grid2Dim)> neighborIntersects2;
602 bool intersectionFound = computeIntersection(currentCandidate1, currentCandidate2,
603 grid1Coords,grid1_element_types, neighborIntersects1,
604 grid2Coords,grid2_element_types, neighborIntersects2);
606 for (
size_t i=0; i<neighborIntersects2.size(); i++)
607 if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
608 seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
611 if (intersectionFound) {
613 for (
size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
615 int neighbor = elementNeighbors1_[currentCandidate1][i];
620 if (isHandled1.find(neighbor) == isHandled1.end()
621 && isCandidate1.find(neighbor) == isCandidate1.end()) {
622 candidates1.push(neighbor);
623 isCandidate1.insert(neighbor);
637 bool seedFound = !candidates2.empty();
638 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
640 int neighbor = elementNeighbors2_[currentCandidate2][i];
646 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
648 isCandidate2[neighbor][0] =
true;
649 candidates2.push(neighbor);
654 if (seedFound || !m_enableFallback)
659 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
661 int neighbor = elementNeighbors2_[currentCandidate2][i];
666 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
673 for (
typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
674 seedIt != isHandled1.end(); ++seedIt) {
676 std::bitset<(1<<grid1Dim)> neighborIntersects1;
677 std::bitset<(1<<grid2Dim)> neighborIntersects2;
678 bool intersectionFound = computeIntersection(*seedIt, neighbor,
679 grid1Coords, grid1_element_types, neighborIntersects1,
680 grid2Coords, grid2_element_types, neighborIntersects2,
684 if (intersectionFound) {
686 Dune::dwarn <<
"Algorithm entered first fallback method and found a new seed in the build algorithm." <<
687 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
696 seed = bruteForceSearch(neighbor,
697 grid1Coords,grid1_element_types,
698 grid2Coords,grid2_element_types);
699 Dune::dwarn <<
"Algorithm entered second fallback method. This probably should not happen." << std::endl;
704 isCandidate2[neighbor] =
true;
711 candidates2.push(neighbor);
712 seeds[neighbor] = seed;
722 if (!seedFound && candidates2.empty()) {
723 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
728template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
729void StandardMerge<T,grid1Dim,grid2Dim,dimworld>::buildBruteForce(
730 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
731 const std::vector<unsigned int>& grid1_elements,
732 const std::vector<Dune::GeometryType>& grid1_element_types,
733 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
734 const std::vector<unsigned int>& grid2_elements,
735 const std::vector<Dune::GeometryType>& grid2_element_types
738 std::bitset<(1<<grid1Dim)> neighborIntersects1;
739 std::bitset<(1<<grid2Dim)> neighborIntersects2;
741 for (
unsigned i = 0; i < grid1_element_types.size(); ++i) {
742 for (
unsigned j = 0; j < grid2_element_types.size(); ++j) {
743 (void) computeIntersection(i, j,
744 grid1Coords, grid1_element_types, neighborIntersects1,
745 grid2Coords, grid2_element_types, neighborIntersects2);
750template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
754 return intersections_.size();