From: oliver Date: Thu, 27 Apr 2006 16:22:54 +0000 (+0000) Subject: Implemented the ABF functionality for the other "interpolate" method and X-Git-Tag: v8.0.0~11847 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d87489ec2508c07abd047086f294315527993f6;p=dealii.git Implemented the ABF functionality for the other "interpolate" method and copied the "initialise_restriction" method from the FE_RT class. git-svn-id: https://svn.dealii.org/trunk@12925 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc index 80f529a48c..c8eb6b843c 100644 --- a/deal.II/deal.II/source/fe/fe_abf.cc +++ b/deal.II/deal.II/source/fe/fe_abf.cc @@ -62,7 +62,8 @@ FE_ABF::FE_ABF (const unsigned int deg) FullMatrix M(n_dofs, n_dofs); FETools::compute_node_matrix(M, *this); - M.print (std::cout); + //TODO: Remove debugging output + // M.print (std::cout); this->inverse_node_matrix.reinit(n_dofs, n_dofs); this->inverse_node_matrix.invert(M); @@ -79,13 +80,14 @@ FE_ABF::FE_ABF (const unsigned int deg) } FETools::compute_embedding_matrices (*this, &this->prolongation[0]); - // initialize_restriction (); + initialize_restriction (); - // TODO std::vector > face_embeddings(1<<(dim-1), FullMatrix(this->dofs_per_face, this->dofs_per_face)); - //FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0); + // TODO: Something goes wrong there. The error of the least squares fit + // is to large ... + // FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0); this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, this->dofs_per_face); unsigned int target_row=0; @@ -326,6 +328,149 @@ FE_ABF::initialize_support_points (const unsigned int deg) #endif +#if deal_II_dimension == 1 + +template +void +FE_ABF::initialize_restriction() +{ + for (unsigned int i=0;i::children_per_cell;++i) + this->restriction[i].reinit(0,0); +} + +#else + +// This function is the same Raviart-Thomas interpolation performed by +// interpolate. Still, we cannot use interpolate, since it was written +// for smooth functions. Thefunctions interpolated here are not +// smooth, maybe even not continuous. Therefore, we must double the +// number of quadrature points in each direction in order to integrate +// only smooth functions. + +// Then again, the interpolated function is chosen such that the +// moments coincide with the function to be interpolated. + +template +void +FE_ABF::initialize_restriction() +{ + QGauss q_base (rt_order+1); + const unsigned int n_face_points = q_base.n_quadrature_points; + // First, compute interpolation on + // subfaces + for (unsigned int face=0;face::faces_per_cell;++face) + { + // 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. + this->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 + // space D_xi Q_k. Here, we cannot + // use the shape functions + std::vector* > polynomials(dim); + for (unsigned int dd=0;dd > > poly(dim); + for (unsigned int d=0;d(poly); + } + + 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 k=0;kdofs_per_cell; ++i_child) + for (unsigned int d=0;dn();++i_weight) + { + this->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 @@ -523,9 +668,33 @@ FE_ABF::interpolate( for (unsigned int d=0;d::faces_per_cell; ++face) + { + double n_orient = (double) GeometryInfo::unit_normal_orientation[face]; + for (unsigned int fp=0; fp < n_face_points; ++fp) + { + // TODO: Check what the face_orientation has to be in 3D + unsigned int k = QProjector::DataSetDescriptor::face (face, false, n_face_points); + for (unsigned int i=0; i::unit_normal_direction[face]+offset); + } + } + // TODO: Check if this "correction" can be removed. + for (unsigned int i=0; i void