From 3ffdd16511d315a6b494851aedecf400eb2bea77 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Thu, 19 Jan 2006 21:51:06 +0000 Subject: [PATCH] restriction implemented git-svn-id: https://svn.dealii.org/trunk@12103 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_raviart_thomas.h | 25 ++- .../deal.II/source/fe/fe_raviart_thomas.cc | 191 ++++++++++++------ 2 files changed, 145 insertions(+), 71 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index 0a2dc17d27..71b3e9af17 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -183,16 +183,6 @@ class FE_RaviartThomas static std::vector get_dpo_vector (const unsigned int degree); - /** - * Compute the vector used for - * the - * @p restriction_is_additive - * field passed to the base - * class's constructor. - */ - static std::vector - get_ria_vector (const unsigned int degree); - /** * Initialize the @p * generalized_support_points @@ -205,6 +195,19 @@ class FE_RaviartThomas */ void initialize_support_points (const unsigned int rt_degree); + /** + * Initialize the interpolation + * from functions on refined mesh + * cells onto the father + * cell. According to the + * philosophy of the + * Raviart-Thomas element, this + * restriction operator preserves + * the divergence of a function + * weakly. + */ + void initialize_restriction (); + /** * Given a set of flags indicating * what quantities are requested 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 2e21ecedb6..d50f499d1c 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -41,11 +41,9 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) deg, FiniteElementData(get_dpo_vector(deg), dim, deg+1, FiniteElementData::Hdiv), - get_ria_vector (deg), - std::vector >( - FiniteElementData(get_dpo_vector(deg), - dim,deg+1).dofs_per_cell, - std::vector(dim,true))), + std::vector(PolynomialsRaviartThomas::compute_n_pols(deg), true), + std::vector >(PolynomialsRaviartThomas::compute_n_pols(deg), + std::vector(dim,true))), rt_order(deg) { Assert (dim >= 2, ExcImpossibleInDim(dim)); @@ -73,10 +71,11 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) for (unsigned int i=0; i::children_per_cell; ++i) { this->prolongation[i].reinit (n_dofs, n_dofs); - this->restriction[i].reinit (0, 0); + this->restriction[i].reinit (n_dofs, n_dofs); } FETools::compute_embedding_matrices (*this, &this->prolongation[0]); + initialize_restriction(); std::vector > face_embeddings(1<<(dim-1), FullMatrix(this->dofs_per_face, @@ -270,6 +269,131 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) #endif +#if deal_II_dimension == 1 + +template +void +FE_RaviartThomas::initialize_restriction() +{ + for (unsigned int i=0;i::children_per_cell;++i) + 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_RaviartThomas::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) + 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]); + restriction[child](face*this->dofs_per_face+i_face, + i_child) = s; + } + } + + 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; + + 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)); + restriction[child](start_cell_dofs+i_weight*dim+d, + i_child) = s; + } + } + + for (unsigned int d=0;d @@ -309,59 +433,6 @@ FE_RaviartThomas::get_dpo_vector (const unsigned int rt_order) -#if deal_II_dimension == 1 - -template <> -std::vector -FE_RaviartThomas<1>::get_ria_vector (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector(); -} - -#endif - - -template -std::vector -FE_RaviartThomas::get_ria_vector (const unsigned int rt_order) -{ - unsigned int dofs_per_cell, dofs_per_face; - switch (dim) - { - case 2: - dofs_per_face = rt_order+1; - dofs_per_cell = 2*(rt_order+1)*(rt_order+2); - break; - case 3: - dofs_per_face = (rt_order+1)*(rt_order+1); - dofs_per_cell = 3*(rt_order+1)*(rt_order+1)*(rt_order+2); - break; - default: - Assert (false, ExcNotImplemented()); - } - Assert (FiniteElementData(get_dpo_vector(rt_order),dim).dofs_per_cell == - dofs_per_cell, - ExcInternalError()); - Assert (FiniteElementData(get_dpo_vector(rt_order),dim).dofs_per_face == - dofs_per_face, - ExcInternalError()); - - // all face dofs need to be - // non-additive, since they have - // continuity requirements. - // however, the interior dofs are - // made additive - std::vector ret_val(dofs_per_cell,false); - for (unsigned int i=GeometryInfo::faces_per_cell*dofs_per_face; - i < dofs_per_cell; ++i) - ret_val[i] = true; - - return ret_val; -} - - - template UpdateFlags FE_RaviartThomas::update_once (const UpdateFlags) const -- 2.39.5