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

interpolation extended to taylorhood example

parent d580cccc
......@@ -18,8 +18,11 @@ target_link_dune_default_libraries("hierarchiciterator")
add_executable("interpolate" interpolate.cc)
target_link_dune_default_libraries("interpolate")
add_executable("taylorhood" taylorhood.cc)
target_link_dune_default_libraries("taylorhood")
if (HAVE_ALBERTA)
add_dune_alberta_flags("interpolate")
add_dune_alberta_flags("taylorhood")
endif ()
find_package(MTL PATHS /opt/development/mtl4)
......
#pragma once
#include <cassert>
#include <tuple>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/functions/common/functionfromcallable.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
namespace Dune
{
template <class Grids>
class Interpolation
{
using HostGrid = typename Grids::HostGrid;
using HostEntity = typename HostGrid::Traits::template Codim<0>::Entity;
using Geometry = typename HostEntity::Geometry;
using LocalCoordinate = typename Geometry::LocalCoordinate;
using GlobalCoordinate = typename Geometry::GlobalCoordinate;
public:
Interpolation(Grids const& grids, std::size_t master = 1)
: grids_(grids)
, master_(master)
{}
template <class Vector, class Basis>
void operator()(Vector& oldVector, Basis const& newBasis) const
{
Vector newVector;
auto oldBasis = newBasis;
oldBasis.update(grids_[0].leafGridView());
(*this)(oldVector, oldBasis, newVector, newBasis);
std::swap(oldVector, newVector);
}
template <class Vector, class Basis>
void operator()(Vector const& oldVector, Basis const& oldBasis,
Vector& newVector, Basis const& newBasis) const
{
auto oldCoeff = Functions::istlVectorBackend(oldVector);
auto newCoeff = Functions::istlVectorBackend(newVector);
newCoeff.resize(newBasis);
for (const auto& e : master_leaf_elements(grids_, master_))
{
// NOTE: entities are not necessarily leaf entities
auto const& oldEntity = e[1 - master_];
auto const& newEntity = e[master_];
forEachLeafNode(newBasis.localView().tree(), [&](auto&&, auto const& treePath)
{
auto oldSubLocalView = Dune::Functions::subspaceBasis(oldBasis, treePath).localView();
auto newSubLocalView = Dune::Functions::subspaceBasis(newBasis, treePath).localView();
interpolate(oldCoeff, oldEntity, oldSubLocalView, newCoeff, newEntity, newSubLocalView);
});
}
}
template <class OldCoeff, class NewCoeff, class LocalView>
void interpolate(OldCoeff oldCoeff, const HostEntity& oldEntity, LocalView& oldLocalView,
NewCoeff newCoeff, const HostEntity& newEntity, LocalView& newLocalView) const
{
// newEntity is regular
if (newEntity.isLeaf()) {
newLocalView.bind(newEntity);
// temporary coefficient vector
auto const& node = newLocalView.tree();
std::vector<double> newDOFs(node.finiteElement().size());
// entity not changed
if (newEntity.level() == oldEntity.level() && newEntity.isLeaf() && oldEntity.isLeaf()) {
oldLocalView.bind(oldEntity);
interpolateToSame(oldCoeff, oldLocalView, newDOFs);
}
// entity was coarsened
else if (newEntity.level() == oldEntity.level() && newEntity.isLeaf() && !oldEntity.isLeaf()) {
interpolateToCoarse(oldCoeff, oldLocalView, oldEntity, newLocalView, newDOFs);
}
// entity was refined
else if (newEntity.level() > oldEntity.level() && newEntity.isLeaf() && oldEntity.isLeaf()) {
oldLocalView.bind(oldEntity);
interpolateToFine(oldCoeff, oldLocalView, newLocalView, newDOFs);
} else {
DUNE_THROW(Dune::Exception, "Unknown case in interpolation.");
}
// copy temporary DOFs to output vector
for (std::size_t i = 0; i < newDOFs.size(); ++i)
newCoeff[newLocalView.index(node.localIndex(i))] = newDOFs[i];
}
// newEntity is iregular
else if (oldEntity.isLeaf()) {
oldLocalView.bind(oldEntity);
int maxLevel = grids_.maxLevel(master_);
auto hItEnd = newEntity.hend(maxLevel);
std::vector<double> newChildDOFs;
for (auto hIt = newEntity.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
if (hIt->isLeaf()) {
auto const& child = *hIt;
newLocalView.bind(child);
auto const& node = newLocalView.tree();
newChildDOFs.resize(node.finiteElement().size());
interpolateToFine(oldCoeff, oldLocalView, newLocalView, newChildDOFs);
// copy childDOFs
for (std::size_t i = 0; i < newChildDOFs.size(); ++i)
newCoeff[newLocalView.index(node.localIndex(i))] = newChildDOFs[i];
}
}
}
// remaining iregular case
else {
int maxLevel = grids_.maxLevel(master_);
auto hItEnd = newEntity.hend(maxLevel);
std::vector<double> newChildDOFs;
for (auto hIt = newEntity.hbegin(maxLevel); hIt != hItEnd; ++hIt) {
if (hIt->isLeaf()) {
auto const& child = *hIt;
newLocalView.bind(child);
auto const& node = newLocalView.tree();
newChildDOFs.resize(node.finiteElement().size());
interpolateToCoarse(oldCoeff, oldLocalView, oldEntity, newLocalView, newChildDOFs);
// copy childDOFs
for (std::size_t i = 0; i < newChildDOFs.size(); ++i)
newCoeff[newLocalView.index(node.localIndex(i))] = newChildDOFs[i];
}
}
}
}
protected:
template <class Coeff, class OldLocalView>
void interpolateToSame(Coeff oldCoeff, OldLocalView const& oldLocalView,
std::vector<double>& newDOFs) const
{
auto const& node = oldLocalView.tree();
std::size_t oldFeSize = node.finiteElement().size();
assert(newDOFs.size() == oldFeSize);
for (std::size_t i = 0; i < newDOFs.size(); ++i)
newDOFs[i] = oldCoeff[oldLocalView.index(node.localIndex(i))];
}
// [[ expects: newLocalView.bound() ]]
template <class Coeff, class OldLocalView, class NewLocalView>
void interpolateToCoarse(Coeff oldCoeff, OldLocalView& oldLocalView, HostEntity const& oldEntity,
NewLocalView const& newLocalView, std::vector<double>& newDOFs) const
{
auto newGeo = newLocalView.element().geometry();
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);
}
// [[ expects: newLocalView.bound() && oldLocalView.bound() ]]
template <class Coeff, class OldLocalView, class NewLocalView>
void interpolateToFine(Coeff oldCoeff, OldLocalView const& oldLocalView,
NewLocalView const& newLocalView, std::vector<double>& newDOFs) const
{
auto oldGeo = oldLocalView.element().geometry();
auto newGeo = newLocalView.element().geometry();
auto evalOld = [&](LocalCoordinate const& x)
{
thread_local std::vector<FieldVector<double,1>> shapeValues;
auto const& node = oldLocalView.tree();
auto const& fe = node.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;
};
auto const& node = newLocalView.tree();
auto const& fe = node.finiteElement();
using FFC = Functions::FunctionFromCallable<double(LocalCoordinate), decltype(evalOld)>;
fe.localInterpolation().interpolate(FFC(evalOld), newDOFs);
}
std::pair<HostEntity,LocalCoordinate> findLeafEntity(HostEntity const& father, GlobalCoordinate const& global) const
{
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)
{
HostEntity 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));
}
private:
Grids const& grids_;
std::size_t master_;
};
} // end namespace Dune
// -*- 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";
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment