AFEPack
|
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