Constraints.hpp 4.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#pragma once

#include <vector>

#include <amdis/linearalgebra/Constraints.hpp>
#include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
#include <amdis/linearalgebra/petsc/VectorBackend.hpp>

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
11
12
  template <class DofMap>
  struct Constraints<PetscMatrix<DofMap>>
13
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
14
15
    using Matrix = PetscMatrix<DofMap>;
    using Vector = PetscVector<DofMap>;
16
17
18
19
20
21
22
23
24
25
26
27
28
29

    /// Dirichlet boundary conditions.
    /**
     * Clear the rows (columns) of the given boundary nodes that are owned by the processor.
     * The columns are only cleared if the matrix is symmetric. When clearing the columns, the scaled
     * solution values are substracted from the rhs.
     *
     * The rhs vector is adapted to the solution values at the dirichlet nodes.
     **/
    template <class BitVector>
    static void dirichletBC(Matrix& mat, Vector& x, Vector& b, BitVector const& nodes, bool setDiagonal = true)
    {
      thread_local std::vector<PetscInt> rows;
      rows.clear();
Praetorius, Simon's avatar
Praetorius, Simon committed
30
      auto const& dofMap = mat.dofMap_;
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
      for (std::size_t i = 0; i < nodes.size(); ++i)
        if (nodes[i])
          rows.push_back(dofMap.global(i));

      // test for symmetry of the matrix
      PetscBool set, flg;
      MatIsSymmetricKnown(mat.matrix(), &set, &flg);

      if (set == PETSC_TRUE && flg == PETSC_TRUE)
        MatZeroRowsColumns(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1.0 : 0.0, x.vector(), b.vector());
      else
        MatZeroRows(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1.0 : 0.0, x.vector(), b.vector());
    }


    /// Periodic boundary conditions
    /**
     * The periodic constraints are not implemented in a symmetric fashion
     **/
    template <class BitVector, class Associations>
    static void periodicBC(Matrix& mat, Vector& x, Vector& b, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
    {
      thread_local std::vector<PetscInt> rows, associated_rows;
      rows.clear();
      associated_rows.clear();
Praetorius, Simon's avatar
Praetorius, Simon committed
56
      auto const& dofMap = mat.dofMap_;
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
      for (PetscInt i = 0; i < PetscInt(left.size()); ++i) {
        if (left[i]) {
          // get global row index
          PetscInt row_local[2] = {i, at(left2right,i)};
          PetscInt row[2] = {dofMap.global(row_local[0]), dofMap.global(row_local[1])};

          rows.push_back(row[0]);
          associated_rows.push_back(row[1]);

          // copy row to associated row
          PetscInt ncols;
          PetscInt const* cols;
          PetscScalar const* vals;
          MatGetRow(mat.matrix(), row[0], &ncols, &cols, &vals);
          MatSetValues(mat.matrix(), 1, &row[1], ncols, cols, vals, ADD_VALUES);
          MatRestoreRow(mat.matrix(), row[0], &ncols, &cols, &vals);
        }
      }

      MatAssemblyBegin(mat.matrix(), MAT_FLUSH_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
      MatAssemblyEnd(mat.matrix(), MAT_FLUSH_ASSEMBLY);

      // clear the periodic rows and add a 1 on the diagonal
      MatSetOption(mat.matrix(), MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
      MatZeroRows(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1 : 0, PETSC_NULL, PETSC_NULL);

      // Insert the -1 value at the off-diagonal associated columns
      // NOTE: the following is probably very slow
      for (std::size_t i = 0; i < rows.size(); ++i)
        MatSetValue(mat.matrix(), rows[i], associated_rows[i], -1, ADD_VALUES);
      MatAssemblyBegin(mat.matrix(), MAT_FINAL_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
      MatAssemblyEnd(mat.matrix(), MAT_FINAL_ASSEMBLY);

      std::vector<PetscScalar> data(rows.size());

      // copy periodic rhs values to associated rows
      VecGetValues(b.vector(), rows.size(), rows.data(), data.data());
      VecSetValues(b.vector(), rows.size(), associated_rows.data(), data.data(), ADD_VALUES);
      VecAssemblyBegin(b.vector());
      VecAssemblyEnd(b.vector());

      // set the periodic rhs components to 0
      data.assign(rows.size(), 0);
      VecSetValues(b.vector(), rows.size(), rows.data(), data.data(), INSERT_VALUES);
      VecAssemblyBegin(b.vector());
      VecAssemblyEnd(b.vector());

      // symmetrize the solution vector
      VecGetValues(x.vector(), rows.size(), rows.data(), data.data());
      VecSetValues(x.vector(), rows.size(), associated_rows.data(), data.data(), INSERT_VALUES);
      VecAssemblyBegin(x.vector());
      VecAssemblyEnd(x.vector());
    }

  private:
    template <class Associations>
    static PetscInt at(Associations const& m, std::size_t idx)
    {
      auto it = m.find(idx);
      assert(it != m.end());
      return it->second;
    }
  };

} // end namespace AMDiS