dune-pdelab  2.7-git
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/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/pdelab.hh>
template<typename GV, typename RF>
{
public:
static constexpr bool permeabilityIsConstantPerCell()
{
return true;
}
A (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
{
for (std::size_t i=0; i<Traits::dimDomain; i++)
for (std::size_t j=0; j<Traits::dimDomain; j++)
I[i][j] = (i==j) ? 1 : 0;
return I;
}
b (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
{
typename Traits::RangeType v(0.0);
return v;
}
c (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
{
return 0.0;
}
f (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
{
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
{
}
g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal) const
{
return 0.0;
}
j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& xlocal) const
{
return 0.0;
}
};
// [Define function]
template <int dim>
: 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;
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)
{
try{
// 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};
typedef Dune::YaspGrid<dim> SquareGrid;
SquareGrid sgrid(L,N);
// [Mapping grid]
typedef typename Dune::GeometryGrid<SquareGrid,GridTransformation> Grid;
Grid grid(sgrid,gTrafo);
// define parameters
typedef double NumberType;
// need a grid in order to test grid functions
constexpr unsigned int degree = 1;
constexpr std::size_t nonzeros = Dune::power(2*degree+1,dim);
// make problem parameters
Problem problem;
BCType bctype(grid.leafGridView(),problem);
// Make FEM space
typedef typename Grid::ctype DF;
FEM fem(grid.leafGridView());
// make function space with constraints
typedef Dune::PDELab::GridFunctionSpace<Grid::LeafGridView,FEM,
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
X x(gfs,0.0);
G g(grid.leafGridView(),problem);
// 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));
vtkwriter.write("mesh",Dune::VTK::appendedraw);
return 0;
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch (...){
std::cerr << "Unknown exception thrown!" << std::endl;
return 1;
}
}
int main(int argc, char **argv)
[Define function]
Definition: recipe-geometry-grid.cc:147
static const int dim
Definition: adaptivity.hh:84
const Entity & e
Definition: localfunctionspace.hh:121
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
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
vtk::OutputCollector< VTKWriter, vtk::DGFTreeCommonData< GFS, X, Predicate > > addSolutionToVTKWriter(VTKWriter &vtk_writer, const GFS &gfs, const X &x, const NameGenerator &name_generator=vtk::defaultNameScheme(), const Predicate &predicate=Predicate())
Definition: vtk.hh:653
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:188
Definition: istl/descriptors.hh:48
Sequential conjugate gradient solver preconditioned with AMG smoothed by SSOR.
Definition: seqistlsolverbackend.hh:856
Dirichlet Constraints construction.
Definition: conforming.hh:38
A grid function space.
Definition: gridfunctionspace.hh:186
Definition: vtk.hh:30
Standard grid operator implementation.
Definition: gridoperator.hh:36
Definition: convectiondiffusionfem.hh:47
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
GV::Intersection IntersectionType
Definition: convectiondiffusionparameter.hh:106
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
Type
Definition: convectiondiffusionparameter.hh:113
@ Dirichlet
Definition: convectiondiffusionparameter.hh:113
Definition: convectiondiffusionparameter.hh:221
Definition: convectiondiffusionparameter.hh:325
Definition: linearproblem.hh:45
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
Dune::PDELab::ConvectionDiffusionParameterTraits< GV, RF > Traits
Definition: recipe-geometry-grid.cc:48
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
Definition: recipe-geometry-grid.cc:125
Base ::DomainVector DomainVector
Definition: recipe-geometry-grid.cc:130
Base ::RangeVector RangeVector
Definition: recipe-geometry-grid.cc:131
void evaluate(const DomainVector &x, RangeVector &y) const
Definition: recipe-geometry-grid.cc:135
GridTransformation()
Definition: recipe-geometry-grid.cc:133