From: kanschat Date: Fri, 17 Sep 2010 17:07:25 +0000 (+0000) Subject: new vector valued DG elements X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d5ae6e9e3ff628c3c4c465f5f9b780116231eb5;p=dealii-svn.git new vector valued DG elements git-svn-id: https://svn.dealii.org/trunk@22018 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_dg_vector.h b/deal.II/deal.II/include/fe/fe_dg_vector.h new file mode 100644 index 0000000000..0400d52374 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_dg_vector.h @@ -0,0 +1,186 @@ +//--------------------------------------------------------------------------- +// $Id: fe_raviart_thomas.h 18166 2009-01-09 21:21:16Z kanschat $ +// +// Copyright (C) 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__fe_dg_vector_h +#define __deal2__fe_dg_vector_h + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +template class MappingQ; + + +/** + * DG elements based on vector valued polynomials. + * + * @ingroup fe + * @author Guido Kanschat + * @date 2010 + */ +template +class FE_DGVector + : + public FE_PolyTensor +{ + public: + /** + * Constructor for the vector + * element of degree @p p. + */ + FE_DGVector (const unsigned int p, MappingType m); + public: + + FiniteElement* clone() const; + + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * FE_RaviartThomas(degree), with + * @p dim and @p degree + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + + + /** + * Check whether a shape function + * may be non-zero on a face. + * + * Returns always + * @p true. + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + virtual void interpolate(std::vector& local_dofs, + const std::vector& values) const; + virtual void interpolate(std::vector& local_dofs, + const std::vector >& values, + unsigned int offset = 0) const; + virtual void interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const; + virtual unsigned int memory_consumption () const; + + private: + /** + * Only for internal use. Its + * full name is + * @p get_dofs_per_object_vector + * function and it creates the + * @p dofs_per_object vector that is + * needed within the constructor to + * be passed to the constructor of + * @p FiniteElementData. + */ + static std::vector + get_dpo_vector (const unsigned int degree); + + /** + * Initialize the @p + * generalized_support_points + * field of the FiniteElement + * class and fill the tables with + * #interior_weights. Called + * from the constructor. + */ + void initialize_support_points (const unsigned int 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 (); + + /** + * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. + */ + class InternalData : public FiniteElement::InternalDataBase + { + public: + /** + * Array with shape function + * values in quadrature + * points. There is one row + * for each shape function, + * containing values for each + * quadrature point. Since + * the shape functions are + * vector-valued (with as + * many components as there + * are space dimensions), the + * value is a tensor. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. The transformation + * to the real space cell is + * then simply done by + * multiplication with the + * Jacobian of the mapping. + */ + std::vector > > shape_values; + + /** + * Array with shape function + * gradients in quadrature + * points. There is one + * row for each shape + * function, containing + * values for each quadrature + * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. + */ + std::vector > > shape_gradients; + }; + Table<3, double> interior_weights; +}; + + +/* -------------- declaration of explicit specializations ------------- */ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/deal.II/include/fe/fe_dg_vector.templates.h b/deal.II/deal.II/include/fe/fe_dg_vector.templates.h new file mode 100644 index 0000000000..46441f3788 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_dg_vector.templates.h @@ -0,0 +1,106 @@ +//--------------------------------------------------------------------------- +// $Id: fe_poly.templates.h 21954 2010-09-13 20:49:06Z kanschat $ +// +// Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +//TODO:[GK] deg+1 is wrong here and should bew fixed after FiniteElementData was cleaned up + +template +FE_DGVector::FE_DGVector ( + const unsigned int deg, MappingType map) + : + FE_PolyTensor( + deg, + FiniteElementData( + get_dpo_vector(deg), dim, deg+1, FiniteElementData::L2, 1), + std::vector(POLY::compute_n_pols(deg), true), + std::vector >(POLY::compute_n_pols(deg), + std::vector(dim,true))) +{ + this->mapping_type = map; + const unsigned int polynomial_degree = this->tensor_degree(); + + QGauss quadrature(polynomial_degree+1); + this->generalized_support_points = quadrature.get_points(); + + this->reinit_restriction_and_prolongation_matrices(true, true); + FETools::compute_projection_matrices (*this, this->restriction, true); + FETools::compute_embedding_matrices (*this, this->prolongation, true); +} + + +template +FiniteElement * +FE_DGVector::clone() const +{ + return new FE_DGVector(*this); +} + + +template +std::string +FE_DGVector::get_name() const +{ + std::ostringstream namebuf; + namebuf << "FE_DGVector_" << this->poly_space.name() + << "<" << dim << ">(" << this->degree-1 << ")"; + return namebuf.str(); +} + + +template +std::vector +FE_DGVector::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(dim+1); + dpo[dim] = POLY::compute_n_pols(deg); + + return dpo; +} + + +template +bool +FE_DGVector::has_support_on_face ( + const unsigned int, + const unsigned int) const +{ + return true; +} + + +template +void +FE_DGVector::interpolate( + std::vector&, + const std::vector&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +void +FE_DGVector::interpolate( + std::vector& /*local_dofs*/, + const std::vector >& /*values*/, + unsigned int /*offset*/) const +{ + Assert(false, ExcNotImplemented()); +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index 10860141a4..acb9b373f9 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -477,8 +477,10 @@ class FETools * to use this function mostly. */ template - static void compute_projection_matrices(const FiniteElement &fe, - std::vector > >& matrices); + static void compute_projection_matrices( + const FiniteElement &fe, + std::vector > >& matrices, + const bool isotropic_only = false); /** * Projects scalar data defined in diff --git a/deal.II/deal.II/source/fe/fe_dg_vector.cc b/deal.II/deal.II/source/fe/fe_dg_vector.cc new file mode 100644 index 0000000000..dac4f7c51a --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_dg_vector.cc @@ -0,0 +1,27 @@ +//--------------------------------------------------------------------------- +// $Id: fe_dgp.cc 17866 2008-12-05 22:27:44Z bangerth $ +// +// Copyright (C) 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class FE_DGVector, deal_II_dimension>; +template class FE_DGVector, deal_II_dimension>; +template class FE_DGVector, deal_II_dimension>; +template class FE_DGVector, deal_II_dimension>; + +DEAL_II_NAMESPACE_CLOSE + 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 58cf547741..d9a6f3dfdf 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, 2006, 2007, 2008, 2009 by the deal.II authors +// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -457,8 +457,9 @@ FE_RaviartThomas::get_dpo_vector (const unsigned int deg) template bool -FE_RaviartThomas::has_support_on_face (const unsigned int shape_index, - const unsigned int face_index) const +FE_RaviartThomas::has_support_on_face ( + const unsigned int shape_index, + const unsigned int face_index) const { Assert (shape_index < this->dofs_per_cell, ExcIndexRange (shape_index, 0, this->dofs_per_cell)); diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 7a9d4320c1..76082357d0 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -933,7 +933,7 @@ FETools::compute_face_embedding_matrices(const FiniteElement& fe, template <> void FETools::compute_projection_matrices(const FiniteElement<1,2>&, - std::vector > >& ) + std::vector > >&, bool) { Assert(false, ExcNotImplemented()); } @@ -942,7 +942,7 @@ FETools::compute_projection_matrices(const FiniteElement<1,2>&, template <> void FETools::compute_projection_matrices(const FiniteElement<2,3>&, - std::vector > >& ) + std::vector > >&, bool) { Assert(false, ExcNotImplemented()); } @@ -953,7 +953,8 @@ FETools::compute_projection_matrices(const FiniteElement<2,3>&, template void FETools::compute_projection_matrices(const FiniteElement& fe, - std::vector > >& matrices) + std::vector > >& matrices, + const bool isotropic_only) { const unsigned int n = fe.dofs_per_cell; const unsigned int nd = fe.n_components(); @@ -1005,9 +1006,12 @@ FETools::compute_projection_matrices(const FiniteElement& fe, mass.gauss_jordan(); } - // loop over all possible refinement cases - for (unsigned int ref_case = RefinementCase::cut_x; - ref_case < RefinementCase::isotropic_refinement+1; ++ref_case) + // loop over all possible + // refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) { const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); @@ -2257,7 +2261,7 @@ void FETools::compute_face_embedding_matrices template void FETools::compute_projection_matrices -(const FiniteElement &, std::vector > >&); +(const FiniteElement &, std::vector > >&, bool); template void FETools::interpolate diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0b75000b1e..ceadb44286 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -267,6 +267,13 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.

deal.II

    + +
  1. New: FE_DGVector implements discontinuous elements based on + vector valued polynomial spaces. +
    + (GK 2010/09/17) +

  2. +
  3. Fixed: The methods VectorTools::interpolate_boundary_values and