Commit 8901e45a authored by Praetorius, Simon's avatar Praetorius, Simon

Merge branch 'feature/entity_hierarchic_iterator'

parents db80a50b e8152f8a
......@@ -34,7 +34,7 @@ namespace Dune
using Super = std::vector<HostGridEntity>;
/// Constructor from std::vector
using Super::vector;
using std::vector<HostGridEntity>::vector;
/// Return a local geometry of source in 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 @@
#include <dune/grid/common/gridenums.hh>
#include "mmentity.hh"
#include "mmhierarchiciterator.hh"
namespace Dune
{
namespace tag
{
struct with_gridview {};
}
/** \brief Iterator over all entities of a given codimension and level of a grid.
* \ingroup MultiMesh
*/
......@@ -55,45 +51,8 @@ namespace Dune
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;
};
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;
};
using EntityTest = typename MultiMeshHierarchicIterator<HostGrid>::EntityTest;
using EntityStack = typename MultiMeshHierarchicIterator<HostGrid>::EntityStack;
public:
/// Constructor. Stores a pointer to the grid
......@@ -142,9 +101,10 @@ namespace Dune
: MultiMeshIterator{multiMesh, -1, endDummy}
{}
/// Construct an iterator from n gridViews of possibly different type
template <class... GridViews,
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)
, maxLevel_{gridViews.grid().maxLevel()...}
, multiEntity_(sizeof...(GridViews))
......@@ -162,11 +122,13 @@ namespace Dune
}
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>()...}
{}
#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>
MultiMeshIterator (const std::vector<std::variant<GV...>>& gridViews)
: incrementAllowed_(gridViews.size(), true)
......@@ -258,11 +220,14 @@ namespace Dune
}
// 3. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(i,child);
child = dereference(i)) {
auto 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);
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
......@@ -282,16 +247,19 @@ namespace Dune
auto& macroIt = macroIterators_[i];
auto const& macroEnd = macroEndIterators_[i];
assert( entityStack.empty() );
assert(entityStack.empty());
if (macroIt == macroEnd)
return;
// 2. go down in tree until leaf entity
for (auto child = dereference(i); !entityTest(i,child);
child = dereference(i)) {
// 1. go down in tree until leaf entity
auto 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);
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
......@@ -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:
std::vector<bool> incrementAllowed_;
std::vector<EntityTest> contains_;
......@@ -322,16 +285,14 @@ namespace Dune
std::vector<EntityStack> entityStacks_;
};
template <class> class MultiMesh;
template <class... GridViews>
inline auto multi_elements(GridViews const&... gridViews)
{
using GridView0 = std::tuple_element_t<0,std::tuple<GridViews...>>;
using Iterator = MultiMeshIterator<0,All_Partition,typename GridView0::Grid>;
using Range = IteratorRange<Iterator>;
return Range{ Iterator{tag::with_gridview{}, gridViews...},
Iterator{tag::with_gridview{}, true, gridViews...} };
return Range{ Iterator{gridViews...},
Iterator{true, gridViews...} };
}
} // end namespace Dune
......
......@@ -57,7 +57,7 @@ namespace Dune
/// Map entity to index.
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_);
}
......@@ -71,7 +71,7 @@ namespace Dune
/// Map a subentity to an index.
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
{
return std::visit([&e,i,codim](auto const* is) { return is->template subIndex<cc>(e,i,codim); }, indexSets_);
......
......@@ -60,7 +60,7 @@ namespace Dune
using LevelIterator = typename Partition<All_Partition>::LevelIterator;
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
......
......@@ -16,7 +16,7 @@ namespace Dune
template <class MultiMesh, class... Bases>
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
MultiBasis (MultiMesh const& mm, Bases&&... bases)
......@@ -27,7 +27,7 @@ namespace Dune
{}
/// Return the leafGridView of the MultiMesh
const MultiGridView& gridView () const { return multiGridView_; }
const MultiLeafGridView& gridView () const { return multiGridView_; }
/// Return the I'th basis
template <std::size_t I>
......@@ -63,7 +63,7 @@ namespace Dune
}
private:
MultiGridView multiGridView_;
MultiLeafGridView multiGridView_;
std::tuple<Bases...> bases_;
std::tuple<typename Bases::LocalView...> localViews_;
......
......@@ -13,6 +13,9 @@ target_link_dune_default_libraries("uggrid")
add_executable("multigridview" multigridview.cc)
target_link_dune_default_libraries("multigridview")
add_executable("hierarchiciterator" hierarchiciterator.cc)
target_link_dune_default_libraries("hierarchiciterator")
find_package(MTL PATHS /opt/development/mtl4)
if (MTL_FOUND)
set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
......@@ -25,12 +28,18 @@ if (MTL_FOUND)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
endif ()
add_executable("phasefield" phasefield.cc)
target_link_dune_default_libraries("phasefield")
if (HAVE_ALBERTA)
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "phasefield")
endif (HAVE_ALBERTA)
target_include_directories("phasefield" PRIVATE ${MTL_INCLUDE_DIRS})
target_compile_definitions("phasefield" PRIVATE ${MTL_COMPILE_DEFINITIONS})
target_compile_options("phasefield" PRIVATE -Wno-deprecated-declarations)
set(MTL_TARGETS "")
list(APPEND MTL_TARGETS "phasefield" "phasefield2" "phasefield3" "phasefield4")
foreach(target ${MTL_TARGETS})
add_executable(${target} ${target}.cc)
target_link_dune_default_libraries(${target})
if (HAVE_ALBERTA)
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 ()
\ 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));
});
writeFile(phase, fine_basis, "phase");
// assemble a sequence of systems for finer and finer grids
for (int l = 2; l < maxLevel; ++l) {
auto coarse_basis = makeBasis(grid.levelGridView(l), lagrange<1>());
VectorType rhs(coarse_basis.dimension());
MatrixType matrix(coarse_basis.dimension(), coarse_basis.dimension());
assembleSystem(coarse_basis, fine_basis, phase, matrix, rhs, eps);
VectorType u(coarse_basis.dimension());
u = 0.0;
umfpack_solve(matrix, u, rhs);
writeFile(u, coarse_basis, "u_" + std::to_string(l));
}
}
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#pragma once
#include <vector>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
#include <boost/numeric/mtl/mtl.hpp>
using namespace Dune;
template <class LocalView, class FineBasis>
void getLocalSystem(const LocalView& localView, const FineBasis& fine_basis,
mtl::dense_vector<double> const& phase,
mtl::dense2D<double>& elementMatrix,
mtl::dense_vector<double>& elementVector,
double eps)
{
using namespace Dune::Indices;
using MultiElement = MultiEntity<typename LocalView::GridView::Grid>;
const int dim = MultiElement::dimension;
elementMatrix.change_dim(localView.size(), localView.size());
elementVector.change_dim(localView.size());
set_to_zero(elementMatrix); // fills the entire matrix with zeroes
elementVector = 0;
// geometry of the basis element
auto geometry = localView.element().geometry();
const auto& localFE = localView.tree().finiteElement();
auto phaseLocalView = fine_basis.localView();
auto phaseIndexSet = fine_basis.localIndexSet();
for (auto const& e : childs(localView.element(), fine_basis.gridView()))
{
phaseLocalView.bind(e);
phaseIndexSet.bind(phaseLocalView);
const auto& phaseLocalFE = phaseLocalView.tree().finiteElement();
auto quadGeometry = e.geometry();
const auto& quad = QuadratureRules<double, dim>::rule(quadGeometry.type(), 2*localFE.localBasis().order());
for (const auto& quadPoint : quad)
{
auto local = MultiElement::local(e, localView.element(), quadPoint.position());
const auto jacobian = geometry.jacobianInverseTransposed(local);
const double dx = quadGeometry.integrationElement(quadPoint.position());
static std::vector<FieldMatrix<double,1,dim> > referenceGradients;
localFE.localBasis().evaluateJacobian(local, referenceGradients);
static std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
for (std::size_t i = 0; i < gradients.size(); ++i)
jacobian.mv(referenceGradients[i][0], gradients[i]);
static std::vector<FieldVector<double,1> > shapeValues;
localFE.localBasis().evaluateFunction(local, shapeValues);
static std::vector<FieldVector<double,1> > phaseShapeValues;
phaseLocalFE.localBasis().evaluateFunction(quadPoint.position(), phaseShapeValues);
double phase_atqp = 0.0;
for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) {
auto row = phaseIndexSet.index(phaseLocalView.tree().localIndex(i));
phase_atqp += phase[row[0]] * phaseShapeValues[i];
}