]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use a different approach: Build a matrix that contains the vector Laplacian with...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 12 Aug 2008 19:15:57 +0000 (19:15 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 12 Aug 2008 19:15:57 +0000 (19:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@16523 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-31/step-31.cc

index addc5e237f46aab3ebc84714a2957426ef2f4ad3..65af9ab698eb192e90be4c3f92a691cb6b501c43 100644 (file)
@@ -113,6 +113,7 @@ class BoussinesqFlowProblem
 
   private:
     void setup_dofs (const bool setup_matrices);
+    void assemble_preconditioner ();
     void assemble_system ();
     void assemble_rhs_T ();
     double get_maximal_velocity () const;
@@ -130,6 +131,7 @@ class BoussinesqFlowProblem
 
     BlockSparsityPattern      sparsity_pattern;
     BlockSparseMatrix<double> system_matrix;
+    BlockSparseMatrix<double> preconditioner_matrix;
 
     double time_step;
     unsigned int timestep_number;
@@ -733,6 +735,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
   if (setup_matrices == true)
     {
       system_matrix.clear ();
+      preconditioner_matrix.clear ();
 
       BlockCompressedSetSparsityPattern csp (3,3);
  
@@ -747,12 +750,15 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
       csp.block(2,2).reinit (n_T, n_T);
       
       csp.collect_sizes ();
-      
+
+                                      // TODO: only build entries
+                                      // that we really need
       DoFTools::make_sparsity_pattern (dof_handler, csp);
       hanging_node_constraints.condense (csp);
       sparsity_pattern.copy_from (csp);
 
       system_matrix.reinit (sparsity_pattern);
+      preconditioner_matrix.reinit (sparsity_pattern);
     }
 
                                 // As last action in this function,
@@ -785,6 +791,74 @@ void BoussinesqFlowProblem<dim>::setup_dofs (const bool setup_matrices)
 
 
 
+template <int dim>
+double scalar_product (const Tensor<2,dim> &t1,
+                      const Tensor<2,dim> &t2)
+{
+  double s = 0;
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      s += t1[i][j] * t2[i][j];
+  return s;
+}
+
+
+template <int dim>
+void
+BoussinesqFlowProblem<dim>::assemble_preconditioner ()
+{
+  preconditioner_matrix = 0;
+
+  QGauss<dim>   quadrature_formula(degree+2);
+  QGauss<dim-1> face_quadrature_formula(degree+2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                          update_JxW_values |
+                          update_gradients);
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+
+  const unsigned int   n_q_points      = quadrature_formula.size();
+
+  FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+  std::vector<Tensor<2,dim> > phi_grad_u (dofs_per_cell);
+
+  const FEValuesExtractors::Vector velocities (0);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      local_matrix = 0;
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+       {
+         for (unsigned int k=0; k<dofs_per_cell; ++k)
+           phi_grad_u[k] = fe_values[velocities].gradient(k,q);
+
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             local_matrix(i,j) += scalar_product (phi_grad_u[i], phi_grad_u[j])
+                                  * fe_values.JxW(q);
+       }
+
+      cell->get_dof_indices (local_dof_indices);
+
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       for (unsigned int j=0; j<dofs_per_cell; ++j)
+         preconditioner_matrix.add (local_dof_indices[i],
+                                    local_dof_indices[j],
+                                    local_matrix(i,j));
+    }
+  
+  hanging_node_constraints.condense (preconditioner_matrix);
+}
+
+
+
                                 // @sect4{BoussinesqFlowProblem::assemble_system}
                                 // 
                                 // The assembly of the Boussinesq 
@@ -1173,10 +1247,12 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
       std::cout << "   Rebuilding preconditioner..." << std::flush;
 
-     A_preconditioner
+      assemble_preconditioner ();      
+      
+      A_preconditioner
        = boost::shared_ptr<typename InnerPreconditioner<dim>::type>
                (new typename InnerPreconditioner<dim>::type());
-      A_preconditioner->initialize (system_matrix.block(0,0),
+      A_preconditioner->initialize (preconditioner_matrix.block(0,0),
                typename InnerPreconditioner<dim>::type::AdditionalData());
 
       Mp_preconditioner

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.