...
 
Commits (21)
---
before_script:
- source ~/toolchain
- export CMAKE_FLAGS="-DCMAKE_C_COMPILER='$CC' -DCMAKE_CXX_COMPILER='$CXX'"
dune-2.6 gcc-6-14:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-9-gcc-6-14
script: duneci-standard-test
dune-2.6 gcc-8-17:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-gcc-8-17
script: duneci-standard-test
dune-2.6 clang-6-17:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-clang-6-libcpp-17
script: duneci-standard-test
dune-git gcc-7-14:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-7-14
script: duneci-standard-test
dune-git gcc-8-17:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-8-noassert-17
script: duneci-standard-test
dune-git clang-6-17:
image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-clang-6-libcpp-17
script: duneci-standard-test
\ No newline at end of file
......@@ -7,5 +7,5 @@ Module: dune-multimesh
Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-common dune-grid
Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry dune-istl
Depends: dune-grid
Suggests: dune-alugrid dune-foamgrid dune-functions dune-istl dune-spgrid dune-uggrid
......@@ -19,20 +19,21 @@ namespace Dune
private:
using ctype = typename HostGrid::ctype;
/// Entity in the host grid
using HostGridEntity = typename HostGrid::Traits::template Codim<0>::Entity;
public:
enum { dimension = HostGrid::dimension };
enum { codimension = 0 };
enum { dimension = HostGridEntity::dimension };
enum { mydimension = HostGridEntity::mydimension };
/// The type of a local geometry
using LocalGeometry = MultiMeshLocalGeometry<dimension, dimension, HostGrid>;
using EntitySeed = MultiMeshEntitySeed<0, HostGrid>;
/// Entity in the host grid
using HostGridEntity = typename HostGrid::Traits::template Codim<0>::Entity;
/// Containertype
using Super = std::vector<HostGridEntity>;
/// Seed representing the vector of entities in the MultiMesh
using EntitySeed = MultiEntitySeed<0, HostGrid>;
public:
/// Constructor from std::vector
using std::vector<HostGridEntity>::vector;
......@@ -62,6 +63,39 @@ namespace Dune
return localGeometry(source, target).global(sourceLocal);
}
/// Return the partitionType of any entity in the multiEntity
// [[ expects: all partitionType are the same ]]
PartitionType partitionType () const
{
return max().partitionType();
}
/// Return the geometry of the entity with maximal level
typename HostGridEntity::Geometry geometry () const
{
return max().geometry();
}
/// Return the type of the entity with maximal level
GeometryType type () const
{
return max().type();
}
/// The entities are always regular, since irregular entities are not allowed in multimesh
bool isRegular() const { return true; }
/// Return the maximal level
int level () const
{
return max().level();
}
/// Return whether the entity with minimal level has boundary intersections
bool hasBoundaryIntersections () const
{
return min().hasBoundaryIntersections();
}
/// \brief Return the entity seed which contains sufficient information
/// to generate the entity again and uses as little memory as possible.
......@@ -69,9 +103,15 @@ namespace Dune
* The MultiMeshEntitySeed contains the HostGridEntitySeed and the index of
* the grid it is extracted from.
**/
EntitySeed seed (std::size_t entity_i) const
MultiMeshEntitySeed<0,HostGrid> seed (std::size_t idx) const
{
return {(*this)[idx], idx};
}
/// Return an entity seed to restore the multi-entity.
EntitySeed seed () const
{
return {(*this)[entity_i], entity_i};
return {*this};
}
/// Return the entity with maximal level
......
......@@ -10,16 +10,16 @@ namespace Dune
* \ingroup MultiMesh
*
*/
template <int codim, class GridImp>
template <int codim, class HostGrid>
class MultiMeshEntitySeed
{
protected:
/// Entity type of the hostgrid
typedef typename GridImp::HostGridType::Traits::template Codim<codim>::Entity HostEntity;
using HostEntity = typename HostGrid::Traits::template Codim<codim>::Entity;
/// EntitySeed type of the hostgrid
typedef typename GridImp::HostGridType::Traits::template Codim<codim>::EntitySeed HostEntitySeed;
using HostEntitySeed = typename HostGrid::Traits::template Codim<codim>::EntitySeed;
public:
enum { codimension = codim };
......@@ -61,6 +61,55 @@ namespace Dune
std::size_t gridIdx_;
};
template <class HG>
class MultiEntity;
// An entity seed for a vector of entities
template <int codim, class HostGrid>
class MultiEntitySeed
{
protected:
/// Entity type of the hostgrid
using HostEntity = typename HostGrid::Traits::template Codim<codim>::Entity;
/// EntitySeed type of the hostgrid
using HostEntitySeed = typename HostGrid::Traits::template Codim<codim>::EntitySeed;
public:
enum { codimension = codim };
/// Construct an empty (i.e. isValid() == false) seed.
MultiEntitySeed()
{}
/// \brief Create EntitySeed from hostgrid Entity
template <class MultiEntity>
MultiEntitySeed (const MultiEntity& multiEntity)
{
for (auto const& e : multiEntity)
hostEntitySeeds_.emplace_back(e.seed());
}
/// Get stored HostEntitySeed
const std::vector<HostEntitySeed>& hostEntitySeeds () const
{
return hostEntitySeeds_;
}
/// Check whether it is safe to create an Entity from this Seed
bool isValid () const
{
return hostEntitySeeds_.empty() ? false :
std::all_of(hostEntitySeeds_.begin(), hostEntitySeeds_.end(),
[](auto const& seed) { return seed.isValid(); });
}
private:
std::vector<HostEntitySeed> hostEntitySeeds_;
};
} // end namespace Dune
#endif // DUNE_MULTIMESH_ENTITYSEED_HH
......@@ -8,6 +8,8 @@
#include <vector>
#include <dune/common/hybridutilities.hh>
#include <dune/common/to_unique_ptr.hh>
#include <dune/common/version.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/grid/common/gridfactory.hh>
......@@ -147,6 +149,7 @@ namespace Dune
* Create n copies of the grid and store it in unique_pointers
* inside the multimesh grid.
**/
#if DUNE_VERSION_LT(DUNE_GRID,2,7)
virtual Grid* createGrid () override
{
Grid* multimesh = new Grid{};
......@@ -154,6 +157,15 @@ namespace Dune
multimesh->grids_.emplace_back(gridFactory->createGrid());
return multimesh;
}
#else
virtual ToUniquePtr<Grid> createGrid () override
{
auto multimesh = makeToUnique<Grid>();
for (auto& gridFactory : gridFactories_)
multimesh->grids_.emplace_back(gridFactory->createGrid());
return std::move(multimesh);
}
#endif
private:
std::vector<std::unique_ptr<GridFactory<HostGrid>>> gridFactories_;
......
......@@ -56,7 +56,7 @@ namespace Dune
using Grid = typename Traits::Grid;
/// Type of the corresponding GridType hosted by the MultiMesh
using HostGrid = typename GridImp::HostGridType;
using HostGrid = typename GridImp::HostGrid;
using IndexSet = typename Traits::IndexSet;
using Intersection = typename Traits::Intersection;
......@@ -84,7 +84,7 @@ namespace Dune
const HostGrid& grid (std::size_t i) const
{
return multiMesh_->grid(i);
return (*multiMesh_)[i];
}
/// Obtain the level-indexSet
......@@ -255,7 +255,7 @@ namespace Dune
const HostGrid& grid (std::size_t i) const
{
return multiMesh_[i];
return (*multiMesh_)[i];
}
/// Obtain the level-indexSet
......
......@@ -5,10 +5,11 @@
#include <numeric>
#include <stack>
#include <variant>
#include <dune/common/iteratorfacades.hh>
#include <dune/common/iteratorrange.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/grid/common/gridenums.hh>
namespace Dune
{
......@@ -30,7 +31,8 @@ namespace Dune
struct EntityStackEntry
{
template <class Entity>
template <class Entity,
std::enable_if_t<std::is_same<std::decay_t<Entity>, HostEntity>::value, int> = 0>
explicit EntityStackEntry (Entity&& entity)
: it(entity.hbegin(entity.level()+1))
, end(entity.hend(entity.level()+1))
......@@ -117,11 +119,17 @@ namespace Dune
// 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_ );
}
// 4. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack_.empty()) {
entityStack_.pop();
child = dereference();
}
assert(contains_(child) && "No valid child element found in gridView");
}
......@@ -150,11 +158,17 @@ namespace Dune
// 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_ );
}
// 2. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack_.empty()) {
entityStack_.pop();
child = dereference();
}
assert(contains_(child) && "No valid child element found in gridView");
}
......@@ -167,7 +181,7 @@ namespace Dune
};
template <class Entity, class GridView>
inline auto childs(const Entity& entity, const GridView& gridView)
inline auto childs (const Entity& entity, const GridView& gridView)
{
using Iterator = MultiMeshHierarchicIterator<typename GridView::Grid>;
return IteratorRange<Iterator>{ Iterator{entity, gridView}, Iterator{true} };
......
......@@ -6,6 +6,9 @@
#include <numeric>
#include <stack>
#include <type_traits>
#if DUNE_HAVE_CXX_VARIANT
#include <variant>
#endif
#include <dune/common/std/type_traits.hh>
#include <dune/grid/common/gridenums.hh>
......@@ -55,8 +58,8 @@ namespace Dune
// go to first leaf entity on all grids
entityStacks_.reserve(multiMesh->size());
for (std::size_t i = 0; i < multiMesh->size(); ++i) {
maxLevel_[i] = (*multiMesh)[i].maxLevel();
entityStacks_.emplace_back((*multiMesh)[i].maxLevel());
maxLevel_[i] = multiMesh->maxLevel(i);
entityStacks_.emplace_back(multiMesh->maxLevel(i));
}
}
......@@ -152,13 +155,18 @@ namespace Dune
// 3. go down in tree until leaf entity
auto child = dereference(i);
for (; !this->levelReached(i, child); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child);
assert(int(entityStack.size()) <= maxLevel_[i]);
}
// 4. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack.empty()) {
entityStack.pop();
child = dereference(i);
}
return entityStack.size();
// assert(contains_[i](child) && "No valid child element found in gridView");
}
using Super::increment;
......@@ -167,12 +175,12 @@ namespace Dune
{
std::size_t size = entityStacks_[i].size();
return std::accumulate(entityStacks_.begin(), entityStacks_.end(), true,
[i,size](bool allowed, auto const& entityStack) {
[size](bool allowed, auto const& entityStack) {
return allowed && (entityStack.size() <= size || entityStack.finished(size));
});
}
// got to first leaf entity on grid i
// go to first entity on grid i on the desired level
int initialIncrement (std::size_t i)
{
auto& entityStack = entityStacks_[i];
......@@ -186,16 +194,22 @@ namespace Dune
// 1. go down in tree until desired level is reached
auto child = dereference(i);
for (; !this->levelReached(i, child); child = dereference(i)) {
assert(child.isRegular() && "No irregular elements allowed in multi-mesh traversal");
entityStack.emplace(child);
assert(int(entityStack.size()) <= maxLevel_[i]);
}
// 2. go up in tree again to the first regular entity, since
// irregular element can not be traversed in a multi-mesh sense
while (!child.isRegular() && !entityStack.empty()) {
entityStack.pop();
child = dereference(i);
}
return entityStack.size();
// assert(contains_[i](child) && "No valid child element found in gridView");
}
using Super::initialIncrement;
// Return the current entity the hierarchic iterator or macro iterator pointing to
HostEntity dereference (std::size_t i) const
{
if (entityStacks_[i].empty()) {
......
......@@ -97,6 +97,8 @@ namespace Dune
using Super = GridDefaultImplementation<HG::dimension, HG::dimensionworld,
typename HG::ctype, MultiMeshFamily<HG> >;
using Self = MultiMesh;
public:
using HostGrid = HG;
......@@ -109,6 +111,8 @@ namespace Dune
/// The type used to store coordinates, inherited from the HostGrid
using ctype = typename HostGrid::ctype;
enum { dimension = HostGrid::dimension };
enum { dimensionworld = HostGrid::dimensionworld };
/// \brief Constructor, stores n copies of the hostgrid
/**
......@@ -127,7 +131,7 @@ namespace Dune
MultiMesh () = default;
/// Return the number of grids handled by this MultiMesh
std::size_t size() const
std::size_t size () const
{
return grids_.size();
}
......@@ -157,6 +161,11 @@ namespace Dune
[](int level, auto const& grid) { return std::max(level, grid->maxLevel()); });
}
/// Return maximum level in the `idx`th grid
int maxLevel (std::size_t idx) const
{
return grids_[idx]->maxLevel();
}
/// Iterator to first entity of given codim on level
template <int codim, PartitionIteratorType PiType = All_Partition>
......@@ -303,12 +312,23 @@ namespace Dune
/// Create Entity from EntitySeed
template <class EntitySeed>
typename Traits::template Codim<EntitySeed::codimension>::Entity
entity(const EntitySeed& seed) const
entity (const EntitySeed& seed) const
{
std::size_t gridIdx = this->getRealImplementation(seed).gridIndex();
return grids_[gridIdx]->entity(this->getRealImplementation(seed).hostEntitySeed());
}
/// Create a MultiEntity from a MultiEntitySeed
template <int codim>
MultiEntity<HostGrid> entity (const MultiEntitySeed<codim, Self>& seed) const
{
assert(seed.isValid());
MultiEntity<HostGrid> multiEntity;
auto const& hostEntitySeeds = seed.hostEntitySeeds();
for (std::size_t i = 0; i < hostEntitySeeds.size(); ++i)
multiEntity.emplace_back(grids_[i]->entity(hostEntitySeeds[i]));
return multiEntity;
}
/** @name Grid Refinement Methods */
/** @{ */
......@@ -320,6 +340,24 @@ namespace Dune
grid->globalRefine(refCount);
}
/// Mark all entities in a multiEntitity for refinement or coarsening
bool mark (int refCount, const MultiEntity<HostGrid>& entities)
{
bool b = false;
for (std::size_t i = 0; i < grids_.size(); ++i)
b = grids_[i]->mark(refCount, entities[i]) || b;
return b;
}
/// Return the marks of the the entities in a multiEntity
std::vector<int> getMark (const MultiEntity<HostGrid>& entities) const
{
std::vector<int> marks(grids_.size());
for (std::size_t i = 0; i < grids_.size(); ++i)
marks[i] = grids_[i]->getMark(entities[i]);
return marks;
}
/// Returns true, if at least one entity is marked for adaption in any of the
/// handled grids
bool preAdapt()
......@@ -344,29 +382,39 @@ namespace Dune
/** @} */
/** @name Grid Parallelization Methods */
/** @{ */
/// Size of the overlap on the leaf level
unsigned int overlapSize (int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[codim](unsigned int overlap, auto const& grid) {
return std::max(overlap, grid->overlapSize(codim)); });
}
/// Size of the ghost cell layer on the leaf level
unsigned int ghostSize (int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[codim](unsigned int ghost, auto const& grid) {
return std::max(ghost, grid->ghostSize(codim)); });
}
/// Size of the overlap on a given level
unsigned int overlapSize (int level, int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[level,codim](unsigned int overlap, auto const& grid) {
return std::max(overlap, grid->overlapSize(level,codim)); });
}
/// Size of the ghost cell layer on a given level
unsigned int ghostSize (int level, int codim) const
{
return 0u; // TODO: implement
return std::accumulate(grids_.begin(), grids_.end(), 0u,
[level,codim](unsigned int ghost, auto const& grid) {
return std::max(ghost, grid->ghostSize(level,codim)); });
}
/// dummy collective communication
......@@ -376,6 +424,22 @@ namespace Dune
return comm(0);
}
/// Rebalance all grids the same way based on an average load calculation.
bool loadBalance ()
{
assert(false && "Needs to be implemented.");
return false;
}
/// Rebalance all grids the same way based on an average load calculation.
template<class DataHandle>
bool loadBalance (DataHandle& data)
{
assert(false && "Needs to be implemented.");
return false;
}
/** @} */
public: // implementation details
......@@ -413,7 +477,7 @@ namespace Dune
template <class HostGrid, int codim>
struct hasEntityIterator<MultiMesh<HostGrid>, codim>
{
static const bool v = hasEntityIterator<HostGrid, codim>::v;
static const bool v = (codim == 0);
};
/** \brief has conforming level grids when host grid has
......
add_executable("testiterator" testiterator.cc)
target_link_dune_default_libraries("testiterator")
if (HAVE_ALBERTA)
add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 "testiterator")
endif (HAVE_ALBERTA)
add_executable("localrefinement" localrefinement.cc)
target_link_dune_default_libraries("localrefinement")
......@@ -10,42 +5,47 @@ target_link_dune_default_libraries("localrefinement")
add_executable("uggrid" uggrid.cc)
target_link_dune_default_libraries("uggrid")
add_executable("multigridview" multigridview.cc)
target_link_dune_default_libraries("multigridview")
add_executable("profiling" profiling.cc)
target_link_dune_default_libraries("profiling")
add_executable("benchmark" benchmark.cc)
target_link_dune_default_libraries("benchmark")
add_executable("hierarchiciterator" hierarchiciterator.cc)
target_link_dune_default_libraries("hierarchiciterator")
add_executable("testiterator" testiterator.cc)
target_link_dune_default_libraries("testiterator")
add_executable("interpolate" interpolate.cc)
target_link_dune_default_libraries("interpolate")
if (HAVE_ALBERTA)
add_dune_alberta_flags("interpolate")
endif ()
find_package(MTL PATHS /opt/development/mtl4)
if (MTL_FOUND)
set(CXX_ELEVEN_FEATURE_LIST "MOVE" "AUTO" "RANGEDFOR" "INITLIST" "STATICASSERT" "DEFAULTIMPL")
set(MTL_COMPILE_DEFINITIONS "")
foreach(feature ${CXX_ELEVEN_FEATURE_LIST})
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_WITH_${feature}")
endforeach()
if (HAVE_UMFPACK OR ENABLE_SUITESPARSE OR SuiteSparse_FOUND)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
endif ()
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()
add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 "testiterator")
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "benchmark")
endif (HAVE_ALBERTA)
if (DUNE_HAVE_CXX_VARIANT)
add_executable("multigridview" multigridview.cc)
target_link_dune_default_libraries("multigridview")
add_executable("hierarchiciterator" hierarchiciterator.cc)
target_link_dune_default_libraries("hierarchiciterator")
endif (DUNE_HAVE_CXX_VARIANT)
if (dune-functions_FOUND AND dune-alugrid_FOUND AND HAVE_ALBERTA)
add_executable("interpolate" interpolate.cc)
target_link_dune_default_libraries("interpolate")
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "interpolate")
add_executable("taylorhood" taylorhood.cc)
target_link_dune_default_libraries("taylorhood")
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 "taylorhood")
find_package(MTL QUIET PATHS /opt/development/mtl4)
if (MTL_FOUND)
foreach(target "phasefield" "phasefield2" "phasefield3" "phasefield4")
add_executable(${target} ${target}.cc)
target_link_dune_default_libraries(${target})
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${target})
target_include_directories(${target} PRIVATE ${MTL_INCLUDE_DIRS})
target_compile_definitions(${target} PRIVATE ${MTL_COMPILE_DEFINITIONS})
endforeach()
endif (MTL_FOUND)
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/filledarray.hh>
#include <dune/common/timer.hh>
#include <dune/grid/onedgrid.hh>
#include <dune/grid/yaspgrid.hh>
#if HAVE_ALBERTA
#include <dune/grid/albertagrid.hh>
#endif
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
#if HAVE_DUNE_FOAMGRID
#include <dune/foamgrid/foamgrid.hh>
#endif
#if HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif
#include <dune/multimesh/multimesh.hh>
#include <dune/multimesh/mmgridfactory.hh>
#include <dune/multimesh/utility/structuredgridbuilder.hh>
using namespace Dune;
template <class Grid>
std::size_t traverse(Grid const& grid)
{
std::size_t num = 0;
double vol = 0.0;
for (auto const& entities : elements(grid.leafGridView())) {
++num;
vol += entities.geometry().volume();
}
std::cout << "volume = " << vol << "\n";
return num;
}
template <class HostGrid>
void simplex_grid(std::string name, int refCount)
{
const int dow = HostGrid::dimensionworld;
using Factory = StructuredGridBuilder<MultiMesh<HostGrid> >;
GridFactory<MultiMesh<HostGrid> > gridFactory(2);
FieldVector<double,dow> lower; lower = -1.5;
FieldVector<double,dow> upper; upper = 1.5;
auto num_elements = filledArray<dow,unsigned int>(2);
Factory::createSimplexGrid(gridFactory, lower, upper, num_elements);
auto gridPtr = gridFactory.createGrid();
auto& grids = *gridPtr;
grids.globalRefine(refCount);
std::cout << name << "\n";
Timer t;
std::size_t num_multimesh = traverse(grids);
double time_multimesh = t.elapsed();
t.reset();
std::size_t num_singlemesh = traverse(grids[0]);
double time_singlemesh = t.elapsed();
std::cout << " [size(multimesh), size(singlemesh)]: [" << num_multimesh << ", " << num_singlemesh << "]\n";
std::cout << " multimesh: " << time_multimesh << "sec\n";
std::cout << " singlemesh: " << time_singlemesh << "sec\n";
std::cout << " ratio: " << (time_multimesh / time_singlemesh) << "\n";
}
template <class HostGrid>
void cubes_grid(std::string name, int refCount)
{
const int dow = HostGrid::dimensionworld;
using Factory = StructuredGridBuilder<MultiMesh<HostGrid> >;
GridFactory<MultiMesh<HostGrid> > gridFactory(2);
FieldVector<double,dow> lower; lower = -1.5;
FieldVector<double,dow> upper; upper = 1.5;
auto num_elements = filledArray<dow,unsigned int>(2);
Factory::createCubeGrid(gridFactory, lower, upper, num_elements);
auto gridPtr = gridFactory.createGrid();
auto& grids = *gridPtr;
grids.globalRefine(refCount);
std::cout << name << "\n";
Timer t;
std::size_t num_multimesh = traverse(grids);
double time_multimesh = t.elapsed();
t.reset();
std::size_t num_singlemesh = traverse(grids[0]);
double time_singlemesh = t.elapsed();
std::cout << " [size(multimesh), size(singlemesh)]: [" << num_multimesh << ", " << num_singlemesh << "]\n";
std::cout << " multimesh: " << time_multimesh << "sec\n";
std::cout << " singlemesh: " << time_singlemesh << "sec\n";
std::cout << " ratio: " << (time_multimesh / time_singlemesh) << "\n";
}
template <class HostGrid, class... Args>
void structured_grid(std::string name, int refCount, Args&&... args)
{
MultiMesh<HostGrid> grids(2, std::forward<Args>(args)...);
grids.globalRefine(refCount);
std::cout << name << "\n";
Timer t;
std::size_t num_multimesh = traverse(grids);
double time_multimesh = t.elapsed();
t.reset();
std::size_t num_singlemesh = traverse(grids[0]);
double time_singlemesh = t.elapsed();
std::cout << " [size(multimesh), size(singlemesh)]: [" << num_multimesh << ", " << num_singlemesh << "]\n";
std::cout << " multimesh: " << time_multimesh << "sec\n";
std::cout << " singlemesh: " << time_singlemesh << "sec\n";
std::cout << " ratio: " << (time_multimesh / time_singlemesh) << "\n";
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
FieldVector<double,1> lower1{-1.5};
FieldVector<double,2> lower2{-1.5,-1.5};
FieldVector<double,3> lower3{-1.5,-1.5,-1.5};
FieldVector<double,1> upper1{1.5};
FieldVector<double,2> upper2{1.5,1.5};
FieldVector<double,3> upper3{1.5,1.5,1.5};
structured_grid<OneDGrid>("OneDGrid", 16, 2, lower1[0], upper1[0]);
structured_grid<YaspGrid<1,EquidistantOffsetCoordinates<double,1>>>("YaspGrid<1>", 16, lower1, upper1, std::array<int,1>{2});
structured_grid<YaspGrid<2,EquidistantOffsetCoordinates<double,2>>>("YaspGrid<2>", 8, lower2, upper2, std::array<int,2>{2,2});
structured_grid<YaspGrid<3,EquidistantOffsetCoordinates<double,3>>>("YaspGrid<3>", 6, lower3, upper3, std::array<int,3>{2,2,2});
#if HAVE_ALBERTA
// simplex_grid<AlbertaGrid<1,1>>("AlbertaGrid<1>", 8);
// simplex_grid<AlbertaGrid<2,2>>("AlbertaGrid<2>", 6);
simplex_grid<AlbertaGrid<ALBERTA_DIM,ALBERTA_DIM>>("AlbertaGrid<" + std::to_string(ALBERTA_DIM) + ">", 4);
#endif
#if HAVE_UG
// simplex_grid<UGGrid<1>>("UGGrid<1>", 8);
simplex_grid<UGGrid<2>>("UGGrid<2>", 8);
simplex_grid<UGGrid<3>>("UGGrid<3>", 4);
cubes_grid<UGGrid<2>>("UGGrid<2,cube>", 8);
cubes_grid<UGGrid<3>>("UGGrid<3,cube>", 4);
#endif
#if HAVE_DUNE_ALUGRID
simplex_grid<Dune::ALUGrid<2, 2, simplex, conforming> >("ALUGrid<2,conforming>", 8);
// simplex_grid<Dune::ALUGrid<3, 3, simplex, conforming> >("ALUGrid<3,conforming>", 4);
simplex_grid<Dune::ALUGrid<2, 2, simplex, nonconforming> >("ALUGrid<2,nonconforming>", 6);
simplex_grid<Dune::ALUGrid<3, 3, simplex, nonconforming> >("ALUGrid<3,nonconforming>", 4);
cubes_grid<Dune::ALUGrid<2, 2, cube, nonconforming> >("ALUGrid<2,cube,nonconforming>", 6);
cubes_grid<Dune::ALUGrid<3, 3, cube, nonconforming> >("ALUGrid<3,cube,nonconforming>", 4);
#endif
#if HAVE_DUNE_SPGRID
structured_grid<SPGrid<double,1>>("SPGrid<1>", 16, lower1, upper1, std::array<int,1>{2});
structured_grid<SPGrid<double,2>>("SPGrid<2>", 8, lower2, upper2, std::array<int,2>{2,2});
structured_grid<SPGrid<double,3>>("SPGrid<3>", 6, lower3, upper3, std::array<int,3>{2,2,2});
#endif
}
......@@ -4,6 +4,8 @@
# include "config.h"
#endif
#if DUNE_HAVE_CXX_VARIANT
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
......@@ -53,6 +55,9 @@ int main(int argc, char** argv)
grid.globalRefine(4);
printGrid(grid.levelGridView(2), grid.leafGridView());
#endif
}
#else // DUNE_HAVE_CXX_VARIANT
int main() {}
#endif // DUNE_HAVE_CXX_VARIANT
\ No newline at end of file
......@@ -29,6 +29,7 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/albertagrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
......@@ -40,6 +41,8 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include "interpolation.hh"
using namespace Dune;
using namespace Dune::Indices;
......@@ -59,126 +62,17 @@ template <class Grids, class SignedDistFct>
void markElements(Grids& grids, SignedDistFct&& signedDistFct, double eps = 0.1)
{
for (auto const& elems : elements(grids.leafGridView())) {
auto geo = elems[0].geometry();
int level = elems[0].level();
auto geo = elems.geometry();
int level = elems.level();
int newLevel = 0;
for (int j = 0; j < geo.corners(); ++j)
newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 12 : 0);
newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 8 : 0);
int m = level < newLevel ? 1 : -1;
grids[0].mark(m, elems[0]);
grids[1].mark(m, elems[1]);
grids.mark(m, elems);
}
}
template <class Element, class GlobalCoordinate>
auto findLeafEntity(Element const& father, GlobalCoordinate const& global)
{
const int childLevel = father.level()+1;
// loop over all child Entities
const auto end = father.hend(childLevel);
for (auto it = father.hbegin(childLevel); it != end; ++it)
{
Element child = *it;
auto geo = child.geometry();
auto local = geo.local(global);
if (referenceElement(geo).checkInside(local))
{
// return if we found the leaf, else search through the child entites
if (child.isLeaf())
return std::make_pair(std::move(child), local);
else
return findLeafEntity(child, global);
}
}
assert(father.isLeaf());
return std::make_pair(father, father.geometry().local(global));
}
template <class Grids, class OldBasis, class NewBasis>
void mmInterpolate(Grids const& grids,
VectorType const& oldCoeff, OldBasis const& oldBasis,
VectorType& newCoeff, NewBasis const& newBasis)
{
using namespace Dune::Indices;
using LocalCoordinate = typename Grids::Traits::template Codim<0>::Geometry::LocalCoordinate;
auto oldLocalView = oldBasis.localView();
auto newLocalView = newBasis.localView();
newCoeff.resize(newBasis.dimension());
for (const auto& e : master_leaf_elements(grids, 1))
{
auto const& oldEntity = e[0];
auto const& newEntity = e[1]; // always a leaf entity
assert(newEntity.isLeaf());
newLocalView.bind(newEntity);
auto newGeo = newEntity.geometry();
auto const& node = newLocalView.tree();
auto const& fe = node.finiteElement();
std::size_t newFeSize = fe.size();
std::vector<double> newDOFs(newFeSize);
// entity not changed
if (newEntity.level() == oldEntity.level() && oldEntity.isLeaf()) {
oldLocalView.bind(oldEntity);
std::size_t oldFeSize = oldLocalView.tree().finiteElement().size();
assert(newFeSize == oldFeSize);
for (std::size_t i = 0; i < newFeSize; ++i)
newDOFs[i] = oldCoeff[oldLocalView.index(oldLocalView.tree().localIndex(i))];
}
// entity was coarsened
if (newEntity.level() == oldEntity.level() && !oldEntity.isLeaf()) {
auto evalLeaf = [&](LocalCoordinate const& x) {
thread_local std::vector<FieldVector<double,1>> shapeValues;
auto leafLocal = findLeafEntity(oldEntity, newGeo.global(x));
oldLocalView.bind(leafLocal.first);
auto const& node = oldLocalView.tree();
auto const& fe = node.finiteElement();
fe.localBasis().evaluateFunction(leafLocal.second, shapeValues);
double y = 0;
for (std::size_t i = 0; i < shapeValues.size(); ++i)
y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
return y;
};
auto const& node = newLocalView.tree();
auto const& fe = node.finiteElement();
using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalLeaf)>;
fe.localInterpolation().interpolate(FFC(evalLeaf), newDOFs);
}
// entity was refined
else if (newEntity.level() >= oldEntity.level()) {
oldLocalView.bind(oldEntity);
auto oldGeo = oldEntity.geometry();
auto evalOld = [&](LocalCoordinate const& x)
{
thread_local std::vector<FieldVector<double,1>> shapeValues;
auto const& fe = oldLocalView.tree().finiteElement();
fe.localBasis().evaluateFunction(oldGeo.local(newGeo.global(x)), shapeValues);
double y = 0;
for (std::size_t i = 0; i < shapeValues.size(); ++i)
y += shapeValues[i] * oldCoeff[oldLocalView.index(node.localIndex(i))];
return y;
};
using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalOld)>;
fe.localInterpolation().interpolate(FFC(evalOld), newDOFs);
}
for (std::size_t i = 0; i < newFeSize; ++i)
newCoeff[newLocalView.index(node.localIndex(i))] = newDOFs[i];
}
}
int main(int argc, char** argv)
{
......@@ -188,7 +82,8 @@ int main(int argc, char** argv)
std::string filename = argv[1];
// using HostGrid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
using HostGrid = Dune::AlbertaGrid<2,2>;
// using HostGrid = Dune::AlbertaGrid<2,2>;
using HostGrid = Dune::UGGrid<2>;
using WorldVector = typename HostGrid::template Codim<0>::Geometry::GlobalCoordinate;
using Grid = MultiMesh<HostGrid>;
AlbertaReader<Grid> albertaReader;
......@@ -197,7 +92,7 @@ int main(int argc, char** argv)
std::unique_ptr<Grid> gridPtr(factory.createGrid());
auto& grids = *gridPtr;
grids.globalRefine(6);
grids.globalRefine(3);
double eps = 0.1;
WorldVector shift(0.0);
......@@ -231,6 +126,8 @@ int main(int argc, char** argv)
double timeInterpolate = 0;
double timeWriteFile = 0;
Interpolation<Grid> interpolation(grids);
// assemble a sequence of systems for finer and finer grids
for (int i = 1; i < 20; ++i) {
std::cout << i << "\n";
......@@ -254,7 +151,7 @@ int main(int argc, char** argv)
timeBasisUpdate += t.elapsed();
t.reset();
mmInterpolate(grids, oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid
interpolation(oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid
timeInterpolate += t.elapsed();
t.reset();
......
This diff is collapsed.
......@@ -4,6 +4,8 @@
# include "config.h"
#endif
#if DUNE_HAVE_CXX_VARIANT
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
......@@ -73,3 +75,7 @@ int main(int argc, char** argv)
#endif
}
#else // DUNE_HAVE_CXX_VARIANT
int main() {}
#endif // DUNE_HAVE_CXX_VARIANT
......@@ -39,7 +39,7 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include "phasefield3.hh"
#include "phasefield4.hh"
using namespace Dune;
using namespace Dune::Indices;
......
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/timer.hh>
#include <dune/grid/onedgrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/multimesh/multimesh.hh>
using namespace Dune;
template <class Grid>
void traverse(Grid const& grid)
{
double area = 0;
int num_boundary_elements = 0;
for (auto const& entities : elements(grid.leafGridView())) {
area += entities.geometry().volume();
num_boundary_elements += entities.hasBoundaryIntersections() ? 1 : 0;
}
std::cout << " boundaryIntersections: " << num_boundary_elements << "\n";
std::cout << " area: " << area << "\n";
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
using HostGrid = YaspGrid<2,EquidistantOffsetCoordinates<double,2>>;
FieldVector<double,2> lower2{-1.5,-1.5};
FieldVector<double,2> upper2{1.5,1.5};
MultiMesh<HostGrid> grids(2, lower2, upper2, std::array<int,2>{32,32});
std::cout << "YaspGrid<2>\n";
Timer t;
traverse(grids);
std::cout << " time:" << t.elapsed() << "sec\n";
}
// -*- 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/common/timer.hh>
#include <dune/alugrid/grid.hh>
#include <dune/multimesh/mmhierarchiciterator.hh>
#include <dune/multimesh/mmgridfactory.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/albertagrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/istl/bvector.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/taylorhoodbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include "interpolation.hh"
using namespace Dune;
using namespace Dune::Indices;
template <class DiscreteFunction>
void writeFile(DiscreteFunction const& x, std::string filename, VTK::FieldInfo info)
{
auto gv = x.basis().gridView();
VTKWriter<decltype(gv)> vtkWriter(gv);
vtkWriter.addVertexData(x, info);
vtkWriter.write(filename);
}
template <class Grids, class SignedDistFct>
void markElements(Grids& grids, SignedDistFct&& signedDistFct, double eps = 0.1)
{
for (auto const& elems : elements(grids.leafGridView())) {
auto geo = elems.geometry();
int level = elems.level();
int newLevel = 0;
for (int j = 0; j < geo.corners(); ++j)
newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 6 : 3);
int m = level < newLevel ? 1 : -1;
grids.mark(m, elems);
}
}
int main(int argc, char** argv)
{
MPIHelper::instance(argc, argv);
assert(argc > 1);
std::string filename = argv[1];
// using HostGrid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
using HostGrid = Dune::AlbertaGrid<2,2>;
//using HostGrid = Dune::UGGrid<2>;
using WorldVector = typename HostGrid::template Codim<0>::Geometry::GlobalCoordinate;
using Grid = MultiMesh<HostGrid>;
AlbertaReader<Grid> albertaReader;
GridFactory<Grid> factory(2);
albertaReader.readGrid(filename, factory);
std::unique_ptr<Grid> gridPtr(factory.createGrid());
auto& grids = *gridPtr;
grids.globalRefine(3);
double eps = 0.1;
WorldVector shift(0.0);
auto signedDistFct = [&shift](WorldVector const& x) { return (x - shift).two_norm() - 0.7; };
// initial refinement
for (int i = 0; i < 8; ++i) {
markElements(grids, signedDistFct, eps);
grids.preAdapt();
grids.adapt();
grids.postAdapt();
}
using namespace Functions::BasisBuilder;
auto basis = makeBasis(grids.leafGridView(1), taylorHood()); // new grid
using VectorType = BlockVector<BlockVector<FieldVector<double,1>>>;
// interpolate phase-field to fine grid
VectorType coefficients;
Functions::interpolate(Functions::subspaceBasis(basis, _0), coefficients, [](auto const& x)
{
return FieldVector<double,2>{std::sin(x[0]), std::cos(x[1])};
});
Functions::interpolate(Functions::subspaceBasis(basis, _1), coefficients, [](auto const& x)
{
return std::sin(x[0])*std::cos(x[1]);
});
using VelocityRange = FieldVector<double,2>;
using PressureRange = double;
VTK::FieldInfo velocityInfo("u", VTK::FieldInfo::Type::vector, 2);
auto velocity
= Functions::makeDiscreteGlobalBasisFunction<VelocityRange>(
Functions::subspaceBasis(basis, _0), coefficients);
VTK::FieldInfo pressureInfo("p", VTK::FieldInfo::Type::scalar, 1);
auto pressure
= Functions::makeDiscreteGlobalBasisFunction<PressureRange>(
Functions::subspaceBasis(basis, _1), coefficients);
writeFile(velocity, "velocity_0", velocityInfo);
writeFile(pressure, "pressure_0", pressureInfo);
Timer t;
double timeMarkElements = 0;
double timeAdapt = 0;
double timeBasisUpdate = 0;
double timeInterpolate = 0;
double timeWriteFile = 0;
Interpolation<Grid> interpolation(grids);
// assemble a sequence of systems for finer and finer grids
for (int i = 1; i < 20; ++i) {
std::cout << i << "\n";
shift[0] += 0.01;
t.reset();
markElements(grids, signedDistFct, eps);
timeMarkElements += t.elapsed();
t.reset();
grids[1].preAdapt();
grids[1].adapt();
grids[1].postAdapt();
timeAdapt += t.elapsed();
t.reset();
basis.update(grids.leafGridView(1));
timeBasisUpdate += t.elapsed();
t.reset();
interpolation(coefficients, basis); // interpolate from old grid to new grid
timeInterpolate += t.elapsed();
t.reset();
writeFile(velocity, "velocity_" + std::to_string(i), velocityInfo);
writeFile(pressure, "pressure_" + std::to_string(i), pressureInfo);
timeWriteFile += t.elapsed();
t.reset();
grids[0].preAdapt();
grids[0].adapt();
grids[0].postAdapt();
timeAdapt += t.elapsed();
}
std::cout << "timings:\n";
std::cout << " time(markElements) = " << timeMarkElements << "\n";
std::cout << " time(adapt) = " << timeAdapt << "\n";
std::cout << " time(basisUpdate) = " << timeBasisUpdate << "\n";
std::cout << " time(interpolate) = " << timeInterpolate << "\n";
std::cout << " time(writeFile) = " << timeWriteFile << "\n";
}
dune_add_test(SOURCES testvolume.cc)
dune_add_test(SOURCES testnumentities.cc)
dune_add_test(SOURCES testtransform.cc)
dune_add_test(SOURCES testgridviews.cc)
\ No newline at end of file
if (DUNE_HAVE_CXX_VARIANT)
dune_add_test(SOURCES testgridviews.cc)
endif (DUNE_HAVE_CXX_VARIANT)
\ No newline at end of file