Dear all,

First of all, thanks for the awesome library. I am starting to learn 
deal.II, and it is great! The documentation has been extremely helpful so 
far.

I am trying to solve the time dependent elastic equation, and there we 
naturally have a mass matrix-like for the velocities. To build this mass 
matrix, I adapted the function Diffusion::assemble_system from step-52. 
However, in my problem the mass matrix has zero determinant, therefore, it 
is not invertible. Is there some mistake in the construction of the mass 
matrix, or something I am missing? 

I have sent attached a minimal source code that reproduces this error, the 
output of a single run, and a file containing the mass matrix.

I am running deal.II 9.0.0 installed via candi in Ubuntu 18.04.

Kind regards,
Bob

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
 2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  6.94444e-03  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  5.55556e-02  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  2.77778e-02  
1.11111e-01  1.11111e-01  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  2.77778e-02  
1.11111e-01  1.11111e-01  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  
0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02 
Number of active cells:       4
Number of degrees of freedom: 18
Preparing mass matrix

MASS MATRIX
Mat Object: 1 MPI processes
  type: mpiaij
 2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  6.94444e-03  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  5.55556e-02  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  2.77778e-02  
1.11111e-01  1.11111e-01  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  2.77778e-02  
1.11111e-01  1.11111e-01  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  
0.00000e+00  0.00000e+00  6.94444e-03  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
2.77778e-02  2.77778e-02  0.00000e+00  0.00000e+00  0.00000e+00  6.94444e-03  
0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  5.55556e-02  0.00000e+00  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02  0.00000e+00 
 0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  
0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00  2.77778e-02 

Testing mass matrix
RHS: 
[Proc0 0-17] 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 
1.000e+00 1.000e+00 1.000e+00 1.000e+00 
Solution: 
[Proc0 0-17] inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf 
inf inf 
Process 0 is out!
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/numerics/vector_tools.h>

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <vector>

using namespace dealii;

template <int dim>
class MassProblem {
    public:
        MassProblem();
        ~MassProblem();

        void run ();
    private:

        void create_slab();
        void prepare_mass_matrix();
        void setup_slab_boundary_conditions();
        void setup_system();
        template<class Tri>
        void prepare_boundary_ids_distributed(Tri& tri) const;

        parallel::shared::Triangulation<dim>  triangulation;
        /*! Data structure with the information of our interpolating functions.
        */
        FESystem<dim>                         fe;
        /*! Object describing how DoFs are connected to the interpolating
         * functions. */
        DoFHandler<dim>                       dof_handler;
        /*! Contains the information of how (and where) to perform
         * the integrations in our system */
        QGauss<dim>                           quadrature_formula;
        /*! Describes the location of the non-null elements in our tangential
         * matrix. */
        SparsityPattern                       sparsity_pattern;
        /*! Matrix describing how the DoFs are connected in order to fulfill not
         * only the BCs, but also how to handle with hanging nodes. */
        ConstraintMatrix                      hanging_node_constraints;
        /*! Set of DoFs which are under the responsibility of the current
         * process. */
        IndexSet                              locally_owned_dofs;
        /*! Set of DoFs which are necessary in order to perform calculations with
         * the DoFs of this proces. */
        IndexSet                              locally_relevant_dofs;
        /*! Partition with which process is responsible for which DoF. This is
         * necessary to create the SparsityPattern */
        std::vector<types::global_dof_index>  local_dofs_per_process;

        /*! PETSc matrix containing the mass matrix of our system */
        PETScWrappers::MPI::SparseMatrix  mass_matrix;

        //MPI info
        /*! Holds which MPI Communicator we are using. Normally, it is
         * MPI_COMM_WORLD. */
        MPI_Comm            mpi_communicator;
        /*! Holds the rank of the current process. */
        const unsigned int  n_mpi_processes;
        /*! Details how many processes in total we have. */
        const unsigned int  this_mpi_process;

        /*! Special ostream that will only write on the screen if we
         * are on the root process. */
        ConditionalOStream pcout;

        /*! Size of the box */
        double size;
};

// Constructor
template <int dim>
MassProblem<dim>::MassProblem ():
    triangulation (MPI_COMM_WORLD, Triangulation<dim, dim >::smoothing_on_refinement, false),
    fe (FE_Q<dim>(1), dim),
    dof_handler (triangulation),
    quadrature_formula(2),
    mpi_communicator(MPI_COMM_WORLD),
    n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
    this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
    pcout (std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)),
    size(1.0)
{
}

template <int dim>
MassProblem<dim>::~MassProblem (){
    dof_handler.clear ();
    triangulation.clear();
    std::cout << "Process " << this_mpi_process << " is out!" << std::endl;
}

template <int dim>
void MassProblem<dim>::setup_slab_boundary_conditions(){

    for(int i = 0; i<dim; i++) {
        int face_l = i*2;
        int face_r = i*2+1;
        const FEValuesExtractors::Scalar i_displacement(i);

        VectorTools:: interpolate_boundary_values(dof_handler,
                face_l,
                ZeroFunction<dim>(dim),
                hanging_node_constraints,
                fe.component_mask(i_displacement));
        VectorTools:: interpolate_boundary_values(dof_handler,
                face_r,
                ZeroFunction<dim>(dim),
                hanging_node_constraints,
                fe.component_mask(i_displacement));
    }

/*         { */
/*             const FEValuesExtractors::Scalar i_displacement(dim-1); */
/*             VectorTools:: interpolate_boundary_values(dof_handler, */
/*                     2*(dim-1), */
/*                     ZeroFunction<dim>(dim), */
/*                     hanging_node_constraints); */
/*         } */

}

//Setup
template <int dim>
void MassProblem<dim>::setup_system(){

    dof_handler.clear();
    dof_handler.distribute_dofs(fe);

    locally_owned_dofs = dof_handler.locally_owned_dofs();
    DoFTools::extract_locally_relevant_dofs (dof_handler,locally_relevant_dofs);
    local_dofs_per_process = dof_handler.n_locally_owned_dofs_per_processor();

    hanging_node_constraints.clear ();
    hanging_node_constraints.reinit(locally_relevant_dofs);
    DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints);
    setup_slab_boundary_conditions();
    hanging_node_constraints.close ();

    DynamicSparsityPattern sparsity_pattern (locally_relevant_dofs);
    DoFTools::make_sparsity_pattern (dof_handler,
            sparsity_pattern,
            hanging_node_constraints, /*keep constrained dofs*/ false);
    SparsityTools::distribute_sparsity_pattern (sparsity_pattern,
            local_dofs_per_process,
            mpi_communicator,
            locally_relevant_dofs);
    mass_matrix.reinit (locally_owned_dofs,
            locally_owned_dofs,
            sparsity_pattern,
            mpi_communicator);

    pcout << "Number of active cells:       " << triangulation.n_active_cells() << std::endl;
    pcout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
}

template <int dim>
template <class Tri>
void MassProblem<dim>::prepare_boundary_ids_distributed(Tri& tri) const {

    for(auto cell: tri.cell_iterators()) {
        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f){
            if (cell->face(f)->at_boundary()){

                const Point<dim> face_center = cell->face(f)->center();

                for(unsigned i = 0; i<dim-1; i++) {
                    unsigned int face_l = 2*i;
                    unsigned int face_r = 2*i + 1;
                    if (face_center(i) < 1e-10) {
                        cell->face(f)->set_boundary_id (face_l); //x==0
                    }
                    if (face_center(i) >= size - 1e-10) {
                        cell->face(f)->set_boundary_id (face_r); //x==L
                    }
                }
            }
        }
    }
}

template <int dim>
void MassProblem<dim>::create_slab() {

    Point<dim> lower;
    Point<dim> upper;
    for(unsigned int i = 0; i<dim; i++) {
        upper(i) = size;
    }
    GridGenerator::hyper_rectangle(triangulation,
            lower, upper,
            true);

    triangulation.refine_global(1);
    prepare_boundary_ids_distributed(triangulation);

}


template <int dim>
void MassProblem<dim>::prepare_mass_matrix(){

    pcout << "Preparing mass matrix" << std::endl << std::endl;

    mass_matrix   = 0.;
    FEValues<dim> fe_values (fe, quadrature_formula,
            update_values | update_quadrature_points | update_JxW_values);

    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    unsigned int n_q_points = quadrature_formula.size();


    FullMatrix<double> cell_mass_matrix(dofs_per_cell, dofs_per_cell);

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
    for (auto cell: dof_handler.active_cell_iterators()){
        if(cell->is_locally_owned()){
            cell_mass_matrix = 0.;
            fe_values.reinit(cell);

            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
                for (unsigned int i = 0; i < dofs_per_cell; ++i)
                    for (unsigned int j = 0; j < dofs_per_cell; ++j) {
                        cell_mass_matrix(i, j) +=
                            fe_values.shape_value(i, q_point) *
                            fe_values.shape_value(j, q_point) *
                            fe_values.JxW(q_point);
                    }
            }
            cell->get_dof_indices(local_dof_indices);
            hanging_node_constraints.distribute_local_to_global(cell_mass_matrix,
                    local_dof_indices,
                    mass_matrix);
        }
    }

    mass_matrix.compress(VectorOperation::add);

    pcout << "MASS MATRIX" << std::endl;
    mass_matrix.write_ascii(PETSC_VIEWER_ASCII_DENSE);
}

template <int dim>
void MassProblem<dim>::run (){
    create_slab();
    setup_system();
    prepare_mass_matrix();

    pcout << std::endl;
    pcout << "Testing mass matrix" << std::endl;

    /*! PETSC vector containing some solution to help us test the mass_matrix
    */
    PETScWrappers::MPI::Vector sol;
    sol.reinit(locally_owned_dofs, MPI_COMM_WORLD);
    sol = 0;

    /*! PETSC vector containing some rhs to help us test the mass_matrix
    */
    PETScWrappers::MPI::Vector rhs;
    rhs.reinit(locally_owned_dofs, MPI_COMM_WORLD);
    rhs = 1;

    pcout << "RHS: " << std::endl;
    rhs.print(std::cout);

    SolverControl solver_control;
    PETScWrappers::SolverPreOnly solver(solver_control, mpi_communicator);
    PETScWrappers::PreconditionLU preconditioner(mass_matrix);

    solver.solve(mass_matrix,
            sol,
            rhs,
            preconditioner);

    pcout << "Solution: " << std::endl;
    sol.print(std::cout);
}

/*! Main function */
int main (int argc, char **argv){

    try{
        Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
        MassProblem<2> problem;
        problem.run();
    } catch (std::exception &exc) {

        std::cerr << std::endl << std::endl
            << "------------------------------------------------"
            << std::endl;
        std::cerr << "Exception on processing: " << std::endl
            << exc.what() << std::endl
            << "Aborting!" << std::endl
            << "--------------------------------------------------"
            << std::endl;

        return 1;
    } catch (...){
        std::cerr << "Unknown exception!" << std::endl
            << "Aborting!" << std::endl
            << "--------------------------------------------------"
            << std::endl;
        return 1;
    }

    return 0;
}

Reply via email to