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