]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Undo previous change, at least for those files that I didn't want to be affected
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 04:05:33 +0000 (04:05 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 04:05:33 +0000 (04:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@14008 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-23/Makefile
deal.II/examples/step-26/Makefile
deal.II/examples/step-26/step-26.cc
deal.II/examples/step-4/Makefile
deal.II/examples/step-4/step-4.cc
deal.II/examples/step-5/Makefile
deal.II/examples/step-5/step-5.cc

index b28e6264d988275fa1cd60b125e68dac0e7829d6..f291d79ad6ee6bbdef723be0389a98dd16fbbdf0 100644 (file)
@@ -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
index 31309f018cdf0b1d400b326df9bdc4f0b65ed3b5..7f74cb964592722e24ff719f737444a7ff26746f 100644 (file)
@@ -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
index aa363f6256a50c9526e285a67bee38db349fe3c4..0cffe5165b0420ac2e47fb55103718e160b1dbcc 100644 (file)
 #include <grid/tria_boundary.h>
 
 
-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<Point<3> > &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<3> > 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<Point<3> > &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<Point<3> > &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<dim>::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<dim>::active_cell_iterator
             cell = triangulation.begin_active();
@@ -609,9 +597,7 @@ void LaplaceProblem<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::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;
index 5380806d90163c84bab06d28597209e43d7ea950..366b8c588d31939549016ff5f6a9dcf19f7c97d9 100644 (file)
@@ -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
index ea78dc8f901e02ea5b526ca9c2d690b62f00ec06..d6ff4ec32a652a889f1aa36b8622249ef2de2974 100644 (file)
 #include <base/logstream.h>
 
 
-#include <algorithm>
-#include <numeric>
-#include <grid/tria_boundary.h>
-
-
-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 <tt>n=points.size()</tt>
-                                     * 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<Point<3> > &points) const;
-
-                                    /**
-                                     * Gives <tt>n=points.size()=m*m</tt>
-                                     * points that splits the
-                                     * p{StraightBoundary} quad into
-                                     * <tt>(m+1)(m+1)</tt> 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<Point<3> > &points) const;
-  private:
-    std::vector<Point<3> > 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<point_list.size(); ++i)
-      point_list[i] -= c_o_m;
-  }
-
-                                  // next do a least squares fit to the
-                                  // function ax+by. this leads to the
-                                  // following equations:
-  
-                                  // min f(a,b) = sum_i (zi-a xi - b yi)^2 / 2
-                                  //
-                                  // f_a = sum_i (zi - a xi - b yi) xi = 0
-                                  // f_b = sum_i (zi - a xi - b yi) yi = 0
-                                  //
-                                  // f_a = (sum_i zi xi) - (sum xi^2) a - (sum xi yi) b = 0
-                                  // f_a = (sum_i zi yi) - (sum xi yi) a - (sum yi^2) b = 0
-  {
-    double A[2][2] = {{0,0},{0,0}};
-    double B[2] = {0,0};
-
-    for (unsigned int i=0; i<point_list.size(); ++i)
-      {
-       A[0][0] += point_list[i][0] * point_list[i][0];
-       A[0][1] += point_list[i][0] * point_list[i][1];
-       A[1][1] += point_list[i][1] * point_list[i][1];
-
-       B[0] += point_list[i][0] * point_list[i][2];
-       B[1] += point_list[i][1] * point_list[i][2];
-      }
-
-    const double det = A[0][0]*A[1][1]-2*A[0][1];
-    const double a = (A[1][1] * B[0] - A[0][1] * B[1]) / det;
-    const double b = (A[0][0] * B[1] - A[0][1] * B[0]) / det;
-
-
-                                    // with this information, we can rotate
-                                    // the points so that the corresponding
-                                    // least-squares fit would be the x-y
-                                    // plane
-    const Point<2> 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<point_list.size(); ++i)
-      {
-                                        // we can do that by, for each point,
-                                        // first subtract the points in the
-                                        // plane:
-       point_list[i][2] -= a*point_list[i][0] + b*point_list[i][1];
-
-                                      // we made a mistake here, though:
-                                      // we've shrunk the plan in the
-                                      // direction parallel to the
-                                      // gradient. we will have to correct
-                                      // for this:
-       const Point<2> 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<Point<3> >::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<Point<3> > &points) const
-{
-  StraightBoundary<3>::get_intermediate_points_on_line (line,
-                                                       points);
-  for (unsigned int i=0; i<points.size(); ++i)
-    points[i] = closest_point(points[i]);
-}
-
-
-
-void
-PointDefinedSurface::
-get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
-                                std::vector<Point<3> > &points) const
-{
-  StraightBoundary<3>::get_intermediate_points_on_quad (quad,
-                                                       points);
-  for (unsigned int i=0; i<points.size(); ++i)
-    points[i] = closest_point(points[i]);
-}
-
-
-
-PointDefinedSurface pds("unique-points");
-
-
-
-
-
-
-
-
                                  // @sect3{The <code>LaplaceProblem</code> class template}
 
                                 // This is again the same
@@ -584,48 +324,16 @@ LaplaceProblem<dim>::LaplaceProblem () :
 template <int dim>
 void LaplaceProblem<dim>::make_grid_and_dofs ()
 {
-  GridGenerator::hyper_cube (triangulation, -30, 30);
-
-  for (unsigned int f=0; f<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<dim>::active_cell_iterator
-            cell = triangulation.begin_active();
-          cell != triangulation.end(); ++cell)
-       for (unsigned int f=0; f<GeometryInfo<dim>::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<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::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<dim>::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 ();
index 053a398931afcf40fdf20e69d754df08e2a11646..5cecab44d4351ffa1ed4b96339d2f4fb1a32234b 100644 (file)
@@ -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)
 
index 4fc5c06aeb71a24800aaacc26497d6665f8c2f1e..0caaae616b1b307952f3f293537f964695e4aa92 100644 (file)
@@ -844,7 +844,7 @@ void LaplaceProblem<dim>::run ()
                                           // data in the file:
          GridIn<dim> 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<dim>::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<dim>::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
+                                           // <code>grid_in.read_xxx</code>
+                                           // function. (See the
+                                           // documentation of the
+                                           // <code>GridIn</code> 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<dim> boundary;
+         triangulation.set_boundary (0, boundary);
        }
 
                                       // Now that we have a mesh for
@@ -921,9 +955,9 @@ void LaplaceProblem<dim>::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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.