From: bangerth Date: Tue, 17 Oct 2006 03:59:47 +0000 (+0000) Subject: Move the presently unfinished step-21 to step-27, to make room for the IMPES tutorial... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=91e4633d6f18299ff3d5b926a7a7799324298ee8;p=dealii-svn.git Move the presently unfinished step-21 to step-27, to make room for the IMPES tutorial program that naturally follows step-20 git-svn-id: https://svn.dealii.org/trunk@14007 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-23/Makefile b/deal.II/examples/step-23/Makefile index f291d79ad6..b28e6264d9 100644 --- a/deal.II/examples/step-23/Makefile +++ b/deal.II/examples/step-23/Makefile @@ -14,7 +14,7 @@ target = $(basename $(shell echo step-*.cc)) # run-time checking of parameters and internal states is performed, so # you should set this value to `on' while you develop your program, # and to `off' when running production computations. -debug-mode = on +debug-mode = off # As third field, we need to give the path to the top-level deal.II diff --git a/deal.II/examples/step-26/Makefile b/deal.II/examples/step-26/Makefile index 7f74cb9645..31309f018c 100644 --- a/deal.II/examples/step-26/Makefile +++ b/deal.II/examples/step-26/Makefile @@ -14,7 +14,7 @@ target = $(basename $(shell echo step-*.cc)) # run-time checking of parameters and internal states is performed, so # you should set this value to `on' while you develop your program, # and to `off' when running production computations. -debug-mode = off +debug-mode = on # As third field, we need to give the path to the top-level deal.II diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc index 0cffe5165b..aa363f6256 100644 --- a/deal.II/examples/step-26/step-26.cc +++ b/deal.II/examples/step-26/step-26.cc @@ -55,12 +55,13 @@ #include -class PointDefinedSurface : public StraightBoundary<3> +class PointCloudSurface : public StraightBoundary<3> { public: - PointDefinedSurface (const std::string &filename); - - Point<3> closest_point (const Point<3> &p) const; + /** + * Constructor. + */ + PointCloudSurface (const std::string &filename); /** * Let the new point be the @@ -125,12 +126,23 @@ class PointDefinedSurface : public StraightBoundary<3> virtual void get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, std::vector > &points) const; + + /** + * A function that, given a point @p p, + * returns the closest point on the + * surface defined by the input file. For + * the time being, we simply return the + * closest point in the point cloud, + * rather than doing any sort of + * interpolation. + */ + Point<3> closest_point (const Point<3> &p) const; private: std::vector > point_list; }; -PointDefinedSurface::PointDefinedSurface (const std::string &filename) +PointCloudSurface::PointCloudSurface (const std::string &filename) { // first read in all the points { @@ -237,7 +249,7 @@ PointDefinedSurface::PointDefinedSurface (const std::string &filename) Point<3> -PointDefinedSurface::closest_point (const Point<3> &p) const +PointCloudSurface::closest_point (const Point<3> &p) const { double distance = p.distance (point_list[0]); Point<3> point = point_list[0]; @@ -258,7 +270,7 @@ PointDefinedSurface::closest_point (const Point<3> &p) const Point<3> -PointDefinedSurface:: +PointCloudSurface:: get_new_point_on_line (const Triangulation<3>::line_iterator &line) const { return closest_point (StraightBoundary<3>::get_new_point_on_line (line)); @@ -267,7 +279,7 @@ get_new_point_on_line (const Triangulation<3>::line_iterator &line) const Point<3> -PointDefinedSurface:: +PointCloudSurface:: get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const { return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad)); @@ -276,7 +288,7 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const void -PointDefinedSurface:: +PointCloudSurface:: get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, std::vector > &points) const { @@ -289,7 +301,7 @@ get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, void -PointDefinedSurface:: +PointCloudSurface:: get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, std::vector > &points) const { @@ -301,7 +313,7 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, -PointDefinedSurface pds("surface-points"); +PointCloudSurface pds("surface-points"); @@ -492,7 +504,7 @@ void LaplaceProblem::make_grid_and_dofs () triangulation.begin()->vertex(v)[1], 0)); - for (unsigned int i=0; i<7; ++i) + for (unsigned int i=0; i<4; ++i) { for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); @@ -597,7 +609,9 @@ void LaplaceProblem::assemble_system () template void LaplaceProblem::solve () { - SolverControl solver_control (1000, 1e-12); + // NEW + SolverControl solver_control (dof_handler.n_dofs(), + 1e-12*system_rhs.l2_norm()); SolverCG<> cg (solver_control); PreconditionSSOR<> preconditioner; diff --git a/deal.II/examples/step-21/Makefile b/deal.II/examples/step-27/Makefile similarity index 100% rename from deal.II/examples/step-21/Makefile rename to deal.II/examples/step-27/Makefile diff --git a/deal.II/examples/step-21/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox similarity index 100% rename from deal.II/examples/step-21/doc/intro.dox rename to deal.II/examples/step-27/doc/intro.dox diff --git a/deal.II/examples/step-21/doc/results.dox b/deal.II/examples/step-27/doc/results.dox similarity index 100% rename from deal.II/examples/step-21/doc/results.dox rename to deal.II/examples/step-27/doc/results.dox diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-27/step-21.cc similarity index 100% rename from deal.II/examples/step-21/step-21.cc rename to deal.II/examples/step-27/step-21.cc diff --git a/deal.II/examples/step-4/Makefile b/deal.II/examples/step-4/Makefile index 366b8c588d..5380806d90 100644 --- a/deal.II/examples/step-4/Makefile +++ b/deal.II/examples/step-4/Makefile @@ -14,7 +14,7 @@ target = $(basename $(shell echo step-*.cc)) # run-time checking of parameters and internal states is performed, so # you should set this value to `on' while you develop your program, # and to `off' when running production computations. -debug-mode = on +debug-mode = off # As third field, we need to give the path to the top-level deal.II diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc index d6ff4ec32a..ea78dc8f90 100644 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@ -50,6 +50,266 @@ #include +#include +#include +#include + + +class PointDefinedSurface : public StraightBoundary<3> +{ + public: + PointDefinedSurface (const std::string &filename); + + Point<3> closest_point (const Point<3> &p) const; + + /** + * Let the new point be the + * arithmetic mean of the two + * vertices of the line. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. + */ + virtual Point<3> + get_new_point_on_line (const Triangulation<3>::line_iterator &line) const; + + /** + * Let the new point be the + * arithmetic mean of the four + * vertices of this quad and the + * four midpoints of the lines, + * which are already created at + * the time of calling this + * function. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class for more + * information. + */ + virtual Point<3> + get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const; + + /** + * Gives n=points.size() + * points that splits the + * p{StraightBoundary} line into + * p{n+1} partitions of equal + * lengths. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, + std::vector > &points) const; + + /** + * Gives n=points.size()=m*m + * points that splits the + * p{StraightBoundary} quad into + * (m+1)(m+1) subquads of equal + * size. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, + std::vector > &points) const; + private: + std::vector > point_list; +}; + + +PointDefinedSurface::PointDefinedSurface (const std::string &filename) +{ + // first read in all the points + { + std::ifstream in (filename.c_str()); + AssertThrow (in, ExcIO()); + + while (in) + { + Point<3> p; + in >> p; + point_list.push_back (p); + } + + AssertThrow (point_list.size() > 1, ExcIO()); + } + + // next fit a linear model through the data + // cloud to rectify it in a local + // coordinate system + // + // the first step is to move the center of + // mass of the points to the origin + { + const Point<3> c_o_m = std::accumulate (point_list.begin(), + point_list.end(), + Point<3>()) / + point_list.size(); + for (unsigned int i=0; i gradient_direction + = Point<2>(a,b) / std::sqrt(a*a+b*b); + const Point<2> orthogonal_direction + = Point<2>(-b,a) / std::sqrt(a*a+b*b); + + const double stretch_factor = std::sqrt(1.+a*a+b*b); + + for (unsigned int i=0; i xy (point_list[i][0], + point_list[i][1]); + const double grad_distance = xy * gradient_direction; + const double orth_distance = xy * orthogonal_direction; + + // we then have to stretch the points + // in the gradient direction. the + // stretch factor is defined above + // (zero if the original plane was + // already the xy plane, infinity if + // it was vertical) + const Point<2> new_xy + = (grad_distance * stretch_factor * gradient_direction + + orth_distance * orthogonal_direction); + point_list[i][0] = new_xy[0]; + point_list[i][1] = new_xy[1]; + } + } +} + + +Point<3> +PointDefinedSurface::closest_point (const Point<3> &p) const +{ + double distance = p.distance (point_list[0]); + Point<3> point = point_list[0]; + + for (std::vector >::const_iterator i=point_list.begin(); + i != point_list.end(); ++i) + { + const double d = p.distance (*i); + if (d < distance) + { + distance = d; + point = *i; + } + } + + return point; +} + + +Point<3> +PointDefinedSurface:: +get_new_point_on_line (const Triangulation<3>::line_iterator &line) const +{ + return closest_point (StraightBoundary<3>::get_new_point_on_line (line)); +} + + + +Point<3> +PointDefinedSurface:: +get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const +{ + return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad)); +} + + + +void +PointDefinedSurface:: +get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, + std::vector > &points) const +{ + StraightBoundary<3>::get_intermediate_points_on_line (line, + points); + for (unsigned int i=0; i::quad_iterator &quad, + std::vector > &points) const +{ + StraightBoundary<3>::get_intermediate_points_on_quad (quad, + points); + for (unsigned int i=0; iLaplaceProblem class template} // This is again the same @@ -324,16 +584,48 @@ LaplaceProblem::LaplaceProblem () : template void LaplaceProblem::make_grid_and_dofs () { - GridGenerator::hyper_cube (triangulation, -1, 1); - triangulation.refine_global (4); + GridGenerator::hyper_cube (triangulation, -30, 30); + + for (unsigned int f=0; f::faces_per_cell; ++f) + if (triangulation.begin()->face(f)->center()[2] > 15) + { + triangulation.begin()->face(f)->set_boundary_indicator (1); + for (unsigned int i=0; i::lines_per_face; ++i) + triangulation.begin()->face(f)->line(i)->set_boundary_indicator (1); + break; + } + triangulation.set_boundary (1, pds); - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl - << " Total number of cells: " - << triangulation.n_cells() - << std::endl; + + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (triangulation.begin()->vertex(v)[2] > 0) + triangulation.begin()->vertex(v) + = pds.closest_point (Point<3>(triangulation.begin()->vertex(v)[0], + triangulation.begin()->vertex(v)[1], + 0)); + + for (unsigned int i=0; i<7; ++i) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->boundary_indicator() == 1) + cell->set_refine_flag (); + + triangulation.execute_coarsening_and_refinement (); + + std::cout << i << std::endl + << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " + << triangulation.n_cells() + << std::endl; + } + + dof_handler.distribute_dofs (fe); std::cout << " Number of degrees of freedom: " @@ -566,19 +858,14 @@ void LaplaceProblem::assemble_system () template void LaplaceProblem::solve () { - SolverControl solver_control (1000, 1e-12); - SolverCG<> cg (solver_control); - cg.solve (system_matrix, solution, system_rhs, - PreconditionIdentity()); - - // We have made one addition, - // though: since we suppress output - // from the linear solvers, we have - // to print the number of - // iterations by hand. - std::cout << " " << solver_control.last_step() - << " CG iterations needed to obtain convergence." - << std::endl; +// SolverControl solver_control (1000, 1e-12); +// SolverCG<> cg (solver_control); + +// PreconditionSSOR<> preconditioner; +// preconditioner.initialize(system_matrix, 1.2); + +// cg.solve (system_matrix, solution, system_rhs, +// preconditioner); } @@ -724,11 +1011,6 @@ void LaplaceProblem::run () int main () { deallog.depth_console (0); - { - LaplaceProblem<2> laplace_problem_2d; - laplace_problem_2d.run (); - } - { LaplaceProblem<3> laplace_problem_3d; laplace_problem_3d.run (); diff --git a/deal.II/examples/step-5/Makefile b/deal.II/examples/step-5/Makefile index 5cecab44d4..053a398931 100644 --- a/deal.II/examples/step-5/Makefile +++ b/deal.II/examples/step-5/Makefile @@ -57,10 +57,10 @@ include $D/common/Make.global_options # # You may need to augment the lists of libraries when compiling your # program for other dimensions, or when using third party libraries -libs.g = $(lib-deal2-2d.g) \ +libs.g = $(lib-deal2-3d.g) \ $(lib-lac.g) \ $(lib-base.g) -libs.o = $(lib-deal2-2d.o) \ +libs.o = $(lib-deal2-3d.o) \ $(lib-lac.o) \ $(lib-base.o) diff --git a/deal.II/examples/step-5/step-5.cc b/deal.II/examples/step-5/step-5.cc index 0caaae616b..4fc5c06aeb 100644 --- a/deal.II/examples/step-5/step-5.cc +++ b/deal.II/examples/step-5/step-5.cc @@ -844,7 +844,7 @@ void LaplaceProblem::run () // data in the file: GridIn grid_in; grid_in.attach_triangulation (triangulation); - std::ifstream input_file("circle-grid.inp"); + std::ifstream input_file("pig.inp"); // We would now like to // read the file. However, // the input file is only @@ -860,7 +860,7 @@ void LaplaceProblem::run () // dimensions, but rather // kill the whole program // if we are not in 2D: - Assert (dim==2, ExcInternalError()); +// Assert (dim==2, ExcInternalError()); // ExcInternalError is a // globally defined // exception, which may be @@ -907,40 +907,6 @@ void LaplaceProblem::run () // visualization program), // for example: grid_in.read_ucd (input_file); - // If you like to use - // another input format, - // you have to use an other - // grid_in.read_xxx - // function. (See the - // documentation of the - // GridIn class to find - // out what input formats - // are presently - // supported.) - - // The grid in the file - // describes a - // circle. Therefore we - // have to use a boundary - // object which tells the - // triangulation where to - // put new points on the - // boundary when the grid - // is refined. This works - // in the same way as in - // the first example. Note - // that the - // HyperBallBoundary - // constructor takes two - // parameters, the center - // of the ball and the - // radius, but that their - // default (the origin and - // 1.0) are the ones which - // we would like to use - // here. - static const HyperBallBoundary boundary; - triangulation.set_boundary (0, boundary); } // Now that we have a mesh for @@ -955,9 +921,9 @@ void LaplaceProblem::run () << triangulation.n_cells() << std::endl; - setup_system (); - assemble_system (); - solve (); +// setup_system (); +// assemble_system (); +// solve (); output_results (cycle); } } @@ -973,8 +939,8 @@ int main () { deallog.depth_console (0); - LaplaceProblem<2> laplace_problem_2d; - laplace_problem_2d.run (); + LaplaceProblem<3> laplace_problem_3d; + laplace_problem_3d.run (); // Finally, we have promised to // trigger an exception in the