AFEPack
Operator.project.templates.h
浏览该文件的文档。
00001 
00011 #ifndef _Operator_project_templates_h_
00012 #define _Operator_project_templates_h_
00013 
00014 #include "Operator.h"
00015 
00016 template <class value_type, int DIM>
00017 void Operator::L2Project(const FEMFunction<value_type, DIM>& f0, 
00018                          FEMFunction<value_type, DIM>& f1, 
00019                          Method method, 
00020                          int algebric_accuracy)
00021 {
00022   if (&(f0.femSpace().mesh()) == &(f1.femSpace().mesh())) {
00023     switch (method) {
00024     case MASS_ACCUMULATION: {
00025       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00026       unsigned int n_dof = fem_space.n_dof();
00027       Vector<double> mass_accumulation(n_dof);
00028       f1.Vector<value_type>::operator=(0.0);
00029       typename FEMSpace<value_type,DIM>::ConstElementIterator 
00030         the_element0 = f0.femSpace().beginElement(),
00031         the_element = fem_space.beginElement();
00032       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00033       for (;the_element != end_element;the_element ++, the_element0 ++) {
00034         double volume = the_element->templateElement().volume();
00035         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00036         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00037         const std::vector<int>& element_dof = the_element->dof();
00038         unsigned int n_element_dof = element_dof.size();
00039         int n_quadrature_point = quad_info.n_quadraturePoint();
00040         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00041         std::vector<value_type> f0_value = f0.value(q_point, *the_element0);
00042         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00043         for (int l = 0;l < n_quadrature_point;l ++) {
00044           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00045           for (unsigned int j = 0;j < n_element_dof;j ++) {
00046             f1(element_dof[j]) += Jxw*f0_value[l]*basis_value[j][l];
00047             mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
00048           }
00049         }
00050       }
00051       for (unsigned int i = 0;i < n_dof;i ++) 
00052         f1(i) /= mass_accumulation(i);
00053       break;
00054     }
00055     case LEAST_SQUARE: {
00056       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00057       f1.Vector<value_type>::operator=(0.0);
00058       unsigned int n_dof = fem_space.n_dof();
00059       MassMatrix<DIM, value_type> mass_matrix(fem_space);
00060       mass_matrix.algebricAccuracy() = algebric_accuracy;
00061       mass_matrix.build();
00062       Vector<double> rhs(n_dof);
00063       typename FEMSpace<value_type,DIM>::ConstElementIterator 
00064         the_element0 = f0.femSpace().beginElement(),
00065         the_element = fem_space.beginElement();
00066       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00067       for (;the_element != end_element;the_element ++, the_element0 ++) {
00068         double volume = the_element->templateElement().volume();
00069         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00070         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00071         const std::vector<int>& element_dof = the_element->dof();
00072         unsigned int n_element_dof = element_dof.size();
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<value_type> f0_value = f0.value(q_point, *the_element0);
00076         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00077         for (int l = 0;l < n_quadrature_point;l ++) {
00078           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00079           for (unsigned int j = 0;j < n_element_dof;j ++) {
00080             rhs(element_dof[j]) += Jxw*f0_value[l]*basis_value[j][l];
00081           }
00082         }
00083       }
00084       AMGSolver solver(mass_matrix);
00085       solver.solve(f1, rhs);
00086       break;
00087     }
00088     case LOCAL_LEAST_SQUARE: {
00089       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00090       unsigned int n_dof = fem_space.n_dof();
00091       std::vector<int> counter(n_dof, 0);
00092       f1.Vector<value_type>::operator=(0.0);
00093       typename FEMSpace<value_type,DIM>::ConstElementIterator 
00094         the_element0 = f0.femSpace().beginElement(),
00095         the_element = fem_space.beginElement();
00096       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00097       for (;the_element != end_element;the_element ++, the_element0 ++) {
00098         double volume = the_element->templateElement().volume();
00099         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00100         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00101         const std::vector<int>& element_dof = the_element->dof();
00102         unsigned int n_element_dof = element_dof.size();
00103         FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
00104         Vector<double> local_rhs(n_element_dof);
00105         Vector<double> local_f1(n_element_dof);
00106         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00107         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00108         std::vector<value_type> f0_value = f0.value(q_point, *the_element0);
00109         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00110         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00111           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00112           for (unsigned int j = 0;j < n_element_dof;j ++) {
00113             for (unsigned int k = 0;k < n_element_dof;k ++) {
00114               local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00115             }
00116             local_rhs(j) += Jxw*f0_value[l]*basis_value[j][l];
00117           }
00118         }
00119         local_mass_matrix.gauss_jordan();
00120         local_mass_matrix.vmult(local_f1, local_rhs);
00121         for (unsigned int i = 0;i < n_element_dof;i ++) {
00122           f1(element_dof[i]) += local_f1(i);
00123           counter[element_dof[i]] ++;
00124         }
00125       }
00126       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00127       break;
00128     }
00129     default: {
00130       Assert(false, ExcNotImplemented());
00131     }
00132     }
00133     return;
00134   }
00135   // otherwise, the two FEM functions are on different meshes
00136   switch (method) {
00137   case MASS_ACCUMULATION: {
00138     const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00139     FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00140     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00141     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00142     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00143     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00144     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00145       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00146                 << std::endl;
00147       Assert(false, ExcInternalError());
00148     }
00149     unsigned int n_dof = fem_space1.n_dof();
00150     Vector<double> mass_accumulation(n_dof);
00151     f1.Vector<value_type>::operator=(0.0);
00152     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00153     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00154     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00155     for (;the_pair != end_pair;++ the_pair) {
00156       const HElement<DIM>& h_element0 = the_pair(0);
00157       const HElement<DIM>& h_element1 = the_pair(1);
00158       const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00159       const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00160       Assert(element0.index() == h_element0.index, ExcInternalError());
00161       Assert(element1.index() == h_element1.index, ExcInternalError());
00162       const std::vector<int>& element_dof1 = element1.dof();
00163       unsigned int n_element_dof1 = element_dof1.size();
00164       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00165         double volume = element1.templateElement().volume();
00166         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00167         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00168         int n_quadrature_point = quad_info.n_quadraturePoint();
00169         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00170         std::vector<value_type> f0_value = f0.value(q_point, element0);
00171         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00172         for (int l = 0;l < n_quadrature_point;l ++) {
00173           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00174           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00175             f1(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00176             mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
00177           }
00178         }
00179       }
00180       else {
00181         double volume = element0.templateElement().volume();
00182         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00183         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00184         int n_quadrature_point = quad_info.n_quadraturePoint();
00185         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00186         std::vector<value_type> f0_value = f0.value(q_point, element0);
00187         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00188         for (int l = 0;l < n_quadrature_point;l ++) {
00189           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00190           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00191             f1(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00192             mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
00193           }
00194         }
00195       }
00196     }
00197     for (unsigned int i = 0;i < n_dof;i ++) 
00198       f1(i) /= mass_accumulation(i);
00199     break;
00200   }
00201   case LEAST_SQUARE: {
00202     const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00203     FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00204     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00205     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00206     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00207     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00208     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00209       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00210                 << std::endl;
00211       Assert(false, ExcInternalError());
00212     }
00213     f1.Vector<value_type>::operator=(0.0);
00214     MassMatrix<DIM, value_type> mass_matrix(fem_space1);
00215     mass_matrix.algebricAccuracy() = algebric_accuracy;
00216     mass_matrix.build();
00217     unsigned int n_dof = fem_space1.n_dof();
00218     Vector<double> rhs(n_dof);
00219     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00220     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00221     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00222     for (;the_pair != end_pair;++ the_pair) {
00223       const HElement<DIM>& h_element0 = the_pair(0);
00224       const HElement<DIM>& h_element1 = the_pair(1);
00225       const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00226       const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00227       Assert(element0.index() == h_element0.index, ExcInternalError());
00228       Assert(element1.index() == h_element1.index, ExcInternalError());
00229       const std::vector<int>& element_dof1 = element1.dof();
00230       unsigned int n_element_dof1 = element_dof1.size();
00231       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00232         double volume = element1.templateElement().volume();
00233         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00234         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00235         int n_quadrature_point = quad_info.n_quadraturePoint();
00236         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00237         std::vector<value_type> f0_value = f0.value(q_point, element0);
00238         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00239         for (int l = 0;l < n_quadrature_point;l ++) {
00240           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00241           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00242             rhs(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00243           }
00244         }
00245       }
00246       else {
00247         double volume = element0.templateElement().volume();
00248         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00249         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00250         int n_quadrature_point = quad_info.n_quadraturePoint();
00251         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00252         std::vector<value_type> f0_value = f0.value(q_point, element0);
00253         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00254         for (int l = 0;l < n_quadrature_point;l ++) {
00255           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00256           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00257             rhs(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00258           }
00259         }
00260       }
00261     }
00262     AMGSolver solver(mass_matrix);
00263     solver.solve(f1, rhs);
00264     break;
00265   }
00266   case LOCAL_LEAST_SQUARE: {
00267     const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00268     FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00269     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00270     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00271     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00272     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00273     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00274       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00275                 << std::endl;
00276       Assert(false, ExcInternalError());
00277     }
00278     unsigned int n_dof = fem_space1.n_dof();
00279     std::vector<int> counter(n_dof, 0);
00280     f1.Vector<value_type>::operator=(0.0);
00281     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00282     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00283     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00284     int the_element_index;
00285     unsigned int n_the_element_dof;
00286     { // only used to make the variables local
00287       const HElement<DIM>& the_h_element = the_pair(1);
00288       the_element_index = the_h_element.index;
00289       const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
00290       Assert(the_element.index() == the_element_index, ExcInternalError());
00291       n_the_element_dof = the_element.n_dof();
00292     }
00293     Vector<double> local_rhs(n_the_element_dof);
00294     for (;the_pair != end_pair;) {
00295       const HElement<DIM>& h_element0 = the_pair(0);
00296       const HElement<DIM>& h_element1 = the_pair(1);
00297       const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00298       const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00299       Assert(element0.index() == h_element0.index, ExcInternalError());
00300       Assert(element1.index() == h_element1.index, ExcInternalError());
00301       const std::vector<int>& element_dof1 = element1.dof();
00302       unsigned int n_element_dof1 = element_dof1.size();
00303       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00304         double volume = element1.templateElement().volume();
00305         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00306         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00307         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00308         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00309         std::vector<value_type> f0_value = f0.value(q_point, element0);
00310         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00311         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00312           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00313           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00314             local_rhs(j) += Jxw*f0_value[l]*basis_value1[j][l];
00315           }
00316         }
00317       }
00318       else {
00319         double volume = element0.templateElement().volume();
00320         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00321         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00322         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00323         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00324         std::vector<value_type> f0_value = f0.value(q_point, element0);
00325         std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00326         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00327           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00328           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00329             local_rhs(j) += Jxw*f0_value[l]*basis_value1[j][l];
00330           }
00331         }
00332       }
00333       ++ the_pair;
00334       if (the_pair == end_pair || the_pair(1).index != the_element_index) {
00335         const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
00336         const std::vector<int>& the_element_dof = the_element.dof();
00337         FullMatrix<double> local_mass_matrix(n_the_element_dof, n_the_element_dof);
00338         Vector<value_type> local_f1(n_the_element_dof);
00339         double volume = the_element.templateElement().volume();
00340         const QuadratureInfo<DIM>& quad_info = the_element.findQuadratureInfo(algebric_accuracy);
00341         std::vector<double> jacobian = the_element.local_to_global_jacobian(quad_info.quadraturePoint());
00342         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00343         std::vector<afepack::Point<DIM> > q_point = the_element.local_to_global(quad_info.quadraturePoint());
00344         std::vector<std::vector<value_type> > basis_value = the_element.basis_function_value(q_point);
00345         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00346           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00347           for (unsigned int j = 0;j < n_the_element_dof;j ++) {
00348             for (unsigned int k = 0;k < n_the_element_dof;k ++) {
00349               local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00350             }
00351           }
00352         }
00353         local_mass_matrix.gauss_jordan();
00354         local_mass_matrix.vmult(local_f1, local_rhs);
00355         for (unsigned int i = 0;i < n_the_element_dof;i ++) {
00356           f1(the_element_dof[i]) += local_f1(i);
00357           counter[the_element_dof[i]] ++;
00358         }
00359         if (the_pair != end_pair) {
00360           the_element_index = the_pair(1).index;
00361           Element<value_type,DIM>& the_new_element = fem_space1.element(the_element_index);
00362           n_the_element_dof = the_new_element.n_dof();
00363           local_rhs.reinit(n_the_element_dof);
00364         }
00365       }
00366     }
00367     for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00368     break;
00369   }
00370   default: {
00371     Assert(false, ExcNotImplemented());
00372   }
00373   }
00374 }
00375 
00376 
00377 
00378 template <class value_type, int DIM> 
00379   void Operator::L2Project(value_type (*f0)(const double *), 
00380                            FEMFunction<value_type, DIM>& f1, 
00381                            Method method, 
00382                            int algebric_accuracy)
00383 {
00384   switch (method) {
00385   case MASS_ACCUMULATION: {
00386     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00387     unsigned int n_dof = fem_space.n_dof();
00388     Vector<double> mass_accumulation(n_dof);
00389     f1.Vector<value_type>::operator=(0.0);
00390     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00391     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00392     for (;the_element != end_element;the_element ++) {
00393       double volume = the_element->templateElement().volume();
00394       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00395       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00396       const std::vector<int>& element_dof = the_element->dof();
00397       unsigned int n_element_dof = element_dof.size();
00398       int n_quadrature_point = quad_info.n_quadraturePoint();
00399       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00400       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00401       for (int l = 0;l < n_quadrature_point;l ++) {
00402         value_type f0_value = (*f0)(q_point[l]);
00403         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00404         for (unsigned int j = 0;j < n_element_dof;j ++) {
00405           f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00406           mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
00407         }
00408       }
00409     }
00410     for (unsigned int i = 0;i < n_dof;i ++) 
00411       f1(i) /= mass_accumulation(i);
00412     break;
00413   }
00414   case LEAST_SQUARE: {
00415     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00416     f1.Vector<value_type>::operator=(0.0);
00417     unsigned int n_dof = fem_space.n_dof();
00418     MassMatrix<DIM, value_type> mass_matrix(fem_space);
00419     mass_matrix.algebricAccuracy() = algebric_accuracy;
00420     mass_matrix.build();
00421     Vector<double> rhs(n_dof);
00422     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00423     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00424     for (;the_element != end_element;the_element ++) {
00425       double volume = the_element->templateElement().volume();
00426       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00427       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00428       const std::vector<int>& element_dof = the_element->dof();
00429       unsigned int n_element_dof = element_dof.size();
00430       int n_quadrature_point = quad_info.n_quadraturePoint();
00431       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00432       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00433       for (int l = 0;l < n_quadrature_point;l ++) {
00434         value_type f0_value = (*f0)(q_point[l]);
00435         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00436         for (unsigned int j = 0;j < n_element_dof;j ++) {
00437           rhs(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00438         }
00439       }
00440     }
00441     AMGSolver solver(mass_matrix);
00442     solver.solve(f1, rhs);
00443     break;
00444   }
00445   case LOCAL_LEAST_SQUARE: {
00446     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00447     unsigned int n_dof = fem_space.n_dof();
00448     std::vector<int> counter(n_dof, 0);
00449     f1.Vector<value_type>::operator=(0.0);
00450     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00451     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00452     for (;the_element != end_element;the_element ++) {
00453       double volume = the_element->templateElement().volume();
00454       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00455       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00456       const std::vector<int>& element_dof = the_element->dof();
00457       unsigned int n_element_dof = element_dof.size();
00458       FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
00459       Vector<double> local_rhs(n_element_dof);
00460       Vector<double> local_f1(n_element_dof);
00461       unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00462       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00463       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00464       for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00465         value_type f0_value = (*f0)(q_point[l]);
00466         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00467         for (unsigned int j = 0;j < n_element_dof;j ++) {
00468           for (unsigned int k = 0;k < n_element_dof;k ++) {
00469             local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00470           }
00471           local_rhs(j) += Jxw*f0_value*basis_value[j][l];
00472         }
00473       }
00474       local_mass_matrix.gauss_jordan();
00475       local_mass_matrix.vmult(local_f1, local_rhs);
00476       for (unsigned int i = 0;i < n_element_dof;i ++) {
00477         f1(element_dof[i]) += local_f1(i);
00478         counter[element_dof[i]] ++;
00479       }
00480     }
00481     for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00482     break;
00483   }
00484   default: {
00485     Assert(false, ExcNotImplemented());
00486   }
00487   }
00488 }
00489 
00490 
00491 
00492 template <class value_type, int DIM> 
00493   void Operator::L2Project(value_type (*f0)(const afepack::Point<DIM>&), 
00494                            FEMFunction<value_type, DIM>& f1, 
00495                            Method method, 
00496                            int algebric_accuracy)
00497 {
00498   switch (method) {
00499   case MASS_ACCUMULATION: {
00500     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00501     unsigned int n_dof = fem_space.n_dof();
00502     Vector<double> mass_accumulation(n_dof);
00503     f1.Vector<value_type>::operator=(0.0);
00504     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00505     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00506     for (;the_element != end_element;the_element ++) {
00507       double volume = the_element->templateElement().volume();
00508       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00509       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00510       const std::vector<int>& element_dof = the_element->dof();
00511       unsigned int n_element_dof = element_dof.size();
00512       int n_quadrature_point = quad_info.n_quadraturePoint();
00513       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00514       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00515       for (int l = 0;l < n_quadrature_point;l ++) {
00516         value_type f0_value = (*f0)(q_point[l]);
00517         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00518         for (unsigned int j = 0;j < n_element_dof;j ++) {
00519           f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00520           mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
00521         }
00522       }
00523     }
00524     for (unsigned int i = 0;i < n_dof;i ++) 
00525       f1(i) /= mass_accumulation(i);
00526     break;
00527   }
00528   case LEAST_SQUARE: {
00529     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00530     f1.Vector<value_type>::operator=(0.0);
00531     unsigned int n_dof = fem_space.n_dof();
00532     MassMatrix<DIM, value_type> mass_matrix(fem_space);
00533     mass_matrix.algebricAccuracy() = algebric_accuracy;
00534     mass_matrix.build();
00535     Vector<double> rhs(n_dof);
00536     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00537     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00538     for (;the_element != end_element;the_element ++) {
00539       double volume = the_element->templateElement().volume();
00540       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00541       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00542       const std::vector<int>& element_dof = the_element->dof();
00543       unsigned int n_element_dof = element_dof.size();
00544       int n_quadrature_point = quad_info.n_quadraturePoint();
00545       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00546       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00547       for (int l = 0;l < n_quadrature_point;l ++) {
00548         value_type f0_value = (*f0)(q_point[l]);
00549         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00550         for (unsigned int j = 0;j < n_element_dof;j ++) {
00551           rhs(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00552         }
00553       }
00554     }
00555     AMGSolver solver(mass_matrix);
00556     solver.solve(f1, rhs);
00557     break;
00558   }
00559   case LOCAL_LEAST_SQUARE: {
00560     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00561     unsigned int n_dof = fem_space.n_dof();
00562     std::vector<int> counter(n_dof, 0);
00563     f1.Vector<value_type>::operator=(0.0);
00564     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00565     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00566     for (;the_element != end_element;the_element ++) {
00567       double volume = the_element->templateElement().volume();
00568       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00569       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00570       const std::vector<int>& element_dof = the_element->dof();
00571       unsigned int n_element_dof = element_dof.size();
00572       FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
00573       Vector<double> local_rhs(n_element_dof);
00574       Vector<double> local_f1(n_element_dof);
00575       unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00576       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00577       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00578       for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00579         value_type f0_value = (*f0)(q_point[l]);
00580         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00581         for (unsigned int j = 0;j < n_element_dof;j ++) {
00582           for (unsigned int k = 0;k < n_element_dof;k ++) {
00583             local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00584           }
00585           local_rhs(j) += Jxw*f0_value*basis_value[j][l];
00586         }
00587       }
00588       local_mass_matrix.gauss_jordan();
00589       local_mass_matrix.vmult(local_f1, local_rhs);
00590       for (unsigned int i = 0;i < n_element_dof;i ++) {
00591         f1(element_dof[i]) += local_f1(i);
00592         counter[element_dof[i]] ++;
00593       }
00594     }
00595     for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00596     break;
00597   }
00598   default: {
00599     Assert(false, ExcNotImplemented());
00600   }
00601   }
00602 }
00603 
00604 
00605 
00606 template <class value_type, int DIM> 
00607   void Operator::L2Project(const Function<value_type>& f0, 
00608                            FEMFunction<value_type, DIM>& f1, 
00609                            Method method, 
00610                            int algebric_accuracy)
00611 {
00612   switch (method) {
00613   case MASS_ACCUMULATION: {
00614     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00615     unsigned int n_dof = fem_space.n_dof();
00616     Vector<double> mass_accumulation(n_dof);
00617     f1.Vector<value_type>::operator=(0.0);
00618     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00619     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00620     for (;the_element != end_element;the_element ++) {
00621       double volume = the_element->templateElement().volume();
00622       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00623       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00624       const std::vector<int>& element_dof = the_element->dof();
00625       unsigned int n_element_dof = element_dof.size();
00626       int n_quadrature_point = quad_info.n_quadraturePoint();
00627       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00628       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00629       for (int l = 0;l < n_quadrature_point;l ++) {
00630         value_type f0_value = f0.value(q_point[l]);
00631         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00632         for (unsigned int j = 0;j < n_element_dof;j ++) {
00633           f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00634           mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
00635         }
00636       }
00637     }
00638     for (unsigned int i = 0;i < n_dof;i ++) 
00639       f1(i) /= mass_accumulation(i);
00640     break;
00641   }
00642   case LEAST_SQUARE: {
00643     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00644     f1.Vector<value_type>::operator=(0.0);
00645     unsigned int n_dof = fem_space.n_dof();
00646     MassMatrix<DIM, value_type> mass_matrix(fem_space);
00647     mass_matrix.algebricAccuracy() = algebric_accuracy;
00648     mass_matrix.build();
00649     Vector<double> rhs(n_dof);
00650     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00651     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00652     for (;the_element != end_element;the_element ++) {
00653       double volume = the_element->templateElement().volume();
00654       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00655       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00656       const std::vector<int>& element_dof = the_element->dof();
00657       unsigned int n_element_dof = element_dof.size();
00658       int n_quadrature_point = quad_info.n_quadraturePoint();
00659       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00660       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00661       for (int l = 0;l < n_quadrature_point;l ++) {
00662         value_type f0_value = f0.value(q_point[l]);
00663         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00664         for (unsigned int j = 0;j < n_element_dof;j ++) {
00665           rhs(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00666         }
00667       }
00668     }
00669     AMGSolver solver(mass_matrix);
00670     solver.solve(f1, rhs);
00671     break;
00672   }
00673   case LOCAL_LEAST_SQUARE: {
00674     FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00675     unsigned int n_dof = fem_space.n_dof();
00676     std::vector<int> counter(n_dof, 0);
00677     f1.Vector<value_type>::operator=(0.0);
00678     typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00679     typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00680     for (;the_element != end_element;the_element ++) {
00681       double volume = the_element->templateElement().volume();
00682       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00683       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00684       const std::vector<int>& element_dof = the_element->dof();
00685       unsigned int n_element_dof = element_dof.size();
00686       FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
00687       Vector<double> local_rhs(n_element_dof);
00688       Vector<double> local_f1(n_element_dof);
00689       unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00690       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00691       std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00692       for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00693         value_type f0_value = f0.value(q_point[l]);
00694         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00695         for (unsigned int j = 0;j < n_element_dof;j ++) {
00696           for (unsigned int k = 0;k < n_element_dof;k ++) {
00697             local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00698           }
00699           local_rhs(j) += Jxw*f0_value*basis_value[j][l];
00700         }
00701       }
00702       local_mass_matrix.gauss_jordan();
00703       local_mass_matrix.vmult(local_f1, local_rhs);
00704       for (unsigned int i = 0;i < n_element_dof;i ++) {
00705         f1(element_dof[i]) += local_f1(i);
00706         counter[element_dof[i]] ++;
00707       }
00708     }
00709     for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00710     break;
00711   }
00712   default: {
00713     Assert(false, ExcNotImplemented());
00714   }
00715   }
00716 }
00717 
00718 
00719 
00720 template <class value_type, int DIM>
00721   void Operator::L2Project(value_type (*f)(const value_type&), 
00722                            const FEMFunction<value_type,DIM>& f0,
00723                            FEMFunction<value_type,DIM>& f1, 
00724                            Method method, 
00725                            int algebric_accuracy)
00726 {
00727   if (&f0.femSpace() == &f1.femSpace()) {
00728     switch (method) {
00729     case MASS_ACCUMULATION: {
00730       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00731       unsigned int n_dof = fem_space.n_dof();
00732       Vector<double> mass_accumulation(n_dof);
00733       f1.Vector<value_type>::operator=(0.0);
00734       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00735       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00736       for (;the_element != end_element;the_element ++) {
00737         double volume = the_element->templateElement().volume();
00738         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00739         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00740         const std::vector<int>& element_dof = the_element->dof();
00741         unsigned int n_element_dof = element_dof.size();
00742         int n_quadrature_point = quad_info.n_quadraturePoint();
00743         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00744         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00745         std::vector<value_type> f0_value = f0.value(q_point, *the_element);
00746         for (int l = 0;l < n_quadrature_point;l ++) {
00747           value_type f_value = (*f)(f0_value[l]);
00748           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00749           for (unsigned int j = 0;j < n_element_dof;j ++) {
00750             f1(element_dof[j]) += Jxw*f_value*basis_value[j][l];
00751             mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
00752           }
00753         }
00754       }
00755       for (unsigned int i = 0;i < n_dof;i ++)
00756         f1(i) /= mass_accumulation(i);
00757       break;
00758     }
00759     case LEAST_SQUARE: {
00760       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00761       unsigned int n_dof = fem_space.n_dof();
00762       f1.Vector<value_type>::operator=(0.0);
00763       MassMatrix<DIM, value_type> mass_matrix(fem_space);
00764       mass_matrix.algebricAccuracy() = algebric_accuracy;
00765       mass_matrix.build();
00766       Vector<double> rhs(n_dof);
00767       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00768       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00769       for (;the_element != end_element;the_element ++) {
00770         double volume = the_element->templateElement().volume();
00771         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00772         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00773         const std::vector<int>& element_dof = the_element->dof();
00774         unsigned int n_element_dof = element_dof.size();
00775         int n_quadrature_point = quad_info.n_quadraturePoint();
00776         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00777         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00778         std::vector<value_type> f0_value = f0.value(q_point, *the_element);
00779         for (int l = 0;l < n_quadrature_point;l ++) {
00780           value_type f_value = (*f)(f0_value[l]);
00781           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00782           for (unsigned int j = 0;j < n_element_dof;j ++) {
00783             rhs(element_dof[j]) += Jxw*f_value*basis_value[j][l];
00784           }
00785         }
00786       }
00787       AMGSolver solver(mass_matrix);
00788       solver.solve(f1, rhs);
00789       break;
00790     }
00791     case LOCAL_LEAST_SQUARE: {
00792       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00793       unsigned int n_dof = fem_space.n_dof();
00794       std::vector<int> counter(n_dof, 0);
00795       f1.Vector<value_type>::operator=(0.0);
00796       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00797       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00798       for (;the_element != end_element;the_element ++) {
00799         double volume = the_element->templateElement().volume();
00800         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00801         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00802         const std::vector<int>& element_dof = the_element->dof();
00803         unsigned int n_element_dof = element_dof.size();
00804         FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
00805         Vector<double> local_rhs(n_element_dof);
00806         Vector<double> local_f1(n_element_dof);
00807         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
00808         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00809         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
00810         std::vector<value_type> f0_value = f0.value(q_point, *the_element);
00811         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
00812           value_type f_value = (*f)(f0_value[l]);
00813           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00814           for (unsigned int j = 0;j < n_element_dof;j ++) {
00815             for (unsigned int k = 0;k < n_element_dof;k ++) {
00816               local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
00817             }
00818             local_rhs(j) += Jxw*f_value*basis_value[j][l];
00819           }
00820         }
00821         local_mass_matrix.gauss_jordan();
00822         local_mass_matrix.vmult(local_f1, local_rhs);
00823         for (unsigned int i = 0;i < n_element_dof;i ++) {
00824           f1(element_dof[i]) += local_f1(i);
00825           counter[element_dof[i]] ++;
00826         }
00827       }
00828       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
00829       break;
00830     }
00831     default: {
00832       Assert(false, ExcNotImplemented());
00833     }
00834     }
00835   }
00836   else {
00837     switch (method) {
00838     case MASS_ACCUMULATION: {
00839       const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00840       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00841       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00842       const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00843       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00844       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00845       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00846         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00847                   << std::endl;
00848         Assert(false, ExcInternalError());
00849       }
00850       unsigned int n_dof = fem_space1.n_dof();
00851       Vector<double> mass_accumulation(n_dof);
00852       f1.Vector<value_type>::operator=(0.0);
00853       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00854       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00855       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00856       for (;the_pair != end_pair;++ the_pair) {
00857         const HElement<DIM>& h_element0 = the_pair(0);
00858         const HElement<DIM>& h_element1 = the_pair(1);
00859         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00860         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00861         Assert(element0.index() == h_element0.index, ExcInternalError());
00862         Assert(element1.index() == h_element1.index, ExcInternalError());
00863         const std::vector<int>& element_dof1 = element1.dof();
00864         unsigned int n_element_dof1 = element_dof1.size();
00865         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00866           double volume = element1.templateElement().volume();
00867           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00868           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00869           int n_quadrature_point = quad_info.n_quadraturePoint();
00870           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00871           std::vector<value_type> f0_value = f0.value(q_point, element0);
00872           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00873           for (int l = 0;l < n_quadrature_point;l ++) {
00874             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00875             value_type f_value = (*f)(f0_value[l]);
00876             for (unsigned int j = 0;j < n_element_dof1;j ++) {
00877               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00878               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
00879             }
00880           }
00881         }
00882         else {
00883           double volume = element0.templateElement().volume();
00884           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00885           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00886           int n_quadrature_point = quad_info.n_quadraturePoint();
00887           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00888           std::vector<value_type> f0_value = f0.value(q_point, element0);
00889           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00890           for (int l = 0;l < n_quadrature_point;l ++) {
00891             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00892             value_type f_value = (*f)(f0_value[l]);
00893             for (unsigned int j = 0;j < n_element_dof1;j ++) {
00894               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00895               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
00896             }
00897           }
00898         }
00899       }
00900       for (unsigned int i = 0;i < n_dof;i ++)
00901         f1(i) /= mass_accumulation(i);
00902       break;
00903     }
00904     case LEAST_SQUARE: {
00905       const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00906       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00907       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00908       const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00909       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00910       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00911       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00912         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00913                   << std::endl;
00914         Assert(false, ExcInternalError());
00915       }
00916       MassMatrix<DIM, value_type> mass_matrix(fem_space1);
00917       mass_matrix.algebricAccuracy() = algebric_accuracy;
00918       mass_matrix.build();
00919       unsigned int n_dof = fem_space1.n_dof();
00920       Vector<double> rhs(n_dof);
00921       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00922       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00923       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00924       for (;the_pair != end_pair;++ the_pair) {
00925         const HElement<DIM>& h_element0 = the_pair(0);
00926         const HElement<DIM>& h_element1 = the_pair(1);
00927         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00928         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00929         Assert(element0.index() == h_element0.index, ExcInternalError());
00930         Assert(element1.index() == h_element1.index, ExcInternalError());
00931         const std::vector<int>& element_dof1 = element1.dof();
00932         unsigned int n_element_dof1 = element_dof1.size();
00933         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00934           double volume = element1.templateElement().volume();
00935           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00936           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00937           int n_quadrature_point = quad_info.n_quadraturePoint();
00938           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00939           std::vector<value_type> f0_value = f0.value(q_point, element0);
00940           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00941           for (int l = 0;l < n_quadrature_point;l ++) {
00942             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00943             value_type f_value = (*f)(f0_value[l]);
00944             for (unsigned int j = 0;j < n_element_dof1;j ++) {
00945               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00946             }
00947           }
00948         }
00949         else {
00950           double volume = element0.templateElement().volume();
00951           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00952           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00953           int n_quadrature_point = quad_info.n_quadraturePoint();
00954           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00955           std::vector<value_type> f0_value = f0.value(q_point, element0);
00956           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
00957           for (int l = 0;l < n_quadrature_point;l ++) {
00958             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00959             value_type f_value = (*f)(f0_value[l]);
00960             for (unsigned int j = 0;j < n_element_dof1;j ++) {
00961               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00962             }
00963           }
00964         }
00965       }
00966       AMGSolver solver(mass_matrix);
00967       solver.solve(f1, rhs);
00968       break;
00969     }
00970     case LOCAL_LEAST_SQUARE: {
00971       const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00972       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00973       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00974       const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00975       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00976       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00977       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00978         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00979                   << std::endl;
00980         Assert(false, ExcInternalError());
00981       }
00982       unsigned int n_dof = fem_space1.n_dof();
00983       std::vector<int> counter(n_dof, 0);
00984       f1.Vector<value_type>::operator=(0.0);
00985       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00986       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00987       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00988       int the_element_index;
00989       unsigned int n_the_element_dof;
00990       { // only used to make the variables local
00991         const HElement<DIM>& the_h_element = the_pair(1);
00992         the_element_index = the_h_element.index;
00993         const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
00994         Assert(the_element.index() == the_element_index, ExcInternalError());
00995         n_the_element_dof = the_element.n_dof();
00996       }
00997       Vector<double> local_rhs(n_the_element_dof);
00998       for (;the_pair != end_pair;) {
00999         const HElement<DIM>& h_element0 = the_pair(0);
01000         const HElement<DIM>& h_element1 = the_pair(1);
01001         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01002         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01003         Assert(element0.index() == h_element0.index, ExcInternalError());
01004         Assert(element1.index() == h_element1.index, ExcInternalError());
01005         const std::vector<int>& element_dof1 = element1.dof();
01006         unsigned int n_element_dof1 = element_dof1.size();
01007         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01008           double volume = element1.templateElement().volume();
01009           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01010           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01011           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01012           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01013           std::vector<value_type> f0_value = f0.value(q_point, element0);
01014           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01015           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01016             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01017             value_type f_value = (*f)(f0_value[l]);
01018             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01019               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01020             }
01021           }
01022         }
01023         else {
01024           double volume = element0.templateElement().volume();
01025           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01026           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01027           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01028           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01029           std::vector<value_type> f0_value = f0.value(q_point, element0);
01030           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01031           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01032             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01033             value_type f_value = (*f)(f0_value[l]);
01034             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01035               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01036             }
01037           }
01038         }
01039         ++ the_pair;
01040         if (the_pair == end_pair || the_pair(1).index != the_element_index) {
01041           const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01042           const std::vector<int>& the_element_dof = the_element.dof();
01043           FullMatrix<double> local_mass_matrix(n_the_element_dof, n_the_element_dof);
01044           Vector<value_type> local_f1(n_the_element_dof);
01045           double volume = element1.templateElement().volume();
01046           const QuadratureInfo<DIM>& quad_info = the_element.findQuadratureInfo(algebric_accuracy);
01047           std::vector<double> jacobian = the_element.local_to_global_jacobian(quad_info.quadraturePoint());
01048           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01049           std::vector<afepack::Point<DIM> > q_point = the_element.local_to_global(quad_info.quadraturePoint());
01050           std::vector<std::vector<value_type> > basis_value = the_element.basis_function_value(q_point);
01051           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01052             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01053             for (unsigned int j = 0;j < n_the_element_dof;j ++) {
01054               for (unsigned int k = 0;k < n_the_element_dof;k ++) {
01055                 local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
01056               }
01057             }
01058           }
01059           local_mass_matrix.gauss_jordan();
01060           local_mass_matrix.vmult(local_f1, local_rhs);
01061           for (unsigned int i = 0;i < n_the_element_dof;i ++) {
01062             f1(the_element_dof[i]) += local_f1(i);
01063             counter[the_element_dof[i]] ++;
01064           }
01065           if (the_pair != end_pair) {
01066             the_element_index = the_pair(1).index;
01067             const Element<value_type,DIM>& the_new_element = fem_space1.element(the_element_index);
01068             n_the_element_dof = the_new_element.n_dof();
01069             local_rhs.reinit(n_the_element_dof);
01070           }
01071         }
01072       }
01073       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
01074       break;
01075     }
01076     default: {
01077       Assert(false, ExcNotImplemented());
01078     }
01079     }
01080   }
01081 }
01082 
01083 
01084 
01085 template <class value_type, int DIM>
01086   void Operator::L2Project(value_type (*f)(const value_type&,const value_type&), 
01087                            const FEMFunction<value_type,DIM>& f00, 
01088                            const FEMFunction<value_type,DIM>& f01,
01089                            FEMFunction<value_type,DIM>& f1, 
01090                            Method method, 
01091                            int algebric_accuracy)
01092 {
01093   if ((&f00.femSpace() != &f01.femSpace())
01094       && (&f00.femSpace() != &f1.femSpace())
01095       && (&f01.femSpace() != &f1.femSpace())) {
01096     std::cerr << "The three FEM functions are on three different finite element spaces."
01097               << std::endl;
01098     abort();
01099   }
01100   else if ((&f00.femSpace() == &f01.femSpace()) &&
01101            (&f00.femSpace() == &f1.femSpace())) {
01102     switch (method) {
01103     case MASS_ACCUMULATION: {
01104       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
01105       unsigned int n_dof = fem_space.n_dof();
01106       Vector<double> mass_accumulation(n_dof);
01107       f1.Vector<value_type>::operator=(0.0);
01108       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
01109       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
01110       for (;the_element != end_element;the_element ++) {
01111         double volume = the_element->templateElement().volume();
01112         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
01113         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
01114         const std::vector<int>& element_dof = the_element->dof();
01115         unsigned int n_element_dof = element_dof.size();
01116         int n_quadrature_point = quad_info.n_quadraturePoint();
01117         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
01118         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
01119         std::vector<value_type> f00_value = f00.value(q_point, *the_element);
01120         std::vector<value_type> f01_value = f01.value(q_point, *the_element);
01121         for (int l = 0;l < n_quadrature_point;l ++) {
01122           value_type f_value = (*f)(f00_value[l], f01_value[l]);
01123           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01124           for (unsigned int j = 0;j < n_element_dof;j ++) {
01125             f1(element_dof[j]) += Jxw*f_value*basis_value[j][l];
01126             mass_accumulation(element_dof[j]) += Jxw*basis_value[j][l];
01127           }
01128         }
01129       }
01130       for (unsigned int i = 0;i < n_dof;i ++)
01131         f1(i) /= mass_accumulation(i);
01132       break;
01133     }
01134     case LEAST_SQUARE: {
01135       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
01136       unsigned int n_dof = fem_space.n_dof();
01137       f1.Vector<value_type>::operator=(0.0);
01138       MassMatrix<DIM, value_type> mass_matrix(fem_space);
01139       mass_matrix.algebricAccuracy() = algebric_accuracy;
01140       mass_matrix.build();
01141       Vector<double> rhs(n_dof);
01142       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
01143       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
01144       for (;the_element != end_element;the_element ++) {
01145         double volume = the_element->templateElement().volume();
01146         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
01147         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
01148         const std::vector<int>& element_dof = the_element->dof();
01149         unsigned int n_element_dof = element_dof.size();
01150         int n_quadrature_point = quad_info.n_quadraturePoint();
01151         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
01152         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
01153         std::vector<value_type> f00_value = f00.value(q_point, *the_element);
01154         std::vector<value_type> f01_value = f01.value(q_point, *the_element);
01155         for (int l = 0;l < n_quadrature_point;l ++) {
01156           value_type f_value = (*f)(f00_value[l], f01_value[l]);
01157           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01158           for (unsigned int j = 0;j < n_element_dof;j ++) {
01159             rhs(element_dof[j]) += Jxw*f_value*basis_value[j][l];
01160           }
01161         }
01162       }
01163       AMGSolver solver(mass_matrix);
01164       solver.solve(f1, rhs);
01165       break;
01166     }
01167     case LOCAL_LEAST_SQUARE: {
01168       FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
01169       unsigned int n_dof = fem_space.n_dof();
01170       std::vector<int> counter(n_dof, 0);
01171       f1.Vector<value_type>::operator=(0.0);
01172       typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
01173       typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
01174       for (;the_element != end_element;the_element ++) {
01175         double volume = the_element->templateElement().volume();
01176         const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
01177         std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
01178         const std::vector<int>& element_dof = the_element->dof();
01179         unsigned int n_element_dof = element_dof.size();
01180         FullMatrix<double> local_mass_matrix(n_element_dof, n_element_dof);
01181         Vector<double> local_rhs(n_element_dof);
01182         Vector<double> local_f1(n_element_dof);
01183         unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01184         std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
01185         std::vector<std::vector<value_type> > basis_value = the_element->basis_function_value(q_point);
01186         std::vector<value_type> f00_value = f00.value(q_point, *the_element);
01187         std::vector<value_type> f01_value = f01.value(q_point, *the_element);
01188         for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01189           value_type f_value = (*f)(f00_value[l], f01_value[l]);
01190           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01191           for (unsigned int j = 0;j < n_element_dof;j ++) {
01192             for (unsigned int k = 0;k < n_element_dof;k ++) {
01193               local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
01194             }
01195             local_rhs(j) += Jxw*f_value*basis_value[j][l];
01196           }
01197         }
01198         local_mass_matrix.gauss_jordan();
01199         local_mass_matrix.vmult(local_f1, local_rhs);
01200         for (unsigned int i = 0;i < n_element_dof;i ++) {
01201           f1(element_dof[i]) += local_f1(i);
01202           counter[element_dof[i]] ++;
01203         }
01204       }
01205       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
01206       break;
01207     }
01208     default: {
01209       Assert(false, ExcNotImplemented());
01210     }
01211     }
01212   }
01213   else if ((&f00.femSpace() == &f01.femSpace())) {
01214     switch (method) {
01215     case MASS_ACCUMULATION: {
01216       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01217       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01218       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01219       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01220       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01221       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01222       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01223         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01224                   << std::endl;
01225         Assert(false, ExcInternalError());
01226       }
01227       unsigned int n_dof = fem_space1.n_dof();
01228       Vector<double> mass_accumulation(n_dof);
01229       f1.Vector<value_type>::operator=(0.0);
01230       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01231       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01232       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01233       for (;the_pair != end_pair;++ the_pair) {
01234         const HElement<DIM>& h_element0 = the_pair(0);
01235         const HElement<DIM>& h_element1 = the_pair(1);
01236         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01237         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01238         Assert(element0.index() == h_element0.index, ExcInternalError());
01239         Assert(element1.index() == h_element1.index, ExcInternalError());
01240         const std::vector<int>& element_dof1 = element1.dof();
01241         unsigned int n_element_dof1 = element_dof1.size();
01242         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01243           double volume = element1.templateElement().volume();
01244           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01245           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01246           int n_quadrature_point = quad_info.n_quadraturePoint();
01247           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01248           std::vector<value_type> f00_value = f00.value(q_point, element0);
01249           std::vector<value_type> f01_value = f01.value(q_point, element0);
01250           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01251           for (int l = 0;l < n_quadrature_point;l ++) {
01252             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01253             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01254             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01255               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01256               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01257             }
01258           }
01259         }
01260         else {
01261           double volume = element0.templateElement().volume();
01262           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01263           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01264           int n_quadrature_point = quad_info.n_quadraturePoint();
01265           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01266           std::vector<value_type> f00_value = f00.value(q_point, element0);
01267           std::vector<value_type> f01_value = f01.value(q_point, element0);
01268           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01269           for (int l = 0;l < n_quadrature_point;l ++) {
01270             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01271             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01272             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01273               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01274               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01275             }
01276           }
01277         }
01278       }
01279       for (unsigned int i = 0;i < n_dof;i ++)
01280         f1(i) /= mass_accumulation(i);
01281       break;
01282     }
01283     case LEAST_SQUARE: {
01284       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01285       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01286       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01287       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01288       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01289       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01290       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01291         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01292                   << std::endl;
01293         Assert(false, ExcInternalError());
01294       }
01295       f1.Vector<value_type>::operator=(0.0);
01296       MassMatrix<DIM, value_type> mass_matrix(fem_space1);
01297       mass_matrix.algebricAccuracy() = algebric_accuracy;
01298       mass_matrix.build();
01299       unsigned int n_dof = fem_space1.n_dof();
01300       Vector<double> rhs(n_dof);
01301       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01302       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01303       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01304       for (;the_pair != end_pair;++ the_pair) {
01305         const HElement<DIM>& h_element0 = the_pair(0);
01306         const HElement<DIM>& h_element1 = the_pair(1);
01307         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01308         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01309         Assert(element0.index() == h_element0.index, ExcInternalError());
01310         Assert(element1.index() == h_element1.index, ExcInternalError());
01311         const std::vector<int>& element_dof1 = element1.dof();
01312         unsigned int n_element_dof1 = element_dof1.size();
01313         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01314           double volume = element1.templateElement().volume();
01315           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01316           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01317           int n_quadrature_point = quad_info.n_quadraturePoint();
01318           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01319           std::vector<value_type> f00_value = f00.value(q_point, element0);
01320           std::vector<value_type> f01_value = f01.value(q_point, element0);
01321           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01322           for (int l = 0;l < n_quadrature_point;l ++) {
01323             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01324             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01325             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01326               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01327             }
01328           }
01329         }
01330         else {
01331           double volume = element0.templateElement().volume();
01332           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01333           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01334           int n_quadrature_point = quad_info.n_quadraturePoint();
01335           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01336           std::vector<value_type> f00_value = f00.value(q_point, element0);
01337           std::vector<value_type> f01_value = f01.value(q_point, element0);
01338           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01339           for (int l = 0;l < n_quadrature_point;l ++) {
01340             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01341             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01342             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01343               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01344             }
01345           }
01346         }
01347       }
01348       AMGSolver solver(mass_matrix);
01349       solver.solve(f1, rhs);
01350       break;
01351     }
01352     case LOCAL_LEAST_SQUARE: {
01353       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01354       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01355       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01356       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01357       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01358       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01359       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01360         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01361                   << std::endl;
01362         Assert(false, ExcInternalError());
01363       }
01364       unsigned int n_dof = fem_space1.n_dof();
01365       std::vector<int> counter(n_dof, 0);
01366       f1.Vector<value_type>::operator=(0.0);
01367       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01368       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01369       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01370       int the_element_index;
01371       unsigned int n_the_element_dof;
01372       { // only used to make the variables local
01373         const HElement<DIM>& the_h_element = the_pair(1);
01374         the_element_index = the_h_element.index;
01375         const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01376         Assert(the_element.index() == the_element_index, ExcInternalError());
01377         n_the_element_dof = the_element.n_dof();
01378       }
01379       Vector<double> local_rhs(n_the_element_dof);
01380       for (;the_pair != end_pair;) {
01381         const HElement<DIM>& h_element0 = the_pair(0);
01382         const HElement<DIM>& h_element1 = the_pair(1);
01383         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01384         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01385         Assert(element0.index() == h_element0.index, ExcInternalError());
01386         Assert(element1.index() == h_element1.index, ExcInternalError());
01387         const std::vector<int>& element_dof1 = element1.dof();
01388         unsigned int n_element_dof1 = element_dof1.size();
01389         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01390           double volume = element1.templateElement().volume();
01391           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01392           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01393           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01394           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01395           std::vector<value_type> f00_value = f00.value(q_point, element0);
01396           std::vector<value_type> f01_value = f01.value(q_point, element0);
01397           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01398           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01399             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01400             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01401             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01402               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01403             }
01404           }
01405         }
01406         else {
01407           double volume = element0.templateElement().volume();
01408           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01409           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01410           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01411           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01412           std::vector<value_type> f00_value = f00.value(q_point, element0);
01413           std::vector<value_type> f01_value = f01.value(q_point, element0);
01414           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01415           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01416             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01417             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01418             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01419               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01420             }
01421           }
01422         }
01423         ++ the_pair;
01424         if (the_pair == end_pair || the_pair(1).index != the_element_index) {
01425           const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01426           const std::vector<int>& the_element_dof = the_element.dof();
01427           FullMatrix<double> local_mass_matrix(n_the_element_dof, n_the_element_dof);
01428           Vector<value_type> local_f1(n_the_element_dof);
01429           double volume = element1.templateElement().volume();
01430           const QuadratureInfo<DIM>& quad_info = the_element.findQuadratureInfo(algebric_accuracy);
01431           std::vector<double> jacobian = the_element.local_to_global_jacobian(quad_info.quadraturePoint());
01432           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01433           std::vector<afepack::Point<DIM> > q_point = the_element.local_to_global(quad_info.quadraturePoint());
01434           std::vector<std::vector<value_type> > basis_value = the_element.basis_function_value(q_point);
01435           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01436             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01437             for (unsigned int j = 0;j < n_the_element_dof;j ++) {
01438               for (unsigned int k = 0;k < n_the_element_dof;k ++) {
01439                 local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
01440               }
01441             }
01442           }
01443           local_mass_matrix.gauss_jordan();
01444           local_mass_matrix.vmult(local_f1, local_rhs);
01445           for (unsigned int i = 0;i < n_the_element_dof;i ++) {
01446             f1(the_element_dof[i]) += local_f1(i);
01447             counter[the_element_dof[i]] ++;
01448           }
01449           if (the_pair != end_pair) {
01450             the_element_index = the_pair(1).index;
01451             Element<value_type,DIM>& the_new_element = fem_space1.element(the_element_index);
01452             n_the_element_dof = the_new_element.n_dof();
01453             local_rhs.reinit(n_the_element_dof);
01454           }
01455         }
01456       }
01457       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
01458       break;
01459     }
01460     default: {
01461       Assert(false, ExcNotImplemented());
01462     }
01463     }
01464   }
01465   else if (&f00.femSpace() == &f1.femSpace()) {
01466     switch (method) {
01467     case MASS_ACCUMULATION: {
01468       const FEMSpace<value_type,DIM>& fem_space0 = f01.femSpace();
01469       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01470       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01471       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01472       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01473       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01474       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01475         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01476                   << std::endl;
01477         Assert(false, ExcInternalError());
01478       }
01479       unsigned int n_dof = fem_space1.n_dof();
01480       Vector<double> mass_accumulation(n_dof);
01481       f1.Vector<value_type>::operator=(0.0);
01482       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01483       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01484       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01485       for (;the_pair != end_pair;++ the_pair) {
01486         const HElement<DIM>& h_element0 = the_pair(0);
01487         const HElement<DIM>& h_element1 = the_pair(1);
01488         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01489         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01490         Assert(element0.index() == h_element0.index, ExcInternalError());
01491         Assert(element1.index() == h_element1.index, ExcInternalError());
01492         const std::vector<int>& element_dof1 = element1.dof();
01493         unsigned int n_element_dof1 = element_dof1.size();
01494         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01495           double volume = element1.templateElement().volume();
01496           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01497           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01498           int n_quadrature_point = quad_info.n_quadraturePoint();
01499           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01500           std::vector<value_type> f00_value = f00.value(q_point, element1);
01501           std::vector<value_type> f01_value = f01.value(q_point, element0);
01502           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01503           for (int l = 0;l < n_quadrature_point;l ++) {
01504             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01505             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01506             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01507               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01508               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01509             }
01510           }
01511         }
01512         else {
01513           double volume = element0.templateElement().volume();
01514           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01515           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01516           int n_quadrature_point = quad_info.n_quadraturePoint();
01517           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01518           std::vector<value_type> f00_value = f00.value(q_point, element1);
01519           std::vector<value_type> f01_value = f01.value(q_point, element0);
01520           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01521           for (int l = 0;l < n_quadrature_point;l ++) {
01522             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01523             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01524             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01525               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01526               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01527             }
01528           }
01529         }
01530       }
01531       for (unsigned int i = 0;i < n_dof;i ++)
01532         f1(i) /= mass_accumulation(i);
01533       break;
01534     }
01535     case LEAST_SQUARE: {
01536       const FEMSpace<value_type,DIM>& fem_space0 = f01.femSpace();
01537       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01538       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01539       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01540       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01541       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01542       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01543         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01544                   << std::endl;
01545         Assert(false, ExcInternalError());
01546       }
01547       f1.Vector<value_type>::operator=(0.0);
01548       MassMatrix<DIM, value_type> mass_matrix(fem_space1);
01549       mass_matrix.algebricAccuracy() = algebric_accuracy;
01550       mass_matrix.build();
01551       unsigned int n_dof = fem_space1.n_dof();
01552       Vector<double> rhs(n_dof);
01553       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01554       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01555       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01556       for (;the_pair != end_pair;++ the_pair) {
01557         const HElement<DIM>& h_element0 = the_pair(0);
01558         const HElement<DIM>& h_element1 = the_pair(1);
01559         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01560         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01561         Assert(element0.index() == h_element0.index, ExcInternalError());
01562         Assert(element1.index() == h_element1.index, ExcInternalError());
01563         const std::vector<int>& element_dof1 = element1.dof();
01564         unsigned int n_element_dof1 = element_dof1.size();
01565         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01566           double volume = element1.templateElement().volume();
01567           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01568           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01569           int n_quadrature_point = quad_info.n_quadraturePoint();
01570           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01571           std::vector<value_type> f00_value = f00.value(q_point, element1);
01572           std::vector<value_type> f01_value = f01.value(q_point, element0);
01573           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01574           for (int l = 0;l < n_quadrature_point;l ++) {
01575             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01576             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01577             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01578               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01579             }
01580           }
01581         }
01582         else {
01583           double volume = element0.templateElement().volume();
01584           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01585           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01586           int n_quadrature_point = quad_info.n_quadraturePoint();
01587           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01588           std::vector<value_type> f00_value = f00.value(q_point, element1);
01589           std::vector<value_type> f01_value = f01.value(q_point, element0);
01590           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01591           for (int l = 0;l < n_quadrature_point;l ++) {
01592             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01593             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01594             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01595               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01596             }
01597           }
01598         }
01599       }
01600       AMGSolver solver(mass_matrix);
01601       solver.solve(f1, rhs);
01602       break;
01603     }
01604     case LOCAL_LEAST_SQUARE: {
01605       const FEMSpace<value_type,DIM>& fem_space0 = f01.femSpace();
01606       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01607       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01608       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01609       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01610       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01611       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01612         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01613                   << std::endl;
01614         Assert(false, ExcInternalError());
01615       }
01616       unsigned int n_dof = fem_space1.n_dof();
01617       std::vector<int> counter(n_dof, 0);
01618       f1.Vector<value_type>::operator=(0.0);
01619       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01620       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01621       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01622       int the_element_index;
01623       unsigned int n_the_element_dof;
01624       { // only used to make the variables local
01625         const HElement<DIM>& the_h_element = the_pair(1);
01626         the_element_index = the_h_element.index;
01627         const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01628         Assert(the_element.index() == the_element_index, ExcInternalError());
01629         n_the_element_dof = the_element.n_dof();
01630       }
01631       Vector<double> local_rhs(n_the_element_dof);
01632       for (;the_pair != end_pair;) {
01633         const HElement<DIM>& h_element0 = the_pair(0);
01634         const HElement<DIM>& h_element1 = the_pair(1);
01635         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01636         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01637         Assert(element0.index() == h_element0.index, ExcInternalError());
01638         Assert(element1.index() == h_element1.index, ExcInternalError());
01639         const std::vector<int>& element_dof1 = element1.dof();
01640         unsigned int n_element_dof1 = element_dof1.size();
01641         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01642           double volume = element1.templateElement().volume();
01643           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01644           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01645           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01646           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01647           std::vector<value_type> f00_value = f00.value(q_point, element1);
01648           std::vector<value_type> f01_value = f01.value(q_point, element0);
01649           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01650           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01651             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01652             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01653             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01654               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01655             }
01656           }
01657         }
01658         else {
01659           double volume = element0.templateElement().volume();
01660           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01661           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01662           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01663           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01664           std::vector<value_type> f00_value = f00.value(q_point, element1);
01665           std::vector<value_type> f01_value = f01.value(q_point, element0);
01666           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01667           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01668             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01669             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01670             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01671               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01672             }
01673           }
01674         }
01675         ++ the_pair;
01676         if (the_pair == end_pair || the_pair(1).index != the_element_index) {
01677           const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01678           const std::vector<int>& the_element_dof = the_element.dof();
01679           FullMatrix<double> local_mass_matrix(n_the_element_dof, n_the_element_dof);
01680           Vector<value_type> local_f1(n_the_element_dof);
01681           double volume = element1.templateElement().volume();
01682           const QuadratureInfo<DIM>& quad_info = the_element.findQuadratureInfo(algebric_accuracy);
01683           std::vector<double> jacobian = the_element.local_to_global_jacobian(quad_info.quadraturePoint());
01684           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01685           std::vector<afepack::Point<DIM> > q_point = the_element.local_to_global(quad_info.quadraturePoint());
01686           std::vector<std::vector<value_type> > basis_value = the_element.basis_function_value(q_point);
01687           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01688             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01689             for (unsigned int j = 0;j < n_the_element_dof;j ++) {
01690               for (unsigned int k = 0;k < n_the_element_dof;k ++) {
01691                 local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
01692               }
01693             }
01694           }
01695           local_mass_matrix.gauss_jordan();
01696           local_mass_matrix.vmult(local_f1, local_rhs);
01697           for (unsigned int i = 0;i < n_the_element_dof;i ++) {
01698             f1(the_element_dof[i]) += local_f1(i);
01699             counter[the_element_dof[i]] ++;
01700           }
01701           if (the_pair != end_pair) {
01702             the_element_index = the_pair(1).index;
01703             const Element<value_type,DIM>& the_new_element = fem_space1.element(the_element_index);
01704             n_the_element_dof = the_new_element.n_dof();
01705             local_rhs.reinit(n_the_element_dof);
01706           }
01707         }
01708       }
01709       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
01710       break;
01711     }
01712     default: {
01713       Assert(false, ExcNotImplemented());
01714     }
01715     }
01716   }
01717   else if (&f01.femSpace() == &f1.femSpace()) {
01718     switch (method) {
01719     case MASS_ACCUMULATION: {
01720       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01721       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01722       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01723       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01724       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01725       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01726       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01727         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01728                   << std::endl;
01729         Assert(false, ExcInternalError());
01730       }
01731       unsigned int n_dof = fem_space1.n_dof();
01732       Vector<double> mass_accumulation(n_dof);
01733       f1.Vector<value_type>::operator=(0.0);
01734       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01735       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01736       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01737       for (;the_pair != end_pair;++ the_pair) {
01738         const HElement<DIM>& h_element0 = the_pair(0);
01739         const HElement<DIM>& h_element1 = the_pair(1);
01740         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01741         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01742         Assert(element0.index() == h_element0.index, ExcInternalError());
01743         Assert(element1.index() == h_element1.index, ExcInternalError());
01744         const std::vector<int>& element_dof1 = element1.dof();
01745         unsigned int n_element_dof1 = element_dof1.size();
01746         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01747           double volume = element1.templateElement().volume();
01748           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01749           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01750           int n_quadrature_point = quad_info.n_quadraturePoint();
01751           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01752           std::vector<value_type> f00_value = f00.value(q_point, element0);
01753           std::vector<value_type> f01_value = f01.value(q_point, element1);
01754           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01755           for (int l = 0;l < n_quadrature_point;l ++) {
01756             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01757             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01758             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01759               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01760               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01761             }
01762           }
01763         }
01764         else {
01765           double volume = element0.templateElement().volume();
01766           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01767           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01768           int n_quadrature_point = quad_info.n_quadraturePoint();
01769           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01770           std::vector<value_type> f00_value = f00.value(q_point, element0);
01771           std::vector<value_type> f01_value = f01.value(q_point, element1);
01772           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01773           for (int l = 0;l < n_quadrature_point;l ++) {
01774             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01775             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01776             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01777               f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01778               mass_accumulation(element_dof1[j]) += Jxw*basis_value1[j][l];
01779             }
01780           }
01781         }
01782       }
01783       for (unsigned int i = 0;i < n_dof;i ++)
01784         f1(i) /= mass_accumulation(i);
01785       break;
01786     }
01787     case LEAST_SQUARE: {
01788       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01789       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01790       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01791       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01792       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01793       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01794       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01795         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01796                   << std::endl;
01797         Assert(false, ExcInternalError());
01798       }
01799       f1.Vector<value_type>::operator=(0.0);
01800       MassMatrix<DIM, value_type> mass_matrix(fem_space1);
01801       mass_matrix.algebricAccuracy() = algebric_accuracy;
01802       mass_matrix.build();
01803       unsigned int n_dof = fem_space1.n_dof();
01804       Vector<double> rhs(n_dof);
01805       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01806       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01807       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01808       for (;the_pair != end_pair;++ the_pair) {
01809         const HElement<DIM>& h_element0 = the_pair(0);
01810         const HElement<DIM>& h_element1 = the_pair(1);
01811         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01812         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01813         Assert(element0.index() == h_element0.index, ExcInternalError());
01814         Assert(element1.index() == h_element1.index, ExcInternalError());
01815         const std::vector<int>& element_dof1 = element1.dof();
01816         unsigned int n_element_dof1 = element_dof1.size();
01817         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01818           double volume = element1.templateElement().volume();
01819           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01820           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01821           int n_quadrature_point = quad_info.n_quadraturePoint();
01822           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01823           std::vector<value_type> f00_value = f00.value(q_point, element0);
01824           std::vector<value_type> f01_value = f01.value(q_point, element1);
01825           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01826           for (int l = 0;l < n_quadrature_point;l ++) {
01827             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01828             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01829             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01830               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01831             }
01832           }
01833         }
01834         else {
01835           double volume = element0.templateElement().volume();
01836           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01837           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01838           int n_quadrature_point = quad_info.n_quadraturePoint();
01839           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01840           std::vector<value_type> f00_value = f00.value(q_point, element0);
01841           std::vector<value_type> f01_value = f01.value(q_point, element1);
01842           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01843           for (int l = 0;l < n_quadrature_point;l ++) {
01844             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01845             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01846             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01847               rhs(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
01848             }
01849           }
01850         }
01851       }
01852       AMGSolver solver(mass_matrix);
01853       solver.solve(f1, rhs);
01854       break;
01855     }
01856     case LOCAL_LEAST_SQUARE: {
01857       const FEMSpace<value_type,DIM>& fem_space0 = f00.femSpace();
01858       FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
01859       const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
01860       RegularMesh<DIM>& regular_mesh1 = static_cast<RegularMesh<DIM> &>(fem_space1.mesh());
01861       IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
01862       IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
01863       if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
01864         std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
01865                   << std::endl;
01866         Assert(false, ExcInternalError());
01867       }
01868       unsigned int n_dof = fem_space1.n_dof();
01869       std::vector<int> counter(n_dof, 0);
01870       f1.Vector<value_type>::operator=(0.0);
01871       IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
01872       ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
01873       ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
01874       int the_element_index;
01875       unsigned int n_the_element_dof;
01876       { // only used to make the variables local
01877         const HElement<DIM>& the_h_element = the_pair(1);
01878         the_element_index = the_h_element.index;
01879         const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01880         Assert(the_element.index() == the_element_index, ExcInternalError());
01881         n_the_element_dof = the_element.n_dof();
01882       }
01883       Vector<double> local_rhs(n_the_element_dof);
01884       for (;the_pair != end_pair;) {
01885         const HElement<DIM>& h_element0 = the_pair(0);
01886         const HElement<DIM>& h_element1 = the_pair(1);
01887         const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
01888         const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
01889         Assert(element0.index() == h_element0.index, ExcInternalError());
01890         Assert(element1.index() == h_element1.index, ExcInternalError());
01891         const std::vector<int>& element_dof1 = element1.dof();
01892         unsigned int n_element_dof1 = element_dof1.size();
01893         if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
01894           double volume = element1.templateElement().volume();
01895           const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
01896           std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
01897           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01898           std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
01899           std::vector<value_type> f00_value = f00.value(q_point, element0);
01900           std::vector<value_type> f01_value = f01.value(q_point, element1);
01901           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01902           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01903             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01904             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01905             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01906               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01907             }
01908           }
01909         }
01910         else {
01911           double volume = element0.templateElement().volume();
01912           const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
01913           std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
01914           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01915           std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
01916           std::vector<value_type> f00_value = f00.value(q_point, element0);
01917           std::vector<value_type> f01_value = f01.value(q_point, element1);
01918           std::vector<std::vector<value_type> > basis_value1 = element1.basis_function_value(q_point);
01919           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01920             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01921             value_type f_value = (*f)(f00_value[l], f01_value[l]);
01922             for (unsigned int j = 0;j < n_element_dof1;j ++) {
01923               local_rhs(j) += Jxw*f_value*basis_value1[j][l];
01924             }
01925           }
01926         }
01927         ++ the_pair;
01928         if (the_pair == end_pair || the_pair(1).index != the_element_index) {
01929           const Element<value_type,DIM>& the_element = fem_space1.element(the_element_index);
01930           const std::vector<int>& the_element_dof = the_element.dof();
01931           FullMatrix<double> local_mass_matrix(n_the_element_dof, n_the_element_dof);
01932           Vector<value_type> local_f1(n_the_element_dof);
01933           double volume = element1.templateElement().volume();
01934           const QuadratureInfo<DIM>& quad_info = the_element.findQuadratureInfo(algebric_accuracy);
01935           std::vector<double> jacobian = the_element.local_to_global_jacobian(quad_info.quadraturePoint());
01936           unsigned int n_quadrature_point = quad_info.n_quadraturePoint();
01937           std::vector<afepack::Point<DIM> > q_point = the_element.local_to_global(quad_info.quadraturePoint());
01938           std::vector<std::vector<value_type> > basis_value = the_element.basis_function_value(q_point);
01939           for (unsigned int l = 0;l < n_quadrature_point;l ++) {
01940             double Jxw = quad_info.weight(l)*jacobian[l]*volume;
01941             for (unsigned int j = 0;j < n_the_element_dof;j ++) {
01942               for (unsigned int k = 0;k < n_the_element_dof;k ++) {
01943                 local_mass_matrix(j, k) += Jxw*basis_value[j][l]*basis_value[k][l];
01944               }
01945             }
01946           }
01947           local_mass_matrix.gauss_jordan();
01948           local_mass_matrix.vmult(local_f1, local_rhs);
01949           for (unsigned int i = 0;i < n_the_element_dof;i ++) {
01950             f1(the_element_dof[i]) += local_f1(i);
01951             counter[the_element_dof[i]] ++;
01952           }
01953           if (the_pair != end_pair) {
01954             the_element_index = the_pair(1).index;
01955             const Element<value_type,DIM>& the_new_element = fem_space1.element(the_element_index);
01956             n_the_element_dof = the_new_element.n_dof();
01957             local_rhs.reinit(n_the_element_dof);
01958           }
01959         }
01960       }
01961       for (unsigned int i = 0;i < n_dof;i ++) f1(i) /= counter[i];
01962       break;
01963     }
01964     default: {
01965       Assert(false, ExcNotImplemented());
01966     }
01967     }
01968   }
01969   else { // something must be wrong
01970     Assert(false, ExcInternalError());
01971   }
01972 }
01973 
01974 #endif
01975