DUNE PDELab (2.8)
recipe-geometry-grid.cc
See explanation at Transforming a cartesian mesh
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/geometrygrid/grid.hh>
#include <dune/pdelab.hh>
template<typename GV, typename RF>
{
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
public:
{
return true;
}
{
typename Traits::PermTensorType I;
return I;
}
{
typename Traits::RangeType v(0.0);
return v;
}
{
return 0.0;
}
typename Traits::RangeFieldType
{
return e.geometry().global(xlocal)[0] < 0.7 ? 5.0 : 1.0;
}
/* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet for Dirichlet boundary conditions
* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Neumann for flux boundary conditions
* return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Outflow for outflow boundary conditions
*/
BCType
bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
}
typename Traits::RangeFieldType
{
return 0.0;
}
typename Traits::RangeFieldType
j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
typename Traits::RangeFieldType
o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
};
// [Define function]
template <int dim>
class GridTransformation
: public Dune :: AnalyticalCoordFunction< double, dim, dim, GridTransformation <dim> >{
typedef GridTransformation This;
typedef Dune :: AnalyticalCoordFunction< double, dim, dim, This > Base;
public:
typedef typename Base :: DomainVector DomainVector;
typedef typename Base :: RangeVector RangeVector;
GridTransformation(){}
void evaluate(const DomainVector &x, RangeVector &y) const{
y = x;
if(x[0] < 0.8)
y[1] = (1.0 + 5.0/4.0 * (sin(M_PI/18.0) - 1.0) * x[0]) * (x[1] - 1.0);
else
y[1] = sin((x[0] - 0.6)/3.6 * M_PI) * (x[1] - 1.0);
if(x[0] > 3.8)
y[0] += 0.5*(x[0] - 3.8) * (1.0 - pow(x[1] - 1.0, 2.0));
}
};
int main(int argc, char** argv)
{
// Initialize Mpi
Dune::MPIHelper::instance(argc, argv);
// [Setting up grid]
const unsigned int dim = 2;
Dune::FieldVector<double,dim> L = {4.0,2.0};
std::array<int,dim> N ={64,32};
SquareGrid sgrid(L,N);
// [Mapping grid]
typedef GridTransformation<dim> GridTransformation;
GridTransformation gTrafo;
Grid grid(sgrid,gTrafo);
// define parameters
typedef double NumberType;
// need a grid in order to test grid functions
// make problem parameters
typedef GenericEllipticProblem<typename Grid::LeafGridView,NumberType> Problem;
Problem problem;
BCType bctype(grid.leafGridView(),problem);
// Make FEM space
typedef typename Grid::ctype DF;
FEM fem(grid.leafGridView());
// make function space with constraints
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed,1>> GFS;
GFS gfs(grid.leafGridView(),fem); gfs.name("solution");
// Make constraints
typedef typename GFS::template ConstraintsContainer<NumberType>::Type CC;
CC cc;
Dune::PDELab::constraints(bctype,gfs,cc);
// Set up local operator
LOP lop(problem);
// Set up grid operator
auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));
// Define solution
typedef Dune::PDELab::Backend::Vector<GFS,NumberType> X;
X x(gfs,0.0);
G g(grid.leafGridView(),problem);
Dune::PDELab::interpolate(g,gfs,x);
// Solve Poisson equation using AMG
LS ls(100,3);
SLP slp(go,ls,x,1e-10);
slp.apply(); // here all the work is done!
// Plot the mesh
Dune::SubsamplingVTKWriter<decltype(grid.leafGridView())> vtkwriter(grid.leafGridView(),Dune::refinementLevels(0));
Dune::PDELab::addSolutionToVTKWriter(vtkwriter,gfs,x);
}
static DUNE_EXPORT MPIHelper & instance(int &argc, char **&argv)
Get the singleton instance of the helper.
Definition: mpihelper.hh:240
Dirichlet Constraints construction.
Definition: conforming.hh:38
Definition: convectiondiffusionparameter.hh:221
Definition: convectiondiffusionparameter.hh:325
Definition: convectiondiffusionfem.hh:47
Sequential conjugate gradient solver preconditioned with AMG smoothed by SSOR.
Definition: seqistlsolverbackend.hh:856
Definition: qkfem.hh:20
Solve linear problems using a residual formulation.
Definition: linearproblem.hh:61
Writer for the output of subsampled grid functions in the vtk format.
Definition: subsamplingvtkwriter.hh:38
Definition: recipe-geometry-grid.cc:44
Traits::RangeFieldType o(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
outflow boundary condition
Definition: recipe-geometry-grid.cc:116
Traits::PermTensorType A(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
tensor diffusion coefficient
Definition: recipe-geometry-grid.cc:58
BCType bctype(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
boundary condition type function
Definition: recipe-geometry-grid.cc:95
static constexpr bool permeabilityIsConstantPerCell()
tensor diffusion constant per cell? return false if you want more than one evaluation of A per cell.
Definition: recipe-geometry-grid.cc:51
Traits::RangeType b(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
velocity field
Definition: recipe-geometry-grid.cc:69
Traits::RangeFieldType c(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
sink term
Definition: recipe-geometry-grid.cc:77
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
Dirichlet boundary condition value.
Definition: recipe-geometry-grid.cc:102
Traits::RangeFieldType j(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
flux boundary condition
Definition: recipe-geometry-grid.cc:109
Traits::RangeFieldType f(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
source term
Definition: recipe-geometry-grid.cc:84
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:76
RefinementIntervals refinementLevels(int levels)
Creates a RefinementIntervals object.
Definition: base.cc:115
Helpers for dealing with MPI.
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
Traits class for convection diffusion parameters.
Definition: convectiondiffusionparameter.hh:76
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: convectiondiffusionparameter.hh:93
@ dimDomain
dimension of the domain
Definition: convectiondiffusionparameter.hh:83
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: convectiondiffusionparameter.hh:102
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: convectiondiffusionparameter.hh:99
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: convectiondiffusionparameter.hh:105
RF RangeFieldType
Export type for range field.
Definition: convectiondiffusionparameter.hh:96
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: convectiondiffusionparameter.hh:90
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:188
Provides subsampled file i/o for the visualization toolkit.
![Creative Commons License](https://i.creativecommons.org/l/by-sa/4.0/80x15.png)