AFEPack
Operator.discretize.templates.h
浏览该文件的文档。
00001 
00011 #ifndef _Operator_discretize_templates_h_
00012 #define _Operator_discretize_templates_h_
00013 
00014 #include "Operator.h"
00015 
00016 template <class value_type, int DIM>
00017   void Operator::L2Discretize(const FEMFunction<value_type, DIM>& f0, 
00018                               Vector<double>& f1, 
00019                               int algebric_accuracy)
00020 {
00021   const FEMSpace<value_type,DIM>& fem_space = f0.femSpace();
00022   typename FEMSpace<value_type,DIM>::ConstElementIterator 
00023     the_element = fem_space.beginElement(),
00024     end_element = fem_space.endElement();
00025   unsigned int n_dof = fem_space.n_dof();
00026   if (f1.size() != n_dof)
00027     f1.reinit(n_dof);
00028   else
00029     f1.Vector<value_type>::operator=(0.0);
00030   for (;the_element != end_element;++ the_element) {
00031     double volume = the_element->templateElement().volume();
00032     const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00033     std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00034     int n_quadrature_point = quad_info.n_quadraturePoint();
00035     std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00036     std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00037     std::vector<double> f0_value = f0.value(q_point, *the_element);
00038     const std::vector<int>& element_dof = the_element->dof();
00039     unsigned int n_element_dof = element_dof.size();
00040     for (int l = 0;l < n_quadrature_point;l ++) {
00041       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00042       for (unsigned int j = 0;j < n_element_dof;j ++) {
00043         f1(element_dof[j]) += Jxw*f0_value[l]*basis_value[j][l];
00044       }
00045     }
00046   }
00047 }
00048 
00049 
00050 
00051 template <class value_type, int DIM>
00052 void Operator::L2Discretize(const FEMFunction<value_type, DIM>& f0, 
00053                             const FEMSpace<value_type, DIM>& fem_space1, 
00054                             Vector<double>& f1, 
00055                             int algebric_accuracy)
00056 {
00057   const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00058   const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00059   const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00060   IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00061   IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00062   if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00063     std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00064               << std::endl;
00065     Assert(false, ExcInternalError());
00066   }
00067   unsigned int n_dof = fem_space1.n_dof();
00068   if (f1.size() != n_dof)
00069     f1.reinit(n_dof);
00070   else
00071     f1.Vector<value_type>::operator=(0.0);
00072   IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00073   ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00074   ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00075   for (;the_pair != end_pair;++ the_pair) {
00076     const HElement<DIM>& h_element0 = the_pair(0);
00077     const HElement<DIM>& h_element1 = the_pair(1);
00078     const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00079     const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00080     Assert(element0.index() == h_element0.index, ExcInternalError());
00081     Assert(element1.index() == h_element1.index, ExcInternalError());
00082     const std::vector<int>& element_dof1 = element1.dof();
00083     unsigned int n_element_dof1 = element_dof1.size();
00084     if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00085       double volume = element1.templateElement().volume();
00086       const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00087       std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00088       int n_quadrature_point = quad_info.n_quadraturePoint();
00089       std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00090       std::vector<double> f0_value = f0.value(q_point, element0);
00091       std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00092       for (int l = 0;l < n_quadrature_point;l ++) {
00093         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00094         for (unsigned int j = 0;j < n_element_dof1;j ++) {
00095           f1(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00096         }
00097       }
00098     }
00099     else {
00100       double volume = element0.templateElement().volume();
00101       const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00102       std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00103       int n_quadrature_point = quad_info.n_quadraturePoint();
00104       std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00105       std::vector<double> f0_value = f0.value(q_point, element0);
00106       std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00107       for (int l = 0;l < n_quadrature_point;l ++) {
00108         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00109         for (unsigned int j = 0;j < n_element_dof1;j ++) {
00110           f1(element_dof1[j]) += Jxw*f0_value[l]*basis_value1[j][l];
00111         }
00112       }
00113     }
00114   }
00115 }
00116 
00117 
00118 
00119 template <class value_type, int DIM> 
00120 void Operator::L2Discretize(value_type (*f0)(const double *), 
00121                             const FEMSpace<value_type, DIM>& fem_space, 
00122                             Vector<double>& f1, 
00123                             int algebric_accuracy)
00124 {
00125   unsigned int n_dof = fem_space.n_dof();
00126   if (f1.size() != n_dof)
00127     f1.reinit(n_dof);
00128   else
00129     f1.Vector<value_type>::operator=(0.0);
00130   typename FEMSpace<value_type,DIM>::ConstElementIterator 
00131     the_element = fem_space.beginElement(),
00132     end_element = fem_space.endElement();
00133   for (;the_element != end_element;the_element ++) {
00134     double volume = the_element->templateElement().volume();
00135     const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00136     std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00137     int n_quadrature_point = quad_info.n_quadraturePoint();
00138     std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00139     std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00140     const std::vector<int>& element_dof = the_element->dof();
00141     unsigned int n_element_dof = element_dof.size();
00142     for (int l = 0;l < n_quadrature_point;l ++) {
00143       value_type f0_value = (*f0)(q_point[l]);
00144       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00145       for (unsigned int j = 0;j < n_element_dof;j ++) {
00146         f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00147       }
00148     }
00149   }
00150 }
00151 
00152 
00153 
00154 template <class value_type, int DIM> 
00155 void Operator::L2Discretize(value_type (*f0)(const afepack::Point<DIM>&), 
00156                             const FEMSpace<value_type, DIM>& fem_space, 
00157                             Vector<double>& f1, 
00158                             int algebric_accuracy)
00159 {
00160   unsigned int n_dof = fem_space.n_dof();
00161   if (f1.size() != n_dof)
00162     f1.reinit(n_dof);
00163   else
00164     f1.Vector<value_type>::operator=(0.0);
00165   typename FEMSpace<value_type,DIM>::ConstElementIterator 
00166     the_element = fem_space.beginElement(),
00167     end_element = fem_space.endElement();
00168   for (;the_element != end_element;the_element ++) {
00169     double volume = the_element->templateElement().volume();
00170     const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00171     std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00172     int n_quadrature_point = quad_info.n_quadraturePoint();
00173     std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00174     std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00175     const std::vector<int>& element_dof = the_element->dof();
00176     unsigned int n_element_dof = element_dof.size();
00177     for (int l = 0;l < n_quadrature_point;l ++) {
00178       value_type f0_value = (*f0)(q_point[l]);
00179       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00180       for (unsigned int j = 0;j < n_element_dof;j ++) {
00181         f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00182       }
00183     }
00184   }
00185 }
00186 
00187 
00188 
00189 
00190 template <class value_type, int DIM> 
00191 void Operator::L2Discretize(const Function<value_type>& f0, 
00192                             const FEMSpace<value_type, DIM>& fem_space, 
00193                             Vector<double>& f1, 
00194                             int algebric_accuracy)
00195 {
00196   unsigned int n_dof = fem_space.n_dof();
00197   if (f1.size() != n_dof)
00198     f1.reinit(n_dof);
00199   else
00200     f1.Vector<value_type>::operator=(0.0);
00201   typename FEMSpace<value_type,DIM>::ConstElementIterator 
00202     the_element = fem_space.beginElement(),
00203     end_element = fem_space.endElement();
00204   for (;the_element != end_element;the_element ++) {
00205     double volume = the_element->templateElement().volume();
00206     const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00207     std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00208     int n_quadrature_point = quad_info.n_quadraturePoint();
00209     std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00210     std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00211     const std::vector<int>& element_dof = the_element->dof();
00212     unsigned int n_element_dof = element_dof.size();
00213     for (int l = 0;l < n_quadrature_point;l ++) {
00214       value_type f0_value = f0.value(q_point[l]);
00215       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00216       for (unsigned int j = 0;j < n_element_dof;j ++) {
00217         f1(element_dof[j]) += Jxw*f0_value*basis_value[j][l];
00218       }
00219     }
00220   }
00221 }
00222 
00223 
00224 
00225 
00226 template <class value_type, int DIM>
00227 void Operator::L2Discretize(value_type (*f)(const value_type&), 
00228                             const FEMFunction<value_type, DIM>& f0,
00229                             Vector<double>& f1, 
00230                             int algebric_accuracy)
00231 {
00232   const FEMSpace<value_type,DIM>& fem_space = f0.femSpace();
00233   typename FEMSpace<value_type,DIM>::ConstElementIterator 
00234     the_element = fem_space.beginElement(),
00235     end_element = fem_space.endElement();
00236   unsigned int n_dof = fem_space.n_dof();
00237   if (f1.size() != n_dof)
00238     f1.reinit(n_dof);
00239   else
00240     f1.Vector<value_type>::operator=(0.0);
00241   for (;the_element != end_element;++ the_element) {
00242     double volume = the_element->templateElement().volume();
00243     const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00244     std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00245     int n_quadrature_point = quad_info.n_quadraturePoint();
00246     std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00247     std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00248     std::vector<double> f0_value = f0.value(q_point, *the_element);
00249     const std::vector<int>& element_dof = the_element->dof();
00250     unsigned int n_element_dof = element_dof.size();
00251     for (int l = 0;l < n_quadrature_point;l ++) {
00252       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00253       value_type f_value = (*f)(f0_value[l]);
00254       for (unsigned int j = 0;j < n_element_dof;j ++) {
00255         f1(element_dof[j]) += Jxw*f_value*basis_value[j][l];
00256       }
00257     }
00258   }
00259 }
00260 
00261 
00262 
00263 template <class value_type, int DIM>
00264 void Operator::L2Discretize(value_type (*f)(const value_type&), 
00265                             const FEMFunction<value_type, DIM>& f0,
00266                             const FEMSpace<value_type, DIM>& fem_space1, 
00267                             Vector<double>& f1, 
00268                             int algebric_accuracy)
00269 {
00270   const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00271   unsigned int n_dof = fem_space1.n_dof();
00272   if (f1.size() != n_dof)
00273     f1.reinit(n_dof);
00274   else
00275     f1.Vector<value_type>::operator=(0.0);
00276   if (&fem_space0 == &fem_space1 || &(fem_space0.mesh()) == &(fem_space1.mesh())) {
00277     typename FEMSpace<value_type,DIM>::ConstElementIterator 
00278       the_element = fem_space0.beginElement(),
00279       end_element = fem_space0.endElement();
00280     typename FEMSpace<value_type,DIM>::ConstElementIterator 
00281       the_element1 = fem_space1.beginElement();
00282     for (;the_element != end_element;++ the_element, ++ the_element1) {
00283       double volume = the_element->templateElement().volume();
00284       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00285       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00286       int n_quadrature_point = quad_info.n_quadraturePoint();
00287       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00288       std::vector<std::vector<double> > basis_value = the_element1->basis_function_value(q_point);
00289       std::vector<double> f0_value = f0.value(q_point, *the_element);
00290       const std::vector<int>& element_dof = the_element1->dof();
00291       unsigned int n_element_dof = element_dof.size();
00292       for (int l = 0;l < n_quadrature_point;l ++) {
00293         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00294         value_type f_value = (*f)(f0_value[l]);
00295         for (unsigned int j = 0;j < n_element_dof;j ++) {
00296           f1(element_dof[j]) += Jxw*f_value*basis_value[j][l];
00297         }
00298       }
00299     }
00300     return;
00301   }
00302   const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00303   const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00304   IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00305   IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00306   if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00307     std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00308               << std::endl;
00309     Assert(false, ExcInternalError());
00310   }
00311   IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00312   ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00313   ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00314   for (;the_pair != end_pair;++ the_pair) {
00315     const HElement<DIM>& h_element0 = the_pair(0);
00316     const HElement<DIM>& h_element1 = the_pair(1);
00317     const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00318     const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00319     Assert(element0.index() == h_element0.index, ExcInternalError());
00320     Assert(element1.index() == h_element1.index, ExcInternalError());
00321     const std::vector<int>& element_dof1 = element1.dof();
00322     unsigned int n_element_dof1 = element_dof1.size();
00323     if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00324       double volume = element1.templateElement().volume();
00325       const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00326       std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00327       int n_quadrature_point = quad_info.n_quadraturePoint();
00328       std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00329       std::vector<double> f0_value = f0.value(q_point, element0);
00330       std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00331       for (int l = 0;l < n_quadrature_point;l ++) {
00332         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00333         value_type f_value = (*f)(f0_value[l]);
00334         for (unsigned int j = 0;j < n_element_dof1;j ++) {
00335           f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00336         }
00337       }
00338     }
00339     else {
00340       double volume = element0.templateElement().volume();
00341       const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00342       std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00343       int n_quadrature_point = quad_info.n_quadraturePoint();
00344       std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00345       std::vector<double> f0_value = f0.value(q_point, element0);
00346       std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00347       for (int l = 0;l < n_quadrature_point;l ++) {
00348         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00349         value_type f_value = (*f)(f0_value[l]);
00350         for (unsigned int j = 0;j < n_element_dof1;j ++) {
00351           f1(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00352         }
00353       }
00354     }
00355   }
00356 }
00357 
00358 
00359 
00360 template <class value_type, int DIM>
00361 void Operator::L2Discretize(value_type (*f)(const value_type&, const value_type&),
00362                             const FEMFunction<value_type, DIM>& f0, 
00363                             const FEMFunction<value_type, DIM>& f1,
00364                             const FEMSpace<value_type, DIM>& fem_space2, 
00365                             Vector<double>& f2, 
00366                             int algebric_accuracy)
00367 {
00368   const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00369   const FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00370   if ((&fem_space0 != &fem_space1) &&
00371       (&fem_space0 != &fem_space2) &&
00372       (&fem_space1 != &fem_space2)) {
00373     std::cerr << "The three FEM functions are on three different finite element spaces."
00374               << std::endl;
00375     abort();
00376   }
00377   unsigned int n_dof = fem_space2.n_dof();
00378   if (f2.size() != n_dof)
00379     f2.reinit(n_dof);
00380   else
00381     f2.Vector<value_type>::operator=(0.0);
00382   if ((&fem_space0 == &fem_space1) && (&fem_space0 == &fem_space2)) {
00383     typename FEMSpace<value_type,DIM>::ConstElementIterator 
00384       the_element = fem_space0.beginElement(),
00385       end_element = fem_space0.endElement();
00386     for (;the_element != end_element;++ the_element) {
00387       double volume = the_element->templateElement().volume();
00388       const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(algebric_accuracy);
00389       std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
00390       int n_quadrature_point = quad_info.n_quadraturePoint();
00391       std::vector<afepack::Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
00392       std::vector<std::vector<double> > basis_value = the_element->basis_function_value(q_point);
00393       std::vector<double> f0_value = f0.value(q_point, *the_element);
00394       std::vector<double> f1_value = f1.value(q_point, *the_element);
00395       const std::vector<int>& element_dof = the_element->dof();
00396       unsigned int n_element_dof = element_dof.size();
00397       for (int l = 0;l < n_quadrature_point;l ++) {
00398         double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00399         value_type f_value = (*f)(f0_value[l], f1_value[l]);
00400         for (unsigned int j = 0;j < n_element_dof;j ++) {
00401           f2(element_dof[j]) += Jxw*f_value*basis_value[j][l];
00402         }
00403       }
00404     }
00405   }
00406   else if (&fem_space0 == &fem_space1) {
00407     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00408     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space2.mesh());
00409     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00410     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00411     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00412       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00413                 << std::endl;
00414       Assert(false, ExcInternalError());
00415     }
00416     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00417     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00418     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00419     for (;the_pair != end_pair;++ the_pair) {
00420       const HElement<DIM>& h_element0 = the_pair(0);
00421       const HElement<DIM>& h_element1 = the_pair(1);
00422       const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00423       const Element<value_type,DIM>& element1 = fem_space2.element(h_element1.index);
00424       Assert(element0.index() == h_element0.index, ExcInternalError());
00425       Assert(element1.index() == h_element1.index, ExcInternalError());
00426       const std::vector<int>& element_dof1 = element1.dof();
00427       unsigned int n_element_dof1 = element_dof1.size();
00428       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00429         double volume = element1.templateElement().volume();
00430         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00431         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00432         int n_quadrature_point = quad_info.n_quadraturePoint();
00433         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00434         std::vector<double> f0_value = f0.value(q_point, element0);
00435         std::vector<double> f1_value = f1.value(q_point, element0);
00436         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00437         for (int l = 0;l < n_quadrature_point;l ++) {
00438           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00439           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00440           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00441             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00442           }
00443         }
00444       }
00445       else {
00446         double volume = element0.templateElement().volume();
00447         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00448         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00449         int n_quadrature_point = quad_info.n_quadraturePoint();
00450         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00451         std::vector<double> f0_value = f0.value(q_point, element0);
00452         std::vector<double> f1_value = f1.value(q_point, element0);
00453         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00454         for (int l = 0;l < n_quadrature_point;l ++) {
00455           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00456           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00457           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00458             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00459           }
00460         }
00461       }
00462     }
00463   }
00464   else if (&fem_space0 == &fem_space2) {
00465     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00466     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space2.mesh());
00467     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00468     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00469     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00470       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00471                 << std::endl;
00472       Assert(false, ExcInternalError());
00473     }
00474     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00475     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00476     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00477     for (;the_pair != end_pair;++ the_pair) {
00478       const HElement<DIM>& h_element0 = the_pair(0);
00479       const HElement<DIM>& h_element1 = the_pair(1);
00480       const Element<value_type,DIM>& element0 = fem_space1.element(h_element0.index);
00481       const Element<value_type,DIM>& element1 = fem_space2.element(h_element1.index);
00482       Assert(element0.index() == h_element0.index, ExcInternalError());
00483       Assert(element1.index() == h_element1.index, ExcInternalError());
00484       const std::vector<int>& element_dof1 = element1.dof();
00485       unsigned int n_element_dof1 = element_dof1.size();
00486       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00487         double volume = element1.templateElement().volume();
00488         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00489         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00490         int n_quadrature_point = quad_info.n_quadraturePoint();
00491         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00492         std::vector<double> f0_value = f0.value(q_point, element1);
00493         std::vector<double> f1_value = f1.value(q_point, element0);
00494         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00495         for (int l = 0;l < n_quadrature_point;l ++) {
00496           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00497           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00498           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00499             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00500           }
00501         }
00502       }
00503       else {
00504         double volume = element0.templateElement().volume();
00505         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00506         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00507         int n_quadrature_point = quad_info.n_quadraturePoint();
00508         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00509         std::vector<double> f0_value = f0.value(q_point, element1);
00510         std::vector<double> f1_value = f1.value(q_point, element0);
00511         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00512         for (int l = 0;l < n_quadrature_point;l ++) {
00513           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00514           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00515           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00516             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00517           }
00518         }
00519       }
00520     }
00521   }
00522   else if (&fem_space1 == &fem_space2) {
00523     const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00524     const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space2.mesh());
00525     IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00526     IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00527     if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00528       std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00529                 << std::endl;
00530       Assert(false, ExcInternalError());
00531     }
00532     IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00533     ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00534     ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00535     for (;the_pair != end_pair;++ the_pair) {
00536       const HElement<DIM>& h_element0 = the_pair(0);
00537       const HElement<DIM>& h_element1 = the_pair(1);
00538       const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00539       const Element<value_type,DIM>& element1 = fem_space2.element(h_element1.index);
00540       Assert(element0.index() == h_element0.index, ExcInternalError());
00541       Assert(element1.index() == h_element1.index, ExcInternalError());
00542       const std::vector<int>& element_dof1 = element1.dof();
00543       unsigned int n_element_dof1 = element_dof1.size();
00544       if (the_pair.state() == ActiveElementPairIterator<DIM>::GREAT_THAN) {
00545         double volume = element1.templateElement().volume();
00546         const QuadratureInfo<DIM>& quad_info = element1.findQuadratureInfo(algebric_accuracy);
00547         std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00548         int n_quadrature_point = quad_info.n_quadraturePoint();
00549         std::vector<afepack::Point<DIM> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00550         std::vector<double> f0_value = f0.value(q_point, element0);
00551         std::vector<double> f1_value = f1.value(q_point, element1);
00552         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00553         for (int l = 0;l < n_quadrature_point;l ++) {
00554           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00555           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00556           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00557             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00558           }
00559         }
00560       }
00561       else {
00562         double volume = element0.templateElement().volume();
00563         const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(algebric_accuracy);
00564         std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00565         int n_quadrature_point = quad_info.n_quadraturePoint();
00566         std::vector<afepack::Point<DIM> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00567         std::vector<double> f0_value = f0.value(q_point, element0);
00568         std::vector<double> f1_value = f1.value(q_point, element1);
00569         std::vector<std::vector<double> > basis_value1 = element1.basis_function_value(q_point);
00570         for (int l = 0;l < n_quadrature_point;l ++) {
00571           double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00572           value_type f_value = (*f)(f0_value[l], f1_value[l]);
00573           for (unsigned int j = 0;j < n_element_dof1;j ++) {
00574             f2(element_dof1[j]) += Jxw*f_value*basis_value1[j][l];
00575           }
00576         }
00577       }
00578     }
00579   }
00580   else { // something must be wrong
00581     Assert(false, ExcInternalError());
00582   }
00583 }
00584 
00585 
00586 
00587 #endif
00588