]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move the presently unfinished step-21 to step-27, to make room for the IMPES tutorial...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 03:59:47 +0000 (03:59 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 17 Oct 2006 03:59:47 +0000 (03:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@14007 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-27/Makefile [moved from deal.II/examples/step-21/Makefile with 100% similarity]
deal.II/examples/step-27/doc/intro.dox [moved from deal.II/examples/step-21/doc/intro.dox with 100% similarity]
deal.II/examples/step-27/doc/results.dox [moved from deal.II/examples/step-21/doc/results.dox with 100% similarity]
deal.II/examples/step-27/step-21.cc [moved from deal.II/examples/step-21/step-21.cc with 100% similarity]
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 f291d79ad6ee6bbdef723be0389a98dd16fbbdf0..b28e6264d988275fa1cd60b125e68dac0e7829d6 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 7f74cb964592722e24ff719f737444a7ff26746f..31309f018cdf0b1d400b326df9bdc4f0b65ed3b5 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 0cffe5165b0420ac2e47fb55103718e160b1dbcc..aa363f6256a50c9526e285a67bee38db349fe3c4 100644 (file)
 #include <grid/tria_boundary.h>
 
 
-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<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;
 };
 
 
-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<Point<3> > &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<Point<3> > &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<dim>::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<dim>::active_cell_iterator
             cell = triangulation.begin_active();
@@ -597,7 +609,9 @@ void LaplaceProblem<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::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;
index 366b8c588d31939549016ff5f6a9dcf19f7c97d9..5380806d90163c84bab06d28597209e43d7ea950 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 d6ff4ec32a652a889f1aa36b8622249ef2de2974..ea78dc8f901e02ea5b526ca9c2d690b62f00ec06 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
@@ -324,16 +584,48 @@ LaplaceProblem<dim>::LaplaceProblem () :
 template <int dim>
 void LaplaceProblem<dim>::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<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);
   
-  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<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;
 
+    }
+  
+  
   dof_handler.distribute_dofs (fe);
 
   std::cout << "   Number of degrees of freedom: "
@@ -566,19 +858,14 @@ void LaplaceProblem<dim>::assemble_system ()
 template <int dim>
 void LaplaceProblem<dim>::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<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 5cecab44d4351ffa1ed4b96339d2f4fb1a32234b..053a398931afcf40fdf20e69d754df08e2a11646 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-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)
 
index 0caaae616b1b307952f3f293537f964695e4aa92..4fc5c06aeb71a24800aaacc26497d6665f8c2f1e 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("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<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,40 +907,6 @@ 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
@@ -955,9 +921,9 @@ void LaplaceProblem<dim>::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

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.