@@ -114,6 +114,40 @@ class SelfProduct : public SInt {
114114 OmegahMeshField &omf;
115115};
116116
117+ template <typename EstimationT, typename EpsStarT, typename ErrorNormT,
118+ typename WeightsT, typename DerJacobiansT>
119+ MeshField::Real
120+ error_par_reduce_impl (EstimationT &estimation, EpsStarT eps_star_atPts,
121+ ErrorNormT errorNorm, WeightsT w, DerJacobiansT dV,
122+ const size_t numPtsPerElem) {
123+ const auto meshDim = static_cast <double >(estimation.mesh .dim ());
124+ const auto orderP = static_cast <double >(estimation.recovered_order );
125+ const auto &eps = estimation.eps ;
126+ MeshField::Real r = 0 ;
127+ Kokkos::parallel_reduce (
128+ " eval" , estimation.mesh .nelems (),
129+ KOKKOS_LAMBDA (const int &elm, MeshField::Real &r_local) {
130+ const auto first = elm * numPtsPerElem;
131+ const auto last = first + numPtsPerElem;
132+ const auto eps_elm = eps[elm];
133+ MeshField::Real sum = 0 ;
134+ for (auto pt = first; pt < last; pt++) {
135+ const auto eps_star_Pt = eps_star_atPts (pt, 0 );
136+ const auto diff = eps_elm - eps_star_Pt;
137+ const auto wPt = w (pt);
138+ const auto dVPt = dV (pt);
139+ sum += (diff * diff) * wPt * dVPt;
140+ }
141+ assert (sum > 0 );
142+ r_local += Kokkos::pow (Kokkos::sqrt (sum),
143+ ((2 * meshDim) / (2 * orderP + meshDim)));
144+ errorNorm (elm) =
145+ Kokkos::pow (Kokkos::sqrt (sum), -(2 / (2 * orderP + meshDim)));
146+ },
147+ r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$
148+ return r;
149+ }
150+
117151/* computes the integral over the element of the
118152 sum of the squared differences between the
119153 original and recovered fields */
@@ -138,28 +172,8 @@ class Error : public SInt {
138172 double orderP = estimation.recovered_order ;
139173
140174 const auto &eps = estimation.eps ;
141- r = 0 ;
142- Kokkos::parallel_reduce (
143- " eval" , estimation.mesh .nelems (),
144- KOKKOS_CLASS_LAMBDA (const int &elm, MeshField::Real &r_local) {
145- const auto first = elm * numPtsPerElem;
146- const auto last = first + numPtsPerElem;
147- const auto eps_elm = eps[elm];
148- MeshField::Real sum = 0 ;
149- for (auto pt = first; pt < last; pt++) {
150- const auto eps_star_Pt = eps_star_atPts (pt, 0 );
151- const auto diff = eps_elm - eps_star_Pt;
152- const auto wPt = w (pt);
153- const auto dVPt = dV (pt);
154- sum += (diff * diff) * wPt * dVPt;
155- }
156- assert (sum > 0 );
157- r_local += Kokkos::pow (Kokkos::sqrt (sum),
158- ((2 * meshDim) / (2 * orderP + meshDim)));
159- errorNorm (elm) =
160- Kokkos::pow (Kokkos::sqrt (sum), -(2 / (2 * orderP + meshDim)));
161- },
162- r); // $\sum_{i=1}^n \|e_\epsilon\|^{\frac{2d}{2p+d}}$
175+ r = error_par_reduce_impl (estimation, eps_star_atPts, errorNorm, w, dV,
176+ numPtsPerElem);
163177 }
164178 EstimationT &estimation;
165179 OmegahMeshField &omf;
0 commit comments