From: Guido Kanschat Date: Sun, 18 Jul 2010 05:02:58 +0000 (+0000) Subject: BDM elements X-Git-Tag: v8.0.0~5807 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=84e8a38fc478a94b60806b526be563d4109bc875;p=dealii.git BDM elements git-svn-id: https://svn.dealii.org/trunk@21515 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_bdm.h b/deal.II/deal.II/include/fe/fe_bdm.h new file mode 100644 index 0000000000..228c09fb1d --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_bdm.h @@ -0,0 +1,124 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 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_bdm_h +#define __deal2__fe_bdm_h + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +/** + * The Brezzi-Douglas-Marini element. + * + *

Degrees of freedom

+ * + * @todo This is for 2D only. + * + * @todo Transformation works only for uniform, Cartesian meshes + * + * The BDM element of order @p p has p+1 degrees of freedom on + * each face. These are implemented as the function values in the + * p+1 Gauss points on each face. + * + * Additionally, for order greater or equal 2, we have additional + * p(p-1), the number of vector valued polynomials in + * Pp, interior degrees of freedom. These are the + * vector function values in the first p(p-1)/2 of the + * p2 Gauss points in the cell. + */ +template +class FE_BDM + : + public FE_PolyTensor, dim> +{ + public: + /** + * Constructor for the BDM + * element of degree @p p. + */ + FE_BDM (const unsigned int p); + + /** + * Return a string that uniquely + * identifies a finite + * element. This class returns + * FE_BDM(degree), with + * @p dim and @p degree + * replaced by appropriate + * values. + */ + virtual std::string get_name () const; + + virtual FiniteElement* clone () 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; + 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); + + /** + * 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 + * FiniteElement::generalized_support_points + * and FiniteElement::generalized_face_support_points + * fields. Called from the + * constructor. + */ + void initialize_support_points (const unsigned int rt_degree); + /** + * The values in the interior + * support points of the + * polynomials needed as test + * functions. The outer vector is + * indexed by quadrature points, + * the inner by the test + * function. + */ + std::vector > test_values; +}; + +#endif diff --git a/deal.II/deal.II/source/fe/fe_bdm.cc b/deal.II/deal.II/source/fe/fe_bdm.cc new file mode 100644 index 0000000000..1822c2788a --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_bdm.cc @@ -0,0 +1,375 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// 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 +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +template +FE_BDM::FE_BDM (const unsigned int deg) + : + FE_PolyTensor, dim> ( + deg, + FiniteElementData(get_dpo_vector(deg), + dim, deg+1, FiniteElementData::Hdiv, 1), + get_ria_vector (deg), + std::vector >(PolynomialsBDM::compute_n_pols(deg), + std::vector(dim,true))) +{ + Assert (dim >= 2, ExcImpossibleInDim(dim)); + const unsigned int n_dofs = this->dofs_per_cell; + + this->mapping_type = mapping_piola; + // 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 + //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); + +// std::cout << std::endl; +// M.print_formatted(std::cout, 2, true); + + this->inverse_node_matrix.reinit(n_dofs, n_dofs); + this->inverse_node_matrix.invert(M); + // From now on, the shape functions + // will be the correct ones, not + // the raw shape functions anymore. + + this->reinit_restriction_and_prolongation_matrices(true, true); + FETools::compute_embedding_matrices (*this, this->prolongation, true); + +// FullMatrix face_embeddings[GeometryInfo::subfaces_per_face]; +// for (unsigned int i=0; i::subfaces_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); +// this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, +// this->dofs_per_face); +// unsigned int target_row=0; +// for (unsigned int d=0;d::subfaces_per_face;++d) +// for (unsigned int i=0;iinterface_constraints(target_row,j) = face_embeddings[d](i,j); +// ++target_row; +// } +} + + + +template +std::string +FE_BDM::get_name () const +{ + // note that the + // FETools::get_fe_from_name + // function depends on the + // particular format of the string + // 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_BDM<" << dim << ">(" << this->degree-1 << ")"; + + return namebuf.str(); +} + + +template +FiniteElement * +FE_BDM::clone() const +{ + return new FE_BDM(*this); +} + + + +template +void +FE_BDM::interpolate( + std::vector&, + const std::vector&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +void +FE_BDM::interpolate( + std::vector&, + const std::vector >&, + unsigned int) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template +void +FE_BDM::interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const +{ + AssertDimension (values.size(), dim); + Assert (values[0].size() == this->generalized_support_points.size(), + ExcDimensionMismatch(values.size(), this->generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + // First do interpolation on + // faces. There, the component + // evaluated depends on the face + // direction and orientation. + unsigned int fbase = 0; + unsigned int f=0; + for (;f::faces_per_cell; + ++f, fbase+=this->dofs_per_face) + { + for (unsigned int i=0;idofs_per_face;++i) + { + local_dofs[fbase+i] = values[GeometryInfo::unit_normal_direction[f]][fbase+i]; + } + } + + // Done for BDM1 + if (fbase == this->dofs_per_cell) return; + + // What's missing are the interior + // degrees of freedom. In each + // point, we take all components of + // the solution. + Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError()); + + // Here, the number of the point + // and of the shape function + // coincides. This will change + // below, since we have more + // support points than test + // functions in the interior. + const unsigned int pbase = fbase; + for (unsigned int d=0;ddofs_per_cell, ExcInternalError()); +} + + + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_BDM<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_BDM::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; d2) + { + interior_dofs *= deg-2; + interior_dofs /= 3; + } + + std::vector dpo(dim+1); + dpo[dim-1] = dofs_per_face; + dpo[dim] = interior_dofs; + + return dpo; +} + + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_BDM<1>::get_ria_vector (const unsigned int) +{ + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); +} + +#endif + + +template +std::vector +FE_BDM::get_ria_vector (const unsigned int deg) +{ + Assert(dim==2, ExcNotImplemented()); + const unsigned int dofs_per_cell = PolynomialsBDM::compute_n_pols(deg); + const unsigned int dofs_per_face = deg+1; + // all 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 +void +FE_BDM::initialize_support_points (const unsigned int deg) +{ + // interior point in 1d + unsigned int npoints = deg; + // interior point in 2d + if (dim >= 2) + { + npoints *= deg; +// npoints /= 2; + } + // interior point in 2d + if (dim >= 3) + { + npoints *= deg; +// npoints /= 3; + } + npoints += GeometryInfo::faces_per_cell * this->dofs_per_face; + + this->generalized_support_points.resize (npoints); + 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 BDM 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<=1) return; + // Although the polynomial space is + // only P_{k-2}, we use the tensor + // product points for Q_{k-2} + QGauss quadrature(deg); + + // Remember where interior points start + const unsigned int ibase=current; +// for (unsigned int k=0;kgeneralized_support_points[current] = quadrature.point(current-ibase); + ++current; + } + Assert(current == npoints, ExcInternalError()); + + // Finaly, compute the values of + // the test functios in the + // interior quadrature points + PolynomialsP poly(deg-2); + + test_values.resize(quadrature.size()); + std::vector > dummy1; + std::vector > dummy2; + + for (unsigned int k=0;k; + diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c7d3e62f95..98c4ba1618 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -150,6 +150,12 @@ inconvenience this causes.

deal.II

    + +
  1. New: Brezzi-Douglas-Marini elements of arbitrary order in FE_BDM. +
    + (GK 2010/07/19) +

    +
  2. Fixed: The FEValues::get_cell() function was unusable from user code