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