From: bangerth Date: Tue, 7 Aug 2007 17:58:47 +0000 (+0000) Subject: Comment a bit further X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f2431f0fdadd303e366f39a65673c81bb4ea3288;p=dealii-svn.git Comment a bit further git-svn-id: https://svn.dealii.org/trunk@14916 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox index d0c2c8be48..7ea6fbd59c 100644 --- a/deal.II/examples/step-27/doc/intro.dox +++ b/deal.II/examples/step-27/doc/intro.dox @@ -238,7 +238,7 @@ i.e. the interesting part of the loop over all cells would look like this: const FEValues &fe_values = hp_fe_values.get_present_fe_values (); ... // assemble local contributions and copy them into global object - } + } @@ -509,6 +509,8 @@ x)$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$.

Complications with linear systems for hp discretizations

+

Creating the sparsity pattern

+ One of the problems with $hp$ methods is that the high polynomial degree of shape functions together with the large number of constrained degrees of freedom leads to matrices with large numbers of nonzero entries in some @@ -568,3 +570,113 @@ class has a slight advantage for lower order elements in 3d), but for high order and $hp$ elements in 3d, the CompressedSetSparsityPattern has a definite edge. We will therefore use it later when we build the sparsity pattern in this tutorial program. + + +

Eliminating constrained degrees of freedom

+ +A second problem particular to $hp$ methods arises because we have so many +constrained degrees of freedom: typically up to about one third of all degrees +of freedom (in 3d) are constrained because they either belong to cells with +hanging nodes or because they are on cells adjacent to cells with a higher or +lower polynomial degree. This is much more than the fraction of constrained +degrees of freedom in non-$hp$ mode. + +In previous programs, for example in step-6, we have dealt with eliminating +these constrained degrees of freedom by first building an object of type +ConstraintMatrix that stores the constraints, and then "condensing" away these +degrees of freedom first from the sparsity pattern. The same scheme is used +for the matrix and right hand side, which are also both first built then +condensed. + +Dealing with sparsity patterns and matrices in this way turns out to be +inefficient because the effort necessary is at least ${\cal O}(N \log N)$ in +the number of unknowns, whereas an ideal finite element program would of +course only have algorithms that are linear in the number of unknowns. The +solution to this problem is to use the ConstraintMatrix object already at the +time of creating the sparsity pattern, or while copying local contributions +into the global matrix object (or global vector). + +So, instead of the code snippet (taken from @ref step_6 "step-6") + +
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  hanging_node_constraints.condense (sparsity_pattern);
+
+ +we now use the following: + +
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
+                                   hanging_node_constraints);
+
+ + +In a similar vein, we have to slightly modify the way we copy local +contributions into global matrices and vectors. In previous tutorial programs, +we have always used a process like this: + +
+  typename hp::DoFHandler::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      ...  // assemble local contributions them into global object
+
+      cell->get_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i
+
+We now replace copying and later condensing into one step, as already shown in
+@ref step_17 "step-17" and @ref step_18 "step-18":
+
+
+  typename hp::DoFHandler::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      ...  // assemble local contributions them into global object
+
+      cell->get_dof_indices (local_dof_indices);
+
+      hanging_node_constraints
+	.distribute_local_to_global (cell_matrix,
+				     local_dof_indices,
+				     system_matrix);
+	  
+      hanging_node_constraints
+	.distribute_local_to_global (cell_rhs,
+				     local_dof_indices,
+				     system_rhs);
+    }
+

+
+Essentially, what the
+hanging_node_constraints.distribute_local_to_global calls do is
+to implement the same loop as before, but whenever we hit a constrained degree
+of freedom, the function does the right thing and already condenses it away.
+
+Using this technique of eliminating constrained nodes already when
+transferring local contributions into the global objects, we avoid the problem
+of having to go back later and change these objects. Timing these operations
+shows that this makes the overall algorithms faster.
diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc
index b62aa07440..bb6052088d 100644
--- a/deal.II/examples/step-27/step-27.cc
+++ b/deal.II/examples/step-27/step-27.cc
@@ -217,8 +217,14 @@ LaplaceProblem::~LaplaceProblem ()
 				 // directly build the sparsity pattern, but
 				 // first create an intermediate object that
 				 // we later copy into the right data
-				 // structure. This is as explained in the
-				 // introduction of this program.
+				 // structure. In another slight deviation, we
+				 // do not first build the sparsity pattern
+				 // then condense away constrained degrees of
+				 // freedom, but pass the constraint matrix
+				 // object directly to the function that
+				 // builds the sparsity pattern. Both of these
+				 // steps are explained in the introduction of
+				 // this program.
 				 //
 				 // The second change, maybe hidden in plain
 				 // sight, is that the dof_handler variable
@@ -251,6 +257,23 @@ void LaplaceProblem::setup_system ()
 
 
 				 // @sect4{LaplaceProblem::assemble_system}
+
+				 // This is the function that assembles the
+				 // global matrix and right hand side vector
+				 // from the local contributions of each
+				 // cell. Its main working is as has been
+				 // described in many of the tutorial programs
+				 // before. The significant deviations are the
+				 // ones necessary for $hp$ finite element
+				 // methods. In particular, that we need to
+				 // use a collection of FEValues object
+				 // (implemented through the hp::FEValues
+				 // class), and that we have to eliminate
+				 // constrained degrees of freedom already
+				 // when copying local contributions into
+				 // global objects. Both of these are
+				 // explained in detail in the introduction of
+				 // this program.
 template 
 void LaplaceProblem::assemble_system () 
 {
@@ -311,6 +334,11 @@ void LaplaceProblem::assemble_system ()
 				     system_rhs);
     }
 
+				   // After the steps already discussed above,
+				   // all we still have to do is to treat
+				   // Dirichlet boundary values
+				   // correctly. This works in exactly the
+				   // same way as for non-$hp$ objects:
   std::map boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
 					    0,
@@ -322,6 +350,16 @@ void LaplaceProblem::assemble_system ()
 				      system_rhs);
 }
 
+
+
+				 // @sect4{LaplaceProblem::solve}
+
+				 // The function solving the linear system is
+				 // entirely unchanged from previous
+				 // examples. We simply try to reduce the
+				 // initial residual (which equals the $l_2$
+				 // norm of the right hand side) by a certain
+				 // factor:
 template 
 void LaplaceProblem::solve () 
 {
@@ -340,175 +378,7 @@ void LaplaceProblem::solve ()
 
 
 
-template 
-void
-LaplaceProblem::
-estimate_smoothness (Vector &smoothness_indicators) const
-{
-  const unsigned int N = max_degree;
-
-				   // form all the Fourier vectors
-				   // that we want to
-				   // consider. exclude k=0 to avoid
-				   // problems with |k|^{-mu} and also
-				   // logarithms of |k|
-  std::vector > k_vectors;
-  std::vector   k_vectors_magnitude;
-  switch (dim)
-    {
-      case 2:
-      {
-	for (unsigned int i=0; i(deal_II_numbers::PI * i,
-						deal_II_numbers::PI * j));
-		k_vectors_magnitude.push_back (i*i+j*j);
-	      }
-	
-	break;
-      }
-
-      case 3:
-      {
-	for (unsigned int i=0; i(deal_II_numbers::PI * i,
-						  deal_II_numbers::PI * j,
-						  deal_II_numbers::PI * k));
-		  k_vectors_magnitude.push_back (i*i+j*j+k*k);
-	      }
-	
-	break;
-      }
-      
-      default:
-	    Assert (false, ExcNotImplemented());
-    }
-
-  const unsigned n_fourier_modes = k_vectors.size();
-  std::vector ln_k (n_fourier_modes);
-  for (unsigned int i=0; i      base_quadrature (2);
-  QIterated quadrature (base_quadrature, N);
-  
-  std::vector > >
-    fourier_transform_matrices (fe_collection.size());
-  for (unsigned int fe=0; fe sum = 0;
-	    for (unsigned int q=0; q x_q = quadrature.point(q);
-		sum += std::exp(std::complex(0,1) *
-				(k_vectors[k] * x_q)) *
-		       fe_collection[fe].shape_value(i,x_q) *
-		       quadrature.weight(q);
-	      }
-	    fourier_transform_matrices[fe](k,i)
-	      = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
-	  }
-    }
-
-				   // the next thing is to loop over
-				   // all cells and do our work there,
-				   // i.e. to locally do the Fourier
-				   // transform and estimate the decay
-				   // coefficient
-  std::vector > fourier_coefficients (n_fourier_modes);
-  Vector                     local_dof_values;
-  
-  typename hp::DoFHandler::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (unsigned int index=0; cell!=endc; ++cell, ++index)
-    {
-      local_dof_values.reinit (cell->get_fe().dofs_per_cell);
-      cell->get_dof_values (solution, local_dof_values);
-
-				       // first compute the Fourier
-				       // transform of the local
-				       // solution
-      std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
-      for (unsigned int f=0; fget_fe().dofs_per_cell; ++i)
-	  fourier_coefficients[f] += 
-	    fourier_transform_matrices[cell->active_fe_index()](f,i)
-	    *
-	    local_dof_values(i);
-
-				       // enter the Fourier
-				       // coefficients into a map with
-				       // the k-magnitudes, to make
-				       // sure that we get only the
-				       // largest magnitude for each
-				       // value of |k|
-      std::map k_to_max_U_map;
-      for (unsigned int f=0; f
 void LaplaceProblem::postprocess (const unsigned int cycle)
@@ -595,6 +465,7 @@ void LaplaceProblem::postprocess (const unsigned int cycle)
 }
 
 
+				 // @sect4{LaplaceProblem::create_coarse_grid}
 template <>
 void LaplaceProblem<2>::create_coarse_grid ()
 {
@@ -668,6 +539,7 @@ void LaplaceProblem<2>::create_coarse_grid ()
 
 
 
+				 // @sect4{LaplaceProblem::run}
 template 
 void LaplaceProblem::run () 
 {
@@ -700,6 +572,186 @@ void LaplaceProblem::run ()
     }
 }
 
+
+
+template 
+void
+LaplaceProblem::
+estimate_smoothness (Vector &smoothness_indicators) const
+{
+  const unsigned int N = max_degree;
+
+				   // form all the Fourier vectors
+				   // that we want to
+				   // consider. exclude k=0 to avoid
+				   // problems with |k|^{-mu} and also
+				   // logarithms of |k|
+  std::vector > k_vectors;
+  std::vector   k_vectors_magnitude;
+  switch (dim)
+    {
+      case 2:
+      {
+	for (unsigned int i=0; i(deal_II_numbers::PI * i,
+						deal_II_numbers::PI * j));
+		k_vectors_magnitude.push_back (i*i+j*j);
+	      }
+	
+	break;
+      }
+
+      case 3:
+      {
+	for (unsigned int i=0; i(deal_II_numbers::PI * i,
+						  deal_II_numbers::PI * j,
+						  deal_II_numbers::PI * k));
+		  k_vectors_magnitude.push_back (i*i+j*j+k*k);
+	      }
+	
+	break;
+      }
+      
+      default:
+	    Assert (false, ExcNotImplemented());
+    }
+
+  const unsigned n_fourier_modes = k_vectors.size();
+  std::vector ln_k (n_fourier_modes);
+  for (unsigned int i=0; i      base_quadrature (2);
+  QIterated quadrature (base_quadrature, N);
+  
+  std::vector > >
+    fourier_transform_matrices (fe_collection.size());
+  for (unsigned int fe=0; fe sum = 0;
+	    for (unsigned int q=0; q x_q = quadrature.point(q);
+		sum += std::exp(std::complex(0,1) *
+				(k_vectors[k] * x_q)) *
+		       fe_collection[fe].shape_value(i,x_q) *
+		       quadrature.weight(q);
+	      }
+	    fourier_transform_matrices[fe](k,i)
+	      = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
+	  }
+    }
+
+				   // the next thing is to loop over
+				   // all cells and do our work there,
+				   // i.e. to locally do the Fourier
+				   // transform and estimate the decay
+				   // coefficient
+  std::vector > fourier_coefficients (n_fourier_modes);
+  Vector                     local_dof_values;
+  
+  typename hp::DoFHandler::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (unsigned int index=0; cell!=endc; ++cell, ++index)
+    {
+      local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+      cell->get_dof_values (solution, local_dof_values);
+
+				       // first compute the Fourier
+				       // transform of the local
+				       // solution
+      std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
+      for (unsigned int f=0; fget_fe().dofs_per_cell; ++i)
+	  fourier_coefficients[f] += 
+	    fourier_transform_matrices[cell->active_fe_index()](f,i)
+	    *
+	    local_dof_values(i);
+
+				       // enter the Fourier
+				       // coefficients into a map with
+				       // the k-magnitudes, to make
+				       // sure that we get only the
+				       // largest magnitude for each
+				       // value of |k|
+      std::map k_to_max_U_map;
+      for (unsigned int f=0; ftry block and catch whatever
+				 // exceptions are thrown, thereby producing
+				 // meaningful output if anything should go
+				 // wrong:
 int main () 
 {
   try