interpolate_noalias with PETSc-Backend
In Code below interpolate_noalias
doesn't work in case of using the PETSc-Backend. The issue is in line 55. Using the function in line 27 or 37 works.
The code runs without an error but in the Result you will see that the interpolation is not acting on the DOFs correctly.
With the following changes the code works:
- Compiling amdis with dune-istl
- using
interpolate
in line 55
With the following changes the code fails also:
- using first oder third order dune-curvedgrid
- refine the grid
-
tag::average{}
instead oftag::assign{}
#include <config.h>
#include <iostream>
#include <amdis/AMDiS.hpp>
#include <amdis/LocalOperators.hpp>
#include <amdis/ProblemStat.hpp>
#include <dune/vtk/pvdwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/curvedgrid/curvedgrid.hh>
#include <dune/foamgrid/foamgrid.hh>
using namespace AMDiS;
using namespace Dune::Functions::BasisFactory;
int main(int argc, char** argv)
{
Environment env(argc, argv);
using gridType = Dune::FoamGrid<2,3> ;
auto refGridPtr = MeshCreator<gridType>::create("mesh");
DOFVector surfaceDofVec{refGridPtr->leafGridView(), power<3>(lagrange(2))};
valueOf(surfaceDofVec).interpolate_noalias([](auto const& x){return x / x.two_norm(); } , tag::assign{});
Dune::CurvedGrid grid{*refGridPtr, valueOf(surfaceDofVec)};
ProblemStat curv("curv",grid, lagrange(1) );
curv.initialize(INIT_ALL);
AdaptInfo adaptInfo("adapt");
// create a DOFVector representing the surface normal
DOFVector nDofVec{curv.gridView(), power<3>(lagrange(2))};
valueOf(nDofVec).interpolate_noalias([](const auto& x) { return x / x.two_norm(); } , tag::assign{});
DOFVector surfaceUpdate{curv.gridView() , power<3>(lagrange(2))};
DOFVector surfaceOldDofVec{curv.gridView() , power<3>(lagrange(2))};
// visualization of the vector field
Dune::Vtk::LagrangeDataCollector dataCollector{curv.gridView(), 2};
using Writer = decltype(Dune::VtkUnstructuredGridWriter{dataCollector});
Dune::PvdWriter<Writer> writer{dataCollector};
writer.addPointData(valueOf(nDofVec), Dune::Vtk::FieldInfo{"n", .size=3, .rangeType=Dune::Vtk::RangeTypes::VECTOR});
for (adaptInfo.setTime(adaptInfo.startTime()) ; !adaptInfo.reachedEndTime() ; adaptInfo.setTime(adaptInfo.time()+adaptInfo.timestep()) , adaptInfo.incTimestepNumber()){
writer.writeTimestep(adaptInfo.time() , "/home/veitz/output/nematodynamic/meanCurvatureFlow/meanCurvatureFlow.pvd" );
surfaceOldDofVec.impl().vector() = surfaceDofVec.impl().vector() ;
// The code line below make issues in case of interpolate_noalias
valueOf(surfaceUpdate).interpolate_noalias(valueOf(surfaceOldDofVec) + adaptInfo.timestep() * X() , tag::assign{});
//valueOf(surfaceUpdate).interpolate(valueOf(surfaceOldDofVec) + adaptInfo.timestep() * X() , tag::assign{});
// ---------------------------------------------------------
surfaceDofVec.impl().vector() = surfaceUpdate.impl().vector() ;
}
return 0;
}