Commit b54b6417 authored by Praetorius, Simon's avatar Praetorius, Simon

Merge branch 'feature/entity_hierarchic_iterator' into 'master'

Feature/entity hierarchic iterator

See merge request !3
parents db80a50b e8152f8a
Pipeline #1448 failed with stages
in 1 minute and 9 seconds
...@@ -34,7 +34,7 @@ namespace Dune ...@@ -34,7 +34,7 @@ namespace Dune
using Super = std::vector<HostGridEntity>; using Super = std::vector<HostGridEntity>;
/// Constructor from std::vector /// Constructor from std::vector
using Super::vector; using std::vector<HostGridEntity>::vector;
/// Return a local geometry of source in target /// Return a local geometry of source in target
static LocalGeometry localGeometry (HostGridEntity const& source, HostGridEntity const& target) static LocalGeometry localGeometry (HostGridEntity const& source, HostGridEntity const& target)
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH
#define DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH
#include <numeric>
#include <stack>
#include <variant>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/std/type_traits.hh>
namespace Dune
{
/** \brief Iterator over all entities of a given codimension and level of a grid.
* \ingroup MultiMesh
*/
template <class HostGrid>
class MultiMeshHierarchicIterator
: public ForwardIteratorFacade<MultiMeshHierarchicIterator<HostGrid>,
typename HostGrid::template Codim<0>::Entity,
typename HostGrid::template Codim<0>::Entity>
{
private:
template <int codim, PartitionIteratorType pitype, class HostGrid_>
friend class MultiMeshIterator;
using HostEntity = typename HostGrid::template Codim<0>::Entity;
using EntityTest = std::function<bool(HostEntity)>;
struct EntityStackEntry
{
template <class Entity>
explicit EntityStackEntry (Entity&& entity)
: it(entity.hbegin(entity.level()+1))
, end(entity.hend(entity.level()+1))
{}
typename HostEntity::HierarchicIterator it;
typename HostEntity::HierarchicIterator end;
friend bool operator== (const EntityStackEntry& lhs, const EntityStackEntry& rhs)
{
return lhs.it == rhs.it;
}
};
class EntityStack
: public std::stack<EntityStackEntry, std::vector<EntityStackEntry>>
{
using Super = std::stack<EntityStackEntry, std::vector<EntityStackEntry>>;
public:
explicit EntityStack (int maxLevel)
{
c.reserve(maxLevel);
}
// return true if all levels >= l are finished, i.e. hierarchic iterators it+1 == end
bool finished (std::size_t l = 0) const
{
bool f = true;
for (auto tmp = c[l].it; f && l < c.size(); ++l) {
tmp = c[l].it;
f = f && (++tmp) == c[l].end;
}
return f;
}
friend bool operator== (const EntityStack& lhs, const EntityStack& rhs)
{
return lhs.size() == rhs.size() && lhs.c == rhs.c;
}
protected:
using Super::c;
};
public:
/// Constructor. Stores a pointer to the entity
template <class Entity, class GridView>
MultiMeshHierarchicIterator (const Entity& entity, const GridView& gridView)
: entity_{&entity}
, contains_{[gv=gridView](const HostEntity& e) { return gv.contains(e); }}
, maxLevel_{gridView.grid().maxLevel()}
, entityStack_{EntityStack{gridView.grid().maxLevel()}}
{
initialIncrement();
}
/// Constructor which create the end iterator
/**
* \param endDummy Here only to distinguish it from the other constructor
*/
MultiMeshHierarchicIterator (bool endDummy)
: entityStack_{0}
{}
/// prefix increment
void increment ()
{
// 1. go up in tree or to next entity on current level until we can go down again
while (!entityStack_.empty()) {
auto& top = entityStack_.top();
++top.it;
if (top.it == top.end) {
entityStack_.pop();
} else {
break;
}
}
// 2. if entityStack is empty, iteration is finished
if (entityStack_.empty())
return;
// 3. go down in tree until leaf entity
auto child = dereference();
for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack_.emplace(child);
assert( entityStack_.size() <= maxLevel_ );
}
assert(contains_(child) && "No valid child element found in gridView");
}
/// dereferencing
decltype(auto) dereference () const
{
if (entityStack_.empty()) {
return *entity_;
} else {
assert(entityStack_.top().it != entityStack_.top().end);
return *entityStack_.top().it;
}
}
/// equality
bool equals (const MultiMeshHierarchicIterator& that) const
{
return entityStack_ == that.entityStack_;
}
protected:
// got to first leaf entity on gridView
void initialIncrement ()
{
// 1. go down in tree until leaf entity
auto child = dereference();
for (; !(contains_(child) || child.isLeaf()); child = dereference()) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack_.emplace(child);
assert( entityStack_.size() <= maxLevel_ );
}
assert(contains_(child) && "No valid child element found in gridView");
}
private:
const HostEntity* entity_ = nullptr;
EntityTest contains_;
int maxLevel_;
EntityStack entityStack_;
};
template <class Entity, class GridView>
inline auto childs(const Entity& entity, const GridView& gridView)
{
using Iterator = MultiMeshHierarchicIterator<typename GridView::Grid>;
return IteratorRange<Iterator>{ Iterator{entity, gridView}, Iterator{true} };
}
} // end namespace Dune
#endif // DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH
...@@ -13,14 +13,10 @@ ...@@ -13,14 +13,10 @@
#include <dune/grid/common/gridenums.hh> #include <dune/grid/common/gridenums.hh>
#include "mmentity.hh" #include "mmentity.hh"
#include "mmhierarchiciterator.hh"
namespace Dune namespace Dune
{ {
namespace tag
{
struct with_gridview {};
}
/** \brief Iterator over all entities of a given codimension and level of a grid. /** \brief Iterator over all entities of a given codimension and level of a grid.
* \ingroup MultiMesh * \ingroup MultiMesh
*/ */
...@@ -55,45 +51,8 @@ namespace Dune ...@@ -55,45 +51,8 @@ namespace Dune
using HostEntity = typename HostGrid::template Codim<0>::Entity; using HostEntity = typename HostGrid::template Codim<0>::Entity;
using EntityTest = std::function<bool(HostEntity)>; using EntityTest = typename MultiMeshHierarchicIterator<HostGrid>::EntityTest;
using EntityStack = typename MultiMeshHierarchicIterator<HostGrid>::EntityStack;
struct EntityStackEntry
{
template <class Entity>
explicit EntityStackEntry (Entity&& entity)
: it(entity.hbegin(entity.level()+1))
, end(entity.hend(entity.level()+1))
{}
typename HostEntity::HierarchicIterator it;
typename HostEntity::HierarchicIterator end;
};
class EntityStack
: public std::stack<EntityStackEntry, std::vector<EntityStackEntry>>
{
using Super = std::stack<EntityStackEntry, std::vector<EntityStackEntry>>;
public:
explicit EntityStack (int maxLevel)
{
c.reserve(maxLevel);
}
// return true if all levels >= l are finished, i.e. hierarchic iterators it+1 == end
bool finished (std::size_t l = 0) const
{
bool f = true;
for (auto tmp = c[l].it; f && l < c.size(); ++l) {
tmp = c[l].it;
f = f && (++tmp) == c[l].end;
}
return f;
}
protected:
using Super::c;
};
public: public:
/// Constructor. Stores a pointer to the grid /// Constructor. Stores a pointer to the grid
...@@ -142,9 +101,10 @@ namespace Dune ...@@ -142,9 +101,10 @@ namespace Dune
: MultiMeshIterator{multiMesh, -1, endDummy} : MultiMeshIterator{multiMesh, -1, endDummy}
{} {}
/// Construct an iterator from n gridViews of possibly different type
template <class... GridViews, template <class... GridViews,
std::enable_if_t<not Std::disjunction<std::is_same<GridViews,bool>...>::value, int> = 0> std::enable_if_t<not Std::disjunction<std::is_same<GridViews,bool>...>::value, int> = 0>
MultiMeshIterator (tag::with_gridview, GridViews const&... gridViews) MultiMeshIterator (GridViews const&... gridViews)
: incrementAllowed_(sizeof...(GridViews), true) : incrementAllowed_(sizeof...(GridViews), true)
, maxLevel_{gridViews.grid().maxLevel()...} , maxLevel_{gridViews.grid().maxLevel()...}
, multiEntity_(sizeof...(GridViews)) , multiEntity_(sizeof...(GridViews))
...@@ -162,11 +122,13 @@ namespace Dune ...@@ -162,11 +122,13 @@ namespace Dune
} }
template <class... GridViews> template <class... GridViews>
MultiMeshIterator (tag::with_gridview, bool endDummy, GridViews const&... gridViews) MultiMeshIterator (bool endDummy, GridViews const&... gridViews)
: macroIterators_{gridViews.grid().levelGridView(0).template end<0,pitype>()...} : macroIterators_{gridViews.grid().levelGridView(0).template end<0,pitype>()...}
{} {}
#if DUNE_HAVE_CXX_VARIANT #if DUNE_HAVE_CXX_VARIANT
/// Construct an iterator from a vector of GridViews, where each GridView can be
/// either a leaf- or a level-gridView, thus we need to use a variant.
template <class... GV> template <class... GV>
MultiMeshIterator (const std::vector<std::variant<GV...>>& gridViews) MultiMeshIterator (const std::vector<std::variant<GV...>>& gridViews)
: incrementAllowed_(gridViews.size(), true) : incrementAllowed_(gridViews.size(), true)
...@@ -258,11 +220,14 @@ namespace Dune ...@@ -258,11 +220,14 @@ namespace Dune
} }
// 3. go down in tree until leaf entity // 3. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(i,child); auto child = dereference(i);
child = dereference(i)) { for (; !(contains_[i](child) || child.isLeaf()); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child); entityStack.emplace(child);
assert( entityStack.size() <= maxLevel_[i] ); assert(entityStack.size() <= maxLevel_[i]);
} }
// assert(contains_[i](child) && "No valid child element found in gridView");
} }
/// Return true, if all stacks with size > stack[i].size are finished /// Return true, if all stacks with size > stack[i].size are finished
...@@ -282,16 +247,19 @@ namespace Dune ...@@ -282,16 +247,19 @@ namespace Dune
auto& macroIt = macroIterators_[i]; auto& macroIt = macroIterators_[i];
auto const& macroEnd = macroEndIterators_[i]; auto const& macroEnd = macroEndIterators_[i];
assert( entityStack.empty() ); assert(entityStack.empty());
if (macroIt == macroEnd) if (macroIt == macroEnd)
return; return;
// 2. go down in tree until leaf entity // 1. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(i,child); auto child = dereference(i);
child = dereference(i)) { for (; !(contains_[i](child) || child.isLeaf()); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child); entityStack.emplace(child);
assert( entityStack.size() <= maxLevel_[i] ); assert(entityStack.size() <= maxLevel_[i]);
} }
// assert(contains_[i](child) && "No valid child element found in gridView");
} }
HostEntity dereference (std::size_t i) const HostEntity dereference (std::size_t i) const
...@@ -305,11 +273,6 @@ namespace Dune ...@@ -305,11 +273,6 @@ namespace Dune
} }
} }
bool entityTest (std::size_t i, HostEntity const& entity) const
{
return contains_[i](entity) || entity.isLeaf() || !entity.isRegular();
}
private: private:
std::vector<bool> incrementAllowed_; std::vector<bool> incrementAllowed_;
std::vector<EntityTest> contains_; std::vector<EntityTest> contains_;
...@@ -322,16 +285,14 @@ namespace Dune ...@@ -322,16 +285,14 @@ namespace Dune
std::vector<EntityStack> entityStacks_; std::vector<EntityStack> entityStacks_;
}; };
template <class> class MultiMesh;
template <class... GridViews> template <class... GridViews>
inline auto multi_elements(GridViews const&... gridViews) inline auto multi_elements(GridViews const&... gridViews)
{ {
using GridView0 = std::tuple_element_t<0,std::tuple<GridViews...>>; using GridView0 = std::tuple_element_t<0,std::tuple<GridViews...>>;
using Iterator = MultiMeshIterator<0,All_Partition,typename GridView0::Grid>; using Iterator = MultiMeshIterator<0,All_Partition,typename GridView0::Grid>;
using Range = IteratorRange<Iterator>; using Range = IteratorRange<Iterator>;
return Range{ Iterator{tag::with_gridview{}, gridViews...}, return Range{ Iterator{gridViews...},
Iterator{tag::with_gridview{}, true, gridViews...} }; Iterator{true, gridViews...} };
} }
} // end namespace Dune } // end namespace Dune
......
...@@ -57,7 +57,7 @@ namespace Dune ...@@ -57,7 +57,7 @@ namespace Dune
/// Map entity to index. /// Map entity to index.
template <int cc> template <int cc>
IndexType index (const typename LeafIndexSet::Traits::template Codim<cc>::Entity& e) const IndexType index (const typename LeafIndexSet::template Codim<cc>::Entity& e) const
{ {
return std::visit([&e](auto const* is) { return is->template index<cc>(e); }, indexSets_); return std::visit([&e](auto const* is) { return is->template index<cc>(e); }, indexSets_);
} }
...@@ -71,7 +71,7 @@ namespace Dune ...@@ -71,7 +71,7 @@ namespace Dune
/// Map a subentity to an index. /// Map a subentity to an index.
template <int cc> template <int cc>
IndexType subIndex (const typename LeafIndexSet::Traits::template Codim<cc>::Entity& e, IndexType subIndex (const typename LeafIndexSet::template Codim<cc>::Entity& e,
int i, unsigned int codim) const int i, unsigned int codim) const
{ {
return std::visit([&e,i,codim](auto const* is) { return is->template subIndex<cc>(e,i,codim); }, indexSets_); return std::visit([&e,i,codim](auto const* is) { return is->template subIndex<cc>(e,i,codim); }, indexSets_);
......
...@@ -60,7 +60,7 @@ namespace Dune ...@@ -60,7 +60,7 @@ namespace Dune
using LevelIterator = typename Partition<All_Partition>::LevelIterator; using LevelIterator = typename Partition<All_Partition>::LevelIterator;
private: private:
friend class HostGrid::GridFamily::Traits::template Codim<cd>::Entity; // friend class HostGrid::GridFamily::Traits::template Codim<cd>::Entity;
}; };
/// The type of view for leaf grid /// The type of view for leaf grid
......
...@@ -16,7 +16,7 @@ namespace Dune ...@@ -16,7 +16,7 @@ namespace Dune
template <class MultiMesh, class... Bases> template <class MultiMesh, class... Bases>
struct MultiBasis struct MultiBasis
{ {
using MultiGridView = typename MultiMesh::LeafGridView; using MultiLeafGridView = typename MultiMesh::LeafGridView;
/// Constructor. Stores a leafgridView of the MultiMesh and LocalView/LocalIndexSet of the bases /// Constructor. Stores a leafgridView of the MultiMesh and LocalView/LocalIndexSet of the bases
MultiBasis (MultiMesh const& mm, Bases&&... bases) MultiBasis (MultiMesh const& mm, Bases&&... bases)
...@@ -27,7 +27,7 @@ namespace Dune ...@@ -27,7 +27,7 @@ namespace Dune
{} {}
/// Return the leafGridView of the MultiMesh /// Return the leafGridView of the MultiMesh
const MultiGridView& gridView () const { return multiGridView_; } const MultiLeafGridView& gridView () const { return multiGridView_; }
/// Return the I'th basis /// Return the I'th basis
template <std::size_t I> template <std::size_t I>
...@@ -63,7 +63,7 @@ namespace Dune ...@@ -63,7 +63,7 @@ namespace Dune
} }
private: private:
MultiGridView multiGridView_; MultiLeafGridView multiGridView_;
std::tuple<Bases...> bases_; std::tuple<Bases...> bases_;
std::tuple<typename Bases::LocalView...> localViews_; std::tuple<typename Bases::LocalView...> localViews_;
......
...@@ -13,6 +13,9 @@ target_link_dune_default_libraries("uggrid") ...@@ -13,6 +13,9 @@ target_link_dune_default_libraries("uggrid")
add_executable("multigridview" multigridview.cc) add_executable("multigridview" multigridview.cc)
target_link_dune_default_libraries("multigridview") target_link_dune_default_libraries("multigridview")
add_executable("hierarchiciterator" hierarchiciterator.cc)
target_link_dune_default_libraries("hierarchiciterator")
find_package(MTL PATHS /opt/development/mtl4) find_package(MTL PATHS /opt/development/mtl4)
if (MTL_FOUND) if (MTL_FOUND)
set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL") set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
...@@ -25,12 +28,18 @@ if (MTL_FOUND) ...@@ -25,12 +28,18 @@ if (MTL_FOUND)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK") list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
endif () endif ()
add_executable("phasefield" phasefield.cc) set(MTL_TARGETS "")
target_link_dune_default_libraries("phasefield") list(APPEND MTL_TARGETS "phasefield" "phasefield2" "phasefield3" "phasefield4")
if (HAVE_ALBERTA)
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "phasefield") foreach(target ${MTL_TARGETS})
endif (HAVE_ALBERTA) add_executable(${target} ${target}.cc)
target_include_directories("phasefield" PRIVATE ${MTL_INCLUDE_DIRS}) target_link_dune_default_libraries(${target})
target_compile_definitions("phasefield" PRIVATE ${MTL_COMPILE_DEFINITIONS}) if (HAVE_ALBERTA)
target_compile_options("phasefield" PRIVATE -Wno-deprecated-declarations) add_dune_alberta_flags(${target})
endif (HAVE_ALBERTA)
target_include_directories(${target} PRIVATE ${MTL_INCLUDE_DIRS})
target_compile_definitions(${target} PRIVATE ${MTL_COMPILE_DEFINITIONS})
target_compile_options(${target} PRIVATE -Wno-deprecated-declarations)
endforeach()
endif () endif ()
\ No newline at end of file
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/common/timer.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
#include <dune/multimesh/mmhierarchiciterator.hh>
#include <dune/multimesh/multigridview.hh>
#include <dune/multimesh/utility/structuredgridbuilder.hh>
using namespace Dune;
template <class CoarseGridView, class FineGridView>
void printGrid(const CoarseGridView& coarse, const FineGridView& fine)
{
volatile std::size_t n = 0;
Dune::Timer t;
for (auto const& coarse_entity : elements(coarse)) {
for (auto const& fine_entity : childs(coarse_entity, fine)) {
n += fine.indexSet().index(fine_entity);
std::cout << "indices = [" << coarse.indexSet().index(coarse_entity) << ", "
<< fine.indexSet().index(fine_entity) << "]\n";
}
}
std::cout << n << "\n";
std::cout << "time: " << t.elapsed() << "\n";
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
#if HAVE_DUNE_ALUGRID
FieldVector<double,2> lower_left = {-1.5, -1.5};
FieldVector<double,2> bbox = {1.5, 1.5};
std::array<unsigned int,2> num_elements = {2, 2};
using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
using Factory = StructuredGridBuilder<Grid>;
auto gridPtr = Factory::createSimplexGrid(lower_left, bbox, num_elements);
auto& grid = *gridPtr;
grid.globalRefine(4);
printGrid(grid.levelGridView(2), grid.leafGridView());
#endif
}
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#if ! HAVE_DUNE_FUNCTIONS
#error "Require dune-functions!"
#endif
#if ! HAVE_DUNE_ALUGRID
#error "Require dune-alugrid!"
#endif
#if ! HAVE_ALBERTA
#error "Require Alberta!"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/alugrid/grid.hh>
#include <dune/multimesh/mmhierarchiciterator.hh>
#include <dune/multimesh/mmgridfactory.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include "phasefield2.hh"
using namespace Dune;
using namespace Dune::Indices;
struct MTLVector
{
using value_type = double;
MTLVector(mtl::dense_vector<double>& vec)
: vec_(vec) {}
void resize(std::size_t s) { vec_.change_dim(s); }
std::size_t size() const { return mtl::size(vec_); }
decltype(auto) operator[](std::size_t i) { return vec_[i]; }
decltype(auto) operator[](std::size_t i) const { return vec_[i]; }
decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) { return vec_[i[0]]; }
decltype(auto) operator[](ReservedVector<long unsigned int, 1> i) const { return vec_[i[0]]; }
mtl::dense_vector<double>& vec_;
};
template <class VectorType, class Basis>
void writeFile(VectorType& x, Basis const& basis, std::string filename)
{
auto gv = basis.gridView();
auto xWrapper = MTLVector(x);
auto xFct = Functions::makeDiscreteGlobalBasisFunction<double>(basis, TypeTree::hybridTreePath(), xWrapper);
VTKWriter<decltype(gv)> vtkWriter(gv);
vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write(filename);
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
assert(argc > 1);
std::string filename = argv[1];
using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
AlbertaReader<Grid> albertaReader;
GridFactory<Grid> factory;
albertaReader.readGrid(filename, factory);
std::unique_ptr<Grid> gridPtr(factory.createGrid());
auto& grid = *gridPtr;
double eps = 0.05;
auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; };
grid.globalRefine(10);
auto maxLevel = grid.maxLevel();
// refine grid[1] along phase-field interface
for (int i = 0; i < 4; ++i) {
std::cout << "refine " << i << "...\n";
for (auto const& elem : elements(grid.leafGridView())) {
auto geo = elem.geometry();
bool refine = false;
for (int j = 0; j < geo.corners(); ++j)
refine = refine || std::abs(signedDistFct(geo.corner(j))) < std::max(1,4-i)*eps;
if (refine)
grid.mark(1, elem);
}
grid.preAdapt();
grid.adapt();
grid.postAdapt();
}
using namespace Functions::BasisBuilder;
auto fine_basis = makeBasis(grid.leafGridView(), lagrange<1>());
using VectorType = mtl::dense_vector<double>;
using MatrixType = mtl::compressed2D<double>;
// interpolate phase-field to fine grid
VectorType phase(fine_basis.dimension());
auto phaseWrapper = MTLVector(phase);
interpolate(fine_basis, phaseWrapper, [eps,d=signedDistFct](auto const& x)
{
return 0.5*(1.0 - std::tanh(3.0*d(x)/eps));
});