]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Example update.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 3 Apr 1998 17:43:11 +0000 (17:43 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 3 Apr 1998 17:43:11 +0000 (17:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@133 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/poisson/equation.cc
deal.II/deal.II/Attic/examples/poisson/poisson.h
deal.II/deal.II/Attic/examples/poisson/poisson.prm
deal.II/deal.II/Attic/examples/poisson/problem.cc
tests/big-tests/poisson/equation.cc
tests/big-tests/poisson/poisson.h
tests/big-tests/poisson/poisson.prm
tests/big-tests/poisson/problem.cc

index c995c110e085b6312c3cc82493620343c30a58cf..405d59704fa18a0c05ed61d2f7a256993a8da1d9 100644 (file)
@@ -4,24 +4,6 @@
 
 
 
-template <int dim>
-double RightHandSide<dim>::operator () (const Point<dim> &p) const {
-  const double pi = 3.1415926536;
-  switch (dim) 
-    {
-      case 1:
-           return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
-      case 2:
-           return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) +
-                   2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) +
-                   5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) -
-                   2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi);
-      default:
-           return 0;
-    };
-};
-
-
 
 
 
@@ -83,9 +65,6 @@ void PoissonEquation<dim>::assemble (dVector             &,
 
 
 
-template class RightHandSide<1>;
-template class RightHandSide<2>;
-
 
 template class PoissonEquation<1>;
 template class PoissonEquation<2>;
index d41e6bf43209fb416abd946c0c4f87a972e9fb0f..74224e3837c0dab3b9fe686ab842631f1754793a 100644 (file)
@@ -34,21 +34,6 @@ extern "C" {
 
 
 
-template <int dim>
-class RightHandSide : public Function<dim> {
-  public:
-                                    /**
-                                     * Return the value of the function
-                                     * at the given point.
-                                     */
-    virtual double operator () (const Point<dim> &p) const;
-};
-
-
-
-
-
-
 template <int dim>
 class PoissonEquation :  public Equation<dim> {
   public:
@@ -89,14 +74,18 @@ class PoissonProblem : public ProblemBase<dim>,
 
 
     bool make_grid (ParameterHandler &prm);
-    void make_ball_grid ();
-    void make_curved_line_grid ();
     void make_zoom_in_grid ();
     void make_random_grid ();
 
+    bool set_right_hand_side (ParameterHandler &prm);
+    bool set_boundary_values (ParameterHandler &prm);
+    
   protected:
     Triangulation<dim> *tria;
     DoFHandler<dim>    *dof;
+    
+    Function<dim>      *rhs;
+    Function<dim>      *boundary_values;
 };
 
 
index fb9f151be991efdf994f2e8e0bf3968ff585a7c2..2bebf9b6d429f0e7ba30d7a39614c7feebef0e03 100644 (file)
@@ -1,3 +1,5 @@
 set Test run          = { zoom in | ball | curved line | random }
-set Global refinement = {{ 2 | 4 | 6 | 0 }}
+set Global refinement = {{ 2 | 5 | 6 | 0 }}
+set Right hand side   = {{ zero | zero | trigpoly | constant }}
+set Boundary values   = {{ sine | sine | zero | zero }}
 set Output file       = gnuplot.{{ zoom_in | ball | curved_line | random }}
index 77fcbcf35d51241242e51e4c1a16d4a93eaec89d..fcc3a74ee3e3fea22e003cfb9530da0c23617211 100644 (file)
@@ -5,6 +5,52 @@
 
 
 
+template <int dim>
+class BoundaryValuesSine : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const {
+      return cos(2*3.1415926536*p(0))*cos(2*3.1415926536*p(1));
+    };
+    
+
+                                    /**
+                                     * Set #values# to the point values
+                                     * of the function at the #points#.
+                                     * It is assumed that #values# be
+                                     * empty.
+                                     */
+    virtual void value_list (const vector<Point<dim> > &points,
+                            vector<double>            &values) const {
+      for (unsigned int i=0; i<points.size(); ++i) 
+       values.push_back (cos(2*3.1415926536*points[i](0)) *
+                         cos(2*3.1415926536*points[i](1)));
+    };
+};
+
+
+
+
+template <int dim>
+class RHSTrigPoly : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const;
+};
+
+
+
+
+
+
+
+
 
 template <int dim>
 class CurvedLine :
@@ -34,17 +80,57 @@ class CurvedLine :
 
 
 
+template <int dim>
+double RHSTrigPoly<dim>::operator () (const Point<dim> &p) const {
+  const double pi = 3.1415926536;
+  switch (dim) 
+    {
+      case 1:
+           return p(0)*p(0)*cos(2*pi*p(0));
+      case 2:
+           return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) +
+                   2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) +
+                   5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) -
+                   2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi);
+      default:
+           return 0;
+    };
+};
+
+
+
+
+
 template <int dim>
 PoissonProblem<dim>::PoissonProblem () :
-               tria(0), dof(0) {};
+               tria(0), dof(0), rhs(0), boundary_values(0) {};
 
 
 
 
 template <int dim>
 void PoissonProblem<dim>::clear () {
-  if (tria != 0) delete tria;
-  if (dof  != 0) delete dof;
+  if (tria != 0) {
+    delete tria;
+    tria = 0;
+  };
+  
+  if (dof != 0) {
+    delete dof;
+    dof = 0;
+  };
+
+  if (rhs != 0) 
+    {
+      delete rhs;
+      rhs = 0;
+    };
+
+  if (boundary_values != 0) 
+    {
+      delete boundary_values;
+      boundary_values = 0;
+    };
 
   ProblemBase<dim>::clear ();
 };
@@ -67,12 +153,17 @@ void PoissonProblem<dim>::create_new (const unsigned int) {
 template <int dim>
 void PoissonProblem<dim>::declare_parameters (ParameterHandler &prm) {
   if (dim>=2)
-    prm.declare_entry ("Test run", "zoom in", "zoom in\\|ball\\|curved line\\|random");
+    prm.declare_entry ("Test run", "zoom in",
+                      "tensor\\|zoom in\\|ball\\|curved line\\|random");
   else
-    prm.declare_entry ("Test run", "zoom in", "zoom in\\|random");
+    prm.declare_entry ("Test run", "zoom in", "tensor\\|zoom in\\|random");
 
   prm.declare_entry ("Global refinement", "0",
                     ParameterHandler::RegularExpressions::Integer);
+  prm.declare_entry ("Right hand side", "zero",
+                    "zero\\|constant\\|trigpoly");
+  prm.declare_entry ("Boundary values", "zero",
+                    "zero\\|sine");
   prm.declare_entry ("Output file", "gnuplot.1");
 };
 
@@ -90,11 +181,13 @@ bool PoissonProblem<dim>::make_grid (ParameterHandler &prm) {
       if (test=="curved line") test_case = 3;
       else
        if (test=="random") test_case = 4;
-       else 
-         {
-           cerr << "This test seems not to be implemented!" << endl;
-           return false;
-         };
+       else
+         if (test=="tensor") test_case = 5;
+         else
+           {
+             cerr << "This test seems not to be implemented!" << endl;
+             return false;
+           };
 
   switch (test_case) 
     {
@@ -102,14 +195,29 @@ bool PoissonProblem<dim>::make_grid (ParameterHandler &prm) {
            make_zoom_in_grid ();
            break;
       case 2:
-           make_ball_grid ();
+                                            // make ball grid around origin with
+                                            // unit radius
+      {
+           static const Point<dim> origin;
+           static const HyperBallBoundary<dim> boundary(origin, 1.);
+           tria->create_hyper_ball (origin, 1.);
+           tria->set_boundary (&boundary);
            break;
+      };
       case 3:
-           make_curved_line_grid ();
+                                            // set the boundary function
+      {
+           static const CurvedLine<dim> boundary;
+           tria->create_hypercube ();
+           tria->set_boundary (&boundary);
            break;
+      };
       case 4:
            make_random_grid ();
            break;
+      case 5:
+           tria->create_hypercube ();
+           break;
       default:
            return false;
     };
@@ -152,30 +260,6 @@ void PoissonProblem<dim>::make_zoom_in_grid () {
 
 
 
-template <int dim>
-void PoissonProblem<dim>::make_ball_grid () {
-                                  // make ball grid around origin with
-                                  // unit radius
-  const Point<dim> origin;
-  static const HyperBallBoundary<dim> boundary(origin, 1.);
-  
-  tria->create_hyper_ball (origin, 1.);
-  tria->set_boundary (&boundary);
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::make_curved_line_grid () {
-                                  // set the boundary function
-  static const CurvedLine<dim> boundary;
-
-  tria->create_hypercube ();
-  tria->set_boundary (&boundary);
-};
-
-
 
 template <int dim>
 void PoissonProblem<dim>::make_random_grid () {
@@ -203,9 +287,52 @@ void PoissonProblem<dim>::make_random_grid () {
       tria->execute_refinement ();
     };
 };
+  
+
+
+
+template <int dim>
+bool PoissonProblem<dim>::set_right_hand_side (ParameterHandler &prm) {
+  String rhs_name = prm.get ("Right hand side");
+
+  if (rhs_name == "zero")
+    rhs = new ZeroFunction<dim>();
+  else
+    if (rhs_name == "constant")
+      rhs = new ConstantFunction<dim>(1.);
+    else
+      if (rhs_name == "trigpoly")
+       rhs = new RHSTrigPoly<dim>();
+      else
+       return false;
+
+  if (rhs != 0)
+    return true;
+  else
+    return false;
+};
 
 
 
+template <int dim>
+bool PoissonProblem<dim>::set_boundary_values (ParameterHandler &prm) {
+  String bv_name = prm.get ("Boundary values");
+
+  if (bv_name == "zero")
+    boundary_values = new ZeroFunction<dim> ();
+  else
+    if (bv_name == "sine")
+      boundary_values = new BoundaryValuesSine<dim> ();
+    else
+      return false;
+
+  if (boundary_values != 0)
+    return true;
+  else
+    return false;
+};
+
+
 
 
 template <int dim>
@@ -217,10 +344,14 @@ void PoissonProblem<dim>::run (ParameterHandler &prm) {
   if (!make_grid (prm))
     return;
   
+  if (!set_right_hand_side (prm))
+    return;
 
+  if (!set_boundary_values (prm))
+    return;
+  
   FELinear<dim>                   fe;
-  static const RightHandSide<dim> rhs;
-  PoissonEquation<dim>            equation (rhs);
+  PoissonEquation<dim>            equation (*rhs);
   QGauss4<dim>                    quadrature;
   
   cout << "    Distributing dofs... "; 
@@ -233,8 +364,7 @@ void PoissonProblem<dim>::run (ParameterHandler &prm) {
   update_flags.jacobians = update_flags.JxW_values = true;
   
   ProblemBase<dim>::DirichletBC dirichlet_bc;
-  ZeroFunction<dim> zero;
-  dirichlet_bc[0] = &zero;
+  dirichlet_bc[0] = boundary_values;
   assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
 
   cout << "    Solving..." << endl;
index c995c110e085b6312c3cc82493620343c30a58cf..405d59704fa18a0c05ed61d2f7a256993a8da1d9 100644 (file)
@@ -4,24 +4,6 @@
 
 
 
-template <int dim>
-double RightHandSide<dim>::operator () (const Point<dim> &p) const {
-  const double pi = 3.1415926536;
-  switch (dim) 
-    {
-      case 1:
-           return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3;
-      case 2:
-           return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) +
-                   2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) +
-                   5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) -
-                   2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi);
-      default:
-           return 0;
-    };
-};
-
-
 
 
 
@@ -83,9 +65,6 @@ void PoissonEquation<dim>::assemble (dVector             &,
 
 
 
-template class RightHandSide<1>;
-template class RightHandSide<2>;
-
 
 template class PoissonEquation<1>;
 template class PoissonEquation<2>;
index d41e6bf43209fb416abd946c0c4f87a972e9fb0f..74224e3837c0dab3b9fe686ab842631f1754793a 100644 (file)
@@ -34,21 +34,6 @@ extern "C" {
 
 
 
-template <int dim>
-class RightHandSide : public Function<dim> {
-  public:
-                                    /**
-                                     * Return the value of the function
-                                     * at the given point.
-                                     */
-    virtual double operator () (const Point<dim> &p) const;
-};
-
-
-
-
-
-
 template <int dim>
 class PoissonEquation :  public Equation<dim> {
   public:
@@ -89,14 +74,18 @@ class PoissonProblem : public ProblemBase<dim>,
 
 
     bool make_grid (ParameterHandler &prm);
-    void make_ball_grid ();
-    void make_curved_line_grid ();
     void make_zoom_in_grid ();
     void make_random_grid ();
 
+    bool set_right_hand_side (ParameterHandler &prm);
+    bool set_boundary_values (ParameterHandler &prm);
+    
   protected:
     Triangulation<dim> *tria;
     DoFHandler<dim>    *dof;
+    
+    Function<dim>      *rhs;
+    Function<dim>      *boundary_values;
 };
 
 
index fb9f151be991efdf994f2e8e0bf3968ff585a7c2..2bebf9b6d429f0e7ba30d7a39614c7feebef0e03 100644 (file)
@@ -1,3 +1,5 @@
 set Test run          = { zoom in | ball | curved line | random }
-set Global refinement = {{ 2 | 4 | 6 | 0 }}
+set Global refinement = {{ 2 | 5 | 6 | 0 }}
+set Right hand side   = {{ zero | zero | trigpoly | constant }}
+set Boundary values   = {{ sine | sine | zero | zero }}
 set Output file       = gnuplot.{{ zoom_in | ball | curved_line | random }}
index 77fcbcf35d51241242e51e4c1a16d4a93eaec89d..fcc3a74ee3e3fea22e003cfb9530da0c23617211 100644 (file)
@@ -5,6 +5,52 @@
 
 
 
+template <int dim>
+class BoundaryValuesSine : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const {
+      return cos(2*3.1415926536*p(0))*cos(2*3.1415926536*p(1));
+    };
+    
+
+                                    /**
+                                     * Set #values# to the point values
+                                     * of the function at the #points#.
+                                     * It is assumed that #values# be
+                                     * empty.
+                                     */
+    virtual void value_list (const vector<Point<dim> > &points,
+                            vector<double>            &values) const {
+      for (unsigned int i=0; i<points.size(); ++i) 
+       values.push_back (cos(2*3.1415926536*points[i](0)) *
+                         cos(2*3.1415926536*points[i](1)));
+    };
+};
+
+
+
+
+template <int dim>
+class RHSTrigPoly : public Function<dim> {
+  public:
+                                    /**
+                                     * Return the value of the function
+                                     * at the given point.
+                                     */
+    virtual double operator () (const Point<dim> &p) const;
+};
+
+
+
+
+
+
+
+
 
 template <int dim>
 class CurvedLine :
@@ -34,17 +80,57 @@ class CurvedLine :
 
 
 
+template <int dim>
+double RHSTrigPoly<dim>::operator () (const Point<dim> &p) const {
+  const double pi = 3.1415926536;
+  switch (dim) 
+    {
+      case 1:
+           return p(0)*p(0)*cos(2*pi*p(0));
+      case 2:
+           return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) +
+                   2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) +
+                   5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) -
+                   2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi);
+      default:
+           return 0;
+    };
+};
+
+
+
+
+
 template <int dim>
 PoissonProblem<dim>::PoissonProblem () :
-               tria(0), dof(0) {};
+               tria(0), dof(0), rhs(0), boundary_values(0) {};
 
 
 
 
 template <int dim>
 void PoissonProblem<dim>::clear () {
-  if (tria != 0) delete tria;
-  if (dof  != 0) delete dof;
+  if (tria != 0) {
+    delete tria;
+    tria = 0;
+  };
+  
+  if (dof != 0) {
+    delete dof;
+    dof = 0;
+  };
+
+  if (rhs != 0) 
+    {
+      delete rhs;
+      rhs = 0;
+    };
+
+  if (boundary_values != 0) 
+    {
+      delete boundary_values;
+      boundary_values = 0;
+    };
 
   ProblemBase<dim>::clear ();
 };
@@ -67,12 +153,17 @@ void PoissonProblem<dim>::create_new (const unsigned int) {
 template <int dim>
 void PoissonProblem<dim>::declare_parameters (ParameterHandler &prm) {
   if (dim>=2)
-    prm.declare_entry ("Test run", "zoom in", "zoom in\\|ball\\|curved line\\|random");
+    prm.declare_entry ("Test run", "zoom in",
+                      "tensor\\|zoom in\\|ball\\|curved line\\|random");
   else
-    prm.declare_entry ("Test run", "zoom in", "zoom in\\|random");
+    prm.declare_entry ("Test run", "zoom in", "tensor\\|zoom in\\|random");
 
   prm.declare_entry ("Global refinement", "0",
                     ParameterHandler::RegularExpressions::Integer);
+  prm.declare_entry ("Right hand side", "zero",
+                    "zero\\|constant\\|trigpoly");
+  prm.declare_entry ("Boundary values", "zero",
+                    "zero\\|sine");
   prm.declare_entry ("Output file", "gnuplot.1");
 };
 
@@ -90,11 +181,13 @@ bool PoissonProblem<dim>::make_grid (ParameterHandler &prm) {
       if (test=="curved line") test_case = 3;
       else
        if (test=="random") test_case = 4;
-       else 
-         {
-           cerr << "This test seems not to be implemented!" << endl;
-           return false;
-         };
+       else
+         if (test=="tensor") test_case = 5;
+         else
+           {
+             cerr << "This test seems not to be implemented!" << endl;
+             return false;
+           };
 
   switch (test_case) 
     {
@@ -102,14 +195,29 @@ bool PoissonProblem<dim>::make_grid (ParameterHandler &prm) {
            make_zoom_in_grid ();
            break;
       case 2:
-           make_ball_grid ();
+                                            // make ball grid around origin with
+                                            // unit radius
+      {
+           static const Point<dim> origin;
+           static const HyperBallBoundary<dim> boundary(origin, 1.);
+           tria->create_hyper_ball (origin, 1.);
+           tria->set_boundary (&boundary);
            break;
+      };
       case 3:
-           make_curved_line_grid ();
+                                            // set the boundary function
+      {
+           static const CurvedLine<dim> boundary;
+           tria->create_hypercube ();
+           tria->set_boundary (&boundary);
            break;
+      };
       case 4:
            make_random_grid ();
            break;
+      case 5:
+           tria->create_hypercube ();
+           break;
       default:
            return false;
     };
@@ -152,30 +260,6 @@ void PoissonProblem<dim>::make_zoom_in_grid () {
 
 
 
-template <int dim>
-void PoissonProblem<dim>::make_ball_grid () {
-                                  // make ball grid around origin with
-                                  // unit radius
-  const Point<dim> origin;
-  static const HyperBallBoundary<dim> boundary(origin, 1.);
-  
-  tria->create_hyper_ball (origin, 1.);
-  tria->set_boundary (&boundary);
-};
-
-
-
-
-template <int dim>
-void PoissonProblem<dim>::make_curved_line_grid () {
-                                  // set the boundary function
-  static const CurvedLine<dim> boundary;
-
-  tria->create_hypercube ();
-  tria->set_boundary (&boundary);
-};
-
-
 
 template <int dim>
 void PoissonProblem<dim>::make_random_grid () {
@@ -203,9 +287,52 @@ void PoissonProblem<dim>::make_random_grid () {
       tria->execute_refinement ();
     };
 };
+  
+
+
+
+template <int dim>
+bool PoissonProblem<dim>::set_right_hand_side (ParameterHandler &prm) {
+  String rhs_name = prm.get ("Right hand side");
+
+  if (rhs_name == "zero")
+    rhs = new ZeroFunction<dim>();
+  else
+    if (rhs_name == "constant")
+      rhs = new ConstantFunction<dim>(1.);
+    else
+      if (rhs_name == "trigpoly")
+       rhs = new RHSTrigPoly<dim>();
+      else
+       return false;
+
+  if (rhs != 0)
+    return true;
+  else
+    return false;
+};
 
 
 
+template <int dim>
+bool PoissonProblem<dim>::set_boundary_values (ParameterHandler &prm) {
+  String bv_name = prm.get ("Boundary values");
+
+  if (bv_name == "zero")
+    boundary_values = new ZeroFunction<dim> ();
+  else
+    if (bv_name == "sine")
+      boundary_values = new BoundaryValuesSine<dim> ();
+    else
+      return false;
+
+  if (boundary_values != 0)
+    return true;
+  else
+    return false;
+};
+
+
 
 
 template <int dim>
@@ -217,10 +344,14 @@ void PoissonProblem<dim>::run (ParameterHandler &prm) {
   if (!make_grid (prm))
     return;
   
+  if (!set_right_hand_side (prm))
+    return;
 
+  if (!set_boundary_values (prm))
+    return;
+  
   FELinear<dim>                   fe;
-  static const RightHandSide<dim> rhs;
-  PoissonEquation<dim>            equation (rhs);
+  PoissonEquation<dim>            equation (*rhs);
   QGauss4<dim>                    quadrature;
   
   cout << "    Distributing dofs... "; 
@@ -233,8 +364,7 @@ void PoissonProblem<dim>::run (ParameterHandler &prm) {
   update_flags.jacobians = update_flags.JxW_values = true;
   
   ProblemBase<dim>::DirichletBC dirichlet_bc;
-  ZeroFunction<dim> zero;
-  dirichlet_bc[0] = &zero;
+  dirichlet_bc[0] = boundary_values;
   assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
 
   cout << "    Solving..." << endl;

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.