PLaSK library

Loading...
Searching...
No Matches
matrix.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2023 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK_COMMON_FEM_MATRIX_HPP
15#define PLASK_COMMON_FEM_MATRIX_HPP
16
17#include <plask/plask.hpp>
18
19 namespace plask {
20
21 struct FemMatrix {
22 const size_t rank;
23 const size_t size;
24 double* data;
25 const Solver* solver;
26
27 FemMatrix(const Solver* solver, size_t rank, size_t size)
28 : rank(rank), size(size), data(aligned_malloc<double>(size)), solver(solver) {
29 clear();
30 }
31
32 FemMatrix(const FemMatrix&) = delete;
33
34 virtual ~FemMatrix() { aligned_free<double>(data); }
35
41 virtual double& operator()(size_t r, size_t c) = 0;
42
44 virtual void clear() {
45 std::fill_n(data, size, 0.);
46 }
47
52 virtual void factorize() {}
53
61 virtual void solverhs(DataVector<double>& B, DataVector<double>& X) = 0;
62
69 void solve(DataVector<double>& B, DataVector<double>& X) {
70 factorize();
71 solverhs(B, X);
72 }
73
80 void solve(DataVector<double>& B) {
81 solve(B, B);
82 }
83
89 virtual void mult(const DataVector<const double>& vector, DataVector<double>& result) = 0;
90
96 virtual void addmult(const DataVector<const double>& vector, DataVector<double>& result) = 0;
97
104 virtual void setBC(DataVector<double>& B, size_t r, double val) = 0;
105
111 template <typename BoundaryConditonsT> void applyBC(const BoundaryConditonsT& bconds, DataVector<double>& B) {
112 // boundary conditions of the first kind
113 for (auto cond : bconds) {
114 for (auto r : cond.place) {
115 setBC(B, r, cond.value);
116 }
117 }
118 }
119
120 virtual std::string describe() const {
121 return format("rank={}, size={}", rank, size);
122 }
123};
124
125 struct BandMatrix : FemMatrix {
126 const size_t ld;
127 const size_t kd;
128
129 BandMatrix(const Solver* solver, size_t rank, size_t kd, size_t ld)
130 : FemMatrix(solver, rank, rank * (ld + 1)), ld(ld), kd(kd) {}
131
132 void setBC(DataVector<double>& B, size_t r, double val) override {
133 B[r] = val;
134 (*this)(r, r) = 1.;
135 size_t start = (r > kd) ? r - kd : 0;
136 size_t end = (r + kd < rank) ? r + kd + 1 : rank;
137 for (size_t c = start; c < r; ++c) {
138 B[c] -= (*this)(r, c) * val;
139 (*this)(r, c) = 0.;
140 }
141 for (size_t c = r + 1; c < end; ++c) {
142 B[c] -= (*this)(r, c) * val;
143 (*this)(r, c) = 0.;
144 }
145 }
146
147 std::string describe() const override {
148 return format("rank={}, bands={}, size={}", rank, kd+1, size);
149 }
150};
151
152
153} // namespace plask
154
155#endif