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