From: bangerth Date: Fri, 23 May 2008 04:28:02 +0000 (+0000) Subject: Move computing refinement indicators into the EulerEquation class. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=177e03eab32d92570ff39261c3be19d6326f0da0;p=dealii-svn.git Move computing refinement indicators into the EulerEquation class. git-svn-id: https://svn.dealii.org/trunk@16172 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 3fcaa6f51e..4215306156 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -103,21 +103,23 @@ using namespace dealii; // @sect3{Euler equation specifics} // Here we define the flux function for this - // particular system of conservation laws, - // the Euler equations for gas dynamics. We - // group all this into a structure that - // defines everything that has to do with the - // flux. All members of this structures are - // static, i.e. the structure has no actual - // state specified by instance member - // variables. The better way to do this, - // rather than a structure with all static - // members would be to use a namespace -- but - // namespaces can't be templatized and we - // want some of the member variables of the - // structure to depend on the space + // particular system of conservation laws, as + // well as pretty much everything else that's + // specific to the Euler equations for gas + // dynamics, for reasons discussed in the + // introduction. We group all this into a + // structure that defines everything that has + // to do with the flux. All members of this + // structures are static, i.e. the structure + // has no actual state specified by instance + // member variables. The better way to do + // this, rather than a structure with all + // static members would be to use a namespace + // -- but namespaces can't be templatized and + // we want some of the member variables of + // the structure to depend on the space // dimension, which we in our usual way - // introduce using a template parameter: + // introduce using a template parameter. template struct EulerEquations { @@ -356,11 +358,11 @@ struct EulerEquations // introduction: template static - void numerical_normal_flux(const Point &normal, - const InputVector &Wplus, - const InputVector &Wminus, - const double alpha, - Sacado::Fad::DFad (&normal_flux)[n_components]) + void numerical_normal_flux (const Point &normal, + const InputVector &Wplus, + const InputVector &Wminus, + const double alpha, + Sacado::Fad::DFad (&normal_flux)[n_components]) { Sacado::Fad::DFad iflux[n_components][dim]; Sacado::Fad::DFad oflux[n_components][dim]; @@ -590,6 +592,64 @@ struct EulerEquations } + // @sect4{EulerEquations::compute_refinement_indicators} + + // In this class, we also want to specify + // how to refine the mesh. The class + // ConservationLaw that will + // use all the information we provide + // here in the EulerEquation + // class is pretty agnostic about the + // particular conservation law it solves: + // as doesn't even really care how many + // components a solution vector + // has. Consequently, it can't know what + // a reasonable refinement indicator + // would be. On the other hand, here we + // do, or at least we can come up with a + // reasonable choice: we simply look at + // the gradient of the density, and + // compute + // $\eta_K=\log\left(1+|\nabla\rho(x_K)|\right)$, + // where $x_K$ is the center of cell $K$. + // + // There are certainly a number of + // equally reasonable refinement + // indicators, but this one does, and it + // is easy to compute: + static + void + compute_refinement_indicators (const DoFHandler &dof_handler, + const Mapping &mapping, + const Vector &solution, + Vector &refinement_indicators) + { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + + const QMidpoint quadrature_formula; + const UpdateFlags update_flags = update_gradients; + FEValues fe_v (mapping, dof_handler.get_fe(), + quadrature_formula, update_flags); + + std::vector > > + dU (1, std::vector >(n_components)); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + { + fe_v.reinit(cell); + fe_v.get_function_grads (solution, dU); + + refinement_indicators(cell_no) + = std::log(1+ + std::sqrt(dU[0][density_component] * + dU[0][density_component])); + } + } + // @sect4{EulerEquations::Postprocessor} @@ -2771,55 +2831,22 @@ ConservationLaw::solve (Vector &newton_update) // @sect4{ConservationLaw::compute_refinement_indicators} - // Loop and assign a value for refinement. We - // simply use the density squared, which selects - // shocks with some success. + // This function is real simple: We don't + // pretend that we know here what a good + // refinement indicator would be. Rather, we + // assume that the EulerEquation + // class would know about this, and so we + // simply defer to the respective function + // we've implemented there: template void ConservationLaw:: compute_refinement_indicators (Vector &refinement_indicators) const -{ - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - UpdateFlags update_flags = update_values - | update_gradients - | update_q_points - | update_JxW_values; - - QGauss quadrature_formula(1); - unsigned int n_q_points = quadrature_formula.n_quadrature_points; - - - FEValues fe_v ( - mapping, fe, quadrature_formula, update_flags); - - std::vector > U(n_q_points, - Vector(EulerEquations::n_components)); - std::vector > > dU(n_q_points, - std::vector >(EulerEquations::n_components)); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) - { - fe_v.reinit(cell); - - fe_v.get_function_values(predictor, U); - fe_v.get_function_grads(predictor, dU); - - refinement_indicators(cell_no) = 0; - for (unsigned int q = 0; q < n_q_points; q++) { - double ng = 0; - for (unsigned int d = 0; d < dim; d++) - ng += dU[q][EulerEquations::density_component][d] * - dU[q][EulerEquations::density_component][d]; - - refinement_indicators(cell_no) += std::log(1+std::sqrt(ng)); - - } - refinement_indicators(cell_no) /= n_q_points; - } +{ + EulerEquations::compute_refinement_indicators (dof_handler, + mapping, + predictor, + refinement_indicators); }