From: bangerth Date: Mon, 12 May 2008 18:54:02 +0000 (+0000) Subject: Rename a few things. Probably won't compile, but will fix in a minute on a system... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=be15e577c9283f059e7821515cfadb61443a6740;p=dealii-svn.git Rename a few things. Probably won't compile, but will fix in a minute on a system with Trilinos. git-svn-id: https://svn.dealii.org/trunk@16077 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 2528bb9b41..9726a40b87 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -112,9 +112,10 @@ using namespace dealii; // structure to depend on the space // dimension, which we in our usual way // introduce using a template parameter: -namespace EulerEquations +template +struct EulerEquations { - // First a few inline functions that + // First a few variables that // describe the various components of our // solution vector in a generic way. This // includes the number of components in the @@ -127,134 +128,128 @@ namespace EulerEquations // vector of the first momentum component, // the density component, and the energy // density component. Note that all these - // numbers depend on the space dimension; + // %numbers depend on the space dimension; // defining them in a generic way (rather // than by implicit convention) makes our // code more flexible and makes it easier // to later extend it, for example by // adding more components to the equations. - template - inline - unsigned int n_components () - { - return dim + 2; - } - - template - inline - unsigned int first_momentum_component () - { - return 0; - } - - template - inline - unsigned int density_component () - { - return dim; - } - - template - inline - unsigned int energy_component () - { - return dim+1; - } - - - // Next, we define the gas constant. This - // value is representative of a gas that - // consists of molecules composed of two - // atoms, such as air which consists up to - // small traces almost entirely of $N_2$ + static const unsigned int n_components = dim + 2; + static const unsigned int first_momentum_component = 0; + static const unsigned int density_component = dim; + static const unsigned int energy_component = dim+1; + + // Next, we define the gas + // constant. We will set it to 1.4 + // in its definition immediately + // following the declaration of + // this class (unlike integer + // variables, like the ones above, + // static const floating point + // member variables cannot be + // initialized within the class + // declaration in C++). This value + // of 1.4 is representative of a + // gas that consists of molecules + // composed of two atoms, such as + // air which consists up to small + // traces almost entirely of $N_2$ // and $O_2$. - const double gas_gamma = 1.4; -} + static const double gas_gamma; + + // We define the flux function $F(W)$ as one large + // matrix. Each row of this matrix + // represents a scalar conservation law for + // the component in that row. We templatize + // the numerical type of the flux function so + // that we may use the automatic + // differentiation type here. The flux + // functions are defined in terms of the + // conserved variables $\rho w_0, \dots, \rho + // w_{d-1}, \rho, E$, so they do not look + // exactly like the Euler equations one is + // used to seeing. We evaluate the flux at a + // single quadrature point. + template + static + void flux_matrix(number (&flux)[n_components][dim], + const std::vector &W) + { + + // Pressure is a dependent variable: $p = + // (\gas_gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. + number rho_normVsqr; + for (unsigned int d=0; d -void Flux(std::vector > &flux, - const Point &/*point*/, - const std::vector &W) -{ - // Pressure is a dependent variable: $p = - // (\gas_gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. - number rho_normVsqr; - for (unsigned int d0 = 0; d0 < dim; d0++) - rho_normVsqr += W[d0]*W[d0]; - // Since W are $\rho v$, we get a $\rho^2$ in the - // numerator, so dividing a $\rho$ out gives the desired $ \rho |v|^2$. - rho_normVsqr /= W[density_component()]; - - number pressure = (gas_gamma-1.0)*(W[energy_component()] - number(0.5)*(rho_normVsqr)); - - // We compute the momentum terms. We divide by the - // density here to get $v_i \rho v_j$ - for (unsigned int d = 0; d < dim; d++) { - for (unsigned int d1 = 0; d1 < dim; d1++) { - flux[d][d1] = W[d]*W[d1]/W[density_component()]; + // On the boundaries of the domain and across hanging nodes we use + // a numerical flux function to enforce boundary conditions. This routine + // is the basic Lax-Friedrich's flux with a stabilization parameter + // $\alpha$. + template + void LFNumFlux(std::vector > > &nflux, + const std::vector > &points, + const std::vector > &normals, + const std::vector > &Wplus, + const std::vector > &Wminus, + double alpha) + { + const unsigned int n_q_points = points.size(); + + // We evaluate the flux at each of the quadrature points. + for (unsigned int q = 0; q < n_q_points; q++) + { + Sacado::Fad::DFad iflux[n_components][dim]; + Sacado::Fad::DFad oflux[n_components][dim]; + + flux_matrix(iflux, Wplus[q]); + flux_matrix(oflux, Wminus[q]); + + for (unsigned int di=0; di()][d] = W[d]; - // And, lastly, conservation of energy. - flux[energy_component()][d] = W[d]/W[density_component()]* - (W[energy_component()] + pressure); // energy - } -} +}; + - // On the boundaries of the domain and across `hanging nodes` we use - // a numerical flux function to enforce boundary conditions. This routine - // is the basic Lax-Friedrich's flux with a stabilization parameter - // $\alpha$. -template -void LFNumFlux( - std::vector > > &nflux, - const std::vector > &points, - const std::vector > &normals, - const std::vector > &Wplus, - const std::vector > &Wminus, - double alpha) -{ - const unsigned int n_q_points = points.size(); +template +const double EulerEquations::gas_gamma = 1.4; - // We evaluate the flux at each of the quadrature points. - for (unsigned int q = 0; q < n_q_points; q++) { - std::vector > > iflux(n_components(), - std::vector >(dim, 0)); - std::vector > > oflux(n_components(), - std::vector >(dim, 0)); - Flux(iflux, points[q], Wplus[q]); - Flux(oflux, points[q], Wminus[q]); - for (unsigned int di = 0; di < n_components(); di++) { - nflux[q][di] = 0; - for (unsigned int d = 0; d < dim; d++) { - nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d); - } - nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]); - } - } -} // @sect3{Initial and side condition parsing} // For the initial condition we use the expression parser function @@ -670,13 +665,12 @@ void ConsLaw::assemble_cell_term( // this could be done in a better way, since this // could be a rather large object, but for now it // seems to work just fine. - std::vector > > > flux(n_q_points, - std::vector > >(n_components(), - std::vector >(dim, 0))); + typedef Sacado::Fad::DFad FluxMatrix[EulerEquations::n_components][dim]; + std::vector flux(n_q_points); - for (unsigned int q=0; q < n_q_points; ++q) { - Flux(flux[q], fe_v.get_quadrature_points()[q], Wcn[q]); - } + for (unsigned int q=0; q < n_q_points; ++q) + flux_matrix(flux[q], Wcn[q]); + // We now have all of the function values/grads/fluxes, // so perform the assembly. We have an outer loop