]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extend test to 3D, beware of interplay of hanging nodes and boundary values.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 08:38:10 +0000 (08:38 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 15 Sep 2010 08:38:10 +0000 (08:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@21972 0785d39b-7218-0410-832d-ea1e28bc413d

tests/deal.II/inhomogeneous_constraints.cc
tests/deal.II/inhomogeneous_constraints/cmp/generic

index 26caf8a8cac11ec37481b9c997ba75cb9ecae2e4..a7965967a2bb6e8eaa0550e102e2373212acc2cf 100644 (file)
@@ -97,6 +97,29 @@ class LaplaceProblem
 
 
 
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+  public:
+    BoundaryValues () : Function<dim> () {}
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>   &p,
+                          const unsigned int  /*component*/) const
+{
+  double sum = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    sum += std::sin(deal_II_numbers::PI*p[d]);
+  return sum;
+}
+
+
 template <int dim>
 class RightHandSide : public Function<dim>
 {
@@ -126,12 +149,20 @@ LaplaceProblem<dim>::LaplaceProblem ()
                dof_handler (triangulation),
                max_degree (5)
 {
-  for (unsigned int degree=2; degree<=max_degree; ++degree)
-    {
-      fe_collection.push_back (FE_Q<dim>(degree));
-      quadrature_collection.push_back (QGauss<dim>(degree+1));
-      face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
-    }
+  if (dim == 2)
+    for (unsigned int degree=2; degree<=max_degree; ++degree)
+      {
+       fe_collection.push_back (FE_Q<dim>(degree));
+       quadrature_collection.push_back (QGauss<dim>(degree+1));
+       face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
+  else
+    for (unsigned int degree=1; degree<max_degree-1; ++degree)
+      {
+       fe_collection.push_back (FE_Q<dim>(degree));
+       quadrature_collection.push_back (QGauss<dim>(degree+1));
+       face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+      }
 }
 
 
@@ -154,29 +185,36 @@ void LaplaceProblem<dim>::setup_system ()
   hanging_nodes_only.clear ();
   test_all_constraints.clear ();
 
-                                  // add boundary conditions as
-                                  // inhomogeneous constraints here. In
-                                  // contrast to step-27, we choose a
-                                  // constant function with value 1 here.
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                          hanging_nodes_only);
+  hanging_nodes_only.close ();
+
+                                  // add boundary conditions as inhomogeneous
+                                  // constraints here, do it after having
+                                  // added the hanging node constraints in
+                                  // order to be consistent and skip dofs
+                                  // that are already constrained (i.e., are
+                                  // hanging nodes on the boundary in 3D). In
+                                  // contrast to step-27, we choose a sine
+                                  // function.
   {
+    test_all_constraints.merge (hanging_nodes_only);
     std::map<unsigned int,double> boundary_values;
     VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
-                                             ConstantFunction<dim>(1.),
+                                             BoundaryValues<dim>(),
                                              boundary_values);
     std::map<unsigned int,double>::const_iterator boundary_value = boundary_values.begin();
     for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
       {
-       test_all_constraints.add_line(boundary_value->first);
-       test_all_constraints.set_inhomogeneity (boundary_value->first, 
-                                               boundary_value->second);
+       if (!test_all_constraints.is_constrained(boundary_value->first))
+         {
+           test_all_constraints.add_line(boundary_value->first);
+           test_all_constraints.set_inhomogeneity (boundary_value->first, 
+                                                   boundary_value->second);
+         }
       }
   }
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          hanging_nodes_only);
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          test_all_constraints);
-  hanging_nodes_only.close ();
   test_all_constraints.close ();
 
   CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
@@ -302,7 +340,7 @@ void LaplaceProblem<dim>::assemble_reference ()
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            0,
-                                           ConstantFunction<dim>(1.),
+                                           BoundaryValues<dim>(),
                                            boundary_values);
   MatrixTools::apply_boundary_values (boundary_values,
                                      reference_matrix,
@@ -473,12 +511,10 @@ void LaplaceProblem<dim>::solve ()
   cg.solve (reference_matrix, solution, reference_rhs,
            preconditioner);
 
-  // test also distribute function
   Vector<double> solution_test (solution);
-
   hanging_nodes_only.distribute (solution);
 
-  // test also distribute function
+                               // test also distribute function
   test_all_constraints.distribute(solution_test);
   solution_test -= solution;
   deallog << "Distribute error: " << solution_test.l2_norm () << std::endl;
@@ -616,6 +652,14 @@ void LaplaceProblem<2>::create_coarse_grid ()
 }
 
 
+template <>
+void LaplaceProblem<3>::create_coarse_grid ()
+{
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (1);
+}
+
+
 
 template <int dim>
 void LaplaceProblem<dim>::run () 
@@ -628,14 +672,17 @@ void LaplaceProblem<dim>::run ()
       setup_system ();
 
       deallog << std::endl << std::endl
-             << "   Number of active cells:       "
+             << "   Number of active cells:             "
              << triangulation.n_active_cells()
              << std::endl
-             << "   Number of degrees of freedom: "
+             << "   Number of degrees of freedom:       "
              << dof_handler.n_dofs()
              << std::endl
-             << "   Number of constraints       : "
+             << "   Number of hanging node constraints: "
              << hanging_nodes_only.n_constraints()
+             << std::endl
+             << "   Total number of constraints:        "
+             << test_all_constraints.n_constraints()
              << std::endl;
 
       assemble_reference ();
@@ -648,7 +695,8 @@ void LaplaceProblem<dim>::run ()
 }
 
 
-                                  // this function is copied verbatim from step-27
+                                  // this function is copied verbatim from
+                                  // step-27
 template <int dim>
 void
 LaplaceProblem<dim>::
@@ -796,6 +844,18 @@ int main ()
   deallog.depth_console(0);
   deallog.threshold_double(1.e-10);
 
-  LaplaceProblem<2> laplace_problem;
-  laplace_problem.run ();
+  {
+    deallog.push("2d");
+    LaplaceProblem<2> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
+
+  {
+    deallog.push("3d");
+    LaplaceProblem<3> laplace_problem;
+    laplace_problem.run ();
+    deallog.pop();
+  }
 }
+
index ff2340e4a32cebe791fc5887b460a25ad511902a..2a78e7a212f4cfda3b18489d968312d7c9dba7d3 100644 (file)
@@ -1,46 +1,97 @@
 
-DEAL::
-DEAL::
-DEAL::   Number of active cells:       192
-DEAL::   Number of degrees of freedom: 864
-DEAL::   Number of constraints       : 0
-DEAL::  Reference matrix nonzeros: 12672, actually: 8480
-DEAL::  Test matrix 1 nonzeros: 12672, actually: 8480
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL::  Test matrix 2 nonzeros: 12672, actually: 8480
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL:cg::Starting value 20.
-DEAL:cg::Convergence step 25 value 2.5e-07
-DEAL::Distribute error: 0
-DEAL::
-DEAL::
-DEAL::   Number of active cells:       309
-DEAL::   Number of degrees of freedom: 1575
-DEAL::   Number of constraints       : 187
-DEAL::  Reference matrix nonzeros: 25901, actually: 16531
-DEAL::  Test matrix 1 nonzeros: 25901, actually: 16531
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL::  Test matrix 2 nonzeros: 25901, actually: 16531
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL:cg::Starting value 22.
-DEAL:cg::Convergence step 42 value 3.2e-07
-DEAL::Distribute error: 0
-DEAL::
-DEAL::
-DEAL::   Number of active cells:       450
-DEAL::   Number of degrees of freedom: 2692
-DEAL::   Number of constraints       : 475
-DEAL::  Reference matrix nonzeros: 51910, actually: 32730
-DEAL::  Test matrix 1 nonzeros: 51910, actually: 32730
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL::  Test matrix 2 nonzeros: 51910, actually: 32730
-DEAL::  Matrix difference norm: 0
-DEAL::  RHS difference norm: 0
-DEAL:cg::Starting value 25.
-DEAL:cg::Convergence step 61 value 3.8e-07
-DEAL::Distribute error: 0
+DEAL:2d::
+DEAL:2d::
+DEAL:2d::   Number of active cells:             192
+DEAL:2d::   Number of degrees of freedom:       864
+DEAL:2d::   Number of hanging node constraints: 0
+DEAL:2d::   Total number of constraints:        192
+DEAL:2d::  Reference matrix nonzeros: 12672, actually: 8480
+DEAL:2d::  Test matrix 1 nonzeros: 12672, actually: 8480
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d::  Test matrix 2 nonzeros: 12672, actually: 8480
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d:cg::Starting value 18.
+DEAL:2d:cg::Convergence step 26 value 2.2e-07
+DEAL:2d::Distribute error: 0
+DEAL:2d::
+DEAL:2d::
+DEAL:2d::   Number of active cells:             237
+DEAL:2d::   Number of degrees of freedom:       1385
+DEAL:2d::   Number of hanging node constraints: 162
+DEAL:2d::   Total number of constraints:        380
+DEAL:2d::  Reference matrix nonzeros: 24829, actually: 15897
+DEAL:2d::  Test matrix 1 nonzeros: 24829, actually: 15897
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d::  Test matrix 2 nonzeros: 24829, actually: 15897
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d:cg::Starting value 22.
+DEAL:2d:cg::Convergence step 38 value 2.3e-07
+DEAL:2d::Distribute error: 0
+DEAL:2d::
+DEAL:2d::
+DEAL:2d::   Number of active cells:             261
+DEAL:2d::   Number of degrees of freedom:       2021
+DEAL:2d::   Number of hanging node constraints: 327
+DEAL:2d::   Total number of constraints:        580
+DEAL:2d::  Reference matrix nonzeros: 48091, actually: 30641
+DEAL:2d::  Test matrix 1 nonzeros: 48091, actually: 30641
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d::  Test matrix 2 nonzeros: 48091, actually: 30641
+DEAL:2d::  Matrix difference norm: 0
+DEAL:2d::  RHS difference norm: 0
+DEAL:2d:cg::Starting value 30.
+DEAL:2d:cg::Convergence step 52 value 3.9e-07
+DEAL:2d::Distribute error: 0
+DEAL:3d::
+DEAL:3d::
+DEAL:3d::   Number of active cells:             8
+DEAL:3d::   Number of degrees of freedom:       27
+DEAL:3d::   Number of hanging node constraints: 0
+DEAL:3d::   Total number of constraints:        26
+DEAL:3d::  Reference matrix nonzeros: 343, actually: 27
+DEAL:3d::  Test matrix 1 nonzeros: 343, actually: 27
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d::  Test matrix 2 nonzeros: 343, actually: 27
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d:cg::Starting value 1.4
+DEAL:3d:cg::Convergence step 1 value 0
+DEAL:3d::Distribute error: 0
+DEAL:3d::
+DEAL:3d::
+DEAL:3d::   Number of active cells:             43
+DEAL:3d::   Number of degrees of freedom:       145
+DEAL:3d::   Number of hanging node constraints: 40
+DEAL:3d::   Total number of constraints:        128
+DEAL:3d::  Reference matrix nonzeros: 3663, actually: 209
+DEAL:3d::  Test matrix 1 nonzeros: 3663, actually: 209
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d::  Test matrix 2 nonzeros: 3663, actually: 209
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d:cg::Starting value 5.8
+DEAL:3d:cg::Convergence step 6 value 4.5e-09
+DEAL:3d::Distribute error: 0
+DEAL:3d::
+DEAL:3d::
+DEAL:3d::   Number of active cells:             99
+DEAL:3d::   Number of degrees of freedom:       561
+DEAL:3d::   Number of hanging node constraints: 197
+DEAL:3d::   Total number of constraints:        422
+DEAL:3d::  Reference matrix nonzeros: 20835, actually: 3425
+DEAL:3d::  Test matrix 1 nonzeros: 20835, actually: 3425
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d::  Test matrix 2 nonzeros: 20835, actually: 3425
+DEAL:3d::  Matrix difference norm: 0
+DEAL:3d::  RHS difference norm: 0
+DEAL:3d:cg::Starting value 4.8
+DEAL:3d:cg::Convergence step 10 value 4.2e-08
+DEAL:3d::Distribute error: 0

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.