From d22c8dbe74c14ba6695f03f82ad1881c8073ec76 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 24 Apr 2002 08:00:35 +0000 Subject: [PATCH] More text. Grid output. git-svn-id: https://svn.dealii.org/trunk@5722 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-14/step-14.cc | 619 ++++++++++++++++++++++------ 1 file changed, 488 insertions(+), 131 deletions(-) diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index fbeaa43fb7..9f9105be06 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -1,5 +1,3 @@ -// TODO: bv for primal problem wrong! -// check Galerkin orthogonality /* $Id$ */ /* Author: Wolfgang Bangerth, ETH Zurich, 2002 */ @@ -27,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -155,6 +154,50 @@ namespace Evaluation << 1.594915543-point_value << std::endl; }; + + template + class GridOutput : public EvaluationBase + { + 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; + }; + + + template + GridOutput:: + GridOutput (const std::string &output_name_base) + : + output_name_base (output_name_base) + {}; + + + template + void + GridOutput::operator () (const DoFHandler &dof_handler, + const Vector &/*solution*/) const + { +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream filename; +#else + std::ostrstream filename; +#endif + filename << output_name_base << "-" + << refinement_cycle + << ".eps" + << std::ends; +#ifdef HAVE_STD_STRINGSTREAM + std::ofstream out (filename.str().c_str()); +#else + std::ofstream out (filename.str()); +#endif + + GridOut().write_eps (dof_handler.get_tria(), out); + }; }; @@ -685,52 +728,75 @@ namespace LaplaceSolver // with the same data sets as in the // previous one, but as it may so // happen that someone wants to run - // the program with a different - // solution and right hand side - // function, we show a simple - // technique to do exactly that. For - // more clarity, we furthermore pack - // everything that has to do with - // equation data into a namespace of - // its own. + // the program with different + // boundary values and right hand side + // functions, or on a different grid, + // we show a simple technique to do + // exactly that. For more clarity, we + // furthermore pack everything that + // has to do with equation data into + // a namespace of its own. // - // Basically, the idea is as follows: - // let us have a structure for each - // set of data, in which we have two + // The underlying assumption is that + // this is a research program, and + // that there we often have a number + // of test cases that consist of a + // domain, a right hand side, + // boundary values, possibly a + // specified coefficient, and a + // number of other parameters. They + // often vary all at the same time + // when shifting from one example to + // another. To make handling such + // sets of problem description + // parameters simple is the goal of + // the following. + // + // Basically, the idea is this: let + // us have a structure for each set + // of data, in which we pack + // everything that describes a test + // case: here, these are two // subclasses, one called - // ``Solution'' for the exact - // solution (and also used as right - // hand side), and one called - // ``RightHandSide''. Since the - // solution of the previous example - // program looked like curved ridges, - // we use this name here for the - // enclosing class. Note that the - // names of the two inner class have - // to be the same for all enclosing - // classes, and also that we have - // attached the dimension template - // argument to the enclosing class - // rather than to the inner ones, to - // make further processing simpler. - // (From a language viewpoint, a - // namespace would be better to - // encapsulate these inner classes, - // rather than a structure. However, - // namespaces cannot be given as - // template arguments, so we use a - // structure to allow a second object - // to select from within its given + // ``BoundaryValues'' for the + // boundary values of the exact + // solution, and one called + // ``RightHandSide'', and then a way + // to generate the coarse grid. Since + // the solution of the previous + // example program looked like curved + // ridges, we use this name here for + // the enclosing class. Note that the + // names of the two inner classes + // have to be the same for all + // enclosing test case classes, and + // also that we have attached the + // dimension template argument to the + // enclosing class rather than to the + // inner ones, to make further + // processing simpler. (From a + // language viewpoint, a namespace + // would be better to encapsulate + // these inner classes, rather than a + // structure. However, namespaces + // cannot be given as template + // arguments, so we use a structure + // to allow a second object to select + // from within its given // argument. The enclosing structure, - // of course, has no members apart - // from the classes it declares, and - // will never be instantiated.) + // of course, has no member variables + // apart from the classes it + // declares, and a static function to + // generate the coarse mesh; it will + // in general never be instantiated.) // - // The idea is then the following: we - // can generate objects for - // solution/boundary values and right - // hand side by simply giving the - // name of the outer class as a + // The idea is then the following + // (this is the right time to also + // take a brief look at the code + // below): we can generate objects + // for boundary values and + // right hand side by simply giving + // the name of the outer class as a // template argument to a class which // we call here ``Data::SetUp'', and // it then creates objects for the @@ -741,17 +807,17 @@ namespace LaplaceSolver // ``Data::SetUp'', // and everything we need to know // about the solution would be static - // member variables of that object. + // member variables and functions of + // that object. // // This approach might seem like // overkill in this case, but will // become very handy once a certain // set up is not only characterized - // by a solution (or Dirichlet - // boundary values) and a right hand - // side function, but in addition by - // material properties, Neumann - // values, different boundary + // by Dirichlet boundary values and a + // right hand side function, but in + // addition by material properties, + // Neumann values, different boundary // descriptors, etc. In that case, // the ``SetUp'' class might consist // of a dozen or more objects, and @@ -799,10 +865,13 @@ namespace Data struct SetUpBase { virtual - const Function & get_solution () const = 0; + const Function & get_boundary_values () const = 0; virtual const Function & get_right_hand_side () const = 0; + + virtual + void create_coarse_grid (Triangulation &coarse_grid) const = 0; }; @@ -825,31 +894,35 @@ namespace Data SetUp () {}; virtual - const Function & get_solution () const; + const Function & get_boundary_values () const; virtual const Function & get_right_hand_side () const; + + virtual + void create_coarse_grid (Triangulation &coarse_grid) const; + private: - static const typename Traits::Solution solution; - static const typename Traits::RightHandSide right_hand_side; + static const typename Traits::BoundaryValues boundary_values; + static const typename Traits::RightHandSide right_hand_side; }; // We have to provide definitions // for the static member variables // of the above class: template - const typename Traits::Solution SetUp::solution; + const typename Traits::BoundaryValues SetUp::boundary_values; template - const typename Traits::RightHandSide SetUp::right_hand_side; + const typename Traits::RightHandSide SetUp::right_hand_side; // And definitions of the member // functions: template const Function & - SetUp::get_solution () const + SetUp::get_boundary_values () const { - return solution; + return boundary_values; }; @@ -859,21 +932,32 @@ namespace Data { return right_hand_side; }; + + + template + void + SetUp:: + create_coarse_grid (Triangulation &coarse_grid) const + { + Traits::create_coarse_grid (coarse_grid); + }; // @sect4{The CurvedRidges class} // The class that is used to - // describe the solution and right - // hand side of the ``curved - // ridge'' problem is like so: + // describe the boundary values and + // right hand side of the ``curved + // ridge'' problem already used in + // the step-13 example program is + // then like so: template struct CurvedRidges { - class Solution : public Function + class BoundaryValues : public Function { public: - Solution () : Function () {}; + BoundaryValues () : Function () {}; virtual double value (const Point &p, const unsigned int component) const; @@ -888,13 +972,18 @@ namespace Data virtual double value (const Point &p, const unsigned int component) const; }; + + static + void + create_coarse_grid (Triangulation &coarse_grid); }; template double - CurvedRidges::Solution::value (const Point &p, - const unsigned int /*component*/) const + CurvedRidges::BoundaryValues:: + value (const Point &p, + const unsigned int /*component*/) const { double q = p(0); for (unsigned int i=1; i + void + CurvedRidges:: + create_coarse_grid (Triangulation &coarse_grid) + { + GridGenerator::hyper_cube (coarse_grid, -1, 1); + coarse_grid.refine_global (2); + }; + + // @sect4{The Exercise_2_3 class} + + // This example program was written + // while giving practical courses + // for a lecture on adaptive finite + // element methods and duality + // based error estimates. For these + // courses, we had one exercise, + // which required to solve the + // Laplace equation with constant + // right hand side on a square + // domain with a square hole in the + // center, and zero boundary + // values. Since the implementation + // of the properties of this + // problem is so particularly + // simple here, lets do it. As the + // number of the exercise was 2.3, + // we take the liberty to retain + // this name for the class as well. + template + struct Exercise_2_3 + { + // We need a class to denote + // the boundary values of the + // problem. In this case, this + // is simple: it's the zero + // function, so don't even + // declare a class, just a + // typedef: + typedef ZeroFunction BoundaryValues; + + // Second, a class that denotes + // the right hand side. Since + // they are constant, just + // subclass the corresponding + // class of the library and be + // done: + class RightHandSide : public ConstantFunction + { + public: + RightHandSide () : ConstantFunction (1.) {}; + }; + + // Finally a function to + // generate the coarse + // grid. This is somewhat more + // complicated here, see + // immediately below. + static + void + create_coarse_grid (Triangulation &coarse_grid); + }; + // As stated above, the grid for + // this example is the square + // [-1,1]^2 with the square + // [-1/2,1/2]^2 as hole in it. We + // create the coarse grid as 3 + // times 3 cells with the middle + // one missing. + // + // Of course, the example has an + // extension to 3d, but since this + // function cannot be written in a + // dimension independent way we + // choose not to implement this + // here, but rather only specialize + // the template for dim=2. If you + // compile the program for 3d, + // you'll get a message from the + // linker that this function is not + // implemented for 3d, and needs to + // be provided. + // + // For the creation of this + // geometry, the library has no + // predefined method. In this case, + // the geometry is still simple + // enough to do the creation by + // hand, rather than using a mesh + // generator. + template <> + void + Exercise_2_3<2>:: + create_coarse_grid (Triangulation<2> &coarse_grid) + { + // First define the space + // dimension, to allow those + // parts of the function that are + // actually dimension independent + // to use this variable. That + // makes it simpler if you later + // takes this as a starting point + // to implement the 3d version. + const unsigned int dim = 2; + + // Then have a list of + // vertices. Here, they are 24 (5 + // times 5, with the middle one + // omitted). It is probably best + // to draw a sketch here. Note + // that we leave the number of + // vertices open at first, but + // then let the compiler compute + // this number afterwards. This + // reduces the possibility of + // having the dimension to large + // and leaving the last ones + // uninitialized. + static const Point<2> vertices_1[] + = { Point<2> (-1., -1.), + Point<2> (-1./2, -1.), + Point<2> (0., -1.), + Point<2> (+1./2, -1.), + Point<2> (+1, -1.), + + Point<2> (-1., -1./2.), + Point<2> (-1./2, -1./2.), + Point<2> (0., -1./2.), + Point<2> (+1./2, -1./2.), + Point<2> (+1, -1./2.), + + Point<2> (-1., 0.), + Point<2> (-1./2, 0.), + Point<2> (+1./2, 0.), + Point<2> (+1, 0.), + + Point<2> (-1., 1./2.), + Point<2> (-1./2, 1./2.), + Point<2> (0., 1./2.), + Point<2> (+1./2, 1./2.), + Point<2> (+1, 1./2.), + + Point<2> (-1., 1.), + Point<2> (-1./2, 1.), + Point<2> (0., 1.), + Point<2> (+1./2, 1.), + Point<2> (+1, 1.) }; + const unsigned int + n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]); + + // From this static list of + // vertices, we generate an STL + // vector of the vertices, as + // this is the data type the + // library wants to see. + const std::vector > vertices (&vertices_1[0], + &vertices_1[n_vertices]); + + // Next, we have to define the + // cells and the vertices they + // contain. Here, we have 8 + // vertices, but leave the number + // open and let it be computed + // afterwards: + static const int cell_vertices[][GeometryInfo::vertices_per_cell] + = {{0, 1, 6,5}, + {1, 2, 7, 6}, + {2, 3, 8, 7}, + {3, 4, 9, 8}, + {5, 6, 11, 10}, + {8, 9, 13, 12}, + {10, 11, 15, 14}, + {12, 13, 18, 17}, + {14, 15, 20, 19}, + {15, 16, 21, 20}, + {16, 17, 22, 21}, + {17, 18, 23, 22}}; + const unsigned int + n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]); + + // Again, we generate a C++ + // vector type from this, but + // this time by looping over the + // cells (yes, this is + // boring). Additionally, we set + // the material indicator to zero + // for all the cells: + std::vector > cells (n_cells, CellData()); + for (unsigned int i=0; i::vertices_per_cell; + ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; + + // Finally pass all this + // information to the library to + // generate a triangulation. The + // last parameter may be used to + // pass information about + // non-zero boundary indicators + // at certain faces of the + // triangulation to the library, + // but we don't want that here, + // so we give an empty object: + coarse_grid.create_triangulation (vertices, + cells, + SubCellData()); + + // And since we want that the + // evaluation point (3/4,3/4) in + // this example is a grid point, + // we refine once globally: + coarse_grid.refine_global (1); + }; +}; + + // @sect4{Discussion} + // + // As you have now read through this + // framework, you may be wondering + // why we have not chosen to + // implement the classes implementing + // a certain setup (like the + // ``CurvedRidges'' class) directly + // as classes derived from + // ``Data::SetUpBase''. Indeed, we + // could have done very well so. The + // only reason is that then we would + // have to have member variables for + // the solution and right hand side + // classes in the ``CurvedRidges'' + // class, as well as member functions + // overloading the abstract functions + // of the base class giving access to + // these member variables. The + // ``SetUp'' class has the sole + // reason to relieve us from the need + // to reiterate these member + // variables and functions that would + // be necessary in all such + // classes. In some way, the template + // mechanism here only provides a way + // to have default implementations + // for a number of functions that + // depend on external quantities and + // can thus not be provided using + // normal virtual functions, at least + // not without the help of templates. + // + // However, there might be good + // reasons to actually implement + // classes derived from + // ``Data::SetUpBase'', for example + // if the solution or right hand side + // classes require constructors that + // take arguments, which the + // ``Data::SetUpBase'' class cannot + // provide. In that case, subclassing + // is a worthwhile strategy. Other + // possibilities for special cases + // are to derive from + // ``Data::SetUp'' where + // ``SomeSetUp'' denotes a class, or + // even to explicitly specialize + // ``Data::SetUp''. The + // latter allows to transparently use + // the way the ``SetUp'' class is + // used for other set-ups, but with + // special actions taken for special + // arguments. + // + // A final observation favoring the + // approach taken here is the + // following: we have found numerous + // times that when starting a + // project, the number of parameters + // (usually boundary values, right + // hand side, coarse grid, just as + // here) was small, and the number of + // test cases was small as well. One + // then starts out by handcoding them + // into a number of ``switch'' + // statements. Over time, projects + // grow, and so does the number of + // test cases. The number of + // ``switch'' statements grows with + // that, and their length as well, + // and one starts to find ways to + // consider impossible examples where + // domains, boundary values, and + // right hand sides do not fit + // together any more, and starts + // loosing the overview over the + // whole structure. Encapsulating + // everything belonging to a certain + // test case into a structure of its + // own has proven worthwhile for + // this, as it keeps everything that + // belongs to one test case in one + // place. Furthermore, it allows to + // put these things all in one or + // more files that are only devoted + // to test cases and their data, + // without having to bring their + // actual implementation into contact + // with the rest of the program. + namespace DualFunctional @@ -1503,17 +1900,25 @@ namespace LaplaceSolver // primal finite element // space. Fortunately, the // library provides functions for - // these two actions. + // these two actions. (In + // general, for transformations + // between different finite + // elements, the ``FETools'' + // namespace provides a number of + // functions.) Vector primal_solution (DualSolver::dof_handler.n_dofs()); FETools::interpolate (PrimalSolver::dof_handler, PrimalSolver::solution, DualSolver::dof_handler, primal_solution); Vector dual_weights (DualSolver::dof_handler.n_dofs()); - FETools::interpolation_difference (DualSolver::dof_handler, - DualSolver::solution, - *PrimalSolver::fe, - dual_weights); +// FETools::interpolation_difference (DualSolver::dof_handler, +// DualSolver::solution, +// *PrimalSolver::fe, +// dual_weights); + dual_weights = DualSolver::solution; + abort (); // check Galerkin orthogonality, also for hanging nodes! + // Then we set up a map between // face iterators and their jump @@ -2197,73 +2602,22 @@ run_simulation (LaplaceSolver::Base &solver, -void -create_triangulation (Triangulation<2> &tria) -{ - const Point<2> - vertices[16] = { Point<2> (-1., -1.), - Point<2> (-1./3, -1.), - Point<2> (+1./3, -1.), - Point<2> (+1, -1.), - Point<2> (-1., -1./3.), - Point<2> (-1./3, -1./3.), - Point<2> (+1./3, -1./3.), - Point<2> (+1, -1./3.), - Point<2> (-1., 1./3.), - Point<2> (-1./3, 1./3.), - Point<2> (+1./3, 1./3.), - Point<2> (+1, 1./3.), - Point<2> (-1., 1.), - Point<2> (-1./3, 1.), - Point<2> (+1./3, 1.), - Point<2> (+1, 1.) }; - - const int cell_vertices[8][4] = {{0, 1, 5, 4}, - {1, 2, 6, 5}, - {2, 3, 7, 6}, - {4, 5, 9, 8}, - {6, 7, 11, 10}, - {8,9,13,12}, - {9,10,14,13}, - {10,11,15,14}}; - - std::vector > cells (8, CellData<2>()); - - for (unsigned int i=0; i<8; ++i) - { - for (unsigned int j=0; j<4; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - }; - - tria.create_triangulation (std::vector >(&vertices[0], &vertices[16]), - cells, - SubCellData()); // no boundary information -}; - - template -void solve_problem (const std::string &solver_name) +void solve_problem () { - const std::string header = "Running tests with \"" + solver_name + - "\" refinement criterion:"; - std::cout << header << std::endl - << std::string (header.size(), '-') << std::endl; - Triangulation triangulation (Triangulation::maximum_smoothing); -// create_triangulation (triangulation); - GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (5); const FE_Q primal_fe(1); const FE_Q dual_fe(2); const QGauss4 quadrature; const QGauss4 face_quadrature; const Data::SetUpBase *data = - new Data::SetUp,dim> (); + new Data::SetUp,dim> (); - const Point evaluation_point(0.5,0.5); + data->create_coarse_grid (triangulation); + + const Point evaluation_point(3./4.,3./4.); const DualFunctional::PointValueEvaluation dual_functional (evaluation_point); @@ -2274,15 +2628,18 @@ void solve_problem (const std::string &solver_name) quadrature, face_quadrature, data->get_right_hand_side(), - data->get_solution(), + data->get_boundary_values(), dual_functional); TableHandler results_table; Evaluation::PointValueEvaluation - postprocessor1 (Point(0.5,0.5), results_table); + postprocessor1 (Point(3./4.,3./4.), results_table); + Evaluation::GridOutput + postprocessor2 ("grid"); std::list *> postprocessor_list; postprocessor_list.push_back (&postprocessor1); + postprocessor_list.push_back (&postprocessor2); run_simulation (*solver, postprocessor_list); @@ -2300,7 +2657,7 @@ int main () { deallog.depth_console (0); - solve_problem<2> ("global"); + solve_problem<2> (); // solve_problem<2> ("kelly"); } catch (std::exception &exc) -- 2.39.5