AFEPack
|
00001 00011 #ifndef _Functional_templates_h_ 00012 #define _Functional_templates_h_ 00013 00014 #include "Functional.h" 00015 00016 template <class value_type, int DIM> 00017 value_type Functional::L1Norm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00018 { 00019 value_type norm = 0; 00020 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00021 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00022 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00023 for (;the_element != end_element;the_element ++) { 00024 double volume = the_element->templateElement().volume(); 00025 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00026 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00027 int n_quadrature_point = quad_info.n_quadraturePoint(); 00028 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00029 std::vector<double> f_value = f.value(q_point, *the_element); 00030 for (int l = 0;l < n_quadrature_point;l ++) { 00031 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00032 norm += Jxw*fabs(f_value[l]); 00033 } 00034 } 00035 return norm; 00036 }; 00037 00038 00039 00040 template <class value_type, int DIM> 00041 value_type Functional::L2Norm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00042 { 00043 value_type norm = 0; 00044 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00045 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00046 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00047 for (;the_element != end_element;the_element ++) { 00048 double volume = the_element->templateElement().volume(); 00049 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00050 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00051 int n_quadrature_point = quad_info.n_quadraturePoint(); 00052 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00053 std::vector<double> f_value = f.value(q_point, *the_element); 00054 for (int l = 0;l < n_quadrature_point;l ++) { 00055 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00056 norm += Jxw*f_value[l]*f_value[l]; 00057 } 00058 } 00059 return sqrt(fabs(norm)); 00060 }; 00061 00062 00063 00064 template <class value_type, int DIM> 00065 value_type Functional::L0Norm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00066 { 00067 value_type norm = 0; 00068 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00069 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00070 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00071 for (;the_element != end_element;the_element ++) { 00072 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00073 int n_quadrature_point = quad_info.n_quadraturePoint(); 00074 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00075 std::vector<double> f_value = f.value(q_point, *the_element); 00076 for (int l = 0;l < n_quadrature_point;l ++) { 00077 f_value[l] = fabs(f_value[l]); 00078 if (f_value[l] > norm) norm = f_value[l]; 00079 } 00080 } 00081 return norm; 00082 }; 00083 00084 00085 00086 template <class value_type, int DIM> 00087 value_type Functional::LpNorm(FEMFunction<value_type, DIM>& f, double p, int algebric_accuracy) 00088 { 00089 value_type norm = 0; 00090 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00091 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00092 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00093 for (;the_element != end_element;the_element ++) { 00094 double volume = the_element->templateElement().volume(); 00095 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00096 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00097 int n_quadrature_point = quad_info.n_quadraturePoint(); 00098 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00099 std::vector<double> f_value = f.value(q_point, *the_element); 00100 for (int l = 0;l < n_quadrature_point;l ++) { 00101 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00102 norm += Jxw*pow(f_value[l], p); 00103 } 00104 } 00105 return pow(norm, 1./p); 00106 }; 00107 00108 00109 00110 template <class value_type, int DIM> 00111 value_type Functional::W11Seminorm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00112 { 00113 int i; 00114 value_type norm, n[DIM]; 00115 for (i = 0;i < DIM;i ++) n[i] = 0; 00116 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00117 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00118 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00119 for (;the_element != end_element;the_element ++) { 00120 double volume = the_element->templateElement().volume(); 00121 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00122 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00123 int n_quadrature_point = quad_info.n_quadraturePoint(); 00124 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00125 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00126 for (int l = 0;l < n_quadrature_point;l ++) { 00127 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00128 for (i = 0;i < DIM;i ++) 00129 n[i] += Jxw*fabs(f_gradient[l][i]); 00130 } 00131 } 00132 for (norm = 0, i = 0;i < DIM;i ++) norm += n[i]; 00133 return norm; 00134 }; 00135 00136 00137 00138 template <class value_type, int DIM> 00139 value_type Functional::H1Seminorm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00140 { 00141 value_type norm = 0; 00142 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00143 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00144 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00145 for (;the_element != end_element;the_element ++) { 00146 double volume = the_element->templateElement().volume(); 00147 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00148 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00149 int n_quadrature_point = quad_info.n_quadraturePoint(); 00150 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00151 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00152 for (int l = 0;l < n_quadrature_point;l ++) { 00153 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00154 for (int i = 0;i < DIM;i ++) 00155 norm += Jxw*f_gradient[l][i]*f_gradient[l][i]; 00156 } 00157 } 00158 return sqrt(fabs(norm)); 00159 }; 00160 00161 00162 00163 template <class value_type, int DIM> 00164 value_type Functional::W1pSeminorm(FEMFunction<value_type, DIM>& f, double p, int algebric_accuracy) 00165 { 00166 value_type norm = 0; 00167 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00168 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00169 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00170 for (;the_element != end_element;the_element ++) { 00171 double volume = the_element->templateElement().volume(); 00172 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00173 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00174 int n_quadrature_point = quad_info.n_quadraturePoint(); 00175 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00176 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00177 for (int l = 0;l < n_quadrature_point;l ++) { 00178 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00179 for (int i = 0;i < DIM;i ++) 00180 norm += Jxw*pow(f_gradient[l][i], p); 00181 } 00182 } 00183 return pow(norm, 1./p); 00184 }; 00185 00186 00187 00188 template <class value_type, int DIM> 00189 value_type Functional::W10Seminorm(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00190 { 00191 int i; 00192 value_type norm, n[DIM]; 00193 for (i = 0;i < DIM;i ++) n[i] = 0; 00194 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00195 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00196 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00197 for (;the_element != end_element;the_element ++) { 00198 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00199 int n_quadrature_point = quad_info.n_quadraturePoint(); 00200 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00201 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00202 for (int l = 0;l < n_quadrature_point;l ++) { 00203 for (i = 0;i < DIM;i ++) { 00204 f_gradient[l][i] = fabs(f_gradient[l][i]); 00205 if (n[i] < f_gradient[l][i]) 00206 n[i] = f_gradient[l][i]; 00207 } 00208 } 00209 } 00210 for (norm = 0, i = 0;i < DIM;i ++) norm += n[i]; 00211 return norm; 00212 }; 00213 00214 00215 00216 template <class value_type, int DIM> 00217 value_type Functional::meanValue(FEMFunction<value_type, DIM>& f, int algebric_accuracy) 00218 { 00219 value_type val = 0; 00220 double vol = 0; 00221 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00222 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00223 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00224 for (;the_element != end_element;the_element ++) { 00225 double volume = the_element->templateElement().volume(); 00226 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00227 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00228 int n_quadrature_point = quad_info.n_quadraturePoint(); 00229 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00230 std::vector<double> f_value = f.value(q_point, *the_element); 00231 for (int l = 0;l < n_quadrature_point;l ++) { 00232 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00233 val += Jxw*f_value[l]; 00234 vol += Jxw; 00235 } 00236 } 00237 return val/vol; 00238 }; 00239 00240 00241 00242 template <class value_type, int DIM> 00243 value_type Functional::meanValue(const Function<value_type>& f, 00244 FEMSpace<value_type,DIM>& fem_space, 00245 int algebric_accuracy) 00246 { 00247 value_type val = 0; 00248 double vol = 0; 00249 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00250 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00251 for (;the_element != end_element;the_element ++) { 00252 double volume = the_element->templateElement().volume(); 00253 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00254 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00255 int n_quadrature_point = quad_info.n_quadraturePoint(); 00256 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00257 for (int l = 0;l < n_quadrature_point;l ++) { 00258 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00259 val += Jxw*f.value(q_point[l]); 00260 vol += Jxw; 00261 } 00262 } 00263 return val/vol; 00264 }; 00265 00266 00267 template <class value_type, int DIM> 00268 value_type Functional::L1Error(FEMFunction<value_type, DIM>& f, 00269 const Function<value_type>& f1, 00270 int algebric_accuracy) 00271 { 00272 value_type error = 0; 00273 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00274 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00275 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00276 for (;the_element != end_element;the_element ++) { 00277 double volume = the_element->templateElement().volume(); 00278 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00279 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00280 int n_quadrature_point = quad_info.n_quadraturePoint(); 00281 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00282 std::vector<double> f_value = f.value(q_point, *the_element); 00283 for (int l = 0;l < n_quadrature_point;l ++) { 00284 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00285 double df_value = f1.value(q_point[l]) - f_value[l]; 00286 error += Jxw*fabs(df_value); 00287 } 00288 } 00289 return error; 00290 }; 00291 00292 00293 00294 template <class value_type, int DIM> 00295 value_type Functional::L2Error(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, int algebric_accuracy) 00296 { 00297 value_type error = 0; 00298 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00299 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00300 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00301 for (;the_element != end_element;the_element ++) { 00302 double volume = the_element->templateElement().volume(); 00303 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00304 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00305 int n_quadrature_point = quad_info.n_quadraturePoint(); 00306 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00307 std::vector<double> f_value = f.value(q_point, *the_element); 00308 for (int l = 0;l < n_quadrature_point;l ++) { 00309 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00310 double df_value = f1.value(q_point[l]) - f_value[l]; 00311 error += Jxw*df_value*df_value; 00312 } 00313 } 00314 return sqrt(fabs(error)); 00315 }; 00316 00317 00318 template <class value_type, int DIM> 00319 value_type Functional::L0Error(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, int algebric_accuracy) 00320 { 00321 value_type error = 0; 00322 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00323 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00324 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00325 for (;the_element != end_element;the_element ++) { 00326 double volume = the_element->templateElement().volume(); 00327 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00328 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00329 int n_quadrature_point = quad_info.n_quadraturePoint(); 00330 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00331 std::vector<double> f_value = f.value(q_point, *the_element); 00332 for (int l = 0;l < n_quadrature_point;l ++) { 00333 double df_value = f1.value(q_point[l]) - f_value[l]; 00334 df_value = fabs(df_value); 00335 if (df_value > error) error = df_value; 00336 } 00337 } 00338 return error; 00339 }; 00340 00341 00342 template <class value_type, int DIM> 00343 value_type Functional::LpError(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, double p, int algebric_accuracy) 00344 { 00345 value_type error = 0; 00346 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00347 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00348 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00349 for (;the_element != end_element;the_element ++) { 00350 double volume = the_element->templateElement().volume(); 00351 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00352 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00353 int n_quadrature_point = quad_info.n_quadraturePoint(); 00354 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00355 std::vector<double> f_value = f.value(q_point, *the_element); 00356 for (int l = 0;l < n_quadrature_point;l ++) { 00357 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00358 double df_value = f1.value(q_point[l]) - f_value[l]; 00359 error += Jxw*pow(fabs(df_value), p); 00360 } 00361 } 00362 return pow(error, 1.0/p); 00363 }; 00364 00365 00366 00367 template <class value_type, int DIM> 00368 value_type Functional::W11SemiError(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, int algebric_accuracy) 00369 { 00370 int i; 00371 value_type error, n[DIM]; 00372 for (i = 0;i < DIM;i ++) n[i] = 0; 00373 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00374 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00375 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00376 for (;the_element != end_element;the_element ++) { 00377 double volume = the_element->templateElement().volume(); 00378 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00379 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00380 int n_quadrature_point = quad_info.n_quadraturePoint(); 00381 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00382 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00383 for (int l = 0;l < n_quadrature_point;l ++) { 00384 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00385 std::vector<value_type> f1_gradient = f1.gradient(q_point[l]); 00386 for (i = 0;i < DIM;i ++) { 00387 n[i] += Jxw*fabs(f_gradient[l][i] - f1_gradient[i]); 00388 } 00389 } 00390 } 00391 for (error = 0, i = 0;i < DIM;i ++) error += n[i]; 00392 return error; 00393 }; 00394 00395 00396 00397 template <class value_type, int DIM> 00398 value_type Functional::H1SemiError(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, int algebric_accuracy) 00399 { 00400 int i; 00401 value_type error, n[DIM]; 00402 for (i = 0;i < DIM;i ++) n[i] = 0; 00403 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00404 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00405 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00406 for (;the_element != end_element;the_element ++) { 00407 double volume = the_element->templateElement().volume(); 00408 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00409 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00410 int n_quadrature_point = quad_info.n_quadraturePoint(); 00411 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00412 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00413 for (int l = 0;l < n_quadrature_point;l ++) { 00414 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00415 std::vector<value_type> f1_gradient = f1.gradient(q_point[l]); 00416 for (i = 0;i < DIM;i ++) { 00417 double df_value = f_gradient[l][i] - f1_gradient[i]; 00418 n[i] += Jxw*df_value*df_value; 00419 } 00420 } 00421 } 00422 for (error = 0, i = 0;i < DIM;i ++) error += n[i]; 00423 return sqrt(fabs(error)); 00424 }; 00425 00426 00427 00428 template <class value_type, int DIM> 00429 value_type Functional::W1pSemiError(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, double p, int algebric_accuracy) 00430 { 00431 int i; 00432 value_type error, n[DIM]; 00433 for (i = 0;i < DIM;i ++) n[i] = 0; 00434 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00435 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00436 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00437 for (;the_element != end_element;the_element ++) { 00438 double volume = the_element->templateElement().volume(); 00439 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00440 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00441 int n_quadrature_point = quad_info.n_quadraturePoint(); 00442 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00443 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00444 for (int l = 0;l < n_quadrature_point;l ++) { 00445 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00446 std::vector<value_type> f1_gradient = f1.gradient(q_point[l]); 00447 for (i = 0;i < DIM;i ++) { 00448 double df_value = fabs(f_gradient[l][i] - f1_gradient[i]); 00449 n[i] += Jxw*pow(df_value, p); 00450 } 00451 } 00452 } 00453 for (error = 0, i = 0;i < DIM;i ++) error += n[i]; 00454 return pow(error, 1./p); 00455 }; 00456 00457 00458 00459 template <class value_type, int DIM> 00460 value_type Functional::W10SemiError(FEMFunction<value_type, DIM>& f, const Function<value_type>& f1, int algebric_accuracy) 00461 { 00462 int i; 00463 value_type error, n[DIM]; 00464 for (i = 0;i < DIM;i ++) n[i] = 0; 00465 FEMSpace<value_type,DIM>& fem_space = f.femSpace(); 00466 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00467 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00468 for (;the_element != end_element;the_element ++) { 00469 double volume = the_element->templateElement().volume(); 00470 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00471 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00472 int n_quadrature_point = quad_info.n_quadraturePoint(); 00473 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00474 std::vector<std::vector<double> > f_gradient = f.gradient(q_point, *the_element); 00475 for (int l = 0;l < n_quadrature_point;l ++) { 00476 std::vector<value_type> f1_gradient = f1.gradient(q_point[l]); 00477 for (i = 0;i < DIM;i ++) { 00478 double df_value = fabs(f_gradient[l][i] - f1_gradient[i]); 00479 if (n[i] < df_value) n[i] = df_value; 00480 } 00481 } 00482 } 00483 for (error = 0, i = 0;i < DIM;i ++) error += n[i]; 00484 return error; 00485 }; 00486 00487 00488 template <class value_type, int DIM> 00489 value_type Functional::L1Norm(const Function<value_type>& f, 00490 FEMSpace<value_type,DIM>& fem_space, 00491 int algebric_accuracy) 00492 { 00493 value_type norm = 0; 00494 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00495 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00496 for (;the_element != end_element;the_element ++) { 00497 double volume = the_element->templateElement().volume(); 00498 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00499 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00500 int n_quadrature_point = quad_info.n_quadraturePoint(); 00501 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00502 for (int l = 0;l < n_quadrature_point;l ++) { 00503 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00504 double f_value = f.value(q_point[l]); 00505 norm += Jxw*fabs(f_value); 00506 } 00507 } 00508 return norm; 00509 }; 00510 00511 00512 00513 template <class value_type, int DIM> 00514 value_type Functional::L2Norm(const Function<value_type>& f, 00515 FEMSpace<value_type,DIM>& fem_space, 00516 int algebric_accuracy) 00517 { 00518 value_type norm = 0; 00519 typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement(); 00520 typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement(); 00521 for (;the_element != end_element;the_element ++) { 00522 double volume = the_element->templateElement().volume(); 00523 const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy); 00524 std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint()); 00525 int n_quadrature_point = quad_info.n_quadraturePoint(); 00526 std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint()); 00527 for (int l = 0;l < n_quadrature_point;l ++) { 00528 double Jxw = quad_info.weight(l)*jacobian[l]*volume; 00529 double f_value = f.value(q_point[l]); 00530 norm += Jxw*f_value*f_value; 00531 } 00532 } 00533 return sqrt(fabs(norm)); 00534 }; 00535 00536 00537 00538 #endif // _Functional_templates_h_ 00539