[deal.II] Trilinos and Neumann boundary

[Markus Buerg]

Hello,

I am trying to solve Laplace with Neumann boundary conditions (step-11). I 
implemented it the following way:

//includes and class definitions

BoundaryValues::BoundaryValues () : Function<2> () {
}

double BoundaryValues::value (const Point<2> &p, const unsigned int) const {
   return 1;
}

double RightHandSide::value (const Point<2> &p, const unsigned int) const {
   return -2;
}

int main (int argc, char *argv[]) {
   Utilities::System::MPI_InitFinalize mpi_initialization (argc, argv);
   LaplaceProblem problem;

   problem.run ();
   return 0;
}

LaplaceProblem::LaplaceProblem () : fe (2), dof_handler (triangulation), 
trilinos_communicator (Utilities::Trilinos::comm_world()) {
}

RightHandSide::RightHandSide () : Function<2> () {
}

void LaplaceProblem::assemble_system () {
//define variables

   for (DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active 
(); cell != dof_handler.end (); ++cell)
      if (cell->subdomain_id () == Utilities::Trilinos::get_this_mpi_process 
(trilinos_communicator)) {
         cell_matrix = 0;
         cell_rhs = 0;
         fe_values.reinit (cell);

         for (unsigned int q_point = 0; q_point < quadrature.size (); 
++q_point) {
            //assemble local_matrix

               cell_rhs (i) += JxW * right_hand_side.value 
(fe_values.quadrature_point (q_point)) * fe_values.shape_value (i, q_point);
            }
         }

         if (cell->at_boundary ())
            for (unsigned int face = 0; face < 
GeometryInfo<2>::faces_per_cell; ++face)
               if (cell->at_boundary (face)) {
                  fe_face_values.reinit (cell, face);

                  for (unsigned int q_point = 0; q_point < 
face_quadrature.size (); ++q_point) {
                     JxW = fe_face_values.JxW (q_point);

                     for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
                        cell_rhs (i) += JxW * boundary_values.value 
(fe_face_values.quadrature_point (q_point)) * fe_face_values.shape_value (i, 
q_point);
                  }
               }

         cell->get_dof_indices (cell_dof_indices);
         constraints.distribute_local_to_global (cell_matrix, 
cell_dof_indices, system_matrix);
         constraints.distribute_local_to_global (cell_rhs, cell_dof_indices, 
system_rhs);
      }
}

void LaplaceProblem::create_grid () {
   GridGenerator::hyper_ball (triangulation);

   {
      static const HyperBallBoundary<2> boundary;

      triangulation.set_boundary (0, boundary);
   }

   triangulation.refine_global (5);
   GridTools::partition_triangulation 
(Utilities::Trilinos::get_n_mpi_processes (trilinos_communicator), 
triangulation);
}

void LaplaceProblem::estimate_error () {
//estimator
}

void LaplaceProblem::refine_grid () {
//p-refinement
}

void LaplaceProblem::run () {
   create_grid ();
   setup_system ();
   assemble_system ();
   solve ();
   estimated_error_per_cell.reinit (Epetra_Map (triangulation.n_active_cells 
(), GridTools::count_cells_with_subdomain_association (triangulation, 
Utilities::Trilinos::get_this_mpi_process (trilinos_communicator)), 0, 
trilinos_communicator));
   estimate_error ();

   while (estimated_error_per_cell.l1_norm () > 1e-8) {
      refine_grid ();
      triangulation.prepare_coarsening_and_refinement ();
      triangulation.execute_coarsening_and_refinement ();
      setup_system ();
      assemble_system ();
      solve ();
      estimated_error_per_cell.reinit (Epetra_Map 
(triangulation.n_active_cells (), 
GridTools::count_cells_with_subdomain_association (triangulation, 
Utilities::Trilinos::get_this_mpi_process (trilinos_communicator)), 0, 
trilinos_communicator));
      estimate_error ();
   }
}

void LaplaceProblem::setup_system () {
   dof_handler.distribute_dofs (fe);
   DoFRenumbering::subdomain_wise (dof_handler);
   constraints.clear ();
   DoFTools::make_hanging_node_constraints (dof_handler, constraints);

   {
      std::vector<bool> boundary_dofs (dof_handler.n_dofs (), false);

      DoFTools::extract_boundary_dofs (dof_handler, std::vector<bool> (1, 
true), boundary_dofs);

      const unsigned int first_boundary_dof = std::distance 
(boundary_dofs.begin (), std::find (boundary_dofs.begin (), boundary_dofs.end 
(), true));

      constraints.add_line (first_boundary_dof);

      for (unsigned int i = first_boundary_dof + 1; i < dof_handler.n_dofs (); 
++i)
         if (boundary_dofs[i])
            constraints.add_entry (first_boundary_dof, i, -1);
   }

   constraints.close ();

   Epetra_Map map (dof_handler.n_dofs (), 
DoFTools::count_dofs_with_subdomain_association (dof_handler, 
Utilities::Trilinos::get_this_mpi_process (trilinos_communicator)), 0, 
trilinos_communicator);

   sparsity_pattern.reinit (map);
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern, 
constraints, false);
   sparsity_pattern.compress ();
   system_matrix.reinit (sparsity_pattern);
   solution.reinit (map);
   system_rhs.reinit (map);
}

void LaplaceProblem::solve () {
   SolverControl solver_control (system_rhs.size (), 1e-8 *  
system_rhs.l2_norm ());
   SolverCG<TrilinosWrappers::MPI::Vector> cg (solver_control);
   TrilinosWrappers::PreconditionSSOR preconditioner;

   preconditioner.initialize (system_matrix);
   cg.solve (system_matrix, solution, system_rhs, preconditioner);
   constraints.distribute (solution);
}

Unfortunately the cg-method does not converge. What do I have to change to get 
this program running?

Thanks,
Markus
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii