Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library 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 GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : // libMesh includes
20 : #include "libmesh/system_norm.h"
21 : #include "libmesh/enum_norm_type.h"
22 : #include "libmesh/int_range.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 35454 : SystemNorm::SystemNorm() :
28 35454 : _norms(1, DISCRETE_L2), _weights(1, 1.0), _weights_sq(1, 1.0)
29 : {
30 35454 : }
31 :
32 :
33 40293 : SystemNorm::SystemNorm(const FEMNormType & t) :
34 40293 : _norms(1, t), _weights(1, 1.0), _weights_sq(1, 1.0)
35 : {
36 40293 : }
37 :
38 :
39 420 : SystemNorm::SystemNorm(std::vector<FEMNormType> norms) :
40 432 : _norms(std::move(norms)), _weights(1, 1.0), _weights_sq(1, 1.0)
41 : {
42 420 : if (_norms.empty())
43 420 : _norms.push_back(DISCRETE_L2);
44 420 : }
45 :
46 :
47 118084 : SystemNorm::SystemNorm(std::vector<FEMNormType> norms,
48 118084 : std::vector<Real> & weights) :
49 114718 : _norms(std::move(norms)), _weights(weights),
50 121450 : _weights_sq(_weights.size(), 0.0)
51 : {
52 118084 : if (_norms.empty())
53 0 : _norms.push_back(DISCRETE_L2);
54 :
55 118084 : if (_weights.empty())
56 : {
57 0 : _weights.push_back(1.0);
58 0 : _weights_sq.push_back(1.0);
59 : }
60 : else
61 269768 : for (auto i : index_range(_weights))
62 156010 : _weights_sq[i] = _weights[i] * _weights[i];
63 118084 : }
64 :
65 0 : SystemNorm::SystemNorm(std::vector<FEMNormType> norms,
66 0 : std::vector<std::vector<Real>> & weights):
67 0 : _norms(std::move(norms)),
68 0 : _weights(weights.size()),
69 0 : _weights_sq(weights.size()),
70 0 : _off_diagonal_weights(weights)
71 : {
72 0 : if (_norms.empty())
73 0 : _norms.push_back(DISCRETE_L2);
74 :
75 0 : if (_weights.empty())
76 : {
77 0 : _weights.push_back(1.0);
78 0 : _weights_sq.push_back(1.0);
79 : }
80 : else
81 : {
82 : // Loop over the entries of the user provided matrix and store its entries in
83 : // the _off_diagonal_weights or _diagonal_weights
84 0 : for (auto i : index_range(_off_diagonal_weights))
85 : {
86 0 : if (_off_diagonal_weights[i].size() > i)
87 : {
88 0 : _weights[i] = _off_diagonal_weights[i][i];
89 0 : _off_diagonal_weights[i][i] = 0;
90 : }
91 : else
92 0 : _weights[i] = 1.0;
93 : }
94 0 : for (auto i : index_range(_weights))
95 0 : _weights_sq[i] = _weights[i] * _weights[i];
96 : }
97 0 : }
98 :
99 125981 : bool SystemNorm::is_discrete() const
100 : {
101 3968 : libmesh_assert (!_norms.empty());
102 :
103 125981 : if (_norms[0] == DISCRETE_L1 ||
104 240622 : _norms[0] == DISCRETE_L2 ||
105 3644 : _norms[0] == DISCRETE_L_INF)
106 11340 : return true;
107 :
108 3644 : return false;
109 : }
110 :
111 :
112 873667 : FEMNormType SystemNorm::type(unsigned int var) const
113 : {
114 109239 : libmesh_assert (!_norms.empty());
115 :
116 955524 : std::size_t i = (var < _norms.size()) ? var : _norms.size() - 1;
117 :
118 873667 : return _norms[i];
119 : }
120 :
121 :
122 :
123 13396 : void SystemNorm::set_type(unsigned int var, const FEMNormType & t)
124 : {
125 456 : libmesh_assert (!_norms.empty());
126 :
127 13852 : if (var >= _norms.size())
128 8400 : _norms.resize(var+1, t);
129 :
130 13852 : _norms[var] = t;
131 13396 : }
132 :
133 :
134 7755092 : Real SystemNorm::weight(unsigned int var) const
135 : {
136 874399 : libmesh_assert (!_weights.empty());
137 :
138 8629303 : return (var < _weights.size()) ? _weights[var] : 1.0;
139 : }
140 :
141 :
142 5600 : void SystemNorm::set_weight(unsigned int var, Real w)
143 : {
144 160 : libmesh_assert (!_weights.empty());
145 :
146 5760 : if (var >= _weights.size())
147 : {
148 4200 : _weights.resize(var+1, 1.0);
149 4200 : _weights_sq.resize(var+1, 1.0);
150 : }
151 :
152 5760 : _weights[var] = w;
153 5760 : _weights_sq[var] = w*w;
154 5600 : }
155 :
156 16800 : void SystemNorm::set_off_diagonal_weight(unsigned int i,
157 : unsigned int j,
158 : Real w)
159 : {
160 480 : libmesh_assert (!_weights.empty());
161 :
162 17280 : if (i >= _off_diagonal_weights.size())
163 : {
164 5600 : _off_diagonal_weights.resize(i+1);
165 : }
166 :
167 17760 : if (j >= _off_diagonal_weights[i].size())
168 : {
169 16800 : _off_diagonal_weights[i].resize(j+1, 0.);
170 : }
171 :
172 17760 : _off_diagonal_weights[i][j] = w;
173 :
174 16800 : }
175 :
176 :
177 274294 : Real SystemNorm::weight_sq(unsigned int var) const
178 : {
179 21913 : libmesh_assert (!_weights_sq.empty());
180 :
181 295943 : return (var < _weights_sq.size()) ? _weights_sq[var] : 1.0;
182 : }
183 :
184 :
185 58576 : Real SystemNorm::calculate_norm(const std::vector<Real> & v1,
186 : const std::vector<Real> & v2)
187 : {
188 : // The vectors are assumed to both be vectors of the (same number
189 : // of) components
190 3200 : std::size_t vsize = v1.size();
191 1600 : libmesh_assert_equal_to (vsize, v2.size());
192 :
193 : // We'll support implicitly defining weights, but if the user sets
194 : // more weights than he uses then something's probably wrong
195 3200 : std::size_t diagsize = this->_weights.size();
196 1600 : libmesh_assert_greater_equal (vsize, diagsize);
197 :
198 : // Initialize the variable val
199 1600 : Real val = 0.;
200 :
201 : // Loop over all the components of the system with explicit
202 : // weights
203 292880 : for (std::size_t i = 0; i != diagsize; i++)
204 : {
205 247104 : val += this->_weights[i] * v1[i] * v2[i];
206 : }
207 : // Loop over all the components of the system with implicit
208 : // weights
209 58576 : for (std::size_t i = diagsize; i < vsize; i++)
210 : {
211 0 : val += v1[i] * v2[i];
212 : }
213 :
214 : // Loop over the components of the system
215 3200 : std::size_t nrows = this->_off_diagonal_weights.size();
216 1600 : libmesh_assert_less_equal (vsize, nrows);
217 :
218 292880 : for (std::size_t i = 0; i != nrows; i++)
219 : {
220 12800 : std::size_t ncols = this->_off_diagonal_weights[i].size();
221 1112944 : for (std::size_t j=0; j != ncols; j++)
222 : {
223 : // The diagonal weights here were set to zero in the
224 : // constructor.
225 926640 : val += this->_off_diagonal_weights[i][j] * v1[i] * v2[j];
226 : }
227 : }
228 :
229 58576 : return(val);
230 : }
231 :
232 0 : Real SystemNorm::calculate_norm(const std::vector<Real> & v1)
233 : {
234 0 : return this->calculate_norm(v1,v1);
235 : }
236 :
237 2498 : bool SystemNorm::is_identity()
238 : {
239 216 : std::size_t nrows = this->_off_diagonal_weights.size();
240 :
241 : // If any of the off-diagonal elements is not 0, then we are in the non-identity case
242 4248 : for (std::size_t i = 0; i != nrows; i++)
243 : {
244 160 : std::size_t ncols = this->_off_diagonal_weights[i].size();
245 12250 : for (std::size_t j = 0; j != ncols; j++)
246 : {
247 10500 : if (_off_diagonal_weights[i][j] != 0)
248 : {
249 30 : return(false);
250 : }
251 : }
252 : }
253 :
254 : // If any of the diagonal elements is not 1, then we are in the non-identity case
255 156 : nrows = this->_weights.size();
256 3246 : for (std::size_t i = 0; i != nrows; i++)
257 2148 : if (_weights[i] != 1)
258 10 : return(false);
259 :
260 : // If all the off-diagonals elements are 0, and diagonal elements 1, then we are in an identity case
261 68 : return(true);
262 : }
263 :
264 : }
|
