Commit f7da1fc5 authored by Andreas Wagner's avatar Andreas Wagner
Browse files

Merge branch 'wagnandr/facedoffunction' into 'master'

Wagnandr/facedoffunction

See merge request !446
parents e2ca9545 ab90e5bd
Pipeline #33856 passed with stages
in 197 minutes and 35 seconds
......@@ -17,18 +17,17 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <core/timing/Timer.h>
#include <core/Environment.h>
#include <core/timing/Timer.h>
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
using walberla::real_t;
using walberla::uint_t;
......
......@@ -24,8 +24,8 @@
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/P1StokesOperator.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesProlongation.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesRestriction.hpp"
......@@ -181,7 +181,7 @@ int main( int argc, char* argv[] )
{
WALBERLA_LOG_PROGRESS_ON_ROOT( "Solving Stokes system..." )
f_dg->interpolateExtended( expr_f, {c_old.get()}, maxLevel );
f_dg->interpolate( expr_f, { *c_old }, maxLevel );
f->uvw[0].integrateDG( *f_dg, *n_x, maxLevel, hyteg::All );
f->uvw[1].integrateDG( *f_dg, *n_y, maxLevel, hyteg::All );
......
......@@ -77,7 +77,7 @@ void VTKDGDoFWriter::write( const VTKOutput& mgr, std::ostream& output, const ui
template < typename value_t >
void VTKDGDoFWriter::writeScalarFunction( std::ostream& output,
const DGFunction< value_t >& function,
const FaceDoFFunction< value_t >& function,
const std::shared_ptr< PrimitiveStorage >& storage,
const uint_t& level,
bool write2D,
......
......@@ -33,7 +33,7 @@ class VTKDGDoFWriter
private:
template < typename value_t >
static void writeScalarFunction( std::ostream& output,
const DGFunction< value_t >& function,
const FaceDoFFunction< value_t >& function,
const std::shared_ptr< PrimitiveStorage >& storage,
const uint_t& level,
bool write2D,
......
......@@ -24,10 +24,10 @@
#include "hyteg/Levelinfo.hpp"
#include "hyteg/celldofspace/CellDoFIndexing.hpp"
#include "hyteg/communication/Syncing.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/edgedofspace/EdgeDoFFunction.hpp"
#include "hyteg/edgedofspace/EdgeDoFIndexing.hpp"
#include "hyteg/edgedofspace/EdgeDoFMacroCell.hpp"
#include "hyteg/facedofspace/FaceDoFFunction.hpp"
#include "hyteg/facedofspace/FaceDoFIndexing.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/VertexDoFFunction.hpp"
......@@ -100,7 +100,7 @@ void VTKOutput::add( const GenericFunction< value_t >& function )
break;
case functionTraits::DG_FUNCTION:
matchFound = tryUnwrapAndAdd< FunctionWrapper< DGFunction< value_t > > >( function );
matchFound = tryUnwrapAndAdd< FunctionWrapper< FaceDoFFunction< value_t > > >( function );
break;
default:
......
......@@ -28,8 +28,8 @@
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/P2P1TaylorHoodFunction.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/edgedofspace/EdgeDoFFunction.hpp"
#include "hyteg/facedofspace/FaceDoFFunction.hpp"
#include "hyteg/functions/BlockFunction.hpp"
#include "hyteg/functions/FunctionMultiStore.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
......@@ -107,7 +107,7 @@ class VTKOutput
}
template < typename value_t >
inline void add( const DGFunction< value_t >& function )
inline void add( const FaceDoFFunction< value_t >& function )
{
dgFunctions_.push_back( function );
}
......@@ -233,7 +233,7 @@ class VTKOutput
FunctionMultiStore< P2VectorFunction > p2VecFunctions_;
FunctionMultiStore< EdgeDoFFunction > edgeDoFFunctions_;
FunctionMultiStore< DGFunction > dgFunctions_;
FunctionMultiStore< FaceDoFFunction > dgFunctions_;
std::shared_ptr< PrimitiveStorage > storage_;
......
/*
* Copyright (c) 2017-2019 Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "hyteg/memory/FunctionMemory.hpp"
#include "hyteg/primitives/all.hpp"
#include "DGMemory.hpp"
namespace hyteg {
template< typename ValueType >
class VertexDGFunctionMemoryDataHandling : public FunctionMemoryDataHandling< FunctionMemory < ValueType >, Vertex >
{
public:
VertexDGFunctionMemoryDataHandling( const uint_t & minLevel, const uint_t & maxLevel )
: minLevel_( minLevel ),
maxLevel_( maxLevel )
{}
std::shared_ptr< FunctionMemory< ValueType > > initialize( const Vertex * const vertex ) const override;
private:
uint_t minLevel_;
uint_t maxLevel_;
};
template< typename ValueType >
class EdgeDGFunctionMemoryDataHandling : public FunctionMemoryDataHandling< FunctionMemory < ValueType >, Edge >
{
public:
EdgeDGFunctionMemoryDataHandling( const uint_t & minLevel, const uint_t & maxLevel )
: minLevel_( minLevel ),
maxLevel_( maxLevel )
{}
std::shared_ptr< FunctionMemory< ValueType > > initialize( const Edge * const vertex ) const override;
private:
uint_t minLevel_;
uint_t maxLevel_;
};
template< typename ValueType >
class FaceDGFunctionMemoryDataHandling : public FunctionMemoryDataHandling< FunctionMemory < ValueType >, Face >
{
public:
FaceDGFunctionMemoryDataHandling( const uint_t & minLevel, const uint_t & maxLevel )
: minLevel_( minLevel ),
maxLevel_( maxLevel )
{}
std::shared_ptr<FunctionMemory<ValueType>> initialize( const Face *const face) const override;
private:
uint_t minLevel_;
uint_t maxLevel_;
};
template< typename ValueType >
std::shared_ptr< FunctionMemory< ValueType > > VertexDGFunctionMemoryDataHandling< ValueType >::initialize( const Vertex * const vertex ) const
{
return std::make_shared< FunctionMemory< ValueType > >( DGVertexFunctionMemorySize, *vertex, minLevel_, maxLevel_ );
}
template< typename ValueType >
std::shared_ptr< FunctionMemory< ValueType > > EdgeDGFunctionMemoryDataHandling< ValueType >::initialize( const Edge * const edge ) const
{
return std::make_shared< FunctionMemory< ValueType > >( DGEdgeFunctionMemorySize, *edge, minLevel_, maxLevel_ );
}
template< typename ValueType >
std::shared_ptr< FunctionMemory< ValueType > > FaceDGFunctionMemoryDataHandling< ValueType >::initialize( const Face * const face ) const
{
return std::make_shared< FunctionMemory< ValueType > >( DGFaceFunctionMemorySize, *face, minLevel_, maxLevel_ );
}
}
This diff is collapsed.
/*
* Copyright (c) 2017-2019 Boerge Struempfel, Daniel Drzisga, Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/dgfunctionspace/DGFunctionMacroEdge.hpp"
#include "hyteg/dgfunctionspace/DGFunctionMacroFace.hpp"
#include "hyteg/dgfunctionspace/DGFunctionMacroVertex.hpp"
namespace hyteg {
template < typename ValueType >
void DGFunction< ValueType >::projectP1( P1Function< ValueType >& src, uint_t level, DoFType flag, UpdateType updateType )
{
this->startTiming( "projectP1" );
src.template startCommunication< Edge, Vertex >( level );
src.template startCommunication< Face, Edge >( level );
src.template endCommunication< Edge, Vertex >( level );
for ( auto& it : this->getStorage()->getVertices() )
{
Vertex& vertex = *it.second;
const DoFType vertexBC = this->getBoundaryCondition().getBoundaryType( vertex.getMeshBoundaryFlag() );
if ( testFlag( vertexBC, flag ) )
{
dgfunction::macrovertex::projectP1< ValueType >(
level, vertex, this->getStorage(), src.getVertexDataID(), this->getVertexDataID(), updateType );
}
}
this->template startCommunication< Vertex, Edge >( level );
src.template endCommunication< Face, Edge >( level );
for ( auto& it : this->getStorage()->getEdges() )
{
Edge& edge = *it.second;
const DoFType edgeBC = this->getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if ( testFlag( edgeBC, flag ) )
{
dgfunction::macroedge::projectP1< ValueType >(
level, edge, this->getStorage(), src.getEdgeDataID(), this->getEdgeDataID(), updateType );
}
}
this->template endCommunication< Vertex, Edge >( level );
this->template startCommunication< Edge, Face >( level );
for ( auto& it : this->getStorage()->getFaces() )
{
Face& face = *it.second;
const DoFType faceBC = this->getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
if ( testFlag( faceBC, flag ) )
{
dgfunction::macroface::projectP1< ValueType >(
level, face, this->getStorage(), src.getFaceDataID(), this->getFaceDataID(), updateType );
}
}
this->template endCommunication< Edge, Face >( level );
this->stopTiming( "projectP1" );
}
template class DGFunction< real_t >;
template class DGFunction< int32_t >;
template class DGFunction< int64_t >;
} // namespace hyteg
This diff is collapsed.
/*
* Copyright (c) 2017-2019 Daniel Drzisga, Dominik Thoennes, Marcus Mohr, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
#include "hyteg/facedofspace/FaceDoFIndexing.hpp"
#include "hyteg/p1functionspace/VertexDoFIndexing.hpp"
namespace hyteg {
namespace dgfunction {
namespace macroedge {
template < typename ValueType >
inline void upwind( const uint_t& Level,
Edge& edge,
const std::shared_ptr< PrimitiveStorage >& storage,
const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& srcId,
const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& dstId,
const std::array< PrimitiveDataID< FunctionMemory< ValueType >, Edge >, 2 >& velocityIds,
UpdateType updateType )
{
auto src = edge.getData( srcId )->getPointer( Level );
auto dst = edge.getData( dstId )->getPointer( Level );
auto u = edge.getData( velocityIds[0] )->getPointer( Level );
auto v = edge.getData( velocityIds[1] )->getPointer( Level );
size_t rowsize = levelinfo::num_microvertices_per_edge( Level );
ValueType tmp;
Point2D u_0, u_1, u_2;
real_t un_0, un_1, un_2;
real_t c_up_0, c_up_1, c_up_2;
// first face (south)
{
Face* face = storage->getFace( edge.neighborFaces()[0] );
real_t faceArea = std::pow( 4.0, -walberla::real_c( Level ) ) * face->area;
real_t faceAreaInv = 1.0 / faceArea;
auto oppositeVertex = face->get_vertex_opposite_to_edge( edge.getID() );
uint_t v0 = face->vertex_index( edge.getVertexID0() );
uint_t v1 = face->vertex_index( edge.getVertexID1() );
uint_t v2 = face->vertex_index( oppositeVertex );
// get edge directions
auto d0 = ( face->coords[v1] - face->coords[v0] ) / walberla::real_c( rowsize - 1 );
auto d1 = ( face->coords[v2] - face->coords[v1] ) / walberla::real_c( rowsize - 1 );
auto d2 = ( face->coords[v0] - face->coords[v2] ) / walberla::real_c( rowsize - 1 );
// compute edge lengths
real_t d0Length = d0.norm();
real_t d1Length = d1.norm();
real_t d2Length = d2.norm();
// compute normals
auto n_0 = d0.normal2D() / d0Length;
auto n_1 = d1.normal2D() / d1Length;
auto n_2 = d2.normal2D() / d2Length;
real_t faceOrientation = math::faceOrientation2D( face->coords[v0], face->coords[v1], face->coords[v2] );
n_0 *= faceOrientation;
n_1 *= faceOrientation;
n_2 *= faceOrientation;
for ( uint_t i = 1; i < rowsize - 2; ++i )
{
u_0[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] );
u_0[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] );
u_1[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_SE )] );
u_1[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_SE )] );
u_2[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_SE )] );
u_2[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_SE )] );
// CENTER <-- CELL_GRAY_SE (i)
// NORTH <-- CELL_GRAY_NE (i)
// WEST <-- CELL_BLUE_SE (i)
// EAST <-- CELL_BLUE_SE (i+1)
un_0 = d0Length * u_0.dot( n_0 );
un_1 = d1Length * u_1.dot( n_1 );
un_2 = d2Length * u_2.dot( n_2 );
if ( un_0 >= 0 )
{
c_up_0 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )];
}
else
{
c_up_0 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )];
}
if ( un_1 >= 0 )
{
c_up_1 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )];
}
else
{
c_up_1 = src[facedof::macroedge::indexFaceFromVertex( Level, i + 1, stencilDirection::CELL_BLUE_SE )];
}
if ( un_2 >= 0 )
{
c_up_2 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )];
}
else
{
c_up_2 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_BLUE_SE )];
}
tmp = un_0 * c_up_0 + un_1 * c_up_1 + un_2 * c_up_2;
tmp *= faceAreaInv;
if ( updateType == Replace )
{
dst[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )] = tmp;
}
else if ( updateType == Add )
{
dst[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )] += tmp;
}
}
}
// second face (north)
if ( edge.getNumNeighborFaces() == 2 )
{
Face* face = storage->getFace( edge.neighborFaces()[1] );
real_t faceArea = std::pow( 4.0, -walberla::real_c( Level ) ) * face->area;
real_t faceAreaInv = 1.0 / faceArea;
auto oppositeVertex = face->get_vertex_opposite_to_edge( edge.getID() );
uint_t v0 = face->vertex_index( edge.getVertexID0() );
uint_t v1 = face->vertex_index( edge.getVertexID1() );
uint_t v2 = face->vertex_index( oppositeVertex );
// get edge directions
auto d0 = ( face->coords[v1] - face->coords[v0] ) / walberla::real_c( rowsize - 1 );
auto d1 = ( face->coords[v2] - face->coords[v1] ) / walberla::real_c( rowsize - 1 );
auto d2 = ( face->coords[v0] - face->coords[v2] ) / walberla::real_c( rowsize - 1 );
// compute edge lengths
real_t d0Length = d0.norm();
real_t d1Length = d1.norm();
real_t d2Length = d2.norm();
// compute normals
auto n_0 = d0.normal2D() / d0Length;
auto n_1 = d1.normal2D() / d1Length;
auto n_2 = d2.normal2D() / d2Length;
real_t faceOrientation = math::faceOrientation2D( face->coords[v0], face->coords[v1], face->coords[v2] );
n_0 *= faceOrientation;
n_1 *= faceOrientation;
n_2 *= faceOrientation;
for ( uint_t i = 1; i < rowsize - 2; ++i )
{
u_0[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] );
u_0[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] );
u_1[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_N )] );
u_1[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_E )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_N )] );
u_2[0] = 0.5 * ( u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
u[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_N )] );
u_2[1] = 0.5 * ( v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_C )] +
v[vertexdof::macroedge::indexFromVertex( Level, i, stencilDirection::VERTEX_N )] );
// CENTER <-- CELL_GRAY_NE (i)
// SOUTH <-- CELL_GRAY_SE (i)
// WEST <-- CELL_BLUE_NW (i+1)
// EAST <-- CELL_BLUE_NW (i)
un_0 = d0Length * u_0.dot( n_0 );
un_1 = d1Length * u_1.dot( n_1 );
un_2 = d2Length * u_2.dot( n_2 );
if ( un_0 >= 0 )
{
c_up_0 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )];
}
else
{
c_up_0 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_SE )];
}
if ( un_1 >= 0 )
{
c_up_1 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )];
}
else
{
c_up_1 = src[facedof::macroedge::indexFaceFromVertex( Level, i + 1, stencilDirection::CELL_BLUE_NW )];
}
if ( un_2 >= 0 )
{
c_up_2 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )];
}
else
{
c_up_2 = src[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_BLUE_NW )];
}
tmp = un_0 * c_up_0 + un_1 * c_up_1 + un_2 * c_up_2;
tmp *= faceAreaInv;
if ( updateType == Replace )
{
dst[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )] = tmp;
}
else if ( updateType == Add )
{
dst[facedof::macroedge::indexFaceFromVertex( Level, i, stencilDirection::CELL_GRAY_NE )] += tmp;
}
}
}
}
template < typename ValueType >
inline void projectP1( const uint_t& Level,
Edge& edge,
const std::shared_ptr< PrimitiveStorage >& storage,
const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& srcId,
const PrimitiveDataID< FunctionMemory< ValueType >, Edge >& dstId,
UpdateType updateType )
{
auto src = edge.getData( srcId )->getPointer( Level );
auto dst = edge.getData( dstId )->getPointer( Level );
size_t rowsize = levelinfo::num_microvertices_per_edge( Level );
ValueType tmp;
// first face (south)
{
for ( uint_t i = 1; i < rowsize - 2; ++i )