Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • pdelab/dune-pdelab-tutorials
  • arne.strehlow/dune-pdelab-tutorials
  • ranjeet.kumar/dune-pdelab-tutorials
3 results
Show changes
Commits on Source (5)
Hide whitespace changes
Inline Side-by-side
Showing
with 104 additions and 172 deletions
tutorial07/doc/tutorial07.tex
View file @ 0d03315b
......@@ -12,7 +12,7 @@
\usepackage{eurosym}
\usepackage{graphicx}
\usepackage{color}
\usepackage{listings,lstautogobble}
\usepackage{listings}
\lstset{language=C++, basicstyle=\ttfamily,
keywordstyle=\color{black}\bfseries, tabsize=4, stringstyle=\ttfamily,
commentstyle=\itshape, extendedchars=true, escapeinside={/*@}{@*/},
......@@ -905,8 +905,7 @@ We include our hyperbolic model and problem to solve
\lstinputlisting[firstline=60, lastline=67,
basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg},
autogobble=true]{../src/tutorial07-linearacoustics.cc}
backgroundcolor=\color{listingbg}]{../src/tutorial07-linearacoustics.cc}
Calling Problem and Model constructors and choose proper numerical flux
......@@ -1096,26 +1095,26 @@ and may now loop over the quadrature points.
\lstinputlisting[firstline=124, lastline=136,
basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg},autogobble=true]{../src/hyperbolicdg.hh}
backgroundcolor=\color{listingbg}]{../src/hyperbolicdg.hh}
\lstinputlisting[firstline=224, lastline=235,
basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg},autogobble=true]{../src/hyperbolicdg.hh}
backgroundcolor=\color{listingbg}]{../src/hyperbolicdg.hh}
\lstinputlisting[firstline=303, lastline=311,
basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg},autogobble=true]{../src/hyperbolicdg.hh}
backgroundcolor=\color{listingbg}]{../src/hyperbolicdg.hh}
\lstinputlisting[firstline=339, lastline=350,
basicstyle=\ttfamily\small,
frame=single,
backgroundcolor=\color{listingbg},autogobble=true]{../src/hyperbolicdg.hh}
backgroundcolor=\color{listingbg}]{../src/hyperbolicdg.hh}
\subsubsection*{\lstinline{jacobian_volume} method}
......
tutorial07/src/driver.hh
View file @ 0d03315b
......@@ -41,9 +41,9 @@ void driver (const GV& gv, const FEMDG& femdg, NUMFLUX& numflux, Dune::Parameter
std::cout << "degrees of freedom: " << gfs.globalSize() << std::endl;
// Make instationary grid operator
using LOP = Dune::PDELab::DGLinearHyperbolicSpatialOperator<NUMFLUX,FEMDG>;
using LOP = Dune::PDELab::DGHyperbolicSpatialOperator<NUMFLUX,FEMDG>;
LOP lop(numflux);
using TLOP = Dune::PDELab::DGLinearHyperbolicTemporalOperator<NUMFLUX,FEMDG>;
using TLOP = Dune::PDELab::DGHyperbolicTemporalOperator<NUMFLUX,FEMDG>;
TLOP tlop(numflux);
using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
......
tutorial07/src/hyperbolicdg.hh
View file @ 0d03315b
......@@ -34,13 +34,13 @@ namespace Dune {
\tparam FEM Finite Element Map needed to select the cache
*/
template<typename NUMFLUX, typename FEM>
class DGLinearHyperbolicSpatialOperator :
public NumericalJacobianApplyVolume<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianVolume<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianApplySkeleton<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianSkeleton<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianApplyBoundary<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianBoundary<DGLinearHyperbolicSpatialOperator<NUMFLUX,FEM> >,
class DGHyperbolicSpatialOperator :
public NumericalJacobianApplyVolume<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianVolume<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianApplySkeleton<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianSkeleton<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianApplyBoundary<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public NumericalJacobianBoundary<DGHyperbolicSpatialOperator<NUMFLUX,FEM> >,
public FullSkeletonPattern,
public FullVolumePattern,
public LocalOperatorDefaultFlags,
......@@ -62,7 +62,7 @@ namespace Dune {
enum { doLambdaVolume = true };
// ! constructor
DGLinearHyperbolicSpatialOperator (NUMFLUX& numflux_, int overintegration_=0)
DGHyperbolicSpatialOperator (NUMFLUX& numflux_, int overintegration_=0)
: numflux(numflux_), overintegration(overintegration_), cache(20)
{
}
......@@ -364,8 +364,8 @@ namespace Dune {
* \f}
*/
template<typename NUMFLUX, typename FEM>
class DGLinearHyperbolicTemporalOperator :
public NumericalJacobianApplyVolume<DGLinearHyperbolicTemporalOperator<NUMFLUX,FEM> >,
class DGHyperbolicTemporalOperator :
public NumericalJacobianApplyVolume<DGHyperbolicTemporalOperator<NUMFLUX,FEM> >,
public LocalOperatorDefaultFlags,
public InstationaryLocalOperatorDefaultMethods<typename NUMFLUX::RF>
{
......@@ -380,7 +380,7 @@ namespace Dune {
// residual assembly flags
enum { doAlphaVolume = true };
DGLinearHyperbolicTemporalOperator (NUMFLUX& numflux_, int overintegration_=0)
DGHyperbolicTemporalOperator (NUMFLUX& numflux_, int overintegration_=0)
: numflux(numflux_), overintegration(overintegration_), cache(20)
{}
......
tutorial07/src/linearacoustics.hh
View file @ 0d03315b
......@@ -10,83 +10,11 @@
\param RT matrix to be filled
*/
template<int dim, typename PROBLEM>
class Model ;
template<typename PROBLEM>
class Model<1,PROBLEM>
class Model
{
public:
static constexpr int dim = 1; // space dimension
static constexpr int m = dim+1; // system size
static constexpr int mplus = 1; // number of positive eigenvalues
static constexpr int mminus = 1; // number of negative eigenvalues
static constexpr int mstar = mplus+mminus; // number of nonzero eigenvalues
using RangeField = typename PROBLEM::RangeField;
Model (PROBLEM& p)
: problem(p)
{
}
template<typename E, typename X, typename T2, typename T3>
void eigenvectors (const E& e, const X& x,
const Dune::FieldVector<T2,dim>& n,
Dune::FieldMatrix<T3,m,m>& RT) const
{
auto c = 1.;
RT[0][0] = 1; RT[1][0] = c*n[0];
RT[0][1] = -1; RT[1][1] = c*n[0];
}
template<typename E, typename X, typename RF>
void diagonal (const E& e, const X& x, Dune::FieldMatrix<RF,m,m>& D) const
{
auto c = 1.; //problem.c(e,x);
D[0][0] = c; D[0][1] = 0.0;
D[1][0] = 0.0; D[1][1] = -c ;
}
template<typename RF>
static void max_eigenvalue (const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha)
{
alpha = 1.0;
}
//Flux function
template<typename E, typename X, typename RF>
void flux (const E& e, const X& x,
const Dune::FieldVector<RF,m>& u,
Dune::FieldMatrix<RF,m,dim>& F) const
{
auto c = 1.;//problem.c(e,x);
F[0][0] = u[1] ;
F[1][0] = c*c*u[0];
}
const PROBLEM& problem;
};//1d
template<typename PROBLEM>
class Model<2,PROBLEM>
{
public:
static constexpr int dim = 2; // space dimension
static constexpr int dim = PROBLEM::dim; // space dimension
static constexpr int m = dim+1; // system size
static constexpr int mplus = 1; // number of positive eigenvalues
static constexpr int mminus = 1; // number of negative eigenvalues
......@@ -108,9 +36,16 @@ public:
{
auto c = problem.c(e,x);
RT[0][0] = 1.0/c; RT[0][1] = -1.0/c; RT[0][2] = 0.0;
RT[1][0] = n[0]; RT[1][1] = n[0]; RT[1][2] = -n[1];
RT[2][0] = n[1]; RT[2][1] = n[1]; RT[2][2] = n[0];
//TODO find a way to write eigenvectors independently of dim
if (dim == 1) {
RT[0][0] = 1; RT[1][0] = c*n[0];
RT[0][1] = -1; RT[1][1] = c*n[0];
}
if (dim == 2) {
RT[0][0] = 1.0/c; RT[0][1] = -1.0/c; RT[0][2] = 0.0;
RT[1][0] = n[0]; RT[1][1] = n[0]; RT[1][2] = -n[1];
RT[2][0] = n[1]; RT[2][1] = n[1]; RT[2][2] = n[0];
}
}
/// tex: eigenvectors
......@@ -121,20 +56,25 @@ public:
{
auto c = problem.c(e,x);
D[0][0] = c; D[0][1] = 0.0; D[0][2] = 0.0;
D[1][0] = 0.0; D[1][1] = -c ; D[1][2] = 0.0;
D[2][0] = 0.0; D[2][1] = 0.0; D[2][2] = 0.0;
for (size_t i=0; i<m; i++)
for (size_t j=0; j<m; j++)
D[i][j] = 0.0;
D[0][0] = c;
D[1][1] = -c ;
}
/// tex: diagonal
template<typename RF>
static void max_eigenvalue (const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha)
template<typename E, typename X, typename RF>
void max_eigenvalue (const E& inside, const X& x_inside,
const E& outside, const X& x_outside,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha) const
{
alpha = 1.0;
alpha = std::max( problem.c(inside,x_inside),
problem.c(outside,x_outside) );
}
/// tex: flux
......@@ -146,13 +86,16 @@ public:
{
auto c = problem.c(e,x);
F[0][0] = u[1] ; F[0][1] = u[2];
F[1][0] = c*c*u[0]; F[1][1] = 0.0;
F[2][0] = 0.0 ; F[2][1] = c*c*u[0];
for (size_t i=0; i<dim; i++) {
F[0][i] = u[i+1];
F[i+1][i] = c*c*u[0];
}
}
/// tex: flux
const PROBLEM& problem;
};//2d
};// Model
#endif //ACOUSTICS_MODEL
tutorial07/src/linearacousticsproblem.hh
View file @ 0d03315b
#ifndef ACOUSTICS_RIEMANNPROBLEM
#define ACOUSTICS_RIEMANNPROBLEM
template<const int dim, typename GV, typename NUMBER>
template<typename GV, typename NUMBER>
class Problem
{
public:
......@@ -8,7 +8,7 @@ public:
using RangeField = NUMBER;
//problem specification depends on dimension
//static constexpr int dim = 2;
static constexpr int dim = GV::dimension;
static constexpr int m = dim+1;
using Range = Dune::FieldVector<NUMBER,m>;
......@@ -26,7 +26,7 @@ public:
int material (const E& e, const X& x) const
{
auto xglobal = e.geometry().center();
if (xglobal[1]>0.625)
if (xglobal[0]>1.625)
return 1;
else
return 2;
......@@ -39,20 +39,13 @@ public:
NUMBER c (const E& e, const X& x) const
{
auto xglobal = e.geometry().center();
if (xglobal[1]>0.625)
if (xglobal[0]>1.625)
return 0.33333;
else
return 1.0;
}
/// tex: speedofsound
//! Neumann boundary condition
template<typename I, typename X>
NUMBER j (const I& i, const X& x) const
{
return 0.0;
}
/// tex: bc
//! Boundary condition value - reflecting bc
template<typename I, typename X, typename R>
......@@ -60,8 +53,10 @@ public:
{
Range u(0.0);
u[0] = s[0];
u[1] = -s[1];
u[2] = -s[2];
for (size_t i=0; i<dim; i++)
u[i+1] = -s[i+1];
return u;
}
/// tex: bc
......
tutorial07/src/maxwellproblem.hh
View file @ 0d03315b
......@@ -49,18 +49,12 @@ public:
return 1.0;
}
//! Neumann boundary condition
template<typename I, typename X>
NUMBER j (const I& i, const X& x) const
{
return 0.0;
}
//! Boundary condition value - reflecting bc
//! Boundary condition
template<typename I, typename X, typename R>
Range g (const I& is, const X& x, const R& s) const
{
Range u(0.0);
//reflecting bc
// u[0] = -s[0];
// u[1] = -s[1];
// u[2] = -s[2];
......@@ -95,10 +89,10 @@ public:
auto c1 = std::cos(alpha*x);
auto st = std::sin(alpha*0.0);
auto ct = std::cos(alpha*0.0);
u[0] += 0; // E_x
u[0] += 0; // E_x
u[1] += -c1*st; // E_y
u[2] += s1*ct; // E_z
u[3] += 0; // H_x
u[3] += 0; // H_x
u[4] += c1*st; // H_y
u[5] += s1*ct; // H_z
......
tutorial07/src/numericalflux.hh
View file @ 0d03315b
......@@ -9,7 +9,7 @@ class LLFflux
public:
static constexpr int dim = MODEL::dim;
static constexpr int m = MODEL::Model::m;
static constexpr int m = MODEL::m;
using Model = MODEL;
using RF = typename MODEL::RangeField; // type for computations
......@@ -43,7 +43,7 @@ public:
//max eigenvalue
RF alpha(0.0);
MODEL::max_eigenvalue(u_s,u_n,n_F,alpha);
model().max_eigenvalue(inside,x_inside,outside,x_outside,u_s,u_n,n_F,alpha);
//add diffusion
for (size_t i =0 ; i<m;i++)
......@@ -69,7 +69,7 @@ class FluxVectorSplitting
public:
static constexpr int dim = MODEL::dim;
static constexpr int m = MODEL::Model::m;
static constexpr int m = MODEL::m;
using Model = MODEL;
using RF = typename MODEL::RangeField; // type for computations
......@@ -142,7 +142,7 @@ class VariableFluxVectorSplitting
public:
static constexpr int dim = MODEL::dim;
static constexpr int m = MODEL::Model::m;
static constexpr int m = MODEL::m;
static constexpr int mstar = MODEL::mstar;
using Model = MODEL;
......
tutorial07/src/shallowwater.hh
View file @ 0d03315b
......@@ -5,10 +5,14 @@
*/
template<int dim, typename PROBLEM>
class Model ;
class Model_ ;
//wrapper for dimension specialisation
template<typename PROBLEM>
using Model = Model_<PROBLEM::dim, PROBLEM>;
template<typename PROBLEM>
class Model<1,PROBLEM>
class Model_<1,PROBLEM>
{
public:
......@@ -17,16 +21,18 @@ public:
using RangeField = typename PROBLEM::RangeField;
Model (PROBLEM& p)
Model_ (PROBLEM& p)
: problem(p)
{
}
template<typename RF>
static void max_eigenvalue (const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha)
template<typename E, typename X, typename RF>
void max_eigenvalue (const E& inside, const X& x_inside,
const E& outside, const X& x_outside,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha) const
{
int g = 1;
......@@ -56,7 +62,7 @@ public:
template<typename PROBLEM>
class Model<2,PROBLEM>
class Model_<2,PROBLEM>
{
public:
......@@ -65,16 +71,18 @@ public:
using RangeField = typename PROBLEM::RangeField;
Model (PROBLEM& p)
Model_ (PROBLEM& p)
: problem(p)
{
}
template<typename RF>
static void max_eigenvalue (const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha)
template<typename E, typename X, typename RF>
void max_eigenvalue (const E& inside, const X& x_inside,
const E& outside, const X& x_outside,
const Dune::FieldVector<RF,m>& u_s,
const Dune::FieldVector<RF,m>& u_n,
const Dune::FieldVector<RF,dim>& n_F,
RF& alpha) const
{
int g = 1;
......
tutorial07/src/shallowwaterproblem.hh
View file @ 0d03315b
#ifndef SHALLOWWATER_RIEMANNPROBLEM
#define SHALLOWWATER_RIEMANNPROBLEM
template<const int dim, typename GV, typename NUMBER>
template<typename GV, typename NUMBER>
class Problem
{
public:
......@@ -8,6 +8,7 @@ public:
using RangeField = NUMBER;
//problem specification depends on dimension
static constexpr int dim = GV::dimension;
static constexpr int m = dim+1;
using Range = Dune::FieldVector<NUMBER,m>;
......@@ -17,14 +18,6 @@ public:
{
}
//! Neumann boundary condition
template<typename I, typename X>
NUMBER j (const I& i, const X& x) const
{
return 0.0;
}
//! Boundary condition value - reflecting bc
template<typename I, typename X, typename R>
Range g (const I& is, const X& x, const R& s) const
......
tutorial07/src/tutorial07-linearacoustics.cc
View file @ 0d03315b
......@@ -130,11 +130,11 @@ int main(int argc, char** argv)
GV gv = grid.leafGridView();
//create problem (setting)
using PROBLEM = Problem<1, GV,GV::Grid::ctype>;
using PROBLEM = Problem<GV,GV::Grid::ctype>;
PROBLEM problem;
//create model on a given setting
using MODEL = Model<dim,PROBLEM >;
using MODEL = Model<PROBLEM >;
MODEL model(problem);
//create numerical flux
......@@ -187,11 +187,11 @@ int main(int argc, char** argv)
GV gv=grid.leafGridView();
//create problem (setting)
using PROBLEM = Problem<2,GV,GV::Grid::ctype>;
using PROBLEM = Problem<GV,GV::Grid::ctype>;
PROBLEM problem;
//create model on a given setting
using MODEL = Model<dim,PROBLEM>;
using MODEL = Model<PROBLEM>;
MODEL model(problem);
//create numerical flux
......
tutorial07/src/tutorial07-linearacoustics.ini
View file @ 0d03315b
[grid]
type=structured
manager=yasp
dim=1
L = 2.5
N = 5
dim=2
L = 2.5 1
N = 5 2
refinement = 4
[fem]
......
tutorial07/src/tutorial07-shallowwater.cc
View file @ 0d03315b
......@@ -130,11 +130,11 @@ int main(int argc, char** argv)
GV gv = grid.leafGridView();
//create problem (setting)
using PROBLEM = Problem<1, GV,GV::Grid::ctype>;
using PROBLEM = Problem<GV,GV::Grid::ctype>;
PROBLEM problem;
//create model on a given setting
using MODEL = Model<dim,PROBLEM >;
using MODEL = Model<PROBLEM>;
MODEL model(problem);
//create numerical flux
......@@ -186,11 +186,11 @@ int main(int argc, char** argv)
GV gv=grid.leafGridView();
//create problem (setting)
using PROBLEM = Problem<2, GV,GV::Grid::ctype>;
using PROBLEM = Problem<GV,GV::Grid::ctype>;
PROBLEM problem;
//create model on a given setting
using MODEL = Model<dim,PROBLEM>;
using MODEL = Model<PROBLEM>;
MODEL model(problem);
//create numerical flux
......