diff --git a/misc/examples/sizes.cc b/misc/examples/sizes.cc
index 04b32ddb370f1c91db69e98dbccde6922ec8a1f5..32c7159d1363cddeb2e9149c17694f46caa1f01e 100644
--- a/misc/examples/sizes.cc
+++ b/misc/examples/sizes.cc
@@ -83,36 +83,46 @@ void checkRefinements( GitterType& grid )
// refine grid globally, i.e. mark all elements and then call adapt
template <class GitterType>
-void globalRefine(GitterType& grid, bool global)
+void globalRefine(GitterType& grid, bool global,int step)
{
{
if (global)
{
- // get LeafIterator which iterates over all leaf elements of the grid
- LeafIterator < Gitter::helement_STI > w (grid) ;
+ // get LeafIterator which iterates over all leaf elements of the grid
+ LeafIterator < Gitter::helement_STI > w (grid) ;
- for (w->first () ; ! w->done () ; w->next ())
- {
- // mark element for refinement
- // w->item ().tagForGlobalRefinement ();
-
- typedef typename GitterType :: Objects :: tetra_IMPL tetra_IMPL ;
- // mark element for refinement
- tetra_IMPL* item = ((tetra_IMPL *) &w->item ());
-
- //item->checkTetra( item, item->nChild(), true );
-
- typedef Gitter ::Geometric :: TetraRule TetraRule ;
- item->request ( TetraRule :: bisect );
-
+ for (w->first () ; ! w->done () ; w->next ())
+ {
+ if( w->item ().type() == tetra )
+ {
+ typedef typename GitterType :: Objects :: tetra_IMPL tetra_IMPL ;
+ // mark element for refinement
+ tetra_IMPL* item = ((tetra_IMPL *) &w->item ());
+
+ //item->checkTetra( item, item->nChild(), true );
+
+ typedef Gitter ::Geometric :: TetraRule TetraRule ;
+ item->request ( TetraRule :: bisect );
+ }
+ else
+ {
+ // mark element for refinement
+ w->item ().tagForGlobalRefinement ();
+ }
// break ;
- }
+ }
// adapt grid
grid.adapt ();
+ grid.printsize () ;
}
else
{
+ double t = double(step)/10.;
double center[3] = {0.2,0.2,0.2};
+ double dir[3] = {1.0,0.0,0.0};
+ center[0] += dir[0]*t;
+ center[1] += dir[1]*t;
+ center[2] += dir[2]*t;
double rad=0.6;
grid.refineBall(center,rad,10);
}
@@ -171,16 +181,19 @@ int main (int argc, char ** argv, const char ** envp)
{
MPI_Init(&argc,&argv);
- int mxl = 0;
+ int mxl = 0, glb = 0;
const char* filename = 0 ;
if (argc < 2)
{
filename = "../macrogrids/reference.tetra";
mxl = 1;
- cout << "usage: "<< argv[0] << " <macro grid> <opt: level> \n";
+ glb = 1;
+ cout << "usage: "<< argv[0] << " <macro grid> <opt: level global> \n";
}
else
+ {
filename = argv[ 1 ];
+ }
{
int rank = 0;
@@ -196,6 +209,13 @@ int main (int argc, char ** argv, const char ** envp)
}
else
mxl = atoi(argv[2]);
+ if (argc < 4)
+ {
+ if( rank == 0 )
+ cout << "Default global refinement = "<< glb << " choosen! \n";
+ }
+ else
+ glb = atoi(argv[3]);
std::string macroname( filename );
@@ -218,84 +238,31 @@ int main (int argc, char ** argv, const char ** envp)
//cout << "P[ " << rank << " ] : Grid generated! \n";
grid.printsize();
cout << "---------------------------------------------\n";
- /*
+
grid.printMemUsage();
- for (int i = 0; i < 6; ++i)
- globalRefine(grid, true);
- */
-
- //int bla;
- // cin >> bla;
- for( int i = 0; i <= mxl; ++i )
+ for (int i = 0; i < glb; ++i)
+ globalRefine(grid, true,-1);
+ for (int i = 0; i < glb; ++i)
+ globalRefine(grid, false,0);
+ for( int i = 0; i < mxl; ++i )
{
std::ostringstream ss;
ss << "out-" << ZeroPadNumber(i) << ".vtu";
grid.tovtk( ss.str().c_str() );
- if( i < mxl )
- globalRefine(grid, false);
+ globalRefine(grid, false,i);
}
-
-#if 0
- checkRefinements( grid );
-
- std::ofstream file( "file.out" );
- grid.duneBackup( file );
- file.close();
-
- /*
{
- ObjectStream os ;
- grid.duneBackup( os );
-
- char* buffer = ObjectStream ::allocateBuffer( os.size() );
- os.read( buffer, os.size() );
-
- std::ofstream obj( "obj.out" );
- const size_t size = os.size();
- obj.write( (const char *) &size, sizeof( size_t ) );
- obj.write( buffer, size );
- ObjectStream :: freeBuffer( buffer );
+ std::ostringstream ss;
+ ss << "out-" << ZeroPadNumber(mxl) << ".vtu";
+ grid.tovtk( ss.str().c_str() );
}
-
- globalRefine(grid, mxl);
-
+ globalCoarsening(grid,3*glb);
{
- size_t size = 0;
-
- std::ifstream obj( "obj.out" );
- obj.read( (char *) &size, sizeof( size_t ) );
- char * buffer = ObjectStream ::allocateBuffer( size );
- obj.read( buffer, size );
- //ObjectStream :: freeBuffer( buffer );
-
- ObjectStream is;
- is.clear();
- is.write( buffer, size );
-
- {
- ObjectStream copy( is );
- std::ofstream obj( "check.out" );
- const size_t size = copy.size();
- obj.write( buffer, size );
- //ObjectStream :: freeBuffer( buffer );
- }
- }*/
-
- {
- std::ifstream file( "file.out" );
-#ifdef PARALLEL
- MpAccessMPI a (MPI_COMM_WORLD);
- GitterDunePll grid2( file, a);
-#else
- GitterDuneImpl grid2( file );
-#endif
- grid2.printsize();
+ std::ostringstream ss;
+ ss << "out-" << ZeroPadNumber(mxl+1) << ".vtu";
+ grid.tovtk( ss.str().c_str() );
}
- //levelwalk(grid, mxl);
- // globalCoarsening(grid, mxl);
- //grid.printMemUsage();
- //cin.get();
-#endif
+ return 1;
}
}
diff --git a/src/duneinterface/gitter_dune_pll_impl.cc b/src/duneinterface/gitter_dune_pll_impl.cc
index 05876facfdebfceffaabde48ec35b785a5fe8ca8..3d0a8e5d2ecd4c8b038a1f84d3a7190964b1e9fc 100644
--- a/src/duneinterface/gitter_dune_pll_impl.cc
+++ b/src/duneinterface/gitter_dune_pll_impl.cc
@@ -1395,8 +1395,6 @@ void GitterDunePll :: duneRestore(const char *fileName)
}
void GitterDunePll :: tovtk( const std::string &fn )
{
- const bool showbnd = false;
- const bool showface = true;
const int myrank = mpAccess ().myrank () ;
const int nProc = mpAccess ().psize () ;
@@ -1405,243 +1403,7 @@ void GitterDunePll :: tovtk( const std::string &fn )
// openfile
std::ofstream vtuFile;
- vtuFile.open( ss.str().c_str() );
-
- // header info
- vtuFile << "<?xml version=\"1.0\"?>" << std::endl;
- vtuFile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
- vtuFile << " <UnstructuredGrid>" << std::endl;
-
- // vertex list
- typedef std::vector< double > Vertex;
- typedef std::map< int, std::pair<int,Vertex> > VertexList;
- VertexList vertexList;
-
- int nCells = 0;
- int nBnd = 0;
- int nFaces = 0;
-#if 1
- typedef LeafIterator < Gitter::helement_STI > Iterator;
- Iterator w (*this) ;
- typedef LeafIterator < Gitter::hbndseg_STI > BndIterator;
- BndIterator wbnd (*this) ;
- typedef LeafIterator < Gitter::hface_STI > FaceIterator;
- FaceIterator wface (*this) ;
-#else
- typedef LevelIterator < Gitter::helement_STI > Iterator;
- Iterator w (*this,0) ;
- typedef LevelIterator < Gitter::hbndseg_STI > BndIterator;
- BndIterator wbnd (*this,0) ;
- typedef LevelIterator < Gitter::hface_STI > FaceIterator;
- FaceIterator wface (*this,0) ;
-#endif
- // loop to find vertexList and count cells
- {
- typedef Objects :: tetra_IMPL tetra_IMPL ;
- for (w->first () ; ! w->done () ; w->next ())
- {
- tetra_IMPL* item = ((tetra_IMPL *) &w->item ());
- for (int i=0;i<4;++i)
- {
- Vertex v(3);
- for (int k=0;k<3;++k)
- v[k] = item->myvertex(i)->Point()[k];
- vertexList[ item->myvertex(i)->getIndex() ] = make_pair(-1,v);
- }
- ++nCells;
- }
- if (showbnd)
- {
- for (wbnd->first () ; ! wbnd->done () ; wbnd->next ())
- ++nBnd;
- }
- if (showface)
- {
- for (wface->first () ; ! wface->done () ; wface->next ())
- ++nFaces;
- }
- }
-
- vtuFile << " <Piece NumberOfPoints=\"" << vertexList.size() << "\" "
- << "NumberOfCells=\"" << nCells+nBnd+nFaces << "\">" << std::endl;
-
- // cell data
- {
- vtuFile << " <CellData Scalars=\"cell-id\">" << std::endl;
- vtuFile << " <DataArray type=\"Float32\" Name=\"cell-id\" NumberOfComponents=\"1\">" << std::endl;
- vtuFile << " ";
-
- typedef Objects :: tetra_IMPL tetra_IMPL ;
- for (w->first () ; ! w->done () ; w->next ())
- {
- tetra_IMPL* item = ((tetra_IMPL *) &w->item ());
- // vtuFile << item->getIndex() << " ";
- bool ok = true;
- for (int k=0;k<4;++k)
- ok &= item->myneighbour( k ).first->isRealObject();
- if (!ok)
- {
- std::cout << "Problem: " << item << std::endl;
- for (int k=0;k<4;++k)
- if (!item->myneighbour( k ).first->isRealObject())
- {
- std::cout << item->myhface3(k) << std::endl;
- if ( item->myhface3(k)->nb.front().first->isRealObject() )
- {
- std::cout << item->myhface3(k)->nb.front().first << std::endl;
- std::cout << ((hbndseg3_GEO *) item->myhface3(k)->nb.front().first)->myhface3(0) << std::endl;
- }
- if ( item->myhface3(k)->nb.rear().first->isRealObject() )
- {
- std::cout << item->myhface3(k)->nb.rear().first << std::endl;
- std::cout << ((hbndseg3_GEO *) item->myhface3(k)->nb.rear().first)->myhface3(0) << std::endl;
- }
- std::cout << std::endl;
- }
- std::cout << std::endl;
- }
-
- vtuFile << ((ok)?1:-1) << " ";
- }
-
- vtuFile << std::endl;
- if (showbnd)
- for (wbnd->first () ; ! wbnd->done () ; wbnd->next ())
- {
- hbndseg3_GEO* item = ((hbndseg3_GEO *) &wbnd->item ());
- bool ok = true;
- ok &= item->myhface3(0)->nb.front().first->isRealObject();
- ok &= item->myhface3(0)->nb.rear().first->isRealObject();
- if (!ok)
- {
- if ( item->myhface3(0)->nb.front().first->isRealObject() )
- assert( item->myhface3(0)->nb.front().first == item );
- if ( item->myhface3(0)->nb.rear().first->isRealObject() )
- assert( item->myhface3(0)->nb.rear().first == item );
- std::cout << "Problem: " << item << std::endl;
- std::cout << item->myhface3(0) << std::endl;
- std::cout << item->myhface3(0)->nb.front().first << std::endl;
- std::cout << item->myhface3(0)->nb.rear().first << std::endl;
- std::cout << std::endl;
- }
- vtuFile << ((ok)?1:-1)*item->myhface3(0)->ref << " ";
- }
- if (showface)
- for (wface->first () ; ! wface->done () ; wface->next ())
- {
- hface3_GEO* item = ((hface3_GEO *) &wface->item ());
- bool ok = true;
- ok &= item->nb.front().first->isRealObject();
- ok &= item->nb.rear().first->isRealObject();
- vtuFile << ((ok)?1:-1)*item->ref << " ";
- }
-
- vtuFile << std::endl;
- vtuFile << " </DataArray>" << std::endl;
- vtuFile << " </CellData>" << std::endl;
- }
-
- // points info
- {
- vtuFile << " <Points>" << std::endl;
- vtuFile << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
-
- const VertexList::iterator end = vertexList.end();
- int nr=0;
- for( VertexList::iterator i = vertexList.begin(); i != end; ++i, ++nr )
- {
- vtuFile << " " << (*i).second.second[ 0 ] << " " << (*i).second.second[ 1 ] << " " << (*i).second.second[ 2 ] << std::endl;
- (*i).second.first = nr;
- }
-
- vtuFile << " </DataArray>" << std::endl;
- vtuFile << " </Points>" << std::endl;
- }
-
- // cell info
- {
- vtuFile << " <Cells>" << std::endl;
- // connectivity
- vtuFile << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">" << std::endl;
- vtuFile << " ";
-
- typedef Objects :: tetra_IMPL tetra_IMPL ;
- for (w->first () ; ! w->done () ; w->next ())
- {
- tetra_IMPL* item = ((tetra_IMPL *) &w->item ());
- for (int i=0;i<4;++i)
- {
- vtuFile << " " << vertexList[item->myvertex(i)->getIndex()].first;
- }
- }
- if (showbnd)
- for (wbnd->first () ; ! wbnd->done () ; wbnd->next ())
- {
- hbndseg3_GEO* item = ((hbndseg3_GEO *) &wbnd->item ());
- for (int i=0;i<3;++i)
- {
- vtuFile << " " << vertexList[item->myvertex(0,i)->getIndex()].first;
- }
- }
- if (showface)
- for (wface->first () ; ! wface->done () ; wface->next ())
- {
- hface3_GEO* item = ((hface3_GEO *) &wface->item ());
- for (int i=0;i<3;++i)
- {
- vtuFile << " " << vertexList[item->myvertex(i)->getIndex()].first;
- }
- }
- vtuFile << std::endl;
- vtuFile << " </DataArray>" << std::endl;
-
- // offsets
- vtuFile << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">" << std::endl;
- vtuFile << " ";
-
- for( int i = 0; i < nCells; ++i )
- {
- vtuFile << " " << (i+1)*4;
- }
- for( int i = 0; i < nBnd; ++i )
- {
- vtuFile << " " << nCells*4 + (i+1)*3;
- }
- for( int i = 0; i < nFaces; ++i )
- {
- vtuFile << " " << nCells*4 + nBnd*3 + (i+1)*3;
- }
- vtuFile << std::endl;
-
- vtuFile << " </DataArray>" << std::endl;
-
- // cell type
- vtuFile << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
- vtuFile << " ";
-
- for( int i = 0; i < nCells; ++i )
- {
- vtuFile << " " << 10; // 10 for tetrahedra
- }
- for( int i = 0; i < nBnd; ++i )
- {
- vtuFile << " " << 5; // 5 for triangle
- }
- for( int i = 0; i < nFaces; ++i )
- {
- vtuFile << " " << 5; // 5 for trianglea
- }
- vtuFile << std::endl;
-
- vtuFile << " </DataArray>" << std::endl;
- }
-
- vtuFile << " </Cells>" << std::endl;
- vtuFile << " </Piece>" << std::endl;
- vtuFile << " </UnstructuredGrid>" << std::endl;
- vtuFile << "</VTKFile>" << std::endl;
-
- vtuFile.close();
+ Gitter :: tovtk( ss.str() );
if( myrank == 0 )
{
diff --git a/src/serial/gitter_geo.cc b/src/serial/gitter_geo.cc
index 5cf293f7f3958ca6aea895d65dd0e943a879186d..7127383ab27ebb6146707a1009c31c5394d3044b 100644
--- a/src/serial/gitter_geo.cc
+++ b/src/serial/gitter_geo.cc
@@ -549,9 +549,9 @@ int Gitter :: Geometric :: Tetra :: tagForBallRefinement (const alucoord_t (&cen
}
}
#endif
- return hit ? (request (myrule_t :: bisect), 1) : (request (myrule_t :: nosplit), 0);
- // (level () < limit ? (request (myrule_t :: iso8), 1)
- // : (request (myrule_t :: nosplit), 0)) : (request (myrule_t :: crs), 1) ;
+ if (!hit) return (request (myrule_t :: crs), 1) ;
+ if (level()>limit) return (request (myrule_t :: nosplit), 0);
+ return (request (myrule_t :: bisect), 1);
}
// ###### #####
diff --git a/src/serial/gitter_sti.cc b/src/serial/gitter_sti.cc
index 732416b7a9938c0fe9d89cf8ee31f27bb2e9bb92..5f0fd22179ddb81b16a18b79ca37bfdf97aeaa56 100644
--- a/src/serial/gitter_sti.cc
+++ b/src/serial/gitter_sti.cc
@@ -312,13 +312,22 @@ bool Gitter :: markNonConform()
{
bool x = true ;
leaf_element__macro_element__iterator i (container ()) ;
- // std::cout << "check non conform refinement" << std::endl;
for( i.first(); ! i.done() ; i.next()) { x &= i.item ().markNonConform () ; }
return x;
}
+void Gitter :: markEdgeCoarsening ()
+{
+ edgeCoarseningFlags_.assign( indexManagerStorage().get(2).getMaxIndex(), true );
+ leaf_element__macro_element__iterator i (container ()) ;
+ for( i.first(); ! i.done() ; i.next() )
+ {
+ i.item().markEdgeCoarsening();
+ }
+}
void Gitter :: coarse()
{
+ markEdgeCoarsening();
assert (debugOption (20) ? (cout << "**INFO Gitter :: coarse ()" << endl, 1) : 1) ;
{
AccessIterator < helement_STI > :: Handle i (container ()) ;
@@ -327,10 +336,23 @@ void Gitter :: coarse()
i.item ().coarse () ;
}
}
+ std::ostringstream ss;
+ int filenr = adaptstep*100+nr;
+ ss << "crs-" << ZeroPadNumber(filenr) << ".vtu";
+ tovtk( ss.str() );
+ ++nr;
}
-void Gitter :: tovtk( const std::string &fn )
+template <class element_t, class bndseg>
+void Gitter :: tovtkImpl( const std::string &fn,
+ const int elementVertices,
+ const element_t*, const bndseg* )
{
+ const bool showbnd = false;
+ const bool showface = true;
+
+ const int nFaceVertices = ( elementVertices == 4 ) ? 3 : 4 ;
+
// openfile
std::ofstream vtuFile;
vtuFile.open( fn.c_str() );
@@ -345,31 +367,41 @@ void Gitter :: tovtk( const std::string &fn )
typedef std::map< int, std::pair<int,Vertex> > VertexList;
VertexList vertexList;
- int nCells = 0;
-
- typedef Gitter :: Geometric :: tetra_GEO tetra_GEO ;
typedef LeafIterator < Gitter::helement_STI > Iterator;
- // typedef LevelIterator < Gitter::helement_STI > Iterator;
Iterator w (*this) ;
+ typedef LeafIterator < Gitter::hbndseg_STI > BndIterator;
+ BndIterator wbnd (*this) ;
+ typedef LeafIterator < Gitter::hface_STI > FaceIterator;
+ FaceIterator wface (*this) ;
+
+ const int nCells = w->size();
+ const int nBnd = showbnd ? wbnd->size() : 0;
+ const int nFaces = showface ? wface->size() : 0;
+
+ typedef typename element_t :: myvertex_t myvertex_t ;
+ typedef typename element_t :: myhface_t myhface_t;
// loop to find vertexList and count cells
{
for (w->first () ; ! w->done () ; w->next ())
{
- tetra_GEO* item = ((tetra_GEO *) &w->item ());
- for (int i=0;i<4;++i)
+ element_t* item = ((element_t *) &w->item ());
+ // store vertices
+ for (int i=0; i < elementVertices; ++i )
{
Vertex v(3);
- for (int k=0;k<3;++k)
- v[k] = item->myvertex(i)->Point()[k];
- vertexList[ item->myvertex(i)->getIndex() ] = make_pair(-1,v);
+ const myvertex_t* vx = item->myvertex( i );
+ assert( vx );
+ const alucoord_t (&coord)[ 3 ] = vx->Point();
+ // copy coordinates
+ for (int k=0;k<3;++k) v[k] = coord[ k ];
+ vertexList[ vx->getIndex() ] = make_pair(-1,v);
}
- ++nCells;
}
}
vtuFile << " <Piece NumberOfPoints=\"" << vertexList.size() << "\" "
- << "NumberOfCells=\"" << nCells << "\">" << std::endl;
+ << "NumberOfCells=\"" << nCells+nBnd+nFaces << "\">" << std::endl;
// cell data
{
@@ -379,7 +411,76 @@ void Gitter :: tovtk( const std::string &fn )
for (w->first () ; ! w->done () ; w->next ())
{
- vtuFile << w->item ().getIndex() << " ";
+ element_t* item = ((element_t *) &w->item ());
+ // vtuFile << item->getIndex() << " ";
+ bool ok = true;
+ const int nFaces = item->nFaces();
+ for (int k=0; k < nFaces; ++k )
+ ok &= item->myneighbour( k ).first->isRealObject();
+
+ if (!ok)
+ {
+ std::cout << "Problem: " << item << std::endl;
+ for (int k=0; k<nFaces; ++k)
+ {
+ if (!item->myneighbour( k ).first->isRealObject())
+ {
+ std::cout << item->myhface(k) << std::endl;
+ if ( item->myhface(k)->nb.front().first->isRealObject() )
+ {
+ std::cout << item->myhface(k)->nb.front().first << std::endl;
+ std::cout << ((bndseg *) item->myhface(k)->nb.front().first)->myhface(0) << std::endl;
+ }
+ if ( item->myhface(k)->nb.rear().first->isRealObject() )
+ {
+ std::cout << item->myhface(k)->nb.rear().first << std::endl;
+ std::cout << ((bndseg *) item->myhface(k)->nb.rear().first)->myhface(0) << std::endl;
+ }
+ std::cout << std::endl;
+ }
+ }
+ std::cout << std::endl;
+ }
+
+ vtuFile << ((ok)?1:-1) << " ";
+ }
+
+ vtuFile << std::endl;
+ if (showbnd)
+ {
+ for (wbnd->first () ; ! wbnd->done () ; wbnd->next ())
+ {
+ bndseg* item = ((bndseg *) &wbnd->item ());
+ bool ok = true;
+ ok &= item->myhface(0)->nb.front().first->isRealObject();
+ ok &= item->myhface(0)->nb.rear().first->isRealObject();
+ if (!ok)
+ {
+ if ( item->myhface(0)->nb.front().first->isRealObject() )
+ assert( item->myhface(0)->nb.front().first == item );
+ if ( item->myhface(0)->nb.rear().first->isRealObject() )
+ assert( item->myhface(0)->nb.rear().first == item );
+ std::cout << "Problem: " << item << std::endl;
+ std::cout << item->myhface(0) << std::endl;
+ std::cout << item->myhface(0)->nb.front().first << std::endl;
+ std::cout << item->myhface(0)->nb.rear().first << std::endl;
+ std::cout << std::endl;
+ }
+ vtuFile << ((ok)?1:-1)*item->myhface(0)->ref << " ";
+ }
+ }
+
+ if (showface)
+ {
+ for (wface->first () ; ! wface->done () ; wface->next ())
+ {
+ myhface_t* item = ((myhface_t *) &wface->item ());
+ bool ok = true;
+ ok &= item->nb.front().first->isRealObject();
+ ok &= item->nb.rear().first->isRealObject();
+ // assert(item->ref>0);
+ vtuFile << ((ok)?1:-1)*item->ref << " ";
+ }
}
vtuFile << std::endl;
@@ -413,10 +514,34 @@ void Gitter :: tovtk( const std::string &fn )
for (w->first () ; ! w->done () ; w->next ())
{
- tetra_GEO* item = ((tetra_GEO *) &w->item ());
- for (int i=0;i<4;++i)
+ element_t* item = ((element_t *) &w->item ());
+ for (int i=0; i<elementVertices; ++i)
+ {
+ vtuFile << " " << vertexList[ item->myvertex(i)->getIndex() ].first;
+ }
+ }
+
+ if (showbnd)
+ {
+ for (wbnd->first () ; ! wbnd->done () ; wbnd->next ())
+ {
+ bndseg* item = ((bndseg *) &wbnd->item ());
+ for (int i=0; i<nFaceVertices; ++i)
+ {
+ vtuFile << " " << vertexList[ item->myvertex(0,i)->getIndex() ].first;
+ }
+ }
+ }
+
+ if (showface)
+ {
+ for (wface->first () ; ! wface->done () ; wface->next ())
{
- vtuFile << " " << vertexList[item->myvertex(i)->getIndex()].first;
+ myhface_t* item = ((myhface_t *) &wface->item ());
+ for (int i=0; i<nFaceVertices; ++i)
+ {
+ vtuFile << " " << vertexList[ item->myvertex(i)->getIndex() ].first;
+ }
}
}
vtuFile << std::endl;
@@ -428,7 +553,15 @@ void Gitter :: tovtk( const std::string &fn )
for( int i = 0; i < nCells; ++i )
{
- vtuFile << " " << (i+1)*4;
+ vtuFile << " " << (i+1)* elementVertices;
+ }
+ for( int i = 0; i < nBnd; ++i )
+ {
+ vtuFile << " " << nCells* elementVertices + (i+1)* nFaceVertices;
+ }
+ for( int i = 0; i < nFaces; ++i )
+ {
+ vtuFile << " " << nCells*elementVertices + nBnd*nFaceVertices + (i+1)*nFaceVertices;
}
vtuFile << std::endl;
@@ -438,15 +571,26 @@ void Gitter :: tovtk( const std::string &fn )
vtuFile << " <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">" << std::endl;
vtuFile << " ";
+ // 10 for tetrahedra, 12 for hexahedron
+ const int elemId = ( elementVertices == 4 ) ? 10 : 12 ;
for( int i = 0; i < nCells; ++i )
{
- vtuFile << " " << 10; // 10 for tetrahedra
+ vtuFile << " " << elemId ;
+ }
+ // 5 for triangle, 9 for quadrilateral
+ const int faceId = ( nFaceVertices == 3 ) ? 5 : 9 ;
+ for( int i = 0; i < nBnd; ++i )
+ {
+ vtuFile << " " << faceId;
+ }
+ for( int i = 0; i < nFaces; ++i )
+ {
+ vtuFile << " " << faceId;
}
vtuFile << std::endl;
vtuFile << " </DataArray>" << std::endl;
}
-
vtuFile << " </Cells>" << std::endl;
vtuFile << " </Piece>" << std::endl;
vtuFile << " </UnstructuredGrid>" << std::endl;
@@ -456,6 +600,21 @@ void Gitter :: tovtk( const std::string &fn )
std::cout << "data written to " << fn << std::endl;
}
+void Gitter :: tovtk( const string& filename )
+{
+ typedef LeafIterator < Gitter::helement_STI > Iterator;
+ Iterator w (*this) ;
+ w->first();
+ if( ! w->done() && w->item().type() == hexa )
+ {
+ tovtkImpl( filename, 8, (Geometric :: hexa_GEO *) 0, (Geometric :: hbndseg4_GEO * ) 0 );
+ }
+ else
+ {
+ tovtkImpl( filename, 4, (Geometric :: tetra_GEO *) 0, (Geometric :: tetra_GEO * ) 0 );
+ }
+}
+
bool Gitter :: adapt ()
{
assert (debugOption (20) ? (cout << "**INFO Gitter :: adapt ()" << endl, 1) : 1) ;
@@ -467,11 +626,9 @@ bool Gitter :: adapt ()
do {
refined &= refine ();
// check for conformity
- // if noconform break;
- std::cout << "check non conform refinement" << std::endl;
x = markNonConform();
}
- while (!x); // need something here on required conformity
+ while (!x);
if (!refined) {
cerr << "**WARNUNG (IGNORIERT) Verfeinerung nicht vollst\"andig (warum auch immer)\n" ;
@@ -480,6 +637,7 @@ bool Gitter :: adapt ()
}
int lap = clock () ;
coarse () ;
+ assert ( markNonConform() );
int end = clock () ;
if (debugOption (1)) {
float u1 = (float)(lap - start)/(float)(CLOCKS_PER_SEC) ;
diff --git a/src/serial/gitter_sti.h b/src/serial/gitter_sti.h
index 96a5056526dbe0efe75e9efe6502091bbeaddaab..5cbca57e7989953051018abf38d487eab2ea2bc2 100644
--- a/src/serial/gitter_sti.h
+++ b/src/serial/gitter_sti.h
@@ -573,6 +573,8 @@ public :
virtual void backupIndex (ObjectStream& os ) const { backupIndexErr(); }
virtual void restoreIndex (ObjectStream& is, RestoreInfo& ) { restoreIndexErr(); }
+
+ virtual bool canCoarsen( std::vector<bool> &edgeCoarseningFlags ) const;
} ;
class hface : public FacePllXDefault,
@@ -687,6 +689,7 @@ public :
virtual int segmentIndex (const int) const { return -1; }
virtual bool markNonConform () = 0;
+ virtual void markEdgeCoarsening () = 0;
// mark element for using iso8 rule
virtual int tagForGlobalRefinement () = 0 ;
// mark element for coarsening
@@ -1435,18 +1438,24 @@ public :
// und Prismen sorgf"altig bedacht werden.
// Der Prototyp steht in 'gitter_geo.cc'.
- typedef class Tetra : public helement_STI, public hasFace3, public MyAlloc {
- protected :
-
+ typedef class Tetra
+ : public helement_STI, public hasFace3, public MyAlloc
+ {
+ public :
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface3_GEO myhface3_t ;
+ typedef hedge1_GEO myhedge_t ;
+ typedef hface3_GEO myhface_t ;
+
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface3_t ;
+
typedef TetraRule myrule_t ;
typedef pair < hasFace3 *, int > myneighbour_t ;
- inline Tetra (myhface3_t *, int, myhface3_t *, int,
- myhface3_t *, int, myhface3_t *, int) ;
+ protected :
+ inline Tetra (myhface_t *, int, myhface_t *, int,
+ myhface_t *, int, myhface_t *, int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
public :
@@ -1463,12 +1472,16 @@ public :
static const vector<int> & facesNotOnFace( const int face ) ;
inline virtual ~Tetra () ;
- inline hface3_GEO * myhface3 (int) ;
- inline const hface3_GEO * myhface3 (int) const ;
+ inline myhface_t* myhface (int) ;
+ inline const myhface_t* myhface (int) const ;
+
+ inline myhface_t* myhface3 ( int i ) { return myhface( i ); } // don't use these methods, use the ones without number
+ inline const myhface_t* myhface3 ( int i ) const { return myhface( i ); } // don't use these methods, use the ones without number
+
inline VertexGeo * myvertex (int) ;
inline const VertexGeo * myvertex (int) const ;
- inline hedge1_GEO * myhedge1(int);
- inline const hedge1_GEO * myhedge1(int) const;
+ inline myhedge_t* myhedge1(int);
+ inline const myhedge_t* myhedge1(int) const;
inline VertexGeo * myvertex (int,int) ;
inline const VertexGeo * myvertex (int,int) const ;
inline pair < hasFace3 *, int > myneighbour (int) ;
@@ -1476,8 +1489,8 @@ public :
// Dune extension
// return pair, first = pointer to face, second = twist of face
- inline pair < hface3_GEO *, int > myintersection (int) ;
- inline pair < const hface3_GEO *, int > myintersection (int) const;
+ inline pair < myhface_t*, int > myintersection (int) ;
+ inline pair < const myhface_t*, int > myintersection (int) const;
virtual int nFaces() const { return 4; }
virtual int nEdges() const { return 6; }
@@ -1515,7 +1528,7 @@ public :
inline int originalVertexTwist(int, int) const;
inline int originalEdgeTwist(int, int) const;
private:
- myhface3_t * f [4] ;
+ myhface_t * f [4] ;
signed char s [4] ;
} tetra_GEO ;
@@ -1526,16 +1539,20 @@ public :
public hasFace3,
public MyAlloc
{
- protected :
+ public:
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface3_GEO myhface3_t ;
+ typedef hedge1_GEO myhedge_t ;
+ typedef hface3_GEO myhface_t ;
+
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface3_t ;
+
typedef Hface3Rule myrule_t ;
- public:
+
typedef pair < hasFace3 *, int > myneighbour_t ;
typedef pair < const hasFace3 *, int > const_myneighbour_t;
protected:
- inline Periodic3 (myhface3_t *, int, myhface3_t *, int) ;
+ inline Periodic3 (myhface_t *, int, myhface_t *, int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
@@ -1543,8 +1560,13 @@ public :
using hasFace3 :: accessPllX ;
static const int prototype [2][3] ;
inline virtual ~Periodic3 () ;
- inline hface3_GEO * myhface3 (int) ;
- inline const hface3_GEO * myhface3 (int) const ;
+
+ inline myhface_t* myhface (int) ;
+ inline const myhface_t* myhface (int) const ;
+
+ inline myhface_t* myhface3 ( int i ) { return myhface( i ); } // use the ones without number
+ inline const myhface_t* myhface3 ( int i ) const { return myhface( i ); } // use the ones without number
+
inline VertexGeo * myvertex (int) ;
inline const VertexGeo * myvertex (int) const ;
inline VertexGeo * myvertex (int,int) ;
@@ -1564,6 +1586,7 @@ public :
virtual myrule_t getrule () const = 0 ;
virtual void request (myrule_t) = 0 ;
virtual bool markNonConform () { return true; }
+ virtual void markEdgeCoarsening () { }
int tagForGlobalRefinement () ;
int tagForGlobalCoarsening () ;
int resetRefinementRequest () ;
@@ -1580,7 +1603,7 @@ public :
// return false for vertex projection
virtual bool hasVertexProjection() const { return false; }
private :
- myhface3_t * f [2] ;
+ myhface_t * f [2] ;
signed char s [2] ;
} periodic3_GEO ;
@@ -1591,17 +1614,24 @@ public :
// aus dem Element herausgeschaut wird. Gegensatz zum Tetraeder.
// Der Prototyp steht in 'gitter_geo.cc'
- typedef class Hexa : public helement_STI, public hasFace4, public MyAlloc {
- protected :
+ typedef class Hexa
+ : public helement_STI, public hasFace4, public MyAlloc
+ {
+ public :
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface4_GEO myhface4_t ;
+ typedef hedge1_GEO myhedge_t ;
+ typedef hface4_GEO myhface_t ;
+
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface4_t ;
+
typedef HexaRule myrule_t ;
typedef pair < hasFace4 *, int > myneighbour_t ;
- inline Hexa (myhface4_t *, int, myhface4_t *, int,
- myhface4_t *, int, myhface4_t *, int,
- myhface4_t *, int, myhface4_t *, int) ;
+ protected :
+ inline Hexa (myhface_t *, int, myhface_t *, int,
+ myhface_t *, int, myhface_t *, int,
+ myhface_t *, int, myhface_t *, int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
@@ -1623,12 +1653,18 @@ public :
static const vector<int> & facesNotOnFace( const int face );
inline virtual ~Hexa () ;
- inline hface4_GEO * myhface4 (int) ;
- inline const hface4_GEO * myhface4 (int) const ;
+ inline myhface_t* myhface (int) ;
+ inline const myhface_t* myhface (int) const ;
+
+ // don't use these methods use the ones without numbers
+ inline myhface_t* myhface4 ( int i ) { return myhface( i ); }
+ inline const myhface_t* myhface4 ( int i ) const { return myhface( i ); }
+
inline VertexGeo * myvertex (int) ;
inline const VertexGeo * myvertex (int) const ;
- inline hedge1_GEO * myhedge1(int);
- inline const hedge1_GEO * myhedge1(int) const;
+ inline myhedge_t* myhedge1(int);
+ inline const myhedge_t* myhedge1(int) const;
+
inline VertexGeo * myvertex (int,int) ;
inline const VertexGeo * myvertex (int,int) const ;
inline pair < hasFace4 *, int > myneighbour (int) ;
@@ -1636,8 +1672,8 @@ public :
// Dune extension
// return pair, first = pointer to face, second = twist of face
- inline pair < hface4_GEO *, int > myintersection (int) ;
- inline pair < const hface4_GEO *, int > myintersection (int) const;
+ inline pair < myhface_t*, int > myintersection (int) ;
+ inline pair < const myhface_t*, int > myintersection (int) const;
virtual int nFaces() const { return 6; }
virtual int nEdges() const { return 12; }
@@ -1655,6 +1691,7 @@ public :
virtual myrule_t requestrule () const = 0;
virtual void request (myrule_t) = 0 ;
virtual bool markNonConform () { return true; }
+ virtual void markEdgeCoarsening () { }
int tagForGlobalRefinement () ;
int tagForGlobalCoarsening () ;
int resetRefinementRequest () ;
@@ -1675,7 +1712,7 @@ public :
inline int evalVertexTwist(int, int) const;
inline int evalEdgeTwist(int, int) const;
private:
- myhface4_t * f [6] ;
+ myhface_t * f [6] ;
signed char s [6] ;
} hexa_GEO ;
@@ -1686,16 +1723,20 @@ public :
public hasFace4,
public MyAlloc
{
- protected :
+ public:
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface4_GEO myhface4_t ;
+ typedef hedge1_GEO myhedge_t ;
+ typedef hface4_GEO myhface_t ;
+
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface4_t ;
+
typedef Hface4Rule myrule_t ;
- public:
+
typedef pair < hasFace4 *, int > myneighbour_t ;
typedef pair < const hasFace4 *, int > const_myneighbour_t;
protected:
- inline Periodic4 (myhface4_t *, int, myhface4_t *, int) ;
+ inline Periodic4 (myhface_t *, int, myhface_t *, int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
@@ -1703,8 +1744,11 @@ public :
using hasFace4 :: accessPllX ;
static const int prototype [2][4] ;
inline virtual ~Periodic4 () ;
- inline hface4_GEO * myhface4 (int) ;
- inline const hface4_GEO * myhface4 (int) const ;
+ inline myhface_t* myhface (int) ;
+ inline const myhface_t* myhface (int) const ;
+ inline myhface_t* myhface4 ( int i ) { return myhface( i ); }
+ inline const myhface_t* myhface4 ( int i ) const { return myhface( i ); }
+
inline VertexGeo * myvertex (int) ;
inline const VertexGeo * myvertex (int) const ;
inline VertexGeo * myvertex (int,int) ;
@@ -1728,6 +1772,7 @@ public :
virtual myrule_t getrule () const = 0 ;
virtual void request (myrule_t) = 0 ;
virtual bool markNonConform () { return true; }
+ virtual void markEdgeCoarsening () { }
int tagForGlobalRefinement () ;
int tagForGlobalCoarsening () ;
int resetRefinementRequest () ;
@@ -1740,7 +1785,7 @@ public :
// returns false because only bnd segments have projections
virtual bool hasVertexProjection () const { return false; }
private :
- myhface4_t * f [2] ;
+ myhface_t * f [2] ;
signed char s [2] ;
} periodic4_GEO ;
@@ -1750,17 +1795,20 @@ public :
// auf die Randfl"ache geschaut wird. Das Vierecksrandelement hat die
// entgegengesetzte Konvention.
- typedef class hbndseg3 : public hbndseg_STI, public hasFace3, public MyAlloc {
+ typedef class hbndseg3
+ : public hbndseg_STI, public hasFace3, public MyAlloc
+ {
public :
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface3_GEO myhface3_t ;
+ typedef hedge1_GEO myhedge_t ;
typedef hface3_GEO myhface_t ;
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface3_t ;
typedef Hface3Rule myrule_t ;
typedef hbndseg_STI :: bnd_t bnd_t;
protected :
- inline hbndseg3 (myhface3_t *,int) ;
+ inline hbndseg3 (myhface_t *,int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
inline bool lockedAgainstCoarsening () const { return false ; }
@@ -1770,7 +1818,8 @@ public :
inline myrule_t getrule () const ;
virtual bool refineLikeElement (balrule_t) = 0 ;
inline myvertex_t * myvertex (int,int) const ;
- inline myhface3_t * myhface3 (int) const ;
+ inline myhface_t * myhface3 ( int i ) const { return myhface( i ); } // use the method without number
+ inline myhface_t * myhface (int) const ;
inline int twist (int) const ;
inline hface3_GEO * subface3 (int,int) const ;
@@ -1778,6 +1827,7 @@ public :
virtual bool isperiodic() const { return false; }
virtual bool markNonConform () { return true; }
+ virtual void markEdgeCoarsening () { }
virtual int nChild () const;
// just returns level
virtual int nbLevel() const { return level(); }
@@ -1796,7 +1846,7 @@ public :
ProjectVertex* projection() { return ( this->isGhost() ) ? 0 : _face->myvertex(0)->myGrid()->vertexProjection(); }
private:
- myhface3_t * _face ;
+ myhface_t * _face ;
int _twist ;
public:
} hbndseg3_GEO ;
@@ -1805,14 +1855,17 @@ public :
typedef class hbndseg4 : public hbndseg_STI, public hasFace4, public MyAlloc {
public :
typedef VertexGeo myvertex_t ;
- typedef hedge1_GEO myhedge1_t ;
- typedef hface4_GEO myhface4_t ;
+ typedef hedge1_GEO myhedge_t ;
typedef hface4_GEO myhface_t ;
+
+ typedef myhedge_t myhedge1_t ;
+ typedef myhface_t myhface4_t
+ ;
typedef Hface4Rule myrule_t ;
typedef hbndseg_STI :: bnd_t bnd_t;
protected :
- inline hbndseg4 (myhface4_t *,int) ;
+ inline hbndseg4 (myhface_t *,int) ;
inline int postRefinement () ;
inline int preCoarsening () ;
inline bool lockedAgainstCoarsening () const { return false ; }
@@ -1822,11 +1875,13 @@ public :
inline myrule_t getrule () const ;
virtual bool refineLikeElement (balrule_t) = 0 ;
inline myvertex_t * myvertex (int,int) const ;
- inline myhface4_t * myhface4 (int) const ;
+ inline myhface_t * myhface (int) const ;
+ inline myhface_t * myhface4 ( int i ) const { return myhface( i ); }
inline int twist (int) const ;
inline hface4_GEO * subface4 (int,int) const ;
virtual bool markNonConform () { return true; }
+ virtual void markEdgeCoarsening () { }
virtual bool isboundary() const { return true; }
virtual bool isperiodic() const { return false; }
virtual int nChild () const;
@@ -1842,7 +1897,7 @@ public :
ProjectVertex* projection() { return ( this->isGhost() ) ? 0 : _face->myvertex(0)->myGrid()->vertexProjection(); }
private:
- myhface4_t * _face ;
+ myhface_t * _face ;
int _twist ;
public:
@@ -2026,6 +2081,7 @@ protected :
// methods for refining and coarsening
virtual bool refine () ;
virtual bool markNonConform () ;
+ virtual void markEdgeCoarsening () ;
virtual void coarse () ;
virtual Makrogitter & container () = 0 ;
@@ -2078,6 +2134,10 @@ public :
virtual void tovtk( const std::string &fn) ;
protected:
+ template <class element_t, class bnd_t>
+ void tovtkImpl( const std::string &fn,
+ const int, const element_t*, const bnd_t* ) ;
+
template <class ostream_t>
void backupImpl( ostream_t& );
@@ -2096,6 +2156,10 @@ protected:
friend class LevelIterator < hbndseg_STI > ;
friend class LevelIterator < hedge_STI > ;
friend class LevelIterator < hface_STI > ;
+
+ public:
+ // flags to say if an edge can be coarsened during conform bisection
+ mutable std::vector<bool> edgeCoarseningFlags_;
} ;
// --endGitter
@@ -2542,7 +2606,7 @@ inline ostream& operator<< (ostream& s, const Gitter :: Geometric :: Tetra* tetr
{
s << "T-Face " << i << " (tw = " << tetra->twist( i ) << ") ";
for(int j=0; j<3; ++j)
- s << tetra->myhface3( i )->myvertex( j ) << " " ;
+ s << tetra->myhface( i )->myvertex( j ) << " " ;
s << endl;
}
s << endl;
@@ -2820,6 +2884,18 @@ inline const Gitter :: Geometric :: hedge1 :: myvertex_t * Gitter :: Geometric :
return i == 1 ? v1 : v0 ;
}
+
+
+inline bool Gitter :: hedge :: canCoarsen( std::vector<bool> &edgeCoarseningFlags ) const
+{
+ int idx = this->getIndex();
+ if ( !edgeCoarseningFlags[ idx ] ) return false;
+ if ( this->down() ) return this->down()->canCoarsen( edgeCoarseningFlags );
+ if ( this->next() ) return this->next()->canCoarsen( edgeCoarseningFlags );
+ return true;
+}
+
+
// # # ##### ######
// # # ###### ## #### ###### # # # # # # # ######
// # # # # # # # # # # # # # # #
@@ -3031,8 +3107,8 @@ Gitter :: Geometric :: hface3 :: face3Neighbour :: setPrevRear ( )
inline void Gitter :: Geometric :: hface3 :: face3Neighbour :: operator = (const face3Neighbour & n)
{
- frontList_ = n.frontList_;
- rearList_ = n.rearList_;
+ frontList_.clear();
+ rearList_.clear();
_faceFront = n._faceFront;
_faceRear = n._faceRear;
_numFront = n._numFront;
@@ -3043,6 +3119,7 @@ inline void Gitter :: Geometric :: hface3 :: face3Neighbour :: operator = (const
inline int Gitter :: Geometric :: hface3 :: face3Neighbour :: complete (const face3Neighbour & n)
{
int ret = 0;
+ // assert(0);
if( front() == null )
{
@@ -3109,7 +3186,6 @@ inline void Gitter :: Geometric :: hface3 :: attachElement (const pair < myconne
ref ++ ;
if (t>=0 && nb.frontList_.size()==0)
ref ++ ;
-
if( t < 0 )
nb.setRear( p );
else
@@ -3472,8 +3548,8 @@ inline ostream &operator<< ( ostream &out, const Gitter :: Geometric :: TetraRul
// # ###### # # # # #
inline Gitter :: Geometric :: Tetra ::
-Tetra (myhface3_t * f0, int t0, myhface3_t * f1, int t1,
- myhface3_t * f2, int t2, myhface3_t * f3, int t3)
+Tetra (myhface_t * f0, int t0, myhface_t * f1, int t1,
+ myhface_t * f2, int t2, myhface_t * f3, int t3)
{
(f [0] = f0)->attachElement (pair < hasFace3 *, int > (InternalHasFace3 ()(this), 0),(s [0] = t0)) ;
(f [1] = f1)->attachElement (pair < hasFace3 *, int > (InternalHasFace3 ()(this), 1),(s [1] = t1)) ;
@@ -3495,12 +3571,12 @@ inline int Gitter :: Geometric :: Tetra :: twist (int i) const {
return s [i] ;
}
-inline Gitter :: Geometric :: Tetra :: myhface3_t * Gitter :: Geometric :: Tetra :: myhface3 (int i) {
+inline Gitter :: Geometric :: Tetra :: myhface_t * Gitter :: Geometric :: Tetra :: myhface (int i) {
assert (i < 4) ;
return f [i] ;
}
-inline const Gitter :: Geometric :: Tetra :: myhface3_t * Gitter :: Geometric :: Tetra :: myhface3 (int i) const {
+inline const Gitter :: Geometric :: Tetra :: myhface_t * Gitter :: Geometric :: Tetra :: myhface (int i) const {
assert (i < 4) ;
return f [i] ;
}
@@ -3546,7 +3622,7 @@ inline Gitter :: Geometric :: Tetra :: myhedge1_t * Gitter :: Geometric :: Tetra
typedef Gitter::Geometric::Tetra ThisType;
- return myhface3(ThisType::edgeMap[edge][0])->
+ return myhface(ThisType::edgeMap[edge][0])->
myhedge1(evalEdgeTwist(ThisType::edgeMap[edge][0],ThisType::edgeMap[edge][1]));
}
@@ -3556,17 +3632,17 @@ inline const Gitter :: Geometric :: Tetra :: myhedge1_t * Gitter :: Geometric ::
typedef Gitter::Geometric::Tetra ThisType;
- return myhface3(ThisType::edgeMap[edge][0])->
+ return myhface(ThisType::edgeMap[edge][0])->
myhedge1(evalEdgeTwist(ThisType::edgeMap[edge][0],ThisType::edgeMap[edge][1]));
}
inline Gitter :: Geometric :: Tetra :: myvertex_t * Gitter :: Geometric :: Tetra :: myvertex (int i, int j) {
- return myhface3(i)->myvertex(evalVertexTwist(i, j));
+ return myhface(i)->myvertex(evalVertexTwist(i, j));
}
inline const Gitter :: Geometric :: Tetra :: myvertex_t * Gitter :: Geometric :: Tetra :: myvertex (int i, int j) const {
- return myhface3(i)->myvertex(evalVertexTwist(i, j));
+ return myhface(i)->myvertex(evalVertexTwist(i, j));
}
//- --tetramyvertex
@@ -3582,24 +3658,24 @@ inline const Gitter :: Geometric :: Tetra :: myvertex_t * Gitter :: Geometric ::
inline pair < Gitter :: Geometric :: hasFace3 *, int > Gitter :: Geometric :: Tetra :: myneighbour (int i)
{
- return twist (i) < 0 ? myhface3 (i)->nb.front () : myhface3 (i)->nb.rear ();
+ return twist (i) < 0 ? myhface( i )->nb.front () : myhface( i )->nb.rear ();
}
inline pair < const Gitter :: Geometric :: hasFace3 *, int > Gitter :: Geometric :: Tetra :: myneighbour (int i) const {
- return twist (i) < 0 ? pair < const hasFace3 *, int > (myhface3 (i)->nb.front ().first, myhface3 (i)->nb.front ().second)
- : pair < const hasFace3 *, int > (myhface3 (i)->nb.rear ().first, myhface3 (i)->nb.rear ().second) ;
+ return twist (i) < 0 ? pair < const hasFace3 *, int > (myhface( i )->nb.front ().first, myhface( i )->nb.front ().second)
+ : pair < const hasFace3 *, int > (myhface( i )->nb.rear ().first, myhface( i )->nb.rear ().second) ;
}
inline pair < Gitter :: Geometric :: hface3_GEO *, int > Gitter :: Geometric :: Tetra :: myintersection (int i)
{
// return pair, first = pointer to face, second = twist of face
- return pair< Gitter::Geometric::hface3_GEO *,int> (myhface3 (i) ,twist (i));
+ return pair< Gitter::Geometric::hface3_GEO *,int> (myhface( i ) ,twist (i));
}
inline pair < const Gitter :: Geometric :: hface3_GEO *, int > Gitter :: Geometric :: Tetra :: myintersection (int i) const
{
// return pair, first = pointer to face, second = twist of face
- return pair< const Gitter::Geometric::hface3_GEO * , int > (myhface3 (i) , twist (i) );
+ return pair< const Gitter::Geometric::hface3_GEO * , int > (myhface( i ) , twist (i) );
}
inline int Gitter :: Geometric :: Tetra :: postRefinement () {
@@ -3619,7 +3695,7 @@ inline int Gitter :: Geometric :: Tetra :: preCoarsening () {
// # ###### # # # #### ##### # #### #####
inline Gitter :: Geometric :: Periodic3 ::
-Periodic3 (myhface3_t * f0, int t0, myhface3_t * f1, int t1)
+Periodic3 (myhface_t * f0, int t0, myhface_t * f1, int t1)
{
(f [0] = f0)->attachElement (pair < hasFace3 *, int > (InternalHasFace3 ()(this), 0),(s [0] = t0)) ;
(f [1] = f1)->attachElement (pair < hasFace3 *, int > (InternalHasFace3 ()(this), 1),(s [1] = t1)) ;
@@ -3637,23 +3713,23 @@ inline int Gitter :: Geometric :: Periodic3 :: twist (int i) const {
return s [i] ;
}
-inline Gitter :: Geometric :: Periodic3 :: myhface3_t * Gitter :: Geometric :: Periodic3 :: myhface3 (int i) {
+inline Gitter :: Geometric :: Periodic3 :: myhface_t * Gitter :: Geometric :: Periodic3 :: myhface (int i) {
assert (0 <= i && i < 2) ;
return f [i] ;
}
-inline const Gitter :: Geometric :: Periodic3 :: myhface3_t * Gitter :: Geometric :: Periodic3 :: myhface3 (int i) const {
+inline const Gitter :: Geometric :: Periodic3 :: myhface_t * Gitter :: Geometric :: Periodic3 :: myhface (int i) const {
assert (0 <= i && i < 2) ;
return f [i] ;
}
inline Gitter :: Geometric :: Periodic3 :: myvertex_t * Gitter :: Geometric :: Periodic3 :: myvertex (int i, int j) {
assert (0 <= i && i < 2) ;
- return (twist(i) < 0) ? myhface3(i)->myvertex((7 - j + twist(i)) % 3) : myhface3(i)->myvertex((j + twist(i)) % 3) ;
+ return (twist(i) < 0) ? myhface( i )->myvertex((7 - j + twist(i)) % 3) : myhface( i )->myvertex((j + twist(i)) % 3) ;
}
inline const Gitter :: Geometric :: Periodic3 :: myvertex_t * Gitter :: Geometric :: Periodic3 :: myvertex (int i, int j) const {
- return (twist(i) < 0) ? myhface3(i)->myvertex((7 - j + twist(i)) % 3) : myhface3(i)->myvertex((j + twist(i)) % 3) ;
+ return (twist(i) < 0) ? myhface( i )->myvertex((7 - j + twist(i)) % 3) : myhface( i )->myvertex((j + twist(i)) % 3) ;
}
inline Gitter :: Geometric :: Periodic3 :: myvertex_t * Gitter :: Geometric :: Periodic3 :: myvertex (int i) {
@@ -3675,13 +3751,13 @@ inline const Gitter :: Geometric :: Periodic3 :: myvertex_t * Gitter :: Geometri
inline pair < Gitter :: Geometric :: hasFace3 *, int > Gitter :: Geometric :: Periodic3 :: myneighbour (int i) {
assert (0 <= i && i < 2) ;
- return twist (i) < 0 ? myhface3 (i)->nb.front () : myhface3 (i)->nb.rear () ;
+ return twist (i) < 0 ? myhface( i )->nb.front () : myhface( i )->nb.rear () ;
}
inline pair < const Gitter :: Geometric :: hasFace3 *, int > Gitter :: Geometric :: Periodic3 :: myneighbour (int i) const {
assert (0 <= i && i < 2) ;
- return twist (i) < 0 ? pair < const hasFace3 *, int > (myhface3 (i)->nb.front ().first, myhface3 (i)->nb.front ().second)
- : pair < const hasFace3 *, int > (myhface3 (i)->nb.rear ().first, myhface3 (i)->nb.rear ().second) ;
+ return twist (i) < 0 ? pair < const hasFace3 *, int > (myhface( i )->nb.front ().first, myhface( i )->nb.front ().second)
+ : pair < const hasFace3 *, int > (myhface( i )->nb.rear ().first, myhface( i )->nb.rear ().second) ;
}
inline int Gitter :: Geometric :: Periodic3 :: postRefinement () {
@@ -3703,7 +3779,7 @@ inline int Gitter :: Geometric :: Periodic3 :: preCoarsening () {
// # ###### # # # #### ##### # #### #
inline Gitter :: Geometric :: Periodic4 ::
-Periodic4 (myhface4_t * f0, int t0, myhface4_t * f1, int t1)
+Periodic4 (myhface_t * f0, int t0, myhface_t * f1, int t1)
{
(f [0] = f0)->attachElement (pair < hasFace4 *, int > (InternalHasFace4 ()(this), 0),(s [0] = t0)) ;
(f [1] = f1)->attachElement (pair < hasFace4 *, int > (InternalHasFace4 ()(this), 1),(s [1] = t1)) ;
@@ -3721,24 +3797,24 @@ inline int Gitter :: Geometric :: Periodic4 :: twist (int i) const {
return s [i] ;
}
-inline Gitter :: Geometric :: Periodic4 :: myhface4_t * Gitter :: Geometric :: Periodic4 :: myhface4 (int i) {
+inline Gitter :: Geometric :: Periodic4 :: myhface_t * Gitter :: Geometric :: Periodic4 :: myhface (int i) {
assert (0 <= i && i < 2) ;
return f [i] ;
}
-inline const Gitter :: Geometric :: Periodic4 :: myhface4_t * Gitter :: Geometric :: Periodic4 :: myhface4 (int i) const {
+inline const Gitter :: Geometric :: Periodic4 :: myhface_t * Gitter :: Geometric :: Periodic4 :: myhface (int i) const {
assert (0 <= i && i < 2) ;
return f [i] ;
}
inline Gitter :: Geometric :: Periodic4 :: myvertex_t * Gitter :: Geometric :: Periodic4 :: myvertex (int i, int j) {
assert (0 <= i && i < 2) ;
- return (twist(i) < 0) ? myhface4(i)->myvertex((9 - j + twist(i)) % 4) : myhface4(i)->myvertex((j + twist(i)) % 4) ;
+ return (twist(i) < 0) ? myhface(i)->myvertex((9 - j + twist(i)) % 4) : myhface(i)->myvertex((j + twist(i)) % 4) ;
}
inline const Gitter :: Geometric :: Periodic4 :: myvertex_t * Gitter :: Geometric :: Periodic4 :: myvertex (int i, int j) const {
assert (0 <= i && i < 2) ;
- return (twist(i) < 0) ? myhface4(i)->myvertex((9 - j + twist(i)) % 4) : myhface4(i)->myvertex((j + twist(i)) % 4) ;
+ return (twist(i) < 0) ? myhface(i)->myvertex((9 - j + twist(i)) % 4) : myhface(i)->myvertex((j + twist(i)) % 4) ;
}
inline Gitter :: Geometric :: Periodic4 :: myvertex_t * Gitter :: Geometric :: Periodic4 :: myvertex (int i) { // ok
@@ -3753,13 +3829,13 @@ inline const Gitter :: Geometric :: Periodic4 :: myvertex_t * Gitter :: Geometri
inline pair < Gitter :: Geometric :: hasFace4 *, int > Gitter :: Geometric :: Periodic4 :: myneighbour (int i) {
assert (0 <= i && i < 2) ;
- return twist (i) < 0 ? myhface4 (i)->nb.front () : myhface4 (i)->nb.rear () ;
+ return twist (i) < 0 ? myhface( i )->nb.front () : myhface( i )->nb.rear () ;
}
inline pair < const Gitter :: Geometric :: hasFace4 *, int > Gitter :: Geometric :: Periodic4 :: myneighbour (int i) const {
assert (0 <= i && i < 2) ;
- return twist (i) < 0 ? pair < const hasFace4 *, int > (myhface4 (i)->nb.front ().first, myhface4 (i)->nb.front ().second)
- : pair < const hasFace4 *, int > (myhface4 (i)->nb.rear ().first, myhface4 (i)->nb.rear ().second) ;
+ return twist (i) < 0 ? pair < const hasFace4 *, int > (myhface( i )->nb.front ().first, myhface( i )->nb.front ().second)
+ : pair < const hasFace4 *, int > (myhface( i )->nb.rear ().first, myhface( i )->nb.rear ().second) ;
}
inline int Gitter :: Geometric :: Periodic4 :: postRefinement () {
@@ -3835,9 +3911,9 @@ inline ostream &operator<< ( ostream &out, const Gitter :: Geometric :: HexaRule
// # # ###### # # # #
inline Gitter :: Geometric :: Hexa ::
-Hexa (myhface4_t * f0, int t0, myhface4_t * f1, int t1,
- myhface4_t * f2, int t2, myhface4_t * f3, int t3,
- myhface4_t * f4, int t4, myhface4_t * f5, int t5)
+Hexa (myhface_t * f0, int t0, myhface_t * f1, int t1,
+ myhface_t * f2, int t2, myhface_t * f3, int t3,
+ myhface_t * f4, int t4, myhface_t * f5, int t5)
{
(f [0] = f0)->attachElement (pair < hasFace4 *, int > (InternalHasFace4 ()(this), 0),(s [0] = t0)) ;
(f [1] = f1)->attachElement (pair < hasFace4 *, int > (InternalHasFace4 ()(this), 1),(s [1] = t1)) ;
@@ -3863,24 +3939,24 @@ inline int Gitter :: Geometric :: Hexa :: twist (int i) const {
return s [i] ;
}
-inline Gitter :: Geometric :: Hexa :: myhface4_t * Gitter :: Geometric :: Hexa :: myhface4 (int i) {
+inline Gitter :: Geometric :: Hexa :: myhface_t * Gitter :: Geometric :: Hexa :: myhface (int i) {
assert (i < 6) ;
return f [i] ;
}
-inline const Gitter :: Geometric :: Hexa :: myhface4_t * Gitter :: Geometric :: Hexa :: myhface4 (int i) const {
+inline const Gitter :: Geometric :: Hexa :: myhface_t * Gitter :: Geometric :: Hexa :: myhface (int i) const {
assert (i < 6) ;
return f [i] ;
}
inline Gitter :: Geometric :: Hexa :: myvertex_t * Gitter :: Geometric :: Hexa :: myvertex (int i, int j)
{
- return myhface4(i)->myvertex(evalVertexTwist(i, j));
+ return myhface(i)->myvertex(evalVertexTwist(i, j));
}
inline const Gitter :: Geometric :: Hexa :: myvertex_t * Gitter :: Geometric :: Hexa :: myvertex (int i, int j) const
{
- return myhface4(i)->myvertex(evalVertexTwist(i, j));
+ return myhface(i)->myvertex(evalVertexTwist(i, j));
}
inline Gitter :: Geometric :: Hexa :: myvertex_t *
@@ -3901,7 +3977,7 @@ inline Gitter :: Geometric :: Hexa :: myhedge1_t * Gitter :: Geometric :: Hexa :
assert (0 <= i && i < 12);
typedef Gitter::Geometric::Hexa MyType;
- return myhface4(MyType::edgeMap[i][0])->
+ return myhface(MyType::edgeMap[i][0])->
myhedge1(evalEdgeTwist(MyType::edgeMap[i][0], MyType::edgeMap[i][1]));
}
@@ -3909,27 +3985,27 @@ inline const Gitter :: Geometric :: Hexa :: myhedge1_t * Gitter :: Geometric ::
assert (0 <= i && i < 12);
typedef Gitter::Geometric::Hexa MyType;
- return myhface4(MyType::edgeMap[i][0])->
+ return myhface(MyType::edgeMap[i][0])->
myhedge1(evalEdgeTwist(MyType::edgeMap[i][0], MyType::edgeMap[i][1]));
}
inline pair < Gitter :: Geometric :: hasFace4 *, int > Gitter :: Geometric :: Hexa :: myneighbour (int i) {
- return twist (i) < 0 ? myhface4 (i)->nb.front () : myhface4 (i)->nb.rear () ;
+ return twist (i) < 0 ? myhface( i )->nb.front () : myhface( i )->nb.rear () ;
}
inline pair < const Gitter :: Geometric :: hasFace4 *, int > Gitter :: Geometric :: Hexa :: myneighbour (int i) const {
- return twist (i) < 0 ? pair < const hasFace4 *, int > (myhface4 (i)->nb.front ().first, myhface4 (i)->nb.front ().second)
- : pair < const hasFace4 *, int > (myhface4 (i)->nb.rear ().first, myhface4 (i)->nb.rear ().second) ;
+ return twist (i) < 0 ? pair < const hasFace4 *, int > (myhface( i )->nb.front ().first, myhface( i )->nb.front ().second)
+ : pair < const hasFace4 *, int > (myhface( i )->nb.rear ().first, myhface( i )->nb.rear ().second) ;
}
inline pair<Gitter::Geometric::hface4_GEO*, int>
Gitter::Geometric::Hexa::myintersection(int i) {
- return make_pair(myhface4(i), twist(i));
+ return make_pair(myhface(i), twist(i));
}
inline pair<const Gitter::Geometric::hface4_GEO*, int>
Gitter::Geometric::Hexa::myintersection(int i) const {
- return make_pair(myhface4(i), twist(i));
+ return make_pair(myhface(i), twist(i));
}
inline int Gitter :: Geometric :: Hexa :: postRefinement () {
@@ -3981,7 +4057,7 @@ inline int Gitter :: Geometric :: Hexa :: originalEdgeTwist (int face, int edge)
// # # ##### # # ##### #### ###### #### #####
inline Gitter :: Geometric :: hbndseg3 ::
-hbndseg3 (myhface3_t * a, int b)
+hbndseg3 (myhface_t * a, int b)
: _face( a ),
_twist (b)
{
@@ -3998,7 +4074,7 @@ inline void Gitter :: Geometric :: hbndseg3 :: attachleafs ()
{
this->addleaf();
- myhface3_t & face = *(myhface3(0));
+ myhface_t & face = *(myhface(0));
face.addleaf();
for (int i=0; i<3; ++i)
{
@@ -4011,7 +4087,7 @@ inline void Gitter :: Geometric :: hbndseg3 :: detachleafs ()
{
this->removeleaf();
- myhface3_t & face = *(myhface3(0));
+ myhface_t & face = *(myhface(0));
face.removeleaf();
for (int i=0; i<3; ++i)
{
@@ -4025,7 +4101,7 @@ inline int Gitter :: Geometric :: hbndseg3 :: postRefinement ()
ProjectVertexPair pv( projection(), segmentIndex() );
if ( pv.first )
{
- myhface3(0)->projectVertex( pv );
+ myhface(0)->projectVertex( pv );
}
return 0 ;
}
@@ -4040,21 +4116,21 @@ inline int Gitter :: Geometric :: hbndseg3 :: twist (int i) const {
return _twist ;
}
-inline Gitter :: Geometric :: hbndseg3 :: myhface3_t * Gitter :: Geometric :: hbndseg3 :: myhface3 (int i) const {
+inline Gitter :: Geometric :: hbndseg3 :: myhface_t * Gitter :: Geometric :: hbndseg3 :: myhface (int i) const {
assert (i == 0) ;
return _face ;
}
inline Gitter :: Geometric :: hbndseg3 :: myvertex_t * Gitter :: Geometric :: hbndseg3 :: myvertex (int,int j) const {
- return (twist (0) < 0) ? myhface3 (0)->myvertex ((7 - j + twist (0)) % 3) : myhface3 (0)->myvertex ((j + twist (0)) % 3) ;
+ return (twist (0) < 0) ? myhface(0)->myvertex ((7 - j + twist (0)) % 3) : myhface(0)->myvertex ((j + twist (0)) % 3) ;
}
-inline Gitter :: Geometric :: hbndseg3 :: myhface3_t * Gitter :: Geometric :: hbndseg3 :: subface3 (int,int i) const {
- return myhface3 (0)->subface3 (i) ;
+inline Gitter :: Geometric :: hbndseg3 :: myhface_t * Gitter :: Geometric :: hbndseg3 :: subface3 (int,int i) const {
+ return myhface(0)->subface3 (i) ;
}
inline Gitter :: Geometric :: hbndseg3 :: myrule_t Gitter :: Geometric :: hbndseg3 :: getrule () const {
- return myhface3 (0)->getrule () ;
+ return myhface(0)->getrule () ;
}
inline int Gitter :: Geometric :: hbndseg3 :: nChild () const {
@@ -4071,7 +4147,7 @@ inline int Gitter :: Geometric :: hbndseg3 :: nChild () const {
// # # # # # ## # # # # # # # #
// # # ##### # # ##### #### ###### #### #
-inline Gitter :: Geometric :: hbndseg4 :: hbndseg4 (myhface4_t * a, int b)
+inline Gitter :: Geometric :: hbndseg4 :: hbndseg4 (myhface_t * a, int b)
: _face( a ),
_twist (b)
{
@@ -4089,7 +4165,7 @@ inline void Gitter :: Geometric :: hbndseg4 :: attachleafs ()
assert(this->leafRefCount()==0);
this->addleaf();
- hface4_GEO & face = *(myhface4(0));
+ myhface_t& face = *(myhface(0));
face.addleaf();
for (int i=0; i<4; ++i)
{
@@ -4103,7 +4179,7 @@ inline void Gitter :: Geometric :: hbndseg4 :: detachleafs ()
assert(this->leafRefCount()==1);
this->removeleaf();
- hface4_GEO & face = *(myhface4(0));
+ myhface_t& face = *(myhface(0));
face.removeleaf();
for (int i=0; i<4; ++i)
{
@@ -4117,7 +4193,7 @@ inline int Gitter :: Geometric :: hbndseg4 :: postRefinement ()
ProjectVertexPair pv( projection(), segmentIndex() );
if( pv.first )
{
- myhface4(0)->projectVertex( pv );
+ myhface(0)->projectVertex( pv );
}
return 0 ;
}
@@ -4133,21 +4209,21 @@ inline int Gitter :: Geometric :: hbndseg4 :: twist (int i) const
return _twist ;
}
-inline Gitter :: Geometric :: hbndseg4 :: myhface4_t * Gitter :: Geometric :: hbndseg4 :: myhface4 (int i) const {
+inline Gitter :: Geometric :: hbndseg4 :: myhface_t * Gitter :: Geometric :: hbndseg4 :: myhface (int i) const {
assert (i == 0) ;
return _face ;
}
inline Gitter :: Geometric :: hbndseg4 :: myvertex_t * Gitter :: Geometric :: hbndseg4 :: myvertex (int,int j) const {
- return (twist (0) < 0) ? myhface4 (0)->myvertex ((9 - j + twist (0)) % 4) : myhface4 (0)->myvertex ((j + twist (0)) % 4) ;
+ return (twist (0) < 0) ? myhface(0)->myvertex ((9 - j + twist (0)) % 4) : myhface(0)->myvertex ((j + twist (0)) % 4) ;
}
-inline Gitter :: Geometric :: hbndseg4 :: myhface4_t * Gitter :: Geometric :: hbndseg4 :: subface4 (int,int i) const {
- return myhface4 (0)->subface4 (i) ;
+inline Gitter :: Geometric :: hbndseg4 :: myhface_t * Gitter :: Geometric :: hbndseg4 :: subface4 (int,int i) const {
+ return myhface(0)->subface4 (i) ;
}
inline Gitter :: Geometric :: hbndseg4 :: myrule_t Gitter :: Geometric :: hbndseg4 :: getrule () const {
- return myhface4 (0)->getrule () ;
+ return myhface(0)->getrule () ;
}
inline int Gitter :: Geometric :: hbndseg4 :: nChild () const {
diff --git a/src/serial/gitter_tetra_top.cc b/src/serial/gitter_tetra_top.cc
index d1df8e0910dd5b348676d6af78be1adfbe034ae1..3fdc72f02c3b102f6723c8987d6d5677b1abccfa 100644
--- a/src/serial/gitter_tetra_top.cc
+++ b/src/serial/gitter_tetra_top.cc
@@ -422,7 +422,6 @@ template < class A > bool Hface3Top < A > :: coarse ()
// vollst"andigt, die eventuell durch Elementvergr"oberung
// durcheinander gekommen war. Die Vergr"oberung geht dann
// auf das n"achste Level "uber.
-
if (f->ref) {
if (f->ref == 1) f->nb.complete (this->nb) ;
f->coarse () ;
@@ -2103,6 +2102,7 @@ template < class A > bool TetraTop < A > :: coarse ()
{
if (this->leaf ())
{
+ if (!up()) return false;
assert (_req == myrule_t :: nosplit || _req == myrule_t :: crs) ;
myrule_t w = _req ;
_req = myrule_t :: nosplit ;
@@ -2132,12 +2132,34 @@ template < class A > bool TetraTop < A > :: coarse ()
// not faces that are not leaf
if (x)
{
+ // test marking on refinement edge
+ assert( this->nEdges() == 6 );
+ // int e = getrule()-2; // did not work....
+ for (int e=0;e<6;++e)
+ {
+ myhedge1_t *edge = this->myhedge1(e);
+ if ( edge->down() )
+ {
+ int idx = edge->getIndex();
+ if ( !this->myGrid()->edgeCoarseningFlags_[ idx ] )
+ return false;
+ // we need to make sure that none of the children of this
+ // refinement edge should not be removed
+ // (we could not go up on edges during marking so we go down here)
+ if (! edge->down()->canCoarsen( this->myGrid()->edgeCoarseningFlags_ ) )
+ {
+ return false;
+ }
+ }
+ }
+
this->preCoarsening () ;
this->attachleafs();
delete _inner ;
_inner = 0 ;
+#if 0
// for bisection refinement we have to again
// set the face neighbours, since they have been overwritten
// by the refined element
@@ -2158,6 +2180,7 @@ template < class A > bool TetraTop < A > :: coarse ()
}
}
}
+#endif
// reset refinement rule
_rule = myrule_t :: nosplit ;
diff --git a/src/serial/gitter_tetra_top.h b/src/serial/gitter_tetra_top.h
index 3496dd23c581d12c6780641ba651aced1c3a77b3..e423a744ed011fa225c5087f2451ac7d659486b6 100644
--- a/src/serial/gitter_tetra_top.h
+++ b/src/serial/gitter_tetra_top.h
@@ -403,7 +403,6 @@ template < class A > class TetraTop : public A
assert( this->nEdges() == 6 );
for (int e=0; e < 6; ++e)
{
- // std::cout << this->myhedge1(e) << std::endl;
if( this->myhedge1(e)->down() )
{
this->request ( myrule_t :: bisect );
@@ -412,6 +411,24 @@ template < class A > class TetraTop : public A
}
return true;
}
+ virtual void markEdgeCoarsening ()
+ {
+ assert( this->nEdges() == 6 );
+ if (!this->up()) return;
+ for (int e=0; e<6; ++e)
+ {
+ myhedge1_t *edge = this->up()->myhedge1(e);
+ // if (e == (int)(this->up()->getrule())-2) // does not seem to // work...
+ if (_req == myrule_t :: crs && edge->down()) // the father of a leaf element can only have one non leaf edge
+ { // we leave thing as they are
+ }
+ else
+ {
+ int idx = edge->getIndex();
+ this->myGrid()->edgeCoarseningFlags_[ idx ] = false;
+ }
+ }
+ }
void refineImmediate (myrule_t) ;
inline void append (innertetra_t * h) ;