AMDiS 2.10
The Adaptive Multi-Dimensional Simulation Toolbox
Constraints.hpp
1#pragma once
2
3#include <vector>
4
5#include <amdis/linearalgebra/Constraints.hpp>
6#include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
7#include <amdis/linearalgebra/petsc/VectorBackend.hpp>
8
9namespace AMDiS
10{
11 template <class DofMap>
12 struct Constraints<PetscMatrix<DofMap>>
13 {
17
18
20
23 template <class BitVector, class Associations>
24 static void periodicBC(Matrix& mat, VectorX& x, VectorB& b, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
25 {
26 thread_local std::vector<PetscInt> rows, associated_rows;
27 rows.clear();
28 associated_rows.clear();
29 auto const& dofMap = mat.rowDofMap_;
30 for (PetscInt i = 0; i < PetscInt(left.size()); ++i) {
31 if (left[i]) {
32 // get global row index
33 PetscInt row_local[2] = {i, PetscInt(left2right.at(i))};
34 PetscInt row[2] = {dofMap.global(row_local[0]), dofMap.global(row_local[1])};
35
36 rows.push_back(row[0]);
37 associated_rows.push_back(row[1]);
38
39 // copy row to associated row
40 PetscInt ncols;
41 PetscInt const* cols;
42 PetscScalar const* vals;
43 MatGetRow(mat.matrix(), row[0], &ncols, &cols, &vals);
44 MatSetValues(mat.matrix(), 1, &row[1], ncols, cols, vals, ADD_VALUES);
45 MatRestoreRow(mat.matrix(), row[0], &ncols, &cols, &vals);
46 }
47 }
48
49 MatAssemblyBegin(mat.matrix(), MAT_FLUSH_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
50 MatAssemblyEnd(mat.matrix(), MAT_FLUSH_ASSEMBLY);
51
52 // clear the periodic rows and add a 1 on the diagonal
53 MatSetOption(mat.matrix(), MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
54 MatZeroRows(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1 : 0, PETSC_NULL, PETSC_NULL);
55
56 // Insert the -1 value at the off-diagonal associated columns
57 // NOTE: the following is probably very slow
58 for (std::size_t i = 0; i < rows.size(); ++i)
59 MatSetValue(mat.matrix(), rows[i], associated_rows[i], -1, ADD_VALUES);
60 MatAssemblyBegin(mat.matrix(), MAT_FINAL_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
61 MatAssemblyEnd(mat.matrix(), MAT_FINAL_ASSEMBLY);
62
63 std::vector<PetscScalar> data(rows.size());
64
65 // copy periodic rhs values to associated rows
66 VecGetValues(b.vector(), rows.size(), rows.data(), data.data());
67 VecSetValues(b.vector(), rows.size(), associated_rows.data(), data.data(), ADD_VALUES);
68 VecAssemblyBegin(b.vector());
69 VecAssemblyEnd(b.vector());
70
71 // set the periodic rhs components to 0
72 data.assign(rows.size(), 0);
73 VecSetValues(b.vector(), rows.size(), rows.data(), data.data(), INSERT_VALUES);
74 VecAssemblyBegin(b.vector());
75 VecAssemblyEnd(b.vector());
76
77 // symmetrize the solution vector
78 VecGetValues(x.vector(), rows.size(), rows.data(), data.data());
79 VecSetValues(x.vector(), rows.size(), associated_rows.data(), data.data(), INSERT_VALUES);
80 VecAssemblyBegin(x.vector());
81 VecAssemblyEnd(x.vector());
82 }
83 };
84
85} // end namespace AMDiS
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:24
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:56
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:23
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:90
static void periodicBC(Matrix &mat, VectorX &x, VectorB &b, BitVector const &left, Associations const &left2right, bool setDiagonal=true)
Periodic boundary conditions.
Definition: Constraints.hpp:24
Definition: Constraints.hpp:14