AFEPack
BilinearOperator.templates.h
浏览该文件的文档。
00001 
00011 #ifndef _BilinearOperator_templates_h_
00012 #define _BilinearOperator_templates_h_
00013 
00014 #include "BilinearOperator.h"
00015 
00016 #define TEMPLATE template <int DIM, class value_type0, class value_type1, int DOW, int TDIM0, int TDIM1>
00017 #define THIS BilinearOperator<DIM,value_type0,value_type1,DOW,TDIM0,TDIM1>
00018 
00019 TEMPLATE
00020 THIS::BilinearOperator()
00021 {}
00022 
00023 
00024 
00025 TEMPLATE
00026 THIS::BilinearOperator(FEMSpace<value_type0,DIM,DOW,TDIM0>& sp0, 
00027                        FEMSpace<value_type1,DIM,DOW,TDIM1>& sp1)
00028 {
00029   fem_space0 = &sp0;
00030   fem_space1 = &sp1;
00031 }
00032 
00033 
00034 
00035 TEMPLATE
00036 THIS::~BilinearOperator()
00037 {
00038   SparseMatrix<double>::clear();
00039 }
00040 
00041 
00042 
00043 TEMPLATE
00044 void THIS::reinit(FEMSpace<value_type0,DIM,DOW,TDIM0>& sp0, 
00045                   FEMSpace<value_type1,DIM,DOW,TDIM1>& sp1)
00046 {
00047   fem_space0 = &sp0;
00048   fem_space1 = &sp1;
00049 }
00050 
00051 
00052 
00053 TEMPLATE
00054 void THIS::build()
00055 {
00056   buildSparsityPattern();
00057   buildSparseMatrix();
00058 }
00059 
00060 
00061 
00062 TEMPLATE
00063 void THIS::buildSparsityPattern()
00064 {
00065   buildDofInfo();
00066   sparsity_pattern.reinit(nDof0(), nDof1(), nMaxCouplingDof());
00067   if ((void *)fem_space0 == (void *)fem_space1) {
00068     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00069     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00070     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00071     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00072       getElementPattern(*the_element0, *the_element1);
00073       addElementPattern();
00074     }
00075   }
00076   else if (&(fem_space0->mesh()) == &(fem_space1->mesh())) {
00077     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00078     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00079     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00080     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00081       Assert (the_element0->index() == the_element1->index(), ExcInternalError());
00082       getElementPattern(*the_element0, *the_element1);
00083       addElementPattern();
00084     }
00085   }
00086   else {
00087     Mesh<DIM,DOW>& mesh0 = fem_space0->mesh();
00088     Mesh<DIM,DOW>& mesh1 = fem_space1->mesh();
00089     RegularMesh<DIM,DOW>& regular_mesh0 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh0);
00090     RegularMesh<DIM,DOW>& regular_mesh1 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh1);
00091     IrregularMesh<DIM,DOW>& irregular_mesh0 = regular_mesh0.irregularMesh();
00092     IrregularMesh<DIM,DOW>& irregular_mesh1 = regular_mesh1.irregularMesh();
00093     Assert (&(irregular_mesh0.geometryTree()) == &(irregular_mesh1.geometryTree()), ExcInternalError());
00094     IrregularMeshPair<DIM,DOW> mesh_pair(irregular_mesh0, irregular_mesh1);
00095     ActiveElementPairIterator<DIM,DOW> the_pair = mesh_pair.beginActiveElementPair();
00096     ActiveElementPairIterator<DIM,DOW> end_pair = mesh_pair.endActiveElementPair();
00097     for (;the_pair != end_pair;++ the_pair) {
00098       const HElement<DIM,DOW>& h_element0 = the_pair(0);
00099       const HElement<DIM,DOW>& h_element1 = the_pair(1);
00100       Element<value_type0,DIM,DOW,TDIM0>& element0 = fem_space0->element(h_element0.index);
00101       Element<value_type1,DIM,DOW,TDIM1>& element1 = fem_space1->element(h_element1.index);
00102       Assert (element0.index() == h_element0.index, ExcInternalError());
00103       Assert (element1.index() == h_element1.index, ExcInternalError());
00104       getElementPattern(element0, element1);
00105       addElementPattern();
00106     }
00107   }
00108   sparsity_pattern.compress();
00109 }
00110 
00111 
00112 
00113 TEMPLATE
00114 void THIS::buildDofInfo()
00115 {
00116   nDof0() = fem_space0->n_dof();
00117   nDof1() = fem_space1->n_dof();
00118   std::vector<int> n_coupling_dof(nDof0(), 0);
00119   if ((void *)fem_space0 == (void *)fem_space1) {
00120     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00121     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00122     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00123     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00124       getElementPattern(*the_element0, *the_element1);
00125       int n_element_dof = elementDof0().size();
00126       for (int i = 0;i < n_element_dof;i ++)
00127         n_coupling_dof[elementDof0(i)] += n_element_dof;
00128     }
00129   }
00130   else if (&(fem_space0->mesh()) == &(fem_space1->mesh())) {
00131     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00132     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00133     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00134     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00135       Assert (the_element0->index() == the_element1->index(), ExcInternalError());
00136       getElementPattern(*the_element0, *the_element1);
00137       int n_element_dof0 = elementDof0().size();
00138       int n_element_dof1 = elementDof1().size();
00139       for (int i = 0;i < n_element_dof0;i ++)
00140         n_coupling_dof[elementDof0(i)] += n_element_dof1;
00141     }
00142   }
00143   else {
00144     Mesh<DIM,DOW>& mesh0 = fem_space0->mesh();
00145     Mesh<DIM,DOW>& mesh1 = fem_space1->mesh();
00146     RegularMesh<DIM,DOW>& regular_mesh0 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh0);
00147     RegularMesh<DIM,DOW>& regular_mesh1 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh1);
00148     IrregularMesh<DIM,DOW>& irregular_mesh0 = regular_mesh0.irregularMesh();
00149     IrregularMesh<DIM,DOW>& irregular_mesh1 = regular_mesh1.irregularMesh();
00150     Assert (&(irregular_mesh0.geometryTree()) == &(irregular_mesh0.geometryTree()), ExcInternalError());
00151     IrregularMeshPair<DIM,DOW> mesh_pair(irregular_mesh0, irregular_mesh1);
00152     ActiveElementPairIterator<DIM,DOW> the_pair = mesh_pair.beginActiveElementPair();
00153     ActiveElementPairIterator<DIM,DOW> end_pair = mesh_pair.endActiveElementPair();
00154     for (;the_pair != end_pair;++ the_pair) {
00155       const HElement<DIM,DOW>& h_element0 = the_pair(0);
00156       const HElement<DIM,DOW>& h_element1 = the_pair(1);
00157       Element<value_type0,DIM,DOW,TDIM0>& element0 = fem_space0->element(h_element0.index);
00158       Element<value_type1,DIM,DOW,TDIM1>& element1 = fem_space1->element(h_element1.index);
00159       Assert (element0.index() == h_element0.index, ExcInternalError());
00160       Assert (element1.index() == h_element1.index, ExcInternalError());
00161       getElementPattern(element0, element1);
00162       int n_element_dof0 = elementDof0().size();
00163       int n_element_dof1 = elementDof1().size();
00164       for (int i = 0;i < n_element_dof0;i ++)
00165         n_coupling_dof[elementDof0(i)] += n_element_dof1;
00166     }
00167   }
00168   nMaxCouplingDof() = *max_element(n_coupling_dof.begin(), n_coupling_dof.end());
00169   if (nMaxCouplingDof() > nDof1()) nMaxCouplingDof() = nDof1();
00170 }
00171 
00172 
00173 
00174 TEMPLATE
00175 void THIS::buildSparseMatrix()
00176 {
00177   SparseMatrix<double>::reinit(sparsity_pattern);
00178   if ((void *)fem_space0 == (void *)fem_space1) {
00179     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00180     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00181     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00182     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00183       getElementPattern(*the_element0, *the_element1);
00184       elementMatrix().reinit(elementDof0().size(), elementDof1().size());
00185       getElementMatrix(*the_element0, *the_element1);
00186       addElementMatrix();
00187     }
00188   }
00189   else if (&(fem_space0->mesh()) == &(fem_space1->mesh())) {
00190     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator the_element0 = fem_space0->beginElement();
00191     typename std::vector<Element<value_type0,DIM,DOW,TDIM0> >::iterator end_element0 = fem_space0->endElement();
00192     typename std::vector<Element<value_type1,DIM,DOW,TDIM1> >::iterator the_element1 = fem_space1->beginElement();
00193     for (;the_element0 != end_element0;the_element0 ++, the_element1 ++) {
00194       Assert (the_element0->index() == the_element1->index(), ExcInternalError());
00195       getElementPattern(*the_element0, *the_element1);
00196       elementMatrix().reinit(elementDof0().size(), elementDof1().size());
00197       getElementMatrix(*the_element0, *the_element1);
00198       addElementMatrix();
00199     }
00200   }
00201   else {
00202     Mesh<DIM,DOW>& mesh0 = fem_space0->mesh();
00203     Mesh<DIM,DOW>& mesh1 = fem_space1->mesh();
00204     RegularMesh<DIM,DOW>& regular_mesh0 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh0);
00205     RegularMesh<DIM,DOW>& regular_mesh1 = dynamic_cast<RegularMesh<DIM,DOW>&>(mesh1);
00206     IrregularMesh<DIM,DOW>& irregular_mesh0 = regular_mesh0.irregularMesh();
00207     IrregularMesh<DIM,DOW>& irregular_mesh1 = regular_mesh1.irregularMesh();
00208     Assert (&(irregular_mesh0.geometryTree()) == &(irregular_mesh1.geometryTree()), ExcInternalError());
00209     IrregularMeshPair<DIM,DOW> mesh_pair(irregular_mesh0, irregular_mesh1);
00210     ActiveElementPairIterator<DIM,DOW> the_pair = mesh_pair.beginActiveElementPair();
00211     ActiveElementPairIterator<DIM,DOW> end_pair = mesh_pair.endActiveElementPair();
00212     for (;the_pair != end_pair;++ the_pair) {
00213       const HElement<DIM,DOW>& h_element0 = the_pair(0);
00214       const HElement<DIM,DOW>& h_element1 = the_pair(1);
00215       Element<value_type0,DIM,DOW,TDIM0>& element0 = fem_space0->element(h_element0.index);
00216       Element<value_type1,DIM,DOW,TDIM1>& element1 = fem_space1->element(h_element1.index);
00217       Assert (element0.index() == h_element0.index, ExcInternalError());
00218       Assert (element1.index() == h_element1.index, ExcInternalError());
00219       getElementPattern(element0, element1);
00220       elementMatrix().reinit(elementDof0().size(), elementDof1().size());
00221       getElementMatrix(element0, element1, the_pair.state());
00222       addElementMatrix();
00223     }
00224   }
00225 }
00226 
00227 
00228 
00229 TEMPLATE
00230 void THIS::addElementPattern()
00231 {
00232   int n_element_dof0 = elementDof0().size();
00233   int n_element_dof1 = elementDof1().size();
00234   for (int i = 0;i < n_element_dof0;i ++) {
00235     for (int j = 0;j < n_element_dof1;j ++) {
00236       sparsity_pattern.add(elementDof0(i), elementDof1(j));
00237     }
00238   }
00239 }
00240 
00241 
00242 
00243 TEMPLATE
00244 void THIS::addElementMatrix()
00245 {
00246   int n_element_dof0 = elementDof0().size();
00247   int n_element_dof1 = elementDof1().size();
00248   for (int i = 0;i < n_element_dof0;i ++) {
00249     for (int j = 0;j < n_element_dof1;j ++) {
00250       SparseMatrix<double>::add(elementDof0(i), elementDof1(j), element_matrix(i,j));
00251     }
00252   }
00253 }
00254 
00255 
00256 
00257 TEMPLATE
00258 void THIS::getElementPattern(
00259                              const Element<value_type0,DIM,DOW,TDIM0>& element0, 
00260                              const Element<value_type1,DIM,DOW,TDIM1>& element1)
00261 {
00262   element_dof0 = &(element0.dof());
00263   element_dof1 = &(element1.dof());
00264 }
00265 
00266 #undef THIS
00267 #undef TEMPLATE
00268 
00269 #define TEMPLATE template <int DIM, class value_type0, class value_type1, int DOW, int TDIM0, int TDIM1>
00270 #define THIS L2InnerProduct<DIM,value_type0,value_type1,DOW,TDIM0,TDIM1>
00271 
00272 TEMPLATE
00273 void THIS::reinit(FEMSpace<value_type0,DIM,DOW,TDIM0>& sp0, 
00274                   FEMSpace<value_type1,DIM,DOW,TDIM1>& sp1)
00275 {
00276   BilinearOperator<DIM,value_type0,value_type1,DOW,TDIM0,TDIM1>::reinit(sp0, sp1);      
00277 }
00278 
00279 TEMPLATE
00280 void THIS::getElementMatrix(const Element<value_type0,DIM,DOW,TDIM0>& element0,
00281                             const Element<value_type1,DIM,DOW,TDIM1>& element1,
00282                             const typename ActiveElementPairIterator<DIM>::State state)
00283 {
00284   const std::vector<int>& ele_dof0 = element0.dof();
00285   const std::vector<int>& ele_dof1 = element1.dof();
00286   int n_element_dof0 = ele_dof0.size();
00287   int n_element_dof1 = ele_dof1.size();
00288   if (state == ActiveElementPairIterator<DIM,DOW>::GREAT_THAN) {
00289     double volume = element1.templateElement().volume();
00290     const int& alg_acc = this->algebricAccuracy();
00291     const QuadratureInfo<TDIM1>& quad_info = element1.findQuadratureInfo(alg_acc);
00292     std::vector<double> jacobian = element1.local_to_global_jacobian(quad_info.quadraturePoint());
00293     int n_quadrature_point = quad_info.n_quadraturePoint();
00294     std::vector<afepack::Point<DOW> > q_point = element1.local_to_global(quad_info.quadraturePoint());
00295     std::vector<std::vector<value_type0> > basis_value0 = element0.basis_function_value(q_point);
00296     std::vector<std::vector<value_type1> > basis_value1 = element1.basis_function_value(q_point);
00297     for (int l = 0;l < n_quadrature_point;l ++) {
00298       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00299       for (int j = 0;j < n_element_dof0;j ++) {
00300         for (int k = 0;k < n_element_dof1;k ++) {
00301           this->elementMatrix(j,k) += Jxw*basis_value0[j][l]*basis_value1[k][l];
00302         }
00303       }
00304     }
00305   }
00306   else {
00307     double volume = element0.templateElement().volume();
00308     const int& alg_acc = this->algebricAccuracy();
00309     const QuadratureInfo<TDIM0>& quad_info = element0.findQuadratureInfo(alg_acc);
00310     std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00311     int n_quadrature_point = quad_info.n_quadraturePoint();
00312     std::vector<afepack::Point<DOW> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00313     std::vector<std::vector<value_type0> > basis_value0 = element0.basis_function_value(q_point);
00314     std::vector<std::vector<value_type1> > basis_value1 = element1.basis_function_value(q_point);
00315     for (int l = 0;l < n_quadrature_point;l ++) {
00316       double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00317       for (int j = 0;j < n_element_dof0;j ++) {
00318         for (int k = 0;k < n_element_dof1;k ++) {
00319           this->elementMatrix(j,k) += Jxw*basis_value0[j][l]*basis_value1[k][l];
00320         }
00321       }
00322     }
00323   }
00324 }
00325 
00326 #undef THIS
00327 #undef TEMPLATE
00328 
00329 #define TEMPLATE template <int DIM, class value_type, int DOW, int TDIM>
00330 #define THIS MassMatrix<DIM,value_type,DOW,TDIM>
00331 
00332 TEMPLATE
00333 void THIS::reinit(FEMSpace<typename MassMatrix::value_type,DIM,DOW,TDIM>& sp)
00334 {
00335   BilinearOperator<DIM,typename MassMatrix::value_type,typename MassMatrix::value_type,DOW,TDIM,TDIM>::reinit(sp, sp);
00336 }
00337 
00338 
00339 
00340 TEMPLATE
00341 void THIS::getElementMatrix(const Element<typename MassMatrix::value_type,DIM,DOW,TDIM>& element0,
00342                             const Element<typename MassMatrix::value_type,DIM,DOW,TDIM>& element1,
00343                             const typename ActiveElementPairIterator<DIM>::State state)
00344 {
00345   const std::vector<int>& ele_dof0 = element0.dof();
00346   const std::vector<int>& ele_dof1 = element1.dof();
00347   int n_element_dof0 = ele_dof0.size();
00348   int n_element_dof1 = ele_dof1.size();
00349   double volume = element0.templateElement().volume();
00350   const int& alg_acc = this->algebricAccuracy();
00351   const QuadratureInfo<DIM>& quad_info = element0.findQuadratureInfo(alg_acc);
00352   std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00353   int n_quadrature_point = quad_info.n_quadraturePoint();
00354   std::vector<afepack::Point<DOW> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00355   std::vector<std::vector<typename MassMatrix::value_type> > basis_value = element0.basis_function_value(q_point);
00356   for (int l = 0;l < n_quadrature_point;l ++) {
00357     double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00358     for (int j = 0;j < n_element_dof0;j ++) {
00359       for (int k = 0;k < n_element_dof1;k ++) {
00360         this->elementMatrix(j,k) += Jxw*basis_value[j][l]*basis_value[k][l];
00361       }
00362     }
00363   }
00364 }
00365 
00366 #undef THIS
00367 #undef TEMPLATE
00368 
00369 #define TEMPLATE template <int DIM, class value_type, int DOW, int TDIM>
00370 #define THIS StiffMatrix<DIM,value_type,DOW,TDIM>
00371 
00372 TEMPLATE
00373 void THIS::reinit(FEMSpace<typename StiffMatrix::value_type,DIM,DOW,TDIM>& sp)
00374 {
00375   BilinearOperator<DIM,typename StiffMatrix::value_type,typename StiffMatrix::value_type,DOW,TDIM,TDIM>::reinit(sp, sp);
00376 }
00377 
00378 
00379 
00380 TEMPLATE
00381 void THIS::getElementMatrix(const Element<typename StiffMatrix::value_type,DIM,DOW,TDIM>& element0,
00382                             const Element<typename StiffMatrix::value_type,DIM,DOW,TDIM>& element1,
00383                             const typename ActiveElementPairIterator<DIM>::State state)
00384 {
00385   const std::vector<int>& ele_dof0 = element0.dof();
00386   const std::vector<int>& ele_dof1 = element1.dof();
00387   int n_element_dof0 = ele_dof0.size();
00388   int n_element_dof1 = ele_dof1.size();
00389   double volume = element0.templateElement().volume();
00390   const int& alg_acc = this->algebricAccuracy();  
00391   const QuadratureInfo<TDIM>& quad_info = element0.findQuadratureInfo(alg_acc);
00392   std::vector<double> jacobian = element0.local_to_global_jacobian(quad_info.quadraturePoint());
00393   int n_quadrature_point = quad_info.n_quadraturePoint();
00394   std::vector<afepack::Point<DOW> > q_point = element0.local_to_global(quad_info.quadraturePoint());
00395   std::vector<std::vector<std::vector<typename StiffMatrix::value_type> > > basis_gradient = element0.basis_function_gradient(q_point);
00396   for (int l = 0;l < n_quadrature_point;l ++) {
00397     double Jxw = quad_info.weight(l)*jacobian[l]*volume;
00398     for (int j = 0;j < n_element_dof0;j ++) {
00399       for (int k = 0;k < n_element_dof1;k ++) {
00400         this->elementMatrix(j,k) += Jxw*innerProduct(basis_gradient[j][l],basis_gradient[k][l]);
00401       }
00402     }
00403   }
00404 }
00405 
00406 #undef THIS
00407 #undef TEMPLATE
00408 
00409 #endif
00410 
00411 //
00412 // end of file