FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
error_computer.hpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#pragma once
7
8// includes, FEAT
9#include <kernel/analytic/function.hpp>
10#include <kernel/assembly/asm_traits.hpp>
11#include <kernel/util/dist.hpp>
12#include <kernel/util/tiny_algebra.hpp>
13
14namespace FEAT
15{
16 namespace Assembly
17 {
19 namespace Intern
20 {
21 template<bool h0_, typename AsmTraits_, typename AnaTraits_>
22 struct SecHelperH0
23 {
24 typedef typename AsmTraits_::DataType DataType;
25
26 template<typename FuncEval_, typename LocalVec_>
27 static DataType eval(
28 FuncEval_&,
29 const typename AsmTraits_::TrafoEvalData&,
30 const typename AsmTraits_::SpaceEvalData&,
31 const LocalVec_&,
32 const int
33 )
34 {
35 return DataType(0);
36 }
37 };
38
39 template<typename AsmTraits_, typename AnaTraits_>
40 struct SecHelperH0<true, AsmTraits_, AnaTraits_>
41 {
42 typedef typename AsmTraits_::DataType DataType;
43
44 template<typename FuncEval_, typename LocalVec_>
45 static DataType eval(
46 FuncEval_& func_eval,
47 const typename AsmTraits_::TrafoEvalData& trafo_data,
48 const typename AsmTraits_::SpaceEvalData& space_data,
49 const LocalVec_& lvad,
50 const int num_loc_dofs
51 )
52 {
53 // evaluate function value
54 typename AnaTraits_::ValueType value = func_eval.value(trafo_data.img_point);
55
56 // subtract FE function value
57 for(int i(0); i < num_loc_dofs; ++i)
58 value -= lvad[i] * space_data.phi[i].value;
59
60 // return result
61 return Math::abs(value);
62 }
63 };
64
65 template<bool h1_, typename AsmTraits_, typename AnaTraits_, bool sub_dim_>
66 struct SecHelperH1
67 {
68 typedef typename AsmTraits_::DataType DataType;
69
70 template<typename FuncEval_, typename LocalVec_>
71 static DataType eval(
72 FuncEval_&,
73 const typename AsmTraits_::TrafoEvalData&,
74 const typename AsmTraits_::SpaceEvalData&,
75 const LocalVec_&,
76 const int
77 )
78 {
79 return DataType(0);
80 }
81 };
82
83 // full-dimension variant for H1
84 template<typename AsmTraits_, typename AnaTraits_>
85 struct SecHelperH1<true, AsmTraits_, AnaTraits_, false>
86 {
87 typedef typename AsmTraits_::DataType DataType;
88
89 template<typename FuncEval_, typename LocalVec_>
90 static DataType eval(
91 FuncEval_& func_eval,
92 const typename AsmTraits_::TrafoEvalData& trafo_data,
93 const typename AsmTraits_::SpaceEvalData& space_data,
94 const LocalVec_& lvad,
95 const int num_loc_dofs
96 )
97 {
98 // evaluate function gradient
99 typename AnaTraits_::GradientType grad = func_eval.gradient(trafo_data.img_point);
100
101 // subtract FE function gradient
102 for(int i(0); i < num_loc_dofs; ++i)
103 grad.axpy(-lvad[i], space_data.phi[i].grad);
104
105 // return result
106 return grad.norm_euclid_sqr();
107 }
108 };
109
110 // sub-dimensional variant for H1
111 template<typename AsmTraits_, typename AnaTraits_>
112 struct SecHelperH1<true, AsmTraits_, AnaTraits_, true>
113 {
114 typedef typename AsmTraits_::DataType DataType;
115
116 template<typename FuncEval_, typename LocalVec_>
117 static DataType eval(
118 FuncEval_& func_eval,
119 const typename AsmTraits_::TrafoEvalData& trafo_data,
120 const typename AsmTraits_::SpaceEvalData& space_data,
121 const LocalVec_& lvad,
122 const int num_loc_dofs
123 )
124 {
125 // get the dimensions
126 static constexpr int dom_dim = AsmTraits_::domain_dim;
127
128 // evaluate function gradient
129 typename AnaTraits_::GradientType grad = func_eval.gradient(trafo_data.img_point);
130
131 // transform back onto reference element by applying chain rule
132 Tiny::Vector<DataType, dom_dim> ref_grad;
133 ref_grad.set_vec_mat_mult(grad, trafo_data.jac_mat);
134
135 // subtract FE function reference gradient
136 for(int i(0); i < num_loc_dofs; ++i)
137 ref_grad.axpy(-lvad[i], space_data.phi[i].ref_grad);
138
139 // compute gram matrix of transformation
140 Tiny::Matrix<DataType, dom_dim, dom_dim> gram_mat, gram_inv;
141 gram_mat.set_gram(trafo_data.jac_mat);
142 gram_inv.set_inverse(gram_mat);
143
144 // return result
145 return gram_inv.scalar_product(ref_grad, ref_grad);
146 }
147 };
148
149 template<bool h2_, typename AsmTraits_, typename AnaTraits_>
150 struct SecHelperH2
151 {
152 typedef typename AsmTraits_::DataType DataType;
153
154 template<typename FuncEval_, typename LocalVec_>
155 static DataType eval(
156 FuncEval_&,
157 const typename AsmTraits_::TrafoEvalData&,
158 const typename AsmTraits_::SpaceEvalData&,
159 const LocalVec_&,
160 const int
161 )
162 {
163 return DataType(0);
164 }
165 };
166
167 template<typename AsmTraits_, typename AnaTraits_>
168 struct SecHelperH2<true, AsmTraits_, AnaTraits_>
169 {
170 typedef typename AsmTraits_::DataType DataType;
171
172 template<typename FuncEval_, typename LocalVec_>
173 static DataType eval(
174 FuncEval_& func_eval,
175 const typename AsmTraits_::TrafoEvalData& trafo_data,
176 const typename AsmTraits_::SpaceEvalData& space_data,
177 const LocalVec_& lvad,
178 const int num_loc_dofs
179 )
180 {
181 // evaluate function hessian
182 typename AnaTraits_::HessianType hess = func_eval.hessian(trafo_data.img_point);
183
184 // subtract FE function hessian
185 for(int i(0); i < num_loc_dofs; ++i)
186 hess.axpy(-lvad[i], space_data.phi[i].hess);
187
188 // return result
189 return hess.norm_hessian_sqr();
190 }
191 };
192 } // namespace Intern
194
201 template<typename DataType_>
202 struct ScalarErrorInfo
203 {
205 bool have_h0;
207 bool have_h1;
209 bool have_h2;
211 bool have_l1;
213 bool have_lmax;
214
216 DataType_ norm_h0;
218 DataType_ norm_h1;
220 DataType_ norm_h2;
221
223 DataType_ norm_l1;
225 DataType_ norm_lmax;
226
228 ScalarErrorInfo() :
229 have_h0(false),
230 have_h1(false),
231 have_h2(false),
232 have_l1(false),
233 have_lmax(false),
234 norm_h0(DataType_(0)),
235 norm_h1(DataType_(0)),
236 norm_h2(DataType_(0)),
237 norm_l1(DataType_(0)),
238 norm_lmax(DataType_(0))
239 {
240 }
241
243 template<typename DT2_>
244 ScalarErrorInfo(const ScalarErrorInfo<DT2_>& other) :
245 have_h0(other.have_h0),
246 have_h1(other.have_h1),
247 have_h2(other.have_h2),
248 have_l1(other.have_l1),
249 have_lmax(other.have_lmax),
250 norm_h0(DataType_(other.norm_h0)),
251 norm_h1(DataType_(other.norm_h1)),
252 norm_h2(DataType_(other.norm_h2)),
253 norm_l1(DataType_(other.norm_l1)),
254 norm_lmax(DataType_(other.norm_lmax))
255 {
256 }
257
259 template<typename DT2_>
260 ScalarErrorInfo& operator=(const ScalarErrorInfo<DT2_>& other)
261 {
262 have_h0 = other.have_h0;
263 have_h1 = other.have_h1;
264 have_h2 = other.have_h2;
265 have_l1 = other.have_l1;
266 have_lmax = other.have_lmax;
267 norm_h0 = DataType_(other.norm_h0);
268 norm_h1 = DataType_(other.norm_h1);
269 norm_h2 = DataType_(other.norm_h2);
270 norm_l1 = DataType_(other.norm_l1);
271 norm_lmax = DataType_(other.norm_lmax);
272 }
273
283 void synchronize(const Dist::Comm& comm)
284 {
285 DataType_ verr[4] =
286 {
287 have_h0 ? Math::sqr(norm_h0) : DataType_(0),
288 have_h1 ? Math::sqr(norm_h1) : DataType_(0),
289 have_h2 ? Math::sqr(norm_h2) : DataType_(0),
290 have_l1 ? norm_l1 : DataType_(0)
291 };
292
293 comm.allreduce(verr, verr, std::size_t(4), Dist::op_sum);
294
295 if(have_h0) norm_h0 = Math::sqrt(verr[0]);
296 if(have_h1) norm_h1 = Math::sqrt(verr[1]);
297 if(have_h2) norm_h2 = Math::sqrt(verr[2]);
298 if(have_l1) norm_l1 = verr[3];
299
300 if(have_lmax)
301 {
302 comm.allreduce(&norm_lmax, &norm_lmax, std::size_t(1), Dist::op_max);
303 }
304 }
305
321 String format_string(int precision = 0, std::size_t pad_size = 8u, char pad_char = '.') const
322 {
323 String s;
324 if(have_h0)
325 s += String("H0-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h0, precision) + "\n";
326 if(have_h1)
327 s += String("H1-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h1, precision) + "\n";
328 if(have_h2)
329 s += String("H2-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h2, precision) + "\n";
330 return s;
331 }
332
334 friend std::ostream& operator<<(std::ostream& os, const ScalarErrorInfo& sei)
335 {
336 // print errors
337 if(sei.have_h0)
338 os << "H0-Error: " << stringify_fp_sci(sei.norm_h0) << "\n";
339 if(sei.have_h1)
340 os << "H1-Error: " << stringify_fp_sci(sei.norm_h1) << "\n";
341 if(sei.have_h2)
342 os << "H2-Error: " << stringify_fp_sci(sei.norm_h2) << "\n";
343 return os;
344 }
345 }; // struct ScalarErrorInfo
346
374 template<int max_norm_ = 0, bool sub_dimensional_ = false>
375 class ScalarErrorComputer
376 {
377 static_assert(max_norm_ >= 0, "invalid max_norm_ parameter");
378
379 private:
381
384 static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_mat | TrafoTags::jac_det;
385
389 static constexpr SpaceTags space_config =
390 ((max_norm_ >= 0) ? SpaceTags::value : SpaceTags::none) |
391 ((max_norm_ >= 1) ? (sub_dimensional_ ? SpaceTags::ref_grad : SpaceTags::grad) : SpaceTags::none) |
392 ((max_norm_ >= 2) ? SpaceTags::hess : SpaceTags::none);
394
395 public:
419 template<
420 typename Vector_,
421 typename Function_,
422 typename Space_,
423 typename CubatureFactory_>
424 static ScalarErrorInfo<typename Vector_::DataType> compute(
425 const Vector_& vector,
426 const Function_& function,
427 const Space_& space,
428 const CubatureFactory_& cubature_factory)
429 {
430 // make sure the function has everything we need
431 static_assert(Function_::can_value, "function values are required for H0-Error");
432 static_assert(Function_::can_grad || (max_norm_ < 1), "function gradients are required for H1-Error");
433 static_assert(Function_::can_hess || (max_norm_ < 2), "function hessians are required for H2-Error");
434
435 // validate vector dimensions
436 XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector size");
437
438 // vector type
439 typedef Vector_ VectorType;
440 // analytic function type
441 typedef Function_ FunctionType;
442 // space type
443 typedef Space_ SpaceType;
444 // assembly traits
445 typedef AsmTraits1<typename Vector_::DataType, SpaceType, trafo_config, space_config> AsmTraits;
446 // data type
447 typedef typename AsmTraits::DataType DataType;
448
449 // create the cubature rule
450 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
451
452 // fetch the trafo
453 const typename AsmTraits::TrafoType& trafo = space.get_trafo();
454
455 // define our analytic evaluation traits
456 typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
457
458 // create a function evaluator
459 typename FunctionType::template Evaluator<AnalyticEvalTraits> func_eval(function);
460
461 // create a trafo evaluator
462 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
463
464 // create a space evaluator and evaluation data
465 typename AsmTraits::SpaceEvaluator space_eval(space);
466
467 // create a dof-mapping
468 typename AsmTraits::DofMapping dof_mapping(space);
469
470 // create trafo evaluation data
471 typename AsmTraits::TrafoEvalData trafo_data;
472
473 // create space evaluation data
474 typename AsmTraits::SpaceEvalData space_data;
475
476 // create local vector data
477 typename AsmTraits::template TLocalVector<DataType> lvad;
478
479 // create matrix scatter-axpy
480 typename VectorType::GatherAxpy gather_axpy(vector);
481
482 // initialize result
483 ScalarErrorInfo<DataType> result;
484 result.have_h0 = (max_norm_ >= 0);
485 result.have_h1 = (max_norm_ >= 1);
486 result.have_h2 = (max_norm_ >= 2);
487 result.have_l1 = (max_norm_ >= 0);
488 result.have_lmax = (max_norm_ >= 0);
489
490 // H0/H1/H2 helpers
491 typedef Intern::SecHelperH0<(max_norm_ >= 0), AsmTraits, AnalyticEvalTraits> SecH0;
492 typedef Intern::SecHelperH1<(max_norm_ >= 1), AsmTraits, AnalyticEvalTraits, sub_dimensional_> SecH1;
493 typedef Intern::SecHelperH2<(max_norm_ >= 2), AsmTraits, AnalyticEvalTraits> SecH2;
494
495 // loop over all cells of the mesh
496 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
497 {
498 // format local vector
499 lvad.format();
500
501 // initialize dof-mapping
502 dof_mapping.prepare(cell);
503
504 // incorporate local matrix
505 gather_axpy(lvad, dof_mapping);
506
507 // finish dof-mapping
508 dof_mapping.finish();
509
510 // prepare trafo evaluator
511 trafo_eval.prepare(cell);
512
513 // prepare space evaluator
514 space_eval.prepare(trafo_eval);
515
516 // fetch number of local dofs
517 int num_loc_dofs = space_eval.get_num_local_dofs();
518
519 // loop over all quadrature points and integrate
520 for(int k(0); k < cubature_rule.get_num_points(); ++k)
521 {
522 // compute trafo data
523 trafo_eval(trafo_data, cubature_rule.get_point(k));
524
525 // compute basis function data
526 space_eval(space_data, trafo_data);
527
528 // compute integration weight
529 const DataType omega = trafo_data.jac_det * cubature_rule.get_weight(k);
530
531 // get absolute point error
532 const DataType abs_err = SecH0::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs);
533
534 // update results
535 result.norm_lmax = Math::max(result.norm_lmax, abs_err);
536 result.norm_l1 += omega * abs_err;
537 result.norm_h0 += omega * abs_err * abs_err;
538 result.norm_h1 += omega * SecH1::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs);
539 result.norm_h2 += omega * SecH2::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs);
540
541 // continue with next cubature point
542 }
543
544 // finish evaluators
545 space_eval.finish();
546 trafo_eval.finish();
547
548 // continue with next cell
549 }
550
551 // take square-roots
552 result.norm_h0 = Math::sqrt(result.norm_h0);
553 result.norm_h1 = Math::sqrt(result.norm_h1);
554 result.norm_h2 = Math::sqrt(result.norm_h2);
555
556 // okay, that's it
557 return result;
558 }
559 }; // class ScalarErrorComputer
560
562 namespace Intern
563 {
564 template<bool h0_, typename AsmTraits_, typename AnaTraits_>
565 struct VecHelperH0
566 {
567 typedef typename AsmTraits_::DataType DataType;
568 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
569
570 template<typename FuncEval_, typename LocalVec_>
571 static ResultType eval(
572 FuncEval_&,
573 const typename AsmTraits_::TrafoEvalData&,
574 const typename AsmTraits_::SpaceEvalData&,
575 const LocalVec_&,
576 const int
577 )
578 {
579 return ResultType::null();
580 }
581 };
582
583 template<typename AsmTraits_, typename AnaTraits_>
584 struct VecHelperH0<true, AsmTraits_, AnaTraits_>
585 {
586 typedef typename AsmTraits_::DataType DataType;
587 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
588
589 template<typename FuncEval_, typename LocalVec_>
590 static ResultType eval(
591 FuncEval_& func_eval,
592 const typename AsmTraits_::TrafoEvalData& trafo_data,
593 const typename AsmTraits_::SpaceEvalData& space_data,
594 const LocalVec_& lvad,
595 const int num_loc_dofs
596 )
597 {
598 // evaluate function value
599 typename AnaTraits_::ValueType value = func_eval.value(trafo_data.img_point);
600
601 // subtract FE function value
602 for(int i(0); i < num_loc_dofs; ++i)
603 value.axpy(-space_data.phi[i].value, lvad[i]);
604
605 // compute absolute norm of each component
606 ResultType result;
607 for(int k(0); k < AnaTraits_::image_dim; ++k)
608 result[k] = Math::abs(value[k]);
609
610 return result;
611 }
612 };
613
614 template<bool h1_, typename AsmTraits_, typename AnaTraits_>
615 struct VecHelperH1
616 {
617 typedef typename AsmTraits_::DataType DataType;
618 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
619
620 template<typename FuncEval_, typename LocalVec_>
621 static ResultType eval(
622 FuncEval_&,
623 const typename AsmTraits_::TrafoEvalData&,
624 const typename AsmTraits_::SpaceEvalData&,
625 const LocalVec_&,
626 const int
627 )
628 {
629 return ResultType::null();
630 }
631 };
632
633 template<typename AsmTraits_, typename AnaTraits_>
634 struct VecHelperH1<true, AsmTraits_, AnaTraits_>
635 {
636 typedef typename AsmTraits_::DataType DataType;
637 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
638
639 template<typename FuncEval_, typename LocalVec_>
640 static ResultType eval(
641 FuncEval_& func_eval,
642 const typename AsmTraits_::TrafoEvalData& trafo_data,
643 const typename AsmTraits_::SpaceEvalData& space_data,
644 const LocalVec_& lvad,
645 const int num_loc_dofs
646 )
647 {
648 // evaluate reference function gradient
649 typename AnaTraits_::GradientType grad = func_eval.gradient(trafo_data.img_point);
650
651 // subtract FE function gradient
652 for(int i(0); i < num_loc_dofs; ++i)
653 grad.add_outer_product(lvad[i], space_data.phi[i].grad, -DataType(1));
654
655 // compute norm of each component
656 ResultType result;
657 for(int k(0); k < AnaTraits_::image_dim; ++k)
658 result[k] = grad[k].norm_euclid_sqr();
659
660 return result;
661 }
662 };
663
664 template<bool h2_, typename AsmTraits_, typename AnaTraits_>
665 struct VecHelperH2
666 {
667 typedef typename AsmTraits_::DataType DataType;
668 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
669
670 template<typename FuncEval_, typename LocalVec_>
671 static ResultType eval(
672 FuncEval_&,
673 const typename AsmTraits_::TrafoEvalData&,
674 const typename AsmTraits_::SpaceEvalData&,
675 const LocalVec_&,
676 const int
677 )
678 {
679 return ResultType::null();
680 }
681 };
682
683 template<typename AsmTraits_, typename AnaTraits_>
684 struct VecHelperH2<true, AsmTraits_, AnaTraits_>
685 {
686 typedef typename AsmTraits_::DataType DataType;
687 typedef Tiny::Vector<DataType, AnaTraits_::image_dim> ResultType;
688
689 template<typename FuncEval_, typename LocalVec_>
690 static ResultType eval(
691 FuncEval_& func_eval,
692 const typename AsmTraits_::TrafoEvalData& trafo_data,
693 const typename AsmTraits_::SpaceEvalData& space_data,
694 const LocalVec_& lvad,
695 const int num_loc_dofs
696 )
697 {
698 // evaluate reference function hessian
699 typename AnaTraits_::HessianType hess = func_eval.hessian(trafo_data.img_point);
700
701 // subtract FE function hessian
702 for(int i(0); i < num_loc_dofs; ++i)
703 {
704 hess.add_vec_mat_outer_product(lvad[i], space_data.phi[i].hess, -DataType(1));
705 }
706
707 // compute norm of each component
708 ResultType result;
709 for(int k(0); k < AnaTraits_::image_dim; ++k)
710 result[k] = hess[k].norm_hessian_sqr();
711
712 return result;
713 }
714 };
715 } // namespace Intern
717
724 template<typename DataType_, int dim_>
725 struct VectorErrorInfo
726 {
728 static constexpr int dim = dim_;
729
731 bool have_h0;
733 bool have_h1;
735 bool have_h2;
737 bool have_l1;
739 bool have_lmax;
740
742 DataType_ norm_h0;
744 DataType_ norm_h1;
746 DataType_ norm_h2;
747
749 DataType_ norm_l1;
751 DataType_ norm_lmax;
752
758 Tiny::Vector<DataType_, dim_> norm_h0_comp;
759
765 Tiny::Vector<DataType_, dim_> norm_h1_comp;
766
772 Tiny::Vector<DataType_, dim_> norm_h2_comp;
773
779 Tiny::Vector<DataType_, dim_> norm_l1_comp;
780
786 Tiny::Vector<DataType_, dim_> norm_lmax_comp;
787
789 VectorErrorInfo() :
790 have_h0(false),
791 have_h1(false),
792 have_h2(false),
793 have_l1(false),
794 have_lmax(false),
795 norm_h0(DataType_(0)),
796 norm_h1(DataType_(0)),
797 norm_h2(DataType_(0)),
798 norm_l1(DataType_(0)),
799 norm_lmax(DataType_(0)),
800 norm_h0_comp(DataType_(0)),
801 norm_h1_comp(DataType_(0)),
802 norm_h2_comp(DataType_(0)),
803 norm_l1_comp(DataType_(0)),
804 norm_lmax_comp(DataType_(0))
805 {
806 }
807
809 template<typename DT2_>
810 VectorErrorInfo(const VectorErrorInfo<DT2_, dim_>& other) :
811 have_h0(other.have_h0),
812 have_h1(other.have_h1),
813 have_h2(other.have_h2),
814 have_l1(other.have_l1),
815 have_lmax(other.have_lmax),
816 norm_h0(DataType_(other.norm_h0)),
817 norm_h1(DataType_(other.norm_h1)),
818 norm_h2(DataType_(other.norm_h2)),
819 norm_l1(DataType_(other.norm_l1)),
820 norm_lmax(DataType_(other.norm_lmax)),
821 norm_h0_comp(other.norm_h0_comp),
822 norm_h1_comp(other.norm_h1_comp),
823 norm_h2_comp(other.norm_h2_comp),
824 norm_l1_comp(other.norm_l1_comp),
825 norm_lmax_comp(other.norm_lmax_comp)
826 {
827 }
828
830 template<typename DT2_>
831 VectorErrorInfo& operator=(const VectorErrorInfo<DT2_, dim_>& other)
832 {
833 have_h0 = other.have_h0;
834 have_h1 = other.have_h1;
835 have_h2 = other.have_h2;
836 have_l1 = other.have_l1;
837 have_lmax = other.have_lmax;
838 norm_h0 = DataType_(other.norm_h0);
839 norm_h1 = DataType_(other.norm_h1);
840 norm_h2 = DataType_(other.norm_h2);
841 norm_l1 = DataType_(other.norm_l1);
842 norm_lmax = DataType_(other.norm_lmax);
843 norm_h0_comp = other.norm_h0_comp;
844 norm_h1_comp = other.norm_h1_comp;
845 norm_h2_comp = other.norm_h2_comp;
846 norm_l1_comp = other.norm_l1_comp;
847 norm_lmax_comp = other.norm_lmax_comp;
848 }
849
859 void synchronize(const Dist::Comm& comm)
860 {
861 DataType_ verr[4+4*dim_] =
862 {
863 have_h0 ? Math::sqr(norm_h0) : DataType_(0),
864 have_h1 ? Math::sqr(norm_h1) : DataType_(0),
865 have_h2 ? Math::sqr(norm_h2) : DataType_(0),
866 have_l1 ? norm_l1 : DataType_(0)
867 };
868 for(int i(0); i < dim_; ++i)
869 {
870 verr[4+0*dim_+i] = (have_h0 ? Math::sqr(norm_h0_comp[i]) : DataType_(0));
871 verr[4+1*dim_+i] = (have_h0 ? Math::sqr(norm_h1_comp[i]) : DataType_(0));
872 verr[4+2*dim_+i] = (have_h0 ? Math::sqr(norm_h2_comp[i]) : DataType_(0));
873 verr[4+3*dim_+i] = (have_l1 ? norm_l1_comp[i] : DataType_(0));
874 }
875
876 comm.allreduce(verr, verr, std::size_t(4*dim+4), Dist::op_sum);
877
878 if(have_h0)
879 {
880 norm_h0 = Math::sqrt(verr[0]);
881 for(int i(0); i < dim_; ++i)
882 norm_h0_comp[i] = Math::sqrt(verr[4+i]);
883 }
884 if(have_h1)
885 {
886 norm_h1 = Math::sqrt(verr[1]);
887 for(int i(0); i < dim_; ++i)
888 norm_h1_comp[i] = Math::sqrt(verr[4+dim_+i]);
889 }
890 if(have_h2)
891 {
892 norm_h2 = Math::sqrt(verr[2]);
893 for(int i(0); i < dim_; ++i)
894 norm_h2_comp[i] = Math::sqrt(verr[4+2*dim_+i]);
895 }
896 if(have_l1)
897 {
898 norm_l1 = verr[3];
899 for(int i(0); i < dim_; ++i)
900 norm_l1_comp[i] = verr[4+3*dim_+i];
901 }
902 if(have_lmax)
903 {
904 DataType_ vmax[1+dim_] = { norm_lmax };
905 for(int i(0); i < dim_; ++i)
906 vmax[i+1] = norm_lmax_comp[i];
907 comm.allreduce(vmax, vmax, std::size_t(dim+1), Dist::op_max);
908 norm_lmax = vmax[0];
909 for(int i(0); i < dim_; ++i)
910 norm_lmax_comp[i] = vmax[i+1];
911 }
912 }
913
929 String format_string(int precision = 0, std::size_t pad_size = 8u, char pad_char = '.') const
930 {
931 String s;
932 if(have_h0)
933 {
934 s += String("H0-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h0, precision) + " [";
935 for(int i(0); i < dim_; ++i)
936 (s += " ") += stringify_fp_sci(norm_h0_comp[i], precision);
937 s += " ]\n";
938 }
939 if(have_h1)
940 {
941 s += String("H1-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h1, precision) + " [";
942 for(int i(0); i < dim_; ++i)
943 (s += " ") += stringify_fp_sci(norm_h1_comp[i], precision);
944 s += " ]\n";
945 }
946 if(have_h2)
947 {
948 s += String("H2-Error").pad_back(pad_size, pad_char) + ": " + stringify_fp_sci(norm_h2, precision) + " [";
949 for(int i(0); i < dim_; ++i)
950 (s += " ") += stringify_fp_sci(norm_h2_comp[i], precision);
951 s += " ]\n";
952 }
953 return s;
954 }
955
957 friend std::ostream& operator<<(std::ostream& os, const VectorErrorInfo& sei)
958 {
959 // print errors
960 if(sei.have_h0)
961 {
962 os << "H0-Error: " << stringify_fp_sci(sei.norm_h0) << " [";
963 for(int i(0); i < dim_; ++i)
964 os << " " << stringify_fp_sci(sei.norm_h0_comp[i]);
965 os << " ]\n";
966 }
967 if(sei.have_h1)
968 {
969 os << "H1-Error: " << stringify_fp_sci(sei.norm_h1) << " [";
970 for(int i(0); i < dim_; ++i)
971 os << " " << stringify_fp_sci(sei.norm_h1_comp[i]);
972 os << " ]\n";
973 }
974 if(sei.have_h2)
975 {
976 os << "H2-Error: " << stringify_fp_sci(sei.norm_h2) << " [";
977 for(int i(0); i < dim_; ++i)
978 os << " " << stringify_fp_sci(sei.norm_h2_comp[i]);
979 os << " ]\n";
980 }
981 return os;
982 }
983 }; // struct VectorErrorInfo
984
1005 template<int max_norm_ = 0>
1006 class VectorErrorComputer
1007 {
1008 static_assert(max_norm_ >= 0, "invalid max_norm_ parameter");
1009
1010 private:
1012
1015 static constexpr TrafoTags trafo_config = TrafoTags::img_point | TrafoTags::jac_det;
1016
1020 static constexpr SpaceTags space_config =
1021 ((max_norm_ >= 0) ? SpaceTags::value : SpaceTags::none) |
1022 ((max_norm_ >= 1) ? SpaceTags::grad : SpaceTags::none) |
1023 ((max_norm_ >= 2) ? SpaceTags::hess : SpaceTags::none);
1025
1026 public:
1050 template<
1051 typename Vector_,
1052 typename Function_,
1053 typename Space_,
1054 typename CubatureFactory_>
1055 static VectorErrorInfo<typename Vector_::DataType, Function_::ImageType::scalar_components> compute(
1056 const Vector_& vector,
1057 const Function_& function,
1058 const Space_& space,
1059 const CubatureFactory_& cubature_factory)
1060 {
1061 // make sure the function has everything we need
1062 static_assert(Function_::can_value, "function values are required for H0-Error");
1063 static_assert(Function_::can_grad || (max_norm_ < 1), "function gradients are required for H1-Error");
1064 static_assert(Function_::can_hess || (max_norm_ < 2), "function hessians are required for H2-Error");
1065
1066 // validate vector dimensions
1067 XASSERTM(vector.size() == space.get_num_dofs(), "invalid vector size");
1068
1069 // vector type
1070 typedef Vector_ VectorType;
1071 // analytic function type
1072 typedef Function_ FunctionType;
1073 // space type
1074 typedef Space_ SpaceType;
1076 static constexpr int dim = Function_::ImageType::scalar_components;
1078 typedef AsmTraits1
1079 <
1080 typename Vector_::DataType,
1081 SpaceType,
1082 trafo_config,
1083 space_config
1084 > AsmTraits;
1085 // data type
1086 typedef typename AsmTraits::DataType DataType;
1087
1088 // create the cubature rule
1089 typename AsmTraits::CubatureRuleType cubature_rule(Cubature::ctor_factory, cubature_factory);
1090
1091 // fetch the trafo
1092 const typename AsmTraits::TrafoType& trafo = space.get_trafo();
1093
1094 // define our analytic evaluation traits
1095 typedef Analytic::EvalTraits<DataType, Function_> AnalyticEvalTraits;
1096
1097 // create a function evaluator
1098 typename FunctionType::template Evaluator<AnalyticEvalTraits> func_eval(function);
1099
1100 // create a trafo evaluator
1101 typename AsmTraits::TrafoEvaluator trafo_eval(trafo);
1102
1103 // create a space evaluator and evaluation data
1104 typename AsmTraits::SpaceEvaluator space_eval(space);
1105
1106 // create a dof-mapping
1107 typename AsmTraits::DofMapping dof_mapping(space);
1108
1109 // create trafo evaluation data
1110 typename AsmTraits::TrafoEvalData trafo_data;
1111
1112 // create space evaluation data
1113 typename AsmTraits::SpaceEvalData space_data;
1114
1115 // create local vector data
1116 typedef Tiny::Vector<DataType, dim> VectorValue;
1117
1118 // create local vector data
1119 typename AsmTraits::template TLocalVector<VectorValue> lvad;
1120
1121 // create matrix scatter-axpy
1122 typename VectorType::GatherAxpy gather_axpy(vector);
1123
1124 // initialize result
1125 VectorErrorInfo<DataType, dim> info;
1126 info.have_h0 = (max_norm_ >= 0);
1127 info.have_h1 = (max_norm_ >= 1);
1128 info.have_h2 = (max_norm_ >= 2);
1129 info.have_l1 = (max_norm_ >= 0);
1130 info.have_lmax = (max_norm_ >= 0);
1131
1132 // H0/H1/H2 helpers
1133 typedef Intern::VecHelperH0<(max_norm_ >= 0), AsmTraits, AnalyticEvalTraits> VecH0;
1134 typedef Intern::VecHelperH1<(max_norm_ >= 1), AsmTraits, AnalyticEvalTraits> VecH1;
1135 typedef Intern::VecHelperH2<(max_norm_ >= 2), AsmTraits, AnalyticEvalTraits> VecH2;
1136
1137 // loop over all cells of the mesh
1138 for(typename AsmTraits::CellIterator cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
1139 {
1140 // format local vector
1141 lvad.format();
1142
1143 // initialize dof-mapping
1144 dof_mapping.prepare(cell);
1145
1146 // incorporate local matrix
1147 gather_axpy(lvad, dof_mapping);
1148
1149 // finish dof-mapping
1150 dof_mapping.finish();
1151
1152 // prepare trafo evaluator
1153 trafo_eval.prepare(cell);
1154
1155 // prepare space evaluator
1156 space_eval.prepare(trafo_eval);
1157
1158 // fetch number of local dofs
1159 int num_loc_dofs = space_eval.get_num_local_dofs();
1160
1161 // loop over all quadrature points and integrate
1162 for(int k(0); k < cubature_rule.get_num_points(); ++k)
1163 {
1164 // compute trafo data
1165 trafo_eval(trafo_data, cubature_rule.get_point(k));
1166
1167 // compute basis function data
1168 space_eval(space_data, trafo_data);
1169
1170 // compute integration weight
1171 const DataType omega = trafo_data.jac_det * cubature_rule.get_weight(k);
1172
1173 // get absolute point error
1174 auto abs_err = VecH0::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs);
1175
1176 // update H0/L1/Lmax results
1177 for(int i(0); i < dim; ++i)
1178 {
1179 info.norm_lmax_comp[i] = Math::max(info.norm_lmax_comp[i], abs_err[i]);
1180 info.norm_l1_comp[i] += omega * abs_err[i];
1181 info.norm_h0_comp[i] += omega * Math::sqr(abs_err[i]);
1182 }
1183
1184 // update H1/H2 results
1185 info.norm_h1_comp.axpy(omega, VecH1::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs));
1186 info.norm_h2_comp.axpy(omega, VecH2::eval(func_eval, trafo_data, space_data, lvad, num_loc_dofs));
1187
1188 // continue with next cubature point
1189 }
1190
1191 // finish evaluators
1192 space_eval.finish();
1193 trafo_eval.finish();
1194
1195 // continue with next cell
1196 }
1197
1198 // compute the rest
1199 for(int i(0); i < dim; ++i)
1200 {
1201 info.norm_h0 += info.norm_h0_comp[i];
1202 info.norm_h1 += info.norm_h1_comp[i];
1203 info.norm_h2 += info.norm_h2_comp[i];
1204 info.norm_l1 += info.norm_l1_comp[i];
1205 info.norm_lmax = Math::max(info.norm_lmax, info.norm_lmax_comp[i]);
1206 info.norm_h0_comp[i] = Math::sqrt(info.norm_h0_comp[i]);
1207 info.norm_h1_comp[i] = Math::sqrt(info.norm_h1_comp[i]);
1208 info.norm_h2_comp[i] = Math::sqrt(info.norm_h2_comp[i]);
1209 }
1210
1211 // take square-roots
1212 info.norm_h0 = Math::sqrt(info.norm_h0);
1213 info.norm_h1 = Math::sqrt(info.norm_h1);
1214 info.norm_h2 = Math::sqrt(info.norm_h2);
1215
1216 // okay, that's it
1217 return info;
1218 }
1219 }; // class VectorErrorComputer
1220
1221 } // namespace Assembly
1222} // namespace FEAT
XASSERTM
#define XASSERTM(expr, msg)
Assertion macro definition with custom message.
Definition: assertion.hpp:263
FEAT::Assembly::AsmTraits1
Common single-space assembly traits class template.
Definition: asm_traits.hpp:83
FEAT::Assembly::ScalarErrorComputer
Scalar H0/H1/H2-error computer class template.
Definition: error_computer.hpp:376
FEAT::Assembly::ScalarErrorComputer::compute
static ScalarErrorInfo< typename Vector_::DataType > compute(const Vector_ &vector, const Function_ &function, const Space_ &space, const CubatureFactory_ &cubature_factory)
Computes the H0/H1/H2-error.
Definition: error_computer.hpp:424
FEAT::Assembly::VectorErrorComputer
Vector H0/H1/H2-error computer class template.
Definition: error_computer.hpp:1007
FEAT::Assembly::VectorErrorComputer::compute
static VectorErrorInfo< typename Vector_::DataType, Function_::ImageType::scalar_components > compute(const Vector_ &vector, const Function_ &function, const Space_ &space, const CubatureFactory_ &cubature_factory)
Computes the H0/H1/H2-error.
Definition: error_computer.hpp:1055
FEAT::Dist::Comm
Communicator class.
Definition: dist.hpp:1349
FEAT::Dist::Comm::allreduce
void allreduce(const void *sendbuf, void *recvbuf, std::size_t count, const Datatype &datatype, const Operation &op) const
Blocking All-Reduce.
Definition: dist.cpp:655
FEAT::String
String class implementation.
Definition: string.hpp:46
FEAT::String::pad_back
String pad_back(size_type len, char c=' ') const
Pads the back of the string up to a desired length.
Definition: string.hpp:415
FEAT::Tiny::Vector
Tiny Vector class template.
Definition: tiny_algebra.hpp:52
FEAT::Tiny::Vector::axpy
CUDA_HOST_DEVICE Vector & axpy(DataType alpha, const Vector< T_, n_, snx_ > &x)
Adds another scaled vector onto this vector.
Definition: tiny_algebra.hpp:609
FEAT::Dist::op_max
const Operation op_max(MPI_MAX)
Operation wrapper for MPI_MAX.
Definition: dist.hpp:273
FEAT::Dist::op_sum
const Operation op_sum(MPI_SUM)
Operation wrapper for MPI_SUM.
Definition: dist.hpp:271
FEAT::Math::sqrt
T_ sqrt(T_ x)
Returns the square-root of a value.
Definition: math.hpp:300
FEAT::Math::sqr
T_ sqr(T_ x)
Returns the square of a value.
Definition: math.hpp:95
FEAT::Math::max
T_ max(T_ a, T_ b)
Returns the maximum of two values.
Definition: math.hpp:137
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::stringify_fp_sci
String stringify_fp_sci(DataType_ value, int precision=0, int width=0, bool sign=false)
Prints a floating point value to a string in scientific notation.
Definition: string.hpp:1088
FEAT::SpaceTags
SpaceTags
Space configuration tags enum.
Definition: eval_tags.hpp:97
FEAT::SpaceTags::value
@ value
specifies whether the space should supply basis function values
FEAT::SpaceTags::hess
@ hess
specifies whether the space should supply basis function hessians
FEAT::SpaceTags::ref_grad
@ ref_grad
specifies whether the space should supply reference basis function gradients
FEAT::SpaceTags::grad
@ grad
specifies whether the space should supply basis function gradients
FEAT::TrafoTags
TrafoTags
Trafo configuration tags enum.
Definition: eval_tags.hpp:22
FEAT::TrafoTags::img_point
@ img_point
specifies whether the trafo should supply image point coordinates
FEAT::TrafoTags::jac_mat
@ jac_mat
specifies whether the trafo should supply jacobian matrices
FEAT::TrafoTags::jac_det
@ jac_det
specifies whether the trafo should supply jacobian determinants
FEAT::Assembly::ScalarErrorInfo
Scalar Error information structure.
Definition: error_computer.hpp:203
FEAT::Assembly::ScalarErrorInfo::have_h1
bool have_h1
Specifies whether the H1-error was computed.
Definition: error_computer.hpp:207
FEAT::Assembly::ScalarErrorInfo::norm_h1
DataType_ norm_h1
The computed H1-error.
Definition: error_computer.hpp:218
FEAT::Assembly::ScalarErrorInfo::operator=
ScalarErrorInfo & operator=(const ScalarErrorInfo< DT2_ > &other)
conversion assignment operator
Definition: error_computer.hpp:260
FEAT::Assembly::ScalarErrorInfo::have_lmax
bool have_lmax
Specifies whether the Lmax-error was computed.
Definition: error_computer.hpp:213
FEAT::Assembly::ScalarErrorInfo::have_h2
bool have_h2
Specifies whether the H2-error was computed.
Definition: error_computer.hpp:209
FEAT::Assembly::ScalarErrorInfo::format_string
String format_string(int precision=0, std::size_t pad_size=8u, char pad_char='.') const
Formats the error information as a string.
Definition: error_computer.hpp:321
FEAT::Assembly::ScalarErrorInfo::norm_lmax
DataType_ norm_lmax
The computed Lmax-error.
Definition: error_computer.hpp:225
FEAT::Assembly::ScalarErrorInfo::synchronize
void synchronize(const Dist::Comm &comm)
Synchronizes the error information over a communicator.
Definition: error_computer.hpp:283
FEAT::Assembly::ScalarErrorInfo::ScalarErrorInfo
ScalarErrorInfo(const ScalarErrorInfo< DT2_ > &other)
conversion constructor
Definition: error_computer.hpp:244
FEAT::Assembly::ScalarErrorInfo::have_h0
bool have_h0
Specifies whether the H0-error was computed.
Definition: error_computer.hpp:205
FEAT::Assembly::ScalarErrorInfo::operator<<
friend std::ostream & operator<<(std::ostream &os, const ScalarErrorInfo &sei)
prints the info to an output stream
Definition: error_computer.hpp:334
FEAT::Assembly::ScalarErrorInfo::ScalarErrorInfo
ScalarErrorInfo()
standard constructor
Definition: error_computer.hpp:228
FEAT::Assembly::ScalarErrorInfo::have_l1
bool have_l1
Specifies whether the L1-error was computed.
Definition: error_computer.hpp:211
FEAT::Assembly::ScalarErrorInfo::norm_h2
DataType_ norm_h2
The computed H2-error.
Definition: error_computer.hpp:220
FEAT::Assembly::ScalarErrorInfo::norm_h0
DataType_ norm_h0
The computed H0-error.
Definition: error_computer.hpp:216
FEAT::Assembly::ScalarErrorInfo::norm_l1
DataType_ norm_l1
The computed L1-error.
Definition: error_computer.hpp:223
FEAT::Assembly::VectorErrorInfo
Vector Error information structure.
Definition: error_computer.hpp:726
FEAT::Assembly::VectorErrorInfo::have_lmax
bool have_lmax
Specifies whether the Lmax-error was computed.
Definition: error_computer.hpp:739
FEAT::Assembly::VectorErrorInfo::have_h1
bool have_h1
Specifies whether the H1-error was computed.
Definition: error_computer.hpp:733
FEAT::Assembly::VectorErrorInfo::norm_h1
DataType_ norm_h1
The computed H1-error.
Definition: error_computer.hpp:744
FEAT::Assembly::VectorErrorInfo::norm_h0
DataType_ norm_h0
The computed H0-error.
Definition: error_computer.hpp:742
FEAT::Assembly::VectorErrorInfo::norm_h0_comp
Tiny::Vector< DataType_, dim_ > norm_h0_comp
Vector field components H0-errors.
Definition: error_computer.hpp:758
FEAT::Assembly::VectorErrorInfo::operator=
VectorErrorInfo & operator=(const VectorErrorInfo< DT2_, dim_ > &other)
conversion assignment operator
Definition: error_computer.hpp:831
FEAT::Assembly::VectorErrorInfo::dim
static constexpr int dim
The dimension of the analysed vector field.
Definition: error_computer.hpp:728
FEAT::Assembly::VectorErrorInfo::norm_l1_comp
Tiny::Vector< DataType_, dim_ > norm_l1_comp
Vector field components L1-errors.
Definition: error_computer.hpp:779
FEAT::Assembly::VectorErrorInfo::VectorErrorInfo
VectorErrorInfo()
standard constructor
Definition: error_computer.hpp:789
FEAT::Assembly::VectorErrorInfo::synchronize
void synchronize(const Dist::Comm &comm)
Synchronizes the error information over a communicator.
Definition: error_computer.hpp:859
FEAT::Assembly::VectorErrorInfo::norm_h1_comp
Tiny::Vector< DataType_, dim_ > norm_h1_comp
Vector field components H1-errors.
Definition: error_computer.hpp:765
FEAT::Assembly::VectorErrorInfo::norm_l1
DataType_ norm_l1
The computed L1-error.
Definition: error_computer.hpp:749
FEAT::Assembly::VectorErrorInfo::VectorErrorInfo
VectorErrorInfo(const VectorErrorInfo< DT2_, dim_ > &other)
conversion constructor
Definition: error_computer.hpp:810
FEAT::Assembly::VectorErrorInfo::norm_lmax_comp
Tiny::Vector< DataType_, dim_ > norm_lmax_comp
Vector field components Lmax-errors.
Definition: error_computer.hpp:786
FEAT::Assembly::VectorErrorInfo::have_h0
bool have_h0
Specifies whether the H0-error was computed.
Definition: error_computer.hpp:731
FEAT::Assembly::VectorErrorInfo::format_string
String format_string(int precision=0, std::size_t pad_size=8u, char pad_char='.') const
Formats the error information as a string.
Definition: error_computer.hpp:929
FEAT::Assembly::VectorErrorInfo::operator<<
friend std::ostream & operator<<(std::ostream &os, const VectorErrorInfo &sei)
prints the info to an output stream
Definition: error_computer.hpp:957
FEAT::Assembly::VectorErrorInfo::norm_lmax
DataType_ norm_lmax
The computed Lmax-error.
Definition: error_computer.hpp:751
FEAT::Assembly::VectorErrorInfo::norm_h2_comp
Tiny::Vector< DataType_, dim_ > norm_h2_comp
Vector field components H2-errors.
Definition: error_computer.hpp:772
FEAT::Assembly::VectorErrorInfo::norm_h2
DataType_ norm_h2
The computed H2-error.
Definition: error_computer.hpp:746
FEAT::Assembly::VectorErrorInfo::have_l1
bool have_l1
Specifies whether the L1-error was computed.
Definition: error_computer.hpp:737
FEAT::Assembly::VectorErrorInfo::have_h2
bool have_h2
Specifies whether the H2-error was computed.
Definition: error_computer.hpp:735