From: heister Date: Tue, 20 Nov 2012 23:15:38 +0000 (+0000) Subject: merge indentation manually X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=769ad4450a3673c35c1d69d0baac874ae5bf71ad;p=dealii-svn.git merge indentation manually git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27630 0785d39b-7218-0410-832d-ea1e28bc413d --- 769ad4450a3673c35c1d69d0baac874ae5bf71ad diff --cc deal.II/examples/step-13/step-13.cc index e8182e5bd0,8e527273d0..63e62cc64b --- a/deal.II/examples/step-13/step-13.cc +++ b/deal.II/examples/step-13/step-13.cc @@@ -168,15 -168,15 +168,15 @@@ namespace Step1 template class EvaluationBase { - public: - virtual ~EvaluationBase (); + public: + virtual ~EvaluationBase (); - void set_refinement_cycle (const unsigned int refinement_cycle); + void set_refinement_cycle (const unsigned int refinement_cycle); - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const = 0; - protected: - unsigned int refinement_cycle; + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const = 0; ++ const Vector &solution) const = 0; + protected: + unsigned int refinement_cycle; }; @@@ -252,20 -252,20 +252,20 @@@ template class PointValueEvaluation : public EvaluationBase { - public: - PointValueEvaluation (const Point &evaluation_point, - TableHandler &results_table); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - private: - const Point evaluation_point; - TableHandler &results_table; + public: + PointValueEvaluation (const Point &evaluation_point, + TableHandler &results_table); + + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; ++ const Vector &solution) const; + + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; + TableHandler &results_table; }; @@@ -292,36 -292,36 +292,36 @@@ void PointValueEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { - // First allocate a variable that - // will hold the point - // value. Initialize it with a - // value that is clearly bogus, - // so that if we fail to set it - // to a reasonable value, we will - // note at once. This may not be - // necessary in a function as - // small as this one, since we - // can easily see all possible - // paths of execution here, but - // it proved to be helpful for - // more complex cases, and so we - // employ this strategy here as - // well. + // First allocate a variable that + // will hold the point + // value. Initialize it with a + // value that is clearly bogus, + // so that if we fail to set it + // to a reasonable value, we will + // note at once. This may not be + // necessary in a function as + // small as this one, since we + // can easily see all possible + // paths of execution here, but + // it proved to be helpful for + // more complex cases, and so we + // employ this strategy here as + // well. double point_value = 1e20; - // Then loop over all cells and - // all their vertices, and check - // whether a vertex matches the - // evaluation point. If this is - // the case, then extract the - // point value, set a flag that - // we have found the point of - // interest, and exit the loop. + // Then loop over all cells and + // all their vertices, and check + // whether a vertex matches the + // evaluation point. If this is + // the case, then extract the + // point value, set a flag that + // we have found the point of + // interest, and exit the loop. typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = dof_handler.begin_active(), + endc = dof_handler.end(); bool evaluation_point_found = false; for (; (cell!=endc) && !evaluation_point_found; ++cell) for (unsigned int vertex=0; @@@ -610,15 -610,15 +610,15 @@@ template class SolutionOutput : public EvaluationBase { - public: - SolutionOutput (const std::string &output_name_base, - const typename DataOut::OutputFormat output_format); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - private: - const std::string output_name_base; - const typename DataOut::OutputFormat output_format; + public: + SolutionOutput (const std::string &output_name_base, + const typename DataOut::OutputFormat output_format); + + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; ++ const Vector &solution) const; + private: + const std::string output_name_base; + const typename DataOut::OutputFormat output_format; }; diff --cc deal.II/examples/step-14/step-14.cc index 4a6185ee26,9c58a26932..c1d4613872 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@@ -72,15 -72,15 +72,15 @@@ namespace Step1 template class EvaluationBase { - public: - virtual ~EvaluationBase (); + public: + virtual ~EvaluationBase (); - void set_refinement_cycle (const unsigned int refinement_cycle); + void set_refinement_cycle (const unsigned int refinement_cycle); - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const = 0; - protected: - unsigned int refinement_cycle; + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const = 0; ++ const Vector &solution) const = 0; + protected: + unsigned int refinement_cycle; }; @@@ -102,18 -102,18 +102,18 @@@ template class PointValueEvaluation : public EvaluationBase { - public: - PointValueEvaluation (const Point &evaluation_point); + public: + PointValueEvaluation (const Point &evaluation_point); - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; ++ const Vector &solution) const; - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - private: - const Point evaluation_point; + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; }; @@@ -186,18 -186,18 +186,18 @@@ template class PointXDerivativeEvaluation : public EvaluationBase { - public: - PointXDerivativeEvaluation (const Point &evaluation_point); + public: + PointXDerivativeEvaluation (const Point &evaluation_point); - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; ++ const Vector &solution) const; - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - private: - const Point evaluation_point; + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + private: + const Point evaluation_point; }; @@@ -216,19 -216,19 +216,19 @@@ void PointXDerivativeEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { - // This time initialize the - // return value with something - // useful, since we will have to - // add up a number of - // contributions and take the - // mean value afterwards... + // This time initialize the + // return value with something + // useful, since we will have to + // add up a number of + // contributions and take the + // mean value afterwards... double point_derivative = 0; - // ...then have some objects of - // which the meaning wil become - // clear below... + // ...then have some objects of + // which the meaning wil become + // clear below... QTrapez vertex_quadrature; FEValues fe_values (dof_handler.get_fe(), vertex_quadrature, @@@ -399,13 -399,13 +399,13 @@@ template class GridOutput : public EvaluationBase { - public: - GridOutput (const std::string &output_name_base); + public: + GridOutput (const std::string &output_name_base); - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - private: - const std::string output_name_base; + virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; ++ const Vector &solution) const; + private: + const std::string output_name_base; }; diff --cc deal.II/examples/step-18/step-18.cc index 4b41046f84,e1528caebb..af124c3f5b --- a/deal.II/examples/step-18/step-18.cc +++ b/deal.II/examples/step-18/step-18.cc @@@ -845,11 -853,11 +853,11 @@@ namespace Step1 template inline void - BodyForce::vector_value (const Point & /*p*/, + BodyForce::vector_value (const Point &/*p*/, - Vector &values) const + Vector &values) const { Assert (values.size() == dim, - ExcDimensionMismatch (values.size(), dim)); + ExcDimensionMismatch (values.size(), dim)); const double g = 9.81; const double rho = 7700; @@@ -966,11 -974,11 +974,11 @@@ template void IncrementalBoundaryValues:: - vector_value (const Point & /*p*/, + vector_value (const Point &/*p*/, - Vector &values) const + Vector &values) const { Assert (values.size() == dim, - ExcDimensionMismatch (values.size(), dim)); + ExcDimensionMismatch (values.size(), dim)); values = 0; values(2) = -present_timestep * velocity; diff --cc deal.II/examples/step-20/step-20.cc index 8a3e1093f3,01303b4d7e..fd816b8be1 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@@ -191,15 -191,15 +191,15 @@@ namespace Step2 }; - // And then we also have to define - // these respective functions, of - // course. Given our discussion in - // the introduction of how the - // solution should look like, the - // following computations should be - // straightforward: + // And then we also have to define + // these respective functions, of + // course. Given our discussion in + // the introduction of how the + // solution should look like, the + // following computations should be + // straightforward: template - double RightHandSide::value (const Point & /*p*/, + double RightHandSide::value (const Point &/*p*/, const unsigned int /*component*/) const { return 0; diff --cc deal.II/examples/step-28/step-28.cc index 64b0276c58,a3c477b375..0363c2dad6 --- a/deal.II/examples/step-28/step-28.cc +++ b/deal.II/examples/step-28/step-28.cc @@@ -1592,94 -1592,94 +1592,94 @@@ namespace Step2 template class NeutronDiffusionProblem { + public: + class Parameters + { public: - class Parameters - { - public: - Parameters (); - - static void declare_parameters (ParameterHandler &prm); - void get_parameters (ParameterHandler &prm); - - unsigned int n_groups; - unsigned int n_refinement_cycles; - - unsigned int fe_degree; - - double convergence_tolerance; - }; - - - - NeutronDiffusionProblem (const Parameters ¶meters); - ~NeutronDiffusionProblem (); - - void run (); - - private: - // @sect5{Private member functions} - - // There are not that many member - // functions in this class since - // most of the functionality has - // been moved into the - // EnergyGroup class - // and is simply called from the - // run() member - // function of this class. The - // ones that remain have - // self-explanatory names: - void initialize_problem(); - - void refine_grid (); - - double get_total_fission_source () const; - - - // @sect5{Private member variables} - - // Next, we have a few member - // variables. In particular, - // these are (i) a reference to - // the parameter object (owned by - // the main function of this - // program, and passed to the - // constructor of this class), - // (ii) an object describing the - // material parameters for the - // number of energy groups - // requested in the input file, - // and (iii) the finite element - // to be used by all energy - // groups: - const Parameters ¶meters; - const MaterialData material_data; - FE_Q fe; - - // Furthermore, we have (iv) the - // value of the computed - // eigenvalue at the present - // iteration. This is, in fact, - // the only part of the solution - // that is shared between all - // energy groups -- all other - // parts of the solution, such as - // neutron fluxes are particular - // to one or the other energy - // group, and are therefore - // stored in objects that - // describe a single energy - // group: - double k_eff; - - // Finally, (v), we have an array - // of pointers to the energy - // group objects. The length of - // this array is, of course, - // equal to the number of energy - // groups specified in the - // parameter file. - std::vector*> energy_groups; + Parameters (); + + static void declare_parameters (ParameterHandler &prm); + void get_parameters (ParameterHandler &prm); + + unsigned int n_groups; + unsigned int n_refinement_cycles; + + unsigned int fe_degree; + + double convergence_tolerance; + }; + + + + NeutronDiffusionProblem (const Parameters ¶meters); + ~NeutronDiffusionProblem (); + + void run (); + + private: + // @sect5{Private member functions} + + // There are not that many member + // functions in this class since + // most of the functionality has + // been moved into the + // EnergyGroup class + // and is simply called from the + // run() member + // function of this class. The + // ones that remain have + // self-explanatory names: + void initialize_problem(); + + void refine_grid (); + + double get_total_fission_source () const; + + + // @sect5{Private member variables} + + // Next, we have a few member + // variables. In particular, + // these are (i) a reference to + // the parameter object (owned by + // the main function of this + // program, and passed to the + // constructor of this class), + // (ii) an object describing the + // material parameters for the + // number of energy groups + // requested in the input file, + // and (iii) the finite element + // to be used by all energy + // groups: - const Parameters ¶meters; ++ const Parameters ¶meters; + const MaterialData material_data; + FE_Q fe; + + // Furthermore, we have (iv) the + // value of the computed + // eigenvalue at the present + // iteration. This is, in fact, + // the only part of the solution + // that is shared between all + // energy groups -- all other + // parts of the solution, such as + // neutron fluxes are particular + // to one or the other energy + // group, and are therefore + // stored in objects that + // describe a single energy + // group: + double k_eff; + + // Finally, (v), we have an array + // of pointers to the energy + // group objects. The length of + // this array is, of course, + // equal to the number of energy + // groups specified in the + // parameter file. + std::vector*> energy_groups; }; diff --cc deal.II/examples/step-29/step-29.cc index 183e15034a,4e0c83615e..4e4bff6dd9 --- a/deal.II/examples/step-29/step-29.cc +++ b/deal.II/examples/step-29/step-29.cc @@@ -572,20 -572,20 +572,20 @@@ namespace Step2 - // The constructor takes the - // ParameterHandler object and stores - // it in a reference. It also - // initializes the DoF-Handler and - // the finite element system, which - // consists of two copies of the - // scalar Q1 field, one for $v$ and - // one for $w$: + // The constructor takes the + // ParameterHandler object and stores + // it in a reference. It also + // initializes the DoF-Handler and + // the finite element system, which + // consists of two copies of the + // scalar Q1 field, one for $v$ and + // one for $w$: template - UltrasoundProblem::UltrasoundProblem (ParameterHandler& param) - : - prm(param), - dof_handler(triangulation), - fe(FE_Q(1), 2) - UltrasoundProblem::UltrasoundProblem (ParameterHandler ¶m) ++ UltrasoundProblem::UltrasoundProblem (ParameterHandler ¶m) + : + prm(param), + dof_handler(triangulation), + fe(FE_Q(1), 2) {} diff --cc deal.II/examples/step-31/step-31.cc index 1db324114d,39b6350bcf..033ff91776 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@@ -192,11 -192,11 +192,11 @@@ namespace Step3 template double - TemperatureRightHandSide::value (const Point &p, + TemperatureRightHandSide::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { Assert (component == 0, - ExcMessage ("Invalid operation for a scalar function.")); + ExcMessage ("Invalid operation for a scalar function.")); Assert ((dim==2) || (dim==3), ExcNotImplemented()); @@@ -488,15 -489,15 +489,15 @@@ template BlockSchurPreconditioner:: - BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, + BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, - const InverseMatrix &Mpinv, - const PreconditionerA &Apreconditioner) - : - stokes_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (stokes_matrix->block(1,1).m()) + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner) + : + stokes_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (stokes_matrix->block(1,1).m()) {} @@@ -581,80 -582,80 +582,80 @@@ template class BoussinesqFlowProblem { - public: - BoussinesqFlowProblem (); - void run (); + public: + BoussinesqFlowProblem (); + void run (); + + private: + void setup_dofs (); + void assemble_stokes_preconditioner (); + void build_stokes_preconditioner (); + void assemble_stokes_system (); + void assemble_temperature_system (const double maximal_velocity); + void assemble_temperature_matrix (); + double get_maximal_velocity () const; + std::pair get_extrapolated_temperature_range () const; + void solve (); + void output_results () const; + void refine_mesh (const unsigned int max_grid_level); - private: - void setup_dofs (); - void assemble_stokes_preconditioner (); - void build_stokes_preconditioner (); - void assemble_stokes_system (); - void assemble_temperature_system (const double maximal_velocity); - void assemble_temperature_matrix (); - double get_maximal_velocity () const; - std::pair get_extrapolated_temperature_range () const; - void solve (); - void output_results () const; - void refine_mesh (const unsigned int max_grid_level); - - double - compute_viscosity(const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector &old_temperature_laplacians, - const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector &gamma_values, - const double global_u_infty, - const double global_T_variation, - const double cell_diameter) const; - - - Triangulation triangulation; - double global_Omega_diameter; - - const unsigned int stokes_degree; - FESystem stokes_fe; - DoFHandler stokes_dof_handler; - ConstraintMatrix stokes_constraints; - - std::vector stokes_block_sizes; - TrilinosWrappers::BlockSparseMatrix stokes_matrix; - TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; - - TrilinosWrappers::BlockVector stokes_solution; - TrilinosWrappers::BlockVector old_stokes_solution; - TrilinosWrappers::BlockVector stokes_rhs; - - - const unsigned int temperature_degree; - FE_Q temperature_fe; - DoFHandler temperature_dof_handler; - ConstraintMatrix temperature_constraints; - - TrilinosWrappers::SparseMatrix temperature_mass_matrix; - TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; - TrilinosWrappers::SparseMatrix temperature_matrix; - - TrilinosWrappers::Vector temperature_solution; - TrilinosWrappers::Vector old_temperature_solution; - TrilinosWrappers::Vector old_old_temperature_solution; - TrilinosWrappers::Vector temperature_rhs; - - - double time_step; - double old_time_step; - unsigned int timestep_number; - - std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; - - bool rebuild_stokes_matrix; - bool rebuild_temperature_matrices; - bool rebuild_stokes_preconditioner; + double + compute_viscosity(const std::vector &old_temperature, + const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, ++ const std::vector > &old_temperature_grads, ++ const std::vector > &old_old_temperature_grads, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, ++ const std::vector > &old_velocity_values, ++ const std::vector > &old_old_velocity_values, + const std::vector &gamma_values, + const double global_u_infty, + const double global_T_variation, + const double cell_diameter) const; + + + Triangulation triangulation; + double global_Omega_diameter; + + const unsigned int stokes_degree; + FESystem stokes_fe; + DoFHandler stokes_dof_handler; + ConstraintMatrix stokes_constraints; + + std::vector stokes_block_sizes; + TrilinosWrappers::BlockSparseMatrix stokes_matrix; + TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + + TrilinosWrappers::BlockVector stokes_solution; + TrilinosWrappers::BlockVector old_stokes_solution; + TrilinosWrappers::BlockVector stokes_rhs; + + + const unsigned int temperature_degree; + FE_Q temperature_fe; + DoFHandler temperature_dof_handler; + ConstraintMatrix temperature_constraints; + + TrilinosWrappers::SparseMatrix temperature_mass_matrix; + TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; + TrilinosWrappers::SparseMatrix temperature_matrix; + + TrilinosWrappers::Vector temperature_solution; + TrilinosWrappers::Vector old_temperature_solution; + TrilinosWrappers::Vector old_old_temperature_solution; + TrilinosWrappers::Vector temperature_rhs; + + + double time_step; + double old_time_step; + unsigned int timestep_number; + + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + + bool rebuild_stokes_matrix; + bool rebuild_temperature_matrices; + bool rebuild_stokes_preconditioner; }; @@@ -961,17 -962,17 +962,17 @@@ double BoussinesqFlowProblem:: compute_viscosity (const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector &old_temperature_laplacians, - const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector &gamma_values, - const double global_u_infty, - const double global_T_variation, - const double cell_diameter) const + const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, ++ const std::vector > &old_temperature_grads, ++ const std::vector > &old_old_temperature_grads, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, ++ const std::vector > &old_velocity_values, ++ const std::vector > &old_old_velocity_values, + const std::vector &gamma_values, + const double global_u_infty, + const double global_T_variation, + const double cell_diameter) const { const double beta = 0.015 * dim; const double alpha = 1; diff --cc deal.II/examples/step-32/step-32.cc index f75422b4d8,db3c7123d8..b81ffa873d --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@@ -177,8 -177,8 +177,8 @@@ namespace Step3 template double - TemperatureInitialValues::value (const Point &p, + TemperatureInitialValues::value (const Point &p, - const unsigned int) const + const unsigned int) const { const double r = p.norm(); const double h = R1-R0; @@@ -268,60 -268,60 +268,60 @@@ template class BlockSchurPreconditioner : public Subscriptor { - public: - BlockSchurPreconditioner (const TrilinosWrappers::BlockSparseMatrix &S, - const TrilinosWrappers::BlockSparseMatrix &Spre, - const PreconditionerMp &Mppreconditioner, - const PreconditionerA &Apreconditioner, - const bool do_solve_A) - : - stokes_matrix (&S), - stokes_preconditioner_matrix (&Spre), - mp_preconditioner (Mppreconditioner), - a_preconditioner (Apreconditioner), - do_solve_A (do_solve_A) - {} - - void vmult (TrilinosWrappers::MPI::BlockVector &dst, - const TrilinosWrappers::MPI::BlockVector &src) const - { - TrilinosWrappers::MPI::Vector utmp(src.block(0)); - - { - SolverControl solver_control(5000, 1e-6 * src.block(1).l2_norm()); - - SolverCG solver(solver_control); - - solver.solve(stokes_preconditioner_matrix->block(1,1), - dst.block(1), src.block(1), - mp_preconditioner); - - dst.block(1) *= -1.0; - } - - { - stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); - utmp*=-1.0; - utmp.add(src.block(0)); - } - - if (do_solve_A == true) - { - SolverControl solver_control(5000, utmp.l2_norm()*1e-2); - TrilinosWrappers::SolverCG solver(solver_control); - solver.solve(stokes_matrix->block(0,0), dst.block(0), utmp, - a_preconditioner); - } - else - a_preconditioner.vmult (dst.block(0), utmp); - } - - private: - const SmartPointer stokes_matrix; - const SmartPointer stokes_preconditioner_matrix; - const PreconditionerMp &mp_preconditioner; - const PreconditionerA &a_preconditioner; - const bool do_solve_A; + public: - BlockSchurPreconditioner (const TrilinosWrappers::BlockSparseMatrix &S, - const TrilinosWrappers::BlockSparseMatrix &Spre, ++ BlockSchurPreconditioner (const TrilinosWrappers::BlockSparseMatrix &S, ++ const TrilinosWrappers::BlockSparseMatrix &Spre, + const PreconditionerMp &Mppreconditioner, + const PreconditionerA &Apreconditioner, + const bool do_solve_A) + : + stokes_matrix (&S), + stokes_preconditioner_matrix (&Spre), + mp_preconditioner (Mppreconditioner), + a_preconditioner (Apreconditioner), + do_solve_A (do_solve_A) + {} + + void vmult (TrilinosWrappers::MPI::BlockVector &dst, + const TrilinosWrappers::MPI::BlockVector &src) const + { + TrilinosWrappers::MPI::Vector utmp(src.block(0)); + + { + SolverControl solver_control(5000, 1e-6 * src.block(1).l2_norm()); + + SolverCG solver(solver_control); + + solver.solve(stokes_preconditioner_matrix->block(1,1), + dst.block(1), src.block(1), + mp_preconditioner); + + dst.block(1) *= -1.0; + } + + { + stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); + utmp*=-1.0; + utmp.add(src.block(0)); + } + + if (do_solve_A == true) + { + SolverControl solver_control(5000, utmp.l2_norm()*1e-2); + TrilinosWrappers::SolverCG solver(solver_control); + solver.solve(stokes_matrix->block(0,0), dst.block(0), utmp, + a_preconditioner); + } + else + a_preconditioner.vmult (dst.block(0), utmp); + } + + private: + const SmartPointer stokes_matrix; + const SmartPointer stokes_preconditioner_matrix; + const PreconditionerMp &mp_preconditioner; - const PreconditionerA &a_preconditioner; ++ const PreconditionerA &a_preconditioner; + const bool do_solve_A; }; } @@@ -893,352 -893,352 +893,352 @@@ template class BoussinesqFlowProblem { - public: - struct Parameters; - BoussinesqFlowProblem (Parameters ¶meters); - void run (); + public: + struct Parameters; + BoussinesqFlowProblem (Parameters ¶meters); + void run (); + + private: + void setup_dofs (); + void assemble_stokes_preconditioner (); + void build_stokes_preconditioner (); + void assemble_stokes_system (); + void assemble_temperature_matrix (); + void assemble_temperature_system (const double maximal_velocity); + void project_temperature_field (); + double get_maximal_velocity () const; + double get_cfl_number () const; + double get_entropy_variation (const double average_temperature) const; + std::pair get_extrapolated_temperature_range () const; + void solve (); + void output_results (); + void refine_mesh (const unsigned int max_grid_level); - private: - void setup_dofs (); - void assemble_stokes_preconditioner (); - void build_stokes_preconditioner (); - void assemble_stokes_system (); - void assemble_temperature_matrix (); - void assemble_temperature_system (const double maximal_velocity); - void project_temperature_field (); - double get_maximal_velocity () const; - double get_cfl_number () const; - double get_entropy_variation (const double average_temperature) const; - std::pair get_extrapolated_temperature_range () const; - void solve (); - void output_results (); - void refine_mesh (const unsigned int max_grid_level); - - double - compute_viscosity(const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector &old_temperature_laplacians, - const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector > &old_strain_rates, - const std::vector > &old_old_strain_rates, - const double global_u_infty, - const double global_T_variation, - const double average_temperature, - const double global_entropy_variation, - const double cell_diameter) const; + double + compute_viscosity(const std::vector &old_temperature, + const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, ++ const std::vector > &old_temperature_grads, ++ const std::vector > &old_old_temperature_grads, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector > &old_strain_rates, - const std::vector > &old_old_strain_rates, ++ const std::vector > &old_velocity_values, ++ const std::vector > &old_old_velocity_values, ++ const std::vector > &old_strain_rates, ++ const std::vector > &old_old_strain_rates, + const double global_u_infty, + const double global_T_variation, + const double average_temperature, + const double global_entropy_variation, + const double cell_diameter) const; + + public: + + // The first significant new + // component is the definition + // of a struct for the + // parameters according to the + // discussion in the + // introduction. This structure + // is initialized by reading + // from a parameter file during + // construction of this object. + struct Parameters + { + Parameters (const std::string ¶meter_filename); - public: + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); - // The first significant new - // component is the definition - // of a struct for the - // parameters according to the - // discussion in the - // introduction. This structure - // is initialized by reading - // from a parameter file during - // construction of this object. - struct Parameters - { - Parameters (const std::string ¶meter_filename); + double end_time; - static void declare_parameters (ParameterHandler &prm); - void parse_parameters (ParameterHandler &prm); + unsigned int initial_global_refinement; + unsigned int initial_adaptive_refinement; - double end_time; + bool generate_graphical_output; + unsigned int graphical_output_interval; - unsigned int initial_global_refinement; - unsigned int initial_adaptive_refinement; + unsigned int adaptive_refinement_interval; - bool generate_graphical_output; - unsigned int graphical_output_interval; + double stabilization_alpha; + double stabilization_c_R; + double stabilization_beta; - unsigned int adaptive_refinement_interval; + unsigned int stokes_velocity_degree; + bool use_locally_conservative_discretization; - double stabilization_alpha; - double stabilization_c_R; - double stabilization_beta; + unsigned int temperature_degree; + }; - unsigned int stokes_velocity_degree; - bool use_locally_conservative_discretization; + private: + Parameters ¶meters; + + // The pcout (for + // %parallel + // std::cout) + // object is used to simplify + // writing output: each MPI + // process can use this to + // generate output as usual, + // but since each of these + // processes will (hopefully) + // produce the same output it + // will just be replicated many + // times over; with the + // ConditionalOStream class, + // only the output generated by + // one MPI process will + // actually be printed to + // screen, whereas the output + // by all the other threads + // will simply be forgotten. + ConditionalOStream pcout; + + // The following member + // variables will then again be + // similar to those in step-31 + // (and to other tutorial + // programs). As mentioned in + // the introduction, we fully + // distribute computations, so + // we will have to use the + // parallel::distributed::Triangulation + // class (see step-40) but the + // remainder of these variables + // is rather standard with two + // exceptions: + // + // - The mapping + // variable is used to denote a + // higher-order polynomial + // mapping. As mentioned in the + // introduction, we use this + // mapping when forming + // integrals through quadrature + // for all cells that are + // adjacent to either the inner + // or outer boundaries of our + // domain where the boundary is + // curved. + // + // - In a bit of naming + // confusion, you will notice + // below that some of the + // variables from namespace + // TrilinosWrappers are taken + // from namespace + // TrilinosWrappers::MPI (such + // as the right hand side + // vectors) whereas others are + // not (such as the various + // matrices). For the matrices, + // we happen to use the same + // class names for %parallel + // and sequential data + // structures, i.e., all + // matrices will actually be + // considered %parallel + // below. On the other hand, + // for vectors, only those from + // namespace + // TrilinosWrappers::MPI are + // actually distributed. In + // particular, we will + // frequently have to query + // velocities and temperatures + // at arbitrary quadrature + // points; consequently, rather + // than importing ghost + // information of a vector + // whenever we need access to + // degrees of freedom that are + // relevant locally but owned + // by another processor, we + // solve linear systems in + // %parallel but then + // immediately initialize a + // vector including ghost + // entries of the solution for + // further processing. The + // various + // *_solution + // vectors are therefore filled + // immediately after solving + // their respective linear + // system in %parallel and will + // always contain values for + // all @ref + // GlossLocallyRelevantDof + // "locally relevant degrees of freedom"; + // the fully + // distributed vectors that we + // obtain from the solution + // process and that only ever + // contain the @ref + // GlossLocallyOwnedDof + // "locally owned degrees of freedom" + // are destroyed + // immediately after the + // solution process and after + // we have copied the relevant + // values into the member + // variable vectors. + parallel::distributed::Triangulation triangulation; + double global_Omega_diameter; + + const MappingQ mapping; + + const FESystem stokes_fe; + DoFHandler stokes_dof_handler; + ConstraintMatrix stokes_constraints; + + TrilinosWrappers::BlockSparseMatrix stokes_matrix; + TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; + + TrilinosWrappers::MPI::BlockVector stokes_solution; + TrilinosWrappers::MPI::BlockVector old_stokes_solution; + TrilinosWrappers::MPI::BlockVector stokes_rhs; + + + FE_Q temperature_fe; + DoFHandler temperature_dof_handler; + ConstraintMatrix temperature_constraints; + + TrilinosWrappers::SparseMatrix temperature_mass_matrix; + TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; + TrilinosWrappers::SparseMatrix temperature_matrix; + + TrilinosWrappers::MPI::Vector temperature_solution; + TrilinosWrappers::MPI::Vector old_temperature_solution; + TrilinosWrappers::MPI::Vector old_old_temperature_solution; + TrilinosWrappers::MPI::Vector temperature_rhs; + + + double time_step; + double old_time_step; + unsigned int timestep_number; + + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + std_cxx1x::shared_ptr T_preconditioner; + + bool rebuild_stokes_matrix; + bool rebuild_stokes_preconditioner; + bool rebuild_temperature_matrices; + bool rebuild_temperature_preconditioner; + + // The next member variable, + // computing_timer + // is used to conveniently + // account for compute time + // spent in certain "sections" + // of the code that are + // repeatedly entered. For + // example, we will enter (and + // leave) sections for Stokes + // matrix assembly and would + // like to accumulate the run + // time spent in this section + // over all time steps. Every + // so many time steps as well + // as at the end of the program + // (through the destructor of + // the TimerOutput class) we + // will then produce a nice + // summary of the times spent + // in the different sections + // into which we categorize the + // run-time of this program. + TimerOutput computing_timer; + + // After these member variables + // we have a number of + // auxiliary functions that + // have been broken out of the + // ones listed + // above. Specifically, there + // are first three functions + // that we call from + // setup_dofs and + // then the ones that do the + // assembling of linear + // systems: + void setup_stokes_matrix (const std::vector &stokes_partitioning); + void setup_stokes_preconditioner (const std::vector &stokes_partitioning); + void setup_temperature_matrices (const IndexSet &temperature_partitioning); + + + // Following the @ref + // MTWorkStream + // "task-based parallelization" + // paradigm, + // we split all the assembly + // routines into two parts: a + // first part that can do all + // the calculations on a + // certain cell without taking + // care of other threads, and a + // second part (which is + // writing the local data into + // the global matrices and + // vectors) which can be + // entered by only one thread + // at a time. In order to + // implement that, we provide + // functions for each of those + // two steps for all the four + // assembly routines that we + // use in this program. The + // following eight functions do + // exactly this: + void + local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::StokesPreconditioner &scratch, + Assembly::CopyData::StokesPreconditioner &data); - unsigned int temperature_degree; - }; + void + copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data); - private: - Parameters ¶meters; - - // The pcout (for - // %parallel - // std::cout) - // object is used to simplify - // writing output: each MPI - // process can use this to - // generate output as usual, - // but since each of these - // processes will (hopefully) - // produce the same output it - // will just be replicated many - // times over; with the - // ConditionalOStream class, - // only the output generated by - // one MPI process will - // actually be printed to - // screen, whereas the output - // by all the other threads - // will simply be forgotten. - ConditionalOStream pcout; - - // The following member - // variables will then again be - // similar to those in step-31 - // (and to other tutorial - // programs). As mentioned in - // the introduction, we fully - // distribute computations, so - // we will have to use the - // parallel::distributed::Triangulation - // class (see step-40) but the - // remainder of these variables - // is rather standard with two - // exceptions: - // - // - The mapping - // variable is used to denote a - // higher-order polynomial - // mapping. As mentioned in the - // introduction, we use this - // mapping when forming - // integrals through quadrature - // for all cells that are - // adjacent to either the inner - // or outer boundaries of our - // domain where the boundary is - // curved. - // - // - In a bit of naming - // confusion, you will notice - // below that some of the - // variables from namespace - // TrilinosWrappers are taken - // from namespace - // TrilinosWrappers::MPI (such - // as the right hand side - // vectors) whereas others are - // not (such as the various - // matrices). For the matrices, - // we happen to use the same - // class names for %parallel - // and sequential data - // structures, i.e., all - // matrices will actually be - // considered %parallel - // below. On the other hand, - // for vectors, only those from - // namespace - // TrilinosWrappers::MPI are - // actually distributed. In - // particular, we will - // frequently have to query - // velocities and temperatures - // at arbitrary quadrature - // points; consequently, rather - // than importing ghost - // information of a vector - // whenever we need access to - // degrees of freedom that are - // relevant locally but owned - // by another processor, we - // solve linear systems in - // %parallel but then - // immediately initialize a - // vector including ghost - // entries of the solution for - // further processing. The - // various - // *_solution - // vectors are therefore filled - // immediately after solving - // their respective linear - // system in %parallel and will - // always contain values for - // all @ref - // GlossLocallyRelevantDof - // "locally relevant degrees of freedom"; - // the fully - // distributed vectors that we - // obtain from the solution - // process and that only ever - // contain the @ref - // GlossLocallyOwnedDof - // "locally owned degrees of freedom" - // are destroyed - // immediately after the - // solution process and after - // we have copied the relevant - // values into the member - // variable vectors. - parallel::distributed::Triangulation triangulation; - double global_Omega_diameter; - - const MappingQ mapping; - - const FESystem stokes_fe; - DoFHandler stokes_dof_handler; - ConstraintMatrix stokes_constraints; - - TrilinosWrappers::BlockSparseMatrix stokes_matrix; - TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix; - - TrilinosWrappers::MPI::BlockVector stokes_solution; - TrilinosWrappers::MPI::BlockVector old_stokes_solution; - TrilinosWrappers::MPI::BlockVector stokes_rhs; - - - FE_Q temperature_fe; - DoFHandler temperature_dof_handler; - ConstraintMatrix temperature_constraints; - - TrilinosWrappers::SparseMatrix temperature_mass_matrix; - TrilinosWrappers::SparseMatrix temperature_stiffness_matrix; - TrilinosWrappers::SparseMatrix temperature_matrix; - - TrilinosWrappers::MPI::Vector temperature_solution; - TrilinosWrappers::MPI::Vector old_temperature_solution; - TrilinosWrappers::MPI::Vector old_old_temperature_solution; - TrilinosWrappers::MPI::Vector temperature_rhs; - - - double time_step; - double old_time_step; - unsigned int timestep_number; - - std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; - std_cxx1x::shared_ptr T_preconditioner; - - bool rebuild_stokes_matrix; - bool rebuild_stokes_preconditioner; - bool rebuild_temperature_matrices; - bool rebuild_temperature_preconditioner; - - // The next member variable, - // computing_timer - // is used to conveniently - // account for compute time - // spent in certain "sections" - // of the code that are - // repeatedly entered. For - // example, we will enter (and - // leave) sections for Stokes - // matrix assembly and would - // like to accumulate the run - // time spent in this section - // over all time steps. Every - // so many time steps as well - // as at the end of the program - // (through the destructor of - // the TimerOutput class) we - // will then produce a nice - // summary of the times spent - // in the different sections - // into which we categorize the - // run-time of this program. - TimerOutput computing_timer; - - // After these member variables - // we have a number of - // auxiliary functions that - // have been broken out of the - // ones listed - // above. Specifically, there - // are first three functions - // that we call from - // setup_dofs and - // then the ones that do the - // assembling of linear - // systems: - void setup_stokes_matrix (const std::vector &stokes_partitioning); - void setup_stokes_preconditioner (const std::vector &stokes_partitioning); - void setup_temperature_matrices (const IndexSet &temperature_partitioning); - - - // Following the @ref - // MTWorkStream - // "task-based parallelization" - // paradigm, - // we split all the assembly - // routines into two parts: a - // first part that can do all - // the calculations on a - // certain cell without taking - // care of other threads, and a - // second part (which is - // writing the local data into - // the global matrices and - // vectors) which can be - // entered by only one thread - // at a time. In order to - // implement that, we provide - // functions for each of those - // two steps for all the four - // assembly routines that we - // use in this program. The - // following eight functions do - // exactly this: - void - local_assemble_stokes_preconditioner (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesPreconditioner &scratch, - Assembly::CopyData::StokesPreconditioner &data); - - void - copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner &data); - - - void - local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesSystem &scratch, - Assembly::CopyData::StokesSystem &data); - - void - copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data); - - - void - local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureMatrix &scratch, - Assembly::CopyData::TemperatureMatrix &data); - - void - copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data); - - - - void - local_assemble_temperature_rhs (const std::pair global_T_range, - const double global_max_velocity, - const double global_entropy_variation, - const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureRHS &scratch, - Assembly::CopyData::TemperatureRHS &data); - - void - copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data); - - // Finally, we forward declare - // a member class that we will - // define later on and that - // will be used to compute a - // number of quantities from - // our solution vectors that - // we'd like to put into the - // output files for - // visualization. - class Postprocessor; + + void + local_assemble_stokes_system (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::StokesSystem &scratch, ++ Assembly::Scratch::StokesSystem &scratch, + Assembly::CopyData::StokesSystem &data); + + void + copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem &data); + + + void + local_assemble_temperature_matrix (const typename DoFHandler::active_cell_iterator &cell, - Assembly::Scratch::TemperatureMatrix &scratch, ++ Assembly::Scratch::TemperatureMatrix &scratch, + Assembly::CopyData::TemperatureMatrix &data); + + void + copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix &data); + + + + void + local_assemble_temperature_rhs (const std::pair global_T_range, + const double global_max_velocity, + const double global_entropy_variation, + const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::TemperatureRHS &scratch, + Assembly::CopyData::TemperatureRHS &data); + + void + copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS &data); + + // Finally, we forward declare + // a member class that we will + // define later on and that + // will be used to compute a + // number of quantities from + // our solution vectors that + // we'd like to put into the + // output files for + // visualization. + class Postprocessor; }; @@@ -1967,20 -1968,20 +1968,20 @@@ double BoussinesqFlowProblem:: compute_viscosity (const std::vector &old_temperature, - const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, - const std::vector &old_temperature_laplacians, - const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector > &old_strain_rates, - const std::vector > &old_old_strain_rates, - const double global_u_infty, - const double global_T_variation, - const double average_temperature, - const double global_entropy_variation, - const double cell_diameter) const + const std::vector &old_old_temperature, - const std::vector > &old_temperature_grads, - const std::vector > &old_old_temperature_grads, ++ const std::vector > &old_temperature_grads, ++ const std::vector > &old_old_temperature_grads, + const std::vector &old_temperature_laplacians, + const std::vector &old_old_temperature_laplacians, - const std::vector > &old_velocity_values, - const std::vector > &old_old_velocity_values, - const std::vector > &old_strain_rates, - const std::vector > &old_old_strain_rates, ++ const std::vector > &old_velocity_values, ++ const std::vector > &old_old_velocity_values, ++ const std::vector > &old_strain_rates, ++ const std::vector > &old_old_strain_rates, + const double global_u_infty, + const double global_T_variation, + const double average_temperature, + const double global_entropy_variation, + const double cell_diameter) const { if (global_u_infty == 0) return 5e-3 * cell_diameter; diff --cc deal.II/examples/step-33/step-33.cc index dccad4abe0,e402227987..1721b8b098 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@@ -107,532 -107,533 +107,533 @@@ namespace Step3 template struct EulerEquations { - // @sect4{Component description} - - // 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 - // system (Euler's equations have one entry - // for momenta in each spatial direction, - // plus the energy and density components, - // for a total of dim+2 - // components), as well as functions that - // describe the index within the solution - // vector of the first momentum component, - // the density component, and the energy - // density component. Note that all these - // %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. - 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; - - // When generating graphical - // output way down in this - // program, we need to specify - // the names of the solution - // variables as well as how the - // various components group into - // vector and scalar fields. We - // could describe this there, but - // in order to keep things that - // have to do with the Euler - // equation localized here and - // the rest of the program as - // generic as possible, we - // provide this sort of - // information in the following - // two functions: - static - std::vector - component_names () - { - std::vector names (dim, "momentum"); - names.push_back ("density"); - names.push_back ("energy_density"); + // @sect4{Component description} + + // 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 + // system (Euler's equations have one entry + // for momenta in each spatial direction, + // plus the energy and density components, + // for a total of dim+2 + // components), as well as functions that + // describe the index within the solution + // vector of the first momentum component, + // the density component, and the energy + // density component. Note that all these + // %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. + 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; + + // When generating graphical + // output way down in this + // program, we need to specify + // the names of the solution + // variables as well as how the + // various components group into + // vector and scalar fields. We + // could describe this there, but + // in order to keep things that + // have to do with the Euler + // equation localized here and + // the rest of the program as + // generic as possible, we + // provide this sort of + // information in the following + // two functions: + static + std::vector + component_names () + { + std::vector names (dim, "momentum"); + names.push_back ("density"); + names.push_back ("energy_density"); - return names; - } + return names; + } - static + static + std::vector + component_interpretation () + { std::vector - component_interpretation () - { - std::vector - data_component_interpretation - (dim, DataComponentInterpretation::component_is_part_of_vector); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - data_component_interpretation - .push_back (DataComponentInterpretation::component_is_scalar); - - return data_component_interpretation; - } + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + return data_component_interpretation; + } - // @sect4{Transformations between variables} - - // 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$. - static const double gas_gamma; - - - // In the following, we will need to - // compute the kinetic energy and the - // pressure from a vector of conserved - // variables. This we can do based on the - // energy density and the kinetic energy - // $\frac 12 \rho |\mathbf v|^2 = - // \frac{|\rho \mathbf v|^2}{2\rho}$ - // (note that the independent variables - // contain the momentum components $\rho - // v_i$, not the velocities $v_i$). - // - // There is one slight problem: We will - // need to call the following functions - // with input arguments of type - // std::vector@ and - // Vector@. The - // problem is that the former has an - // access operator - // operator[] whereas the - // latter, for historical reasons, has - // operator(). We wouldn't - // be able to write the function in a - // generic way if we were to use one or - // the other of these. Fortunately, we - // can use the following trick: instead - // of writing v[i] or - // v(i), we can use - // *(v.begin() + i), i.e. we - // generate an iterator that points to - // the ith element, and then - // dereference it. This works for both - // kinds of vectors -- not the prettiest - // solution, but one that works. - template - static - number - compute_kinetic_energy (const InputVector &W) - { - number kinetic_energy = 0; - for (unsigned int d=0; dstd::vector@ and + // Vector@. The + // problem is that the former has an + // access operator + // operator[] whereas the + // latter, for historical reasons, has + // operator(). We wouldn't + // be able to write the function in a + // generic way if we were to use one or + // the other of these. Fortunately, we + // can use the following trick: instead + // of writing v[i] or + // v(i), we can use + // *(v.begin() + i), i.e. we + // generate an iterator that points to + // the ith element, and then + // dereference it. This works for both + // kinds of vectors -- not the prettiest + // solution, but one that works. + template + static + number + compute_kinetic_energy (const InputVector &W) + { + number kinetic_energy = 0; + for (unsigned int d=0; d - static - number - compute_pressure (const InputVector &W) - { - return ((gas_gamma-1.0) * - (*(W.begin() + energy_component) - - compute_kinetic_energy(W))); - } + template + static + number + compute_pressure (const InputVector &W) + { + return ((gas_gamma-1.0) * + (*(W.begin() + energy_component) - + compute_kinetic_energy(W))); + } - // @sect4{EulerEquations::compute_flux_matrix} - - // 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. The - // exact form of this matrix is - // given in the - // introduction. Note that we - // know the size of the matrix: - // it has as many rows as the - // system has components, and - // dim columns; - // rather than using a FullMatrix - // object for such a matrix - // (which has a variable number - // of rows and columns and must - // therefore allocate memory on - // the heap each time such a - // matrix is created), we use a - // rectangular array of numbers - // right away. - // - // We templatize the numerical type of - // the flux function so that we may use - // the automatic differentiation type - // here. Similarly, we will call the - // function with different input vector - // data types, so we templatize on it as - // well: - template - static - void compute_flux_matrix (const InputVector &W, - number (&flux)[n_components][dim]) + // @sect4{EulerEquations::compute_flux_matrix} + + // 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. The + // exact form of this matrix is + // given in the + // introduction. Note that we + // know the size of the matrix: + // it has as many rows as the + // system has components, and + // dim columns; + // rather than using a FullMatrix + // object for such a matrix + // (which has a variable number + // of rows and columns and must + // therefore allocate memory on + // the heap each time such a + // matrix is created), we use a + // rectangular array of numbers + // right away. + // + // We templatize the numerical type of + // the flux function so that we may use + // the automatic differentiation type + // here. Similarly, we will call the + // function with different input vector + // data types, so we templatize on it as + // well: + template + static + void compute_flux_matrix (const InputVector &W, + number (&flux)[n_components][dim]) + { + // First compute the pressure that + // appears in the flux matrix, and + // then compute the first + // dim columns of the + // matrix that correspond to the + // momentum terms: + const number pressure = compute_pressure (W); + + for (unsigned int d=0; ddim columns of the - // matrix that correspond to the - // momentum terms: - const number pressure = compute_pressure (W); - - for (unsigned int d=0; d - static - 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]; - - compute_flux_matrix (Wplus, iflux); - compute_flux_matrix (Wminus, oflux); - for (unsigned int di=0; di + static + 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]; - normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]); - } - } + compute_flux_matrix (Wplus, iflux); + compute_flux_matrix (Wminus, oflux); - // @sect4{EulerEquations::compute_forcing_vector} - - // In the same way as describing the flux - // function $\mathbf F(\mathbf w)$, we - // also need to have a way to describe - // the right hand side forcing term. As - // mentioned in the introduction, we - // consider only gravity here, which - // leads to the specific form $\mathbf - // G(\mathbf w) = \left( - // g_1\rho, g_2\rho, g_3\rho, 0, - // \rho \mathbf g \cdot \mathbf v - // \right)^T$, shown here for - // the 3d case. More specifically, we - // will consider only $\mathbf - // g=(0,0,-1)^T$ in 3d, or $\mathbf - // g=(0,-1)^T$ in 2d. This naturally - // leads to the following function: - template - static - void compute_forcing_vector (const InputVector &W, - number (&forcing)[n_components]) + for (unsigned int di=0; di + static + void compute_forcing_vector (const InputVector &W, + number (&forcing)[n_components]) + { + const double gravity = -1.0; - // @sect4{Dealing with boundary conditions} + for (unsigned int c=0; cWminus will of course be - // modified, so it shouldn't be a - // const argument. Yet it is - // in the implementation below, and needs - // to be in order to allow the code to - // compile. The reason is that we call - // this function at a place where - // Wminus is of type - // Table@<2,Sacado::Fad::DFad@ - // @>, this being 2d table with - // indices representing the quadrature - // point and the vector component, - // respectively. We call this function - // with Wminus[q] as last - // argument; subscripting a 2d table - // yields a temporary accessor object - // representing a 1d vector, just what we - // want here. The problem is that a - // temporary accessor object can't be - // bound to a non-const reference - // argument of a function, as we would - // like here, according to the C++ 1998 - // and 2003 standards (something that - // will be fixed with the next standard - // in the form of rvalue references). We - // get away with making the output - // argument here a constant because it is - // the accessor object that's - // constant, not the table it points to: - // that one can still be written to. The - // hack is unpleasant nevertheless - // because it restricts the kind of data - // types that may be used as template - // argument to this function: a regular - // vector isn't going to do because that - // one can not be written to when marked - // const. With no good - // solution around at the moment, we'll - // go with the pragmatic, even if not - // pretty, solution shown here: - template - static - void - compute_Wminus (const BoundaryKind (&boundary_kind)[n_components], - const Point &normal_vector, - const DataVector &Wplus, - const Vector &boundary_values, - const DataVector &Wminus) - { - for (unsigned int c = 0; c < n_components; c++) - switch (boundary_kind[c]) - { - case inflow_boundary: - { - Wminus[c] = boundary_values(c); - break; - } + // Another thing we have to deal with is + // boundary conditions. To this end, let + // us first define the kinds of boundary + // conditions we currently know how to + // deal with: + enum BoundaryKind + { + inflow_boundary, + outflow_boundary, + no_penetration_boundary, + pressure_boundary + }; - case outflow_boundary: - { - Wminus[c] = Wplus[c]; - break; - } - // Prescribed pressure boundary - // conditions are a bit more - // complicated by the fact that - // even though the pressure is - // prescribed, we really are - // setting the energy component - // here, which will depend on - // velocity and pressure. So - // even though this seems like - // a Dirichlet type boundary - // condition, we get - // sensitivities of energy to - // velocity and density (unless - // these are also prescribed): - case pressure_boundary: - { - const typename DataVector::value_type - density = (boundary_kind[density_component] == - inflow_boundary - ? - boundary_values(density_component) - : - Wplus[density_component]); + // The next part is to actually decide + // what to do at each kind of + // boundary. To this end, remember from + // the introduction that boundary + // conditions are specified by choosing a + // value $\mathbf w^-$ on the outside of + // a boundary given an inhomogeneity + // $\mathbf j$ and possibly the + // solution's value $\mathbf w^+$ on the + // inside. Both are then passed to the + // numerical flux $\mathbf + // H(\mathbf{w}^+, \mathbf{w}^-, + // \mathbf{n})$ to define boundary + // contributions to the bilinear form. + // + // Boundary conditions can in some cases + // be specified for each component of the + // solution vector independently. For + // example, if component $c$ is marked + // for inflow, then $w^-_c = j_c$. If it + // is an outflow, then $w^-_c = + // w^+_c$. These two simple cases are + // handled first in the function below. + // + // There is a little snag that makes this + // function unpleasant from a C++ + // language viewpoint: The output vector + // Wminus will of course be + // modified, so it shouldn't be a + // const argument. Yet it is + // in the implementation below, and needs + // to be in order to allow the code to + // compile. The reason is that we call + // this function at a place where + // Wminus is of type + // Table@<2,Sacado::Fad::DFad@ + // @>, this being 2d table with + // indices representing the quadrature + // point and the vector component, + // respectively. We call this function + // with Wminus[q] as last + // argument; subscripting a 2d table + // yields a temporary accessor object + // representing a 1d vector, just what we + // want here. The problem is that a + // temporary accessor object can't be + // bound to a non-const reference + // argument of a function, as we would + // like here, according to the C++ 1998 + // and 2003 standards (something that + // will be fixed with the next standard + // in the form of rvalue references). We + // get away with making the output + // argument here a constant because it is + // the accessor object that's + // constant, not the table it points to: + // that one can still be written to. The + // hack is unpleasant nevertheless + // because it restricts the kind of data + // types that may be used as template + // argument to this function: a regular + // vector isn't going to do because that + // one can not be written to when marked + // const. With no good + // solution around at the moment, we'll + // go with the pragmatic, even if not + // pretty, solution shown here: + template + static + void + compute_Wminus (const BoundaryKind (&boundary_kind)[n_components], + const Point &normal_vector, + const DataVector &Wplus, + const Vector &boundary_values, + const DataVector &Wminus) + { + for (unsigned int c = 0; c < n_components; c++) + switch (boundary_kind[c]) + { + case inflow_boundary: + { + Wminus[c] = boundary_values(c); + break; + } - typename DataVector::value_type kinetic_energy = 0; - for (unsigned int d=0; d vdotn = 0; + for (unsigned int d = 0; d < dim; d++) + { + vdotn += Wplus[d]*normal_vector[d]; + } - case no_penetration_boundary: - { - // We prescribe the - // velocity (we are dealing with a - // particular component here so - // that the average of the - // velocities is orthogonal to the - // surface normal. This creates - // sensitivies of across the - // velocity components. - Sacado::Fad::DFad vdotn = 0; - for (unsigned int d = 0; d < dim; d++) { - vdotn += Wplus[d]*normal_vector[d]; - } - - Wminus[c] = Wplus[c] - 2.0*vdotn*normal_vector[c]; - break; - } + Wminus[c] = Wplus[c] - 2.0*vdotn*normal_vector[c]; + break; + } - default: - Assert (false, ExcNotImplemented()); - } - } + default: + Assert (false, ExcNotImplemented()); + } + } - // @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); + // @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, ++ 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); + 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)); + 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); + 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])); - } + refinement_indicators(cell_no) + = std::log(1+ + std::sqrt(dU[0][density_component] * + dU[0][density_component])); } + } @@@ -764,26 -765,26 +765,26 @@@ EulerEquations::Postprocessor:: compute_derived_quantities_vector (const std::vector > &uh, const std::vector > > &duh, - const std::vector > > & /*dduh*/, - const std::vector > & /*normals*/, - const std::vector > & /*evaluation_points*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &/*evaluation_points*/, std::vector > &computed_quantities) const { - // At the beginning of the function, let us - // make sure that all variables have the - // correct sizes, so that we can access - // individual vector elements without - // having to wonder whether we might read - // or write invalid elements; we also check - // that the duh vector only - // contains data if we really need it (the - // system knows about this because we say - // so in the - // get_needed_update_flags() - // function below). For the inner vectors, - // we check that at least the first element - // of the outer vector has the correct - // inner size: + // At the beginning of the function, let us + // make sure that all variables have the + // correct sizes, so that we can access + // individual vector elements without + // having to wonder whether we might read + // or write invalid elements; we also check + // that the duh vector only + // contains data if we really need it (the + // system knows about this because we say + // so in the + // get_needed_update_flags() + // function below). For the inner vectors, + // we check that at least the first element + // of the outer vector has the correct + // inner size: const unsigned int n_quadrature_points = uh.size(); if (do_schlieren_plot == true) diff --cc deal.II/examples/step-35/step-35.cc index cd57ab99b7,cb37555fe7..109da8c0dd --- a/deal.II/examples/step-35/step-35.cc +++ b/deal.II/examples/step-35/step-35.cc @@@ -419,267 -419,267 +419,267 @@@ namespace Step3 template class NavierStokesProjection { - public: - NavierStokesProjection (const RunTimeParameters::Data_Storage &data); - - void run (const bool verbose = false, - const unsigned int n_plots = 10); - protected: - RunTimeParameters::MethodFormulation type; - - const unsigned int deg; - const double dt; - const double t_0, T, Re; - - EquationData::Velocity vel_exact; - std::map boundary_values; - std::vector boundary_indicators; - - Triangulation triangulation; - - FE_Q fe_velocity; - FE_Q fe_pressure; - - DoFHandler dof_handler_velocity; - DoFHandler dof_handler_pressure; - - QGauss quadrature_pressure; - QGauss quadrature_velocity; - - SparsityPattern sparsity_pattern_velocity; - SparsityPattern sparsity_pattern_pressure; - SparsityPattern sparsity_pattern_pres_vel; - - SparseMatrix vel_Laplace_plus_Mass; - SparseMatrix vel_it_matrix[dim]; - SparseMatrix vel_Mass; - SparseMatrix vel_Laplace; - SparseMatrix vel_Advection; - SparseMatrix pres_Laplace; - SparseMatrix pres_Mass; - SparseMatrix pres_Diff[dim]; - SparseMatrix pres_iterative; - - Vector pres_n; - Vector pres_n_minus_1; - Vector phi_n; - Vector phi_n_minus_1; - Vector u_n[dim]; - Vector u_n_minus_1[dim]; - Vector u_star[dim]; - Vector force[dim]; - Vector v_tmp; - Vector pres_tmp; - Vector rot_u; - - SparseILU prec_velocity[dim]; - SparseILU prec_pres_Laplace; - SparseDirectUMFPACK prec_mass; - SparseDirectUMFPACK prec_vel_mass; - - DeclException2 (ExcInvalidTimeStep, - double, double, - << " The time step " << arg1 << " is out of range." - << std::endl - << " The permitted range is (0," << arg2 << "]"); - - void create_triangulation_and_dofs (const unsigned int n_refines); - - void initialize(); - - void interpolate_velocity (); - - void diffusion_step (const bool reinit_prec); - - void projection_step (const bool reinit_prec); - - void update_pressure (const bool reinit_prec); - - private: - unsigned int vel_max_its; - unsigned int vel_Krylov_size; - unsigned int vel_off_diagonals; - unsigned int vel_update_prec; - double vel_eps; - double vel_diag_strength; - - void initialize_velocity_matrices(); - - void initialize_pressure_matrices(); - - // The next few structures and functions - // are for doing various things in - // parallel. They follow the scheme laid - // out in @ref threads, using the - // WorkStream class. As explained there, - // this requires us to declare two - // structures for each of the assemblers, - // a per-task data and a scratch data - // structure. These are then handed over - // to functions that assemble local - // contributions and that copy these - // local contributions to the global - // objects. - // - // One of the things that are specific to - // this program is that we don't just - // have a single DoFHandler object that - // represents both the velocities and the - // pressure, but we use individual - // DoFHandler objects for these two kinds - // of variables. We pay for this - // optimization when we want to assemble - // terms that involve both variables, - // such as the divergence of the velocity - // and the gradient of the pressure, - // times the respective test - // functions. When doing so, we can't - // just anymore use a single FEValues - // object, but rather we need two, and - // they need to be initialized with cell - // iterators that point to the same cell - // in the triangulation but different - // DoFHandlers. - // - // To do this in practice, we declare a - // "synchronous" iterator -- an object - // that internally consists of several - // (in our case two) iterators, and each - // time the synchronous iteration is - // moved up one step, each of the - // iterators stored internally is moved - // up one step as well, thereby always - // staying in sync. As it so happens, - // there is a deal.II class that - // facilitates this sort of thing. - typedef std_cxx1x::tuple< typename DoFHandler::active_cell_iterator, - typename DoFHandler::active_cell_iterator - > IteratorTuple; - - typedef SynchronousIterators IteratorPair; - - void initialize_gradient_operator(); - - struct InitGradPerTaskData - { - unsigned int d; - unsigned int vel_dpc; - unsigned int pres_dpc; - FullMatrix local_grad; - std::vector vel_local_dof_indices; - std::vector pres_local_dof_indices; - - InitGradPerTaskData (const unsigned int dd, - const unsigned int vdpc, - const unsigned int pdpc) - : - d(dd), - vel_dpc (vdpc), - pres_dpc (pdpc), - local_grad (vdpc, pdpc), - vel_local_dof_indices (vdpc), - pres_local_dof_indices (pdpc) - {} - }; + public: + NavierStokesProjection (const RunTimeParameters::Data_Storage &data); + + void run (const bool verbose = false, + const unsigned int n_plots = 10); + protected: + RunTimeParameters::MethodFormulation type; + + const unsigned int deg; + const double dt; + const double t_0, T, Re; + + EquationData::Velocity vel_exact; + std::map boundary_values; + std::vector boundary_indicators; + + Triangulation triangulation; + + FE_Q fe_velocity; + FE_Q fe_pressure; + + DoFHandler dof_handler_velocity; + DoFHandler dof_handler_pressure; + + QGauss quadrature_pressure; + QGauss quadrature_velocity; + + SparsityPattern sparsity_pattern_velocity; + SparsityPattern sparsity_pattern_pressure; + SparsityPattern sparsity_pattern_pres_vel; + + SparseMatrix vel_Laplace_plus_Mass; + SparseMatrix vel_it_matrix[dim]; + SparseMatrix vel_Mass; + SparseMatrix vel_Laplace; + SparseMatrix vel_Advection; + SparseMatrix pres_Laplace; + SparseMatrix pres_Mass; + SparseMatrix pres_Diff[dim]; + SparseMatrix pres_iterative; + + Vector pres_n; + Vector pres_n_minus_1; + Vector phi_n; + Vector phi_n_minus_1; + Vector u_n[dim]; + Vector u_n_minus_1[dim]; + Vector u_star[dim]; + Vector force[dim]; + Vector v_tmp; + Vector pres_tmp; + Vector rot_u; + + SparseILU prec_velocity[dim]; + SparseILU prec_pres_Laplace; + SparseDirectUMFPACK prec_mass; + SparseDirectUMFPACK prec_vel_mass; + + DeclException2 (ExcInvalidTimeStep, + double, double, + << " The time step " << arg1 << " is out of range." + << std::endl + << " The permitted range is (0," << arg2 << "]"); + + void create_triangulation_and_dofs (const unsigned int n_refines); + + void initialize(); + + void interpolate_velocity (); + + void diffusion_step (const bool reinit_prec); + + void projection_step (const bool reinit_prec); + + void update_pressure (const bool reinit_prec); + + private: + unsigned int vel_max_its; + unsigned int vel_Krylov_size; + unsigned int vel_off_diagonals; + unsigned int vel_update_prec; + double vel_eps; + double vel_diag_strength; + + void initialize_velocity_matrices(); + + void initialize_pressure_matrices(); + + // The next few structures and functions + // are for doing various things in + // parallel. They follow the scheme laid + // out in @ref threads, using the + // WorkStream class. As explained there, + // this requires us to declare two + // structures for each of the assemblers, + // a per-task data and a scratch data + // structure. These are then handed over + // to functions that assemble local + // contributions and that copy these + // local contributions to the global + // objects. + // + // One of the things that are specific to + // this program is that we don't just + // have a single DoFHandler object that + // represents both the velocities and the + // pressure, but we use individual + // DoFHandler objects for these two kinds + // of variables. We pay for this + // optimization when we want to assemble + // terms that involve both variables, + // such as the divergence of the velocity + // and the gradient of the pressure, + // times the respective test + // functions. When doing so, we can't + // just anymore use a single FEValues + // object, but rather we need two, and + // they need to be initialized with cell + // iterators that point to the same cell + // in the triangulation but different + // DoFHandlers. + // + // To do this in practice, we declare a + // "synchronous" iterator -- an object + // that internally consists of several + // (in our case two) iterators, and each + // time the synchronous iteration is + // moved up one step, each of the + // iterators stored internally is moved + // up one step as well, thereby always + // staying in sync. As it so happens, + // there is a deal.II class that + // facilitates this sort of thing. + typedef std_cxx1x::tuple< typename DoFHandler::active_cell_iterator, + typename DoFHandler::active_cell_iterator + > IteratorTuple; + + typedef SynchronousIterators IteratorPair; + + void initialize_gradient_operator(); + + struct InitGradPerTaskData + { + unsigned int d; + unsigned int vel_dpc; + unsigned int pres_dpc; + FullMatrix local_grad; + std::vector vel_local_dof_indices; + std::vector pres_local_dof_indices; + + InitGradPerTaskData (const unsigned int dd, + const unsigned int vdpc, + const unsigned int pdpc) + : + d(dd), + vel_dpc (vdpc), + pres_dpc (pdpc), + local_grad (vdpc, pdpc), + vel_local_dof_indices (vdpc), + pres_local_dof_indices (pdpc) + {} + }; - struct InitGradScratchData - { - unsigned int nqp; - FEValues fe_val_vel; - FEValues fe_val_pres; - InitGradScratchData (const FE_Q &fe_v, - const FE_Q &fe_p, - const QGauss &quad, - const UpdateFlags flags_v, - const UpdateFlags flags_p) - : - nqp (quad.size()), - fe_val_vel (fe_v, quad, flags_v), - fe_val_pres (fe_p, quad, flags_p) - {} - InitGradScratchData (const InitGradScratchData &data) - : - nqp (data.nqp), - fe_val_vel (data.fe_val_vel.get_fe(), - data.fe_val_vel.get_quadrature(), - data.fe_val_vel.get_update_flags()), - fe_val_pres (data.fe_val_pres.get_fe(), - data.fe_val_pres.get_quadrature(), - data.fe_val_pres.get_update_flags()) - {} - }; + struct InitGradScratchData + { + unsigned int nqp; + FEValues fe_val_vel; + FEValues fe_val_pres; + InitGradScratchData (const FE_Q &fe_v, + const FE_Q &fe_p, + const QGauss &quad, + const UpdateFlags flags_v, + const UpdateFlags flags_p) + : + nqp (quad.size()), + fe_val_vel (fe_v, quad, flags_v), + fe_val_pres (fe_p, quad, flags_p) + {} + InitGradScratchData (const InitGradScratchData &data) + : + nqp (data.nqp), + fe_val_vel (data.fe_val_vel.get_fe(), + data.fe_val_vel.get_quadrature(), + data.fe_val_vel.get_update_flags()), + fe_val_pres (data.fe_val_pres.get_fe(), + data.fe_val_pres.get_quadrature(), + data.fe_val_pres.get_update_flags()) + {} + }; - void assemble_one_cell_of_gradient (const IteratorPair &SI, - InitGradScratchData &scratch, - InitGradPerTaskData &data); - void assemble_one_cell_of_gradient (const IteratorPair &SI, ++ void assemble_one_cell_of_gradient (const IteratorPair &SI, + InitGradScratchData &scratch, + InitGradPerTaskData &data); - void copy_gradient_local_to_global (const InitGradPerTaskData &data); + void copy_gradient_local_to_global (const InitGradPerTaskData &data); - // The same general layout also applies - // to the following classes and functions - // implementing the assembly of the - // advection term: - void assemble_advection_term(); + // The same general layout also applies + // to the following classes and functions + // implementing the assembly of the + // advection term: + void assemble_advection_term(); - struct AdvectionPerTaskData - { - FullMatrix local_advection; - std::vector local_dof_indices; - AdvectionPerTaskData (const unsigned int dpc) - : - local_advection (dpc, dpc), - local_dof_indices (dpc) - {} - }; + struct AdvectionPerTaskData + { + FullMatrix local_advection; + std::vector local_dof_indices; + AdvectionPerTaskData (const unsigned int dpc) + : + local_advection (dpc, dpc), + local_dof_indices (dpc) + {} + }; - struct AdvectionScratchData - { - unsigned int nqp; - unsigned int dpc; - std::vector< Point > u_star_local; - std::vector< Tensor<1,dim> > grad_u_star; - std::vector u_star_tmp; - FEValues fe_val; - AdvectionScratchData (const FE_Q &fe, - const QGauss &quad, - const UpdateFlags flags) - : - nqp (quad.size()), - dpc (fe.dofs_per_cell), - u_star_local (nqp), - grad_u_star (nqp), - u_star_tmp (nqp), - fe_val (fe, quad, flags) - {} - - AdvectionScratchData (const AdvectionScratchData &data) - : - nqp (data.nqp), - dpc (data.dpc), - u_star_local (nqp), - grad_u_star (nqp), - u_star_tmp (nqp), - fe_val (data.fe_val.get_fe(), - data.fe_val.get_quadrature(), - data.fe_val.get_update_flags()) - {} - }; + struct AdvectionScratchData + { + unsigned int nqp; + unsigned int dpc; + std::vector< Point > u_star_local; + std::vector< Tensor<1,dim> > grad_u_star; + std::vector u_star_tmp; + FEValues fe_val; + AdvectionScratchData (const FE_Q &fe, + const QGauss &quad, + const UpdateFlags flags) + : + nqp (quad.size()), + dpc (fe.dofs_per_cell), + u_star_local (nqp), + grad_u_star (nqp), + u_star_tmp (nqp), + fe_val (fe, quad, flags) + {} + + AdvectionScratchData (const AdvectionScratchData &data) + : + nqp (data.nqp), + dpc (data.dpc), + u_star_local (nqp), + grad_u_star (nqp), + u_star_tmp (nqp), + fe_val (data.fe_val.get_fe(), + data.fe_val.get_quadrature(), + data.fe_val.get_update_flags()) + {} + }; - void assemble_one_cell_of_advection (const typename DoFHandler::active_cell_iterator &cell, - AdvectionScratchData &scratch, - AdvectionPerTaskData &data); + void assemble_one_cell_of_advection (const typename DoFHandler::active_cell_iterator &cell, + AdvectionScratchData &scratch, + AdvectionPerTaskData &data); - void copy_advection_local_to_global (const AdvectionPerTaskData &data); + void copy_advection_local_to_global (const AdvectionPerTaskData &data); - // The final few functions implement the - // diffusion solve as well as - // postprocessing the output, including - // computing the curl of the velocity: - void diffusion_component_solve (const unsigned int d); + // The final few functions implement the + // diffusion solve as well as + // postprocessing the output, including + // computing the curl of the velocity: + void diffusion_component_solve (const unsigned int d); - void output_results (const unsigned int step); + void output_results (const unsigned int step); - void assemble_vorticity (const bool reinit_prec); + void assemble_vorticity (const bool reinit_prec); }; diff --cc deal.II/examples/step-37/step-37.cc index a664f42dcc,807b710c8e..080c0d9277 --- a/deal.II/examples/step-37/step-37.cc +++ b/deal.II/examples/step-37/step-37.cc @@@ -389,46 -389,46 +389,46 @@@ namespace Step3 template class LaplaceOperator : public Subscriptor { - public: - LaplaceOperator (); + public: + LaplaceOperator (); - void clear(); + void clear(); - void reinit (const MGDoFHandler &dof_handler, - const ConstraintMatrix &constraints, - const unsigned int level = numbers::invalid_unsigned_int); + void reinit (const MGDoFHandler &dof_handler, - const ConstraintMatrix &constraints, ++ const ConstraintMatrix &constraints, + const unsigned int level = numbers::invalid_unsigned_int); - unsigned int m () const; - unsigned int n () const; + unsigned int m () const; + unsigned int n () const; - void vmult (Vector &dst, - const Vector &src) const; - void Tvmult (Vector &dst, - const Vector &src) const; - void vmult_add (Vector &dst, - const Vector &src) const; - void Tvmult_add (Vector &dst, - const Vector &src) const; + void vmult (Vector &dst, + const Vector &src) const; + void Tvmult (Vector &dst, + const Vector &src) const; + void vmult_add (Vector &dst, + const Vector &src) const; + void Tvmult_add (Vector &dst, + const Vector &src) const; - number el (const unsigned int row, - const unsigned int col) const; - void set_diagonal (const Vector &diagonal); + number el (const unsigned int row, + const unsigned int col) const; + void set_diagonal (const Vector &diagonal); - std::size_t memory_consumption () const; + std::size_t memory_consumption () const; - private: - void local_apply (const MatrixFree &data, - Vector &dst, - const Vector &src, - const std::pair &cell_range) const; + private: + void local_apply (const MatrixFree &data, + Vector &dst, + const Vector &src, + const std::pair &cell_range) const; - void evaluate_coefficient(const Coefficient &function); + void evaluate_coefficient(const Coefficient &function); - MatrixFree data; - AlignedVector > coefficient; + MatrixFree data; + AlignedVector > coefficient; - Vector diagonal_values; - bool diagonal_is_available; + Vector diagonal_values; + bool diagonal_is_available; }; @@@ -572,8 -572,8 +572,8 @@@ template void LaplaceOperator::reinit (const MGDoFHandler &dof_handler, - const ConstraintMatrix &constraints, - const unsigned int level) - const ConstraintMatrix &constraints, ++ const ConstraintMatrix &constraints, + const unsigned int level) { typename MatrixFree::AdditionalData additional_data; additional_data.tasks_parallel_scheme = diff --cc deal.II/examples/step-4/step-4.cc index 6b01f7b7d8,38c6d6b8a6..81fa0ad715 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@@ -437,24 -437,24 +437,24 @@@ void Step4::assemble_system ( FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - // Next, we again have to loop over all - // cells and assemble local contributions. - // Note, that a cell is a quadrilateral in - // two space dimensions, but a hexahedron - // in 3D. In fact, the - // active_cell_iterator data - // type is something different, depending - // on the dimension we are in, but to the - // outside world they look alike and you - // will probably never see a difference - // although the classes that this typedef - // stands for are in fact completely - // unrelated: + // Next, we again have to loop over all + // cells and assemble local contributions. + // Note, that a cell is a quadrilateral in + // two space dimensions, but a hexahedron + // in 3D. In fact, the + // active_cell_iterator data + // type is something different, depending + // on the dimension we are in, but to the + // outside world they look alike and you + // will probably never see a difference + // although the classes that this typedef + // stands for are in fact completely + // unrelated: typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + cell = dof_handler.begin_active(), + endc = dof_handler.end(); for (; cell!=endc; ++cell) { diff --cc deal.II/examples/step-42/step-42.cc index aca98b16e2,7bbab130b1..4f1f7df1f3 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@@ -170,14 -170,14 +170,14 @@@ namespace Step4 MPI_Comm _mpi_communicator, ConditionalOStream _pcout); - void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor, + void plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor, - unsigned int &elast_points, - unsigned int &plast_points, + unsigned int &elast_points, + unsigned int &plast_points, double &yield); - void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, - SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor); + void linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, + SymmetricTensor<4,dim> &stress_strain_tensor, + SymmetricTensor<2,dim> &strain_tensor); inline SymmetricTensor<2,dim> get_strain (const FEValues &fe_values, const unsigned int shape_func, const unsigned int q_point) const; @@@ -260,36 -260,36 +260,36 @@@ } template - void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, - SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor) + void ConstitutiveLaw::linearized_plast_linear_hardening (SymmetricTensor<4,dim> &stress_strain_tensor_linearized, - SymmetricTensor<4,dim> &stress_strain_tensor, - SymmetricTensor<2,dim> &strain_tensor) ++ SymmetricTensor<4,dim> &stress_strain_tensor, ++ SymmetricTensor<2,dim> &strain_tensor) { if (dim == 3) - { - SymmetricTensor<2,dim> stress_tensor; - stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; - double tmp = E/((1+nu)*(1-2*nu)); + { + SymmetricTensor<2,dim> stress_tensor; + stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)*strain_tensor; + double tmp = E/((1+nu)*(1-2*nu)); - stress_strain_tensor = stress_strain_tensor_mu; - stress_strain_tensor_linearized = stress_strain_tensor_mu; + stress_strain_tensor = stress_strain_tensor_mu; + stress_strain_tensor_linearized = stress_strain_tensor_mu; - SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); + SymmetricTensor<2,dim> deviator_stress_tensor = deviator(stress_tensor); - double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); + double deviator_stress_tensor_norm = deviator_stress_tensor.norm (); - double beta = 1.0; - if (deviator_stress_tensor_norm >= sigma_0) - { - beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; - stress_strain_tensor *= beta; - stress_strain_tensor_linearized *= beta; - deviator_stress_tensor /= deviator_stress_tensor_norm; - stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); - } + double beta = 1.0; + if (deviator_stress_tensor_norm >= sigma_0) + { + beta = (sigma_0 + gamma)/deviator_stress_tensor_norm; + stress_strain_tensor *= beta; + stress_strain_tensor_linearized *= beta; + deviator_stress_tensor /= deviator_stress_tensor_norm; + stress_strain_tensor_linearized -= beta*2*mu*outer_product(deviator_stress_tensor, deviator_stress_tensor); + } - stress_strain_tensor += stress_strain_tensor_kappa; - stress_strain_tensor_linearized += stress_strain_tensor_kappa; - } + stress_strain_tensor += stress_strain_tensor_kappa; + stress_strain_tensor_linearized += stress_strain_tensor_kappa; + } } namespace EquationData diff --cc deal.II/examples/step-43/step-43.cc index 6876372425,d24bb00104..ec3d2f2fa3 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@@ -460,15 -460,15 +460,15 @@@ namespace Step4 template BlockSchurPreconditioner:: - BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, + BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, const InverseMatrix &Mpinv, + PreconditionerMp> &Mpinv, const PreconditionerA &Apreconditioner) - : - darcy_matrix (&S), - m_inverse (&Mpinv), - a_preconditioner (Apreconditioner), - tmp (darcy_matrix->block(1,1).m()) + : + darcy_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (darcy_matrix->block(1,1).m()) {} @@@ -548,130 -548,130 +548,130 @@@ template class TwoPhaseFlowProblem { - public: - TwoPhaseFlowProblem (const unsigned int degree); - void run (); - - private: - void setup_dofs (); - void assemble_darcy_preconditioner (); - void build_darcy_preconditioner (); - void assemble_darcy_system (); - void assemble_saturation_system (); - void assemble_saturation_matrix (); - void assemble_saturation_rhs (); - void assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, - const FEValues &darcy_fe_values, - const double global_max_u_F_prime, - const double global_S_variation, - const std::vector &local_dof_indices); - void assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, - const FEFaceValues &darcy_fe_face_values, - const std::vector &local_dof_indices); - void solve (); - void refine_mesh (const unsigned int min_grid_level, - const unsigned int max_grid_level); - void output_results () const; - - // We follow with a number of - // helper functions that are - // used in a variety of places - // throughout the program: - double get_max_u_F_prime () const; - std::pair get_extrapolated_saturation_range () const; - bool determine_whether_to_solve_for_pressure_and_velocity () const; - void project_back_saturation (); - double compute_viscosity (const std::vector &old_saturation, - const std::vector &old_old_saturation, - const std::vector > &old_saturation_grads, - const std::vector > &old_old_saturation_grads, - const std::vector > &present_darcy_values, - const double global_max_u_F_prime, - const double global_S_variation, - const double cell_diameter) const; - - - // This all is followed by the - // member variables, most of - // which are similar to the - // ones in step-31, with the - // exception of the ones that - // pertain to the macro time - // stepping for the - // velocity/pressure system: - Triangulation triangulation; - double global_Omega_diameter; - - const unsigned int degree; - - const unsigned int darcy_degree; - FESystem darcy_fe; - DoFHandler darcy_dof_handler; - ConstraintMatrix darcy_constraints; - - ConstraintMatrix darcy_preconditioner_constraints; - - TrilinosWrappers::BlockSparseMatrix darcy_matrix; - TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix; - - TrilinosWrappers::BlockVector darcy_solution; - TrilinosWrappers::BlockVector darcy_rhs; - - TrilinosWrappers::BlockVector last_computed_darcy_solution; - TrilinosWrappers::BlockVector second_last_computed_darcy_solution; - - - const unsigned int saturation_degree; - FE_Q saturation_fe; - DoFHandler saturation_dof_handler; - ConstraintMatrix saturation_constraints; - - TrilinosWrappers::SparseMatrix saturation_matrix; - - - TrilinosWrappers::Vector saturation_solution; - TrilinosWrappers::Vector old_saturation_solution; - TrilinosWrappers::Vector old_old_saturation_solution; - TrilinosWrappers::Vector saturation_rhs; - - TrilinosWrappers::Vector saturation_matching_last_computed_darcy_solution; - - const double saturation_refinement_threshold; - - double time; - const double end_time; - - double current_macro_time_step; - double old_macro_time_step; - - double time_step; - double old_time_step; - unsigned int timestep_number; - - const double viscosity; - const double porosity; - const double AOS_threshold; - - std_cxx1x::shared_ptr Amg_preconditioner; - std_cxx1x::shared_ptr Mp_preconditioner; - - bool rebuild_saturation_matrix; - - // At the very end we declare a - // variable that denotes the - // material model. Compared to - // step-21, we do this here as - // a member variable since we - // will want to use it in a - // variety of places and so - // having a central place where - // such a variable is declared - // will make it simpler to - // replace one class by another - // (e.g. replace - // RandomMedium::KInverse by - // SingleCurvingCrack::KInverse). - const RandomMedium::KInverse k_inverse; + public: + TwoPhaseFlowProblem (const unsigned int degree); + void run (); + + private: + void setup_dofs (); + void assemble_darcy_preconditioner (); + void build_darcy_preconditioner (); + void assemble_darcy_system (); + void assemble_saturation_system (); + void assemble_saturation_matrix (); + void assemble_saturation_rhs (); + void assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, + const FEValues &darcy_fe_values, + const double global_max_u_F_prime, + const double global_S_variation, + const std::vector &local_dof_indices); + void assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, + const FEFaceValues &darcy_fe_face_values, + const std::vector &local_dof_indices); + void solve (); + void refine_mesh (const unsigned int min_grid_level, + const unsigned int max_grid_level); + void output_results () const; + + // We follow with a number of + // helper functions that are + // used in a variety of places + // throughout the program: + double get_max_u_F_prime () const; + std::pair get_extrapolated_saturation_range () const; + bool determine_whether_to_solve_for_pressure_and_velocity () const; + void project_back_saturation (); + double compute_viscosity (const std::vector &old_saturation, + const std::vector &old_old_saturation, - const std::vector > &old_saturation_grads, - const std::vector > &old_old_saturation_grads, ++ const std::vector > &old_saturation_grads, ++ const std::vector > &old_old_saturation_grads, + const std::vector > &present_darcy_values, + const double global_max_u_F_prime, + const double global_S_variation, + const double cell_diameter) const; + + + // This all is followed by the + // member variables, most of + // which are similar to the + // ones in step-31, with the + // exception of the ones that + // pertain to the macro time + // stepping for the + // velocity/pressure system: + Triangulation triangulation; + double global_Omega_diameter; + + const unsigned int degree; + + const unsigned int darcy_degree; + FESystem darcy_fe; + DoFHandler darcy_dof_handler; + ConstraintMatrix darcy_constraints; + + ConstraintMatrix darcy_preconditioner_constraints; + + TrilinosWrappers::BlockSparseMatrix darcy_matrix; + TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix; + + TrilinosWrappers::BlockVector darcy_solution; + TrilinosWrappers::BlockVector darcy_rhs; + + TrilinosWrappers::BlockVector last_computed_darcy_solution; + TrilinosWrappers::BlockVector second_last_computed_darcy_solution; + + + const unsigned int saturation_degree; + FE_Q saturation_fe; + DoFHandler saturation_dof_handler; + ConstraintMatrix saturation_constraints; + + TrilinosWrappers::SparseMatrix saturation_matrix; + + + TrilinosWrappers::Vector saturation_solution; + TrilinosWrappers::Vector old_saturation_solution; + TrilinosWrappers::Vector old_old_saturation_solution; + TrilinosWrappers::Vector saturation_rhs; + + TrilinosWrappers::Vector saturation_matching_last_computed_darcy_solution; + + const double saturation_refinement_threshold; + + double time; + const double end_time; + + double current_macro_time_step; + double old_macro_time_step; + + double time_step; + double old_time_step; + unsigned int timestep_number; + + const double viscosity; + const double porosity; + const double AOS_threshold; + + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + + bool rebuild_saturation_matrix; + + // At the very end we declare a + // variable that denotes the + // material model. Compared to + // step-21, we do this here as + // a member variable since we + // will want to use it in a + // variety of places and so + // having a central place where + // such a variable is declared + // will make it simpler to + // replace one class by another + // (e.g. replace + // RandomMedium::KInverse by + // SingleCurvingCrack::KInverse). + const RandomMedium::KInverse k_inverse; }; diff --cc deal.II/examples/step-5/step-5.cc index ff03e580f5,fc5f1e582e..0b059c108d --- a/deal.II/examples/step-5/step-5.cc +++ b/deal.II/examples/step-5/step-5.cc @@@ -442,61 -442,61 +442,61 @@@ void Step5::assemble_system ( FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - // Here is one difference: for this - // program, we will again use a - // constant right hand side - // function and zero boundary - // values, but a variable - // coefficient. We have already - // declared the class that - // represents this coefficient - // above, so we only have to - // declare a corresponding object - // here. - // - // Then, below, we will ask the - // coefficient function object - // to compute the values of the - // coefficient at all quadrature - // points on one cell at once. The - // reason for this is that, if you - // look back at how we did this in - // step-4, you will realize that we - // called the function computing - // the right hand side value inside - // nested loops over all degrees of - // freedom and over all quadrature - // points, - // i.e. dofs_per_cell*n_q_points - // times. For the coefficient that - // is used inside the matrix, this - // would actually be - // dofs_per_cell*dofs_per_cell*n_q_points. On - // the other hand, the function - // will of course return the same - // value every time it is called - // with the same quadrature point, - // independently of what shape - // function we presently treat; - // secondly, these are virtual - // function calls, so are rather - // expensive. Obviously, there are - // only n_q_point different values, - // and we shouldn't call the - // function more often than - // that. Or, even better than this, - // compute all of these values at - // once, and get away with a single - // function call per cell. - // - // This is exactly what we are - // going to do. For this, we need - // some space to store the values - // in. We therefore also have to - // declare an array to hold these - // values: + // Here is one difference: for this + // program, we will again use a + // constant right hand side + // function and zero boundary + // values, but a variable + // coefficient. We have already + // declared the class that + // represents this coefficient + // above, so we only have to + // declare a corresponding object + // here. + // + // Then, below, we will ask the + // coefficient function object + // to compute the values of the + // coefficient at all quadrature + // points on one cell at once. The + // reason for this is that, if you + // look back at how we did this in + // step-4, you will realize that we + // called the function computing + // the right hand side value inside + // nested loops over all degrees of + // freedom and over all quadrature + // points, + // i.e. dofs_per_cell*n_q_points + // times. For the coefficient that + // is used inside the matrix, this + // would actually be + // dofs_per_cell*dofs_per_cell*n_q_points. On + // the other hand, the function + // will of course return the same + // value every time it is called + // with the same quadrature point, + // independently of what shape + // function we presently treat; + // secondly, these are virtual + // function calls, so are rather + // expensive. Obviously, there are + // only n_q_point different values, + // and we shouldn't call the + // function more often than + // that. Or, even better than this, + // compute all of these values at + // once, and get away with a single + // function call per cell. + // + // This is exactly what we are + // going to do. For this, we need + // some space to store the values + // in. We therefore also have to + // declare an array to hold these + // values: const Coefficient coefficient; std::vector coefficient_values (n_q_points); diff --cc deal.II/examples/step-9/step-9.cc index 40d12c9b88,ba43fa9e49..3ae105a770 --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@@ -510,129 -510,129 +510,129 @@@ namespace Step - // @sect3{GradientEstimation class declaration} - - // Now, finally, here comes the class - // that will compute the difference - // approximation of the gradient on - // each cell and weighs that with a - // power of the mesh size, as - // described in the introduction. - // This class is a simple version of - // the DerivativeApproximation - // class in the library, that uses - // similar techniques to obtain - // finite difference approximations - // of the gradient of a finite - // element field, or if higher - // derivatives. - // - // The - // class has one public static - // function estimate that is - // called to compute a vector of - // error indicators, and one private - // function that does the actual work - // on an interval of all active - // cells. The latter is called by the - // first one in order to be able to - // do the computations in parallel if - // your computer has more than one - // processor. While the first - // function accepts as parameter a - // vector into which the error - // indicator is written for each - // cell. This vector is passed on to - // the second function that actually - // computes the error indicators on - // some cells, and the respective - // elements of the vector are - // written. By the way, we made it - // somewhat of a convention to use - // vectors of floats for error - // indicators rather than the common - // vectors of doubles, as the - // additional accuracy is not - // necessary for estimated values. - // - // In addition to these two - // functions, the class declares to - // exceptions which are raised when a - // cell has no neighbors in each of - // the space directions (in which - // case the matrix described in the - // introduction would be singular and - // can't be inverted), while the - // other one is used in the more - // common case of invalid parameters - // to a function, namely a vector of - // wrong size. - // - // Two annotations to this class are - // still in order: the first is that - // the class has no non-static member - // functions or variables, so this is - // not really a class, but rather - // serves the purpose of a - // namespace in C++. The reason - // that we chose a class over a - // namespace is that this way we can - // declare functions that are - // private, i.e. visible to the - // outside world but not - // callable. This can be done with - // namespaces as well, if one - // declares some functions in header - // files in the namespace and - // implements these and other - // functions in the implementation - // file. The functions not declared - // in the header file are still in - // the namespace but are not callable - // from outside. However, as we have - // only one file here, it is not - // possible to hide functions in the - // present case. - // - // The second is that the dimension - // template parameter is attached to - // the function rather than to the - // class itself. This way, you don't - // have to specify the template - // parameter yourself as in most - // other cases, but the compiler can - // figure its value out itself from - // the dimension of the DoF handler - // object that one passes as first - // argument. - // - // Finally note that the - // IndexInterval typedef is - // introduced as a convenient - // abbreviation for an otherwise - // lengthy type name. + // @sect3{GradientEstimation class declaration} + + // Now, finally, here comes the class + // that will compute the difference + // approximation of the gradient on + // each cell and weighs that with a + // power of the mesh size, as + // described in the introduction. + // This class is a simple version of + // the DerivativeApproximation + // class in the library, that uses + // similar techniques to obtain + // finite difference approximations + // of the gradient of a finite + // element field, or if higher + // derivatives. + // + // The + // class has one public static + // function estimate that is + // called to compute a vector of + // error indicators, and one private + // function that does the actual work + // on an interval of all active + // cells. The latter is called by the + // first one in order to be able to + // do the computations in parallel if + // your computer has more than one + // processor. While the first + // function accepts as parameter a + // vector into which the error + // indicator is written for each + // cell. This vector is passed on to + // the second function that actually + // computes the error indicators on + // some cells, and the respective + // elements of the vector are + // written. By the way, we made it + // somewhat of a convention to use + // vectors of floats for error + // indicators rather than the common + // vectors of doubles, as the + // additional accuracy is not + // necessary for estimated values. + // + // In addition to these two + // functions, the class declares to + // exceptions which are raised when a + // cell has no neighbors in each of + // the space directions (in which + // case the matrix described in the + // introduction would be singular and + // can't be inverted), while the + // other one is used in the more + // common case of invalid parameters + // to a function, namely a vector of + // wrong size. + // + // Two annotations to this class are + // still in order: the first is that + // the class has no non-static member + // functions or variables, so this is + // not really a class, but rather + // serves the purpose of a + // namespace in C++. The reason + // that we chose a class over a + // namespace is that this way we can + // declare functions that are + // private, i.e. visible to the + // outside world but not + // callable. This can be done with + // namespaces as well, if one + // declares some functions in header + // files in the namespace and + // implements these and other + // functions in the implementation + // file. The functions not declared + // in the header file are still in + // the namespace but are not callable + // from outside. However, as we have + // only one file here, it is not + // possible to hide functions in the + // present case. + // + // The second is that the dimension + // template parameter is attached to + // the function rather than to the + // class itself. This way, you don't + // have to specify the template + // parameter yourself as in most + // other cases, but the compiler can + // figure its value out itself from + // the dimension of the DoF handler + // object that one passes as first + // argument. + // + // Finally note that the + // IndexInterval typedef is + // introduced as a convenient + // abbreviation for an otherwise + // lengthy type name. class GradientEstimation { - public: - template - static void estimate (const DoFHandler &dof, - const Vector &solution, - Vector &error_per_cell); - - DeclException2 (ExcInvalidVectorLength, - int, int, - << "Vector has length " << arg1 << ", but should have " - << arg2); - DeclException0 (ExcInsufficientDirections); - - private: - typedef std::pair IndexInterval; - - template - static void estimate_interval (const DoFHandler &dof, - const Vector &solution, - const IndexInterval &index_interval, - Vector &error_per_cell); + public: + template + static void estimate (const DoFHandler &dof, - const Vector &solution, ++ const Vector &solution, + Vector &error_per_cell); + + DeclException2 (ExcInvalidVectorLength, + int, int, + << "Vector has length " << arg1 << ", but should have " + << arg2); + DeclException0 (ExcInsufficientDirections); + + private: + typedef std::pair IndexInterval; + + template + static void estimate_interval (const DoFHandler &dof, - const Vector &solution, ++ const Vector &solution, + const IndexInterval &index_interval, + Vector &error_per_cell); }; @@@ -1459,24 -1459,24 +1459,24 @@@ template void GradientEstimation::estimate (const DoFHandler &dof_handler, - const Vector &solution, + const Vector &solution, Vector &error_per_cell) { - // Before starting with the work, - // we check that the vector into - // which the results are written, - // has the right size. It is a - // common error that such - // parameters have the wrong size, - // but the resulting damage by not - // catching these errors are very - // subtle as they are usually - // corruption of data somewhere in - // memory. Often, the problems - // emerging from this are not - // reproducible, and we found that - // it is well worth the effort to - // check for such things. + // Before starting with the work, + // we check that the vector into + // which the results are written, + // has the right size. It is a + // common error that such + // parameters have the wrong size, + // but the resulting damage by not + // catching these errors are very + // subtle as they are usually + // corruption of data somewhere in + // memory. Often, the problems + // emerging from this are not + // reproducible, and we found that + // it is well worth the effort to + // check for such things. Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(), ExcInvalidVectorLength (error_per_cell.size(), dof_handler.get_tria().n_active_cells())); diff --cc deal.II/include/deal.II/base/conditional_ostream.h index a4c799f48b,d786f5d02e..753ade963e --- a/deal.II/include/deal.II/base/conditional_ostream.h +++ b/deal.II/include/deal.II/base/conditional_ostream.h @@@ -80,81 -80,81 +80,81 @@@ DEAL_II_NAMESPACE_OPE */ class ConditionalOStream { - public: - /** - * Constructor. Set the stream to which - * we want to write, and the condition - * based on which writes are actually - * forwarded. Per default the condition - * of an object is active. - */ - ConditionalOStream (std::ostream &stream, - const bool active = true); - - /** - * Depending on the - * active flag set the - * condition of this stream to - * active (true) or non-active - * (false). An object of this - * class prints to cout - * if and only if its condition - * is active. - */ - void set_condition (const bool active); - - /** - * Return the condition of the object. - */ - bool is_active() const; - - /** - * Return a reference to the stream - * currently in use. - */ - std::ostream & get_stream () const; - - /** - * Output a constant something through - * this stream. This function must be @p - * const so that member objects of this - * type can also be used from @p const - * member functions of the surrounding - * class. - */ - template - const ConditionalOStream & - operator << (const T &t) const; - - /** - * Treat ostream manipulators. This - * function must be @p const so that - * member objects of this type can also - * be used from @p const member functions - * of the surrounding class. - * - * Note that compilers want to see this - * treated differently from the general - * template above since functions like @p - * std::endl are actually overloaded and - * can't be bound directly to a template - * type. - */ - const ConditionalOStream & - operator<< (std::ostream& (*p) (std::ostream&)) const; - - private: - /** - * Reference to the stream we - * want to write to. - */ - std::ostream &output_stream; - - /** - * Stores the actual condition - * the object is in. - */ - bool active_flag; + public: + /** + * Constructor. Set the stream to which + * we want to write, and the condition + * based on which writes are actually + * forwarded. Per default the condition + * of an object is active. + */ + ConditionalOStream (std::ostream &stream, + const bool active = true); + + /** + * Depending on the + * active flag set the + * condition of this stream to + * active (true) or non-active + * (false). An object of this + * class prints to cout + * if and only if its condition + * is active. + */ + void set_condition (const bool active); + + /** + * Return the condition of the object. + */ + bool is_active() const; + + /** + * Return a reference to the stream + * currently in use. + */ + std::ostream &get_stream () const; + + /** + * Output a constant something through + * this stream. This function must be @p + * const so that member objects of this + * type can also be used from @p const + * member functions of the surrounding + * class. + */ + template + const ConditionalOStream & + operator << (const T &t) const; + + /** + * Treat ostream manipulators. This + * function must be @p const so that + * member objects of this type can also + * be used from @p const member functions + * of the surrounding class. + * + * Note that compilers want to see this + * treated differently from the general + * template above since functions like @p + * std::endl are actually overloaded and + * can't be bound directly to a template + * type. + */ + const ConditionalOStream & + operator<< (std::ostream& (*p) (std::ostream &)) const; + + private: + /** + * Reference to the stream we + * want to write to. + */ - std::ostream &output_stream; ++ std::ostream &output_stream; + + /** + * Stores the actual condition + * the object is in. + */ + bool active_flag; }; diff --cc deal.II/include/deal.II/base/data_out_base.h index b42e689aa9,f1ff018f5c..a8124528fc --- a/deal.II/include/deal.II/base/data_out_base.h +++ b/deal.II/include/deal.II/base/data_out_base.h @@@ -2178,640 -2179,640 +2179,640 @@@ private template class DataOutInterface : private DataOutBase { - public: - /* - * Import a few names that were - * previously in this class and have then - * moved to the base class. Since the - * base class is inherited from - * privately, we need to re-import these - * symbols to make sure that references - * to DataOutInterface::XXX - * remain valid. - */ - using DataOutBase::OutputFormat; - using DataOutBase::default_format; - using DataOutBase::dx; - using DataOutBase::gnuplot; - using DataOutBase::povray; - using DataOutBase::eps; - using DataOutBase::tecplot; - using DataOutBase::tecplot_binary; - using DataOutBase::vtk; - using DataOutBase::vtu; - using DataOutBase::deal_II_intermediate; - using DataOutBase::parse_output_format; - using DataOutBase::get_output_format_names; - using DataOutBase::determine_intermediate_format_dimensions; - - /** - * Constructor. - */ - DataOutInterface (); - - /** - * Destructor. Does nothing, but is - * declared virtual since this class has - * virtual functions. - */ - virtual ~DataOutInterface (); - - /** - * Obtain data through get_patches() - * and write it to out - * in OpenDX format. See - * DataOutBase::write_dx. - */ - void write_dx (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in EPS format. See - * DataOutBase::write_eps. - */ - void write_eps (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in GMV format. See - * DataOutBase::write_gmv. - */ - void write_gmv (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in GNUPLOT format. See - * DataOutBase::write_gnuplot. - */ - void write_gnuplot (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in POVRAY format. See - * DataOutBase::write_povray. - */ - void write_povray (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in Tecplot format. See - * DataOutBase::write_tecplot. - */ - void write_tecplot (std::ostream &out) const; - - /** - * Obtain data through - * get_patches() and write it in - * the Tecplot binary output - * format. Note that the name of - * the output file must be - * specified through the - * TecplotFlags interface. - */ - void write_tecplot_binary (std::ostream &out) const; - - /** - * Obtain data through - * get_patches() and write it to - * out in UCD format for - * AVS. See - * DataOutBase::write_ucd. - */ - void write_ucd (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in Vtk format. See - * DataOutBase::write_vtk. - */ - void write_vtk (std::ostream &out) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in Vtu (VTK's XML) format. See - * DataOutBase::write_vtu. - * - * Some visualization programs, - * such as ParaView, can read - * several separate VTU files to - * parallelize visualization. In - * that case, you need a - * .pvtu file that - * describes which VTU files form - * a group. The - * DataOutInterface::write_pvtu_record() - * function can generate such a - * master record. Likewise, - * DataOutInterface::write_visit_record() - * does the same for VisIt. Finally, - * DataOutInterface::write_pvd_record() - * can be used to group together - * the files that jointly make up - * a time dependent simulation. - */ - void write_vtu (std::ostream &out) const; - - /** - * Collective MPI call to write the - * solution from all participating nodes - * (those in the given communicator) to a - * single compressed .vtu file on a - * shared file system. The communicator - * can be a sub communicator of the one - * used by the computation. This routine - * uses MPI I/O to achieve high - * performance on parallel filesystems. - * Also see - * DataOutInterface::write_vtu(). - */ - void write_vtu_in_parallel (const char* filename, MPI_Comm comm) const; - - /** - * Some visualization programs, such as - * ParaView, can read several separate - * VTU files to parallelize - * visualization. In that case, you need - * a .pvtu file that - * describes which VTU files (written, - * for example, through the write_vtu() - * function) form a group. The current - * function can generate such a master - * record. - * - * The file so written contains a list of - * (scalar or vector) fields whose values - * are described by the individual files - * that comprise the set of parallel VTU - * files along with the names of these - * files. This function gets the names - * and types of fields through the - * get_patches() function of this class - * like all the other write_xxx() - * functions. The second argument to this - * function specifies the names of the - * files that form the parallel set. - * - * @note See DataOutBase::write_vtu for - * writing each piece. Also note that - * only one parallel process needs to - * call the current function, listing the - * names of the files written by all - * parallel processes. - * - * @note The use of this function is - * explained in step-40. - * - * @note In order to tell Paraview to - * group together multiple pvtu - * files that each describe one time - * step of a time dependent simulation, - * see the - * DataOutInterface::write_pvd_record() - * function. - * - * @note At the time of writing, - * the other big VTK-based - * visualization program, VisIt, - * can not read pvtu - * records. However, it can read - * visit records as written by - * the write_visit_record() - * function. - */ - void write_pvtu_record (std::ostream &out, - const std::vector &piece_names) const; - - /** - * In ParaView it is possible to visualize time-dependent - * data tagged with the current - * integration time of a time dependent simulation. To use this - * feature you need a .pvd - * file that describes which VTU or PVTU file - * belongs to which timestep. This function writes a file that - * provides this mapping, i.e., it takes a list of pairs each of - * which indicates a particular time instant and the corresponding - * file that contains the graphical data for this time instant. - * - * A typical use case, in program that computes a time dependent - * solution, would be the following (time and - * time_step are member variables of the class with types - * double and unsigned int, respectively; - * the variable times_and_names is of type - * std::vector@ @>): - * - * @code - * template - * void MyEquation::output_results () const - * { - * DataOut data_out; - * - * data_out.attach_dof_handler (dof_handler); - * data_out.add_data_vector (solution, "U"); - * data_out.build_patches (); - * - * const std::string filename = "solution-" + - * Utilities::int_to_string (timestep_number, 3) + - * ".vtu"; - * std::ofstream output (filename.c_str()); - * data_out.write_vtu (output); - * - * times_and_names.push_back (std::pair (time, filename)); - * std::ofstream pvd_output ("solution.pvd"); - * data_out.write_pvd_record (pvd_output, times_and_names); - * } - * @endcode - * - * @note See DataOutBase::write_vtu or - * DataOutInterface::write_pvtu_record for - * writing solutions at each timestep. - * - * @note The second element of each pair, i.e., the file in which - * the graphical data for each time is stored, may itself be again - * a file that references other files. For example, it could be - * the name for a .pvtu file that references multiple - * parts of a parallel computation. - * - * @author Marco Engelhard, 2012 - */ - void write_pvd_record (std::ostream &out, - const std::vector > ×_and_names) const; - - /** - * This function is the exact - * equivalent of the - * write_pvtu_record() function - * but for the VisIt - * visualization program. See - * there for the purpose of this - * function. - * - * This function is documented - * in the "Creating a master file - * for parallel" section (section 5.7) - * of the "Getting data into VisIt" - * report that can be found here: - * https://wci.llnl.gov/codes/visit/2.0.0/GettingDataIntoVisIt2.0.0.pdf - */ - void write_visit_record (std::ostream &out, - const std::vector &piece_names) const; - - /** - * Obtain data through get_patches() - * and write it to out - * in deal.II intermediate - * format. See - * DataOutBase::write_deal_II_intermediate. - * - * Note that the intermediate - * format is what its name - * suggests: a direct - * representation of internal - * data. It isn't standardized - * and will change whenever we - * change our internal - * representation. You can only - * expect to process files - * written in this format using - * the same version of deal.II - * that was used for writing. - */ - void write_deal_II_intermediate (std::ostream &out) const; - - XDMFEntry create_xdmf_entry (const char *h5_filename, - const double cur_time, - MPI_Comm comm) const; - - void write_xdmf_file (const std::vector &entries, - const char *filename, - MPI_Comm comm) const; - - void write_hdf5_parallel (const char* filename, MPI_Comm comm) const; - /** - * Write data and grid to out - * according to the given data - * format. This function simply - * calls the appropriate - * write_* function. If no - * output format is requested, - * the default_format is - * written. - * - * An error occurs if no format - * is provided and the default - * format is default_format. - */ - void write (std::ostream &out, - const OutputFormat output_format = default_format) const; - - /** - * Set the default format. The - * value set here is used - * anytime, output for format - * default_format is - * requested. - */ - void set_default_format (const OutputFormat default_format); - - /** - * Set the flags to be used for - * output in OpenDX format. - */ - void set_flags (const DXFlags &dx_flags); - - /** - * Set the flags to be used for - * output in UCD format. - */ - void set_flags (const UcdFlags &ucd_flags); - - /** - * Set the flags to be used for - * output in GNUPLOT format. - */ - void set_flags (const GnuplotFlags &gnuplot_flags); - - /** - * Set the flags to be used for - * output in POVRAY format. - */ - void set_flags (const PovrayFlags &povray_flags); - - /** - * Set the flags to be used for - * output in EPS output. - */ - void set_flags (const EpsFlags &eps_flags); - - /** - * Set the flags to be used for - * output in GMV format. - */ - void set_flags (const GmvFlags &gmv_flags); - - /** - * Set the flags to be used for - * output in Tecplot format. - */ - void set_flags (const TecplotFlags &tecplot_flags); - - /** - * Set the flags to be used for - * output in VTK format. - */ - void set_flags (const VtkFlags &vtk_flags); - - /** - * Set the flags to be used for output in - * deal.II intermediate format. - */ - void set_flags (const Deal_II_IntermediateFlags &deal_II_intermediate_flags); - - /** - * A function that returns the same - * string as the respective function in - * the base class does; the only - * exception being that if the parameter - * is omitted, then the value for the - * present default format is returned, - * i.e. the correct suffix for the format - * that was set through - * set_default_format() or - * parse_parameters() before calling this - * function. - */ - std::string - default_suffix (const OutputFormat output_format = default_format) const; - - /** - * Declare parameters for all - * output formats by declaring - * subsections within the - * parameter file for each output - * format and call the respective - * declare_parameters - * functions of the flag classes - * for each output format. - * - * Some of the declared - * subsections may not contain - * entries, if the respective - * format does not export any - * flags. - * - * Note that the top-level - * parameters denoting the number - * of subdivisions per patch and - * the output format are not - * declared, since they are only - * passed to virtual functions - * and are not stored inside - * objects of this type. You have - * to declare them yourself. - */ - static void declare_parameters (ParameterHandler &prm); - - /** - * Read the parameters declared - * in declare_parameters and - * set the flags for the output - * formats accordingly. - * - * The flags thus obtained - * overwrite all previous - * contents of the flag objects - * as default-constructed or set - * by the set_flags() function. - */ - void parse_parameters (ParameterHandler &prm); - - /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this - * object. Since sometimes - * the size of objects can - * not be determined exactly - * (for example: what is the - * memory consumption of an - * STL std::map type with a - * certain number of - * elements?), this is only - * an estimate. however often - * quite close to the true - * value. - */ - std::size_t memory_consumption () const; - - protected: - /** - * This is the abstract function - * through which derived classes - * propagate preprocessed data in - * the form of Patch - * structures (declared in the - * base class DataOutBase) to - * the actual output - * function. You need to overload - * this function to allow the - * output functions to know what - * they shall print. - */ - virtual - const std::vector > & - get_patches () const = 0; - - /** - * Abstract virtual function - * through which the names of - * data sets are obtained by the - * output functions of the base - * class. - */ - virtual - std::vector - get_dataset_names () const = 0; - - /** - * This functions returns - * information about how the - * individual components of - * output files that consist of - * more than one data set are to - * be interpreted. - * - * It returns a list of index - * pairs and corresponding name - * indicating which components of - * the output are to be - * considered vector-valued - * rather than just a collection - * of scalar data. The index - * pairs are inclusive; for - * example, if we have a Stokes - * problem in 2d with components - * (u,v,p), then the - * corresponding vector data - * range should be (0,1), and the - * returned list would consist of - * only a single element with a - * tuple such as (0,1,"velocity"). - * - * Since some of the derived - * classes do not know about - * vector data, this function has - * a default implementation that - * simply returns an empty - * string, meaning that all data - * is to be considered a - * collection of scalar fields. - */ - virtual - std::vector > - get_vector_data_ranges () const; - - /** - * The default number of - * subdivisions for patches. This - * is filled by parse_parameters() - * and should be obeyed by - * build_patches() in derived - * classes. - */ - unsigned int default_subdivisions; + public: + /* + * Import a few names that were + * previously in this class and have then + * moved to the base class. Since the + * base class is inherited from + * privately, we need to re-import these + * symbols to make sure that references + * to DataOutInterface::XXX + * remain valid. + */ + using DataOutBase::OutputFormat; + using DataOutBase::default_format; + using DataOutBase::dx; + using DataOutBase::gnuplot; + using DataOutBase::povray; + using DataOutBase::eps; + using DataOutBase::tecplot; + using DataOutBase::tecplot_binary; + using DataOutBase::vtk; + using DataOutBase::vtu; + using DataOutBase::deal_II_intermediate; + using DataOutBase::parse_output_format; + using DataOutBase::get_output_format_names; + using DataOutBase::determine_intermediate_format_dimensions; + + /** + * Constructor. + */ + DataOutInterface (); + + /** + * Destructor. Does nothing, but is + * declared virtual since this class has + * virtual functions. + */ + virtual ~DataOutInterface (); + + /** + * Obtain data through get_patches() + * and write it to out + * in OpenDX format. See + * DataOutBase::write_dx. + */ + void write_dx (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in EPS format. See + * DataOutBase::write_eps. + */ + void write_eps (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in GMV format. See + * DataOutBase::write_gmv. + */ + void write_gmv (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in GNUPLOT format. See + * DataOutBase::write_gnuplot. + */ + void write_gnuplot (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in POVRAY format. See + * DataOutBase::write_povray. + */ + void write_povray (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in Tecplot format. See + * DataOutBase::write_tecplot. + */ + void write_tecplot (std::ostream &out) const; + + /** + * Obtain data through + * get_patches() and write it in + * the Tecplot binary output + * format. Note that the name of + * the output file must be + * specified through the + * TecplotFlags interface. + */ + void write_tecplot_binary (std::ostream &out) const; + + /** + * Obtain data through + * get_patches() and write it to + * out in UCD format for + * AVS. See + * DataOutBase::write_ucd. + */ + void write_ucd (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in Vtk format. See + * DataOutBase::write_vtk. + */ + void write_vtk (std::ostream &out) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in Vtu (VTK's XML) format. See + * DataOutBase::write_vtu. + * + * Some visualization programs, + * such as ParaView, can read + * several separate VTU files to + * parallelize visualization. In + * that case, you need a + * .pvtu file that + * describes which VTU files form + * a group. The + * DataOutInterface::write_pvtu_record() + * function can generate such a + * master record. Likewise, + * DataOutInterface::write_visit_record() + * does the same for VisIt. Finally, + * DataOutInterface::write_pvd_record() + * can be used to group together + * the files that jointly make up + * a time dependent simulation. + */ + void write_vtu (std::ostream &out) const; + + /** + * Collective MPI call to write the + * solution from all participating nodes + * (those in the given communicator) to a + * single compressed .vtu file on a + * shared file system. The communicator + * can be a sub communicator of the one + * used by the computation. This routine + * uses MPI I/O to achieve high + * performance on parallel filesystems. + * Also see + * DataOutInterface::write_vtu(). + */ + void write_vtu_in_parallel (const char *filename, MPI_Comm comm) const; + + /** + * Some visualization programs, such as + * ParaView, can read several separate + * VTU files to parallelize + * visualization. In that case, you need + * a .pvtu file that + * describes which VTU files (written, + * for example, through the write_vtu() + * function) form a group. The current + * function can generate such a master + * record. + * + * The file so written contains a list of + * (scalar or vector) fields whose values + * are described by the individual files + * that comprise the set of parallel VTU + * files along with the names of these + * files. This function gets the names + * and types of fields through the + * get_patches() function of this class + * like all the other write_xxx() + * functions. The second argument to this + * function specifies the names of the + * files that form the parallel set. + * + * @note See DataOutBase::write_vtu for + * writing each piece. Also note that + * only one parallel process needs to + * call the current function, listing the + * names of the files written by all + * parallel processes. + * + * @note The use of this function is + * explained in step-40. + * + * @note In order to tell Paraview to + * group together multiple pvtu + * files that each describe one time + * step of a time dependent simulation, + * see the + * DataOutInterface::write_pvd_record() + * function. + * + * @note At the time of writing, + * the other big VTK-based + * visualization program, VisIt, + * can not read pvtu + * records. However, it can read + * visit records as written by + * the write_visit_record() + * function. + */ + void write_pvtu_record (std::ostream &out, + const std::vector &piece_names) const; + + /** + * In ParaView it is possible to visualize time-dependent + * data tagged with the current + * integration time of a time dependent simulation. To use this + * feature you need a .pvd + * file that describes which VTU or PVTU file + * belongs to which timestep. This function writes a file that + * provides this mapping, i.e., it takes a list of pairs each of + * which indicates a particular time instant and the corresponding + * file that contains the graphical data for this time instant. + * + * A typical use case, in program that computes a time dependent + * solution, would be the following (time and + * time_step are member variables of the class with types + * double and unsigned int, respectively; + * the variable times_and_names is of type + * std::vector@ @>): + * + * @code + * template + * void MyEquation::output_results () const + * { + * DataOut data_out; + * + * data_out.attach_dof_handler (dof_handler); + * data_out.add_data_vector (solution, "U"); + * data_out.build_patches (); + * + * const std::string filename = "solution-" + + * Utilities::int_to_string (timestep_number, 3) + + * ".vtu"; + * std::ofstream output (filename.c_str()); + * data_out.write_vtu (output); + * + * times_and_names.push_back (std::pair (time, filename)); + * std::ofstream pvd_output ("solution.pvd"); + * data_out.write_pvd_record (pvd_output, times_and_names); + * } + * @endcode + * + * @note See DataOutBase::write_vtu or + * DataOutInterface::write_pvtu_record for + * writing solutions at each timestep. + * + * @note The second element of each pair, i.e., the file in which + * the graphical data for each time is stored, may itself be again + * a file that references other files. For example, it could be + * the name for a .pvtu file that references multiple + * parts of a parallel computation. + * + * @author Marco Engelhard, 2012 + */ + void write_pvd_record (std::ostream &out, - const std::vector > ×_and_names) const; ++ const std::vector > ×_and_names) const; + + /** + * This function is the exact + * equivalent of the + * write_pvtu_record() function + * but for the VisIt + * visualization program. See + * there for the purpose of this + * function. + * + * This function is documented + * in the "Creating a master file + * for parallel" section (section 5.7) + * of the "Getting data into VisIt" + * report that can be found here: + * https://wci.llnl.gov/codes/visit/2.0.0/GettingDataIntoVisIt2.0.0.pdf + */ + void write_visit_record (std::ostream &out, + const std::vector &piece_names) const; + + /** + * Obtain data through get_patches() + * and write it to out + * in deal.II intermediate + * format. See + * DataOutBase::write_deal_II_intermediate. + * + * Note that the intermediate + * format is what its name + * suggests: a direct + * representation of internal + * data. It isn't standardized + * and will change whenever we + * change our internal + * representation. You can only + * expect to process files + * written in this format using + * the same version of deal.II + * that was used for writing. + */ + void write_deal_II_intermediate (std::ostream &out) const; + + XDMFEntry create_xdmf_entry (const char *h5_filename, + const double cur_time, + MPI_Comm comm) const; + + void write_xdmf_file (const std::vector &entries, + const char *filename, + MPI_Comm comm) const; + + void write_hdf5_parallel (const char *filename, MPI_Comm comm) const; + /** + * Write data and grid to out + * according to the given data + * format. This function simply + * calls the appropriate + * write_* function. If no + * output format is requested, + * the default_format is + * written. + * + * An error occurs if no format + * is provided and the default + * format is default_format. + */ + void write (std::ostream &out, + const OutputFormat output_format = default_format) const; + + /** + * Set the default format. The + * value set here is used + * anytime, output for format + * default_format is + * requested. + */ + void set_default_format (const OutputFormat default_format); + + /** + * Set the flags to be used for + * output in OpenDX format. + */ + void set_flags (const DXFlags &dx_flags); + + /** + * Set the flags to be used for + * output in UCD format. + */ + void set_flags (const UcdFlags &ucd_flags); + + /** + * Set the flags to be used for + * output in GNUPLOT format. + */ + void set_flags (const GnuplotFlags &gnuplot_flags); + + /** + * Set the flags to be used for + * output in POVRAY format. + */ + void set_flags (const PovrayFlags &povray_flags); + + /** + * Set the flags to be used for + * output in EPS output. + */ + void set_flags (const EpsFlags &eps_flags); + + /** + * Set the flags to be used for + * output in GMV format. + */ + void set_flags (const GmvFlags &gmv_flags); + + /** + * Set the flags to be used for + * output in Tecplot format. + */ + void set_flags (const TecplotFlags &tecplot_flags); + + /** + * Set the flags to be used for + * output in VTK format. + */ + void set_flags (const VtkFlags &vtk_flags); + + /** + * Set the flags to be used for output in + * deal.II intermediate format. + */ + void set_flags (const Deal_II_IntermediateFlags &deal_II_intermediate_flags); + + /** + * A function that returns the same + * string as the respective function in + * the base class does; the only + * exception being that if the parameter + * is omitted, then the value for the + * present default format is returned, + * i.e. the correct suffix for the format + * that was set through + * set_default_format() or + * parse_parameters() before calling this + * function. + */ + std::string + default_suffix (const OutputFormat output_format = default_format) const; + + /** + * Declare parameters for all + * output formats by declaring + * subsections within the + * parameter file for each output + * format and call the respective + * declare_parameters + * functions of the flag classes + * for each output format. + * + * Some of the declared + * subsections may not contain + * entries, if the respective + * format does not export any + * flags. + * + * Note that the top-level + * parameters denoting the number + * of subdivisions per patch and + * the output format are not + * declared, since they are only + * passed to virtual functions + * and are not stored inside + * objects of this type. You have + * to declare them yourself. + */ + static void declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters declared + * in declare_parameters and + * set the flags for the output + * formats accordingly. + * + * The flags thus obtained + * overwrite all previous + * contents of the flag objects + * as default-constructed or set + * by the set_flags() function. + */ + void parse_parameters (ParameterHandler &prm); + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this + * object. Since sometimes + * the size of objects can + * not be determined exactly + * (for example: what is the + * memory consumption of an + * STL std::map type with a + * certain number of + * elements?), this is only + * an estimate. however often + * quite close to the true + * value. + */ + std::size_t memory_consumption () const; + + protected: + /** + * This is the abstract function + * through which derived classes + * propagate preprocessed data in + * the form of Patch + * structures (declared in the + * base class DataOutBase) to + * the actual output + * function. You need to overload + * this function to allow the + * output functions to know what + * they shall print. + */ + virtual + const std::vector > & + get_patches () const = 0; + + /** + * Abstract virtual function + * through which the names of + * data sets are obtained by the + * output functions of the base + * class. + */ + virtual + std::vector + get_dataset_names () const = 0; + + /** + * This functions returns + * information about how the + * individual components of + * output files that consist of + * more than one data set are to + * be interpreted. + * + * It returns a list of index + * pairs and corresponding name + * indicating which components of + * the output are to be + * considered vector-valued + * rather than just a collection + * of scalar data. The index + * pairs are inclusive; for + * example, if we have a Stokes + * problem in 2d with components + * (u,v,p), then the + * corresponding vector data + * range should be (0,1), and the + * returned list would consist of + * only a single element with a + * tuple such as (0,1,"velocity"). + * + * Since some of the derived + * classes do not know about + * vector data, this function has + * a default implementation that + * simply returns an empty + * string, meaning that all data + * is to be considered a + * collection of scalar fields. + */ + virtual + std::vector > + get_vector_data_ranges () const; + + /** + * The default number of + * subdivisions for patches. This + * is filled by parse_parameters() + * and should be obeyed by + * build_patches() in derived + * classes. + */ + unsigned int default_subdivisions; - private: - /** - * Standard output format. Use - * this format, if output format - * default_format is - * requested. It can be changed - * by the set_format function - * or in a parameter file. - */ - OutputFormat default_fmt; - - /** - * Flags to be used upon output - * of OpenDX data. Can be changed by - * using the set_flags - * function. - */ - DXFlags dx_flags; - - /** - * Flags to be used upon output - * of UCD data. Can be changed by - * using the set_flags - * function. - */ - UcdFlags ucd_flags; - - /** - * Flags to be used upon output - * of GNUPLOT data. Can be - * changed by using the - * set_flags function. - */ - GnuplotFlags gnuplot_flags; - - /** - * Flags to be used upon output - * of POVRAY data. Can be changed - * by using the set_flags - * function. - */ - PovrayFlags povray_flags; - - /** - * Flags to be used upon output - * of EPS data in one space - * dimension. Can be changed by - * using the set_flags - * function. - */ - EpsFlags eps_flags; - - /** - * Flags to be used upon output - * of gmv data in one space - * dimension. Can be changed by - * using the set_flags - * function. - */ - GmvFlags gmv_flags; - - /** - * Flags to be used upon output - * of Tecplot data in one space - * dimension. Can be changed by - * using the set_flags - * function. - */ - TecplotFlags tecplot_flags; - - /** - * Flags to be used upon output - * of vtk data in one space - * dimension. Can be changed by - * using the set_flags - * function. - */ - VtkFlags vtk_flags; - - /** - * Flags to be used upon output of - * deal.II intermediate data in one space - * dimension. Can be changed by using the - * set_flags function. - */ - Deal_II_IntermediateFlags deal_II_intermediate_flags; + private: + /** + * Standard output format. Use + * this format, if output format + * default_format is + * requested. It can be changed + * by the set_format function + * or in a parameter file. + */ + OutputFormat default_fmt; + + /** + * Flags to be used upon output + * of OpenDX data. Can be changed by + * using the set_flags + * function. + */ + DXFlags dx_flags; + + /** + * Flags to be used upon output + * of UCD data. Can be changed by + * using the set_flags + * function. + */ + UcdFlags ucd_flags; + + /** + * Flags to be used upon output + * of GNUPLOT data. Can be + * changed by using the + * set_flags function. + */ + GnuplotFlags gnuplot_flags; + + /** + * Flags to be used upon output + * of POVRAY data. Can be changed + * by using the set_flags + * function. + */ + PovrayFlags povray_flags; + + /** + * Flags to be used upon output + * of EPS data in one space + * dimension. Can be changed by + * using the set_flags + * function. + */ + EpsFlags eps_flags; + + /** + * Flags to be used upon output + * of gmv data in one space + * dimension. Can be changed by + * using the set_flags + * function. + */ + GmvFlags gmv_flags; + + /** + * Flags to be used upon output + * of Tecplot data in one space + * dimension. Can be changed by + * using the set_flags + * function. + */ + TecplotFlags tecplot_flags; + + /** + * Flags to be used upon output + * of vtk data in one space + * dimension. Can be changed by + * using the set_flags + * function. + */ + VtkFlags vtk_flags; + + /** + * Flags to be used upon output of + * deal.II intermediate data in one space + * dimension. Can be changed by using the + * set_flags function. + */ + Deal_II_IntermediateFlags deal_II_intermediate_flags; }; diff --cc deal.II/include/deal.II/base/exceptions.h index 87c3c8fc09,002ddae58f..ed9b1615e1 --- a/deal.II/include/deal.II/base/exceptions.h +++ b/deal.II/include/deal.II/base/exceptions.h @@@ -47,131 -47,131 +47,131 @@@ DEAL_II_NAMESPACE_OPE */ class ExceptionBase : public std::exception { - public: - /** - * Default constructor. - */ - ExceptionBase (); - - /** - * The constructor takes the file in which the - * error happened, the line and the violated - * condition as well as the name of the - * exception class as a char* as arguments. - */ - ExceptionBase (const char* f, const int l, const char *func, - const char* c, const char *e); - - /** - * Copy constructor. - */ - ExceptionBase (const ExceptionBase &exc); - - /** - * Destructor. Empty, but needed - * for the sake of exception - * specification, since the base - * class has this exception - * specification and the - * automatically generated - * destructor would have a - * different one due to member - * objects. - */ - virtual ~ExceptionBase () throw(); - - /** - * Set the file name and line of where the - * exception appeared as well as the violated - * condition and the name of the exception as - * a char pointer. - */ - void set_fields (const char *f, - const int l, - const char *func, - const char *c, - const char *e); - - /** - * Print out the general part of the error - * information. - */ - void print_exc_data (std::ostream &out) const; - - /** - * Print more specific information about the - * exception which occured. Overload this - * function in your own exception classes. - */ - virtual void print_info (std::ostream &out) const; - - - /** - * Function derived from the base class - * which allows to pass information like - * the line and name of the file where the - * exception occurred as well as user - * information. - * - * This function is mainly used - * when using exceptions - * declared by the - * DeclException* - * macros with the - * throw mechanism or - * the AssertThrow - * macro. - */ - virtual const char * what () const throw (); - - /** - * Print a stacktrace, if one has - * been recorded previously, to - * the given stream. - */ - void print_stack_trace (std::ostream &out) const; - - protected: - /** - * Name of the file this exception happen in. - */ - const char *file; - - /** - * Line number in this file. - */ - unsigned int line; - - /** - * Name of the function, pretty printed. - */ - const char *function; - - /** - * The violated condition, as a string. - */ - const char *cond; - - /** - * Name of the exception and call sequence. - */ - const char *exc; - - /** - * A backtrace to the position - * where the problem happened, if - * the system supports this. - */ - char ** stacktrace; - - /** - * The number of stacktrace - * frames that are stored in the - * previous variable. Zero if the - * system does not support stack - * traces. - */ - int n_stacktrace_frames; + public: + /** + * Default constructor. + */ + ExceptionBase (); + + /** + * The constructor takes the file in which the + * error happened, the line and the violated + * condition as well as the name of the + * exception class as a char* as arguments. + */ + ExceptionBase (const char *f, const int l, const char *func, + const char *c, const char *e); + + /** + * Copy constructor. + */ + ExceptionBase (const ExceptionBase &exc); + + /** + * Destructor. Empty, but needed + * for the sake of exception + * specification, since the base + * class has this exception + * specification and the + * automatically generated + * destructor would have a + * different one due to member + * objects. + */ + virtual ~ExceptionBase () throw(); + + /** + * Set the file name and line of where the + * exception appeared as well as the violated + * condition and the name of the exception as + * a char pointer. + */ + void set_fields (const char *f, + const int l, + const char *func, + const char *c, + const char *e); + + /** + * Print out the general part of the error + * information. + */ + void print_exc_data (std::ostream &out) const; + + /** + * Print more specific information about the + * exception which occured. Overload this + * function in your own exception classes. + */ + virtual void print_info (std::ostream &out) const; + + + /** + * Function derived from the base class + * which allows to pass information like + * the line and name of the file where the + * exception occurred as well as user + * information. + * + * This function is mainly used + * when using exceptions + * declared by the + * DeclException* + * macros with the + * throw mechanism or + * the AssertThrow + * macro. + */ + virtual const char *what () const throw (); + + /** + * Print a stacktrace, if one has + * been recorded previously, to + * the given stream. + */ + void print_stack_trace (std::ostream &out) const; + + protected: + /** + * Name of the file this exception happen in. + */ - const char *file; ++ const char *file; + + /** + * Line number in this file. + */ + unsigned int line; + + /** + * Name of the function, pretty printed. + */ - const char *function; ++ const char *function; + + /** + * The violated condition, as a string. + */ - const char *cond; ++ const char *cond; + + /** + * Name of the exception and call sequence. + */ - const char *exc; ++ const char *exc; + + /** + * A backtrace to the position + * where the problem happened, if + * the system supports this. + */ + char **stacktrace; + + /** + * The number of stacktrace + * frames that are stored in the + * previous variable. Zero if the + * system does not support stack + * traces. + */ + int n_stacktrace_frames; }; diff --cc deal.II/include/deal.II/base/index_set.h index b1de672aed,f8666206b6..3731e02d56 --- a/deal.II/include/deal.II/base/index_set.h +++ b/deal.II/include/deal.II/base/index_set.h @@@ -53,354 -53,441 +53,441 @@@ DEAL_II_NAMESPACE_OPE */ class IndexSet { - public: - /** - * Default constructor. - */ - IndexSet (); - - /** - * Constructor that also sets the - * overall size of the index - * range. - */ - explicit IndexSet (const types::global_dof_index size); - - /** - * Remove all indices from this - * index set. The index set retains - * its size, however. - */ - void clear (); - - /** - * Set the maximal size of the - * indices upon which this object - * operates. - * - * This function can only be - * called if the index set does - * not yet contain any elements. - * This can be achieved by calling - * clear(), for example. - */ - void set_size (const unsigned int size); - - /** - * Return the size of the index - * space of which this index set - * is a subset of. - * - * Note that the result is not equal to - * the number of indices within this - * set. The latter information is - * returned by n_elements(). - */ - unsigned int size () const; - - /** - * Add the half-open range - * $[\text{begin},\text{end})$ to - * the set of indices represented - * by this class. - */ - void add_range (const unsigned int begin, - const unsigned int end); - - /** - * Add an individual index to the - * set of indices. - */ - void add_index (const unsigned int index); - - /** - * Add a whole set of indices - * described by dereferencing - * every element of the the - * iterator range - * [begin,end). - */ - template - void add_indices (const ForwardIterator &begin, - const ForwardIterator &end); - - /** - * Add the given IndexSet @p other to the - * current one, constructing the union of - * *this and @p other. - */ - void add_indices(const IndexSet & other); - - /** - * Return whether the specified - * index is an element of the - * index set. - */ - bool is_element (const types::global_dof_index index) const; - - /** - * Return whether the index set - * stored by this object defines - * a contiguous range. This is - * true also if no indices are - * stored at all. - */ - bool is_contiguous () const; - - /** - * Return the number of elements - * stored in this index set. - */ - types::global_dof_index n_elements () const; - - /** - * Return the global index of the local - * index with number @p local_index - * stored in this index set. @p - * local_index obviously needs to be less - * than n_elements(). - */ - types::global_dof_index nth_index_in_set (const unsigned int local_index) const; - - /** - * Return the how-manyth element of this - * set (counted in ascending order) @p - * global_index is. @p global_index needs - * to be less than the size(). This - * function throws an exception if the - * index @p global_index is not actually - * a member of this index set, i.e. if - * is_element(global_index) is false. - */ - unsigned int index_within_set (const types::global_dof_index global_index) const; - - /** - * Each index set can be - * represented as the union of a - * number of contiguous intervals - * of indices, where if necessary - * intervals may only consist of - * individual elements to - * represent isolated members of - * the index set. - * - * This function returns the - * minimal number of such - * intervals that are needed to - * represent the index set under - * consideration. - */ - unsigned int n_intervals () const; - - /** - * Compress the internal - * representation by merging - * individual elements with - * contiguous ranges, etc. This - * function does not have any - * external effect. - */ - void compress () const; - - /** - * Comparison for equality of - * index sets. This operation is - * only allowed if the size of - * the two sets is the same - * (though of course they do not - * have to have the same number - * of indices). - */ - bool operator == (const IndexSet &is) const; - - /** - * Comparison for inequality of - * index sets. This operation is - * only allowed if the size of - * the two sets is the same - * (though of course they do not - * have to have the same number - * of indices). - */ - bool operator != (const IndexSet &is) const; - - /** - * Return the intersection of the - * current index set and the - * argument given, i.e. a set of - * indices that are elements of - * both index sets. The two index - * sets must have the same size - * (though of course they do not - * have to have the same number - * of indices). - */ - IndexSet operator & (const IndexSet &is) const; - - /** - * This command takes an interval - * [begin, end) and returns - * the intersection of the current - * index set with the interval, shifted - * to the range [0, - * end-begin). - */ - IndexSet get_view (const types::global_dof_index begin, - const types::global_dof_index end) const; - - - /** - * Removes all elements contained in @p - * other from this set. In other words, - * if $x$ is the current object and $o$ - * the argument, then we compute $x - * \leftarrow x \backslash o$. - */ - void subtract_set (const IndexSet & other); - - - /** - * Fills the given vector with all - * indices contained in this IndexSet. - */ - void fill_index_vector(std::vector & indices) const; - - /** - * Fill the given vector with either - * zero or one elements, providing - * a binary representation of this - * index set. The given vector is - * assumed to already have the correct - * size. - * - * The given argument is filled with - * integer values zero and one, using - * vector.operator[]. Thus, - * any object that has such an operator - * can be used as long as it allows - * conversion of integers zero and one to - * elements of the vector. Specifically, - * this is the case for classes Vector, - * BlockVector, but also - * std::vector@, - * std::vector@, and - * std::vector@. - */ - template - void fill_binary_vector (Vector &vector) const; - - /** - * Outputs a text representation of this - * IndexSet to the given stream. Used for - * testing. - */ - template - void print(STREAM &out) const; - - /** - * Writes the IndexSet into a text based - * file format, that can be read in again - * using the read() function. - */ - void write(std::ostream & out) const; - - /** - * Constructs the IndexSet from a text - * based representation given by the - * stream @param in written by the - * write() function. - */ - void read(std::istream & in); - - /** - * Writes the IndexSet into a binary, - * compact representation, that can be - * read in again using the block_read() - * function. - */ - void block_write(std::ostream & out) const; - - /** - * Constructs the IndexSet from a binary - * representation given by the stream - * @param in written by the write_block() - * function. - */ - void block_read(std::istream & in); + public: + /** + * Default constructor. + */ + IndexSet (); + + /** + * Constructor that also sets the + * overall size of the index + * range. + */ - explicit IndexSet (const unsigned int size); ++ explicit IndexSet (const types::global_dof_index size); + + /** + * Remove all indices from this + * index set. The index set retains + * its size, however. + */ + void clear (); + + /** + * Set the maximal size of the + * indices upon which this object + * operates. + * + * This function can only be + * called if the index set does + * not yet contain any elements. + * This can be achieved by calling + * clear(), for example. + */ + void set_size (const unsigned int size); + + /** + * Return the size of the index + * space of which this index set + * is a subset of. + * + * Note that the result is not equal to + * the number of indices within this + * set. The latter information is + * returned by n_elements(). + */ + unsigned int size () const; + + /** + * Add the half-open range + * $[\text{begin},\text{end})$ to + * the set of indices represented + * by this class. + */ + void add_range (const unsigned int begin, + const unsigned int end); + + /** + * Add an individual index to the + * set of indices. + */ + void add_index (const unsigned int index); + + /** + * Add a whole set of indices + * described by dereferencing + * every element of the the + * iterator range + * [begin,end). + */ + template + void add_indices (const ForwardIterator &begin, + const ForwardIterator &end); + + /** + * Add the given IndexSet @p other to the + * current one, constructing the union of + * *this and @p other. + */ + void add_indices(const IndexSet &other); + + /** + * Return whether the specified + * index is an element of the + * index set. + */ - bool is_element (const unsigned int index) const; ++ bool is_element (const types::global_dof_index index) const; + + /** + * Return whether the index set + * stored by this object defines + * a contiguous range. This is + * true also if no indices are + * stored at all. + */ + bool is_contiguous () const; + + /** + * Return the number of elements + * stored in this index set. + */ - unsigned int n_elements () const; ++ types::global_dof_index n_elements () const; + + /** + * Return the global index of the local + * index with number @p local_index + * stored in this index set. @p + * local_index obviously needs to be less + * than n_elements(). + */ - unsigned int nth_index_in_set (const unsigned int local_index) const; ++ types::global_dof_index nth_index_in_set (const unsigned int local_index) const; + + /** + * Return the how-manyth element of this + * set (counted in ascending order) @p + * global_index is. @p global_index needs + * to be less than the size(). This + * function throws an exception if the + * index @p global_index is not actually + * a member of this index set, i.e. if + * is_element(global_index) is false. + */ - unsigned int index_within_set (const unsigned int global_index) const; ++ unsigned int index_within_set (const types::global_dof_index global_index) const; + + /** + * Each index set can be + * represented as the union of a + * number of contiguous intervals + * of indices, where if necessary + * intervals may only consist of + * individual elements to + * represent isolated members of + * the index set. + * + * This function returns the + * minimal number of such + * intervals that are needed to + * represent the index set under + * consideration. + */ + unsigned int n_intervals () const; + + /** + * Compress the internal + * representation by merging + * individual elements with + * contiguous ranges, etc. This + * function does not have any + * external effect. + */ + void compress () const; + + /** + * Comparison for equality of + * index sets. This operation is + * only allowed if the size of + * the two sets is the same + * (though of course they do not + * have to have the same number + * of indices). + */ + bool operator == (const IndexSet &is) const; + + /** + * Comparison for inequality of + * index sets. This operation is + * only allowed if the size of + * the two sets is the same + * (though of course they do not + * have to have the same number + * of indices). + */ + bool operator != (const IndexSet &is) const; + + /** + * Return the intersection of the + * current index set and the + * argument given, i.e. a set of + * indices that are elements of + * both index sets. The two index + * sets must have the same size + * (though of course they do not + * have to have the same number + * of indices). + */ + IndexSet operator & (const IndexSet &is) const; + + /** + * This command takes an interval + * [begin, end) and returns + * the intersection of the current + * index set with the interval, shifted + * to the range [0, + * end-begin). + */ - IndexSet get_view (const unsigned int begin, - const unsigned int end) const; ++ IndexSet get_view (const types::global_dof_index begin, ++ const types::global_dof_index end) const; + + + /** + * Removes all elements contained in @p + * other from this set. In other words, + * if $x$ is the current object and $o$ + * the argument, then we compute $x + * \leftarrow x \backslash o$. + */ + void subtract_set (const IndexSet &other); + + + /** + * Fills the given vector with all + * indices contained in this IndexSet. + */ - void fill_index_vector(std::vector &indices) const; ++ void fill_index_vector(std::vector &indices) const; + + /** + * Fill the given vector with either + * zero or one elements, providing + * a binary representation of this + * index set. The given vector is + * assumed to already have the correct + * size. + * + * The given argument is filled with + * integer values zero and one, using + * vector.operator[]. Thus, + * any object that has such an operator + * can be used as long as it allows + * conversion of integers zero and one to + * elements of the vector. Specifically, + * this is the case for classes Vector, + * BlockVector, but also + * std::vector@, + * std::vector@, and + * std::vector@. + */ + template + void fill_binary_vector (Vector &vector) const; + + /** + * Outputs a text representation of this + * IndexSet to the given stream. Used for + * testing. + */ + template + void print(STREAM &out) const; + + /** + * Writes the IndexSet into a text based + * file format, that can be read in again + * using the read() function. + */ + void write(std::ostream &out) const; + + /** + * Constructs the IndexSet from a text + * based representation given by the + * stream @param in written by the + * write() function. + */ + void read(std::istream &in); + + /** + * Writes the IndexSet into a binary, + * compact representation, that can be + * read in again using the block_read() + * function. + */ + void block_write(std::ostream &out) const; + + /** + * Constructs the IndexSet from a binary + * representation given by the stream + * @param in written by the write_block() + * function. + */ + void block_read(std::istream &in); #ifdef DEAL_II_USE_TRILINOS - /** - * Given an MPI communicator, - * create a Trilinos map object - * that represents a distribution - * of vector elements or matrix - * rows in which we will locally - * store those elements or rows - * for which we store the index - * in the current index set, and - * all the other elements/rows - * elsewhere on one of the other - * MPI processes. - * - * The last argument only plays a - * role if the communicator is a - * parallel one, distributing - * computations across multiple - * processors. In that case, if - * the last argument is false, - * then it is assumed that the - * index sets this function is - * called on on all processors - * are mutually exclusive but - * together enumerate each index - * exactly once. In other words, - * if you call this function on - * two processors, then the index - * sets this function is called - * with must together have all - * possible indices from zero to - * size()-1, and no index must - * appear in both index - * sets. This corresponds, for - * example, to the case where we - * want to split the elements of - * vectors into unique subsets to - * be stored on different - * processors -- no element - * should be owned by more than - * one processor, but each - * element must be owned by one. - * - * On the other hand, if the - * second argument is true, then - * the index sets can be - * overlapping, though they still - * need to contain each index - * exactly once on all processors - * taken together. This is a - * useful operation if we want to - * create vectors that not only - * contain the locally owned - * indices, but for example also - * the elements that correspond - * to degrees of freedom located - * on ghost cells. - */ - Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD, - const bool overlapping = false) const; + /** + * Given an MPI communicator, + * create a Trilinos map object + * that represents a distribution + * of vector elements or matrix + * rows in which we will locally + * store those elements or rows + * for which we store the index + * in the current index set, and + * all the other elements/rows + * elsewhere on one of the other + * MPI processes. + * + * The last argument only plays a + * role if the communicator is a + * parallel one, distributing + * computations across multiple + * processors. In that case, if + * the last argument is false, + * then it is assumed that the + * index sets this function is + * called on on all processors + * are mutually exclusive but + * together enumerate each index + * exactly once. In other words, + * if you call this function on + * two processors, then the index + * sets this function is called + * with must together have all + * possible indices from zero to + * size()-1, and no index must + * appear in both index + * sets. This corresponds, for + * example, to the case where we + * want to split the elements of + * vectors into unique subsets to + * be stored on different + * processors -- no element + * should be owned by more than + * one processor, but each + * element must be owned by one. + * + * On the other hand, if the + * second argument is true, then + * the index sets can be + * overlapping, though they still + * need to contain each index + * exactly once on all processors + * taken together. This is a + * useful operation if we want to + * create vectors that not only + * contain the locally owned + * indices, but for example also + * the elements that correspond + * to degrees of freedom located + * on ghost cells. + */ + Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD, + const bool overlapping = false) const; #endif - /** - * Determine an estimate for the memory - * consumption (in bytes) of this - * object. - */ - std::size_t memory_consumption () const; + /** + * Determine an estimate for the memory + * consumption (in bytes) of this + * object. + */ + std::size_t memory_consumption () const; + + DeclException1 (ExcIndexNotPresent, int, + << "The global index " << arg1 + << " is not an element of this set."); + + /** + * Write or read the data of this object to or + * from a stream for the purpose of serialization + */ + template + void serialize (Archive &ar, const unsigned int version); + + private: + /** + * A type that denotes the half + * open index range + * [begin,end). + * + * The nth_index_in_set denotes + * the how many-th index within + * this IndexSet the first + * element of the current range + * is. This information is only + * accurate if + * IndexSet::compress() has been + * called after the last + * insertion. + */ + struct Range + { - unsigned int begin; - unsigned int end; ++ types::global_dof_index begin; ++ types::global_dof_index end; + - unsigned int nth_index_in_set; ++ types::global_dof_index nth_index_in_set; - DeclException1 (ExcIndexNotPresent, int, - << "The global index " << arg1 - << " is not an element of this set."); + /** + * Default constructor. Since there is no useful choice for + * a default constructed interval, this constructor simply + * creates something that resembles an invalid range. We + * need this constructor for serialization purposes, but the + * invalid range should be filled with something read from + * the archive before it is used, so we should hopefully + * never get to see an invalid range in the wild. + **/ + Range (); + + /** + * Constructor. Create a half-open interval with the given indices. + * + * @param i1 Left end point of the interval. + * @param i2 First index greater than the last index of the indicated range. + **/ - Range (const unsigned int i1, - const unsigned int i2); ++ Range (const types::global_dof_index i1, ++ const types::global_dof_index i2); + + friend + inline bool operator< (const Range &range_1, + const Range &range_2) + { + return ((range_1.begin < range_2.begin) + || + ((range_1.begin == range_2.begin) + && + (range_1.end < range_2.end))); + } + + static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y) + { + return x.end < y.end; + } + + static bool nth_index_compare (const IndexSet::Range &x, + const IndexSet::Range &y) + { + return (x.nth_index_in_set+(x.end-x.begin) < + y.nth_index_in_set+(y.end-y.begin)); + } + + friend + inline bool operator== (const Range &range_1, + const Range &range_2) + { + return ((range_1.begin == range_2.begin) + || + (range_1.begin == range_2.begin)); + } + + std::size_t memory_consumption () const + { + return sizeof(Range); + } /** * Write or read the data of this object to or @@@ -572,11 -572,11 +572,11 @@@ IndexSet::Range::Range ( inline -IndexSet::Range::Range (const unsigned int i1, - const unsigned int i2) +IndexSet::Range::Range (const types::global_dof_index i1, + const types::global_dof_index i2) - : - begin(i1), - end(i2) + : + begin(i1), + end(i2) {} @@@ -591,11 -591,11 +591,11 @@@ IndexSet::IndexSet ( inline -IndexSet::IndexSet (const unsigned int size) +IndexSet::IndexSet (const types::global_dof_index size) - : - is_compressed (true), - index_space_size (size), - largest_range (deal_II_numbers::invalid_unsigned_int) + : + is_compressed (true), + index_space_size (size), + largest_range (deal_II_numbers::invalid_unsigned_int) {} @@@ -706,14 -706,14 +706,14 @@@ voi IndexSet::add_indices (const ForwardIterator &begin, const ForwardIterator &end) { - // insert each element of the - // range. if some of them happen to - // be consecutive, merge them to a - // range + // insert each element of the + // range. if some of them happen to + // be consecutive, merge them to a + // range for (ForwardIterator p=begin; p!=end;) { - const unsigned int begin_index = *p; - unsigned int end_index = begin_index + 1; + const types::global_dof_index begin_index = *p; + types::global_dof_index end_index = begin_index + 1; ForwardIterator q = p; ++q; while ((q != end) && (*q == end_index)) @@@ -828,11 -828,11 +828,11 @@@ IndexSet::is_contiguous () cons inline -unsigned int +types::global_dof_index IndexSet::n_elements () const { - // make sure we have - // non-overlapping ranges + // make sure we have + // non-overlapping ranges compress (); unsigned int v = 0; @@@ -858,10 -858,10 +858,10 @@@ inline unsigned int -IndexSet::nth_index_in_set (const unsigned int n) const +IndexSet::nth_index_in_set (const types::global_dof_index n) const { - // to make this call thread-safe, compress() - // must not be called through this function + // to make this call thread-safe, compress() + // must not be called through this function Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed.")); Assert (n < n_elements(), ExcIndexRange (n, 0, n_elements())); @@@ -908,11 -908,11 +908,11 @@@ inline -unsigned int +types::global_dof_index IndexSet::index_within_set (const unsigned int n) const { - // to make this call thread-safe, compress() - // must not be called through this function + // to make this call thread-safe, compress() + // must not be called through this function Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed.")); Assert (is_element(n) == true, ExcIndexNotPresent (n)); Assert (n < size(), ExcIndexRange (n, 0, size())); diff --cc deal.II/include/deal.II/base/logstream.h index e8223227c5,d4d47ef84b..cb50979e38 --- a/deal.II/include/deal.II/base/logstream.h +++ b/deal.II/include/deal.II/base/logstream.h @@@ -85,530 -85,530 +85,530 @@@ DEAL_II_NAMESPACE_OPE */ class LogStream : public Subscriptor { + public: + /** + * A subclass allowing for the + * safe generation and removal of + * prefices. + * + * Somewhere at the beginning of + * a block, create one of these + * objects, and it will appear as + * a prefix in LogStream output + * like @p deallog. At the end of + * the block, the prefix will + * automatically be removed, when + * this object is destroyed. + */ + class Prefix + { public: - /** - * A subclass allowing for the - * safe generation and removal of - * prefices. - * - * Somewhere at the beginning of - * a block, create one of these - * objects, and it will appear as - * a prefix in LogStream output - * like @p deallog. At the end of - * the block, the prefix will - * automatically be removed, when - * this object is destroyed. - */ - class Prefix - { - public: - /** - * Set a new prefix for - * @p deallog, which will be - * removed when the variable - * is destroyed . - */ - Prefix(const std::string& text); - - /** - * Set a new prefix for the - * given stream, which will - * be removed when the - * variable is destroyed . - */ - Prefix(const std::string& text, LogStream& stream); - - /** - * Remove the prefix - * associated with this - * variable. - */ - ~Prefix (); - - private: - SmartPointer stream; - }; - - /** - * Standard constructor, since we - * intend to provide an object - * deallog in the library. Set the - * standard output stream to std::cerr. - */ - LogStream (); - - /** - * Destructor. - */ - ~LogStream(); - - /** - * Enable output to a second - * stream o. - */ - void attach (std::ostream& o); - - /** - * Disable output to the second - * stream. You may want to call - * close on the stream that was - * previously attached to this object. - */ - void detach (); - - /** - * Setup the logstream for - * regression test mode. - * - * This sets the parameters - * #double_threshold, - * #float_threshold, and #offset - * to nonzero values. The exact - * values being used have been - * determined experimentally and - * can be found in the source - * code. - * - * Called with an argument - * false, switches off - * test mode and sets all - * involved parameters to zero. - */ - void test_mode (bool on=true); - - /** - * Gives the default stream (std_out). - */ - std::ostream& get_console (); - - /** - * Gives the file stream. - */ - std::ostream& get_file_stream (); - - /** - * @return true, if file stream - * has already been attached. - */ - bool has_file () const; - - /** - * Reroutes cerr to LogStream. - * Works as a switch, turning - * logging of cerr on - * and off alternatingly with - * every call. - */ - void log_cerr (); - - /** - * Return the prefix string. - */ - const std::string& get_prefix () const; - - /** - * @deprecated Use Prefix instead - * - * Push another prefix on the - * stack. Prefixes are - * automatically separated by a - * colon and there is a double - * colon after the last prefix. - */ - void push (const std::string& text); - - /** - * @deprecated Use Prefix instead - * - * Remove the last prefix. - */ - void pop (); - - /** - * Maximum number of levels to be - * printed on the console. This - * function allows to restrict - * console output to the upmost - * levels of iterations. Only - * output with less than n - * prefixes is printed. By calling - * this function with n=0, no - * console output will be written. - * - * The previous value of this - * parameter is returned. - */ - unsigned int depth_console (const unsigned int n); - - /** - * Maximum number of levels to be - * written to the log file. The - * functionality is the same as - * depth_console, nevertheless, - * this function should be used - * with care, since it may spoile - * the value of a log file. - * - * The previous value of this - * parameter is returned. - */ - unsigned int depth_file (const unsigned int n); - - /** - * Set time printing flag. If this flag - * is true, each output line will - * be prepended by the user time used - * by the running program so far. - * - * The previous value of this - * parameter is returned. - */ - bool log_execution_time (const bool flag); - - /** - * Output time differences - * between consecutive logs. If - * this function is invoked with - * true, the time difference - * between the previous log line - * and the recent one is - * printed. If it is invoked with - * false, the accumulated - * time since start of the - * program is printed (default - * behavior). - * - * The measurement of times is - * not changed by this function, - * just the output. - * - * The previous value of this - * parameter is returned. - */ - bool log_time_differences (const bool flag); - - /** - * Write detailed timing - * information. - * - * - */ - void timestamp(); - - /** - * Log the thread id. - */ - bool log_thread_id (const bool flag); - - /** - * Set a threshold for the - * minimal absolute value of - * double values. All numbers - * with a smaller absolute value - * will be printed as zero. - * - * The default value for this - * threshold is zero, - * i.e. numbers are printed - * according to their real value. - * - * This feature is mostly useful - * for automated tests: there, - * one would like to reproduce - * the exact same solution in - * each run of a - * testsuite. However, subtle - * difference in processor, - * operating system, or compiler - * version can lead to - * differences in the last few - * digits of numbers, due to - * different rounding. While one - * can avoid trouble for most - * numbers when comparing with - * stored results by simply - * limiting the accuracy of - * output, this does not hold for - * numbers very close to zero, - * i.e. zero plus accumulated - * round-off. For these numbers, - * already the first digit is - * tainted by round-off. Using - * the present function, it is - * possible to eliminate this - * source of problems, by simply - * writing zero to the output in - * this case. - */ - void threshold_double(const double t); - /** - * The same as - * threshold_double(), but for - * float values. - */ - void threshold_float(const float t); - - /** - * Output a constant something - * through this stream. - */ - template - LogStream & operator << (const T &t); - - /** - * Output double precision - * numbers through this - * stream. - * - * If they are set, this function - * applies the methods for making - * floating point output - * reproducible as discussed in - * the introduction. - */ - LogStream & operator << (const double t); - - /** - * Output single precision - * numbers through this - * stream. - * - * If they are set, this function - * applies the methods for making - * floating point output - * reproducible as discussed in - * the introduction. - */ - LogStream & operator << (const float t); - - /** - * Treat ostream - * manipulators. This passes on - * the whole thing to the - * template function with the - * exception of the - * std::endl - * manipulator, for which special - * action is performed: write the - * temporary stream buffer - * including a header to the file - * and std::cout and - * empty the buffer. - * - * An overload of this function is needed - * anyway, since the compiler can't bind - * manipulators like @p std::endl - * directly to template arguments @p T - * like in the previous general - * template. This is due to the fact that - * @p std::endl is actually an overloaded - * set of functions for @p std::ostream, - * @p std::wostream, and potentially more - * of this kind. This function is - * therefore necessary to pick one - * element from this overload set. - */ - LogStream & operator<< (std::ostream& (*p) (std::ostream&)); - - /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this - * object. Since sometimes - * the size of objects can - * not be determined exactly - * (for example: what is the - * memory consumption of an - * STL std::map type with a - * certain number of - * elements?), this is only - * an estimate. however often - * quite close to the true - * value. - */ - std::size_t memory_consumption () const; - - /** - * Exception. - */ - DeclException0(ExcNoFileStreamGiven); + /** + * Set a new prefix for + * @p deallog, which will be + * removed when the variable + * is destroyed . + */ + Prefix(const std::string &text); + + /** + * Set a new prefix for the + * given stream, which will + * be removed when the + * variable is destroyed . + */ + Prefix(const std::string &text, LogStream &stream); + + /** + * Remove the prefix + * associated with this + * variable. + */ + ~Prefix (); private: - - /** - * Stack of strings which are printed - * at the beginning of each line to - * allow identification where the - * output was generated. - */ - std::stack prefixes; - - /** - * Default stream, where the output - * is to go to. This stream defaults - * to std::cerr, but can be set to another - * stream through the constructor. - */ - std::ostream *std_out; - - /** - * Pointer to a stream, where a copy of - * the output is to go to. Usually, this - * will be a file stream. - * - * You can set and reset this stream - * by the attach function. - */ - std::ostream *file; - - /** - * Value denoting the number of - * prefixes to be printed to the - * standard output. If more than - * this number of prefixes is - * pushed to the stack, then no - * output will be generated until - * the number of prefixes shrinks - * back below this number. - */ - unsigned int std_depth; - - /** - * Same for the maximum depth of - * prefixes for output to a file. - */ - unsigned int file_depth; - - /** - * Flag for printing execution time. - */ - bool print_utime; - - /** - * Flag for printing time differences. - */ - bool diff_utime; - - /** - * Time of last output line. - */ - double last_time; - - /** - * Threshold for printing double - * values. Every number with - * absolute value less than this - * is printed as zero. - */ - double double_threshold; - - /** - * Threshold for printing float - * values. Every number with - * absolute value less than this - * is printed as zero. - */ - float float_threshold; - - /** - * An offset added to every float - * or double number upon - * output. This is done after the - * number is compared to - * #double_threshold or #float_threshold, - * but before rounding. - * - * This functionality was - * introduced to produce more - * reproducible floating point - * output for regression - * tests. The rationale is, that - * an exact output value is much - * more likely to be 1/8 than - * 0.124997. If we round to two - * digits though, 1/8 becomes - * unreliably either .12 or .13 - * due to machine accuracy. On - * the other hand, if we add a - * something above machine - * accuracy first, we will always - * get .13. - * - * It is safe to leave this - * value equal to zero. For - * regression tests, the function - * test_mode() sets it to a - * reasonable value. - * - * The offset is relative to the - * magnitude of the number. - */ - double offset; - - /** - * Flag for printing thread id. - */ - bool print_thread_id; - - /** - * The value times() returned - * on initialization. - */ - double reference_time_val; - - /** - * The tms structure times() - * filled on initialization. - */ - struct tms reference_tms; - - /** - * Original buffer of - * std::cerr. We store - * the address of that buffer - * when #log_cerr is called, and - * reset it to this value if - * #log_cerr is called a second - * time, or when the destructor - * of this class is run. - */ - std::streambuf *old_cerr; - - /** - * Print head of line. This prints - * optional time information and - * the contents of the prefix stack. - */ - void print_line_head (); - - /** - * Actually do the work of - * writing output. This function - * unifies the work that is - * common to the two - * operator<< functions. - */ - template - void print (const T &t); - /** - * Check if we are on a new line - * and print the header before - * the data. - */ - std::ostringstream& get_stream(); - - /** - * Type of the stream map - */ - typedef std::map > stream_map_type; - - /** - * We generate a stringstream for - * every process that sends log - * messages. - */ - stream_map_type outstreams; + SmartPointer stream; + }; + + /** + * Standard constructor, since we + * intend to provide an object + * deallog in the library. Set the + * standard output stream to std::cerr. + */ + LogStream (); + + /** + * Destructor. + */ + ~LogStream(); + + /** + * Enable output to a second + * stream o. + */ + void attach (std::ostream &o); + + /** + * Disable output to the second + * stream. You may want to call + * close on the stream that was + * previously attached to this object. + */ + void detach (); + + /** + * Setup the logstream for + * regression test mode. + * + * This sets the parameters + * #double_threshold, + * #float_threshold, and #offset + * to nonzero values. The exact + * values being used have been + * determined experimentally and + * can be found in the source + * code. + * + * Called with an argument + * false, switches off + * test mode and sets all + * involved parameters to zero. + */ + void test_mode (bool on=true); + + /** + * Gives the default stream (std_out). + */ + std::ostream &get_console (); + + /** + * Gives the file stream. + */ + std::ostream &get_file_stream (); + + /** + * @return true, if file stream + * has already been attached. + */ + bool has_file () const; + + /** + * Reroutes cerr to LogStream. + * Works as a switch, turning + * logging of cerr on + * and off alternatingly with + * every call. + */ + void log_cerr (); + + /** + * Return the prefix string. + */ + const std::string &get_prefix () const; + + /** + * @deprecated Use Prefix instead + * + * Push another prefix on the + * stack. Prefixes are + * automatically separated by a + * colon and there is a double + * colon after the last prefix. + */ + void push (const std::string &text); + + /** + * @deprecated Use Prefix instead + * + * Remove the last prefix. + */ + void pop (); + + /** + * Maximum number of levels to be + * printed on the console. This + * function allows to restrict + * console output to the upmost + * levels of iterations. Only + * output with less than n + * prefixes is printed. By calling + * this function with n=0, no + * console output will be written. + * + * The previous value of this + * parameter is returned. + */ + unsigned int depth_console (const unsigned int n); + + /** + * Maximum number of levels to be + * written to the log file. The + * functionality is the same as + * depth_console, nevertheless, + * this function should be used + * with care, since it may spoile + * the value of a log file. + * + * The previous value of this + * parameter is returned. + */ + unsigned int depth_file (const unsigned int n); + + /** + * Set time printing flag. If this flag + * is true, each output line will + * be prepended by the user time used + * by the running program so far. + * + * The previous value of this + * parameter is returned. + */ + bool log_execution_time (const bool flag); + + /** + * Output time differences + * between consecutive logs. If + * this function is invoked with + * true, the time difference + * between the previous log line + * and the recent one is + * printed. If it is invoked with + * false, the accumulated + * time since start of the + * program is printed (default + * behavior). + * + * The measurement of times is + * not changed by this function, + * just the output. + * + * The previous value of this + * parameter is returned. + */ + bool log_time_differences (const bool flag); + + /** + * Write detailed timing + * information. + * + * + */ + void timestamp(); + + /** + * Log the thread id. + */ + bool log_thread_id (const bool flag); + + /** + * Set a threshold for the + * minimal absolute value of + * double values. All numbers + * with a smaller absolute value + * will be printed as zero. + * + * The default value for this + * threshold is zero, + * i.e. numbers are printed + * according to their real value. + * + * This feature is mostly useful + * for automated tests: there, + * one would like to reproduce + * the exact same solution in + * each run of a + * testsuite. However, subtle + * difference in processor, + * operating system, or compiler + * version can lead to + * differences in the last few + * digits of numbers, due to + * different rounding. While one + * can avoid trouble for most + * numbers when comparing with + * stored results by simply + * limiting the accuracy of + * output, this does not hold for + * numbers very close to zero, + * i.e. zero plus accumulated + * round-off. For these numbers, + * already the first digit is + * tainted by round-off. Using + * the present function, it is + * possible to eliminate this + * source of problems, by simply + * writing zero to the output in + * this case. + */ + void threshold_double(const double t); + /** + * The same as + * threshold_double(), but for + * float values. + */ + void threshold_float(const float t); + + /** + * Output a constant something + * through this stream. + */ + template + LogStream &operator << (const T &t); + + /** + * Output double precision + * numbers through this + * stream. + * + * If they are set, this function + * applies the methods for making + * floating point output + * reproducible as discussed in + * the introduction. + */ + LogStream &operator << (const double t); + + /** + * Output single precision + * numbers through this + * stream. + * + * If they are set, this function + * applies the methods for making + * floating point output + * reproducible as discussed in + * the introduction. + */ + LogStream &operator << (const float t); + + /** + * Treat ostream + * manipulators. This passes on + * the whole thing to the + * template function with the + * exception of the + * std::endl + * manipulator, for which special + * action is performed: write the + * temporary stream buffer + * including a header to the file + * and std::cout and + * empty the buffer. + * + * An overload of this function is needed + * anyway, since the compiler can't bind + * manipulators like @p std::endl + * directly to template arguments @p T + * like in the previous general + * template. This is due to the fact that + * @p std::endl is actually an overloaded + * set of functions for @p std::ostream, + * @p std::wostream, and potentially more + * of this kind. This function is + * therefore necessary to pick one + * element from this overload set. + */ + LogStream &operator<< (std::ostream& (*p) (std::ostream &)); + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this + * object. Since sometimes + * the size of objects can + * not be determined exactly + * (for example: what is the + * memory consumption of an + * STL std::map type with a + * certain number of + * elements?), this is only + * an estimate. however often + * quite close to the true + * value. + */ + std::size_t memory_consumption () const; + + /** + * Exception. + */ + DeclException0(ExcNoFileStreamGiven); + + private: + + /** + * Stack of strings which are printed + * at the beginning of each line to + * allow identification where the + * output was generated. + */ + std::stack prefixes; + + /** + * Default stream, where the output + * is to go to. This stream defaults + * to std::cerr, but can be set to another + * stream through the constructor. + */ - std::ostream *std_out; ++ std::ostream *std_out; + + /** + * Pointer to a stream, where a copy of + * the output is to go to. Usually, this + * will be a file stream. + * + * You can set and reset this stream + * by the attach function. + */ - std::ostream *file; ++ std::ostream *file; + + /** + * Value denoting the number of + * prefixes to be printed to the + * standard output. If more than + * this number of prefixes is + * pushed to the stack, then no + * output will be generated until + * the number of prefixes shrinks + * back below this number. + */ + unsigned int std_depth; + + /** + * Same for the maximum depth of + * prefixes for output to a file. + */ + unsigned int file_depth; + + /** + * Flag for printing execution time. + */ + bool print_utime; + + /** + * Flag for printing time differences. + */ + bool diff_utime; + + /** + * Time of last output line. + */ + double last_time; + + /** + * Threshold for printing double + * values. Every number with + * absolute value less than this + * is printed as zero. + */ + double double_threshold; + + /** + * Threshold for printing float + * values. Every number with + * absolute value less than this + * is printed as zero. + */ + float float_threshold; + + /** + * An offset added to every float + * or double number upon + * output. This is done after the + * number is compared to + * #double_threshold or #float_threshold, + * but before rounding. + * + * This functionality was + * introduced to produce more + * reproducible floating point + * output for regression + * tests. The rationale is, that + * an exact output value is much + * more likely to be 1/8 than + * 0.124997. If we round to two + * digits though, 1/8 becomes + * unreliably either .12 or .13 + * due to machine accuracy. On + * the other hand, if we add a + * something above machine + * accuracy first, we will always + * get .13. + * + * It is safe to leave this + * value equal to zero. For + * regression tests, the function + * test_mode() sets it to a + * reasonable value. + * + * The offset is relative to the + * magnitude of the number. + */ + double offset; + + /** + * Flag for printing thread id. + */ + bool print_thread_id; + + /** + * The value times() returned + * on initialization. + */ + double reference_time_val; + + /** + * The tms structure times() + * filled on initialization. + */ + struct tms reference_tms; + + /** + * Original buffer of + * std::cerr. We store + * the address of that buffer + * when #log_cerr is called, and + * reset it to this value if + * #log_cerr is called a second + * time, or when the destructor + * of this class is run. + */ + std::streambuf *old_cerr; + + /** + * Print head of line. This prints + * optional time information and + * the contents of the prefix stack. + */ + void print_line_head (); + + /** + * Actually do the work of + * writing output. This function + * unifies the work that is + * common to the two + * operator<< functions. + */ + template + void print (const T &t); + /** + * Check if we are on a new line + * and print the header before + * the data. + */ + std::ostringstream &get_stream(); + + /** + * Type of the stream map + */ + typedef std::map > stream_map_type; + + /** + * We generate a stringstream for + * every process that sends log + * messages. + */ + stream_map_type outstreams; }; diff --cc deal.II/include/deal.II/base/parallel.h index 285322472f,2036053859..0afd44ae30 --- a/deal.II/include/deal.II/base/parallel.h +++ b/deal.II/include/deal.II/base/parallel.h @@@ -341,13 -341,13 +341,13 @@@ namespace paralle namespace internal { #if DEAL_II_USE_MT == 1 - /** - * Take a range argument and call the - * given function with its begin and end. - */ + /** + * Take a range argument and call the + * given function with its begin and end. + */ template void apply_to_subranges (const tbb::blocked_range &range, - const Function &f) + const Function &f) { f (range.begin(), range.end()); } diff --cc deal.II/include/deal.II/base/parameter_handler.h index 9a81d2a33a,8d706cba21..f19849ef96 --- a/deal.II/include/deal.II/base/parameter_handler.h +++ b/deal.II/include/deal.II/base/parameter_handler.h @@@ -54,611 -54,611 +54,611 @@@ class LogStream namespace Patterns { - /** - * Base class to declare common - * interface. The purpose of this - * class is mostly to define the - * interface of patterns, and to - * force derived classes to have a - * clone function. It is thus, - * in the languages of the "Design - * Patterns" book (Gamma et al.), a - * "prototype". - */ + /** + * Base class to declare common + * interface. The purpose of this + * class is mostly to define the + * interface of patterns, and to + * force derived classes to have a + * clone function. It is thus, + * in the languages of the "Design + * Patterns" book (Gamma et al.), a + * "prototype". + */ class PatternBase { - public: - /** - * Make destructor of this and all - * derived classes virtual. - */ - virtual ~PatternBase (); - - /** - * Return true if the given string - * matches the pattern. - */ - virtual bool match (const std::string &test_string) const = 0; - - /** - * Return a string describing the - * pattern. - */ - virtual std::string description () const = 0; - - /** - * Return a pointer to an - * exact copy of the - * object. This is necessary - * since we want to store - * objects of this type in - * containers, were we need - * to copy objects without - * knowledge of their actual - * data type (we only have - * pointers to the base - * class). - * - * Ownership of the objects - * returned by this function - * is passed to the caller of - * this function. - */ - virtual PatternBase * clone () const = 0; - - /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this object. To - * avoid unnecessary - * overhead, we do not force - * derived classes to provide - * this function as a virtual - * overloaded one, but rather - * try to cast the present - * object to one of the known - * derived classes and if - * that fails then take the - * size of this base class - * instead and add 32 byte - * (this value is arbitrary, - * it should account for - * virtual function tables, - * and some possible data - * elements). Since there are - * usually not many thousands - * of objects of this type - * around, and since the - * memory_consumption - * mechanism is used to find - * out where memory in the - * range of many megabytes - * is, this seems like a - * reasonable approximation. - * - * On the other hand, if you - * know that your class - * deviates from this - * assumption significantly, - * you can still overload - * this function. - */ - virtual std::size_t memory_consumption () const; + public: + /** + * Make destructor of this and all + * derived classes virtual. + */ + virtual ~PatternBase (); + + /** + * Return true if the given string + * matches the pattern. + */ + virtual bool match (const std::string &test_string) const = 0; + + /** + * Return a string describing the + * pattern. + */ + virtual std::string description () const = 0; + + /** + * Return a pointer to an + * exact copy of the + * object. This is necessary + * since we want to store + * objects of this type in + * containers, were we need + * to copy objects without + * knowledge of their actual + * data type (we only have + * pointers to the base + * class). + * + * Ownership of the objects + * returned by this function + * is passed to the caller of + * this function. + */ + virtual PatternBase *clone () const = 0; + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this object. To + * avoid unnecessary + * overhead, we do not force + * derived classes to provide + * this function as a virtual + * overloaded one, but rather + * try to cast the present + * object to one of the known + * derived classes and if + * that fails then take the + * size of this base class + * instead and add 32 byte + * (this value is arbitrary, + * it should account for + * virtual function tables, + * and some possible data + * elements). Since there are + * usually not many thousands + * of objects of this type + * around, and since the + * memory_consumption + * mechanism is used to find + * out where memory in the + * range of many megabytes + * is, this seems like a + * reasonable approximation. + * + * On the other hand, if you + * know that your class + * deviates from this + * assumption significantly, + * you can still overload + * this function. + */ + virtual std::size_t memory_consumption () const; }; - /** - * Returns pointer to the correct - * derived class based on description. - */ - PatternBase * pattern_factory (const std::string& description); - - /** - * Test for the string being an - * integer. If bounds are given - * to the constructor, then the - * integer given also needs to be - * within the interval specified - * by these bounds. Note that - * unlike common convention in - * the C++ standard library, both - * bounds of this interval are - * inclusive; the reason is that - * in practice in most cases, one - * needs closed intervals, but - * these can only be realized - * with inclusive bounds for - * non-integer values. We thus - * stay consistent by always - * using closed intervals. - * - * If the upper bound given to - * the constructor is smaller - * than the lower bound, then the - * infinite interval is implied, - * i.e. every integer is allowed. - * - * Giving bounds may be useful if - * for example a value can only - * be positive and less than a - * reasonable upper bound (for - * example the number of - * refinement steps to be - * performed), or in many other - * cases. - */ + /** + * Returns pointer to the correct + * derived class based on description. + */ + PatternBase *pattern_factory (const std::string &description); + + /** + * Test for the string being an + * integer. If bounds are given + * to the constructor, then the + * integer given also needs to be + * within the interval specified + * by these bounds. Note that + * unlike common convention in + * the C++ standard library, both + * bounds of this interval are + * inclusive; the reason is that + * in practice in most cases, one + * needs closed intervals, but + * these can only be realized + * with inclusive bounds for + * non-integer values. We thus + * stay consistent by always + * using closed intervals. + * + * If the upper bound given to + * the constructor is smaller + * than the lower bound, then the + * infinite interval is implied, + * i.e. every integer is allowed. + * + * Giving bounds may be useful if + * for example a value can only + * be positive and less than a + * reasonable upper bound (for + * example the number of + * refinement steps to be + * performed), or in many other + * cases. + */ class Integer : public PatternBase { - public: - /** - * Minimal integer value. If - * the numeric_limits class - * is available use this - * information to obtain the - * extremal values, otherwise - * set it so that this class - * understands that all values - * are allowed. - */ - static const int min_int_value; - - /** - * Maximal integer value. If - * the numeric_limits class - * is available use this - * information to obtain the - * extremal values, otherwise - * set it so that this class - * understands that all values - * are allowed. - */ - static const int max_int_value; - - /** - * Constructor. Bounds can be - * specified within which a - * valid parameter has to - * be. If the upper bound is - * smaller than the lower - * bound, then the infinite - * interval is meant. The - * default values are chosen - * such that no bounds are - * enforced on parameters. - */ - Integer (const int lower_bound = min_int_value, - const int upper_bound = max_int_value); - - /** - * Return true if the - * string is an integer and - * its value is within the - * specified range. - */ - virtual bool match (const std::string &test_string) const; - - /** - * Return a description of - * the pattern that valid - * strings are expected to - * match. If bounds were - * specified to the - * constructor, then include - * them into this - * description. - */ - virtual std::string description () const; - - /** - * Return a copy of the - * present object, which is - * newly allocated on the - * heap. Ownership of that - * object is transferred to - * the caller of this - * function. - */ - virtual PatternBase * clone () const; - - /** - * Creates new object if the start of - * description matches - * description_init. Ownership of that - * object is transferred to the caller - * of this function. - */ - static Integer* create (const std::string& description); - - private: - /** - * Value of the lower - * bound. A number that - * satisfies the @ref match - * operation of this class - * must be equal to this - * value or larger, if the - * bounds of the interval for - * a valid range. - */ - const int lower_bound; - - /** - * Value of the upper - * bound. A number that - * satisfies the @ref match - * operation of this class - * must be equal to this - * value or less, if the - * bounds of the interval for - * a valid range. - */ - const int upper_bound; - - /** - * Initial part of description - */ - static const char* description_init; + public: + /** + * Minimal integer value. If + * the numeric_limits class + * is available use this + * information to obtain the + * extremal values, otherwise + * set it so that this class + * understands that all values + * are allowed. + */ + static const int min_int_value; + + /** + * Maximal integer value. If + * the numeric_limits class + * is available use this + * information to obtain the + * extremal values, otherwise + * set it so that this class + * understands that all values + * are allowed. + */ + static const int max_int_value; + + /** + * Constructor. Bounds can be + * specified within which a + * valid parameter has to + * be. If the upper bound is + * smaller than the lower + * bound, then the infinite + * interval is meant. The + * default values are chosen + * such that no bounds are + * enforced on parameters. + */ + Integer (const int lower_bound = min_int_value, + const int upper_bound = max_int_value); + + /** + * Return true if the + * string is an integer and + * its value is within the + * specified range. + */ + virtual bool match (const std::string &test_string) const; + + /** + * Return a description of + * the pattern that valid + * strings are expected to + * match. If bounds were + * specified to the + * constructor, then include + * them into this + * description. + */ + virtual std::string description () const; + + /** + * Return a copy of the + * present object, which is + * newly allocated on the + * heap. Ownership of that + * object is transferred to + * the caller of this + * function. + */ + virtual PatternBase *clone () const; + + /** + * Creates new object if the start of + * description matches + * description_init. Ownership of that + * object is transferred to the caller + * of this function. + */ + static Integer *create (const std::string &description); + + private: + /** + * Value of the lower + * bound. A number that + * satisfies the @ref match + * operation of this class + * must be equal to this + * value or larger, if the + * bounds of the interval for + * a valid range. + */ + const int lower_bound; + + /** + * Value of the upper + * bound. A number that + * satisfies the @ref match + * operation of this class + * must be equal to this + * value or less, if the + * bounds of the interval for + * a valid range. + */ + const int upper_bound; + + /** + * Initial part of description + */ + static const char *description_init; }; - /** - * Test for the string being a - * double. If bounds are - * given to the constructor, then - * the integer given also needs - * to be within the interval - * specified by these - * bounds. Note that unlike - * common convention in the C++ - * standard library, both bounds - * of this interval are - * inclusive; the reason is that - * in practice in most cases, one - * needs closed intervals, but - * these can only be realized - * with inclusive bounds for - * non-integer values. We thus - * stay consistent by always - * using closed intervals. - * - * If the upper bound given to - * the constructor is smaller - * than the lower bound, then the - * infinite interval is implied, - * i.e. every integer is allowed. - * - * Giving bounds may be useful if - * for example a value can only - * be positive and less than a - * reasonable upper bound (for - * example damping parameters are - * frequently only reasonable if - * between zero and one), or in - * many other cases. - */ + /** + * Test for the string being a + * double. If bounds are + * given to the constructor, then + * the integer given also needs + * to be within the interval + * specified by these + * bounds. Note that unlike + * common convention in the C++ + * standard library, both bounds + * of this interval are + * inclusive; the reason is that + * in practice in most cases, one + * needs closed intervals, but + * these can only be realized + * with inclusive bounds for + * non-integer values. We thus + * stay consistent by always + * using closed intervals. + * + * If the upper bound given to + * the constructor is smaller + * than the lower bound, then the + * infinite interval is implied, + * i.e. every integer is allowed. + * + * Giving bounds may be useful if + * for example a value can only + * be positive and less than a + * reasonable upper bound (for + * example damping parameters are + * frequently only reasonable if + * between zero and one), or in + * many other cases. + */ class Double : public PatternBase { - public: - /** - * Minimal double value. If the - * std::numeric_limits - * class is available use this - * information to obtain the - * extremal values, otherwise - * set it so that this class - * understands that all values - * are allowed. - */ - static const double min_double_value; - - /** - * Maximal double value. If the - * numeric_limits class is - * available use this - * information to obtain the - * extremal values, otherwise - * set it so that this class - * understands that all values - * are allowed. - */ - static const double max_double_value; - - /** - * Constructor. Bounds can be - * specified within which a - * valid parameter has to - * be. If the upper bound is - * smaller than the lower - * bound, then the infinite - * interval is meant. The - * default values are chosen - * such that no bounds are - * enforced on parameters. - */ - Double (const double lower_bound = min_double_value, - const double upper_bound = max_double_value); - - /** - * Return true if the - * string is a number and its - * value is within the - * specified range. - */ - virtual bool match (const std::string &test_string) const; - - /** - * Return a description of - * the pattern that valid - * strings are expected to - * match. If bounds were - * specified to the - * constructor, then include - * them into this - * description. - */ - virtual std::string description () const; - - /** - * Return a copy of the - * present object, which is - * newly allocated on the - * heap. Ownership of that - * object is transferred to - * the caller of this - * function. - */ - virtual PatternBase * clone () const; - - /** - * Creates new object if the start of - * description matches - * description_init. Ownership of that - * object is transferred to the caller - * of this function. - */ - static Double* create (const std::string& description); - - private: - /** - * Value of the lower - * bound. A number that - * satisfies the @ref match - * operation of this class - * must be equal to this - * value or larger, if the - * bounds of the interval for - * a valid range. - */ - const double lower_bound; - - /** - * Value of the upper - * bound. A number that - * satisfies the @ref match - * operation of this class - * must be equal to this - * value or less, if the - * bounds of the interval for - * a valid range. - */ - const double upper_bound; - - /** - * Initial part of description - */ - static const char* description_init; + public: + /** + * Minimal double value. If the + * std::numeric_limits + * class is available use this + * information to obtain the + * extremal values, otherwise + * set it so that this class + * understands that all values + * are allowed. + */ + static const double min_double_value; + + /** + * Maximal double value. If the + * numeric_limits class is + * available use this + * information to obtain the + * extremal values, otherwise + * set it so that this class + * understands that all values + * are allowed. + */ + static const double max_double_value; + + /** + * Constructor. Bounds can be + * specified within which a + * valid parameter has to + * be. If the upper bound is + * smaller than the lower + * bound, then the infinite + * interval is meant. The + * default values are chosen + * such that no bounds are + * enforced on parameters. + */ + Double (const double lower_bound = min_double_value, + const double upper_bound = max_double_value); + + /** + * Return true if the + * string is a number and its + * value is within the + * specified range. + */ + virtual bool match (const std::string &test_string) const; + + /** + * Return a description of + * the pattern that valid + * strings are expected to + * match. If bounds were + * specified to the + * constructor, then include + * them into this + * description. + */ + virtual std::string description () const; + + /** + * Return a copy of the + * present object, which is + * newly allocated on the + * heap. Ownership of that + * object is transferred to + * the caller of this + * function. + */ + virtual PatternBase *clone () const; + + /** + * Creates new object if the start of + * description matches + * description_init. Ownership of that + * object is transferred to the caller + * of this function. + */ + static Double *create (const std::string &description); + + private: + /** + * Value of the lower + * bound. A number that + * satisfies the @ref match + * operation of this class + * must be equal to this + * value or larger, if the + * bounds of the interval for + * a valid range. + */ + const double lower_bound; + + /** + * Value of the upper + * bound. A number that + * satisfies the @ref match + * operation of this class + * must be equal to this + * value or less, if the + * bounds of the interval for + * a valid range. + */ + const double upper_bound; + + /** + * Initial part of description + */ + static const char *description_init; + }; + + /** + * Test for the string being one + * of a sequence of values given + * like a regular expression. For + * example, if the string given + * to the constructor is + * "red|blue|black", then the + * @ref match function returns + * true exactly if the string + * is either "red" or "blue" or + * "black". Spaces around the + * pipe signs do not matter and + * are eliminated. + */ + class Selection : public PatternBase + { + public: + /** + * Constructor. Take the + * given parameter as the + * specification of valid + * strings. + */ + Selection (const std::string &seq); + + /** + * Return true if the + * string is an element of + * the description list + * passed to the constructor. + */ + virtual bool match (const std::string &test_string) const; + + /** + * Return a description of + * the pattern that valid + * strings are expected to + * match. Here, this is the + * list of valid strings + * passed to the constructor. + */ + virtual std::string description () const; + + /** + * Return a copy of the + * present object, which is + * newly allocated on the + * heap. Ownership of that + * object is transferred to + * the caller of this + * function. + */ + virtual PatternBase *clone () const; + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this object. + */ + std::size_t memory_consumption () const; + + /** + * Creates new object if the start of + * description matches + * description_init. Ownership of that + * object is transferred to the caller + * of this function. + */ + static Selection *create (const std::string &description); + + private: + /** + * List of valid strings as + * passed to the + * constructor. We don't make + * this string constant, as + * we process it somewhat in + * the constructor. + */ + std::string sequence; + + /** + * Initial part of description + */ + static const char *description_init; }; - /** - * Test for the string being one - * of a sequence of values given - * like a regular expression. For - * example, if the string given - * to the constructor is - * "red|blue|black", then the - * @ref match function returns - * true exactly if the string - * is either "red" or "blue" or - * "black". Spaces around the - * pipe signs do not matter and - * are eliminated. - */ - class Selection : public PatternBase - { - public: - /** - * Constructor. Take the - * given parameter as the - * specification of valid - * strings. - */ - Selection (const std::string &seq); - - /** - * Return true if the - * string is an element of - * the description list - * passed to the constructor. - */ - virtual bool match (const std::string &test_string) const; - - /** - * Return a description of - * the pattern that valid - * strings are expected to - * match. Here, this is the - * list of valid strings - * passed to the constructor. - */ - virtual std::string description () const; - - /** - * Return a copy of the - * present object, which is - * newly allocated on the - * heap. Ownership of that - * object is transferred to - * the caller of this - * function. - */ - virtual PatternBase * clone () const; - - /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this object. - */ - std::size_t memory_consumption () const; - - /** - * Creates new object if the start of - * description matches - * description_init. Ownership of that - * object is transferred to the caller - * of this function. - */ - static Selection* create (const std::string& description); - - private: - /** - * List of valid strings as - * passed to the - * constructor. We don't make - * this string constant, as - * we process it somewhat in - * the constructor. - */ - std::string sequence; - - /** - * Initial part of description - */ - static const char* description_init; - }; + /** + * This pattern matches a list of + * comma-separated values each of which + * have to match a pattern given to the + * constructor. With two additional + * parameters, the number of elements this + * list has to have can be specified. If + * none is specified, the list may have + * zero or more entries. + */ + class List : public PatternBase + { + public: + /** + * Maximal integer value. If + * the numeric_limits class + * is available use this + * information to obtain the + * extremal values, otherwise + * set it so that this class + * understands that all values + * are allowed. + */ + static const unsigned int max_int_value; + + /** + * Constructor. Take the + * given parameter as the + * specification of valid + * elements of the list. + * + * The two other arguments can + * be used to denote minimal + * and maximal allowable + * lengths of the list. + */ - List (const PatternBase &base_pattern, ++ List (const PatternBase &base_pattern, + const unsigned int min_elements = 0, + const unsigned int max_elements = max_int_value); + + /** + * Destructor. + */ + virtual ~List (); + + /** + * Return true if the + * string is a comma-separated + * list of strings each of + * which match the pattern + * given to the constructor. + */ + virtual bool match (const std::string &test_string) const; + + /** + * Return a description of + * the pattern that valid + * strings are expected to + * match. + */ + virtual std::string description () const; + + /** + * Return a copy of the + * present object, which is + * newly allocated on the + * heap. Ownership of that + * object is transferred to + * the caller of this + * function. + */ + virtual PatternBase *clone () const; + + /** + * Creates new object if the start of + * description matches + * description_init. Ownership of that + * object is transferred to the caller + * of this function. + */ + static List *create (const std::string &description); + + /** + * Determine an estimate for + * the memory consumption (in + * bytes) of this object. + */ + std::size_t memory_consumption () const; + + /** @addtogroup Exceptions + * @{ */ + + /** + * Exception. + */ + DeclException2 (ExcInvalidRange, + int, int, + << "The values " << arg1 << " and " << arg2 + << " do not form a valid range."); + //@} + private: + /** + * Copy of the pattern that + * each element of the list has + * to satisfy. + */ + PatternBase *pattern; + + /** + * Minimum number of elements + * the list must have. + */ + const unsigned int min_elements; + + /** + * Maximum number of elements + * the list must have. + */ + const unsigned int max_elements; - /** - * This pattern matches a list of - * comma-separated values each of which - * have to match a pattern given to the - * constructor. With two additional - * parameters, the number of elements this - * list has to have can be specified. If - * none is specified, the list may have - * zero or more entries. - */ - class List : public PatternBase - { - public: - /** - * Maximal integer value. If - * the numeric_limits class - * is available use this - * information to obtain the - * extremal values, otherwise - * set it so that this class - * understands that all values - * are allowed. - */ - static const unsigned int max_int_value; - - /** - * Constructor. Take the - * given parameter as the - * specification of valid - * elements of the list. - * - * The two other arguments can - * be used to denote minimal - * and maximal allowable - * lengths of the list. - */ - List (const PatternBase &base_pattern, - const unsigned int min_elements = 0, - const unsigned int max_elements = max_int_value); - - /** - * Destructor. - */ - virtual ~List (); - - /** - * Return true if the - * string is a comma-separated - * list of strings each of - * which match the pattern - * given to the constructor. - */ - virtual bool match (const std::string &test_string) const; - - /** - * Return a description of - * the pattern that valid - * strings are expected to - * match. - */ - virtual std::string description () const; - - /** - * Return a copy of the - * present object, which is - * newly allocated on the - * heap. Ownership of that - * object is transferred to - * the caller of this - * function. - */ - virtual PatternBase * clone () const; - - /** - * Creates new object if the start of - * description matches - * description_init. Ownership of that - * object is transferred to the caller - * of this function. - */ - static List* create (const std::string& description); - - /** - * Determine an estimate for - * the memory consumption (in - * bytes) of this object. - */ - std::size_t memory_consumption () const; - - /** @addtogroup Exceptions - * @{ */ - - /** - * Exception. - */ - DeclException2 (ExcInvalidRange, - int, int, - << "The values " << arg1 << " and " << arg2 - << " do not form a valid range."); - //@} - private: - /** - * Copy of the pattern that - * each element of the list has - * to satisfy. - */ - PatternBase *pattern; - - /** - * Minimum number of elements - * the list must have. - */ - const unsigned int min_elements; - - /** - * Maximum number of elements - * the list must have. - */ - const unsigned int max_elements; - - /** - * Initial part of description - */ - static const char* description_init; + /** + * Initial part of description + */ + static const char *description_init; }; diff --cc deal.II/include/deal.II/base/partitioner.h index 101c00544c,9203409bf6..fff826924b --- a/deal.II/include/deal.II/base/partitioner.h +++ b/deal.II/include/deal.II/base/partitioner.h @@@ -470,7 -470,7 +470,7 @@@ namespace Utilitie inline - const IndexSet& Partitioner::ghost_indices() const - const IndexSet &Partitioner::ghost_indices() const ++ const IndexSet &Partitioner::ghost_indices() const { return ghost_indices_data; } diff --cc deal.II/include/deal.II/base/qprojector.h index 7ffa03220a,e63769f99b..c89c41ffe7 --- a/deal.II/include/deal.II/base/qprojector.h +++ b/deal.II/include/deal.II/base/qprojector.h @@@ -72,380 -72,380 +72,380 @@@ DEAL_II_NAMESPACE_OPE template class QProjector { + public: + /** + * Define a typedef for a + * quadrature that acts on an + * object of one dimension + * less. For cells, this would + * then be a face quadrature. + */ + typedef Quadrature SubQuadrature; + + /** + * Compute the quadrature points + * on the cell if the given + * quadrature formula is used on + * face face_no. For further + * details, see the general doc + * for this class. + */ + static void project_to_face (const SubQuadrature &quadrature, + const unsigned int face_no, + std::vector > &q_points); + + /** + * Compute the cell quadrature + * formula corresponding to using + * quadrature on face + * face_no. For further + * details, see the general doc + * for this class. + */ + static Quadrature + project_to_face (const SubQuadrature &quadrature, + const unsigned int face_no); + + /** + * Compute the quadrature points on the + * cell if the given quadrature formula is + * used on face face_no, subface + * number subface_no corresponding + * to RefineCase::Type + * ref_case. The last argument is + * only used in 3D. + * + * @note Only the points are + * transformed. The quadrature + * weights are the same as those + * of the original rule. + */ + static void project_to_subface (const SubQuadrature &quadrature, + const unsigned int face_no, + const unsigned int subface_no, - std::vector > &q_points, ++ std::vector > &q_points, + const RefinementCase &ref_case=RefinementCase::isotropic_refinement); + + /** + * Compute the cell quadrature formula + * corresponding to using + * quadrature on subface + * subface_no of face + * face_no with + * RefinementCase + * ref_case. The last argument is + * only used in 3D. + * + * @note Only the points are + * transformed. The quadrature + * weights are the same as those + * of the original rule. + */ + static Quadrature + project_to_subface (const SubQuadrature &quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const RefinementCase &ref_case=RefinementCase::isotropic_refinement); + + /** + * Take a face quadrature formula + * and generate a cell quadrature + * formula from it where the + * quadrature points of the given + * argument are projected on all + * faces. + * + * The weights of the new rule + * are replications of the + * original weights. Thus, the + * sum of the weights is not one, + * but the number of faces, which + * is the surface of the + * reference cell. + * + * This in particular allows us + * to extract a subset of points + * corresponding to a single face + * and use it as a quadrature on + * this face, as is done in + * FEFaceValues. + * + * @note In 3D, this function + * produces eight sets of + * quadrature points for each + * face, in order to cope + * possibly different + * orientations of the mesh. + */ + static Quadrature + project_to_all_faces (const SubQuadrature &quadrature); + + /** + * Take a face quadrature formula + * and generate a cell quadrature + * formula from it where the + * quadrature points of the given + * argument are projected on all + * subfaces. + * + * Like in project_to_all_faces(), + * the weights of the new rule + * sum up to the number of faces + * (not subfaces), which + * is the surface of the + * reference cell. + * + * This in particular allows us + * to extract a subset of points + * corresponding to a single subface + * and use it as a quadrature on + * this face, as is done in + * FESubfaceValues. + */ + static Quadrature + project_to_all_subfaces (const SubQuadrature &quadrature); + + /** + * Project a given quadrature + * formula to a child of a + * cell. You may want to use this + * function in case you want to + * extend an integral only over + * the area which a potential + * child would occupy. The child + * numbering is the same as the + * children would be numbered + * upon refinement of the cell. + * + * As integration using this + * quadrature formula now only + * extends over a fraction of the + * cell, the weights of the + * resulting object are divided by + * GeometryInfo::children_per_cell. + */ + static + Quadrature - project_to_child (const Quadrature &quadrature, ++ project_to_child (const Quadrature &quadrature, + const unsigned int child_no); + + /** + * Project a quadrature rule to + * all children of a + * cell. Similarly to + * project_to_all_subfaces(), + * this function replicates the + * formula generated by + * project_to_child() for all + * children, such that the + * weights sum up to one, the + * volume of the total cell + * again. + * + * The child + * numbering is the same as the + * children would be numbered + * upon refinement of the cell. + */ + static + Quadrature - project_to_all_children (const Quadrature &quadrature); ++ project_to_all_children (const Quadrature &quadrature); + + /** + * Project the onedimensional + * rule quadrature to + * the straight line connecting + * the points p1 and + * p2. + */ + static + Quadrature + project_to_line(const Quadrature<1> &quadrature, + const Point &p1, + const Point &p2); + + /** + * Since the + * project_to_all_faces() and + * project_to_all_subfaces() + * functions chain together the + * quadrature points and weights + * of all projections of a face + * quadrature formula to the + * faces or subfaces of a cell, + * we need a way to identify + * where the starting index of + * the points and weights for a + * particular face or subface + * is. This class provides this: + * there are static member + * functions that generate + * objects of this type, given + * face or subface indices, and + * you can then use the generated + * object in place of an integer + * that denotes the offset of a + * given dataset. + * + * @author Wolfgang Bangerth, 2003 + */ + class DataSetDescriptor + { public: - /** - * Define a typedef for a - * quadrature that acts on an - * object of one dimension - * less. For cells, this would - * then be a face quadrature. - */ - typedef Quadrature SubQuadrature; - - /** - * Compute the quadrature points - * on the cell if the given - * quadrature formula is used on - * face face_no. For further - * details, see the general doc - * for this class. - */ - static void project_to_face (const SubQuadrature &quadrature, - const unsigned int face_no, - std::vector > &q_points); - - /** - * Compute the cell quadrature - * formula corresponding to using - * quadrature on face - * face_no. For further - * details, see the general doc - * for this class. - */ - static Quadrature - project_to_face (const SubQuadrature &quadrature, - const unsigned int face_no); - - /** - * Compute the quadrature points on the - * cell if the given quadrature formula is - * used on face face_no, subface - * number subface_no corresponding - * to RefineCase::Type - * ref_case. The last argument is - * only used in 3D. - * - * @note Only the points are - * transformed. The quadrature - * weights are the same as those - * of the original rule. - */ - static void project_to_subface (const SubQuadrature &quadrature, - const unsigned int face_no, - const unsigned int subface_no, - std::vector > &q_points, - const RefinementCase &ref_case=RefinementCase::isotropic_refinement); - - /** - * Compute the cell quadrature formula - * corresponding to using - * quadrature on subface - * subface_no of face - * face_no with - * RefinementCase - * ref_case. The last argument is - * only used in 3D. - * - * @note Only the points are - * transformed. The quadrature - * weights are the same as those - * of the original rule. - */ - static Quadrature - project_to_subface (const SubQuadrature &quadrature, - const unsigned int face_no, - const unsigned int subface_no, - const RefinementCase &ref_case=RefinementCase::isotropic_refinement); - - /** - * Take a face quadrature formula - * and generate a cell quadrature - * formula from it where the - * quadrature points of the given - * argument are projected on all - * faces. - * - * The weights of the new rule - * are replications of the - * original weights. Thus, the - * sum of the weights is not one, - * but the number of faces, which - * is the surface of the - * reference cell. - * - * This in particular allows us - * to extract a subset of points - * corresponding to a single face - * and use it as a quadrature on - * this face, as is done in - * FEFaceValues. - * - * @note In 3D, this function - * produces eight sets of - * quadrature points for each - * face, in order to cope - * possibly different - * orientations of the mesh. - */ - static Quadrature - project_to_all_faces (const SubQuadrature &quadrature); - - /** - * Take a face quadrature formula - * and generate a cell quadrature - * formula from it where the - * quadrature points of the given - * argument are projected on all - * subfaces. - * - * Like in project_to_all_faces(), - * the weights of the new rule - * sum up to the number of faces - * (not subfaces), which - * is the surface of the - * reference cell. - * - * This in particular allows us - * to extract a subset of points - * corresponding to a single subface - * and use it as a quadrature on - * this face, as is done in - * FESubfaceValues. - */ - static Quadrature - project_to_all_subfaces (const SubQuadrature &quadrature); - - /** - * Project a given quadrature - * formula to a child of a - * cell. You may want to use this - * function in case you want to - * extend an integral only over - * the area which a potential - * child would occupy. The child - * numbering is the same as the - * children would be numbered - * upon refinement of the cell. - * - * As integration using this - * quadrature formula now only - * extends over a fraction of the - * cell, the weights of the - * resulting object are divided by - * GeometryInfo::children_per_cell. - */ + /** + * Default constructor. This + * doesn't do much except + * generating an invalid + * index, since you didn't + * give a valid descriptor of + * the cell, face, or subface + * you wanted. + */ + DataSetDescriptor (); + + /** + * Static function to + * generate the offset of a + * cell. Since we only have + * one cell per quadrature + * object, this offset is of + * course zero, but we carry + * this function around for + * consistency with the other + * static functions. + */ + static DataSetDescriptor cell (); + + /** + * Static function to generate an + * offset object for a given face of a + * cell with the given face + * orientation, flip and rotation. This + * function of course is only allowed + * if dim>=2, and the face + * orientation, flip and rotation are + * ignored if the space dimension + * equals 2. + * + * The last argument denotes + * the number of quadrature + * points the + * lower-dimensional face + * quadrature formula (the + * one that has been + * projected onto the faces) + * has. + */ static - Quadrature - project_to_child (const Quadrature &quadrature, - const unsigned int child_no); - - /** - * Project a quadrature rule to - * all children of a - * cell. Similarly to - * project_to_all_subfaces(), - * this function replicates the - * formula generated by - * project_to_child() for all - * children, such that the - * weights sum up to one, the - * volume of the total cell - * again. - * - * The child - * numbering is the same as the - * children would be numbered - * upon refinement of the cell. - */ + DataSetDescriptor + face (const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const unsigned int n_quadrature_points); + + /** + * Static function to generate an + * offset object for a given subface of + * a cell with the given face + * orientation, flip and rotation. This + * function of course is only allowed + * if dim>=2, and the face + * orientation, flip and rotation are + * ignored if the space dimension + * equals 2. + * + * The last but one argument denotes + * the number of quadrature + * points the + * lower-dimensional face + * quadrature formula (the + * one that has been + * projected onto the faces) + * has. + * + * Through the last argument + * anisotropic refinement can be + * respected. + */ static - Quadrature - project_to_all_children (const Quadrature &quadrature); - - /** - * Project the onedimensional - * rule quadrature to - * the straight line connecting - * the points p1 and - * p2. - */ - static - Quadrature - project_to_line(const Quadrature<1>& quadrature, - const Point& p1, - const Point& p2); - - /** - * Since the - * project_to_all_faces() and - * project_to_all_subfaces() - * functions chain together the - * quadrature points and weights - * of all projections of a face - * quadrature formula to the - * faces or subfaces of a cell, - * we need a way to identify - * where the starting index of - * the points and weights for a - * particular face or subface - * is. This class provides this: - * there are static member - * functions that generate - * objects of this type, given - * face or subface indices, and - * you can then use the generated - * object in place of an integer - * that denotes the offset of a - * given dataset. - * - * @author Wolfgang Bangerth, 2003 - */ - class DataSetDescriptor - { - public: - /** - * Default constructor. This - * doesn't do much except - * generating an invalid - * index, since you didn't - * give a valid descriptor of - * the cell, face, or subface - * you wanted. - */ - DataSetDescriptor (); - - /** - * Static function to - * generate the offset of a - * cell. Since we only have - * one cell per quadrature - * object, this offset is of - * course zero, but we carry - * this function around for - * consistency with the other - * static functions. - */ - static DataSetDescriptor cell (); - - /** - * Static function to generate an - * offset object for a given face of a - * cell with the given face - * orientation, flip and rotation. This - * function of course is only allowed - * if dim>=2, and the face - * orientation, flip and rotation are - * ignored if the space dimension - * equals 2. - * - * The last argument denotes - * the number of quadrature - * points the - * lower-dimensional face - * quadrature formula (the - * one that has been - * projected onto the faces) - * has. - */ - static - DataSetDescriptor - face (const unsigned int face_no, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, - const unsigned int n_quadrature_points); - - /** - * Static function to generate an - * offset object for a given subface of - * a cell with the given face - * orientation, flip and rotation. This - * function of course is only allowed - * if dim>=2, and the face - * orientation, flip and rotation are - * ignored if the space dimension - * equals 2. - * - * The last but one argument denotes - * the number of quadrature - * points the - * lower-dimensional face - * quadrature formula (the - * one that has been - * projected onto the faces) - * has. - * - * Through the last argument - * anisotropic refinement can be - * respected. - */ - static - DataSetDescriptor - subface (const unsigned int face_no, - const unsigned int subface_no, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, - const unsigned int n_quadrature_points, - const internal::SubfaceCase ref_case=internal::SubfaceCase::case_isotropic); - - /** - * Conversion operator to an - * integer denoting the - * offset of the first - * element of this dataset in - * the set of quadrature - * formulas all projected - * onto faces and - * subfaces. This conversion - * operator allows us to use - * offset descriptor objects - * in place of integer - * offsets. - */ - operator unsigned int () const; - - private: - /** - * Store the integer offset - * for a given cell, face, or - * subface. - */ - const unsigned int dataset_offset; - - /** - * This is the real - * constructor, but it is - * private and thus only - * available to the static - * member functions above. - */ - DataSetDescriptor (const unsigned int dataset_offset); - }; + DataSetDescriptor + subface (const unsigned int face_no, + const unsigned int subface_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const unsigned int n_quadrature_points, + const internal::SubfaceCase ref_case=internal::SubfaceCase::case_isotropic); + + /** + * Conversion operator to an + * integer denoting the + * offset of the first + * element of this dataset in + * the set of quadrature + * formulas all projected + * onto faces and + * subfaces. This conversion + * operator allows us to use + * offset descriptor objects + * in place of integer + * offsets. + */ + operator unsigned int () const; private: - /** - * Given a quadrature object in - * 2d, reflect all quadrature - * points at the main diagonal - * and return them with their - * original weights. - * - * This function is necessary for - * projecting a 2d quadrature - * rule onto the faces of a 3d - * cube, since there we need both - * orientations. - */ - static Quadrature<2> reflect (const Quadrature<2> &q); - - /** - * Given a quadrature object in - * 2d, rotate all quadrature - * points by @p n_times * 90 degrees - * counterclockwise - * and return them with their - * original weights. - * - * This function is necessary for - * projecting a 2d quadrature - * rule onto the faces of a 3d - * cube, since there we need all - * rotations to account for - * face_flip and face_rotation - * of non-standard faces. - */ - static Quadrature<2> rotate (const Quadrature<2> &q, - const unsigned int n_times); + /** + * Store the integer offset + * for a given cell, face, or + * subface. + */ + const unsigned int dataset_offset; + + /** + * This is the real + * constructor, but it is + * private and thus only + * available to the static + * member functions above. + */ + DataSetDescriptor (const unsigned int dataset_offset); + }; + + private: + /** + * Given a quadrature object in + * 2d, reflect all quadrature + * points at the main diagonal + * and return them with their + * original weights. + * + * This function is necessary for + * projecting a 2d quadrature + * rule onto the faces of a 3d + * cube, since there we need both + * orientations. + */ + static Quadrature<2> reflect (const Quadrature<2> &q); + + /** + * Given a quadrature object in + * 2d, rotate all quadrature + * points by @p n_times * 90 degrees + * counterclockwise + * and return them with their + * original weights. + * + * This function is necessary for + * projecting a 2d quadrature + * rule onto the faces of a 3d + * cube, since there we need all + * rotations to account for + * face_flip and face_rotation + * of non-standard faces. + */ + static Quadrature<2> rotate (const Quadrature<2> &q, + const unsigned int n_times); }; /*@}*/ diff --cc deal.II/include/deal.II/base/smartpointer.h index f580e80acc,afb56e0f90..1efb98f93d --- a/deal.II/include/deal.II/base/smartpointer.h +++ b/deal.II/include/deal.II/base/smartpointer.h @@@ -58,194 -58,194 +58,194 @@@ DEAL_II_NAMESPACE_OPE template class SmartPointer { - public: - /** - * Standard constructor for null - * pointer. The id of this - * pointer is set to the name of - * the class P. - */ - SmartPointer (); - - /* - * Copy constructor for - * SmartPointer. We do now - * copy the object subscribed to - * from tt, but subscribe - * ourselves to it again. - */ - template - SmartPointer (const SmartPointer &tt); - - /* - * Copy constructor for - * SmartPointer. We do now - * copy the object subscribed to - * from tt, but subscribe - * ourselves to it again. - */ - SmartPointer (const SmartPointer &tt); - - /** - * Constructor taking a normal - * pointer. If possible, i.e. if - * the pointer is not a null - * pointer, the constructor - * subscribes to the given object - * to lock it, i.e. to prevent - * its destruction before the end - * of its use. - * - * The id is used in the - * call to - * Subscriptor::subscribe(id) and - * by ~SmartPointer() in the call - * to Subscriptor::unsubscribe(). - */ - SmartPointer (T *t, const char* id); - - /** - * Constructor taking a normal - * pointer. If possible, i.e. if - * the pointer is not a null - * pointer, the constructor - * subscribes to the given object - * to lock it, i.e. to prevent - * its destruction before the end - * of its use. The id of this - * pointer is set to the name of - * the class P. - */ - SmartPointer (T *t); - - - /** - * Destructor, removing the - * subscription. - */ - ~SmartPointer(); - - /** - * Assignment operator for normal - * pointers. The pointer - * subscribes to the new object - * automatically and unsubscribes - * to an old one if it exists. It - * will not try to subscribe to a - * null-pointer, but still - * delete the old subscription. - */ - SmartPointer & operator= (T *tt); - - /** - * Assignment operator for - * SmartPointer. The pointer - * subscribes to the new object - * automatically and unsubscribes - * to an old one if it exists. - */ - template - SmartPointer & operator= (const SmartPointer &tt); - - /** - * Assignment operator for - * SmartPointer. The pointer - * subscribes to the new object - * automatically and unsubscribes - * to an old one if it exists. - */ - SmartPointer & operator= (const SmartPointer &tt); - - /** - * Delete the object pointed to - * and set the pointer to zero. - */ - void clear (); - - /** - * Conversion to normal pointer. - */ - operator T* () const; - - /** - * Dereferencing operator. This - * operator throws an - * ExcNotInitialized if the - * pointer is a null pointer. - */ - T& operator * () const; - - /** - * Dereferencing operator. This - * operator throws an - * ExcNotInitialized if the - * pointer is a null pointer. - */ - T * operator -> () const; - - /** - * Exchange the pointers of this - * object and the argument. Since - * both the objects to which is - * pointed are subscribed to - * before and after, we do not - * have to change their - * subscription counters. - * - * Note that this function (with - * two arguments) and the - * respective functions where one - * of the arguments is a pointer - * and the other one is a C-style - * pointer are implemented in - * global namespace. - */ - template - void swap (SmartPointer &tt); - - /** - * Swap pointers between this - * object and the pointer - * given. As this releases the - * object pointed to presently, - * we reduce its subscription - * count by one, and increase it - * at the object which we will - * point to in the future. - * - * Note that we indeed need a - * reference of a pointer, as we - * want to change the pointer - * variable which we are given. - */ - void swap (T *&tt); - - /** - * Return an estimate of the - * amount of memory (in bytes) - * used by this class. Note in - * particular, that this only - * includes the amount of memory - * used by this object, not - * by the object pointed to. - */ - std::size_t memory_consumption () const; - - private: - /** - * Pointer to the object we want - * to subscribt to. Since it is - * often necessary to follow this - * pointer when debugging, we - * have deliberately chosen a - * short name. - */ - T * t; - /** - * The identification for the - * subscriptor. - */ - const char* const id; + public: + /** + * Standard constructor for null + * pointer. The id of this + * pointer is set to the name of + * the class P. + */ + SmartPointer (); + + /* + * Copy constructor for + * SmartPointer. We do now + * copy the object subscribed to + * from tt, but subscribe + * ourselves to it again. + */ + template + SmartPointer (const SmartPointer &tt); + + /* + * Copy constructor for + * SmartPointer. We do now + * copy the object subscribed to + * from tt, but subscribe + * ourselves to it again. + */ + SmartPointer (const SmartPointer &tt); + + /** + * Constructor taking a normal + * pointer. If possible, i.e. if + * the pointer is not a null + * pointer, the constructor + * subscribes to the given object + * to lock it, i.e. to prevent + * its destruction before the end + * of its use. + * + * The id is used in the + * call to + * Subscriptor::subscribe(id) and + * by ~SmartPointer() in the call + * to Subscriptor::unsubscribe(). + */ + SmartPointer (T *t, const char *id); + + /** + * Constructor taking a normal + * pointer. If possible, i.e. if + * the pointer is not a null + * pointer, the constructor + * subscribes to the given object + * to lock it, i.e. to prevent + * its destruction before the end + * of its use. The id of this + * pointer is set to the name of + * the class P. + */ + SmartPointer (T *t); + + + /** + * Destructor, removing the + * subscription. + */ + ~SmartPointer(); + + /** + * Assignment operator for normal + * pointers. The pointer + * subscribes to the new object + * automatically and unsubscribes + * to an old one if it exists. It + * will not try to subscribe to a + * null-pointer, but still + * delete the old subscription. + */ + SmartPointer &operator= (T *tt); + + /** + * Assignment operator for + * SmartPointer. The pointer + * subscribes to the new object + * automatically and unsubscribes + * to an old one if it exists. + */ + template + SmartPointer &operator= (const SmartPointer &tt); + + /** + * Assignment operator for + * SmartPointer. The pointer + * subscribes to the new object + * automatically and unsubscribes + * to an old one if it exists. + */ + SmartPointer &operator= (const SmartPointer &tt); + + /** + * Delete the object pointed to + * and set the pointer to zero. + */ + void clear (); + + /** + * Conversion to normal pointer. + */ + operator T *() const; + + /** + * Dereferencing operator. This + * operator throws an + * ExcNotInitialized if the + * pointer is a null pointer. + */ + T &operator * () const; + + /** + * Dereferencing operator. This + * operator throws an + * ExcNotInitialized if the + * pointer is a null pointer. + */ + T *operator -> () const; + + /** + * Exchange the pointers of this + * object and the argument. Since + * both the objects to which is + * pointed are subscribed to + * before and after, we do not + * have to change their + * subscription counters. + * + * Note that this function (with + * two arguments) and the + * respective functions where one + * of the arguments is a pointer + * and the other one is a C-style + * pointer are implemented in + * global namespace. + */ + template + void swap (SmartPointer &tt); + + /** + * Swap pointers between this + * object and the pointer + * given. As this releases the + * object pointed to presently, + * we reduce its subscription + * count by one, and increase it + * at the object which we will + * point to in the future. + * + * Note that we indeed need a + * reference of a pointer, as we + * want to change the pointer + * variable which we are given. + */ - void swap (T *&tt); ++ void swap (T *&tt); + + /** + * Return an estimate of the + * amount of memory (in bytes) + * used by this class. Note in + * particular, that this only + * includes the amount of memory + * used by this object, not + * by the object pointed to. + */ + std::size_t memory_consumption () const; + + private: + /** + * Pointer to the object we want + * to subscribt to. Since it is + * often necessary to follow this + * pointer when debugging, we + * have deliberately chosen a + * short name. + */ + T *t; + /** + * The identification for the + * subscriptor. + */ + const char *const id; }; diff --cc deal.II/include/deal.II/base/types.h index 51ac444b5f,e635865640..39ff6bf5c8 --- a/deal.II/include/deal.II/base/types.h +++ b/deal.II/include/deal.II/base/types.h @@@ -22,118 -23,76 +22,118 @@@ DEAL_II_NAMESPACE_OPE */ namespace types { - /** - * The type used to denote - * subdomain_ids of cells. - * - * See the @ref GlossSubdomainId - * "glossary" for more information. - * - * There is a special value, - * numbers::invalid_subdomain_id - * that is used to indicate an - * invalid value of this type. - */ + /** + * The type used to denote + * subdomain_ids of cells. + * + * See the @ref GlossSubdomainId + * "glossary" for more information. + * + * There is a special value, + * numbers::invalid_subdomain_id + * that is used to indicate an + * invalid value of this type. + */ typedef unsigned int subdomain_id; - /** - * @deprecated Old name for the typedef above. - */ + /** + * @deprecated Old name for the typedef above. + */ typedef subdomain_id subdomain_id_t; - /** - * @deprecated Use numbers::invalid_subdomain_id - */ + /** + * @deprecated Use numbers::invalid_subdomain_id + */ const unsigned int invalid_subdomain_id = static_cast(-1); - /** - * @deprecated Use numbers::artificial_subdomain_id - */ + /** + * @deprecated Use numbers::artificial_subdomain_id + */ const unsigned int artificial_subdomain_id = static_cast(-2); +//#define DEAL_II_USE_LARGE_INDEX_TYPE +#ifdef DEAL_II_USE_LARGE_INDEX_TYPE - /** - * The type used for global indices of - * degrees of freedom. While in sequential - * computations the 4 billion indices of - * 32-bit unsigned integers is plenty, - * parallel computations using the - * parallel::distributed::Triangulation - * class can overflow this number and we - * need a bigger index space. - * - * The data type always indicates an - * unsigned integer type. - */ - // TODO: we should check that unsigned long long int - // has the same size as uint64_t ++ /** ++ * The type used for global indices of ++ * degrees of freedom. While in sequential ++ * computations the 4 billion indices of ++ * 32-bit unsigned integers is plenty, ++ * parallel computations using the ++ * parallel::distributed::Triangulation ++ * class can overflow this number and we ++ * need a bigger index space. ++ * ++ * The data type always indicates an ++ * unsigned integer type. ++ */ ++ // TODO: we should check that unsigned long long int ++ // has the same size as uint64_t + typedef unsigned long long int global_dof_index; + - /** - * An identifier that denotes the MPI type - * associated with types::global_dof_index. - */ + /** - * The type used to denote global dof - * indices. ++ * An identifier that denotes the MPI type ++ * associated with types::global_dof_index. ++ */ +# define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED_LONG_LONG +#else - /** - * The type used for global indices of - * degrees of freedom. While in sequential - * computations the 4 billion indices of - * 32-bit unsigned integers is plenty, - * parallel computations using the - * parallel::distributed::Triangulation - * class can overflow this number and we - * need a bigger index space. - * - * The data type always indicates an - * unsigned integer type. - */ ++ /** ++ * The type used for global indices of ++ * degrees of freedom. While in sequential ++ * computations the 4 billion indices of ++ * 32-bit unsigned integers is plenty, ++ * parallel computations using the ++ * parallel::distributed::Triangulation ++ * class can overflow this number and we ++ * need a bigger index space. ++ * ++ * The data type always indicates an ++ * unsigned integer type. + */ typedef unsigned int global_dof_index; - /** - * An identifier that denotes the MPI type - * associated with types::global_dof_index. - */ ++ /** ++ * An identifier that denotes the MPI type ++ * associated with types::global_dof_index. ++ */ +# define DEAL_II_DOF_INDEX_MPI_TYPE MPI_UNSIGNED +#endif + - /** - * @deprecated Use numbers::invalid_dof_index - */ + /** + * @deprecated Use numbers::invalid_dof_index + */ const global_dof_index invalid_dof_index = static_cast(-1); + - /** - * The type used to denote boundary indicators associated with every - * piece of the boundary and, in the case of meshes that describe - * manifolds in higher dimensions, associated with every cell. - * - * There is a special value, numbers::internal_face_boundary_id - * that is used to indicate an invalid value of this type and that - * is used as the boundary indicator for faces that are in the interior - * of the domain and therefore not part of any addressable boundary - * component. - */ + /** + * The type used to denote boundary indicators associated with every + * piece of the boundary and, in the case of meshes that describe + * manifolds in higher dimensions, associated with every cell. + * + * There is a special value, numbers::internal_face_boundary_id + * that is used to indicate an invalid value of this type and that + * is used as the boundary indicator for faces that are in the interior + * of the domain and therefore not part of any addressable boundary + * component. + */ typedef unsigned char boundary_id; - /** - * @deprecated Old name for the typedef above. - */ + /** + * @deprecated Old name for the typedef above. + */ typedef boundary_id boundary_id_t; - /** - * The type used to denote material indicators associated with every - * cell. - * - * There is a special value, numbers::invalid_material_id - * that is used to indicate an invalid value of this type. - */ + /** + * The type used to denote material indicators associated with every + * cell. + * + * There is a special value, numbers::invalid_material_id + * that is used to indicate an invalid value of this type. + */ typedef unsigned char material_id; - /** - * @deprecated Old name for the typedef above. - */ + /** + * @deprecated Old name for the typedef above. + */ typedef material_id material_id_t; } @@@ -162,62 -120,63 +162,62 @@@ namespace number static const unsigned int invalid_unsigned_int = static_cast (-1); - /** - * An invalid value for indices of degrees - * of freedom. - */ + /** + * An invalid value for indices of degrees + * of freedom. + */ const types::global_dof_index invalid_dof_index = static_cast(-1); - /** - * Invalid material_id which we - * need in several places as a - * default value. We assume that - * all material_ids lie in the - * range [0, invalid_material_id). - */ + /** + * Invalid material_id which we + * need in several places as a + * default value. We assume that + * all material_ids lie in the + * range [0, invalid_material_id). + */ const types::material_id invalid_material_id = static_cast(-1); - /** - * The number which we reserve for - * internal faces. We assume that - * all boundary_ids lie in the - * range [0, - * internal_face_boundary_id). - */ + /** + * The number which we reserve for + * internal faces. We assume that + * all boundary_ids lie in the + * range [0, + * internal_face_boundary_id). + */ const types::boundary_id internal_face_boundary_id = static_cast(-1); - /** - * A special id for an invalid - * subdomain id. This value may not - * be used as a valid id but is - * used, for example, for default - * arguments to indicate a - * subdomain id that is not to be - * used. - * - * See the @ref GlossSubdomainId - * "glossary" for more information. - */ + /** + * A special id for an invalid + * subdomain id. This value may not + * be used as a valid id but is + * used, for example, for default + * arguments to indicate a + * subdomain id that is not to be + * used. + * + * See the @ref GlossSubdomainId + * "glossary" for more information. + */ const types::subdomain_id invalid_subdomain_id = static_cast(-1); - /** - * The subdomain id assigned to a - * cell whose true subdomain id we - * don't know, for example because - * it resides on a different - * processor on a mesh that is kept - * distributed on many - * processors. Such cells are - * called "artificial". - * - * See the glossary entries on @ref - * GlossSubdomainId "subdomain ids" - * and @ref GlossArtificialCell - * "artificial cells" as well as - * the @ref distributed module for - * more information. - */ + /** + * The subdomain id assigned to a + * cell whose true subdomain id we + * don't know, for example because + * it resides on a different + * processor on a mesh that is kept + * distributed on many + * processors. Such cells are + * called "artificial". + * + * See the glossary entries on @ref + * GlossSubdomainId "subdomain ids" + * and @ref GlossArtificialCell + * "artificial cells" as well as + * the @ref distributed module for + * more information. + */ const types::subdomain_id artificial_subdomain_id = static_cast(-2); - } diff --cc deal.II/include/deal.II/base/utilities.h index 994630cb33,30df1f3884..4ded7f1412 --- a/deal.II/include/deal.II/base/utilities.h +++ b/deal.II/include/deal.II/base/utilities.h @@@ -265,15 -265,15 +265,15 @@@ namespace Utilitie Iterator lower_bound (Iterator first, Iterator last, - const T &val); + const T &val); - /** - * The same function as above, but taking - * an argument that is used to compare - * individual elements of the sequence of - * objects pointed to by the iterators. - */ + /** + * The same function as above, but taking + * an argument that is used to compare + * individual elements of the sequence of + * objects pointed to by the iterators. + */ template Iterator lower_bound (Iterator first, @@@ -602,43 -602,43 +602,43 @@@ void destroy_communicator (Epetra_Comm &communicator); - /** - * Return the number of MPI processes - * there exist in the given communicator - * object. If this is a sequential job, - * it returns 1. - */ + /** + * Return the number of MPI processes + * there exist in the given communicator + * object. If this is a sequential job, + * it returns 1. + */ unsigned int get_n_mpi_processes (const Epetra_Comm &mpi_communicator); - /** - * Return the number of the present MPI - * process in the space of processes - * described by the given - * communicator. This will be a unique - * value for each process between zero - * and (less than) the number of all - * processes (given by - * get_n_mpi_processes()). - */ + /** + * Return the number of the present MPI + * process in the space of processes + * described by the given + * communicator. This will be a unique + * value for each process between zero + * and (less than) the number of all + * processes (given by + * get_n_mpi_processes()). + */ unsigned int get_this_mpi_process (const Epetra_Comm &mpi_communicator); - /** - * Given a Trilinos Epetra map, create a - * new map that has the same subdivision - * of elements to processors but uses the - * given communicator object instead of - * the one stored in the first - * argument. In essence, this means that - * we create a map that communicates - * among the same processors in the same - * way, but using a separate channel. - * - * This function is typically used with a - * communicator that has been obtained by - * the duplicate_communicator() function. - */ + /** + * Given a Trilinos Epetra map, create a + * new map that has the same subdivision + * of elements to processors but uses the + * given communicator object instead of + * the one stored in the first + * argument. In essence, this means that + * we create a map that communicates + * among the same processors in the same + * way, but using a separate channel. + * + * This function is typically used with a + * communicator that has been obtained by + * the duplicate_communicator() function. + */ Epetra_Map - duplicate_map (const Epetra_BlockMap &map, + duplicate_map (const Epetra_BlockMap &map, const Epetra_Comm &comm); } diff --cc deal.II/include/deal.II/dofs/block_info.h index a8a874055a,9738685202..27341acc0b --- a/deal.II/include/deal.II/dofs/block_info.h +++ b/deal.II/include/deal.II/dofs/block_info.h @@@ -88,147 -88,147 +88,147 @@@ namespace h */ class BlockInfo : public Subscriptor { - public: - /** - * @brief Fill the object with values - * describing block structure - * of the DoFHandler. - * - * This function will also clear - * the local() indices. - */ - template - void initialize(const DoFHandler&, bool levels_only = false, bool multigrid = false); - - /** - * @brief Fill the object with values - * describing level block - * structure of the - * MGDoFHandler. If - * levels_only is false, - * the other initialize() is - * called as well. - * - * This function will also clear - * the local() indices. - */ - template - void initialize(const MGDoFHandler&, bool levels_only = false); - - /** - * @brief Initialize block structure - * on cells and compute - * renumbering between cell - * dofs and block cell dofs. - */ - template - void initialize_local(const DoFHandler&); - - /** - * Access the BlockIndices - * structure of the global - * system. - */ - const BlockIndices& global() const; - - /** - * Access BlockIndices for the - * local system on a cell. - */ - const BlockIndices& local() const; - - /** - * Access the BlockIndices - * structure of a level in the - * multilevel hierarchy. - */ - const BlockIndices& level(unsigned int level) const; - - /** - * Return the index after local - * renumbering. - * - * The input of this function is - * an index between zero and the - * number of dofs per cell, - * numbered in local block - * ordering, that is first all - * indices of the first system - * block, then all of the second - * block and so forth. The - * function then outputs the index - * in the standard local - * numbering of DoFAccessor. - */ - unsigned int renumber (const unsigned int i) const; - - /** - * The number of base elements. - */ - unsigned int n_base_elements() const; - - /** - * Return the base element of - * this index. - */ - unsigned int base_element (const unsigned int i) const; - - /** - * Write a summary of the block - * structure to the stream. - */ - template - void - print(OS& stream) const; - - /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. - */ - std::size_t memory_consumption () const; - - /** - * Read or write the data of this object to or - * from a stream for the purpose of serialization - */ - template - void serialize (Archive & ar, - const unsigned int version); - - private: - /** - * @brief The block structure - * of the global system. - */ - BlockIndices bi_global; - /** - * @brief The multilevel block structure. - */ - std::vector levels; - - /** - * @brief The block structure - * of the cell systems. - */ - BlockIndices bi_local; - - /** - * The base element associated - * with each block. - */ - std::vector base_elements; - - /** - * A vector containing the - * renumbering from the - * standard order of degrees of - * freedom on a cell to a - * component wise - * ordering. Filled by - * initialize(). - */ - std::vector local_renumbering; + public: + /** + * @brief Fill the object with values + * describing block structure + * of the DoFHandler. + * + * This function will also clear + * the local() indices. + */ + template - void initialize(const DoFHandler &); ++ void initialize(const DoFHandler &, bool levels_only = false, bool multigrid = false); + + /** + * @brief Fill the object with values + * describing level block + * structure of the + * MGDoFHandler. If + * levels_only is false, + * the other initialize() is + * called as well. + * + * This function will also clear + * the local() indices. + */ + template + void initialize(const MGDoFHandler &, bool levels_only = false); + + /** + * @brief Initialize block structure + * on cells and compute + * renumbering between cell + * dofs and block cell dofs. + */ + template + void initialize_local(const DoFHandler &); + + /** + * Access the BlockIndices + * structure of the global + * system. + */ + const BlockIndices &global() const; + + /** + * Access BlockIndices for the + * local system on a cell. + */ + const BlockIndices &local() const; + + /** + * Access the BlockIndices + * structure of a level in the + * multilevel hierarchy. + */ + const BlockIndices &level(unsigned int level) const; + + /** + * Return the index after local + * renumbering. + * + * The input of this function is + * an index between zero and the + * number of dofs per cell, + * numbered in local block + * ordering, that is first all + * indices of the first system + * block, then all of the second + * block and so forth. The + * function then outputs the index + * in the standard local + * numbering of DoFAccessor. + */ + unsigned int renumber (const unsigned int i) const; + + /** + * The number of base elements. + */ + unsigned int n_base_elements() const; + + /** + * Return the base element of + * this index. + */ + unsigned int base_element (const unsigned int i) const; + + /** + * Write a summary of the block + * structure to the stream. + */ + template + void + print(OS &stream) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + std::size_t memory_consumption () const; + + /** + * Read or write the data of this object to or + * from a stream for the purpose of serialization + */ + template + void serialize (Archive &ar, + const unsigned int version); + + private: + /** + * @brief The block structure + * of the global system. + */ + BlockIndices bi_global; + /** + * @brief The multilevel block structure. + */ + std::vector levels; + + /** + * @brief The block structure + * of the cell systems. + */ + BlockIndices bi_local; + + /** + * The base element associated + * with each block. + */ + std::vector base_elements; + + /** + * A vector containing the + * renumbering from the + * standard order of degrees of + * freedom on a cell to a + * component wise + * ordering. Filled by + * initialize(). + */ + std::vector local_renumbering; }; diff --cc deal.II/include/deal.II/dofs/dof_accessor.h index 13d358b470,a0fd1a8c2c..c87cdde900 --- a/deal.II/include/deal.II/dofs/dof_accessor.h +++ b/deal.II/include/deal.II/dofs/dof_accessor.h @@@ -174,626 -174,616 +174,626 @@@ namespace interna template class DoFAccessor : public dealii::internal::DoFAccessor::Inheritance::BaseClass { - public: - - /** - * A static variable that allows users of - * this class to discover the value of - * the second template argument. - */ - static const unsigned int dimension=DH::dimension; - - /** - * A static variable that allows users of - * this class to discover the value of - * the third template argument. - */ - static const unsigned int space_dimension=DH::space_dimension; - - /** - * Declare a typedef to the base - * class to make accessing some - * of the exception classes - * simpler. - */ - typedef - typename dealii::internal::DoFAccessor::Inheritance::BaseClass - BaseClass; - - /** - * Data type passed by the iterator class. - */ - typedef DH AccessorData; - - /** - * @name Constructors - */ - /** - * @{ - */ - - /** - * Default constructor. Provides - * an accessor that can't be - * used. - */ - DoFAccessor (); - - /** - * Constructor - */ - DoFAccessor (const Triangulation *tria, - const int level, - const int index, - const DH *local_data); - - /** - * Conversion constructor. This - * constructor exists to make certain - * constructs simpler to write in - * dimension independent code. For - * example, it allows assigning a face - * iterator to a line iterator, an - * operation that is useful in 2d but - * doesn't make any sense in 3d. The - * constructor here exists for the - * purpose of making the code conform to - * C++ but it will unconditionally abort; - * in other words, assigning a face - * iterator to a line iterator is better - * put into an if-statement that checks - * that the dimension is two, and assign - * to a quad iterator in 3d (an operator - * that, without this constructor would - * be illegal if we happen to compile for - * 2d). - */ - template - DoFAccessor (const InvalidAccessor &); - - /** - * Another conversion operator - * between objects that don't - * make sense, just like the - * previous one. - */ - template - DoFAccessor (const DoFAccessor &); - - /** - * @} - */ - - /** - * Return a handle on the - * DoFHandler object which we - * are using. - */ - const DH & - get_dof_handler () const; - - /** - * Implement the copy operator needed - * for the iterator classes. - */ - void copy_from (const DoFAccessor &a); - - /** - * Copy operator used by the - * iterator class. Keeps the - * previously set dof handler, - * but sets the object - * coordinates of the TriaAccessor. - */ - void copy_from (const TriaAccessorBase &da); - - /** - * Return an iterator pointing to - * the the parent. - */ - TriaIterator > - parent () const; - - /** - * @name Accessing sub-objects - */ - /** - * @{ - */ - - /** - * Return an iterator pointing to - * the the @p c-th child. - */ - TriaIterator > - child (const unsigned int c) const; - - /** - * Pointer to the @p ith line - * bounding this object. If the - * current object is a line itself, - * then the only valid index is - * @p i equals to zero, and the - * function returns an iterator - * to itself. - */ - typename dealii::internal::DoFHandler::Iterators::line_iterator - line (const unsigned int i) const; - - /** - * Pointer to the @p ith quad - * bounding this object. If the - * current object is a quad itself, - * then the only valid index is - * @p i equals to zero, and the - * function returns an iterator - * to itself. - */ - typename dealii::internal::DoFHandler::Iterators::quad_iterator - quad (const unsigned int i) const; - - /** - * @} - */ - - /** - * @name Accessing the DoF indices of this object - */ - /** - * @{ - */ - - /** - * Return the indices of the dofs of this - * object in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * The vector has to have the - * right size before being passed - * to this function. - * - * This function is most often - * used on active objects (edges, - * faces, cells). It can be used - * on non-active objects as well - * (i.e. objects that have - * children), but only if the - * finite element under - * consideration has degrees of - * freedom exclusively on - * vertices. Otherwise, the - * function doesn't make much - * sense, since for example - * inactive edges do not have - * degrees of freedom associated - * with them at all. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - * - * For cells, there is only a - * single possible finite element - * index (namely the one for that - * cell, returned by - * cell-@>active_fe_index. Consequently, - * the derived DoFCellAccessor - * class has an overloaded - * version of this function that - * calls the present function - * with - * cell-@>active_fe_index - * as last argument. - */ - void get_dof_indices (std::vector &dof_indices, - const unsigned int fe_index = DH::default_fe_index) const; - - void get_mg_dof_indices (const int level, std::vector& dof_indices, const unsigned int fe_index = DH::default_fe_index) const; - - /** - * Global DoF index of the i - * degree associated with the @p vertexth - * vertex of the present cell. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - */ - types::global_dof_index vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const unsigned int fe_index = DH::default_fe_index) const; - - types::global_dof_index mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const unsigned int fe_index = DH::default_fe_index) const; - - /** - * Index of the ith degree - * of freedom of this object. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - * - * @note While the get_dof_indices() - * function returns an array that - * contains the indices of all degrees of - * freedom that somehow live on this - * object (i.e. on the vertices, edges or - * interior of this object), the current - * dof_index() function only considers - * the DoFs that really belong to this - * particular object's interior. In other - * words, as an example, if the current - * object refers to a quad (a cell in 2d, - * a face in 3d) and the finite element - * associated with it is a bilinear one, - * then the get_dof_indices() will return - * an array of size 4 while dof_index() - * will produce an exception because no - * degrees are defined in the interior of - * the face. - */ - types::global_dof_index dof_index (const unsigned int i, - const unsigned int fe_index = DH::default_fe_index) const; - - types::global_dof_index mg_dof_index (const int level, const unsigned int i) const; - - /** - * @} - */ - - /** - * @name Accessing the finite element associated with this object - */ - /** - * @{ - */ - - /** - * Return the number of finite - * elements that are active on a - * given object. - * - * For non-hp DoFHandler objects, - * the answer is of course always - * one. However, for - * hp::DoFHandler objects, this - * isn't the case: If this is a - * cell, the answer is of course - * one. If it is a face, the - * answer may be one or two, - * depending on whether the two - * adjacent cells use the same - * finite element or not. If it - * is an edge in 3d, the possible - * return value may be one or any - * other value larger than that. - */ - unsigned int - n_active_fe_indices () const; - - /** - * Return the @p n-th active fe - * index on this object. For - * cells and all non-hp objects, - * there is only a single active - * fe index, so the argument must - * be equal to zero. For - * lower-dimensional hp objects, - * there are - * n_active_fe_indices() active - * finite elements, and this - * function can be queried for - * their indices. - */ - unsigned int - nth_active_fe_index (const unsigned int n) const; - - /** - * Return true if the finite - * element with given index is - * active on the present - * object. For non-hp DoF - * accessors, this is of course - * the case only if @p fe_index - * equals zero. For cells, it is - * the case if @p fe_index equals - * active_fe_index() of this - * cell. For faces and other - * lower-dimensional objects, - * there may be more than one @p - * fe_index that are active on - * any given object (see - * n_active_fe_indices()). - */ - bool - fe_index_is_active (const unsigned int fe_index) const; - - /** - * Return a reference to the finite - * element used on this object with the - * given @p fe_index. @p fe_index must be - * used on this object, - * i.e. fe_index_is_active(fe_index) - * must return true. - */ - const FiniteElement & - get_fe (const unsigned int fe_index) const; - - /** - * @} - */ - - /** - * Exceptions for child classes - * - * @ingroup Exceptions - */ - DeclException0 (ExcInvalidObject); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcVectorNotEmpty); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcVectorDoesNotMatch); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcMatrixDoesNotMatch); - /** - * A function has been called for - * a cell which should be active, - * but is refined. @ref GlossActive - * - * @ingroup Exceptions - */ - DeclException0 (ExcNotActive); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcCantCompareIterators); - - protected: - - /** - * Store the address of the DoFHandler object - * to be accessed. - */ - DH *dof_handler; - - /** - * Compare for equality. - */ - bool operator == (const DoFAccessor &) const; - - /** - * Compare for inequality. - */ - bool operator != (const DoFAccessor &) const; - - /** - * Reset the DoF handler pointer. - */ - void set_dof_handler (DH *dh); - - /** - * Set the index of the - * ith degree of freedom - * of this object to @p index. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - */ - void set_dof_index (const unsigned int i, - const types::global_dof_index index, - const unsigned int fe_index = DH::default_fe_index) const; - - void set_mg_dof_index (const int level, const unsigned int i, const types::global_dof_index index) const; - - /** - * Set the global index of the i - * degree on the @p vertex-th vertex of - * the present cell to @p index. - * - * The last argument denotes the - * finite element index. For the - * standard ::DoFHandler class, - * this value must be equal to - * its default value since that - * class only supports the same - * finite element on all cells - * anyway. - * - * However, for hp objects - * (i.e. the hp::DoFHandler - * class), different finite - * element objects may be used on - * different cells. On faces - * between two cells, as well as - * vertices, there may therefore - * be two sets of degrees of - * freedom, one for each of the - * finite elements used on the - * adjacent cells. In order to - * specify which set of degrees - * of freedom to work on, the - * last argument is used to - * disambiguate. Finally, if this - * function is called for a cell - * object, there can only be a - * single set of degrees of - * freedom, and fe_index has to - * match the result of - * active_fe_index(). - */ - void set_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const types::global_dof_index index, - const unsigned int fe_index = DH::default_fe_index) const; - - void set_mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const types::global_dof_index index, const unsigned int fe_index = DH::default_fe_index) const; - - /** - * Iterator classes need to be friends - * because they need to access operator== - * and operator!=. - */ - template friend class TriaRawIterator; - - - private: - /** - * Copy operator. This is normally used - * in a context like iterator a,b; - * *a=*b;. Presumably, the intent - * here is to copy the object pointed to - * by @p b to the object pointed to by - * @p a. However, the result of - * dereferencing an iterator is not an - * object but an accessor; consequently, - * this operation is not useful for - * iterators on triangulations. We - * declare this function here private, - * thus it may not be used from outside. - * Furthermore it is not implemented and - * will give a linker error if used - * anyway. - */ - DoFAccessor & - operator = (const DoFAccessor &da); - - /** - * Make the DoFHandler class a friend so - * that it can call the set_xxx() - * functions. - */ - template friend class DoFHandler; - template friend class hp::DoFHandler; - - friend struct dealii::internal::DoFHandler::Policy::Implementation; - friend struct dealii::internal::DoFHandler::Implementation; - friend struct dealii::internal::hp::DoFHandler::Implementation; - friend struct dealii::internal::DoFCellAccessor::Implementation; + public: + + /** + * A static variable that allows users of + * this class to discover the value of + * the second template argument. + */ + static const unsigned int dimension=DH::dimension; + + /** + * A static variable that allows users of + * this class to discover the value of + * the third template argument. + */ + static const unsigned int space_dimension=DH::space_dimension; + + /** + * Declare a typedef to the base + * class to make accessing some + * of the exception classes + * simpler. + */ + typedef + typename dealii::internal::DoFAccessor::Inheritance::BaseClass + BaseClass; + + /** + * Data type passed by the iterator class. + */ + typedef DH AccessorData; + + /** + * @name Constructors + */ + /** + * @{ + */ + + /** + * Default constructor. Provides + * an accessor that can't be + * used. + */ + DoFAccessor (); + + /** + * Constructor + */ + DoFAccessor (const Triangulation *tria, + const int level, + const int index, + const DH *local_data); + + /** + * Conversion constructor. This + * constructor exists to make certain + * constructs simpler to write in + * dimension independent code. For + * example, it allows assigning a face + * iterator to a line iterator, an + * operation that is useful in 2d but + * doesn't make any sense in 3d. The + * constructor here exists for the + * purpose of making the code conform to + * C++ but it will unconditionally abort; + * in other words, assigning a face + * iterator to a line iterator is better + * put into an if-statement that checks + * that the dimension is two, and assign + * to a quad iterator in 3d (an operator + * that, without this constructor would + * be illegal if we happen to compile for + * 2d). + */ + template + DoFAccessor (const InvalidAccessor &); + + /** + * Another conversion operator + * between objects that don't + * make sense, just like the + * previous one. + */ + template + DoFAccessor (const DoFAccessor &); + + /** + * @} + */ + + /** + * Return a handle on the + * DoFHandler object which we + * are using. + */ + const DH & + get_dof_handler () const; + + /** + * Implement the copy operator needed + * for the iterator classes. + */ + void copy_from (const DoFAccessor &a); + + /** + * Copy operator used by the + * iterator class. Keeps the + * previously set dof handler, + * but sets the object + * coordinates of the TriaAccessor. + */ + void copy_from (const TriaAccessorBase &da); + + /** + * Return an iterator pointing to + * the the parent. + */ + TriaIterator > + parent () const; + + /** + * @name Accessing sub-objects + */ + /** + * @{ + */ + + /** + * Return an iterator pointing to + * the the @p c-th child. + */ + TriaIterator > + child (const unsigned int c) const; + + /** + * Pointer to the @p ith line + * bounding this object. If the + * current object is a line itself, + * then the only valid index is + * @p i equals to zero, and the + * function returns an iterator + * to itself. + */ + typename dealii::internal::DoFHandler::Iterators::line_iterator + line (const unsigned int i) const; + + /** + * Pointer to the @p ith quad + * bounding this object. If the + * current object is a quad itself, + * then the only valid index is + * @p i equals to zero, and the + * function returns an iterator + * to itself. + */ + typename dealii::internal::DoFHandler::Iterators::quad_iterator + quad (const unsigned int i) const; + + /** + * @} + */ + + /** + * @name Accessing the DoF indices of this object + */ + /** + * @{ + */ + + /** + * Return the indices of the dofs of this + * object in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * The vector has to have the + * right size before being passed + * to this function. + * + * This function is most often + * used on active objects (edges, + * faces, cells). It can be used + * on non-active objects as well + * (i.e. objects that have + * children), but only if the + * finite element under + * consideration has degrees of + * freedom exclusively on + * vertices. Otherwise, the + * function doesn't make much + * sense, since for example + * inactive edges do not have + * degrees of freedom associated + * with them at all. + * + * The last argument denotes the + * finite element index. For the + * standard ::DoFHandler class, + * this value must be equal to + * its default value since that + * class only supports the same + * finite element on all cells + * anyway. + * + * However, for hp objects + * (i.e. the hp::DoFHandler + * class), different finite + * element objects may be used on + * different cells. On faces + * between two cells, as well as + * vertices, there may therefore + * be two sets of degrees of + * freedom, one for each of the + * finite elements used on the + * adjacent cells. In order to + * specify which set of degrees + * of freedom to work on, the + * last argument is used to + * disambiguate. Finally, if this + * function is called for a cell + * object, there can only be a + * single set of degrees of + * freedom, and fe_index has to + * match the result of + * active_fe_index(). + * + * For cells, there is only a + * single possible finite element + * index (namely the one for that + * cell, returned by + * cell-@>active_fe_index. Consequently, + * the derived DoFCellAccessor + * class has an overloaded + * version of this function that + * calls the present function + * with + * cell-@>active_fe_index + * as last argument. + */ - void get_dof_indices (std::vector &dof_indices, ++ void get_dof_indices (std::vector &dof_indices, + const unsigned int fe_index = DH::default_fe_index) const; + ++ void get_mg_dof_indices (const int level, std::vector &dof_indices, const unsigned int fe_index = DH::default_fe_index) const; ++ + /** + * Global DoF index of the i + * degree associated with the @p vertexth + * vertex of the present cell. + * + * The last argument denotes the + * finite element index. For the + * standard ::DoFHandler class, + * this value must be equal to + * its default value since that + * class only supports the same + * finite element on all cells + * anyway. + * + * However, for hp objects + * (i.e. the hp::DoFHandler + * class), different finite + * element objects may be used on + * different cells. On faces + * between two cells, as well as + * vertices, there may therefore + * be two sets of degrees of + * freedom, one for each of the + * finite elements used on the + * adjacent cells. In order to + * specify which set of degrees + * of freedom to work on, the + * last argument is used to + * disambiguate. Finally, if this + * function is called for a cell + * object, there can only be a + * single set of degrees of + * freedom, and fe_index has to + * match the result of + * active_fe_index(). + */ - unsigned int vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const unsigned int fe_index = DH::default_fe_index) const; ++ types::global_dof_index vertex_dof_index (const unsigned int vertex, ++ const unsigned int i, ++ const unsigned int fe_index = DH::default_fe_index) const; ++ ++ types::global_dof_index mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const unsigned int fe_index = DH::default_fe_index) const; + + /** + * Index of the ith degree + * of freedom of this object. + * + * The last argument denotes the + * finite element index. For the + * standard ::DoFHandler class, + * this value must be equal to + * its default value since that + * class only supports the same + * finite element on all cells + * anyway. + * + * However, for hp objects + * (i.e. the hp::DoFHandler + * class), different finite + * element objects may be used on + * different cells. On faces + * between two cells, as well as + * vertices, there may therefore + * be two sets of degrees of + * freedom, one for each of the + * finite elements used on the + * adjacent cells. In order to + * specify which set of degrees + * of freedom to work on, the + * last argument is used to + * disambiguate. Finally, if this + * function is called for a cell + * object, there can only be a + * single set of degrees of + * freedom, and fe_index has to + * match the result of + * active_fe_index(). + * + * @note While the get_dof_indices() + * function returns an array that + * contains the indices of all degrees of + * freedom that somehow live on this + * object (i.e. on the vertices, edges or + * interior of this object), the current + * dof_index() function only considers + * the DoFs that really belong to this + * particular object's interior. In other + * words, as an example, if the current + * object refers to a quad (a cell in 2d, + * a face in 3d) and the finite element + * associated with it is a bilinear one, + * then the get_dof_indices() will return + * an array of size 4 while dof_index() + * will produce an exception because no + * degrees are defined in the interior of + * the face. + */ - unsigned int dof_index (const unsigned int i, - const unsigned int fe_index = DH::default_fe_index) const; ++ types::global_dof_index dof_index (const unsigned int i, ++ const unsigned int fe_index = DH::default_fe_index) const; ++ ++ types::global_dof_index mg_dof_index (const int level, const unsigned int i) const; + + /** + * @} + */ + + /** + * @name Accessing the finite element associated with this object + */ + /** + * @{ + */ + + /** + * Return the number of finite + * elements that are active on a + * given object. + * + * For non-hp DoFHandler objects, + * the answer is of course always + * one. However, for + * hp::DoFHandler objects, this + * isn't the case: If this is a + * cell, the answer is of course + * one. If it is a face, the + * answer may be one or two, + * depending on whether the two + * adjacent cells use the same + * finite element or not. If it + * is an edge in 3d, the possible + * return value may be one or any + * other value larger than that. + */ + unsigned int + n_active_fe_indices () const; + + /** + * Return the @p n-th active fe + * index on this object. For + * cells and all non-hp objects, + * there is only a single active + * fe index, so the argument must + * be equal to zero. For + * lower-dimensional hp objects, + * there are + * n_active_fe_indices() active + * finite elements, and this + * function can be queried for + * their indices. + */ + unsigned int + nth_active_fe_index (const unsigned int n) const; + + /** + * Return true if the finite + * element with given index is + * active on the present + * object. For non-hp DoF + * accessors, this is of course + * the case only if @p fe_index + * equals zero. For cells, it is + * the case if @p fe_index equals + * active_fe_index() of this + * cell. For faces and other + * lower-dimensional objects, + * there may be more than one @p + * fe_index that are active on + * any given object (see + * n_active_fe_indices()). + */ + bool + fe_index_is_active (const unsigned int fe_index) const; + + /** + * Return a reference to the finite + * element used on this object with the + * given @p fe_index. @p fe_index must be + * used on this object, + * i.e. fe_index_is_active(fe_index) + * must return true. + */ + const FiniteElement & + get_fe (const unsigned int fe_index) const; + + /** + * @} + */ + + /** + * Exceptions for child classes + * + * @ingroup Exceptions + */ + DeclException0 (ExcInvalidObject); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcVectorNotEmpty); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcVectorDoesNotMatch); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcMatrixDoesNotMatch); + /** + * A function has been called for + * a cell which should be active, + * but is refined. @ref GlossActive + * + * @ingroup Exceptions + */ + DeclException0 (ExcNotActive); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcCantCompareIterators); + + protected: + + /** + * Store the address of the DoFHandler object + * to be accessed. + */ + DH *dof_handler; + + /** + * Compare for equality. + */ + bool operator == (const DoFAccessor &) const; + + /** + * Compare for inequality. + */ + bool operator != (const DoFAccessor &) const; + + /** + * Reset the DoF handler pointer. + */ + void set_dof_handler (DH *dh); + + /** + * Set the index of the + * ith degree of freedom + * of this object to @p index. + * + * The last argument denotes the + * finite element index. For the + * standard ::DoFHandler class, + * this value must be equal to + * its default value since that + * class only supports the same + * finite element on all cells + * anyway. + * + * However, for hp objects + * (i.e. the hp::DoFHandler + * class), different finite + * element objects may be used on + * different cells. On faces + * between two cells, as well as + * vertices, there may therefore + * be two sets of degrees of + * freedom, one for each of the + * finite elements used on the + * adjacent cells. In order to + * specify which set of degrees + * of freedom to work on, the + * last argument is used to + * disambiguate. Finally, if this + * function is called for a cell + * object, there can only be a + * single set of degrees of + * freedom, and fe_index has to + * match the result of + * active_fe_index(). + */ + void set_dof_index (const unsigned int i, - const unsigned int index, ++ const types::global_dof_index index, + const unsigned int fe_index = DH::default_fe_index) const; + ++ void set_mg_dof_index (const int level, const unsigned int i, const types::global_dof_index index) const; ++ + /** + * Set the global index of the i + * degree on the @p vertex-th vertex of + * the present cell to @p index. + * + * The last argument denotes the + * finite element index. For the + * standard ::DoFHandler class, + * this value must be equal to + * its default value since that + * class only supports the same + * finite element on all cells + * anyway. + * + * However, for hp objects + * (i.e. the hp::DoFHandler + * class), different finite + * element objects may be used on + * different cells. On faces + * between two cells, as well as + * vertices, there may therefore + * be two sets of degrees of + * freedom, one for each of the + * finite elements used on the + * adjacent cells. In order to + * specify which set of degrees + * of freedom to work on, the + * last argument is used to + * disambiguate. Finally, if this + * function is called for a cell + * object, there can only be a + * single set of degrees of + * freedom, and fe_index has to + * match the result of + * active_fe_index(). + */ + void set_vertex_dof_index (const unsigned int vertex, + const unsigned int i, - const unsigned int index, ++ const types::global_dof_index index, + const unsigned int fe_index = DH::default_fe_index) const; + ++ void set_mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const types::global_dof_index index, const unsigned int fe_index = DH::default_fe_index) const; ++ + /** + * Iterator classes need to be friends + * because they need to access operator== + * and operator!=. + */ + template friend class TriaRawIterator; + + + private: + /** + * Copy operator. This is normally used + * in a context like iterator a,b; + * *a=*b;. Presumably, the intent + * here is to copy the object pointed to + * by @p b to the object pointed to by + * @p a. However, the result of + * dereferencing an iterator is not an + * object but an accessor; consequently, + * this operation is not useful for + * iterators on triangulations. We + * declare this function here private, + * thus it may not be used from outside. + * Furthermore it is not implemented and + * will give a linker error if used + * anyway. + */ + DoFAccessor & + operator = (const DoFAccessor &da); + + /** + * Make the DoFHandler class a friend so + * that it can call the set_xxx() + * functions. + */ + template friend class DoFHandler; + template friend class hp::DoFHandler; + + friend struct dealii::internal::DoFHandler::Policy::Implementation; + friend struct dealii::internal::DoFHandler::Implementation; + friend struct dealii::internal::hp::DoFHandler::Implementation; + friend struct dealii::internal::DoFCellAccessor::Implementation; }; @@@ -811,642 -801,642 +811,642 @@@ template