From 503c05b4c2c63f5a74a5903ec80c57c431adf4e9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 24 Apr 2002 15:58:48 +0000 Subject: [PATCH] More doc. git-svn-id: https://svn.dealii.org/trunk@5730 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-14/step-14.cc | 285 ++++++++++++++++++++++------ 1 file changed, 227 insertions(+), 58 deletions(-) diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index 1e420deb3c..b82d9ef7f4 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -12,6 +12,7 @@ /* further information on this license. */ + // Start out with well known things... #include #include #include @@ -350,6 +351,18 @@ namespace Evaluation // variable: point_derivative += solution_gradients[q_point][0]; ++evaluation_point_hits; + + // Finally break out of + // the innermost loop + // iterating over the + // vertices of the + // present cell, since if + // we have found the + // evaluation point at + // one vertex it cannot + // be at a following + // vertex as well: + break; }; // Now we have looped over all @@ -764,7 +777,7 @@ namespace LaplaceSolver void Solver::LinearSystem::solve (Vector &solution) const { - SolverControl solver_control (10000, 1e-12); //TODO! + SolverControl solver_control (5000, 1e-12); PrimitiveVectorMemory<> vector_memory; SolverCG<> cg (solver_control, vector_memory); @@ -779,14 +792,16 @@ namespace LaplaceSolver - // @sect{The PrimalSolver class} + // @sect4{The PrimalSolver class} // The ``PrimalSolver'' class is // also mostly unchanged except for // overloading the functions // ``solve_problem'', ``n_dofs'', // and ``postprocess'' of the base - // class. These overloaded + // class, and implementing the + // ``output_solution'' + // function. These overloaded // functions do nothing particular // besides calling the functions of // the base class -- that seems @@ -817,9 +832,24 @@ namespace LaplaceSolver // derived class using these // additional functions. // - // Except for the reimplementation - // of these three functions, this - // class is also unchanged. + // Regarding the implementation of + // the ``output_solution'' + // function, we keep the + // ``GlobalRefinement'' and + // ``RefinementKelly'' classes in + // this program, and they can then + // rely on the default + // implementation of this function + // which simply outputs the primal + // solution. The class implementing + // dual weighted error estimators + // will overload this function + // itself, to also output the dual + // solution. + // + // Except for this, the class is + // unchanged with respect to the + // previous example. template class PrimalSolver : public Solver { @@ -832,16 +862,16 @@ namespace LaplaceSolver const Function &boundary_values); virtual - void - solve_problem (); + void solve_problem (); virtual - unsigned int - n_dofs () const; + unsigned int n_dofs () const; virtual - void - postprocess (const Evaluation::EvaluationBase &postprocessor) const; + void postprocess (const Evaluation::EvaluationBase &postprocessor) const; + + virtual + void output_solution () const; protected: const SmartPointer > rhs_function; @@ -890,6 +920,34 @@ namespace LaplaceSolver { return Solver::postprocess(postprocessor); }; + + + template + void + PrimalSolver::output_solution () const + { + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream filename; +#else + std::ostrstream filename; +#endif + filename << "solution-" + << refinement_cycle + << ".gnuplot" + << std::ends; +#ifdef HAVE_STD_STRINGSTREAM + std::ofstream out (filename.str().c_str()); +#else + std::ofstream out (filename.str()); +#endif + + data_out.write (out, DataOut::gnuplot); + }; @@ -940,7 +998,12 @@ namespace LaplaceSolver }; - //TODO!! + // @sect4{The RefinementGlobal and RefinementKelly classes} + + // For the following two classes, + // the same applies as for most of + // the above: the class is taken + // from the previous example as-is: template class RefinementGlobal : public PrimalSolver { @@ -1300,7 +1363,6 @@ namespace Data q += sin(10*p(i)+5*p(0)*p(0)); const double exponential = exp(q); return exponential; -// return 0; // TODO! }; @@ -1328,8 +1390,6 @@ namespace Data t1 = t1*t1; return -u*(t1+t2+t3); -// const double pi = 3.1415926536; -// return 2.*pi*pi*sin(pi*p(0))*sin(pi*p(1)); //TODO!! }; @@ -1647,7 +1707,7 @@ namespace Data // with the rest of the program. - + //TODO namespace DualFunctional { template @@ -1715,13 +1775,120 @@ namespace DualFunctional AssertThrow (evaluation_point_found, ExcEvaluationPointNotFound(evaluation_point)); }; + + + + template + class PointXDerivativeEvaluation : public DualFunctionalBase + { + public: + PointXDerivativeEvaluation (const Point &evaluation_point, + const double tolerance); + + virtual + void + assemble_rhs (const DoFHandler &dof_handler, + Vector &rhs) const; + DeclException1 (ExcEvaluationPointNotFound, + Point, + << "The evaluation point " << arg1 + << " was not found among the vertices of the present grid."); + + protected: + const Point evaluation_point; + const double tolerance; + }; + + + template + PointXDerivativeEvaluation:: + PointXDerivativeEvaluation (const Point &evaluation_point, + const double tolerance) + : + evaluation_point (evaluation_point), + tolerance (tolerance) + {}; + + + template + void + PointXDerivativeEvaluation:: + assemble_rhs (const DoFHandler &dof_handler, + Vector &rhs) const + { + rhs.reinit (dof_handler.n_dofs()); + + QTrapez<1> q_trapez; + QIterated quadrature (q_trapez, 4); + FEValues fe_values (dof_handler.get_fe(), quadrature, + update_gradients | + update_q_points | + update_JxW_values); + const unsigned int n_q_points = fe_values.n_quadrature_points; + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + double total_volume = 0; + + for (; cell!=endc; ++cell) + if (cell->center().distance(evaluation_point) - + cell->diameter()/2 + < + tolerance) + { + fe_values.reinit (cell); + cell_rhs.clear (); + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + for (unsigned int i=0; i class DualSolver : public Solver { @@ -1732,7 +1899,6 @@ namespace LaplaceSolver const Quadrature &face_quadrature, const DualFunctional::DualFunctionalBase &dual_functional); - //TODO!! virtual void solve_problem (); @@ -1807,6 +1973,9 @@ namespace LaplaceSolver }; + // @sect4{The WeightedResidual class} + + //TODO! template class WeightedResidual : public PrimalSolver, public DualSolver @@ -1915,6 +2084,20 @@ namespace LaplaceSolver * memory itself, or synchronise * with other threads. */ + struct CellData + { + FEValues fe_values; + const SmartPointer > right_hand_side; + + std::vector cell_residual; + std::vector rhs_values; + std::vector dual_weights; + typename std::vector > cell_grad_grads; + CellData (const FiniteElement &dof_handler, + const Quadrature &quadrature, + const Function &right_hand_side); + }; + struct FaceData { FEFaceValues fe_face_values_cell; @@ -1929,19 +2112,6 @@ namespace LaplaceSolver const Quadrature &face_quadrature); }; - struct CellData - { - FEValues fe_values; - const SmartPointer > right_hand_side; - - std::vector cell_residual; - std::vector rhs_values; - std::vector dual_weights; - typename std::vector > cell_grad_grads; - CellData (const FiniteElement &dof_handler, - const Quadrature &quadrature, - const Function &right_hand_side); - }; @@ -2014,6 +2184,30 @@ namespace LaplaceSolver + template + WeightedResidual::CellData:: + CellData (const FiniteElement &fe, + const Quadrature &quadrature, + const Function &right_hand_side) + : + fe_values (fe, quadrature, + update_values | + update_second_derivatives | + update_q_points | + update_JxW_values), + right_hand_side (&right_hand_side) + { + const unsigned int n_q_points + = quadrature.n_quadrature_points; + + cell_residual.resize(n_q_points); + rhs_values.resize(n_q_points); + dual_weights.resize(n_q_points); + cell_grad_grads.resize(n_q_points); + }; + + + template WeightedResidual::FaceData:: FaceData (const FiniteElement &fe, @@ -2043,30 +2237,6 @@ namespace LaplaceSolver - template - WeightedResidual::CellData:: - CellData (const FiniteElement &fe, - const Quadrature &quadrature, - const Function &right_hand_side) - : - fe_values (fe, quadrature, - update_values | - update_second_derivatives | - update_q_points | - update_JxW_values), - right_hand_side (&right_hand_side) - { - const unsigned int n_q_points - = quadrature.n_quadrature_points; - - cell_residual.resize(n_q_points); - rhs_values.resize(n_q_points); - dual_weights.resize(n_q_points); - cell_grad_grads.resize(n_q_points); - }; - - - template WeightedResidual:: @@ -2943,8 +3113,8 @@ void solve_problem () data->create_coarse_grid (triangulation); const Point evaluation_point(0.75,0.75); - const DualFunctional::PointValueEvaluation - dual_functional (evaluation_point); + const DualFunctional::PointXDerivativeEvaluation + dual_functional (evaluation_point, 0.01); LaplaceSolver::Base * solver = 0; solver = new LaplaceSolver::WeightedResidual (triangulation, @@ -2986,7 +3156,6 @@ int main () deallog.depth_console (0); solve_problem<2> (); -// solve_problem<2> ("kelly"); } catch (std::exception &exc) { -- 2.39.5