From: kanschat Date: Fri, 9 Jan 2009 20:40:08 +0000 (+0000) Subject: do some cleaning in the RT elements X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=10d3841bb140d6909a2420fa9bce38363d229e9c;p=dealii-svn.git do some cleaning in the RT elements git-svn-id: https://svn.dealii.org/trunk@18164 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 fdb283259d..fde8131c31 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -112,34 +112,6 @@ class FE_RaviartThomas virtual std::string get_name () const; - /** - * Number of base elements in a - * mixed discretization. Here, - * this is of course equal to - * one. - */ - virtual unsigned int n_base_elements () const; - - /** - * Access to base element - * objects. Since this element is - * atomic, base_element(0) is - * @p this, and all other - * indices throw an error. - */ - virtual const FiniteElement & - base_element (const unsigned int index) const; - - /** - * Multiplicity of base element - * @p index. Since this is an - * atomic element, - * element_multiplicity(0) - * returns one, and all other - * indices will throw an error. - */ - virtual unsigned int element_multiplicity (const unsigned int index) const; - /** * Check whether a shape function * may be non-zero on a face. @@ -164,17 +136,6 @@ class FE_RaviartThomas virtual FiniteElement * clone() const; private: - /** - * The order of the - * Raviart-Thomas element. The - * lowest order elements are - * usually referred to as RT0, - * even though their shape - * functions are piecewise - * linears. - */ - const unsigned int rt_order; - /** * Only for internal use. Its * full name is @@ -213,45 +174,6 @@ class FE_RaviartThomas */ void initialize_restriction (); - /** - * Given a set of flags indicating - * what quantities are requested - * from a @p FEValues object, - * return which of these can be - * precomputed once and for - * all. Often, the values of - * shape function at quadrature - * points can be precomputed, for - * example, in which case the - * return value of this function - * would be the logical and of - * the input @p flags and - * @p update_values. - * - * For the present kind of finite - * element, this is exactly the - * case. - */ - virtual UpdateFlags update_once (const UpdateFlags flags) const; - - /** - * This is the opposite to the - * above function: given a set of - * flags indicating what we want - * to know, return which of these - * need to be computed each time - * we visit a new cell. - * - * If for the computation of one - * quantity something else is - * also required (for example, we - * often need the covariant - * transformation when gradients - * need to be computed), include - * this in the result as well. - */ - virtual UpdateFlags update_each (const UpdateFlags flags) const; - /** * Fields of cell-independent data. * @@ -442,8 +364,6 @@ class FE_RaviartThomasNodal virtual FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement &fe_other) const; - virtual UpdateFlags update_once (const UpdateFlags flags) const; - virtual UpdateFlags update_each (const UpdateFlags flags) const; private: /** * Only for internal use. Its diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index 1ce943a7c7..9e1dc328fc 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -860,19 +860,30 @@ template UpdateFlags FE_PolyTensor::update_each (const UpdateFlags flags) const { - const bool values_once = (mapping_type == mapping_none); - UpdateFlags out = update_default; - if (!values_once && (flags & update_values)) - out |= update_values | update_covariant_transformation; - if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation; - if (flags & update_hessians) - out |= update_hessians | update_covariant_transformation; - if (flags & update_cell_normal_vectors) - out |= update_cell_normal_vectors | update_JxW_values; - + switch (mapping_type) + { + case mapping_piola: + { + if (flags & update_values) + out |= update_values | update_piola; + + if (flags & update_gradients) + out |= update_gradients | update_piola | update_covariant_transformation; + + if (flags & update_hessians) + out |= update_hessians | update_piola | update_covariant_transformation; + + break; + } + + default: + { + Assert (false, ExcNotImplemented()); + } + } + return out; } 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 cce73d0fcc..cf158ff4b2 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -35,6 +35,7 @@ DEAL_II_NAMESPACE_OPEN + template FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) : @@ -44,8 +45,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) dim, deg+1, FiniteElementData::Hdiv, 1), std::vector(PolynomialsRaviartThomas::compute_n_pols(deg), true), std::vector >(PolynomialsRaviartThomas::compute_n_pols(deg), - std::vector(dim,true))), - rt_order(deg) + std::vector(dim,true))) { Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; @@ -60,6 +60,12 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) //matrix, generating the correct //basis functions from the raw //ones. + + // We use an auxiliary matrix in + // this function. Therefore, + // inverse_node_matrix is still + // empty and shape_value_component + // returns the 'raw' shape values. FullMatrix M(n_dofs, n_dofs); FETools::compute_node_matrix(M, *this); this->inverse_node_matrix.reinit(n_dofs, n_dofs); @@ -82,9 +88,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i=0; i::max_children_per_face; ++i) face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face); - FETools::compute_face_embedding_matrices(*this, - face_embeddings, - 0, 0); + FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0); this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, this->dofs_per_face); unsigned int target_row=0; @@ -110,14 +114,17 @@ FE_RaviartThomas::get_name () const // this function returns, so they // have to be kept in synch + // note that this->degree is the maximal + // polynomial degree and is thus one higher + // than the argument given to the + // constructor std::ostringstream namebuf; - namebuf << "FE_RaviartThomas<" << dim << ">(" << rt_order << ")"; + namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree-1 << ")"; return namebuf.str(); } - template FiniteElement * FE_RaviartThomas::clone() const @@ -288,7 +295,7 @@ FE_RaviartThomas::initialize_restriction() { const unsigned int iso=RefinementCase::isotropic_refinement-1; - QGauss q_base (rt_order+1); + QGauss q_base (this->degree); const unsigned int n_face_points = q_base.size(); // First, compute interpolation on // subfaces @@ -354,7 +361,7 @@ FE_RaviartThomas::initialize_restriction() } } - if (rt_order==0) return; + if (this->degree == 1) return; // Create Legendre basis for the // space D_xi Q_k. Here, we cannot @@ -364,13 +371,13 @@ FE_RaviartThomas::initialize_restriction() { std::vector > > poly(dim); for (unsigned int d=0;ddegree-1); + poly[dd] = Polynomials::Legendre::generate_complete_basis(this->degree-2); polynomials[dd] = new AnisotropicPolynomials(poly); } - QGauss q_cell(rt_order+1); + QGauss q_cell(this->degree); const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; @@ -444,72 +451,11 @@ FE_RaviartThomas::get_dpo_vector (const unsigned int rt_order) -template -UpdateFlags -FE_RaviartThomas::update_once (const UpdateFlags) const -{ - // even the values have to be - // computed on the real cell, so - // nothing can be done in advance - return update_default; -} - - - -template -UpdateFlags -FE_RaviartThomas::update_each (const UpdateFlags flags) const -{ - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values | update_piola; - - if (flags & update_gradients) - out |= update_gradients | update_piola | update_covariant_transformation; - - if (flags & update_hessians) - out |= update_hessians | update_piola | update_covariant_transformation; - - return out; -} - //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- - - -template -unsigned int -FE_RaviartThomas::n_base_elements () const -{ - return 1; -} - - - -template -const FiniteElement & -FE_RaviartThomas::base_element (const unsigned int index) const -{ - Assert (index==0, ExcIndexRange(index, 0, 1)); - return *this; -} - - - -template -unsigned int -FE_RaviartThomas::element_multiplicity (const unsigned int index) const -{ - Assert (index==0, ExcIndexRange(index, 0, 1)); - return 1; -} - - - template bool FE_RaviartThomas::has_support_on_face (const unsigned int shape_index, @@ -523,9 +469,9 @@ FE_RaviartThomas::has_support_on_face (const unsigned int shape_index, // Return computed values if we // know them easily. Otherwise, it // is always safe to return true. - switch (rt_order) + switch (this->degree) { - case 0: + case 1: { switch (dim) { diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc index 9c0087f7cd..8eb0fa196c 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -25,6 +25,7 @@ #include + DEAL_II_NAMESPACE_OPEN @@ -43,14 +44,12 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) const unsigned int n_dofs = this->dofs_per_cell; this->mapping_type = mapping_raviart_thomas; - // These must be done first, since - // they change the evaluation of - // basis functions - - // Set up the generalized support - // points - initialize_support_points (deg); - //Now compute the inverse node + // First, initialize the + // generalized support points and + // quadrature weights, since they + // are required for interpolation. + initialize_support_points(deg); + // Now compute the inverse node //matrix, generating the correct //basis functions from the raw //ones. @@ -131,6 +130,147 @@ FE_RaviartThomasNodal::clone() const } +//--------------------------------------------------------------------------- +// Auxiliary and internal functions +//--------------------------------------------------------------------------- + + + +template +void +FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) +{ + this->generalized_support_points.resize (this->dofs_per_cell); + this->generalized_face_support_points.resize (this->dofs_per_face); + + // Number of the point being entered + unsigned int current = 0; + + // On the faces, we choose as many + // Gauss points as necessary to + // determine the normal component + // uniquely. This is the deg of + // the Raviart-Thomas element plus + // one. + if (dim>1) + { + QGauss face_points (deg+1); + Assert (face_points.size() == this->dofs_per_face, + ExcInternalError()); + for (unsigned int k=0;kdofs_per_face;++k) + this->generalized_face_support_points[k] = face_points.point(k); + Quadrature faces = QProjector::project_to_all_faces(face_points); + for (unsigned int k=0; + kdofs_per_face*GeometryInfo::faces_per_cell; + ++k) + this->generalized_support_points[k] = faces.point(k+QProjector + ::DataSetDescriptor::face(0, + true, + false, + false, + this->dofs_per_face)); + + current = this->dofs_per_face*GeometryInfo::faces_per_cell; + } + + if (deg==0) return; + // In the interior, we need + // anisotropic Gauss quadratures, + // different for each direction. + QGauss<1> high(deg+1); + QGauss<1> low(deg); + + for (unsigned int d=0;d* quadrature; + if (dim == 1) quadrature = new QAnisotropic(high); + if (dim == 2) quadrature = new QAnisotropic(((d==0) ? low : high), + ((d==1) ? low : high)); + if (dim == 3) quadrature = new QAnisotropic(((d==0) ? low : high), + ((d==1) ? low : high), + ((d==2) ? low : high)); + Assert(dim<=3, ExcNotImplemented()); + + for (unsigned int k=0;ksize();++k) + this->generalized_support_points[current++] = quadrature->point(k); + delete quadrature; + } + Assert (current == this->dofs_per_cell, ExcInternalError()); +} + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(2); + dpo[0] = 1; + dpo[1] = deg; + return dpo; +} + +#endif + + +template +std::vector +FE_RaviartThomasNodal::get_dpo_vector (const unsigned int deg) +{ + // the element is face-based and we have + // (deg+1)^(dim-1) DoFs per face + unsigned int dofs_per_face = 1; + for (unsigned int d=1; d dpo(dim+1); + dpo[dim-1] = dofs_per_face; + dpo[dim] = interior_dofs; + + return dpo; +} + + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int) +{ + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); +} + +#endif + + +template +std::vector +FE_RaviartThomasNodal::get_ria_vector (const unsigned int deg) +{ + const unsigned int dofs_per_cell = PolynomialsRaviartThomas::compute_n_pols(deg); + unsigned int dofs_per_face = deg+1; + for (unsigned int d=2;d 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 void @@ -634,169 +774,6 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, -#if deal_II_dimension == 1 - -template <> -std::vector -FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg) -{ - std::vector dpo(2); - dpo[0] = 1; - dpo[1] = deg; - return dpo; -} - -#endif - - -template -std::vector -FE_RaviartThomasNodal::get_dpo_vector (const unsigned int deg) -{ - // the element is face-based and we have - // (deg+1)^(dim-1) DoFs per face - unsigned int dofs_per_face = 1; - for (unsigned int d=1; d dpo(dim+1); - dpo[dim-1] = dofs_per_face; - dpo[dim] = interior_dofs; - - return dpo; -} - - - -#if deal_II_dimension == 1 - -template <> -std::vector -FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector(); -} - -#endif - - -template -std::vector -FE_RaviartThomasNodal::get_ria_vector (const unsigned int deg) -{ - const unsigned int dofs_per_cell = PolynomialsRaviartThomas::compute_n_pols(deg); - unsigned int dofs_per_face = deg+1; - for (unsigned int d=2;d 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 -void -FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) -{ - this->generalized_support_points.resize (this->dofs_per_cell); - this->generalized_face_support_points.resize (this->dofs_per_face); - - // Number of the point being entered - unsigned int current = 0; - - // On the faces, we choose as many - // Gauss points as necessary to - // determine the normal component - // uniquely. This is the deg of - // the Raviart-Thomas element plus - // one. - if (dim>1) - { - QGauss face_points (deg+1); - Assert (face_points.size() == this->dofs_per_face, - ExcInternalError()); - for (unsigned int k=0;kdofs_per_face;++k) - this->generalized_face_support_points[k] = face_points.point(k); - Quadrature faces = QProjector::project_to_all_faces(face_points); - for (unsigned int k=0; - kdofs_per_face*GeometryInfo::faces_per_cell; - ++k) - this->generalized_support_points[k] = faces.point(k+QProjector - ::DataSetDescriptor::face(0, - true, - false, - false, - this->dofs_per_face)); - - current = this->dofs_per_face*GeometryInfo::faces_per_cell; - } - - if (deg==0) return; - // In the interior, we need - // anisotropic Gauss quadratures, - // different for each direction. - QGauss<1> high(deg+1); - QGauss<1> low(deg); - - for (unsigned int d=0;d* quadrature; - if (dim == 1) quadrature = new QAnisotropic(high); - if (dim == 2) quadrature = new QAnisotropic(((d==0) ? low : high), - ((d==1) ? low : high)); - if (dim == 3) quadrature = new QAnisotropic(((d==0) ? low : high), - ((d==1) ? low : high), - ((d==2) ? low : high)); - Assert(dim<=3, ExcNotImplemented()); - - for (unsigned int k=0;ksize();++k) - this->generalized_support_points[current++] = quadrature->point(k); - delete quadrature; - } - Assert (current == this->dofs_per_cell, ExcInternalError()); -} - - -template -UpdateFlags -FE_RaviartThomasNodal::update_once (const UpdateFlags) const -{ - return update_default; -} - - -template -UpdateFlags -FE_RaviartThomasNodal::update_each (const UpdateFlags flags) const -{ - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values | update_piola; - - if (flags & update_gradients) - out |= update_gradients | update_piola | update_covariant_transformation; - - if (flags & update_hessians) - out |= update_hessians | update_piola | update_covariant_transformation; - - return out; -} - - template class FE_RaviartThomasNodal; DEAL_II_NAMESPACE_CLOSE