AFEPack
Functional.templates.h
浏览该文件的文档。
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