Skip to content
Snippets Groups Projects
Commit 58d23ec0 authored by Piotr Minakowski's avatar Piotr Minakowski
Browse files

code cleanup

parent 6cbb897c
No related branches found
No related tags found
1 merge request!38Piotr's Tutorial 7 on hyperbolic problems
Hide whitespace changes
Inline Side-by-side
tutorial07/src/linearacoustics.hh
+ 23
− 22
View file @ 58d23ec0
......@@ -2,7 +2,7 @@
Linear Acoustics model class
*/
/** \brief provide matrix which contains rowwise the eigenvectors of linear acoustics problem
\tparam dim the space dimension
\tparam dim the space dimension
\param c speed of sound
\param n unit outer normal vector
\param RT matrix to be filled
......@@ -19,50 +19,51 @@ public:
static constexpr int dim = 2;
static constexpr int m = 3;
using RangeField = typename PROBLEM::RangeField;
using RangeField = typename PROBLEM::RangeField;
Model (PROBLEM& p)
: problem(p)
: problem(p)
{
}
template<typename T2, typename T3>
static void eigenvectors (const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,m,m>& RT)
{
int c = 1;
template<typename T2, typename T3>
static void eigenvectors (const Dune::FieldVector<T2,dim>& n, Dune::FieldMatrix<T3,m,m>& RT)
{
int c = 1;
RT[0][0] = 0; RT[0][1] = 1; RT[0][2] = -1;
RT[1][0] =-n[1]; RT[1][1] = c*n[0]; RT[1][2] = c*n[0];
RT[2][0] = n[0]; RT[2][1] = c*n[1]; RT[2][2] = c*n[1];
RT[0][0] = 0; RT[1][0] = -n[1]; RT[2][0] = n[0];
RT[0][1] = 1; RT[1][1] = c*n[0]; RT[2][1] = c*n[1];
RT[0][2] = -1; RT[1][2] = c*n[0]; RT[2][2] = c*n[1];
}
}
//one can also provide eigenvectors inverse
template<typename RF>
static void coefficients (Dune::FieldMatrix<RF,m,m>& A)
static void coefficients (Dune::FieldMatrix<RF,m,m>& A)
{
int c2 = 1;
int c2 = 1;
A[0][0] = 0.0; A[0][1] = 1.0; A[0][2] = 1.0;
A[1][0] = c2; A[1][1] = 0.0; A[1][2] = 0.0;
A[2][0] = c2; A[2][1] = 0.0; A[2][2] = 0.0;
}
}
template<typename RF>
static void diagonal (Dune::FieldMatrix<RF,m,m>& D)
static void diagonal (Dune::FieldMatrix<RF,m,m>& D)
{
int c = 1;
int c = 1;
D[0][0] = 0.0; 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] = -c;
}
}
//Flux function
//template<typename RF,typename EG>
template<typename RF>
static void flux (Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F)
{
static void flux (Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F)
{
//fetch parameters
/*
// Reference to cell
......@@ -74,14 +75,14 @@ public:
auto localcenter = ref_el.position(0,0);
auto c = param.c(cell,localcenter);
*/
int c = 1;
int c = 1;
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];
}
}
PROBLEM& problem;
PROBLEM& problem;
};
/*
......
tutorial07/src/maxwell.hh
+ 64
− 64
View file @ 58d23ec0
......@@ -2,7 +2,7 @@
Maxwell model class
*/
/** \brief provide matrix which contains rowwise the eigenvectors of linear acoustics problem
\tparam dim the space dimension
\tparam dim the space dimension
\param c speed of sound
\param n unit outer normal vector
\param RT matrix to be filled
......@@ -28,62 +28,62 @@ public:
}
template<typename T1, typename T2>
static void eigenvectors (const Dune::FieldVector<T1,dim>& n, Dune::FieldMatrix<T2,m,m>& R)
static void eigenvectors (const Dune::FieldVector<T1,dim>& n, Dune::FieldMatrix<T2,m,m>& R)
{
int eps = 1.0;
int mu = 1.0;
int eps = 1.0;
int mu = 1.0;
auto a=n[0], b=n[1], c=n[2];
Dune::FieldVector<T2,dim> alpha, beta;
if (std::abs(c)<0.5)
{
alpha[0]=a*c; alpha[1]=b*c; alpha[2]=c*c-1;
beta[0]=-b; beta[1]=a; beta[2]=0;
}
else
{
alpha[0]=a*b; alpha[1]=b*b-1; alpha[2]=b*c;
beta[0]=c; beta[1]=0.0; beta[2]=-a;
}
// \lambda_0,1 = s
R[0][0] = alpha[0]; R[0][1] = -beta[0];
R[1][0] = alpha[1]; R[1][1] = -beta[1];
R[2][0] = alpha[2]; R[2][1] = -beta[2];
R[3][0] = beta[0]; R[3][1] = alpha[0];
R[4][0] = beta[1]; R[4][1] = alpha[1];
R[5][0] = beta[2]; R[5][1] = alpha[2];
// \lambda_2,3 = -s
R[0][2] = beta[0]; R[0][3] = alpha[0];
R[1][2] = beta[1]; R[1][3] = alpha[1];
R[2][2] = beta[2]; R[2][3] = alpha[2];
R[3][2] = alpha[0]; R[3][3] = -beta[0];
R[4][2] = alpha[1]; R[4][3] = -beta[1];
R[5][2] = alpha[2]; R[5][3] = -beta[2];
// \lambda_4,5 = 0
R[0][4] = a; R[0][5] = 0;
R[1][4] = b; R[1][5] = 0;
R[2][4] = c; R[2][5] = 0;
R[3][4] = 0; R[3][5] = a;
R[4][4] = 0; R[4][5] = b;
R[5][4] = 0; R[5][5] = c;
// apply scaling
T1 weps=sqrt(eps);
T1 wmu=sqrt(mu);
for (std::size_t i=0; i<3; i++)
for (std::size_t j=0; j<6; j++)
R[i][j] *= weps;
for (std::size_t i=3; i<6; i++)
for (std::size_t j=0; j<6; j++)
R[i][j] *= wmu;
return;
};
Dune::FieldVector<T2,dim> alpha, beta;
if (std::abs(c)<0.5)
{
alpha[0]=a*c; alpha[1]=b*c; alpha[2]=c*c-1;
beta[0]=-b; beta[1]=a; beta[2]=0;
}
else
{
alpha[0]=a*b; alpha[1]=b*b-1; alpha[2]=b*c;
beta[0]=c; beta[1]=0.0; beta[2]=-a;
}
// \lambda_0,1 = s
R[0][0] = alpha[0]; R[0][1] = -beta[0];
R[1][0] = alpha[1]; R[1][1] = -beta[1];
R[2][0] = alpha[2]; R[2][1] = -beta[2];
R[3][0] = beta[0]; R[3][1] = alpha[0];
R[4][0] = beta[1]; R[4][1] = alpha[1];
R[5][0] = beta[2]; R[5][1] = alpha[2];
// \lambda_2,3 = -s
R[0][2] = beta[0]; R[0][3] = alpha[0];
R[1][2] = beta[1]; R[1][3] = alpha[1];
R[2][2] = beta[2]; R[2][3] = alpha[2];
R[3][2] = alpha[0]; R[3][3] = -beta[0];
R[4][2] = alpha[1]; R[4][3] = -beta[1];
R[5][2] = alpha[2]; R[5][3] = -beta[2];
// \lambda_4,5 = 0
R[0][4] = a; R[0][5] = 0;
R[1][4] = b; R[1][5] = 0;
R[2][4] = c; R[2][5] = 0;
R[3][4] = 0; R[3][5] = a;
R[4][4] = 0; R[4][5] = b;
R[5][4] = 0; R[5][5] = c;
// apply scaling
T1 weps=sqrt(eps);
T1 wmu=sqrt(mu);
for (std::size_t i=0; i<3; i++)
for (std::size_t j=0; j<6; j++)
R[i][j] *= weps;
for (std::size_t i=3; i<6; i++)
for (std::size_t j=0; j<6; j++)
R[i][j] *= wmu;
return;
};
//one can also provide eigenvectors inverse
template<typename RF>
......@@ -95,10 +95,10 @@ public:
A[0][0] = 0.0; A[0][1] = 0.0; A[0][2] = 0.0; A[0][3] = 0.0; A[0][4] = 1./mu; A[0][5] =-1./mu;
A[1][0] = 0.0; A[1][1] = 0.0; A[1][2] = 0.0; A[1][3] =-1./mu; A[1][4] = 0.0; A[1][5] = 1./mu;
A[2][0] = 0.0; A[2][1] = 0.0; A[2][2] = 0.0; A[2][3] = 1./mu; A[2][4] =-1./mu; A[2][5] = 1.0;
A[3][0] = 0.0; A[3][1] =-1./mu; A[3][2] = 1./mu; A[3][3] = 0.0; A[3][4] = 1.0; A[3][5] = 1.0;
A[4][0] = 1./mu; A[4][1] = 0.0; A[4][2] =-1./mu; A[4][3] = 0.0; A[4][4] = 1.0; A[4][5] = 1.0;
A[5][0] =-1./mu; A[5][1] = 1./mu; A[5][2] = 0.0; A[5][3] = 0.0; A[5][4] = 1.0; A[5][5] = 1.0;
}
A[3][0] = 0.0; A[3][1] =-1./mu; A[3][2] = 1./mu; A[3][3] = 0.0; A[3][4] = 1.0; A[3][5] = 1.0;
A[4][0] = 1./mu; A[4][1] = 0.0; A[4][2] =-1./mu; A[4][3] = 0.0; A[4][4] = 1.0; A[4][5] = 1.0;
A[5][0] =-1./mu; A[5][1] = 1./mu; A[5][2] = 0.0; A[5][3] = 0.0; A[5][4] = 1.0; A[5][5] = 1.0;
}
template<typename RF>
static void diagonal (Dune::FieldMatrix<RF,m,m>& D)
......@@ -108,15 +108,15 @@ public:
D[0][0] = c; D[0][1] = 0.0; D[0][2] = 0.0; D[0][3] = 0.0; D[0][4] = 0.0; D[0][5] = 0.0;
D[1][0] = 0.0; D[1][1] = c; D[1][2] = 0.0; D[1][3] = 0.0; D[1][4] = 0.0; D[1][5] = 0.0;
D[2][0] = 0.0; D[2][1] = 0.0; D[2][2] = -c; D[2][3] = 0.0; D[2][4] = 0.0; D[2][5] = 0.0;
D[3][0] = 0.0; D[3][1] = 0.0; D[3][2] = 0.0; D[3][3] = -c; D[3][4] = 0.0; D[3][5] = 0.0;
D[4][0] = 0.0; D[4][1] = 0.0; D[4][2] = 0.0; D[4][3] = 0.0; D[4][4] = 0.0; D[4][5] = 0.0;
D[5][0] = 0.0; D[5][1] = 0.0; D[5][2] = 0.0; D[5][3] = 0.0; D[5][4] = 0.0; D[5][5] = 0.0;
}
D[3][0] = 0.0; D[3][1] = 0.0; D[3][2] = 0.0; D[3][3] = -c; D[3][4] = 0.0; D[3][5] = 0.0;
D[4][0] = 0.0; D[4][1] = 0.0; D[4][2] = 0.0; D[4][3] = 0.0; D[4][4] = 0.0; D[4][5] = 0.0;
D[5][0] = 0.0; D[5][1] = 0.0; D[5][2] = 0.0; D[5][3] = 0.0; D[5][4] = 0.0; D[5][5] = 0.0;
}
//Flux function
template<typename RF>
static void flux (Dune::FieldVector<RF,m>& u, Dune::FieldMatrix<RF,m,dim>& F)
{
{
RF mu(1.0);
RF ep(1.0);
......@@ -126,9 +126,9 @@ public:
F[3][0] = 0.0 ; F[3][1] = 1/ep*u[2]; F[3][2] =-1/ep*u[1];
F[4][0] =-1/ep*u[2]; F[4][1] = 0.0; F[4][2] = 1/ep*u[0];
F[5][0] = 1/ep*u[1]; F[5][1] =-1/ep*u[0]; F[5][2] = 0.0;
}
}
PROBLEM& problem;
PROBLEM& problem;
};
......
tutorial07/src/maxwellproblem.hh
+ 7
− 7
View file @ 58d23ec0
......@@ -70,17 +70,17 @@ public:
return rhs;
}
//! initial value
//! initial value
template<typename E, typename X>
Range u0 (const E& e, const X& p) const
{
X xglobal = e.geometry().global(p);
Range u(0.0);
auto x=xglobal[0], y=xglobal[1], z=xglobal[2];
/*
/*
if (x>0.4 && x<0.6 && y>0.4 && y<0.6 && z>0.4 && z<0.6)
{
auto p1 = sin(pi*(x-0.4)/0.2)*sin(pi*(x-0.4)/0.2);
......@@ -93,9 +93,9 @@ public:
u[1] = -0.5*p1prime*q2*r3prime;
u[2] = -0.5*p1prime*q2prime*r3;
}
*/
auto alpha = 1.0;
*/
auto alpha = 1.0;
auto s1 = std::sin(alpha*x);
auto c1 = std::cos(alpha*x);
......
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