From e30e30f107c55b8371f4610d8a655953462a6825 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Thu, 28 Jun 2018 16:12:17 +0200 Subject: [PATCH] MultiMeshHierarchicIterator added to traverse leaf-childs of an element --- dune/multimesh/mmhierarchiciterator.hh | 178 +++++++++++++++++++++++++ dune/multimesh/mmiterator.hh | 85 ++++-------- dune/multimesh/utility/multibasis.hh | 6 +- src/CMakeLists.txt | 25 ++-- src/hierarchiciterator.cc | 58 ++++++++ src/phasefield2.cc | 147 ++++++++++++++++++++ src/phasefield2.hh | 122 +++++++++++++++++ src/phasefield3.cc | 146 ++++++++++++++++++++ src/phasefield3.hh | 125 +++++++++++++++++ src/phasefield4.cc | 146 ++++++++++++++++++++ src/phasefield4.hh | 127 ++++++++++++++++++ 11 files changed, 1092 insertions(+), 73 deletions(-) create mode 100644 dune/multimesh/mmhierarchiciterator.hh create mode 100644 src/hierarchiciterator.cc create mode 100644 src/phasefield2.cc create mode 100644 src/phasefield2.hh create mode 100644 src/phasefield3.cc create mode 100644 src/phasefield3.hh create mode 100644 src/phasefield4.cc create mode 100644 src/phasefield4.hh diff --git a/dune/multimesh/mmhierarchiciterator.hh b/dune/multimesh/mmhierarchiciterator.hh new file mode 100644 index 0000000..7f00e15 --- /dev/null +++ b/dune/multimesh/mmhierarchiciterator.hh @@ -0,0 +1,178 @@ +// -*- 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 +#include +#include + +#include +#include + +namespace Dune +{ + /** \brief Iterator over all entities of a given codimension and level of a grid. + * \ingroup MultiMesh + */ + template + class MultiMeshHierarchicIterator + : public ForwardIteratorFacade, + typename HostGrid::template Codim<0>::Entity, + typename HostGrid::template Codim<0>::Entity> + { + private: + template + friend class MultiMeshIterator; + + using HostEntity = typename HostGrid::template Codim<0>::Entity; + using EntityTest = std::function; + + struct EntityStackEntry + { + template + 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> + { + using Super = std::stack>; + + 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 + 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 + inline auto childs(const Entity& entity, const GridView& gridView) + { + using Iterator = MultiMeshHierarchicIterator; + return IteratorRange{ Iterator{entity, gridView}, Iterator{true} }; + } + +} // end namespace Dune + +#endif // DUNE_MULTIMESH_HIERARCHIC_ITERATOR_HH diff --git a/dune/multimesh/mmiterator.hh b/dune/multimesh/mmiterator.hh index 008a6a9..750f608 100644 --- a/dune/multimesh/mmiterator.hh +++ b/dune/multimesh/mmiterator.hh @@ -13,14 +13,10 @@ #include #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; - - struct EntityStackEntry - { - template - 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> - { - using Super = std::stack>; - - 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::EntityTest; + using EntityStack = typename MultiMeshHierarchicIterator::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 ...>::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 - 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 MultiMeshIterator (const std::vector>& 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(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 incrementAllowed_; std::vector contains_; @@ -322,16 +285,14 @@ namespace Dune std::vector entityStacks_; }; - template class MultiMesh; - template inline auto multi_elements(GridViews const&... gridViews) { using GridView0 = std::tuple_element_t<0,std::tuple>; using Iterator = MultiMeshIterator<0,All_Partition,typename GridView0::Grid>; using Range = IteratorRange; - return Range{ Iterator{tag::with_gridview{}, gridViews...}, - Iterator{tag::with_gridview{}, true, gridViews...} }; + return Range{ Iterator{gridViews...}, + Iterator{true, gridViews...} }; } } // end namespace Dune diff --git a/dune/multimesh/utility/multibasis.hh b/dune/multimesh/utility/multibasis.hh index 5d09085..55f249c 100644 --- a/dune/multimesh/utility/multibasis.hh +++ b/dune/multimesh/utility/multibasis.hh @@ -16,7 +16,7 @@ namespace Dune template 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 @@ -63,7 +63,7 @@ namespace Dune } private: - MultiGridView multiGridView_; + MultiLeafGridView multiGridView_; std::tuple bases_; std::tuple localViews_; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f7a999a..0b4a6fe 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/hierarchiciterator.cc b/src/hierarchiciterator.cc new file mode 100644 index 0000000..7a78bf4 --- /dev/null +++ b/src/hierarchiciterator.cc @@ -0,0 +1,58 @@ +// -*- 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 +#include // An initializer of MPI +#include // We use exceptions +#include + +#if HAVE_DUNE_ALUGRID +#include +#endif + +#include +#include +#include + +using namespace Dune; + +template +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 lower_left = {-1.5, -1.5}; + FieldVector bbox = {1.5, 1.5}; + std::array num_elements = {2, 2}; + using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>; + + using Factory = StructuredGridBuilder; + auto gridPtr = Factory::createSimplexGrid(lower_left, bbox, num_elements); + auto& grid = *gridPtr; + + grid.globalRefine(4); + + printGrid(grid.levelGridView(2), grid.leafGridView()); + +#endif +} diff --git a/src/phasefield2.cc b/src/phasefield2.cc new file mode 100644 index 0000000..1ae5901 --- /dev/null +++ b/src/phasefield2.cc @@ -0,0 +1,147 @@ +// -*- 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 +#include // An initializer of MPI +#include // We use exceptions + +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "phasefield2.hh" + +using namespace Dune; +using namespace Dune::Indices; + +struct MTLVector +{ + using value_type = double; + + MTLVector(mtl::dense_vector& 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 i) { return vec_[i[0]]; } + decltype(auto) operator[](ReservedVector i) const { return vec_[i[0]]; } + + mtl::dense_vector& vec_; +}; + +template +void writeFile(VectorType& x, Basis const& basis, std::string filename) +{ + auto gv = basis.gridView(); + auto xWrapper = MTLVector(x); + auto xFct = Functions::makeDiscreteGlobalBasisFunction(basis, TypeTree::hybridTreePath(), xWrapper); + VTKWriter 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 albertaReader; + GridFactory factory; + albertaReader.readGrid(filename, factory); + + std::unique_ptr 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; + using MatrixType = mtl::compressed2D; + + // 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)); + } +} diff --git a/src/phasefield2.hh b/src/phasefield2.hh new file mode 100644 index 0000000..c5cbb1f --- /dev/null +++ b/src/phasefield2.hh @@ -0,0 +1,122 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#pragma once + +#include + +#include +#include +#include +#include + +using namespace Dune; + +template +void getLocalSystem(const LocalView& localView, const FineBasis& fine_basis, + mtl::dense_vector const& phase, + mtl::dense2D& elementMatrix, + mtl::dense_vector& elementVector, + double eps) +{ + using namespace Dune::Indices; + + using MultiElement = MultiEntity; + 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::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 > referenceGradients; + localFE.localBasis().evaluateJacobian(local, referenceGradients); + + static std::vector > gradients(referenceGradients.size()); + for (std::size_t i = 0; i < gradients.size(); ++i) + jacobian.mv(referenceGradients[i][0], gradients[i]); + + static std::vector > shapeValues; + localFE.localBasis().evaluateFunction(local, shapeValues); + + static std::vector > 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]; + } + + for (size_t i=0; i +void assembleSystem(const Basis& basis, const FineBasis& fine_basis, Vector const& phase, Matrix& matrix, Vector& rhs, double eps) +{ + using namespace Dune::Indices; + + matrix.change_dim(basis.dimension(), basis.dimension()); + matrix = 0; + + rhs.change_dim(basis.dimension()); + rhs = 0; + + auto localView = basis.localView(); + auto localIndexSet = basis.localIndexSet(); + + mtl::mat::inserter > ins(matrix, 20); + for (const auto& element : elements(basis.gridView())) { + localView.bind(element); + localIndexSet.bind(localView); + + mtl::dense2D elementMatrix; + mtl::dense_vector elementVector; + getLocalSystem(localView, fine_basis, phase, elementMatrix, elementVector, eps); + + for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) { + auto row = localIndexSet.index(i); + if (elementVector[i] != 0) + rhs[row[0]] += elementVector[i]; + + for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) { + auto col = localIndexSet.index(j); + if (elementMatrix[i][j] != 0) + ins[row[0]][col[0]] += elementMatrix[i][j]; + } + } + } +} diff --git a/src/phasefield3.cc b/src/phasefield3.cc new file mode 100644 index 0000000..b681cbe --- /dev/null +++ b/src/phasefield3.cc @@ -0,0 +1,146 @@ +// -*- 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 +#include // An initializer of MPI +#include // We use exceptions + +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "phasefield3.hh" + +using namespace Dune; +using namespace Dune::Indices; + +struct MTLVector +{ + using value_type = double; + + MTLVector(mtl::dense_vector& 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 i) { return vec_[i[0]]; } + decltype(auto) operator[](ReservedVector i) const { return vec_[i[0]]; } + + mtl::dense_vector& vec_; +}; + +template +void writeFile(VectorType& x, Basis const& basis, std::string filename) +{ + auto gv = basis.gridView(); + auto xWrapper = MTLVector(x); + auto xFct = Functions::makeDiscreteGlobalBasisFunction(basis, TypeTree::hybridTreePath(), xWrapper); + VTKWriter 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 albertaReader; + GridFactory factory; + albertaReader.readGrid(filename, factory); + + std::unique_ptr gridPtr(factory.createGrid()); + auto& grid = *gridPtr; + + double eps = 0.05; + auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; }; + + grid.globalRefine(6); + + // refine grid[1] along phase-field interface + for (int i = 0; i < 8; ++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,6-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; + using MatrixType = mtl::compressed2D; + + // 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 < 7; ++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)); + } +} diff --git a/src/phasefield3.hh b/src/phasefield3.hh new file mode 100644 index 0000000..c94248d --- /dev/null +++ b/src/phasefield3.hh @@ -0,0 +1,125 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#pragma once + +#include + +#include +#include +#include +#include + +using namespace Dune; + +template +void getLocalSystem(const LocalView& localView, + const FineLocalView& fineLocalView, + const FineLocalIndexSet& fineLocalIndexSet, + const MultiElement& multiElement, + mtl::dense_vector const& phase, + mtl::dense2D& elementMatrix, + mtl::dense_vector& elementVector, + double eps) +{ + using namespace Dune::Indices; + + 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(); + + const auto& phaseLocalFE = fineLocalView.tree().finiteElement(); + + auto const& quadElement = multiElement.max(); + auto quadGeometry = quadElement.geometry(); + const auto& quad = QuadratureRules::rule(quadGeometry.type(), 2*localFE.localBasis().order()); + for (const auto& quadPoint : quad) + { + auto local = MultiElement::local(quadElement, localView.element(), quadPoint.position()); + auto phaseLocal = MultiElement::local(quadElement, fineLocalView.element(), quadPoint.position()); + + const auto jacobian = geometry.jacobianInverseTransposed(local); + const double dx = quadGeometry.integrationElement(quadPoint.position()); + + static std::vector > referenceGradients; + localFE.localBasis().evaluateJacobian(local, referenceGradients); + + static std::vector > gradients(referenceGradients.size()); + for (std::size_t i = 0; i < gradients.size(); ++i) + jacobian.mv(referenceGradients[i][0], gradients[i]); + + static std::vector > shapeValues; + localFE.localBasis().evaluateFunction(local, shapeValues); + + static std::vector > phaseShapeValues; + phaseLocalFE.localBasis().evaluateFunction(phaseLocal, phaseShapeValues); + + double phase_atqp = 0.0; + for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) { + auto row = fineLocalIndexSet.index(fineLocalView.tree().localIndex(i)); + phase_atqp += phase[row[0]] * phaseShapeValues[i]; + } + + for (size_t i=0; i +void assembleSystem(const Basis& coarse_basis, const FineBasis& fine_basis, + const Vector& phase, Matrix& matrix, Vector& rhs, double eps) +{ + using namespace Dune::Indices; + + matrix.change_dim(coarse_basis.dimension(), coarse_basis.dimension()); + matrix = 0; + + rhs.change_dim(coarse_basis.dimension()); + rhs = 0; + + auto localView = coarse_basis.localView(); + auto localIndexSet = coarse_basis.localIndexSet(); + + auto fineLocalView = fine_basis.localView(); + auto fineLocalIndexSet = fine_basis.localIndexSet(); + + mtl::mat::inserter > ins(matrix, 20); + for (const auto& multiElement : multi_elements(coarse_basis.gridView(), fine_basis.gridView())) { + localView.bind(multiElement[0]); + localIndexSet.bind(localView); + + fineLocalView.bind(multiElement[1]); + fineLocalIndexSet.bind(fineLocalView); + + mtl::dense2D elementMatrix; + mtl::dense_vector elementVector; + getLocalSystem(localView, fineLocalView, fineLocalIndexSet, multiElement, phase, elementMatrix, elementVector, eps); + + for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) { + auto row = localIndexSet.index(i); + if (elementVector[i] != 0) + rhs[row[0]] += elementVector[i]; + + for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) { + auto col = localIndexSet.index(j); + if (elementMatrix[i][j] != 0) + ins[row[0]][col[0]] += elementMatrix[i][j]; + } + } + } +} diff --git a/src/phasefield4.cc b/src/phasefield4.cc new file mode 100644 index 0000000..b681cbe --- /dev/null +++ b/src/phasefield4.cc @@ -0,0 +1,146 @@ +// -*- 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 +#include // An initializer of MPI +#include // We use exceptions + +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "phasefield3.hh" + +using namespace Dune; +using namespace Dune::Indices; + +struct MTLVector +{ + using value_type = double; + + MTLVector(mtl::dense_vector& 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 i) { return vec_[i[0]]; } + decltype(auto) operator[](ReservedVector i) const { return vec_[i[0]]; } + + mtl::dense_vector& vec_; +}; + +template +void writeFile(VectorType& x, Basis const& basis, std::string filename) +{ + auto gv = basis.gridView(); + auto xWrapper = MTLVector(x); + auto xFct = Functions::makeDiscreteGlobalBasisFunction(basis, TypeTree::hybridTreePath(), xWrapper); + VTKWriter 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 albertaReader; + GridFactory factory; + albertaReader.readGrid(filename, factory); + + std::unique_ptr gridPtr(factory.createGrid()); + auto& grid = *gridPtr; + + double eps = 0.05; + auto signedDistFct = [](auto const& x) { return x.two_norm() - 1.0; }; + + grid.globalRefine(6); + + // refine grid[1] along phase-field interface + for (int i = 0; i < 8; ++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,6-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; + using MatrixType = mtl::compressed2D; + + // 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 < 7; ++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)); + } +} diff --git a/src/phasefield4.hh b/src/phasefield4.hh new file mode 100644 index 0000000..fe97215 --- /dev/null +++ b/src/phasefield4.hh @@ -0,0 +1,127 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#pragma once + +#include + +#include +#include +#include +#include + +using namespace Dune; + +template +void getLocalSystem(const LocalView& localView, + const FineLocalView& fineLocalView, + const FineLocalIndexSet& fineLocalIndexSet, + const MultiElement& multiElement, + mtl::dense_vector const& phase, + mtl::dense2D& elementMatrix, + mtl::dense_vector& elementVector, + double eps) +{ + using namespace Dune::Indices; + + 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(); + + const auto& phaseLocalFE = fineLocalView.tree().finiteElement(); + + auto const& quadElement = multiElement.max(); + auto quadGeometry = quadElement.geometry(); + const auto& quad = QuadratureRules::rule(quadGeometry.type(), 2*localFE.localBasis().order()); + for (const auto& quadPoint : quad) + { + auto local = MultiElement::local(quadElement, localView.element(), quadPoint.position()); + auto phaseLocal = MultiElement::local(quadElement, fineLocalView.element(), quadPoint.position()); + + const auto jacobian = geometry.jacobianInverseTransposed(local); + const double dx = quadGeometry.integrationElement(quadPoint.position()); + + static std::vector > referenceGradients; + localFE.localBasis().evaluateJacobian(local, referenceGradients); + + static std::vector > gradients(referenceGradients.size()); + for (std::size_t i = 0; i < gradients.size(); ++i) + jacobian.mv(referenceGradients[i][0], gradients[i]); + + static std::vector > shapeValues; + localFE.localBasis().evaluateFunction(local, shapeValues); + + static std::vector > phaseShapeValues; + phaseLocalFE.localBasis().evaluateFunction(phaseLocal, phaseShapeValues); + + double phase_atqp = 0.0; + for (std::size_t i = 0; i < phaseLocalFE.size(); ++i) { + auto row = fineLocalIndexSet.index(fineLocalView.tree().localIndex(i)); + phase_atqp += phase[row[0]] * phaseShapeValues[i]; + } + + for (size_t i=0; i +void assembleSystem(const Basis& coarse_basis, const FineBasis& fine_basis, + const Vector& phase, Matrix& matrix, Vector& rhs, double eps) +{ + using namespace Dune::Indices; + + matrix.change_dim(coarse_basis.dimension(), coarse_basis.dimension()); + matrix = 0; + + rhs.change_dim(coarse_basis.dimension()); + rhs = 0; + + auto localView = coarse_basis.localView(); + auto localIndexSet = coarse_basis.localIndexSet(); + + auto fineLocalView = fine_basis.localView(); + auto fineLocalIndexSet = fine_basis.localIndexSet(); + + MultiGridView multiGridView(coarse_basis.gridView(), fine_basis.gridView()); + + mtl::mat::inserter > ins(matrix, 20); + for (const auto& multiElement : elements(multiGridView)) { + localView.bind(multiElement[0]); + localIndexSet.bind(localView); + + fineLocalView.bind(multiElement[1]); + fineLocalIndexSet.bind(fineLocalView); + + mtl::dense2D elementMatrix; + mtl::dense_vector elementVector; + getLocalSystem(localView, fineLocalView, fineLocalIndexSet, multiElement, phase, elementMatrix, elementVector, eps); + + for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) { + auto row = localIndexSet.index(i); + if (elementVector[i] != 0) + rhs[row[0]] += elementVector[i]; + + for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) { + auto col = localIndexSet.index(j); + if (elementMatrix[i][j] != 0) + ins[row[0]][col[0]] += elementMatrix[i][j]; + } + } + } +} -- 2.24.1