From 31f111d04024b09a6423119010f50b7041316588 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 21 Nov 2018 23:23:29 -0500 Subject: [PATCH] interpolation example added --- dune.module | 2 +- dune/multimesh/mmiterator.hh | 4 +- src/CMakeLists.txt | 6 + src/interpolate.cc | 281 +++++++++++++++++++++++++++++++++++ 4 files changed, 290 insertions(+), 3 deletions(-) create mode 100644 src/interpolate.cc diff --git a/dune.module b/dune.module index 07911dd..eea5b8d 100644 --- a/dune.module +++ b/dune.module @@ -8,4 +8,4 @@ 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 +Suggests: dune-uggrid dune-alugrid dune-foamgrid dune-functions dune-geometry dune-istl diff --git a/dune/multimesh/mmiterator.hh b/dune/multimesh/mmiterator.hh index cc689fc..4065c09 100644 --- a/dune/multimesh/mmiterator.hh +++ b/dune/multimesh/mmiterator.hh @@ -176,7 +176,7 @@ namespace Dune MultiMeshMasterLeafIterator (tag::begin_iterator, const GridImp* multiMesh, std::size_t master) : Super{tag::begin_iterator{}, multiMesh} , incrementAllowed_(multiMesh->size(), true) - , level_(multiMesh->size()) + , level_(multiMesh->size(), 0) , multiEntity_(multiMesh->size()) , master_(master) { @@ -225,7 +225,7 @@ namespace Dune // got down to leaf entity on master grid and until the master grid level for the other grids bool levelReached_impl (std::size_t i, const HostEntity& entity) const { - return (i == master_ && entity.isLeaf()) || entity.level() == level_[master_]; + return entity.isLeaf() || (i != master_ && entity.level() == level_[master_]); } private: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0b4a6fe..01eb4e6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,6 +16,12 @@ target_link_dune_default_libraries("multigridview") add_executable("hierarchiciterator" hierarchiciterator.cc) target_link_dune_default_libraries("hierarchiciterator") +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") diff --git a/src/interpolate.cc b/src/interpolate.cc new file mode 100644 index 0000000..7ad6e7e --- /dev/null +++ b/src/interpolate.cc @@ -0,0 +1,281 @@ +// -*- 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 + +using namespace Dune; +using namespace Dune::Indices; + +using VectorType = BlockVector>; + +template +void writeFile(VectorType const& x, Basis const& basis, std::string filename) +{ + auto gv = basis.gridView(); + auto xFct = Functions::makeDiscreteGlobalBasisFunction(basis, TypeTree::hybridTreePath(), x); + VTKWriter vtkWriter(gv); + vtkWriter.addVertexData(xFct, VTK::FieldInfo("u", VTK::FieldInfo::Type::scalar, 1)); + vtkWriter.write(filename); +} + +template +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(); + int newLevel = 0; + for (int j = 0; j < geo.corners(); ++j) + newLevel = std::max(newLevel, std::abs(signedDistFct(geo.corner(j))) < eps ? 12 : 0); + + int m = level < newLevel ? 1 : -1; + grids[0].mark(m, elems[0]); + grids[1].mark(m, elems[1]); + } +} + +template +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 +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 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> 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; + 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> 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; + 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) +{ + 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 WorldVector = typename HostGrid::template Codim<0>::Geometry::GlobalCoordinate; + using Grid = MultiMesh; + AlbertaReader albertaReader; + GridFactory factory(2); + albertaReader.readGrid(filename, factory); + + std::unique_ptr gridPtr(factory.createGrid()); + auto& grids = *gridPtr; + grids.globalRefine(6); + + 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 basis0 = makeBasis(grids.leafGridView(0), lagrange<2>()); // old grid + auto basis1 = makeBasis(grids.leafGridView(1), lagrange<2>()); // new grid + + // interpolate phase-field to fine grid + VectorType oldPhase(basis0.dimension()); + VectorType phase(basis1.dimension()); + Functions::interpolate(basis1, phase, [eps,d=signedDistFct](auto const& x) + { + return 0.5*(1.0 - std::tanh(2.0*d(x)/eps)); + }); + writeFile(phase, basis1, "phase_0"); + + Timer t; + double timeMarkElements = 0; + double timeAdapt = 0; + double timeBasisUpdate = 0; + double timeInterpolate = 0; + double timeWriteFile = 0; + + // assemble a sequence of systems for finer and finer grids + for (int i = 1; i < 20; ++i) { + std::cout << i << "\n"; + oldPhase.resize(basis0.dimension()); + oldPhase = phase; + + 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(); + basis1.update(grids.leafGridView(1)); + timeBasisUpdate += t.elapsed(); + + t.reset(); + mmInterpolate(grids, oldPhase, basis0, phase, basis1); // interpolate from old grid to new grid + timeInterpolate += t.elapsed(); + + t.reset(); + writeFile(phase, basis1, "phase_" + std::to_string(i)); + timeWriteFile += t.elapsed(); + + t.reset(); + grids[0].preAdapt(); + grids[0].adapt(); + grids[0].postAdapt(); + timeAdapt += t.elapsed(); + + t.reset(); + basis0.update(grids.leafGridView(0)); + timeBasisUpdate += 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"; +} -- GitLab