#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <sys/stat.h>
#include <cmath>
#include <math.h>
#include <vector>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/densematrix.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/typetree/treepath.hh>
template <typename Number, typename LType>
{
Number time = 0.;
Number finaltime = 1.;
Number dt = 0.2;
Number othertime = 0.;
Number grav = -9.81;
Number rho = 1.;
Number viscosity = 1.002e-3;
Number pw_ini = -1.;
Number K_intristic = 1.76e-4;
Number phi = 0.32;
Number alpha = 3.5;
Number n = 3.19;
Number m;
Number C0_ini = 0.;
Number C1_ini = 1.;
const LType& L;
Number inflow_ = 0.001;
Number inflow_C0 = 1.;
Number inflow_C1 = 0.;
public:
: L(L_)
{
m = 1.-1./n;
}
{
time = t_;
}
{
return time;
}
{
return finaltime;
}
{
return dt;
}
{
return L[0];
}
{
return L[1];
}
{
dt = dt_;
}
{
dt *= coef;
}
{
othertime = t_;
}
{
return othertime;
}
Number
Sw (
const Number& pw_)
const
{
return (pw_ < 0.) ? 1./pow(1.+pow(-alpha*pw_,n),m) : 1.;
}
Number
K (
const Number& Sw_)
const
{
}
{
Number pom = 1. - pow( 1.-pow(eSw,1./m) ,m);
Number result=sqrt(eSw)*pom*pom;
return result;
}
{
return phi;
}
{
return rho;
}
{
return grav;
}
Number
getD0 (
const Number& Sw_)
const
{
return D0*phi*Sw_*Sw_*Sw_*cbrt(phi*Sw_);
}
Number
getD1 (
const Number& Sw_)
const
{
return D1*phi*Sw_*Sw_*Sw_*cbrt(phi*Sw_);
}
template <typename E, typename X>
Number
g (
const E&
e,
const X& x)
const
{
auto global =
e.geometry().global(x);
return pw_ini+global[1]*grav*rho;
}
template <typename E, typename X>
Number
g_C0 (
const E&
e,
const X& x)
const
{
return C0_ini;
}
template <typename E, typename X>
Number
g_C1 (
const E&
e,
const X& x)
const
{
return C1_ini;
}
Number
harg(Number a, Number
b)
{
return 0.;
}
template <typename I, typename X>
bool b (
const I& i,
const X& x)
{
return false;
}
template <typename I, typename X>
bool bC0 (
const I& i,
const X& x)
{
return false;
}
template <typename I, typename X>
bool bC1 (
const I& i,
const X& x)
{
return false;
}
template <typename I, typename X>
Number
inflow (
const I& i,
const X& x)
{
return inflow_;
}
template <typename I, typename X>
Number
inflowC0 (
const I& i,
const X& x)
{
return inflow_C0;
}
template <typename I, typename X>
Number
inflowC1 (
const I& i,
const X& x)
{
return inflow_C1;
}
};
{};
template<typename Param>
{
Param& param;
using RF = typename Param::value_type;
public:
: param(param_)
{}
{
param.setTime(t);
}
template<typename IG, typename LFSU, typename X,
typename LFSV, typename R>
const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
R& r_i, R& r_o) const
{
auto cell_inside =
ig.inside();
auto cell_outside =
ig.outside();
auto insidegeo = cell_inside.geometry();
auto outsidegeo = cell_outside.geometry();
auto inside_global = insidegeo.center();
auto outside_global = outsidegeo.center();
inside_global -= outside_global;
auto distance = inside_global.two_norm();
auto facegeo =
ig.geometry();
auto volume = facegeo.volume();
auto pw_i = x_i(lfsu_i,0);
auto Sw_i = param.Sw(pw_i);
auto pw_o = x_o(lfsu_o,0);
auto Sw_o = param.Sw(pw_o);
auto dpwdn = (pw_o-pw_i)/distance;
auto g = param.getgrav()*param.getrho()*inside_global[1]/distance;
auto K_i = param.K(Sw_i);
auto K_o = param.K(Sw_o);
auto K = param.harg(K_i,K_o);
auto q = -K*(dpwdn+g);
r_i.accumulate(lfsv_i,0, q*volume);
r_o.accumulate(lfsv_o,0,-q*volume);
}
template<typename IG, typename LFSU, typename X,
typename LFSV, typename R>
const LFSU& lfsu_i, const X& x_i,
const LFSV& lfsv_i, R& r_i) const
{
const auto& volume =
ig.geometry().volume();
const auto& global =
ig.geometry().center();
{
if ( (global[1]<(0.9*param.getHeight()) ) && (global[1]>(0.5*param.getHeight()) ) )
{
r_i.accumulate(lfsv_i,0,-param.inflow(
ig,x_i)*volume);
}
}
else if (global[0]>=param.getWidth()-1
e-6)
{
auto pw = x_i(lfsu_i,0);
auto K = param.K(param.Sw(pw));
auto outflow = -K*(0.-pw);
r_i.accumulate(lfsv_i,0,std::max(0.,outflow)*volume);
}
}
};
template<typename Param>
{
Param& param;
using RF = typename Param::value_type;
public:
: param(param_)
{}
template<typename EG, typename LFSU, typename X,
typename LFSV, typename R>
void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r) const
{
const auto volume = eg.geometry().volume();
auto pw = x(lfsu,0);
auto Sw = param.Sw(pw);
auto phi = param.getphi();
r.accumulate(lfsv,0, phi*Sw*volume );
}
};
template<typename Param, typename GFSF, typename ZF, typename GFSC, typename ZC, bool convection>
{
Param& param;
using RF = typename Param::value_type;
std::shared_ptr<ZF> dataf;
std::shared_ptr<const GFSF> pgfsf;
mutable LFSF lfsf;
mutable LFSFCache lfsf_cache;
std::shared_ptr<ZC> datac;
std::shared_ptr<const GFSC> pgfsc;
mutable LFSC lfsc;
mutable LFSCCache lfsc_cache;
public:
: param(param_)
, dataf(stackobject_to_shared_ptr(zf_))
, pgfsf(stackobject_to_shared_ptr(gfsf_))
, lfsf(pgfsf)
, lfsf_cache(lfsf)
, datac(stackobject_to_shared_ptr(zc_))
, pgfsc(stackobject_to_shared_ptr(gfsc_))
, lfsc(pgfsc)
, lfsc_cache(lfsc)
{}
{
param.setTime(t);
}
{
param.adjustTimeStep(coef);
}
template<typename IG, typename LFSU, typename X,
typename LFSV, typename R>
const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
R& r_i, R& r_o) const
{
using namespace Dune::TypeTree::Indices;
auto lfsv_i0 = lfsv_i.child(_0);
auto lfsv_i1 = lfsv_i.child(_1);
auto lfsu_i0 = lfsu_i.child(_0);
auto lfsu_i1 = lfsu_i.child(_1);
auto lfsv_o0 = lfsv_o.child(_0);
auto lfsv_o1 = lfsv_o.child(_1);
auto lfsu_o0 = lfsu_o.child(_0);
auto lfsu_o1 = lfsu_o.child(_1);
auto cell_inside =
ig.inside();
typename ZF::template LocalView<LFSFCache> p_view(*dataf);
lfsf.bind(cell_inside);
lfsf_cache.update();
std::vector<RF> pw(lfsf.size());
p_view.bind(lfsf_cache);
p_view.read(pw);
p_view.unbind();
auto cell_outside =
ig.outside();
auto insidegeo = cell_inside.geometry();
auto outsidegeo = cell_outside.geometry();
auto inside_global = insidegeo.center();
auto outside_global = outsidegeo.center();
inside_global -= outside_global;
auto distance = inside_global.two_norm();
auto facegeo =
ig.geometry();
auto volume = facegeo.volume();
auto pw_i = pw[0];
auto Sw_i = param.Sw(pw_i);
lfsf.bind(cell_outside);
lfsf_cache.update();
p_view.bind(lfsf_cache);
p_view.read(pw);
p_view.unbind();
auto pw_o = pw[0];
auto Sw_o = param.Sw(pw_o);
typename ZC::template LocalView<LFSCCache> c_view(*datac);
lfsc.bind(cell_inside);
lfsc_cache.update();
std::vector<RF> C_i(lfsc.size());
c_view.bind(lfsc_cache);
c_view.read(C_i);
c_view.unbind();
lfsc.bind(cell_outside);
lfsc_cache.update();
std::vector<RF> C_o(lfsc.size());
c_view.bind(lfsc_cache);
c_view.read(C_o);
c_view.unbind();
auto C0c_i = convection ? x_i(lfsu_i0,0) : C_i[0];
auto C1c_i = convection ? x_i(lfsu_i1,0) : C_i[1];
auto C0c_o = convection ? x_o(lfsu_o0,0) : C_o[0];
auto C1c_o = convection ? x_o(lfsu_o1,0) : C_o[1];
auto C0d_i = convection ? C_i[0] : x_i(lfsu_i0,0);
auto C1d_i = convection ? C_i[1] : x_i(lfsu_i1,0);
auto C0d_o = convection ? C_o[0] : x_o(lfsu_o0,0);
auto C1d_o = convection ? C_o[1] : x_o(lfsu_o1,0);
auto dpwdn = (pw_o-pw_i)/distance;
auto g = param.getgrav()*param.getrho()*inside_global[1]/distance;
auto K_i = param.K(Sw_i);
auto K_o = param.K(Sw_o);
auto K = param.harg(K_i,K_o);
auto q = -K*(dpwdn+g);
auto dC0dn = (C0d_o-C0d_i)/distance;
auto dC1dn = (C1d_o-C1d_i)/distance;
decltype(C0c_i) C0;
decltype(C1c_i) C1;
if ( q > 0 )
{
C0 = C0c_i;
C1 = C1c_i;
}
else
{
C0 = C0c_o;
C1 = C1c_o;
}
r_i.accumulate(lfsv_i0,0, q*C0*volume);
r_o.accumulate(lfsv_o0,0,-q*C0*volume);
r_i.accumulate(lfsv_i1,0, q*C1*volume);
r_o.accumulate(lfsv_o1,0,-q*C1*volume);
auto D0_i = param.getD0(Sw_i);
auto D0_o = param.getD0(Sw_o);
auto D0 = param.harg(D0_i,D0_o);
auto D1_i = param.getD1(Sw_i);
auto D1_o = param.getD1(Sw_o);
auto D1 = param.harg(D1_i,D1_o);
r_i.accumulate(lfsv_i0,0,-D0*dC0dn*volume);
r_o.accumulate(lfsv_o0,0, D0*dC0dn*volume);
r_i.accumulate(lfsv_i1,0,-D1*dC1dn*volume);
r_o.accumulate(lfsv_o1,0, D1*dC1dn*volume);
}
template<typename IG, typename LFSU, typename X,
typename LFSV, typename R>
const LFSU& lfsu_i, const X& x_i,
const LFSV& lfsv_i, R& r_i) const
{
using namespace Dune::TypeTree::Indices;
auto lfsv_i0 = lfsv_i.child(_0);
auto lfsv_i1 = lfsv_i.child(_1);
auto lfsu_i0 = lfsu_i.child(_0);
auto lfsu_i1 = lfsu_i.child(_1);
const auto& volume =
ig.geometry().volume();
const auto& global =
ig.geometry().center();
const auto& inside_cell =
ig.inside();
{
if (global[1]<0.9*param.getHeight() && global[1]>0.5*param.getHeight())
{
r_i.accumulate(lfsv_i0,0,-param.inflow(
ig,x_i)*param.inflowC0(
ig,x_i)*volume);
r_i.accumulate(lfsv_i1,0,-param.inflow(
ig,x_i)*param.inflowC1(
ig,x_i)*volume);
}
}
else if (global[0]>=param.getWidth()-1
e-6)
{
typename ZF::template LocalView<LFSFCache> p_view(*dataf);
lfsf.bind(inside_cell);
lfsf_cache.update();
std::vector<RF> pw(lfsf.size());
p_view.bind(lfsf_cache);
p_view.read(pw);
p_view.unbind();
auto K = param.K(param.Sw(pw[0]));
auto outflow = -K*(0.-pw[0]);
auto C0 = x_i(lfsv_i0,0);
auto C1 = x_i(lfsv_i1,0);
r_i.accumulate(lfsv_i0,0,std::max(0.,outflow)*C0*volume);
r_i.accumulate(lfsv_i1,0,std::max(0.,outflow)*C1*volume);
}
}
};
template<typename Param, typename GFSF, typename ZF, bool convection>
{
Param& param;
using RF = typename Param::value_type;
std::shared_ptr<ZF> dataf;
std::shared_ptr<ZF> datafold;
std::shared_ptr<const GFSF> pgfsf;
mutable LFSF lfsf;
mutable LFSFCache lfsf_cache;
public:
: param(param_)
, dataf(stackobject_to_shared_ptr(zf_))
, datafold(stackobject_to_shared_ptr(zfold_))
, pgfsf(stackobject_to_shared_ptr(gfsf_))
, lfsf(pgfsf)
, lfsf_cache(lfsf)
{}
template<typename EG, typename LFSU, typename X,
typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{
using namespace Dune::TypeTree::Indices;
const auto& lfsu_0 = lfsu.child(_0);
const auto& lfsu_1 = lfsu.child(_1);
const auto& lfsv_0 = lfsv.child(_0);
const auto& lfsv_1 = lfsv.child(_1);
const auto& inside_cell = eg.entity();
bool new_or_old = param.getTime() > param.getOtherTime();
typename ZF::template LocalView<LFSFCache> p_view(new_or_old ? *dataf : *datafold);
lfsf.bind(inside_cell);
lfsf_cache.update();
std::vector<RF> pw(lfsf.size());
p_view.bind(lfsf_cache);
p_view.read(pw);
p_view.unbind();
auto Sw = param.Sw(pw[0]);
auto C0 = x(lfsu_0,0);
auto C1 = x(lfsu_1,0);
auto phi = param.getphi();
const auto& volume = eg.geometry().volume();
r.accumulate(lfsv_0,0, phi*Sw*C0*volume);
r.accumulate(lfsv_1,0, phi*Sw*C1*volume);
}
};
template<typename GV, typename FEM, typename Param>
void driver (
const GV& gv,
const FEM& fem, Param& param)
{
using RF = typename Param::value_type;
GFS gfs(gv,fem);
gfs.name("pressure");
GFSC gfsc(gfs);
using namespace Dune::TypeTree::Indices;
gfsc.child(_0).name("C0");
gfsc.child(_1).name("C1");
Z z(gfs);
ZC zc(gfsc);
Z zfnew(z);
ZC zcnew(zc);
Z zfdata(z);
ZC zcdata(zc);
ZC zcstep(zc);
ZC zchalfstep(zc);
RF time = 0.0;
auto glambdaf = [&](
const auto&
e,
const auto& x)
{
RF gf;
return gf;
};
auto glambdac = [&] (
const auto&
e,
const auto& x)
{
Dune::FieldVector<RF,2> gc;
return gc;
};
using CF = typename GFS::template ConstraintsContainer<RF>::Type;
using CC = typename GFSC::template ConstraintsContainer<RF>::Type;
CF cf;
CC cc;
MBE mbe((int)10);
LOPF lopf(param);
GOF0 gof0(gfs,cf,gfs,cf,lopf,mbe);
TLOPF tlopf(param);
GOF1 gof1(gfs,cf,gfs,cf,tlopf,mbe);
IGOF igof(gof0,gof1);
constexpr bool convection{true};
LOPC lopc(param,gfs,zfdata,gfsc,zcdata);
GOC0 goc0(gfsc,cc,gfsc,cc,lopc,mbe);
TLOPC tlopc(param,gfs,zfdata,z);
GOC1 goc1(gfsc,cc,gfsc,cc,tlopc,mbe);
IGOC igoc(goc0,goc1);
constexpr bool diffusion{false};
LOPD lopd(param,gfs,zfdata,gfsc,zcdata);
GOD0 god0(gfsc,cc,gfsc,cc,lopd,mbe);
TLOPD tlopd(param,gfs,zfdata,z);
GOD1 god1(gfsc,cc,gfsc,cc,tlopd,mbe);
IGOD igod(god0,god1);
int verbose=0;
if (gfs.gridView().comm().rank()==0) verbose=1;
LSF lsf(gfs,cf,100,5,verbose);
LSC lsc(gfsc,cc,100,5,verbose);
constexpr int newtonMaxIt{10};
PDESOLVERF pdesolverf(igof,lsf);
Dune::ParameterTree newtonParam;
newtonParam["reassemble_threshold"] = "0.0";
newtonParam["verbosity"] = "2";
newtonParam["reduction"] = "1e-8";
newtonParam["min_linear_reduction"] = "1e-4";
newtonParam["terminate.max_iterations"] = std::to_string(newtonMaxIt);
newtonParam["line_search.line_search_max_iterations"] = "10";
pdesolverf.setParameters(newtonParam);
PDESOLVERC pdesolverc(igoc,lsc);
pdesolverc.setParameters(newtonParam);
PDESOLVERD pdesolverd(igod,lsc);
pdesolverd.setParameters(newtonParam);
OSMF osmf(*pmethod,igof,pdesolverf);
osmf.setVerbosityLevel(2);
OSMC osmc(*pmethod,igoc,pdesolverc);
osmc.setVerbosityLevel(2);
OSMD osmd(*pmethod,igod,pdesolverd);
osmd.setVerbosityLevel(2);
GFSC_Sub0 gfsc_sub0(gfsc);
SubC0 subC0(gfsc_sub0,zcstep);
SubC0 subC0new(gfsc_sub0,zcnew);
ErrC0 errC0(subC0,subC0new);
GFSC_Sub1 gfsc_sub1(gfsc);
SubC1 subC1(gfsc_sub1,zcstep);
SubC1 subC1new(gfsc_sub1,zcnew);
ErrC1 errC1(subC1,subC1new);
int subsampling=1;
VTKWRITER vtkwriter(gv,Dune::refinementIntervals(subsampling));
std::string filename="recipe-operator-splitting_parallel";
struct stat stf;
if( stat( filename.c_str(), &stf ) != 0 )
{
int statf = 0;
statf = mkdir(filename.c_str(),S_IRWXU|S_IRWXG|S_IRWXO);
if( statf != 0 && statf != -1)
std::cout << "Error: Cannot create directory "
<< filename << std::endl;
}
VTKSEQUENCEWRITER vtkSequenceWriter(
std::make_shared<VTKWRITER>(vtkwriter),filename,filename,"");
vtkSequenceWriter.write(0.0,Dune::VTK::appendedraw);
RF T = param.getT();
param.setTime(time);
param.setOtherTime(time);
RF lastvtk = 0.;
{
try
{
unsigned int iterations{0};
unsigned int split_iter{0};
bool continue_opsplit{true};
if (verbose>=1)
std::cout << std::endl << std::endl << "flow part:\n";
osmf.apply(time,param.getTimeStep(),z,gf,zfnew);
zfdata = zfnew;
iterations=pdesolverf.result().iterations;
do
{
if (verbose>=1)
std::cout << std::endl << std::endl << "contaminant convection part 1/2:\n";
osmc.apply(time,param.getTimeStep()/2,zc,gc,zcnew);
zcdata = zcnew;
zchalfstep = zcnew;
iterations=std::max(iterations,pdesolverc.result().iterations);
if (verbose>=1)
std::cout << std::endl << "contaminant diffusion part:\n";
osmd.apply(time,param.getTimeStep(),zc,gc,zcnew);
zcdata = zcnew;
iterations=std::max(iterations,pdesolverd.result().iterations);
if (verbose>=1)
std::cout << std::endl << "contaminant convection part 2/2:\n";
osmc.apply(time+param.getTimeStep()/2,param.getTimeStep()/2,zchalfstep,gc,zcnew);
zcdata = zcnew;
iterations=std::max(iterations,pdesolverc.result().iterations);
typename ErrC0::Traits::RangeType spliterrorcont0;
typename ErrC1::Traits::RangeType spliterrorcont1;
RF sperrC0 = spliterrorcont0.one_norm();
sperrC0 = gv.comm().sum(sperrC0);
RF sperrC1 = spliterrorcont1.one_norm();
sperrC1 = gv.comm().sum(sperrC1);
if (sperrC0 < 1
e-17 && sperrC1 < 1e-19 && split_iter > 0)
{
if (verbose>=1)
std::cout << "low spliterrors (" << sperrC0 << ", " << sperrC1 << "), go to next timestep" << std::endl;
continue_opsplit = false;
}
else
{
if (split_iter >= 5)
{
if (verbose>=1)
std::cout << "max operator splitting iteration number reached, current errors: " << sperrC0 << ", " << sperrC1 << ", reseting time step" << std::endl;
continue_opsplit = true;
}
else
{
if (verbose>=1)
std::cout << "going to " << split_iter+1 << ". operator splitting iteration, errors: " << sperrC0 << ", " << sperrC1 << std::endl;
continue_opsplit = true;
}
}
zcstep = zcnew;
zcdata = zcnew;
++split_iter;
} while (continue_opsplit);
zc=zcstep;
z=zfnew;
time+=param.getTimeStep();
param.setOtherTime(time);
if (time-lastvtk > time*0.01)
{
vtkSequenceWriter.write(time,Dune::VTK::appendedraw);
lastvtk = time;
}
if (iterations<(newtonMaxIt*4)/10)
{
param.adjustTimeStep(1.25);
}
else
if ((iterations>(newtonMaxIt*7)/10) && param.getTimeStep()>1
e-4)
{
iterations = 0;
param.adjustTimeStep(0.75);
}
}
catch(Dune::Exception&
e)
{
if (verbose==1)
{
std::cout <<
e.what() << std::endl;
}
if (param.getTimeStep()<1
e-4)
{
if (verbose==1)
{
std::cout << "too little time step" << std::endl;
}
vtkSequenceWriter.write(time+T,Dune::VTK::appendedraw);
throw;
}
else
{
param.adjustTimeStep(0.75);
zfnew = z;
zcnew = zc;
if (verbose==1)
{
std::cout << "Reducing TimeStep size to " << param.getTimeStep() << std::endl;
}
}
}
{
if (verbose==1)
{
std::cout <<
e.what() << std::endl;
}
if (param.getTimeStep()<1
e-4)
{
if (verbose==1)
{
std::cout << "too little time step" << std::endl;
}
vtkSequenceWriter.write(time+T,Dune::VTK::appendedraw);
throw;
}
else
{
param.adjustTimeStep(0.75);
zfnew = z;
zcnew = zc;
if (verbose==1)
{
std::cout << "Reducing TimeStep size to " << param.getTimeStep() << std::endl;
}
}
}
catch(...)
{
throw;
}
}
}
template <int dim>
{
public:
using iTuple = std::array<int,dim>;
void loadbalance (const iTuple& size, int P, iTuple& dims) const
{
for (
int j=0; j<
dim; ++j)
{
int nextP = P;
for (
int i=1; pow(i,
dim-j)<=P; ++i)
{
if(P%i == 0)
{
dims[j] = i;
nextP = P/i;
}
}
P = nextP;
}
auto greater = [](int a,int b) -> bool {return a>b;};
std::sort(dims.begin(),dims.end(),greater);
}
};
int main(
int argc,
char** argv)
{
try{
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
if(Dune::MPIHelper::isFake)
std::cout<< "This is a sequential program." << std::endl;
else
std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
<<" processes!"<<std::endl;
using RF = double;
using Grid = Dune::YaspGrid<dim>;
using DF = Grid::ctype;
using LType = Dune::FieldVector<DF,dim>;
LType L;
L[0] = 0.1;
L[1] = 0.1;
std::array<int,dim> N;
N[0] = 16;
N[1] = 16;
int overlap=1;
int refinement = 1;
std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>(new Grid(L,N,periodic,overlap,Dune::MPIHelper::getCollectiveCommunication(),&yp));
gridp->refineOptions(false);
gridp->globalRefine(refinement);
using GV = Grid::LeafGridView;
GV gv=gridp->leafGridView();
FEM fem(Dune::GeometryTypes::cube(
dim));
Param param(L);
}
catch (Dune::Exception &
e){
std::cerr <<
"Dune reported error: " <<
e << std::endl;
}
catch (...){
std::cerr << "Unknown exception thrown!" << std::endl;
}
}