From: Guido Kanschat Date: Thu, 19 Jan 2006 23:06:56 +0000 (+0000) Subject: accelerate restriction matrices by reordering X-Git-Tag: v8.0.0~12570 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fa0a4225db61fd9da9e30997985e97be3dc55325;p=dealii.git accelerate restriction matrices by reordering git-svn-id: https://svn.dealii.org/trunk@12105 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index 7fd78f79d6..bf2165b382 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -300,56 +300,66 @@ FE_RaviartThomas::initialize_restriction() // First, compute interpolation on // subfaces for (unsigned int face=0;face::faces_per_cell;++face) - for (unsigned int sub=0;sub::subfaces_per_face;++sub) - { - // The shape functions of the - // child cell are evaluated - // in the quadrature points - // of a full face. - Quadrature q_face - = QProjector::project_to_face(q_base, face); - // The weight fuctions for - // the coarse face are - // evaluated on the subface - // only. - Quadrature q_sub - = QProjector::project_to_subface(q_base, face, sub); - const unsigned int child - = GeometryInfo::child_cell_on_face(face, sub); - - // On a certain face, we must - // compute the moments of ALL - // fine level functions with - // the coarse level weight - // functions belonging to - // that face. Due to the - // orthogonalization process - // when building the shape - // functions, these weights - // are equal to the - // corresponding shpe - // functions. - for (unsigned int i_child = 0; i_child < this->dofs_per_cell; ++i_child) - for (unsigned int i_face = 0; i_face < this->dofs_per_face; ++i_face) - { - double s = 0.; - // The quadrature - // weights on the - // subcell are NOT - // transformed, so we - // have to do it here. - for (unsigned int k=0;kshape_value_component(i_child, q_face.point(k), - GeometryInfo::unit_normal_direction[face]) - * this->shape_value_component(face*this->dofs_per_face+i_face, - q_sub.point(k), - GeometryInfo::unit_normal_direction[face]); - this->restriction[child](face*this->dofs_per_face+i_face, - i_child) = s; - } - } - + { + // The shape functions of the + // child cell are evaluated + // in the quadrature points + // of a full face. + Quadrature q_face + = QProjector::project_to_face(q_base, face); + // Store shape values, since the + // evaluation suffers if not + // ordered by point + Table<2,double> cached_values(this->dofs_per_cell, q_face.n_quadrature_points); + for (unsigned int k=0;kdofs_per_cell; ++i) + cached_values(i,k) + = this->shape_value_component(i, q_face.point(k), + GeometryInfo::unit_normal_direction[face]); + + for (unsigned int sub=0;sub::subfaces_per_face;++sub) + { + // The weight fuctions for + // the coarse face are + // evaluated on the subface + // only. + Quadrature q_sub + = QProjector::project_to_subface(q_base, face, sub); + const unsigned int child + = GeometryInfo::child_cell_on_face(face, sub); + + // On a certain face, we must + // compute the moments of ALL + // fine level functions with + // the coarse level weight + // functions belonging to + // that face. Due to the + // orthogonalization process + // when building the shape + // functions, these weights + // are equal to the + // corresponding shpe + // functions. + for (unsigned int k=0;kdofs_per_cell; ++i_child) + for (unsigned int i_face = 0; i_face < this->dofs_per_face; ++i_face) + { + // The quadrature + // weights on the + // subcell are NOT + // transformed, so we + // have to do it here. + restriction[child](face*this->dofs_per_face+i_face, + i_child) + += std::pow(.5, dim-1.) * q_sub.weight(k) + * cached_values(i_child, k) + * this->shape_value_component(face*this->dofs_per_face+i_face, + q_sub.point(k), + GeometryInfo::unit_normal_direction[face]); + } + } + } + if (rt_order==0) return; // Create Legendre basis for the @@ -369,23 +379,31 @@ FE_RaviartThomas::initialize_restriction() QGauss q_cell(rt_order+1); const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + + // Store shape values, since the + // evaluation suffers if not + // ordered by point + Table<3,double> cached_values(this->dofs_per_cell, q_cell.n_quadrature_points, dim); + for (unsigned int k=0;kdofs_per_cell; ++i) + for (unsigned int d=0;dshape_value_component(i, q_cell.point(k), d); for (unsigned int child=0;child::children_per_cell;++child) { Quadrature q_sub = QProjector::project_to_child(q_cell, child); - for (unsigned int i_child = 0; i_child < this->dofs_per_cell; ++i_child) - for (unsigned int d=0;dn();++i_weight) - { - double s = 0.; - for (unsigned int k=0;kshape_value_component(i_child, q_cell.point(k), d) - * polynomials[d]->compute_value(i_weight, q_sub.point(k)); - this->restriction[child](start_cell_dofs+i_weight*dim+d, - i_child) = s; - } + for (unsigned int k=0;kdofs_per_cell; ++i_child) + for (unsigned int d=0;dn();++i_weight) + { + restriction[child](start_cell_dofs+i_weight*dim+d, + i_child) + += q_sub.weight(k) + * cached_values(i_child, k, d) + * polynomials[d]->compute_value(i_weight, q_sub.point(k)); + } } for (unsigned int d=0;d