1 #ifndef ALUGRID_COMPATIBILITY_CHECK_HH_INCLUDED
2 #define ALUGRID_COMPATIBILITY_CHECK_HH_INCLUDED
13 #include <dune/common/timer.hh>
39 template <
class VertexVector>
93 const std::vector<ElementType>& elements,
110 std::cout <<
"Build neighbors took " << timer.elapsed() <<
" sec." << std::endl;
117 bool verbose =
false;
118 unsigned int bndFaces = 0;
119 std::vector<int> nonCompatFacesAtVertex(
nVertices_, 0 );
122 if( face.second[0] == face.second[1] )
126 else if(! checkStrongCompatibility(face, verbose))
129 for(
size_t i =0; i < 3; ++i)
131 ++(nonCompatFacesAtVertex[face.first[i]]);
135 std::cout <<
"NotStrongCompatibleMacroFaces" <<
" InnerFaces " <<
" TotalFaces " <<
"Maximum/Vertex " <<
" Minimum/Vertex " << std::endl;
136 std::cout << result <<
" " <<
neighbours_.size() - bndFaces <<
" " <<
neighbours_.size() <<
" " << *(std::max_element(nonCompatFacesAtVertex.begin(), nonCompatFacesAtVertex.end())) <<
" " << *(std::min_element(nonCompatFacesAtVertex.begin(), nonCompatFacesAtVertex.end())) << std::endl << std::endl;
146 if(!checkFaceCompatibility(face))
return false;
154 applyStandardOrientation();
157 applyStandardOrientation();
166 std::cout << face.first[0] <<
"," << face.first[1] <<
"," << face.first[2] <<
" : " << face.second[0] <<
", " << face.second[1] << std::endl;
175 std::cout <<
"[" << el[0] ;
176 for(
int i=1; i < 4; ++i)
177 std::cout <<
","<< el[i] ;
178 std::cout <<
"] Refinement Edges: ";
179 for(
int i=0; i< 4; ++i)
181 getRefinementEdge(el, i, edge,
types_[index]);
182 std::cout <<
"[" << edge[0] <<
"," << edge[1] <<
"] ";
184 std::cout <<
" Type: " <<
types_[ index ] <<
" V1: ";
185 for (
size_t i = 0; i < 4; ++i) {
187 std::cout << el[i] <<
"," ;
189 std::cout << std::endl;
207 std::list<std::pair<FaceType, EdgeType> > activeFaceList;
208 std::vector<bool> doneElements(
elements_.size(),
false);
210 std::vector< std::pair< double, std::pair< int,int > > > vertexOrder(
nVertices_ , std::make_pair(-1.0, std::make_pair(-1,-2) ) );
211 const double eps = std::numeric_limits< double >::epsilon() * double(
nVertices_) * 10.0;
214 const unsigned int l1 = -1;
216 const int numberOfElements =
elements_.size();
218 for(
int counter = 0; counter < numberOfElements ; ++counter)
229 for(
unsigned int i=0 ; i < 4 ; ++i)
232 vertexOrder[ vtx ].first = i+1;
234 vertexOrder[ vtx ].second.first = i > 0 ? el0[ i-1 ] : -1;
235 vertexOrder[ vtx ].second.second = i < 3 ? el0[ i+1 ] : -2;
240 assert( ! activeFaceList.empty() );
242 auto it = activeFaceList.begin();
243 int el1 = it->second[ 0 ];
244 int el2 = it->second[ 1 ];
246 while( doneElements[ el1 ] && doneElements[ el2 ] )
248 activeFaceList.erase( it++ );
269 assert( it != activeFaceList.end() );
270 el1 = it->second[ 0 ];
271 el2 = it->second[ 1 ];
278 int elIndex = faceElement.second[0];
281 int neighIndex = faceElement.second[1];
282 if(!doneElements[elIndex])
284 assert(doneElements[neighIndex] || counter == 0 );
285 std::swap(elIndex, neighIndex);
288 assert(!doneElements[neighIndex] || counter == 0);
289 doneElements[neighIndex] =
true;
292 unsigned int faceInEl = getFaceIndex(el, faceElement.first);
293 unsigned int nodeInEl = 3 - faceInEl;
294 unsigned int faceInNeigh = getFaceIndex(neigh, faceElement.first);
295 unsigned int nodeInNeigh = 3 - faceInNeigh;
302 if( vertexOrder[ neigh [ nodeInNeigh ] ].first < 0 )
304 int vxIndex = el[ nodeInEl ];
310 auto& vxPair = vertexOrder[ el[ nodeInNeigh ] ];
312 vxIndex = vxPair.second.second;
315 vertexOrder[ neigh [ nodeInNeigh ] ] = std::make_pair( vxPair.first+1.0, std::make_pair( el[ nodeInNeigh ], -2 ) );
316 vxPair.second.second = neigh [ nodeInNeigh ];
325 vxIndex = el[ nodeInNeigh ];
328 auto& vxPair = vertexOrder[ vxIndex ];
329 assert( vxPair.first >= 0 );
332 const int prevIdx = vxPair.second.first ;
333 double prevOrder = ( prevIdx == -1 ) ? 0.0 : vertexOrder[ prevIdx ].first;
334 double newOrder = 0.5 * (vxPair.first + prevOrder);
336 vertexOrder[ neigh [ nodeInNeigh ] ] = std::make_pair( newOrder, std::make_pair( prevIdx, vxIndex ) );
338 vertexOrder[ prevIdx ].second.second = neigh [ nodeInNeigh ];
339 vxPair.second.first = neigh [ nodeInNeigh ];
340 assert( vxPair.first > newOrder );
342 if( (vxPair.first - newOrder) < eps )
345 Dune::Timer restimer;
346 std::cout <<
"Rescale vertex order weights." << std::endl;
348 std::map< double, int > newVertexOrder;
349 const int size = vertexOrder.size();
350 for(
int j=0; j<size; ++j )
352 if( vertexOrder[ j ].first > 0.0 )
354 newVertexOrder.insert( std::make_pair( vertexOrder[ j ].first, j ) );
359 for(
const auto& vx: newVertexOrder )
361 assert( vx.second >=0 && vx.second < size );
362 vertexOrder[ vx.second ].first = count++ ;
365 for(
int j=0; j<size; ++j )
367 if( vertexOrder[ j ].first > 0.0 && vertexOrder[ j ].second.first >= 0 )
369 if( vertexOrder[ j ].first <= vertexOrder[ vertexOrder[ j ].second.first ].first )
374 std::cout <<
"Rescale done, time = " << restimer.elapsed() << std::endl;
382 std::map< double, int > vx;
383 for(
int j=0; j<4; ++j )
385 assert( vertexOrder[ neigh[ j ] ].first > 0 );
386 vx.insert( std::make_pair( vertexOrder[ neigh[ j ] ].first, j ) );
390 for(
const auto& vxItem : vx )
392 newNeigh[ count++ ] = neigh[ vxItem.second ] ;
397 unsigned int type = 0;
399 for(
unsigned int i =0; i < 4; ++i)
405 if( ( contained[0] && contained[1] && contained[2] && contained[3] ) || ( !contained[0] && !contained[1] && !contained[2] && !contained[3] ) )
414 for(
unsigned int i = 0 ; i < 4; ++i)
418 auto it = std::find(V1Part.begin(),V1Part.end(),newNeigh[i]);
423 auto it = std::find(V0Part.begin(),V0Part.end(),newNeigh[i]);
428 for(
unsigned int i = 0; i < 4; ++i)
431 newNeigh[ i ] = V0Part[ i ];
433 newNeigh[ i ] = V1Part[ i - 1 ] ;
435 newNeigh[ i ] = V0Part[ i - type];
438 types_[neighIndex] = type % 3;
443 for(
unsigned int i =0 ; i < 3; ++i)
445 if( newNeigh[i] != neigh[i] )
447 auto neighIt = std::find(neigh.begin() + i,neigh.end(),newNeigh[i]);
448 std::swap(*neighIt,neigh[i]);
449 neighOrientation = ! neighOrientation;
454 for(
unsigned int i = 0; i < 4 ; ++i)
456 getFace(neigh, i, faceElement);
458 if(faceElement.second[0] == faceElement.second[1])
463 activeFaceList.push_back(faceElement);
465 activeFaceList.push_front(faceElement);
469 const int onePercent = numberOfElements / 100 ;
470 if( onePercent > 0 && counter % onePercent == 0 )
472 std::cout <<
"Done: element " << counter <<
" of " << numberOfElements <<
" time used = " << timer.elapsed() << std::endl;
496 std::vector<int> nodePriority;
500 const unsigned int numberOfElements =
elements_.size();
502 for(
unsigned int elIndex =0 ; elIndex < numberOfElements; ++elIndex)
506 int priorityNode = -1;
509 for(
int i = 0; i < 4; ++i)
511 int tmpPrio = nodePriority[el[i]];
515 if( priorityNode < 0 || tmpPrio > nodePriority[el[priorityNode]] )
518 getFace(el,3 - i,face );
520 if(freeFaces.find(face) != freeFaces.end())
523 if(priorityNode > -1)
525 fixNode(el, priorityNode);
527 else if(freeNode > -1)
529 nodePriority[el[freeNode]] = currNodePriority;
530 fixNode(el, freeNode);
535 nodePriority[el[
type1node_]] = currNodePriority;
546 const auto freeFacesEnd = freeFaces.end();
547 const auto freeFaceIt = freeFaces.find(face);
548 if( freeFaceIt != freeFacesEnd)
550 while(!checkFaceCompatibility(faceElement))
554 freeFaceIt->second[1] = elIndex;
556 else if(activeFaces.find(face) != activeFaces.end())
558 while(!checkFaceCompatibility(faceElement))
562 activeFaces.erase(face);
566 freeFaces.insert({{face,{elIndex,elIndex}}});
570 for(
auto&& i : {0,1,2,3})
575 getFace(el,i,faceElement);
577 unsigned int neighborIndex = faceElement.second[0] == elIndex ? faceElement.second[1] : faceElement.second[0];
578 if(freeFaces.find(face) != freeFacesEnd)
580 while(!checkFaceCompatibility(faceElement))
585 else if(activeFaces.find(face) != activeFaces.end())
587 if(!checkFaceCompatibility(faceElement))
589 checkFaceCompatibility(faceElement,
true) ;
592 activeFaces.erase(face);
596 activeFaces.insert({{face,{elIndex,elIndex}}});
616 std::vector<bool>& elementOrientation,
617 std::vector<int>& types,
618 const bool stevenson =
false )
632 for(
int el = 0; el<nElements; ++el )
651 swapRefinementType();
659 swapRefinementType();
664 void swapRefinementType()
667 for(
int el = 0; el<nElements; ++el )
682 void applyStandardOrientation ()
689 std::swap(element[2],element[3]);
698 bool checkFaceCompatibility(std::pair<FaceType, EdgeType> face,
bool verbose =
false)
701 int elIndex = face.second[0];
702 int neighIndex = face.second[1];
703 if(elIndex != neighIndex)
705 getRefinementEdge(
elements_[elIndex], face.first, edge1,
types_[elIndex]);
706 getRefinementEdge(
elements_[neighIndex], face.first, edge2,
types_[neighIndex]);
725 bool checkStrongCompatibility(std::pair<FaceType, EdgeType> face,
bool verbose =
false)
727 int elIndex = face.second[0];
728 int neighIndex = face.second[1];
731 if(elIndex != neighIndex)
733 int elType =
types_[elIndex];
734 int neighType =
types_[neighIndex];
735 unsigned int elFaceIndex = getFaceIndex(
elements_[elIndex], face.first);
736 unsigned int elNodeIndex = 3 - elFaceIndex;
737 unsigned int neighFaceIndex = getFaceIndex(
elements_[neighIndex], face.first);
738 unsigned int neighNodeIndex = 3 - neighFaceIndex;
739 if(elType == neighType)
745 if( elNodeIndex == neighNodeIndex ||
752 if( elType % 3 == (neighType +1) %3 )
754 std::swap(elType, neighType);
755 std::swap(elNodeIndex, neighNodeIndex);
756 std::swap(elFaceIndex, neighFaceIndex);
757 std::swap(elIndex, neighIndex);
761 assert(neighType == 1);
766 assert(neighType == 2);
771 assert(neighType == 0);
776 std::cerr <<
"No other types than 0,1,2 implemented. Aborting" << std::endl;
785 if (verbose && result ==
false )
823 face = {el[1],el[2],el[3]};
826 face = {el[0],el[2],el[3]};
829 face = {el[0],el[1],el[3]};
832 face = {el[0],el[1],el[2]};
835 std::cerr <<
"index " << faceIndex <<
" NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
838 std::sort(face.begin(), face.end());
845 getFace(el, faceIndex, face);
861 edge = {el[0],el[2]};
864 edge = {el[0],el[3]};
867 edge = {el[0],el[3]};
870 edge = {el[1],el[3]};
873 std::cerr <<
"index " << faceIndex <<
" NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
882 edge = {el[0],el[1]};
885 edge = {el[0],el[1]};
888 edge = {el[0],el[2]};
891 edge = {el[1],el[3]};
894 std::cerr <<
"index " << faceIndex <<
" NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
899 else if(type == 1 || type == 2)
906 edge = {el[0],el[2]};
909 edge = {el[0],el[3]};
912 edge = {el[0],el[3]};
915 edge = {el[2],el[3]};
918 std::cerr <<
"index " << faceIndex <<
" NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
927 edge = {el[0],el[1]};
930 edge = {el[0],el[1]};
933 edge = {el[0],el[3]};
936 edge = {el[1],el[3]};
939 std::cerr <<
"index " << faceIndex <<
" NOT IMPLEMENTED FOR TETRAHEDRONS" << std::endl;
946 std::cerr <<
"no other types than 0, 1, 2 implemented." << std::endl;
949 std::sort(edge.begin(),edge.end());
956 getRefinementEdge(el, getFaceIndex(el, face), edge, type);
964 for(
int i =0; i<4 ; ++i)
966 if(!( el[i] == face[0] || el[i] == face[1] || el[i] == face[2] ) )
977 void buildNeighbors()
985 unsigned int index = 0;
989 for(
int i = 0; i< 4; ++i)
991 getFace(el, i, face);
993 if(faceInList == nend)
995 indexPair = {index, index};
996 neighbours_.insert(std::make_pair (face, indexPair ) );
1001 faceInList->second[1] = index;
1016 void calculateV0(
int variante,
int threshold)
1023 std::vector<int> numberOfAdjacentRefEdges(
nVertices_, 0);
1028 numberOfAdjacentRefEdges [ el [
type0nodes_[ 0 ] ] ] ++;
1029 numberOfAdjacentRefEdges [ el [
type0nodes_[ 1 ] ] ] ++;
1033 if(numberOfAdjacentRefEdges[ i ] >= threshold )
1042 std::vector<int> numberOfAdjacentElements(
nVertices_, 0);
1043 std::vector<bool> vertexOnBoundary(
nVertices_,
false);
1047 if(neigh.second[1] == neigh.second[0])
1049 for(
unsigned int i =0; i < 3 ; ++i)
1051 vertexOnBoundary[ neigh.first [ i ] ] =
true;
1058 for(
unsigned int i =0; i < 4; ++i)
1060 numberOfAdjacentElements[ el [ i ] ] ++;
1065 double bound = vertexOnBoundary[ i ] ? threshold / 2. : threshold ;
1066 if(numberOfAdjacentElements[ i ] < bound )
1076 std::default_random_engine generator;
1077 std::uniform_int_distribution<int> distribution(1,6);
1080 int roll = distribution(generator);
1104 std::cout <<
"#V0 #V1" << std::endl << sizeOfV0 <<
" " << sizeOfV1 << std::endl;
Definition: bisectioncompatibility.hh:16
static int & variant()
Definition: bisectioncompatibility.hh:17
static int & threshold()
Definition: bisectioncompatibility.hh:23
static int & useAnnouncedEdge()
Definition: bisectioncompatibility.hh:29
Definition: bisectioncompatibility.hh:41
bool type1Algorithm()
Definition: bisectioncompatibility.hh:483
void printElement(int index)
Definition: bisectioncompatibility.hh:171
BisectionCompatibility(const VertexVectorType &vertices, const std::vector< ElementType > &elements, const bool stevenson)
Definition: bisectioncompatibility.hh:92
std::vector< bool > elementOrientation_
Definition: bisectioncompatibility.hh:58
void printNB()
Definition: bisectioncompatibility.hh:162
std::vector< int > types_
Definition: bisectioncompatibility.hh:67
unsigned int type1face_
Definition: bisectioncompatibility.hh:79
FaceMapType neighbours_
Definition: bisectioncompatibility.hh:60
std::array< unsigned int, 3 > FaceType
Definition: bisectioncompatibility.hh:47
const int variant_
Definition: bisectioncompatibility.hh:85
int stronglyCompatibleFaces()
Definition: bisectioncompatibility.hh:114
const VertexVectorType & vertices_
Definition: bisectioncompatibility.hh:54
const size_t nVertices_
Definition: bisectioncompatibility.hh:62
const int threshold_
Definition: bisectioncompatibility.hh:86
unsigned int type1node_
Definition: bisectioncompatibility.hh:77
VertexVector VertexVectorType
Definition: bisectioncompatibility.hh:45
std::vector< unsigned int > ElementType
Definition: bisectioncompatibility.hh:48
const bool useAnnouncedEdge_
Definition: bisectioncompatibility.hh:87
EdgeType type0nodes_
Definition: bisectioncompatibility.hh:73
std::array< unsigned int, 2 > EdgeType
Definition: bisectioncompatibility.hh:49
void alberta2Stevenson()
Definition: bisectioncompatibility.hh:655
std::pair< FaceType, EdgeType > FaceElementType
Definition: bisectioncompatibility.hh:51
void returnElements(std::vector< ElementType > &elements, std::vector< bool > &elementOrientation, std::vector< int > &types, const bool stevenson=false)
Definition: bisectioncompatibility.hh:615
std::vector< ElementType > elements_
Definition: bisectioncompatibility.hh:57
void stevenson2Alberta()
Definition: bisectioncompatibility.hh:647
std::map< FaceType, EdgeType > FaceMapType
Definition: bisectioncompatibility.hh:50
bool make6CompatibilityCheck()
Definition: bisectioncompatibility.hh:151
EdgeType type0faces_
Definition: bisectioncompatibility.hh:75
std::vector< bool > containedInV0_
Definition: bisectioncompatibility.hh:65
bool stevensonRefinement_
Definition: bisectioncompatibility.hh:70
bool type0Algorithm()
Definition: bisectioncompatibility.hh:195
bool compatibilityCheck()
Definition: bisectioncompatibility.hh:142