8 #ifndef DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
9 #define DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
20 #include <dune/common/fvector.hh>
21 #include <dune/common/bitsetvector.hh>
22 #include <dune/common/stdstreams.hh>
23 #include <dune/common/timer.hh>
24 #include <dune/common/version.hh>
26 #include <dune/geometry/referenceelements.hh>
27 #include <dune/grid/common/grid.hh>
49 template<
class T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
51 :
public Merger<T,grid1Dim,grid2Dim,dimworld>
123 const std::vector<Dune::FieldVector<T,dimworld> >& grid1ElementCorners,
124 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
125 unsigned int grid1Index,
126 const Dune::GeometryType& grid2ElementType,
127 const std::vector<Dune::FieldVector<T,dimworld> >& grid2ElementCorners,
128 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
129 unsigned int grid2Index,
130 std::vector<RemoteSimplicialIntersection>& intersections) = 0;
136 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
137 const std::vector<Dune::GeometryType>& grid1_element_types,
138 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
139 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
140 const std::vector<Dune::GeometryType>& grid2_element_types,
141 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
163 virtual void build(
const std::vector<Dune::FieldVector<T,dimworld> >& grid1_Coords,
164 const std::vector<unsigned int>& grid1_elements,
165 const std::vector<Dune::GeometryType>& grid1_element_types,
166 const std::vector<Dune::FieldVector<T,dimworld> >& grid2_coords,
167 const std::vector<unsigned int>& grid2_elements,
168 const std::vector<Dune::GeometryType>& grid2_element_types
181 purge(intersections_);
182 purge(grid1ElementCorners_);
183 purge(grid2ElementCorners_);
192 static void purge(V & v)
206 unsigned int grid1Parents(
unsigned int idx)
const;
213 unsigned int grid2Parents(
unsigned int idx)
const;
220 unsigned int grid1Parent(
unsigned int idx,
unsigned int parId = 0)
const;
227 unsigned int grid2Parent(
unsigned int idx,
unsigned int parId = 0)
const;
239 Grid1Coords grid1ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
248 Grid2Coords grid2ParentLocal(
unsigned int idx,
unsigned int corner,
unsigned int parId = 0)
const;
254 void generateSeed(std::vector<int>& seeds,
255 Dune::BitSetVector<1>& isHandled2,
256 std::stack<unsigned>& candidates2,
257 const std::vector<Dune::FieldVector<T, dimworld> >& grid1Coords,
258 const std::vector<Dune::GeometryType>& grid1_element_types,
259 const std::vector<Dune::FieldVector<T, dimworld> >& grid2Coords,
260 const std::vector<Dune::GeometryType>& grid2_element_types);
265 int insertIntersections(
unsigned int candidate1,
unsigned int candidate2,std::vector<RemoteSimplicialIntersection>& intersections);
270 int bruteForceSearch(
int candidate1,
271 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
272 const std::vector<Dune::GeometryType>& grid1_element_types,
273 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
274 const std::vector<Dune::GeometryType>& grid2_element_types);
279 int intersectionIndex(
unsigned int grid1Index,
unsigned int grid2Index,
285 template <
int gr
idDim>
286 void computeNeighborsPerElement(
const std::vector<Dune::GeometryType>& gridElementTypes,
287 const std::vector<std::vector<unsigned int> >& gridElementCorners,
288 std::vector<std::vector<int> >& elementNeighbors);
294 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
296 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
297 const std::vector<Dune::GeometryType>& grid1_element_types,
298 std::bitset<(1<<grid1Dim)>& neighborIntersects1,
299 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
300 const std::vector<Dune::GeometryType>& grid2_element_types,
301 std::bitset<(1<<grid2Dim)>& neighborIntersects2,
305 int grid1NumVertices = grid1ElementCorners_[candidate0].size();
306 std::vector<Dune::FieldVector<T,dimworld> > grid1ElementCorners(grid1NumVertices);
307 for (
int i=0; i<grid1NumVertices; i++)
308 grid1ElementCorners[i] = grid1Coords[grid1ElementCorners_[candidate0][i]];
311 int grid2NumVertices = grid2ElementCorners_[candidate1].size();
312 std::vector<Dune::FieldVector<T,dimworld> > grid2ElementCorners(grid2NumVertices);
313 for (
int i=0; i<grid2NumVertices; i++)
314 grid2ElementCorners[i] = grid2Coords[grid2ElementCorners_[candidate1][i]];
320 std::vector<RemoteSimplicialIntersection> intersections(0);
324 neighborIntersects1, candidate0,
325 grid2_element_types[candidate1], grid2ElementCorners,
326 neighborIntersects2, candidate1,
330 if(insert && intersections.size() > 0)
331 insertIntersections(candidate0,candidate1,intersections);
334 return (intersections.size() > 0 || neighborIntersects1.any() || neighborIntersects2.any());
338 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
340 const std::vector<Dune::FieldVector<T,dimworld> >& grid1Coords,
341 const std::vector<Dune::GeometryType>& grid1_element_types,
342 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
343 const std::vector<Dune::GeometryType>& grid2_element_types)
345 std::bitset<(1<<grid1Dim)> neighborIntersects1;
346 std::bitset<(1<<grid2Dim)> neighborIntersects2;
347 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
350 grid1Coords, grid1_element_types, neighborIntersects1,
351 grid2Coords, grid2_element_types, neighborIntersects2,
355 if (intersectionFound)
364 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
365 template<
int gr
idDim>
368 const std::vector<std::vector<unsigned int> >& gridElementCorners,
369 std::vector<std::vector<int> >& elementNeighbors)
371 typedef std::vector<unsigned int> FaceType;
372 typedef std::map<FaceType, std::pair<unsigned int, unsigned int> > FaceSetType;
378 elementNeighbors.resize(gridElementTypes.size());
380 for (
size_t i=0; i<gridElementTypes.size(); i++)
381 #
if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
382 elementNeighbors[i].resize(Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
384 elementNeighbors[i].resize(Dune::GenericReferenceElements<T,gridDim>::general(gridElementTypes[i]).size(1), -1);
387 for (
size_t i=0; i<gridElementTypes.size(); i++) {
389 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
390 const Dune::ReferenceElement<T,gridDim>& refElement = Dune::ReferenceElements<T,gridDim>::general(gridElementTypes[i]);
392 const Dune::GenericReferenceElement<T,gridDim>& refElement = Dune::GenericReferenceElements<T,gridDim>::general(gridElementTypes[i]);
395 for (
size_t j=0; j<(size_t)refElement.size(1); j++) {
399 for (
size_t k=0; k<(size_t)refElement.size(j,1,gridDim); k++)
400 face.push_back(gridElementCorners[i][refElement.subEntity(j,1,k,gridDim)]);
403 std::sort(face.begin(), face.end());
405 typename FaceSetType::iterator faceHandle = faces.find(face);
407 if (faceHandle == faces.end()) {
410 faces.insert(std::make_pair(face, std::make_pair(i,j)));
415 elementNeighbors[i][j] = faceHandle->second.first;
416 elementNeighbors[faceHandle->second.first][faceHandle->second.second] = i;
418 faces.erase(faceHandle);
432 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
434 const std::vector<unsigned int>& grid1_elements,
435 const std::vector<Dune::GeometryType>& grid1_element_types,
436 const std::vector<Dune::FieldVector<T,dimworld> >& grid2Coords,
437 const std::vector<unsigned int>& grid2_elements,
438 const std::vector<Dune::GeometryType>& grid2_element_types
442 std::cout <<
"StandardMerge building merged grid..." << std::endl;
447 intersections_.clear();
457 grid1ElementCorners_.resize(grid1_element_types.size());
459 unsigned int grid1CornerCounter = 0;
461 for (std::size_t i=0; i<grid1_element_types.size(); i++) {
464 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
465 int numVertices = Dune::ReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
467 int numVertices = Dune::GenericReferenceElements<T,grid1Dim>::general(grid1_element_types[i]).size(grid1Dim);
469 grid1ElementCorners_[i].resize(numVertices);
470 for (
int j=0; j<numVertices; j++)
471 grid1ElementCorners_[i][j] = grid1_elements[grid1CornerCounter++];
476 grid2ElementCorners_.resize(grid2_element_types.size());
478 unsigned int grid2CornerCounter = 0;
480 for (std::size_t i=0; i<grid2_element_types.size(); i++) {
483 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
484 int numVertices = Dune::ReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
486 int numVertices = Dune::GenericReferenceElements<T,grid2Dim>::general(grid2_element_types[i]).size(grid2Dim);
488 grid2ElementCorners_[i].resize(numVertices);
489 for (
int j=0; j<numVertices; j++)
490 grid2ElementCorners_[i][j] = grid2_elements[grid2CornerCounter++];
501 std::cout <<
"setup took " << watch.elapsed() <<
" seconds." << std::endl;
507 std::stack<unsigned int> candidates1;
508 std::stack<unsigned int> candidates2;
510 std::vector<int> seeds(grid2_element_types.size(), -1);
518 Dune::BitSetVector<1> isHandled2(grid2_element_types.size());
521 Dune::BitSetVector<1> isCandidate2(grid2_element_types.size());
523 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
529 std::set<unsigned int> isHandled1;
531 std::set<unsigned int> isCandidate1;
533 while (!candidates2.empty()) {
536 unsigned int currentCandidate2 = candidates2.top();
537 int seed = seeds[currentCandidate2];
541 isHandled2[currentCandidate2] =
true;
545 candidates1.push(seed);
548 isCandidate1.clear();
550 while (!candidates1.empty()) {
552 unsigned int currentCandidate1 = candidates1.top();
554 isHandled1.insert(currentCandidate1);
557 std::bitset<(1<<grid1Dim)> neighborIntersects1;
558 std::bitset<(1<<grid2Dim)> neighborIntersects2;
560 grid1Coords,grid1_element_types, neighborIntersects1,
561 grid2Coords,grid2_element_types, neighborIntersects2);
563 for (
size_t i=0; i<neighborIntersects2.size(); i++)
564 if (neighborIntersects2[i] && elementNeighbors2_[currentCandidate2][i] != -1)
565 seeds[elementNeighbors2_[currentCandidate2][i]] = currentCandidate1;
568 if (intersectionFound) {
570 for (
size_t i=0; i<elementNeighbors1_[currentCandidate1].size(); i++) {
572 int neighbor = elementNeighbors1_[currentCandidate1][i];
577 if (isHandled1.find(neighbor) == isHandled1.end()
578 && isCandidate1.find(neighbor) == isCandidate1.end()) {
579 candidates1.push(neighbor);
580 isCandidate1.insert(neighbor);
594 bool seedFound = !candidates2.empty();
595 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
597 int neighbor = elementNeighbors2_[currentCandidate2][i];
603 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0] && seeds[neighbor]>-1) {
605 isCandidate2[neighbor][0] =
true;
606 candidates2.push(neighbor);
616 for (
size_t i=0; i<elementNeighbors2_[currentCandidate2].size(); i++) {
618 int neighbor = elementNeighbors2_[currentCandidate2][i];
623 if (!isHandled2[neighbor][0] && !isCandidate2[neighbor][0]) {
630 for (
typename std::set<unsigned int>::iterator seedIt = isHandled1.begin();
631 seedIt != isHandled1.end(); ++seedIt) {
633 std::bitset<(1<<grid1Dim)> neighborIntersects1;
634 std::bitset<(1<<grid2Dim)> neighborIntersects2;
636 grid1Coords, grid1_element_types, neighborIntersects1,
637 grid2Coords, grid2_element_types, neighborIntersects2,
641 if (intersectionFound) {
643 Dune::dwarn <<
"Algorithm entered first fallback method and found a new seed in the build algorithm." <<
644 "Probably, the neighborIntersects bitsets computed in computeIntersection specialization is wrong." << std::endl;
653 seed = bruteForceSearch(neighbor,
654 grid1Coords,grid1_element_types,
655 grid2Coords,grid2_element_types);
656 Dune::dwarn <<
"Algorithm entered second fallback method. This probably should not happen." << std::endl;
661 isCandidate2[neighbor] =
true;
668 candidates2.push(neighbor);
669 seeds[neighbor] = seed;
679 if (!seedFound && candidates2.empty()) {
680 generateSeed(seeds, isHandled2, candidates2, grid1Coords, grid1_element_types, grid2Coords, grid2_element_types);
685 std::cout <<
"intersection construction took " << watch.elapsed() <<
" seconds." << std::endl;
689 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
693 return intersections_.size();
696 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
700 return (intersections_[idx].grid1Entities_).size();
703 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
707 return (intersections_[idx].grid2Entities_).size();
711 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
715 return intersections_[idx].grid1Entities_[parId];
719 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
723 return intersections_[idx].grid2Entities_[parId];
727 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
731 return intersections_[idx].grid1Local_[parId][corner];
735 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
739 return intersections_[idx].grid2Local_[parId][corner];
742 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
743 void 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)
745 for (std::size_t j=0; j<grid2_element_types.size(); j++) {
747 if (seeds[j] > 0 || isHandled2[j][0])
750 int seed = bruteForceSearch(j,grid1Coords,grid1_element_types,grid2Coords,grid2_element_types);
757 isHandled2[j] =
true;
761 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
763 std::vector<RemoteSimplicialIntersection>& intersections)
765 typedef typename std::vector<RemoteSimplicialIntersection>::size_type size_t;
768 for (
size_t i = 0; i < intersections.size(); ++i) {
770 int index = intersectionIndex(candidate1,candidate2,intersections[i]);
772 if (index >= this->intersections_.size()) {
773 this->intersections_.push_back(intersections[i]);
776 }
else if (index > -1) {
778 for (
size_t j = 0; j < intersections[i].grid1Entities_.size(); ++j) {
779 this->intersections_[index].grid1Entities_.push_back(candidate1);
780 this->intersections_[index].grid1Local_.push_back(intersections[i].grid1Local_[j]);
784 for (
size_t j = 0; j < intersections[i].grid2Entities_.size(); ++j) {
785 this->intersections_[index].grid2Entities_.push_back(candidate2);
786 this->intersections_[index].grid2Local_.push_back(intersections[i].grid2Local_[j]);
791 Dune::dwarn <<
"Computed the same intersection twice!" << std::endl;
797 template<
typename T,
int gr
id1Dim,
int gr
id2Dim,
int dimworld>
799 RemoteSimplicialIntersection& intersection) {
805 std::size_t n_intersections = this->intersections_.size();
808 for (std::size_t i = 0; i < n_intersections; ++i) {
811 for (std::size_t ei = 0; ei < this->intersections_[i].grid1Entities_.size(); ++ei)
813 if (this->intersections_[i].grid1Entities_[ei] == grid1Index)
815 for (std::size_t er = 0; er < intersection.grid1Entities_.size(); ++er)
817 bool found_all =
true;
819 for (std::size_t ci = 0; ci < this->intersections_[i].grid1Local_[ei].size(); ++ci)
821 Dune::FieldVector<T,grid1Dim> ni = this->intersections_[i].grid1Local_[ei][ci];
822 bool found_ni =
false;
823 for (std::size_t cr = 0; cr < intersection.grid1Local_[er].size(); ++cr)
825 Dune::FieldVector<T,grid1Dim> nr = intersection.grid1Local_[er][cr];
827 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
831 found_all = found_all && found_ni;
837 if (found_all && (this->intersections_[i].grid2Entities_[ei] != grid2Index))
846 for (std::size_t ei = 0; ei < this->intersections_[i].grid2Entities_.size(); ++ei)
848 if (this->intersections_[i].grid2Entities_[ei] == grid2Index)
850 for (std::size_t er = 0; er < intersection.grid2Entities_.size(); ++er)
852 bool found_all =
true;
854 for (std::size_t ci = 0; ci < this->intersections_[i].grid2Local_[ei].size(); ++ci)
856 Dune::FieldVector<T,grid2Dim> ni = this->intersections_[i].grid2Local_[ei][ci];
857 bool found_ni =
false;
858 for (std::size_t cr = 0; cr < intersection.grid2Local_[er].size(); ++cr)
860 Dune::FieldVector<T,grid2Dim> nr = intersection.grid2Local_[er][cr];
861 found_ni = found_ni || ((ni-nr).infinity_norm() < eps);
866 found_all = found_all && found_ni;
872 if (found_all && (this->intersections_[i].grid1Entities_[ei] != grid1Index))
881 return n_intersections;
885 #define STANDARD_MERGE_INSTANTIATE(T,A,B,C) \
887 void StandardMerge<T,A,B,C>::build(const std::vector<Dune::FieldVector<T,C> >& grid1Coords, \
888 const std::vector<unsigned int>& grid1_elements, \
889 const std::vector<Dune::GeometryType>& grid1_element_types, \
890 const std::vector<Dune::FieldVector<T,C> >& grid2Coords, \
891 const std::vector<unsigned int>& grid2_elements, \
892 const std::vector<Dune::GeometryType>& grid2_element_types \
898 #undef STANDARD_MERGE_INSTANTIATE
901 #endif // DUNE_GRIDGLUE_MERGING_STANDARDMERGE_HH
std::vector< Dune::array< Dune::FieldVector< T, grid1Dim >, nVertices > > grid1Local_
Definition: standardmerge.hh:106
std::vector< int > grid2Entities_
Definition: standardmerge.hh:114
unsigned int nSimplices() const
get the number of simplices in the merged grid The indices are then in 0..nSimplices()-1 ...
Definition: standardmerge.hh:690
#define STANDARD_MERGE_INSTANTIATE(T, A, B, C)
Definition: standardmerge.hh:885
std::vector< std::vector< unsigned int > > grid2ElementCorners_
Definition: standardmerge.hh:151
std::vector< std::vector< unsigned int > > grid1ElementCorners_
Temporary internal data.
Definition: standardmerge.hh:150
RemoteSimplicialIntersection(int grid1Entity, int grid2Entity)
Constructor for two given entity indices.
Definition: standardmerge.hh:94
Definition: standardmerge.hh:79
std::vector< Dune::array< Dune::FieldVector< T, grid2Dim >, nVertices > > grid2Local_
Definition: standardmerge.hh:109
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:295
void clear()
Definition: standardmerge.hh:178
Definition: standardmerge.hh:76
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid1Coords Grid1Coords
Type used for local coordinates on the grid1 side.
Definition: standardmerge.hh:62
std::vector< std::vector< int > > elementNeighbors1_
Definition: standardmerge.hh:153
Merger< T, grid1Dim, grid2Dim, dimworld >::Grid2Coords Grid2Coords
Type used for local coordinates on the grid2 side.
Definition: standardmerge.hh:65
StandardMerge()
Definition: standardmerge.hh:74
Abstract base for all classes that take extracted grids and build sets of intersections.
Definition: merger.hh:13
unsigned int counter
Counts the number of times the computeIntersection method has been called.
Definition: merger.hh:192
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< RemoteSimplicialIntersection > &intersections)=0
Compute the intersection between two overlapping elements.
T ctype
the numeric type used in this interface
Definition: standardmerge.hh:59
std::vector< std::vector< int > > elementNeighbors2_
Definition: standardmerge.hh:154
Common base class for many merger implementations: produce pairs of entities that may intersect...
Definition: standardmerge.hh:50
virtual 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)
builds the merged grid
Definition: standardmerge.hh:433
RemoteSimplicialIntersection()
Default constructor.
Definition: standardmerge.hh:85
std::vector< int > grid1Entities_
Definition: standardmerge.hh:112
Definition: standardmerge.hh:82
std::vector< RemoteSimplicialIntersection > intersections_
The computed intersections.
Definition: standardmerge.hh:147
Dune::FieldVector< T, dimworld > WorldCoords
the coordinate type used in this interface
Definition: standardmerge.hh:68
bool valid
Definition: standardmerge.hh:72