See explanation at Assembling a linear system from a PDE
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/filledarray.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
template<typename GV, typename RF>
{
public:
{
return true;
}
{
I[i][
j] = (i==
j) ? 1 : 0;
return I;
}
{
return v;
}
{
return 0.0;
}
{
return 0.0;
}
BCType
{
}
{
if (x[0] < 0.5)
return 1.0;
return 0.0;
}
{
return 0.0;
}
{
return 0.0;
}
};
int main(
int argc,
char **argv)
{
try {
Dune::MPIHelper::instance(argc,argv);
typedef double NumberType;
constexpr
unsigned int dim = 1;
constexpr unsigned int degree = 1;
constexpr std::size_t nonzeros = Dune::power(2*degree+1,
dim);
Dune::FieldVector<NumberType,dim> L(1.0);
std::array<int,dim> N(Dune::filledArray<dim,int>(5));
typedef Dune::YaspGrid<dim> Grid;
Grid grid(L,N);
typedef Dune::YaspGrid<dim> GM;
Problem problem;
BCType bctype(grid.leafGridView(),problem);
typedef typename GM::ctype DF;
FEM fem(grid.leafGridView());
> GFS;
GFS gfs(grid.leafGridView(),fem);
typedef typename GFS::template ConstraintsContainer<NumberType>::Type CC;
CC cc;
LOP lop(problem);
auto go = GO(gfs,cc,gfs,cc,lop,MBE(nonzeros));
X x(gfs,0.0);
G g(grid.leafGridView(),problem);
X d(gfs,0.0);
go.residual(x,d);
Dune::printvector(std::cout,
native(d),
"Residual vector",
"");
typedef GO::Jacobian M;
M A(go);
go.jacobian(x,A);
Dune::printmatrix(std::cout,
native(A),
"Stiffness matrix",
"");
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;
}
}
Definition: convectiondiffusionparameter.hh:217
Standard grid operator implementation.
Definition: gridoperator.hh:35
Definition: istl/descriptors.hh:47
Traits::PermTensorType A(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
tensor diffusion coefficient
Definition: recipe-geometry-grid.cc:58
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Dirichlet Constraints construction.
Definition: conforming.hh:36
BCType bctype(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
boundary condition type function
Definition: recipe-geometry-grid.cc:95
Definition: recipe-geometry-grid.cc:43
Definition: convectiondiffusionparameter.hh:320
Traits::RangeFieldType f(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
source term
Definition: recipe-geometry-grid.cc:84
Type
Definition: convectiondiffusionparameter.hh:113
@ Dirichlet
Definition: convectiondiffusionparameter.hh:113
Traits class for convection diffusion parameters.
Definition: convectiondiffusionparameter.hh:75
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: convectiondiffusionparameter.hh:99
Traits::RangeFieldType o(const typename Traits::IntersectionType &is, const typename Traits::IntersectionDomainType &xlocal) const
outflow boundary condition
Definition: recipe-geometry-grid.cc:116
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: convectiondiffusionparameter.hh:102
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
Dirichlet boundary condition value.
Definition: recipe-geometry-grid.cc:102
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:187
Definition: convectiondiffusionfem.hh:39
int main(int argc, char **argv)
Definition: recipe-linear-system-assembly.cc:154
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: convectiondiffusionparameter.hh:90
const Entity & e
Definition: localfunctionspace.hh:121
RF RangeFieldType
Export type for range field.
Definition: convectiondiffusionparameter.hh:96
A grid function space.
Definition: gridfunctionspace.hh:178
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
GV::Intersection IntersectionType
Definition: convectiondiffusionparameter.hh:106
Dune::PDELab::ConvectionDiffusionParameterTraits< GV, RF > Traits
Definition: recipe-geometry-grid.cc:48
@ dimDomain
dimension of the domain
Definition: convectiondiffusionparameter.hh:83
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 c(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
sink term
Definition: recipe-geometry-grid.cc:77
Traits::RangeType b(const typename Traits::ElementType &e, const typename Traits::DomainType &xlocal) const
velocity field
Definition: recipe-geometry-grid.cc:69
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: convectiondiffusionparameter.hh:93
GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: convectiondiffusionparameter.hh:105
static const int dim
Definition: adaptivity.hh:84
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