lac_initializer.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by Luca Heltai
//
// This file is part of the FSI-suite platform, based on the deal.II library.
//
// The FSI-suite platform is free software; you can use it, redistribute it,
// and/or modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3.0 of the License,
// or (at your option) any later version. The full text of the license can be
// found in the file LICENSE at the top level of the FSI-suite platform
// distribution.
//
// ---------------------------------------------------------------------
 
#ifndef fsi_lac_initializer_h
#define fsi_lac_initializer_h
 
// This includes all types we know of.
#include <mpi.h>
 
#include "lac.h"
 
namespace LAC
{
  class BlockInitializer
  {
  public:
    BlockInitializer(
      const std::vector<dealii::types::global_dof_index> &dofs_per_block,
      const std::vector<dealii::IndexSet>                &owned,
      const std::vector<dealii::IndexSet>                &relevant,
      const MPI_Comm                                     &comm = MPI_COMM_WORLD)
      : dofs_per_block(dofs_per_block)
      , owned(owned)
      , relevant(relevant)
      , comm(comm){};
 
    void
    operator()(LATrilinos::BlockVector &v, bool fast = false)
    {
      v.reinit(owned, comm, fast);
    };
 
 
    void
    ghosted(LATrilinos::BlockVector &v, bool fast = false)
    {
      v.reinit(owned, relevant, comm, fast);
    };
 
    void
    operator()(LAPETSc::BlockVector &v, bool fast = false)
    {
      (void)fast;
      v.reinit(owned, comm);
    };
 
 
    void
    ghosted(LAPETSc::BlockVector &v, bool fast = false)
    {
      (void)fast;
      v.reinit(owned, relevant, comm);
    };
 
    void
    operator()(LAdealii::BlockVector &v, bool fast = false)
    {
      v.reinit(dofs_per_block, fast);
    };
 
 
    void
    ghosted(LAdealii::BlockVector &v, bool fast = false)
    {
      v.reinit(dofs_per_block, fast);
    };
 
    template <int dim, int spacedim>
    void
    operator()(dealii::TrilinosWrappers::BlockSparsityPattern     &s,
               const dealii::DoFHandler<dim, spacedim>            &dh,
               const dealii::AffineConstraints<double>            &cm,
               const dealii::Table<2, dealii::DoFTools::Coupling> &coupling)
    {
      s.reinit(owned, owned, relevant, comm);
      dealii::DoFTools::make_sparsity_pattern(
        dh,
        coupling,
        s,
        cm,
        false,
        dealii::Utilities::MPI::this_mpi_process(comm));
      s.compress();
    }
 
    template <int dim, int spacedim>
    void
    operator()(dealii::BlockSparsityPattern                       &s,
               const dealii::DoFHandler<dim, spacedim>            &dh,
               const dealii::AffineConstraints<double>            &cm,
               const dealii::Table<2, dealii::DoFTools::Coupling> &coupling)
    {
      dsp =
        std::make_unique<dealii::BlockDynamicSparsityPattern>(dofs_per_block,
                                                              dofs_per_block);
 
      dealii::DoFTools::make_sparsity_pattern(dh, coupling, *dsp, cm, false);
      dsp->compress();
      s.copy_from(*dsp);
    }
 
    void
    operator()(const LAdealii::BlockSparsityPattern &sparsity,
               LAdealii::BlockSparseMatrix          &matrix)
    {
      matrix.reinit(sparsity);
    };
 
 
    void
    operator()(const LATrilinos::BlockSparsityPattern &sparsity,
               LATrilinos::BlockSparseMatrix          &matrix)
    {
      matrix.reinit(sparsity);
    };
 
 
 
    void
    operator()(const LAPETSc::BlockSparsityPattern &,
               LAPETSc::BlockSparseMatrix &matrix)
    {
      Assert(dsp, dealii::ExcNotInitialized());
      matrix.reinit(owned, *dsp, comm);
    };
 
 
  private:
    std::unique_ptr<dealii::BlockDynamicSparsityPattern> dsp;
 
    const std::vector<dealii::types::global_dof_index> &dofs_per_block;
 
    const std::vector<dealii::IndexSet> &owned;
 
    const std::vector<dealii::IndexSet> &relevant;
 
    const MPI_Comm &comm;
  };
 
  class Initializer
  {
  public:
    Initializer(const dealii::IndexSet &owned_rows,
                const dealii::IndexSet &relevant_rows,
                const MPI_Comm         &comm             = MPI_COMM_WORLD,
                const dealii::IndexSet &owned_columns    = dealii::IndexSet(),
                const dealii::IndexSet &relevant_columns = dealii::IndexSet())
      : owned_rows(owned_rows)
      , relevant_rows(relevant_rows)
      , comm(comm)
      , owned_columns(owned_columns.size() > 0 ? owned_columns : owned_rows)
      , relevant_columns(relevant_columns.size() > 0 ? relevant_columns :
                                                       relevant_rows){};
 
    void
    operator()(LATrilinos::Vector &v, bool fast = false)
    {
      v.reinit(owned_rows, comm, fast);
    };
 
 
    void
    ghosted(LATrilinos::Vector &v, bool fast = false)
    {
      v.reinit(owned_rows, relevant_rows, comm, fast);
    };
 
    void
    operator()(LAPETSc::Vector &v, bool fast = false)
    {
      (void)fast;
      v.reinit(owned_rows, comm);
    };
 
 
    void
    ghosted(LAPETSc::Vector &v, bool fast = false)
    {
      (void)fast;
      v.reinit(owned_rows, relevant_rows, comm);
    };
 
    void
    operator()(LAdealii::Vector &v, bool fast = false)
    {
      v.reinit(owned_rows.size(), fast);
    };
 
 
    void
    ghosted(LAdealii::Vector &v, bool fast = false)
    {
      v.reinit(owned_rows.size(), fast);
    };
 
    void
    operator()(dealii::TrilinosWrappers::SparsityPattern &s)
    {
      s.reinit(owned_rows, owned_columns, comm);
    }
 
    void
    operator()(dealii::SparsityPattern &s)
    {
      s.reinit(owned_rows.size(), owned_columns.size(), 0);
    }
 
    template <int dim, int spacedim>
    void
    operator()(dealii::TrilinosWrappers::SparsityPattern          &s,
               const dealii::DoFHandler<dim, spacedim>            &dh,
               const dealii::AffineConstraints<double>            &cm,
               const dealii::Table<2, dealii::DoFTools::Coupling> &coupling)
    {
      s.reinit(owned_rows, owned_columns, relevant_rows, comm);
      dealii::DoFTools::make_sparsity_pattern(
        dh,
        coupling,
        s,
        cm,
        false,
        dealii::Utilities::MPI::this_mpi_process(comm));
      s.compress();
    }
 
    template <int dim, int spacedim>
    void
    operator()(dealii::SparsityPattern                            &s,
               const dealii::DoFHandler<dim, spacedim>            &dh,
               const dealii::AffineConstraints<double>            &cm,
               const dealii::Table<2, dealii::DoFTools::Coupling> &coupling)
    {
      dsp =
        std::make_unique<dealii::DynamicSparsityPattern>(owned_rows.size(),
                                                         owned_columns.size());
 
      dealii::DoFTools::make_sparsity_pattern(dh, coupling, *dsp, cm, false);
      dsp->compress();
      s.copy_from(*dsp);
    }
 
    void
    operator()(const LAdealii::SparsityPattern &sparsity,
               LAdealii::SparseMatrix          &matrix)
    {
      matrix.reinit(sparsity);
    };
 
 
    void
    operator()(const LATrilinos::SparsityPattern &sparsity,
               LATrilinos::SparseMatrix          &matrix)
    {
      matrix.reinit(sparsity);
    };
 
 
 
    void
    operator()(const LAPETSc::SparsityPattern &sparsity,
               LAPETSc::SparseMatrix          &matrix)
    {
      matrix.reinit(owned_rows, owned_columns, sparsity, comm);
    };
 
 
  private:
    std::unique_ptr<dealii::DynamicSparsityPattern> dsp;
 
    const dealii::IndexSet owned_rows;
    const dealii::IndexSet relevant_rows;
    const MPI_Comm &comm;
 
    const dealii::IndexSet owned_columns;
    const dealii::IndexSet relevant_columns;
  };
} // namespace LAC
#endif