Commit 0639b731 authored by Linus Seelinger's avatar Linus Seelinger

Sync with exercise

parent 99e16405
......@@ -17,7 +17,6 @@
#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/static_assert.hh>
#include<dune/common/timer.hh>
#include<dune/grid/io/file/gmshreader.hh>
......@@ -40,12 +39,9 @@
#include<dune/istl/solvers.hh>
#include<dune/istl/superlu.hh>
#include<dune/pdelab/backend/istlmatrixbackend.hh>
#include<dune/pdelab/backend/istlsolverbackend.hh>
#include<dune/pdelab/backend/istlvectorbackend.hh>
#include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/constraints/constraints.hh>
#include<dune/pdelab/constraints/conforming.hh>
#include<dune/pdelab/finiteelementmap/p0fem.hh>
#include<dune/pdelab/finiteelementmap/pkfem.hh>
......
......@@ -8,7 +8,7 @@ template<class GV> void example02_Q1 (const GV& gv)
typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,Real,1> FEM;
FEM fem(gv);
typedef Dune::PDELab::ConformingDirichletConstraints CONSTRAINTS; // constraints class
typedef Dune::PDELab::ISTLVectorBackend<> VBE;
typedef Dune::PDELab::istl::VectorBackend<> VBE;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CONSTRAINTS,VBE> GFS;
GFS gfs(gv,fem);
BCTypeParam bctype; // boundary condition type
......@@ -20,7 +20,7 @@ template<class GV> void example02_Q1 (const GV& gv)
// <<<3>>> Make grid operator
typedef Example02LocalOperator<BCTypeParam> LOP; // operator including boundary
LOP lop( bctype );
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::istl::BCRSMatrixBackend<std::size_t> MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
......@@ -29,7 +29,7 @@ template<class GV> void example02_Q1 (const GV& gv)
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
GO go(gfs,cc,gfs,cc,lop,MBE(9));
// <<<4>>> Make FE function extending Dirichlet boundary conditions
typedef typename GO::Traits::Domain U;
......@@ -51,6 +51,7 @@ template<class GV> void example02_Q1 (const GV& gv)
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
typedef Dune::PDELab::VTKGridFunctionAdapter<DGF> VTKADAPTER;
vtkwriter.addVertexData(std::make_shared<VTKADAPTER>(udgf,"solution"));
vtkwriter.write("example02_Q1",Dune::VTK::appendedraw);
}
#include<dune/common/fvector.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/constraints/constraintsparameters.hh>
//! \brief Parameter class selecting boundary conditions
class BCTypeParam
......
......@@ -68,12 +68,11 @@ public:
rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
// loop over quadrature points
for (typename Dune::QuadratureRule<DF,dim>::const_iterator
it=rule.begin(); it!=rule.end(); ++it)
for (const auto& ip : rule)
{
// evaluate basis functions on reference element
std::vector<Range> phi(lfsu.size());
lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
// compute u at integration point
RF u=0.0;
......@@ -82,11 +81,11 @@ public:
// evaluate gradient of basis functions on reference element
std::vector<Jacobian> js(lfsu.size());
lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
// transform gradients from reference element to real element
const Dune::FieldMatrix<DF,dimw,dim>
&jac = eg.geometry().jacobianInverseTransposed(it->position());
&jac = eg.geometry().jacobianInverseTransposed(ip.position());
std::vector<Gradient> gradphi(lfsu.size());
for (size_type i=0; i<lfsu.size(); i++)
jac.mv(js[i][0],gradphi[i]);
......@@ -98,12 +97,12 @@ public:
// evaluate parameters;
// Dune::FieldVector<RF,dim>
// globalpos = eg.geometry().global(it->position());
// globalpos = eg.geometry().global(ip.position());
RF f = 0;
RF a = 0;
// integrate grad u * grad phi_i + a*u*phi_i - f phi_i
RF factor = it->weight()*eg.geometry().integrationElement(it->position());
RF factor = ip.weight()*eg.geometry().integrationElement(ip.position());
for (size_type i=0; i<lfsu.size(); ++i)
r.accumulate(lfsu, i, (gradu*gradphi[i] + a*u*phi[i] - f*phi[i]) * factor);
}
......@@ -137,15 +136,14 @@ public:
rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
// loop over quadrature points and integrate normal flux
for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin();
it!=rule.end(); ++it)
for (const auto& ip : rule)
{
// skip rest if we are on Dirichlet boundary
if ( bctype.isDirichlet( ig, it->position() ) )
if ( bctype.isDirichlet( ig, ip.position() ) )
continue;
// position of quadrature point in local coordinates of element
Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
// evaluate basis functions at integration point
std::vector<Range> phi(lfsu_s.size());
......@@ -158,7 +156,7 @@ public:
// evaluate flux boundary condition
Dune::FieldVector<RF,dim>
globalpos = ig.geometry().global(it->position());
globalpos = ig.geometry().global(ip.position());
RF j;
if (globalpos[1]<0.5)
j = 1.0;
......@@ -166,7 +164,7 @@ public:
j = -1.0; // some outflow
// integrate j
RF factor = it->weight()*ig.geometry().integrationElement(it->position());
RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
for (size_type i=0; i<lfsu_s.size(); ++i)
r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
}
......
......@@ -17,7 +17,6 @@
#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/static_assert.hh>
#include<dune/common/timer.hh>
#include<dune/grid/io/file/gmshreader.hh>
......@@ -43,12 +42,9 @@
#include<dune/istl/solvers.hh>
#include<dune/istl/superlu.hh>
#include<dune/pdelab/backend/istlmatrixbackend.hh>
#include<dune/pdelab/backend/istlsolverbackend.hh>
#include<dune/pdelab/backend/istlvectorbackend.hh>
#include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/constraints/constraints.hh>
#include<dune/pdelab/constraints/conforming.hh>
#include<dune/pdelab/finiteelementmap/p0fem.hh>
#include<dune/pdelab/finiteelementmap/pkfem.hh>
......
......@@ -47,8 +47,7 @@ void sdsolver( U& u, const GV& gv, const GFS& gfs, const CC& cc,
// <<<3>>> Make grid operator space
typedef Example02bLocalOperator<BCTypeParam,BETA,Eps> LOP; // operator including boundary
LOP lop( bctype, beta, epsilon, method );
//typedef typename U::Backend::MatrixBackend MBE;
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::istl::BCRSMatrixBackend<std::size_t> MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
......@@ -57,7 +56,7 @@ void sdsolver( U& u, const GV& gv, const GFS& gfs, const CC& cc,
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
GO go(gfs,cc,gfs,cc,lop,MBE(9));
// <<<5>>> Select a linear solver backend
// typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR LS;
......@@ -74,7 +73,8 @@ void sdsolver( U& u, const GV& gv, const GFS& gfs, const CC& cc,
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
typedef Dune::PDELab::VTKGridFunctionAdapter<DGF> VTKADAPTER;
vtkwriter.addVertexData(std::make_shared<VTKADAPTER>(udgf,"solution"));
vtkwriter.write(outputname,Dune::VTK::appendedraw);
std::cout << "Output of 2b: " << methodname << " --> "
<< outputname << ".vtu" << std::endl << std::endl;
......@@ -91,7 +91,7 @@ void example02b_Q1 (const GV& gv, const VEC& dx, Eps epsilon)
typedef Dune::PDELab::QkLocalFiniteElementMap<GV,Coord,Real,1> FEM;
FEM fem(gv);
typedef Dune::PDELab::ConformingDirichletConstraints CONSTRAINTS; // constraints class
typedef Dune::PDELab::ISTLVectorBackend<> VBE;
typedef Dune::PDELab::istl::VectorBackend<> VBE;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CONSTRAINTS,VBE> GFS;
GFS gfs(gv,fem);
......
#include<dune/common/fvector.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/constraints/constraintsparameters.hh>
//! \brief Parameter class selecting boundary conditions
template<class BETA>
......@@ -38,7 +37,7 @@ public:
return false;
typename BETA::Traits::RangeType beta;
betaF.evaluate( *intersection.inside(),
betaF.evaluate( intersection.inside(),
intersection.geometryInInside().global(coord),
beta );
......
......@@ -85,12 +85,11 @@ public:
rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
// loop over quadrature points
for (typename Dune::QuadratureRule<DF,dim>::const_iterator
it=rule.begin(); it!=rule.end(); ++it)
for (const auto& ip : rule)
{
// evaluate basis functions on reference element
std::vector<Range> phi(lfsu.size());
lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
// compute u at integration point
RF u=0.0;
......@@ -99,11 +98,11 @@ public:
// evaluate gradient of basis functions on reference element
std::vector<Jacobian> js(lfsu.size());
lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
// transform gradients from reference element to real element
const Dune::FieldMatrix<DF,dimw,dim>
&jac = eg.geometry().jacobianInverseTransposed(it->position());
&jac = eg.geometry().jacobianInverseTransposed(ip.position());
std::vector<Gradient> gradphi(lfsu.size());
for (size_type i=0; i<lfsu.size(); i++)
jac.mv(js[i][0],gradphi[i]);
......@@ -115,7 +114,7 @@ public:
// evaluate parameters;
// Dune::FieldVector<RF,dim>
// globalpos = eg.geometry().global(it->position());
// globalpos = eg.geometry().global(ip.position());
RF f = 0;
// get the vector-field beta:
......@@ -131,7 +130,7 @@ public:
// Diffusion Term: eps * grad u * grad phi_i
// Convection Term: - u * beta * grad phi_i
// Source Term: - f phi_i
RF factor = it->weight()*eg.geometry().integrationElement(it->position());
RF factor = ip.weight()*eg.geometry().integrationElement(ip.position());
for (size_type i=0; i<lfsu.size(); ++i)
r.accumulate(lfsu, i, factor * ( epsilon * (gradu*gradphi[i])
- u * (beta*gradphi[i])
......@@ -197,22 +196,20 @@ public:
Dune::FieldVector<DF,dim> facecenterinelement =
ig.geometryInInside().global( facecenterlocal );
typedef typename IG::EntityPointer CellEntityPointer;
const CellEntityPointer eg = ig.inside();
const auto& eg = ig.inside();
typename BETA::Traits::RangeType beta;
betaF.evaluate( *eg, facecenterinelement, beta );
betaF.evaluate( eg, facecenterinelement, beta );
// loop over quadrature points and integrate normal flux
for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin();
it!=rule.end(); ++it)
for (const auto& ip : rule)
{
// skip rest if we are on Dirichlet boundary
if ( bctype.isDirichlet( ig, it->position() ) )
if ( bctype.isDirichlet( ig, ip.position() ) )
continue;
// position of quadrature point in local coordinates of element
Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
// evaluate basis functions at integration point
std::vector<Range> phi(lfsu_s.size());
......@@ -224,7 +221,7 @@ public:
// transform gradients from reference element to real element
const Dune::FieldMatrix<DF,dimw,dim>
jac = eg->geometry().jacobianInverseTransposed(local);
jac = eg.geometry().jacobianInverseTransposed(local);
std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu_s.size());
for (size_type i=0; i<lfsu_s.size(); i++)
jac.mv(js[i][0],gradphi[i]);
......@@ -239,9 +236,9 @@ public:
for (size_type i=0; i<lfsu_s.size(); ++i)
u += x_s(lfsu_s,i)*phi[i];
RF factor = it->weight()*ig.geometry().integrationElement(it->position());
Dune::FieldVector<DF, dimw> normal = ig.unitOuterNormal(it->position());
if ( bctype.isNeumann( ig, it->position() ) ) {
RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
Dune::FieldVector<DF, dimw> normal = ig.unitOuterNormal(ip.position());
if ( bctype.isNeumann( ig, ip.position() ) ) {
// evaluate Neuman boundary condition
// prescribe total flux j
// \{ -eps \nabla u + u * beta \} \cdot n = j
......@@ -251,7 +248,7 @@ public:
r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
}
if ( bctype.isOutflow( ig, it->position() ) ) {
if ( bctype.isOutflow( ig, ip.position() ) ) {
// evaluate outflow boundary condition
// prescribe diffusive flux j_eps
// - eps * \nabla u \cdot n = j_eps
......
......@@ -17,7 +17,6 @@
#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/static_assert.hh>
#include<dune/common/timer.hh>
#include<dune/grid/io/file/gmshreader.hh>
......@@ -43,12 +42,9 @@
#include<dune/istl/solvers.hh>
#include<dune/istl/superlu.hh>
#include<dune/pdelab/backend/istlmatrixbackend.hh>
#include<dune/pdelab/backend/istlsolverbackend.hh>
#include<dune/pdelab/backend/istlvectorbackend.hh>
#include<dune/pdelab/backend/istl.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/vtkexport.hh>
#include<dune/pdelab/constraints/constraints.hh>
#include<dune/pdelab/constraints/conforming.hh>
#include<dune/pdelab/finiteelementmap/p0fem.hh>
#include<dune/pdelab/finiteelementmap/pkfem.hh>
......
......@@ -9,7 +9,7 @@ void example02c_P1 (const GV& gv, const std::string& vtkprefix)
typedef Dune::PDELab::PkLocalFiniteElementMap<GV,Coord,Real,1> FEM;
FEM fem(gv);
typedef Dune::PDELab::ConformingDirichletConstraints CONSTRAINTS; // constraints class
typedef Dune::PDELab::ISTLVectorBackend<> VBE;
typedef Dune::PDELab::istl::VectorBackend<> VBE;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CONSTRAINTS,VBE> GFS;
GFS gfs(gv,fem);
BCTypeParam bctype; // boundary condition type
......@@ -22,8 +22,7 @@ void example02c_P1 (const GV& gv, const std::string& vtkprefix)
// <<<3>>> Make grid operator
typedef Example02LocalOperator<BCTypeParam> LOP; // operator including boundary
LOP lop( bctype );
//typedef VBE::MatrixBackend MBE;
typedef Dune::PDELab::ISTLMatrixBackend MBE;
typedef Dune::PDELab::istl::BCRSMatrixBackend<> MBE;
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
......@@ -32,7 +31,7 @@ void example02c_P1 (const GV& gv, const std::string& vtkprefix)
Real,Real,Real, /* field types for domain, range and jacobian */
CC,CC /* constraints transformation for ansatz and test space */
> GO;
GO go(gfs,cc,gfs,cc,lop);
GO go(gfs,cc,gfs,cc,lop,MBE(5));
// <<<4>>> Make FE function extending Dirichlet boundary conditions
typedef typename GO::Traits::Domain U;
......@@ -54,6 +53,7 @@ void example02c_P1 (const GV& gv, const std::string& vtkprefix)
typedef Dune::PDELab::DiscreteGridFunction<GFS,U> DGF;
DGF udgf(gfs,u);
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
vtkwriter.addVertexData(new Dune::PDELab::VTKGridFunctionAdapter<DGF>(udgf,"solution"));
typedef Dune::PDELab::VTKGridFunctionAdapter<DGF> VTKADAPTER;
vtkwriter.addVertexData(std::make_shared<VTKADAPTER>(udgf,"solution"));
vtkwriter.write(vtkprefix,Dune::VTK::appendedraw);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment