AFEPack
Operator.interpolate.templates.h
浏览该文件的文档。
00001 
00011 #ifndef _Operator_interpolate_templates_h_
00012 #define _Operator_interpolate_templates_h_
00013 
00014 #include "Operator.h"
00015 
00016 template <class value_type, int DIM> 
00017 void Operator::L2Interpolate(const FEMFunction<value_type, DIM>& f0, FEMFunction<value_type, DIM>& f1)
00018 {
00019   const FEMSpace<value_type,DIM>& fem_space0 = f0.femSpace();
00020   FEMSpace<value_type,DIM>& fem_space1 = f1.femSpace();
00021   const RegularMesh<DIM>& regular_mesh0 = static_cast<const RegularMesh<DIM> &>(fem_space0.mesh());
00022   const RegularMesh<DIM>& regular_mesh1 = static_cast<const RegularMesh<DIM> &>(fem_space1.mesh());
00023   IrregularMesh<DIM>& irregular_mesh0 = regular_mesh0.irregularMesh();
00024   IrregularMesh<DIM>& irregular_mesh1 = regular_mesh1.irregularMesh();
00025   if (&(irregular_mesh0.geometryTree()) != &(irregular_mesh1.geometryTree())) {
00026     std::cerr << "The two FEM functions are even not on the same hierarchy geometry tree."
00027               << std::endl;
00028     Assert(false, ExcInternalError());
00029   }
00030   std::vector<bool> flag(f1.size(), false);
00031   f1.Vector<value_type>::operator=(0.0);
00032   IrregularMeshPair<DIM> mesh_pair(irregular_mesh0, irregular_mesh1);
00033   ActiveElementPairIterator<DIM> the_pair = mesh_pair.beginActiveElementPair();
00034   ActiveElementPairIterator<DIM> end_pair = mesh_pair.endActiveElementPair();
00035   for (;the_pair != end_pair;++ the_pair) {
00036     const HElement<DIM>& h_element0 = the_pair(0);
00037     const HElement<DIM>& h_element1 = the_pair(1);
00038     const Element<value_type,DIM>& element0 = fem_space0.element(h_element0.index);
00039     const Element<value_type,DIM>& element1 = fem_space1.element(h_element1.index);
00040     Assert(element0.index() == h_element0.index, ExcInternalError());
00041     Assert(element1.index() == h_element1.index, ExcInternalError());
00042     const std::vector<int>& element_dof1 = element1.dof();
00043     unsigned int n_element_dof1 = element_dof1.size();
00044     if (the_pair.state() != ActiveElementPairIterator<DIM>::LESS_THAN) {
00045       for (unsigned int i = 0;i < n_element_dof1;i ++) {
00046         unsigned int j = element_dof1[i];
00047         const afepack::Point<DIM>& interp_point = fem_space1.dofInfo(j).interp_point;
00048         f1(j) = f0.value(interp_point, element0);
00049         flag[j] = true;
00050       }
00051     }
00052     else {
00053       for (unsigned int i = 0;i < n_element_dof1;i ++) {
00054         unsigned int j = element_dof1[i];
00055         if (flag[j]) continue;
00056         const afepack::Point<DIM>& interp_point = fem_space1.dofInfo(j).interp_point;
00057         if (h_element0.isIncludePoint(interp_point)) {
00058           f1(j) = f0.value(interp_point, element0);
00059           flag[j] = true;
00060         }
00061       }
00062     }
00063   }
00064 }
00065 
00066 
00067 
00068 template <class value_type, int DIM> 
00069 void Operator::L2Interpolate(value_type (*f0)(const double *), FEMFunction<value_type, DIM>& f1)
00070 {
00071   FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00072   typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00073   typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00074   for (;the_element != end_element;the_element ++) {
00075     const std::vector<int>& element_dof = the_element->dof();
00076     unsigned int n_element_dof = element_dof.size();
00077     for (unsigned int i = 0;i < n_element_dof;i ++) {
00078       unsigned int j = element_dof[i];
00079       const afepack::Point<DIM>& interp_point = fem_space.dofInfo(j).interp_point;
00080       f1(j) = f0(interp_point);
00081     }
00082   }
00083 }
00084 
00085 
00086 
00087 template <class value_type, int DIM> 
00088 void Operator::L2Interpolate(value_type (*f0)(const afepack::Point<DIM>&), FEMFunction<value_type, DIM>& f1)
00089 {
00090   FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00091   typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00092   typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00093   for (;the_element != end_element;the_element ++) {
00094     const std::vector<int>& element_dof = the_element->dof();
00095     unsigned int n_element_dof = element_dof.size();
00096     for (unsigned int i = 0;i < n_element_dof;i ++) {
00097       unsigned int j = element_dof[i];
00098       const afepack::Point<DIM>& interp_point = fem_space.dofInfo(j).interp_point;
00099       f1(j) = f0(interp_point);
00100     }
00101   }
00102 }
00103 
00104 
00105 
00106 template <class value_type, int DIM> 
00107 void Operator::L2Interpolate(const Function<value_type>& f0, FEMFunction<value_type, DIM>& f1)
00108 {
00109   FEMSpace<value_type,DIM>& fem_space = f1.femSpace();
00110   typename FEMSpace<value_type,DIM>::ElementIterator the_element = fem_space.beginElement();
00111   typename FEMSpace<value_type,DIM>::ElementIterator end_element = fem_space.endElement();
00112   for (;the_element != end_element;the_element ++) {
00113     const std::vector<int>& element_dof = the_element->dof();
00114     unsigned int n_element_dof = element_dof.size();
00115     for (unsigned int i = 0;i < n_element_dof;i ++) {
00116       unsigned int j = element_dof[i];
00117       const afepack::Point<DIM>& interp_point = fem_space.dofInfo(j).interp_point;
00118       f1(j) = f0.value(interp_point);
00119     }
00120   }
00121 }
00122 
00123 #endif
00124