5#include <dune/common/ftraits.hh>
6#include <dune/common/rangeutilities.hh>
7#include <dune/common/shared_ptr.hh>
8#include <dune/functions/common/multiindex.hh>
10#include <amdis/Output.hpp>
11#include <amdis/algorithm/ForEach.hpp>
12#include <amdis/algorithm/InnerProduct.hpp>
13#include <amdis/algorithm/Transform.hpp>
14#include <amdis/common/FakeContainer.hpp>
15#include <amdis/operations/Assigner.hpp>
16#include <amdis/typetree/MultiIndex.hpp>
21 template <
class DofMap>
25 using AssignMode = std::conditional_t<std::is_same_v<A,Assigner::plus_assign>,
26 std::integral_constant<::InsertMode, ADD_VALUES>,
27 std::integral_constant<::InsertMode, INSERT_VALUES>>;
47 : dofMap_(other.dofMap_)
54 : dofMap_(other.dofMap_)
55 , initialized_(other.initialized_)
57 if (other.initialized_) {
58 VecDuplicate(other.vector_, &vector_);
61 Vec local, otherLocal;
62 VecGhostGetLocalForm(vector_, &local);
63 VecGhostGetLocalForm(other.vector_, &otherLocal);
64 VecCopy(otherLocal, local);
65 VecGhostRestoreLocalForm(vector_, &local);
66 VecGhostRestoreLocalForm(other.vector_, &otherLocal);
104 return dofMap_.localSize();
110 return dofMap_.globalSize();
115 template <
class Basis>
116 void init(Basis
const&,
bool clear)
124 info(3,
" create new vector needed {} seconds", t.elapsed());
138 VecAssemblyBegin(vector_);
139 VecAssemblyEnd(vector_);
140 info(3,
" finish vector assembly needed {} seconds", t.elapsed());
148 VecGhostUpdateBegin(vector_, INSERT_VALUES, SCATTER_FORWARD);
149 VecGhostUpdateEnd(vector_, INSERT_VALUES, SCATTER_FORWARD);
150 info(3,
" synchronize vector needed {} seconds", t.elapsed());
155 template <
class MultiIndex>
159 gather(std::array<MultiIndex,1>{idx}, &y);
165 template <
class MultiIndex,
class Assign>
168 assert(initialized_);
169 VecSetValue(vector_, dofMap_.global(flatMultiIndex(dof)), value, AssignMode<Assign>::value);
174 template <
class IndexRange,
class OutputIterator>
175 void gather(IndexRange
const& localInd, OutputIterator buffer)
const
177 forEach(localInd, [&buffer](
auto&&,
auto const& value) { *buffer = value; ++buffer; });
186 template <
class IndexRange,
class LocalValues,
class MaskRange,
class Assign>
187 void scatter(IndexRange
const& localInd, LocalValues
const& localVal, MaskRange
const& mask, Assign)
189 if constexpr(not std::is_same_v<typename LocalValues::value_type, PetscScalar>)
191 std::vector<PetscScalar> newLocalVal(localVal.begin(), localVal.end());
192 scatter(localInd, newLocalVal, mask, Assign{});
196 assert(initialized_);
197 assert(localInd.size() == localVal.size());
200 std::vector<PetscInt> globalInd;
201 globalInd.reserve(localInd.size());
203 auto ind_it = std::begin(localInd);
204 auto mask_it = std::begin(mask);
205 for (; ind_it != std::end(localInd); ++mask_it, ++ind_it)
206 globalInd.push_back(*mask_it ? dofMap_.global(flatMultiIndex(*ind_it)) : -1);
209 VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
211 VecSetValues(vector_, globalInd.size(), globalInd.data(), localVal.data(), AssignMode<Assign>::value);
212 VecSetOption(vector_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_FALSE);
217 template <
class IndexRange,
class Func>
218 void forEach(IndexRange
const& localInd, Func&& f)
const
221 Vec localGhostedVec =
nullptr;
222 VecGhostGetLocalForm(vector_, &localGhostedVec);
223 assert(localGhostedVec !=
nullptr);
227 VecGetArray(localGhostedVec, &ptr);
229 for (
auto const& dof : localInd)
232 if (dofMap_.owner(i)) {
234 f(dof, ptr[dofMap_.dofToLocal(i)]);
238 f(dof, ptr[dofMap_.localSize() + dofMap_.dofToGhost(i)]);
242 VecRestoreArray(localGhostedVec, &ptr);
243 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
248 template <
class Func>
252 PetscScalar
const* ptr;
253 VecGetArrayRead(vector_, &ptr);
256 f(ptr[dofMap_.dofToLocal(i)]);
258 VecRestoreArrayRead(vector_, &ptr);
264 template <
class Func>
269 VecGetArray(vector_, &ptr);
272 f(ptr[dofMap_.dofToLocal(i)]);
274 VecRestoreArray(vector_, &ptr);
279 template <
class Operation,
class... VecIns>
283 std::array<PetscScalar
const*,
sizeof...(VecIns)> ptrsIn{};
286 VecGetArray(vecOut.vector_, &ptrOut);
287 Ranges::forIndices<
sizeof...(VecIns)>([&](
auto ii) {
288 VecGetArrayRead(std::get<ii>(std::tie(vecIn...)).vector_, &ptrsIn[ii]);
292 auto idx = vecOut.dofMap_.dofToLocal(i);
293 Ranges::applyIndices<
sizeof...(VecIns)>([&](
auto... ii) {
294 ptrOut[idx] = op(ptrsIn[ii][idx]...);
299 VecRestoreArray(vecOut.vector_, &ptrOut);
300 Ranges::forIndices<
sizeof...(VecIns)>([&](
auto ii) {
301 VecRestoreArrayRead(std::get<ii>(std::tie(vecIn...)).vector_, &ptrsIn[ii]);
307 template <
class T,
class BinOp1,
class BinOp2>
310 PetscScalar
const* ptr1;
311 PetscScalar
const* ptr2;
314 VecGetArrayRead(in1.vector_, &ptr1);
315 VecGetArrayRead(in2.vector_, &ptr2);
318 auto idx = in1.dofMap_.dofToLocal(i);
319 init = op1(std::move(
init), op2(ptr1[idx], ptr2[idx]));
323 VecRestoreArrayRead(in1.vector_, &ptr1);
324 VecRestoreArrayRead(in2.vector_, &ptr2);
331 void set(PetscScalar value)
333 Vec localGhostedVec =
nullptr;
334 VecGhostGetLocalForm(vector_, &localGhostedVec);
335 assert(localGhostedVec !=
nullptr);
337 VecSet(localGhostedVec, value);
338 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
345 Vec localGhostedVec =
nullptr;
346 VecGhostGetLocalForm(vector_, &localGhostedVec);
347 assert(localGhostedVec !=
nullptr);
349 VecZeroEntries(localGhostedVec);
350 VecGhostRestoreLocalForm(vector_, &localGhostedVec);
355 assert(&lhs.dofMap_ == &rhs.dofMap_);
356 std::swap(lhs.vector_, rhs.vector_);
357 std::swap(lhs.initialized_, rhs.initialized_);
361 MPI_Comm comm()
const
363 return dofMap_.comm();
370 VecCreateGhost(comm(),
371 dofMap_.localSize(), dofMap_.globalSize(), dofMap_.ghostSize(), dofMap_.ghostIndices().data(),
379 VecDestroy(&vector_);
381 initialized_ =
false;
386 DofMap
const& dofMap_;
389 mutable Vec vector_ =
nullptr;
390 bool initialized_ =
false;
396 template <
class DofMap>
399 template <
class Vec,
class UnaryFunction>
400 static void impl(Vec&& vec, UnaryFunction f)
406 template <
class DofMap>
409 template <
class Operation,
class... VecIns>
416 template <
class DofMap>
419 template <
class T,
class BinOp1,
class BinOp2>
A container-like data-structure not storing anything and with empty implementations in many container...
Definition: FakeContainer.hpp:36
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:23
PetscVector & operator=(PetscVector const &other)
Copy assignment operator.
Definition: VectorBackend.hpp:84
void insert(MultiIndex const &dof, PetscScalar value, Assign)
Access the entry i of the vector with write-access.
Definition: VectorBackend.hpp:166
void forEachLocal(Func &&f)
Perform the f(value) on the local elements of this vector.
Definition: VectorBackend.hpp:265
PetscVector & operator=(PetscVector &&other)
Move assignment operator.
Definition: VectorBackend.hpp:77
::Vec BaseVector
The type of the base vector.
Definition: VectorBackend.hpp:31
void forEachLocal(Func &&f) const
Perform the f(value) on the local elements of this vector.
Definition: VectorBackend.hpp:249
void setZero()
Zero all entries, including the ghost entries.
Definition: VectorBackend.hpp:343
void synchronize()
Update the ghost regions with values from the owning process.
Definition: VectorBackend.hpp:145
void finish()
Finish assembly. Must be called before vector can be used in a KSP.
Definition: VectorBackend.hpp:135
PetscVector(PetscVector &&other)
Move constructor.
Definition: VectorBackend.hpp:46
PetscScalar at(MultiIndex const &idx) const
Access the entry i of the vector with read-only-access.
Definition: VectorBackend.hpp:156
static T innerProductLocal(PetscVector const &in1, PetscVector const &in2, T init, BinOp1 op1, BinOp2 op2)
Perform the op1(init, op2(value1, value2)) on the local elements of [in1] and [in2].
Definition: VectorBackend.hpp:308
PetscVector(PetscVector const &other)
Copy constructor.
Definition: VectorBackend.hpp:53
~PetscVector()
Destructor.
Definition: VectorBackend.hpp:71
BaseVector & vector()
Return the data-vector vector_.
Definition: VectorBackend.hpp:96
std::size_t globalSize() const
Return the number of entries in the global vector.
Definition: VectorBackend.hpp:108
static void transformLocal(PetscVector &vecOut, Operation op, VecIns const &... vecIn)
Perform the valueOut = op(valueIn...) on the local elements of [vecOut] and [vecIn....
Definition: VectorBackend.hpp:280
void init(Basis const &, bool clear)
Resize the vector_ to the size s.
Definition: VectorBackend.hpp:116
void scatter(IndexRange const &localInd, LocalValues const &localVal, MaskRange const &mask, Assign)
Store all values given by localVal with indices localInd in the vector.
Definition: VectorBackend.hpp:187
PetscVector(DofMap const &dofMap)
Constructor. Constructs new BaseVector.
Definition: VectorBackend.hpp:41
std::size_t localSize() const
Return the number of entries in the local part of the vector.
Definition: VectorBackend.hpp:102
PetscInt size_type
The index/size - type.
Definition: VectorBackend.hpp:37
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:90
void set(PetscScalar value)
Set all entries to value, including the ghost entries.
Definition: VectorBackend.hpp:331
void gather(IndexRange const &localInd, OutputIterator buffer) const
Collect all values from this vector to indices given by localInd and store it into the output buffer.
Definition: VectorBackend.hpp:175
PetscScalar value_type
The type of the elements of the DOFVector.
Definition: VectorBackend.hpp:34
void forEach(IndexRange const &localInd, Func &&f) const
Apply the functor f to each vector entry to local index in localInd range.
Definition: VectorBackend.hpp:218
constexpr bool MultiIndex
A multi-index type.
Definition: Concepts.hpp:149
Definition: ForEach.hpp:39
General implementation of recursive inner-product.
Definition: InnerProduct.hpp:58