Dear Wolfgang,
> It is hard to imagine situations in which the mass matrix would be
singular.
> It is a positive definite form that gives rise to the mass matrix and so
it
> really shouldn't be singular at all. Can you show the code again with
which
> you build it?
It seems that my mesh is neither singular nor degenerate. I wonder if the
problem is that I am solving a vector valued problem. To build this mass
matrix, I adapted the function Diffusion::assemble_system from step-52.
You will find the code attached, as well as the output of one run. If I
isolate the mass matrix built and calculate its determinant in Python, the
result is indeed zero.
I am really at a loss here, so any help is appreciated.
Bests,
Bob
On Mon, 15 Apr 2019 at 16:26, Wolfgang Bangerth <bange...@colostate.edu>
wrote:
> On 4/14/19 11:59 PM, Robert Spartus wrote:
> >
> > Thanks for the insightful discussion on the integrating issue. Wolfgang,
> I
> > guess your last argument is the same as you gave in one of your
> fantastic
> > lectures?
>
> Yes. (Also, thanks for the compliment :-) )
>
>
> > Incidentally, do you have any ideas on how to solve the singularity of
> the
> > mass matrix in this vector-valued problem?
>
> It is hard to imagine situations in which the mass matrix would be
> singular.
> It is a positive definite form that gives rise to the mass matrix and so
> it
> really shouldn't be singular at all. Can you show the code again with
> which
> you build it?
>
> The only situation where the mass matrix may become singular is if you
> have a
> degenerate mesh with elements of zero or negative volume.
>
> Best
> WB
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> 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.
>
--
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
#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_mesh();
void prepare_mass_matrix();
void setup_slab_boundary_conditions();
void setup_system();
template<class Tri>
void prepare_boundary_ids_distributed(Tri& tri) const;
/*! Degree of the approximation polynomial */
int fe_degree;
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 ():
fe_degree(1),
triangulation (MPI_COMM_WORLD, Triangulation<dim, dim >::smoothing_on_refinement, false),
fe (FE_Q<dim>(fe_degree), dim),
dof_handler (triangulation),
quadrature_formula(fe_degree + 1),
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();
}
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));
}
}
//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_mesh() {
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_mesh();
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;
}
