DUNE-FEM (unstable)

lfematrix.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
6#define DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
7
8#include <cassert>
9#include <vector>
10
11#include "field.hh"
12
13namespace Dune
14{
15
16 template< class F >
17 class [[deprecated("The usage of LFEMatrix is discouraged. Use DynamicMatrix instead")]] LFEMatrix
18 {
19 typedef std::vector< F > Row;
20 typedef std::vector<Row> RealMatrix;
21
22 public:
23 typedef F Field;
24
25 operator const RealMatrix & () const
26 {
27 return matrix_;
28 }
29
30 operator RealMatrix & ()
31 {
32 return matrix_;
33 }
34
35 template <class Vector>
36 void row( const unsigned int row, Vector &vec ) const
37 {
38 assert(row<rows());
39 for (int i=0; i<cols(); ++i)
40 field_cast(matrix_[row][i], vec[i]);
41 }
42
43 const Field &operator() ( const unsigned int row, const unsigned int col ) const
44 {
45 assert(row<rows());
46 assert(col<cols());
47 return matrix_[ row ][ col ];
48 }
49
50 Field &operator() ( const unsigned int row, const unsigned int col )
51 {
52 assert(row<rows());
53 assert(col<cols());
54 return matrix_[ row ][ col ];
55 }
56
57 unsigned int rows () const
58 {
59 return rows_;
60 }
61
62 unsigned int cols () const
63 {
64 return cols_;
65 }
66
67 const Field *rowPtr ( const unsigned int row ) const
68 {
69 assert(row<rows());
70 return &(matrix_[row][0]);
71 }
72
73 Field *rowPtr ( const unsigned int row )
74 {
75 assert(row<rows());
76 return &(matrix_[row][0]);
77 }
78
79 void resize ( const unsigned int rows, const unsigned int cols )
80 {
81 matrix_.resize(rows);
82 for (unsigned int i=0; i<rows; ++i)
83 matrix_[i].resize(cols);
84 rows_ = rows;
85 cols_ = cols;
86 }
87
88 bool invert ()
89 {
90 using std::abs;
91 assert( rows() == cols() );
92 std::vector<unsigned int> p(rows());
93 for (unsigned int j=0; j<rows(); ++j)
94 p[j] = j;
95 for (unsigned int j=0; j<rows(); ++j)
96 {
97 // pivot search
98 unsigned int r = j;
99 Field max = abs( (*this)(j,j) );
100 for (unsigned int i=j+1; i<rows(); ++i)
101 {
102 if ( abs( (*this)(i,j) ) > max )
103 {
104 max = abs( (*this)(i,j) );
105 r = i;
106 }
107 }
108 if (max == Zero<Field>())
109 return false;
110 // row swap
111 if (r > j)
112 {
113 for (unsigned int k=0; k<cols(); ++k)
114 std::swap( (*this)(j,k), (*this)(r,k) );
115 std::swap( p[j], p[r] );
116 }
117 // transformation
118 Field hr = Unity<Field>()/(*this)(j,j);
119 for (unsigned int i=0; i<rows(); ++i)
120 (*this)(i,j) *= hr;
121 (*this)(j,j) = hr;
122 for (unsigned int k=0; k<cols(); ++k)
123 {
124 if (k==j) continue;
125 for (unsigned int i=0; i<rows(); ++i)
126 {
127 if (i==j) continue;
128 (*this)(i,k) -= (*this)(i,j)*(*this)(j,k);
129 }
130 (*this)(j,k) *= -hr;
131 }
132 }
133 // column exchange
134 Row hv(rows());
135 for (unsigned int i=0; i<rows(); ++i)
136 {
137 for (unsigned int k=0; k<rows(); ++k)
138 hv[ p[k] ] = (*this)(i,k);
139 for (unsigned int k=0; k<rows(); ++k)
140 (*this)(i,k) = hv[k];
141 }
142 return true;
143 }
144
145 private:
146 RealMatrix matrix_;
147 unsigned int cols_,rows_;
148 };
149
150 template< class Field >
151 inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field> &mat)
152 {
153 for (unsigned int r=0; r<mat.rows(); ++r)
154 {
155 out << field_cast<double>(mat(r,0));
156 for (unsigned int c=1; c<mat.cols(); ++c)
157 {
158 out << " , " << field_cast<double>(mat(r,c));
159 }
160 out << std::endl;
161 }
162 return out;
163 }
164}
165
166#endif // #ifndef DUNE_LOCALFUNCTIONS_UTILITY_LFEMATRIX_HH
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Dune namespace.
Definition: alignedallocator.hh:13
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)