petscvectorbackend.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_ASSEMBLER_BACKENDS_PETSCVECTORBACKEND_HH
#define DUNE_ASSEMBLER_BACKENDS_PETSCVECTORBACKEND_HH
 
#include <cassert>
#include <vector>
 
#include <dune/assembler/backends/flatvector.hh>
#include <dune/common/typetraits.hh>
 
#include <petscvec.h>
 
namespace Dune::Assembler {
 
class PETScVectorBackend
{
public:
  using Vector = ::Vec;
  using value_type = PetscScalar;
  using size_type = PetscInt;
  using LocalVector = FlatVector<value_type>;
 
public:
  explicit PETScVectorBackend (Vector& vector) noexcept
    : vector_(&vector)
  {}
 
  template <class Basis>
  void resize (const Basis& basis)
  {
    assert(vector_);
 
    if (*vector_ == nullptr)
      VecCreate(PETSC_COMM_SELF, vector_);
 
    size_type n = basis.dimension();
    VecSetSizes(*vector_, n, n);
    VecSetType(*vector_, VECSEQ);
  }
 
  void setZero ()
  {
    VecZeroEntries(*vector_);
  }
 
  void finish ()
  {
    assert(vector_);
    VecAssemblyBegin(*vector_);
    VecAssemblyEnd(*vector_);
  }
 
  template <class LocalView>
  void scatter (const LocalView& localView, const LocalVector& localVector)
  {
    assert(vector_);
 
    indices_.resize(localView.size());
    for (std::size_t i = 0; i < localView.size(); ++i)
      indices_[i] = localView.index(i);
 
    VecSetValues(*vector_, indices_.size(), indices_.data(),
                 localVector.container_data(), ADD_VALUES);
  }
 
 
  struct Setter
  {
    operator value_type() const
    {
      value_type value{};
      VecGetValues(vector, 1, &idx, &value);
      return value;
    }
 
    Setter& operator= (value_type value)
    {
      VecSetValue(vector, idx, value, INSERT_VALUES);
      return *this;
    }
 
    Setter& operator+= (value_type value)
    {
      VecSetValue(vector, idx, value, ADD_VALUES);
      return *this;
    }
 
    Setter& operator-= (value_type value)
    {
      VecSetValue(vector, idx, -value, ADD_VALUES);
      return *this;
    }
 
    Vector& vector;
    size_type idx;
  };
 
  Setter operator[] (size_type idx)
  {
    assert(vector_);
    return Setter{*vector_, idx};
  }
 
  value_type operator[] (size_type idx) const
  {
    value_type value{};
    VecGetValues(*vector_, 1, &idx, &value);
    return value;
  }
 
  PETScVectorBackend& operator= (value_type value)
  {
    assert(vector_);
    VecSet(*vector_, value);
    return *this;
  }
 
private:
  Vector* vector_ = nullptr;
 
  std::vector<size_type> indices_{};
};
 
} // end namespace Dune::Assembler
 
template <>
struct Dune::AutonomousValueType<Dune::Assembler::PETScVectorBackend::Setter>
{
  using type = typename Dune::Assembler::PETScVectorBackend::value_type;
};
 
#endif // DUNE_ASSEMBLER_BACKENDS_PETSCVECTORBACKEND_HH