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