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