Commit 1f53050e authored by Dominik Thoennes's avatar Dominik Thoennes
Browse files

Merge branch 'p2_restrict' into 'master'

Add P2 Functions

See merge request terraneo/tinyhhg!135
parents 5e0a907b c28c7c08
Pipeline #8233 passed with stages
in 76 minutes and 37 seconds
---
Language: Cpp
AccessModifierOffset: -2
AlignAfterOpenBracket: Align
AlignConsecutiveAssignments: true
AlignConsecutiveDeclarations: true
AlignEscapedNewlines: Left
AlignOperands: true
AlignTrailingComments: true
AllowAllParametersOfDeclarationOnNextLine: false
AllowShortBlocksOnASingleLine: true
AllowShortCaseLabelsOnASingleLine: false
AllowShortFunctionsOnASingleLine: Inline
AllowShortIfStatementsOnASingleLine: false
AllowShortLoopsOnASingleLine: false
AlwaysBreakAfterDefinitionReturnType: None
AlwaysBreakAfterReturnType: None
AlwaysBreakBeforeMultilineStrings: false
AlwaysBreakTemplateDeclarations: true
BinPackArguments: false
BinPackParameters: false
BraceWrapping:
AfterClass: true
AfterControlStatement: true
AfterEnum: true
AfterFunction: true
AfterNamespace: false
AfterObjCDeclaration: true
AfterStruct: true
AfterUnion: true
AfterExternBlock: true
BeforeCatch: false
BeforeElse: false
IndentBraces: false
SplitEmptyFunction: false
SplitEmptyRecord: false
SplitEmptyNamespace: false
BreakBeforeBinaryOperators: None
BreakBeforeBraces: Custom
BreakBeforeInheritanceComma: false
BreakBeforeTernaryOperators: false
BreakConstructorInitializersBeforeComma: false
BreakConstructorInitializers: BeforeComma
BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: false
ColumnLimit: 130
CommentPragmas: '^ IWYU pragma:'
CompactNamespaces: false
ConstructorInitializerAllOnOneLineOrOnePerLine: false
ConstructorInitializerIndentWidth: 0
ContinuationIndentWidth: 4
Cpp11BracedListStyle: true
DerivePointerAlignment: false
DisableFormat: false
ExperimentalAutoDetectBinPacking: false
FixNamespaceComments: true
ForEachMacros:
- foreach
- Q_FOREACH
- BOOST_FOREACH
IncludeBlocks: Regroup
IncludeCategories:
- Regex: '^<'
Priority: 1
- Regex: '^"core/'
Priority: 2
- Regex: '^"tinyhhg_core'
Priority: 3
- Regex: '.*'
Priority: 4
IncludeIsMainRegex: '(Test)?$'
IndentCaseLabels: false
IndentPPDirectives: None
IndentWidth: 3
IndentWrappedFunctionNames: true
JavaScriptQuotes: Leave
JavaScriptWrapImports: true
KeepEmptyLinesAtTheStartOfBlocks: false
MacroBlockBegin: ''
MacroBlockEnd: ''
MaxEmptyLinesToKeep: 1
NamespaceIndentation: None
ObjCBlockIndentWidth: 2
ObjCSpaceAfterProperty: false
ObjCSpaceBeforeProtocolList: true
PenaltyBreakAssignment: 2
PenaltyBreakBeforeFirstCallParameter: 19
PenaltyBreakComment: 300
PenaltyBreakFirstLessLess: 120
PenaltyBreakString: 1000
PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 60
PointerAlignment: Left
RawStringFormats:
- Delimiter: pb
Language: TextProto
BasedOnStyle: google
ReflowComments: false
SortIncludes: true
SortUsingDeclarations: true
SpaceAfterCStyleCast: true
SpaceAfterTemplateKeyword: true
SpaceBeforeAssignmentOperators: true
SpaceBeforeParens: Never
SpaceInEmptyParentheses: false
SpacesBeforeTrailingComments: 1
SpacesInAngles: true
SpacesInContainerLiterals: true
SpacesInCStyleCastParentheses: false
SpacesInParentheses: true
SpacesInSquareBrackets: false
Standard: Cpp11
TabWidth: 8
UseTab: Never
...
......@@ -72,6 +72,10 @@ waLBerla_add_executable( NAME gmg_P2
FILES gmg_P2.cpp
DEPENDS tinyhhg_core)
waLBerla_add_executable( NAME gmg_P2_h_refinement
FILES gmg_P2_h_refinement.cpp
DEPENDS tinyhhg_core)
if( HHG_BUILD_WITH_PETSC )
waLBerla_add_executable( NAME stokes_mini_petsc
FILES stokes_mini_petsc.cpp
......
#include "core/timing/Timer.h"
#include "core/math/Random.h"
#include "core/Environment.h"
#include "core/config/Config.h"
#include "core/logging/Logging.h"
#include "tinyhhg_core/mesh/MeshInfo.hpp"
#include "tinyhhg_core/primitivestorage/PrimitiveStorage.hpp"
#include "tinyhhg_core/primitivestorage/SetupPrimitiveStorage.hpp"
#include "tinyhhg_core/primitivestorage/loadbalancing/SimpleBalancer.hpp"
#include "tinyhhg_core/p2functionspace/P2Function.hpp"
#include "tinyhhg_core/p2functionspace/P2ConstantOperator.hpp"
#include "tinyhhg_core/solvers/gmultigrid.hpp"
#include "tinyhhg_core/solvers/cgsolver.hpp"
#include "tinyhhg_core/format.hpp"
#include "tinyhhg_core/vtkwriter.hpp"
#include "tinyhhg_core/misc/ExactStencilWeights.hpp"
using walberla::real_t;
using walberla::uint_t;
using walberla::uint_c;
using namespace hhg;
int main(int argc, char* argv[]) {
walberla::Environment walberlaEnv(argc, argv);
walberla::logging::Logging::instance()->setLogLevel(walberla::logging::Logging::PROGRESS);
walberla::MPIManager::instance()->useWorldComm();
walberla::shared_ptr<walberla::config::Config> cfg(new walberla::config::Config);
cfg->readParameterFile("../data/param/gmg_P2.prm");
walberla::Config::BlockHandle parameters = cfg->getOneBlock("Parameters");
const uint_t minLevel = parameters.getParameter<uint_t>("minLevel");
const uint_t maxLevel = parameters.getParameter<uint_t>("maxLevel");
const uint_t max_outer_iter = parameters.getParameter<uint_t>("max_outer_iter");
const uint_t max_cg_iter = parameters.getParameter<uint_t>("max_cg_iter");
const real_t mg_tolerance = parameters.getParameter<real_t>("mg_tolerance");
const real_t coarse_tolerance = parameters.getParameter<real_t>("coarse_tolerance");
MeshInfo meshInfo = MeshInfo::fromGmshFile( parameters.getParameter<std::string>("mesh") );
SetupPrimitiveStorage setupStorage( meshInfo, uint_c ( walberla::mpi::MPIManager::instance()->numProcesses() ) );
hhg::loadbalancing::roundRobin( setupStorage );
std::shared_ptr< walberla::WcTimingTree > timingTree( new walberla::WcTimingTree() );
std::shared_ptr<PrimitiveStorage> storage = std::make_shared<PrimitiveStorage>(setupStorage, timingTree);
hhg::P2Function< real_t > r("r", storage, minLevel, maxLevel);
hhg::P2Function< real_t > f("f", storage, minLevel, maxLevel);
hhg::P2Function< real_t > u("u", storage, minLevel, maxLevel);
hhg::P2Function< real_t > Lu("Lu", storage, minLevel, maxLevel);
hhg::P2Function< real_t > tmp("tmp", storage, minLevel, maxLevel);
hhg::P2Function< real_t > u_exact("u_exact", storage, minLevel, maxLevel);
hhg::P2Function< real_t > err("err", storage, minLevel, maxLevel);
hhg::P2Function< real_t > npoints_helper("npoints_helper", storage, minLevel, maxLevel);
std::function<real_t(const hhg::Point3D &)> exact = [](const hhg::Point3D &x) { return sin(x[0])*sinh(x[1]); };
std::function<real_t(const hhg::Point3D &)> rhs = [](const hhg::Point3D & ) { return 0; };
std::function<real_t(const hhg::Point3D &)> zero = [](const hhg::Point3D & ) { return 0.0; };
std::function<real_t(const hhg::Point3D &)> ones = [](const hhg::Point3D & ) { return 1.0; };
walberla::math::seedRandomGenerator(0);
std::function<real_t(const Point3D &)> rand = [](const Point3D &) { return walberla::math::realRandom(0.0, 20.0); };
WALBERLA_LOG_INFO_ON_ROOT("Interpolating u");
u.interpolate(rand, maxLevel, hhg::Inner);
u.interpolate(exact, maxLevel, hhg::DirichletBoundary);
u_exact.interpolate(exact, maxLevel);
// WALBERLA_LOG_INFO_ON_ROOT("Interpolating and integrating rhs");
// npoints_helper.interpolate(rhs, maxLevel);
// M.apply(npoints_helper, f, maxLevel, hhg::All);
WALBERLA_LOG_INFO_ON_ROOT("Setting up stiffness operator");
auto start = walberla::timing::getWcTime();
hhg::P2ConstantLaplaceOperator L(storage, minLevel, maxLevel );
auto end = walberla::timing::getWcTime();
real_t setupTime = end - start;
npoints_helper.interpolate(ones, maxLevel);
real_t npoints = npoints_helper.dot(npoints_helper, maxLevel);
typedef hhg::CGSolver<hhg::P2Function<real_t>, hhg::P2ConstantLaplaceOperator> CoarseSolver;
auto coarseLaplaceSolver = std::make_shared<CoarseSolver>(storage, minLevel, minLevel);
typedef GMultigridSolver<hhg::P2Function<real_t>, hhg::P2ConstantLaplaceOperator, CoarseSolver> LaplaceSover;
LaplaceSover laplaceSolver(storage, coarseLaplaceSolver, minLevel, maxLevel);
if (parameters.getParameter<bool>("useExactWeights")) {
WALBERLA_LOG_INFO("WARNING: works only on tri_1el mesh");
auto weights = hhg::stencilWeights::tri_1el();
for (uint_t i = minLevel; i <= maxLevel; ++i) {
real_t *vToV = storage->getFace(PrimitiveID(6))->getData(L.getVertexToVertexOpr().getFaceStencilID())->getPointer(i);
for (uint_t j = 0; j < 7; ++j) {
vToV[j] = weights.vertexToVertexStencil[j];
}
real_t *eToV = storage->getFace(PrimitiveID(6))->getData(L.getEdgeToVertexOpr().getFaceStencilID())->getPointer(i);
for (uint_t j = 0; j < 12; ++j) {
eToV[j] = weights.edgeToVertexStencil[j];
}
real_t *vToE = storage->getFace(PrimitiveID(6))->getData(L.getVertexToEdgeOpr().getFaceStencilID())->getPointer(i);
for (uint_t j = 0; j < 12; ++j) {
vToE[j] = weights.vertexToEdgeStencil[j];
}
real_t *eToE = storage->getFace(PrimitiveID(6))->getData(L.getEdgeToEdgeOpr().getFaceStencilID())->getPointer(i);
for (uint_t j = 0; j < 15; ++j) {
eToE[j] = weights.edgeToEdgeStencil[j];
}
}
}
WALBERLA_LOG_INFO_ON_ROOT("Starting V cycles");
WALBERLA_LOG_INFO_ON_ROOT(hhg::format("%6s|%10s|%10s|%10s|%10s|%10s","iter","abs_res","rel_res","conv","L2-error","Time"));
real_t rel_res = 1.0;
L.apply(u, Lu, maxLevel, hhg::Inner);
r.assign({1.0, -1.0}, {&f, &Lu}, maxLevel, hhg::Inner);
real_t begin_res = std::sqrt(r.dot(r, maxLevel, hhg::Inner));
real_t abs_res_old = begin_res;
err.assign({1.0, -1.0}, {&u, &u_exact}, maxLevel);
real_t discr_l2_err = std::sqrt(err.dot(err, maxLevel) / npoints);
//WALBERLA_LOG_INFO_ON_ROOT(fmt::format("{:3d} {:e} {:e} {:e} {:e} -", 0, begin_res, rel_res, begin_res/abs_res_old, discr_l2_err));
WALBERLA_LOG_INFO_ON_ROOT(hhg::format("%6d|%10.3e|%10.3e|%10.3e|%10.3e|%10.3e", 0, begin_res, rel_res, begin_res/abs_res_old, discr_l2_err,0))
real_t solveTime = real_c(0.0);
real_t averageConvergenceRate = real_c(0.0);
const uint_t convergenceStartIter = 3;
uint_t i = 0;
for (; i < max_outer_iter; ++i)
{
start = walberla::timing::getWcTime();
// Apply P2 geometric multigrid solver
laplaceSolver.solve(L, u, f, r, maxLevel, coarse_tolerance, max_cg_iter, hhg::Inner, LaplaceSover::CycleType::VCYCLE, false);
end = walberla::timing::getWcTime();
L.apply(u, Lu, maxLevel, hhg::Inner);
r.assign({1.0, -1.0}, { &f, &Lu }, maxLevel, hhg::Inner);
real_t abs_res = std::sqrt(r.dot(r, maxLevel, hhg::Inner));
rel_res = abs_res / begin_res;
err.assign({1.0, -1.0}, { &u, &u_exact }, maxLevel);
discr_l2_err = std::sqrt(err.dot(err, maxLevel) / npoints);
//WALBERLA_LOG_INFO_ON_ROOT(fmt::format("{:3d} {:e} {:e} {:e} {:e} {:e}", i+1, abs_res, rel_res, abs_res/abs_res_old, discr_l2_err, end-start));
WALBERLA_LOG_INFO_ON_ROOT(hhg::format("%6d|%10.3e|%10.3e|%10.3e|%10.3e|%10.3e", i+1, abs_res, rel_res, abs_res/abs_res_old, discr_l2_err,end - start))
solveTime += end-start;
if (i >= convergenceStartIter) {
averageConvergenceRate += abs_res/abs_res_old;
}
abs_res_old = abs_res;
if (rel_res < mg_tolerance)
{
break;
}
}
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "MinLevel: " << minLevel);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "MaxLevel: " << maxLevel);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "DoFs on MaxLevel: " << (uint_t) npoints);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "Setup time: " << std::defaultfloat << setupTime);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "Solve time: " << std::defaultfloat << solveTime);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "Time to solution: " << std::defaultfloat << setupTime + solveTime);
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "Avg. convergence rate: " << std::scientific << averageConvergenceRate / real_c(i-convergenceStartIter));
WALBERLA_LOG_INFO_ON_ROOT(std::setw(25) << "L^2 error: " << std::scientific << discr_l2_err);
if (parameters.getParameter<bool>("vtkOutput")) {
VTKOutput vtkOutput("../output", "gmg_P2_h_refinement");
vtkOutput.add(&u);
vtkOutput.add(&u_exact);
vtkOutput.add(&f);
vtkOutput.add(&r);
vtkOutput.add(&err);
vtkOutput.add(&npoints_helper);
vtkOutput.write(maxLevel);
}
if (parameters.getParameter<bool>("printTiming")) {
walberla::WcTimingTree tt = timingTree->getReduced();
WALBERLA_LOG_INFO_ON_ROOT(tt);
}
return 0;
}
......@@ -14,7 +14,7 @@ $Elements
2 15 2 1 2 2
3 15 2 1 3 3
4 15 2 1 4 4
5 1 2 1 1 1 2
5 1 2 0 1 1 2
6 1 2 1 2 2 4
7 1 2 1 2 4 1
8 1 2 1 3 2 3
......
Parameters
{
mesh ../data/meshes/tri_1el.msh;
mesh ../data/meshes/quad_16el.msh;
level 3;
maxiter 10000;
tolerance 1e-10;
......
......@@ -2,10 +2,12 @@ Parameters
{
mesh ../data/meshes/tri_1el.msh;
minLevel 2;
maxLevel 7;
max_outer_iter 100;
maxLevel 8;
max_outer_iter 20;
max_cg_iter 10000;
mg_tolerance 1e-10;
coarse_tolerance 1e-10;
mg_tolerance 1e-16;
coarse_tolerance 1e-16;
vtkOutput false;
printTiming false;
useExactWeights false;
}
Parameters
{
mesh ../../data/meshes/quad_4el.msh;
level 2;
maxiter 100;
vtkOutput false;
printTiming false;
}
Parameters
{
mesh ../../data/meshes/tri_1el.msh;
level 2;
maxiter 100;
vtkOutput false;
printTiming false;
}
......@@ -212,6 +212,10 @@ constexpr std::array<stencilDirection ,5> neighborsFromHorizontalEdge =
sD::EDGE_DI_N, sD::EDGE_VE_NW
}};
constexpr std::array<stencilDirection ,4> neighborsFromHorizontalEdgeWithoutCenter =
{{ sD::EDGE_DI_S, sD::EDGE_VE_SE,
sD::EDGE_DI_N, sD::EDGE_VE_NW
}};
inline constexpr uint_t indexFromDiagonalEdge( const uint_t & level, const uint_t & col, const uint_t & row, const stencilDirection & dir )
{
......@@ -239,6 +243,11 @@ constexpr std::array<stencilDirection ,5> neighborsFromDiagonalEdge =
sD::EDGE_HO_N, sD::EDGE_VE_W
}};
constexpr std::array<stencilDirection ,4> neighborsFromDiagonalEdgeWithoutCenter =
{{ sD::EDGE_HO_S, sD::EDGE_VE_E,
sD::EDGE_HO_N, sD::EDGE_VE_W
}};
inline constexpr uint_t indexFromVerticalEdge( const uint_t & level, const uint_t & col, const uint_t & row, const stencilDirection & dir )
{
switch( dir )
......@@ -265,6 +274,11 @@ constexpr std::array<stencilDirection ,5> neighborsFromVerticalEdge =
sD::EDGE_HO_NW, sD::EDGE_DI_W
}};
constexpr std::array<stencilDirection ,4> neighborsFromVerticalEdgeWithoutCenter =
{{ sD::EDGE_HO_SE, sD::EDGE_DI_E,
sD::EDGE_HO_NW, sD::EDGE_DI_W
}};
inline constexpr uint_t indexFromVertex( const uint_t & level, const uint_t & col, const uint_t & row, const stencilDirection & dir )
{
// first neighbor == south
......@@ -470,5 +484,22 @@ constexpr inline uint_t stencilIndexFromVerticalEdge(const stencilDirection dir)
}
}
inline bool isHorizontalEdgeOnBoundary(const uint_t level, const hhg::indexing::Index& idx){
/// level is only needed in the diagonal case
WALBERLA_UNUSED( level );
return ( idx.row() == 0 );
}
inline bool isVerticalEdgeOnBoundary(const uint_t level, const hhg::indexing::Index& idx){
/// level is only needed in the diagonal case
WALBERLA_UNUSED( level );
return ( idx.col() == 0 );
}
inline bool isDiagonalEdgeOnBoundary(const uint_t level, const hhg::indexing::Index& idx){
return ( (idx.col() + idx.row()) == (hhg::levelinfo::num_microedges_per_edge( level ) - 1) );
}
} // namespace edgedof
} // namespace hhg
#pragma once
#include <vector>
#include "tinyhhg_core/primitives/Edge.hpp"
#include "tinyhhg_core/primitives/Face.hpp"
#include "tinyhhg_core/levelinfo.hpp"
......@@ -8,6 +9,8 @@
#include "tinyhhg_core/FunctionMemory.hpp"
#include "tinyhhg_core/StencilMemory.hpp"
#include "core/math/KahanSummation.h"
namespace hhg {
namespace edgedof {
namespace macroedge {
......@@ -116,7 +119,7 @@ inline real_t dot( const uint_t & Level, Edge & edge,
auto lhsData = edge.getData( lhsId )->getPointer( Level );
auto rhsData = edge.getData( rhsId )->getPointer( Level );
real_t scalarProduct = real_c( 0 );
walberla::math::KahanAccumulator< ValueType > scalarProduct;
for ( const auto & it : edgedof::macroedge::Iterator( Level ) )
{
......@@ -124,7 +127,7 @@ inline real_t dot( const uint_t & Level, Edge & edge,
scalarProduct += lhsData[ idx ] * rhsData[ idx ];
}
return scalarProduct;
return scalarProduct.get();
}
......@@ -186,7 +189,7 @@ inline void apply(const uint_t & Level, Edge &edge,
template< typename ValueType >
inline void printFunctionMemory(const uint_t & Level, Edge& edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId){
inline void printFunctionMemory(const uint_t & Level, const Edge& edge, const PrimitiveDataID<FunctionMemory< ValueType >, Edge> &dstId){
ValueType* edgeMemory = edge.getData(dstId)->getPointer( Level );
using namespace std;
cout << setfill('=') << setw(100) << "" << endl;
......
......@@ -221,7 +221,7 @@ inline real_t dot( const uint_t & Level, Face & face,
auto lhsData = face.getData( lhsId )->getPointer( Level );
auto rhsData = face.getData( rhsId )->getPointer( Level );
real_t scalarProduct = real_c( 0 );
walberla::math::KahanAccumulator< ValueType > scalarProduct;
for ( const auto & it : edgedof::macroface::Iterator( Level, 0 ) )
{
......@@ -247,7 +247,7 @@ inline real_t dot( const uint_t & Level, Face & face,
}
}
return scalarProduct;
return scalarProduct.get();
}
......
#pragma once
namespace hhg{
namespace stencilWeights{
class tri_1el {
public:
real_t vertexToVertexStencil[7];
real_t edgeToVertexStencil[12];
real_t vertexToEdgeStencil[12];
real_t edgeToEdgeStencil[15];
tri_1el() {
vertexToVertexStencil[0] = 1.0 / 3.0;
vertexToVertexStencil[1] = 0.0;
vertexToVertexStencil[2] = 1.0 / 3.0;
vertexToVertexStencil[3] = 4.0;
vertexToVertexStencil[4] = 1.0 / 3.0;
vertexToVertexStencil[5] = 0.0;
vertexToVertexStencil[6] = 1.0 / 3.0;
edgeToVertexStencil[0] = -(1.0 + 1.0 / 3.0);
edgeToVertexStencil[1] = 0.0;
edgeToVertexStencil[2] = -(1.0 + 1.0 / 3.0);
edgeToVertexStencil[3] = 0.0;
edgeToVertexStencil[4] = 0.0;
edgeToVertexStencil[5] = 0.0;
edgeToVertexStencil[6] = -(1.0 + 1.0 / 3.0);
edgeToVertexStencil[7] = 0.0;
edgeToVertexStencil[8] = -(1.0 + 1.0 / 3.0);
edgeToVertexStencil[9] = 0.0;
edgeToVertexStencil[10] = 0.0;
edgeToVertexStencil[11] = 0.0;
vertexToEdgeStencil[0] = - (1.0 + 1.0/3.0);
vertexToEdgeStencil[1] = - (1.0 + 1.0/3.0);
vertexToEdgeStencil[2] = 0.0;
vertexToEdgeStencil[3] = 0.0;
vertexToEdgeStencil[4] = 0.0;
vertexToEdgeStencil[5] = 0.0;
vertexToEdgeStencil[6] = 0.0;
vertexToEdgeStencil[7] = 0.0;
vertexToEdgeStencil[8] = - (1.0 + 1.0/3.0);
vertexToEdgeStencil[9] = 0.0;
vertexToEdgeStencil[10] = - (1.0 + 1.0/3.0);;
vertexToEdgeStencil[11] = 0.0;
edgeToEdgeStencil[0] = (5.0 + 1.0/3.0);
edgeToEdgeStencil[1] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[2] = 0.0;
edgeToEdgeStencil[3] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[4] = 0.0;
edgeToEdgeStencil[5] = (5.0 + 1.0/3.0);
edgeToEdgeStencil[6] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[7] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[8] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[9] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[10] = (5.0 + 1.0/3.0);
edgeToEdgeStencil[11] = 0.0;
edgeToEdgeStencil[12] = - (1.0 + 1.0/3.0);
edgeToEdgeStencil[13] = 0.0;
edgeToEdgeStencil[14] = - (1.0 + 1.0/3.0);
}
};
}/// namespace stencilWeights
}/// namespace hhg
\ No newline at end of file
......@@ -324,8 +324,8 @@ public:
class BorderIterator : public hhg::indexing::FaceBorderIterator
{
public:
BorderIterator( const uint_t & level, const hhg::indexing::FaceBorderDirection & direction, const uint_t & offsetToCenter = 0 ) :
FaceBorderIterator( levelinfo::num_microvertices_per_edge( level ), direction, offsetToCenter )
BorderIterator( const uint_t & level, const hhg::indexing::FaceBorderDirection & direction, const uint_t & offsetToCenter = 0, const uint_t & offsetFromVertices = 0) :
FaceBorderIterator( levelinfo::num_microvertices_per_edge( level ), direction, offsetToCenter, offsetFromVertices )
{}
};
......@@ -555,5 +555,17 @@ constexpr inline uint_t stencilIndexFromBlueFace( const stencilDirection & dir )
}
}
inline bool isVertexOnBoundary(const uint_t level, const hhg::indexing::Index& idx){
if (idx.row() == 0 ){
return true;
} else if( idx.col() == 0 ){
return true;
} else if( (idx.row() + idx.col()) == ( hhg::levelinfo::num_microvertices_per_edge( level ) - 1 ) ){
return true;
} else {
return false;
}