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 : #ifndef LIBMESH_EXACT_SOLUTION_H
19 : #define LIBMESH_EXACT_SOLUTION_H
20 :
21 :
22 : // Local Includes
23 : #include "libmesh/libmesh_common.h" // for Number
24 :
25 : // C++ includes
26 : #include <map>
27 : #include <vector>
28 : #include <memory>
29 : #include <set>
30 :
31 : namespace libMesh
32 : {
33 :
34 :
35 : // Forward Declarations
36 : class Point;
37 : class EquationSystems;
38 : class Parameters;
39 : class Mesh;
40 : template <typename Output> class FunctionBase;
41 : template <typename Output> class FEMFunctionBase;
42 : enum FEMNormType : int;
43 :
44 : // Is there any way to simplify this?
45 : // All we need are Tensor and Gradient. - RHS
46 : template <typename T> class TensorValue;
47 : template <typename T> class VectorValue;
48 : typedef TensorValue<Number> NumberTensorValue;
49 : typedef NumberTensorValue Tensor;
50 : typedef VectorValue<Number> NumberVectorValue;
51 : typedef NumberVectorValue Gradient;
52 :
53 : /**
54 : * This class handles the computation of the L2 and/or H1 error for
55 : * the Systems in the EquationSystems object which is passed to it.
56 : *
57 : * \note For this to be useful, the user must attach at least one, and
58 : * possibly two, functions which can compute the exact solution and
59 : * its derivative for any component of any system. These are the \p
60 : * exact_value and \p exact_deriv functions below.
61 : *
62 : * \author Benjamin S. Kirk
63 : * \author John W. Peterson (modifications for libmesh)
64 : * \date 2004
65 : */
66 606 : class ExactSolution
67 : {
68 :
69 : public:
70 : /**
71 : * Constructor. The ExactSolution object
72 : * must be initialized with an EquationSystems
73 : * object.
74 : */
75 : explicit
76 : ExactSolution (const EquationSystems & es);
77 :
78 : /**
79 : * The copy constructor and copy/move assignment operators are
80 : * deleted. This class has containers of unique_ptrs so it can't be
81 : * default (shallow) copied, and it has a const reference so it
82 : * can't be assigned to after creation.
83 : */
84 : ExactSolution(const ExactSolution &) = delete;
85 : ExactSolution & operator= (const ExactSolution &) = delete;
86 : ExactSolution & operator= (ExactSolution &&) = delete;
87 :
88 : /**
89 : * Move constructor and destructor are defaulted out-of-line (in the
90 : * C file) to play nicely with our forward declarations.
91 : */
92 : ExactSolution(ExactSolution &&);
93 : ~ExactSolution();
94 :
95 : /**
96 : * The user can indicate that elements in certain subdomains should be
97 : * excluded from the error calculation by passing in a set of subdomain
98 : * ids to ignore. By default, all subdomains are considered in the error
99 : * calculation.
100 : */
101 : void set_excluded_subdomains(const std::set<subdomain_id_type> & excluded);
102 :
103 : /**
104 : * Attach function similar to system.h which
105 : * allows the user to attach a second EquationSystems
106 : * object with a reference fine grid solution.
107 : */
108 : void attach_reference_solution (const EquationSystems * es_fine);
109 :
110 : /**
111 : * Clone and attach arbitrary functors which compute the exact
112 : * values of the EquationSystems' solutions at any point.
113 : */
114 : void attach_exact_values (const std::vector<FunctionBase<Number> *> & f);
115 :
116 : /**
117 : * Clone and attach arbitrary functors which compute the exact
118 : * values of the EquationSystems' solutions at any point.
119 : */
120 : void attach_exact_values (const std::vector<FEMFunctionBase<Number> *> & f);
121 :
122 : /**
123 : * Clone and attach an arbitrary functor which computes the exact
124 : * value of the system \p sys_num solution at any point.
125 : */
126 : void attach_exact_value (unsigned int sys_num,
127 : FunctionBase<Number> * f);
128 :
129 : /**
130 : * Clone and attach an arbitrary functor which computes the exact
131 : * value of the system \p sys_num solution at any point.
132 : */
133 : void attach_exact_value (unsigned int sys_num,
134 : FEMFunctionBase<Number> * f);
135 :
136 : /**
137 : * Attach an arbitrary function which computes the exact value of
138 : * the solution at any point.
139 : */
140 : typedef Number (*ValueFunctionPointer)(const Point & p,
141 : const Parameters & Parameters,
142 : const std::string & sys_name,
143 : const std::string & unknown_name);
144 : void attach_exact_value (ValueFunctionPointer fptr);
145 :
146 : /**
147 : * Clone and attach arbitrary functors which compute the exact
148 : * gradients of the EquationSystems' solutions at any point.
149 : */
150 : void attach_exact_derivs (const std::vector<FunctionBase<Gradient> *> & g);
151 :
152 : /**
153 : * Clone and attach arbitrary functors which compute the exact
154 : * gradients of the EquationSystems' solutions at any point.
155 : */
156 : void attach_exact_derivs (const std::vector<FEMFunctionBase<Gradient> *> & g);
157 :
158 : /**
159 : * Clone and attach an arbitrary functor which computes the exact
160 : * gradient of the system \p sys_num solution at any point.
161 : */
162 : void attach_exact_deriv (unsigned int sys_num,
163 : FunctionBase<Gradient> * g);
164 :
165 : /**
166 : * Clone and attach an arbitrary functor which computes the exact
167 : * gradient of the system \p sys_num solution at any point.
168 : */
169 : void attach_exact_deriv (unsigned int sys_num,
170 : FEMFunctionBase<Gradient> * g);
171 :
172 : /**
173 : * Attach an arbitrary function which computes the exact gradient of
174 : * the solution at any point.
175 : */
176 : typedef Gradient (*GradientFunctionPointer)(const Point & p,
177 : const Parameters & parameters,
178 : const std::string & sys_name,
179 : const std::string & unknown_name);
180 : void attach_exact_deriv (GradientFunctionPointer gptr);
181 :
182 : /**
183 : * Clone and attach arbitrary functors which compute the exact
184 : * second derivatives of the EquationSystems' solutions at any point.
185 : */
186 : void attach_exact_hessians (std::vector<FunctionBase<Tensor> *> h);
187 : /**
188 : * Clone and attach arbitrary functors which compute the exact
189 : * second derivatives of the EquationSystems' solutions at any point.
190 : */
191 : void attach_exact_hessians (std::vector<FEMFunctionBase<Tensor> *> h);
192 :
193 : /**
194 : * Clone and attach an arbitrary functor which computes the exact
195 : * second derivatives of the system \p sys_num solution at any point.
196 : */
197 : void attach_exact_hessian (unsigned int sys_num,
198 : FunctionBase<Tensor> * h);
199 :
200 : /**
201 : * Clone and attach an arbitrary functor which computes the exact
202 : * second derivatives of the system \p sys_num solution at any point.
203 : */
204 : void attach_exact_hessian (unsigned int sys_num,
205 : FEMFunctionBase<Tensor> * h);
206 :
207 : /**
208 : * Attach an arbitrary function which computes the exact second
209 : * derivatives of the solution at any point.
210 : */
211 : typedef Tensor (*HessianFunctionPointer)(const Point & p,
212 : const Parameters & parameters,
213 : const std::string & sys_name,
214 : const std::string & unknown_name);
215 : void attach_exact_hessian (HessianFunctionPointer hptr);
216 :
217 : /**
218 : * Increases or decreases the order of the quadrature rule used for numerical
219 : * integration. The default \p extraorder is 1, because properly
220 : * integrating L2 error requires integrating the squares of terms
221 : * with order p+1, and 2p+2 is 1 higher than what we default to
222 : * using for reasonable mass matrix integration.
223 : */
224 0 : void extra_quadrature_order (const int extraorder)
225 0 : { _extra_order = extraorder; }
226 :
227 : /**
228 : * Computes and stores the error in the solution value e = u-u_h,
229 : * the gradient grad(e) = grad(u) - grad(u_h), and possibly the hessian
230 : * grad(grad(e)) = grad(grad(u)) - grad(grad(u_h)). Does not return
231 : * any value. For that you need to call the l2_error(), h1_error()
232 : * or h2_error() functions respectively.
233 : */
234 : void compute_error(std::string_view sys_name,
235 : std::string_view unknown_name);
236 :
237 : /**
238 : * \returns The integrated L2 error for the system \p sys_name for the
239 : * unknown \p unknown_name.
240 : *
241 : * \note No error computations are actually performed, you must call
242 : * \p compute_error() for that.
243 : */
244 : Real l2_error(std::string_view sys_name,
245 : std::string_view unknown_name);
246 :
247 : /**
248 : * \returns The integrated L1 error for the system \p sys_name for
249 : * the unknown \p unknown_name.
250 : *
251 : * \note No error computations are actually performed, you must call
252 : * \p compute_error() for that.
253 : */
254 : Real l1_error(std::string_view sys_name,
255 : std::string_view unknown_name);
256 :
257 : /**
258 : * \returns The L_INF error for the system \p sys_name for
259 : * the unknown \p unknown_name.
260 : *
261 : * \note No error computations are actually performed, you must call
262 : * compute_error() for that.
263 : *
264 : * \note The result (as for the other norms as well) is not exact,
265 : * but an approximation based on the chosen quadrature rule: to
266 : * compute it, we take the max of the absolute value of the error
267 : * over all the quadrature points.
268 : */
269 : Real l_inf_error(std::string_view sys_name,
270 : std::string_view unknown_name);
271 :
272 : /**
273 : * \returns The H1 error for the system \p sys_name for the unknown
274 : * \p unknown_name.
275 : *
276 : * \note No error computations are actually performed, you must call
277 : * \p compute_error() for that.
278 : */
279 : Real h1_error(std::string_view sys_name,
280 : std::string_view unknown_name);
281 :
282 : /**
283 : * \returns The H(curl) error for the system \p sys_name for the
284 : * unknown \p unknown_name.
285 : *
286 : * \note No error computations are actually performed, you must call
287 : * \p compute_error() for that.
288 : *
289 : * \note This is only valid for vector-valued elements. An error is
290 : * thrown if requested for scalar-valued elements.
291 : */
292 : Real hcurl_error(std::string_view sys_name,
293 : std::string_view unknown_name);
294 :
295 : /**
296 : * \returns The H(div) error for the system \p sys_name for the
297 : * unknown \p unknown_name.
298 : *
299 : * \note No error computations are actually performed, you must call
300 : * \p compute_error() for that.
301 : *
302 : * \note This is only valid for vector-valued elements. An error is
303 : * thrown if requested for scalar-valued elements.
304 : */
305 : Real hdiv_error(std::string_view sys_name,
306 : std::string_view unknown_name);
307 :
308 : /**
309 : * \returns The H2 error for the system \p sys_name for the unknown
310 : * \p unknown_name.
311 : *
312 : * \note No error computations are actually performed, you must call
313 : * \p compute_error() for that.
314 : */
315 : Real h2_error(std::string_view sys_name,
316 : std::string_view unknown_name);
317 :
318 : /**
319 : * \returns The error in the requested norm for the system \p
320 : * sys_name for the unknown \p unknown_name.
321 : *
322 : * \note No error computations are actually performed, you must call
323 : * \p compute_error() for that.
324 : *
325 : * \note The result is not exact, but an approximation based on the
326 : * chosen quadrature rule.
327 : */
328 : Real error_norm(std::string_view sys_name,
329 : std::string_view unknown_name,
330 : const FEMNormType & norm);
331 : private:
332 :
333 : /**
334 : * This function computes the error (in the solution and its first
335 : * derivative) for a single unknown in a single system. It is a
336 : * private function since it is used by the implementation when
337 : * solving for several unknowns in several systems.
338 : */
339 : template<typename OutputShape>
340 : void _compute_error(std::string_view sys_name,
341 : std::string_view unknown_name,
342 : std::vector<Real> & error_vals);
343 :
344 : /**
345 : * This function is responsible for checking the validity of the \p
346 : * sys_name and \p unknown_name inputs.
347 : *
348 : * \returns A reference to the proper vector for storing the values.
349 : */
350 : std::vector<Real> & _check_inputs(std::string_view sys_name,
351 : std::string_view unknown_name);
352 :
353 : /**
354 : * User-provided functors which compute the exact value of the
355 : * solution for each system.
356 : */
357 : std::vector<std::unique_ptr<FEMFunctionBase<Number>>> _exact_values;
358 :
359 : /**
360 : * User-provided functors which compute the exact derivative of the
361 : * solution for each system.
362 : */
363 : std::vector<std::unique_ptr<FEMFunctionBase<Gradient>>> _exact_derivs;
364 :
365 : /**
366 : * User-provided functors which compute the exact hessians of the
367 : * solution for each system.
368 : */
369 : std::vector<std::unique_ptr<FEMFunctionBase<Tensor>>> _exact_hessians;
370 :
371 : /**
372 : * Data structure which stores the errors:
373 : * ||e|| = ||u - u_h||
374 : * ||grad(e)|| = ||grad(u) - grad(u_h)||
375 : * for each unknown in a single system.
376 : * The name of the unknown is
377 : * the key for the map.
378 : */
379 : typedef std::map<std::string, std::vector<Real>, std::less<>>
380 : SystemErrorMap;
381 :
382 : /**
383 : * A map of SystemErrorMaps, which contains entries
384 : * for each system in the EquationSystems object.
385 : * This is required, since it is possible for two
386 : * systems to have unknowns with the *same name*.
387 : */
388 : std::map<std::string, SystemErrorMap, std::less<>> _errors;
389 :
390 : /**
391 : * Constant reference to the \p EquationSystems object
392 : * used for the simulation.
393 : */
394 : const EquationSystems & _equation_systems;
395 :
396 : /**
397 : * Constant pointer to the \p EquationSystems object
398 : * containing the fine grid solution.
399 : */
400 : const EquationSystems * _equation_systems_fine;
401 :
402 : /**
403 : * Extra order to use for quadrature rule
404 : */
405 : int _extra_order;
406 :
407 : /**
408 : * Elements in a subdomain from this set are skipped during the
409 : * error computation.
410 : */
411 : std::set<subdomain_id_type> _excluded_subdomains;
412 : };
413 :
414 :
415 :
416 : } // namespace libMesh
417 :
418 :
419 : #endif // LIBMESH_EXACT_SOLUTION_H
|
