AMDiS  0.1 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
9 namespace AMDiS
10 {
11  template <class DofMap>
12  struct Constraints<PetscMatrix<DofMap>>
13  {
16
18
25  template <class BitVector>
26  static void dirichletBC(Matrix& mat, Vector& x, Vector& b, BitVector const& nodes, bool setDiagonal = true)
27  {
29  rows.clear();
30  auto const& dofMap = mat.dofMap_;
31  for (std::size_t i = 0; i < nodes.size(); ++i)
32  if (nodes[i])
33  rows.push_back(dofMap.global(i));
34
35  // test for symmetry of the matrix
36  PetscBool set, flg;
37  MatIsSymmetricKnown(mat.matrix(), &set, &flg);
38
39  if (set == PETSC_TRUE && flg == PETSC_TRUE)
40  MatZeroRowsColumns(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1.0 : 0.0, x.vector(), b.vector());
41  else
42  MatZeroRows(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1.0 : 0.0, x.vector(), b.vector());
43  }
44
45
47
50  template <class BitVector, class Associations>
51  static void periodicBC(Matrix& mat, Vector& x, Vector& b, BitVector const& left, Associations const& left2right, bool setDiagonal = true)
52  {
54  rows.clear();
55  associated_rows.clear();
56  auto const& dofMap = mat.dofMap_;
57  for (PetscInt i = 0; i < PetscInt(left.size()); ++i) {
58  if (left[i]) {
59  // get global row index
60  PetscInt row_local[2] = {i, at(left2right,i)};
61  PetscInt row[2] = {dofMap.global(row_local[0]), dofMap.global(row_local[1])};
62
63  rows.push_back(row[0]);
64  associated_rows.push_back(row[1]);
65
66  // copy row to associated row
67  PetscInt ncols;
68  PetscInt const* cols;
69  PetscScalar const* vals;
70  MatGetRow(mat.matrix(), row[0], &ncols, &cols, &vals);
71  MatSetValues(mat.matrix(), 1, &row[1], ncols, cols, vals, ADD_VALUES);
72  MatRestoreRow(mat.matrix(), row[0], &ncols, &cols, &vals);
73  }
74  }
75
76  MatAssemblyBegin(mat.matrix(), MAT_FLUSH_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
77  MatAssemblyEnd(mat.matrix(), MAT_FLUSH_ASSEMBLY);
78
79  // clear the periodic rows and add a 1 on the diagonal
80  MatSetOption(mat.matrix(), MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
81  MatZeroRows(mat.matrix(), rows.size(), rows.data(), setDiagonal ? 1 : 0, PETSC_NULL, PETSC_NULL);
82
83  // Insert the -1 value at the off-diagonal associated columns
84  // NOTE: the following is probably very slow
85  for (std::size_t i = 0; i < rows.size(); ++i)
86  MatSetValue(mat.matrix(), rows[i], associated_rows[i], -1, ADD_VALUES);
87  MatAssemblyBegin(mat.matrix(), MAT_FINAL_ASSEMBLY); // Maybe MAT_FLUSH_ASSEMBLY
88  MatAssemblyEnd(mat.matrix(), MAT_FINAL_ASSEMBLY);
89
90  std::vector<PetscScalar> data(rows.size());
91
92  // copy periodic rhs values to associated rows
93  VecGetValues(b.vector(), rows.size(), rows.data(), data.data());
94  VecSetValues(b.vector(), rows.size(), associated_rows.data(), data.data(), ADD_VALUES);
95  VecAssemblyBegin(b.vector());
96  VecAssemblyEnd(b.vector());
97
98  // set the periodic rhs components to 0
99  data.assign(rows.size(), 0);
100  VecSetValues(b.vector(), rows.size(), rows.data(), data.data(), INSERT_VALUES);
101  VecAssemblyBegin(b.vector());
102  VecAssemblyEnd(b.vector());
103
104  // symmetrize the solution vector
105  VecGetValues(x.vector(), rows.size(), rows.data(), data.data());
106  VecSetValues(x.vector(), rows.size(), associated_rows.data(), data.data(), INSERT_VALUES);
107  VecAssemblyBegin(x.vector());
108  VecAssemblyEnd(x.vector());
109  }
110
111  private:
112  template <class Associations>
113  static PetscInt at(Associations const& m, std::size_t idx)
114  {
115  auto it = m.find(idx);
116  assert(it != m.end());
117  return it->second;
118  }
119  };
120
121 } // end namespace AMDiS
BaseVector const & vector() const
Return the data-vector vector_.
Definition: VectorBackend.hpp:81
Contains all classes needed for solving linear and non linear equation systems.
static void periodicBC(Matrix &mat, Vector &x, Vector &b, BitVector const &left, Associations const &left2right, bool setDiagonal=true)
Periodic boundary conditions.
Definition: Constraints.hpp:51
The basic container that stores a base matrix.
Definition: MatrixBackend.hpp:20
BaseMatrix & matrix()
Return a reference to the data-matrix matrix.
Definition: MatrixBackend.hpp:53
Definition: Constraints.hpp:13
static void dirichletBC(Matrix &mat, Vector &x, Vector &b, BitVector const &nodes, bool setDiagonal=true)
Dirichlet boundary conditions.
Definition: Constraints.hpp:26
The basic container that stores a base vector data.
Definition: VectorBackend.hpp:19