From: bangerth Date: Tue, 17 Oct 2006 04:05:33 +0000 (+0000) Subject: Undo previous change, at least for those files that I didn't want to be affected X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f52279d39b373f17e4b5f72e39e45e0e1483ef08;p=dealii-svn.git Undo previous change, at least for those files that I didn't want to be affected git-svn-id: https://svn.dealii.org/trunk@14008 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-23/Makefile b/deal.II/examples/step-23/Makefile index b28e6264d9..f291d79ad6 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 = 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/Makefile b/deal.II/examples/step-26/Makefile index 31309f018c..7f74cb9645 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 = 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/step-26.cc b/deal.II/examples/step-26/step-26.cc index aa363f6256..0cffe5165b 100644 --- a/deal.II/examples/step-26/step-26.cc +++ b/deal.II/examples/step-26/step-26.cc @@ -55,13 +55,12 @@ #include -class PointCloudSurface : public StraightBoundary<3> +class PointDefinedSurface : public StraightBoundary<3> { public: - /** - * Constructor. - */ - PointCloudSurface (const std::string &filename); + PointDefinedSurface (const std::string &filename); + + Point<3> closest_point (const Point<3> &p) const; /** * Let the new point be the @@ -126,23 +125,12 @@ class PointCloudSurface : 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; }; -PointCloudSurface::PointCloudSurface (const std::string &filename) +PointDefinedSurface::PointDefinedSurface (const std::string &filename) { // first read in all the points { @@ -249,7 +237,7 @@ PointCloudSurface::PointCloudSurface (const std::string &filename) Point<3> -PointCloudSurface::closest_point (const Point<3> &p) const +PointDefinedSurface::closest_point (const Point<3> &p) const { double distance = p.distance (point_list[0]); Point<3> point = point_list[0]; @@ -270,7 +258,7 @@ PointCloudSurface::closest_point (const Point<3> &p) const Point<3> -PointCloudSurface:: +PointDefinedSurface:: get_new_point_on_line (const Triangulation<3>::line_iterator &line) const { return closest_point (StraightBoundary<3>::get_new_point_on_line (line)); @@ -279,7 +267,7 @@ get_new_point_on_line (const Triangulation<3>::line_iterator &line) const Point<3> -PointCloudSurface:: +PointDefinedSurface:: get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const { return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad)); @@ -288,7 +276,7 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const void -PointCloudSurface:: +PointDefinedSurface:: get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, std::vector > &points) const { @@ -301,7 +289,7 @@ get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line, void -PointCloudSurface:: +PointDefinedSurface:: get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, std::vector > &points) const { @@ -313,7 +301,7 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, -PointCloudSurface pds("surface-points"); +PointDefinedSurface pds("surface-points"); @@ -504,7 +492,7 @@ void LaplaceProblem::make_grid_and_dofs () triangulation.begin()->vertex(v)[1], 0)); - for (unsigned int i=0; i<4; ++i) + for (unsigned int i=0; i<7; ++i) { for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); @@ -609,9 +597,7 @@ void LaplaceProblem::assemble_system () template void LaplaceProblem::solve () { - // NEW - SolverControl solver_control (dof_handler.n_dofs(), - 1e-12*system_rhs.l2_norm()); + SolverControl solver_control (1000, 1e-12); SolverCG<> cg (solver_control); PreconditionSSOR<> preconditioner; diff --git a/deal.II/examples/step-4/Makefile b/deal.II/examples/step-4/Makefile index 5380806d90..366b8c588d 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 = 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-4/step-4.cc b/deal.II/examples/step-4/step-4.cc index ea78dc8f90..d6ff4ec32a 100644 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@ -50,266 +50,6 @@ #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 @@ -584,48 +324,16 @@ LaplaceProblem::LaplaceProblem () : template void LaplaceProblem::make_grid_and_dofs () { - 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); + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); - - 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; + std::cout << " 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: " @@ -858,14 +566,19 @@ void LaplaceProblem::assemble_system () template void LaplaceProblem::solve () { -// 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); + 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; } @@ -1011,6 +724,11 @@ 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 053a398931..5cecab44d4 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-3d.g) \ +libs.g = $(lib-deal2-2d.g) \ $(lib-lac.g) \ $(lib-base.g) -libs.o = $(lib-deal2-3d.o) \ +libs.o = $(lib-deal2-2d.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 4fc5c06aeb..0caaae616b 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("pig.inp"); + std::ifstream input_file("circle-grid.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,6 +907,40 @@ 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 @@ -921,9 +955,9 @@ void LaplaceProblem::run () << triangulation.n_cells() << std::endl; -// setup_system (); -// assemble_system (); -// solve (); + setup_system (); + assemble_system (); + solve (); output_results (cycle); } } @@ -939,8 +973,8 @@ int main () { deallog.depth_console (0); - LaplaceProblem<3> laplace_problem_3d; - laplace_problem_3d.run (); + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); // Finally, we have promised to // trigger an exception in the