--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/config.h>
+#include <base/table.h>
+#include <base/polynomials_raviart_thomas.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <base/geometry_info.h>
+#include <fe/fe.h>
+#include <fe/fe_poly_tensor.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int spacedim> class MappingQ;
+
+
+/**
+ * DG elements based on vector valued polynomials.
+ *
+ * @ingroup fe
+ * @author Guido Kanschat
+ * @date 2010
+ */
+template <class POLY, int dim, int spacedim=dim>
+class FE_DGVector
+ :
+ public FE_PolyTensor<POLY, dim, spacedim>
+{
+ public:
+ /**
+ * Constructor for the vector
+ * element of degree @p p.
+ */
+ FE_DGVector (const unsigned int p, MappingType m);
+ public:
+
+ FiniteElement<dim, spacedim>* clone() const;
+
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * <tt>FE_RaviartThomas<dim>(degree)</tt>, 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<double>& local_dofs,
+ const std::vector<double>& values) const;
+ virtual void interpolate(std::vector<double>& local_dofs,
+ const std::vector<Vector<double> >& values,
+ unsigned int offset = 0) const;
+ virtual void interpolate(
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& 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<unsigned int>
+ 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<dim>::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<std::vector<Tensor<1,dim> > > 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<std::vector<Tensor<2,dim> > > shape_gradients;
+ };
+ Table<3, double> interior_weights;
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <fe/fe_dg_vector.h>
+#include <fe/fe_tools.h>
+#include <base/quadrature_lib.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+//TODO:[GK] deg+1 is wrong here and should bew fixed after FiniteElementData was cleaned up
+
+template <class POLY, int dim, int spacedim>
+FE_DGVector<POLY,dim,spacedim>::FE_DGVector (
+ const unsigned int deg, MappingType map)
+ :
+ FE_PolyTensor<POLY, dim, spacedim>(
+ deg,
+ FiniteElementData<dim>(
+ get_dpo_vector(deg), dim, deg+1, FiniteElementData<dim>::L2, 1),
+ std::vector<bool>(POLY::compute_n_pols(deg), true),
+ std::vector<std::vector<bool> >(POLY::compute_n_pols(deg),
+ std::vector<bool>(dim,true)))
+{
+ this->mapping_type = map;
+ const unsigned int polynomial_degree = this->tensor_degree();
+
+ QGauss<dim> 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 <class POLY, int dim, int spacedim>
+FiniteElement<dim, spacedim> *
+FE_DGVector<POLY,dim,spacedim>::clone() const
+{
+ return new FE_DGVector<POLY, dim, spacedim>(*this);
+}
+
+
+template <class POLY, int dim, int spacedim>
+std::string
+FE_DGVector<POLY,dim,spacedim>::get_name() const
+{
+ std::ostringstream namebuf;
+ namebuf << "FE_DGVector_" << this->poly_space.name()
+ << "<" << dim << ">(" << this->degree-1 << ")";
+ return namebuf.str();
+}
+
+
+template <class POLY, int dim, int spacedim>
+std::vector<unsigned int>
+FE_DGVector<POLY,dim,spacedim>::get_dpo_vector (const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1);
+ dpo[dim] = POLY::compute_n_pols(deg);
+
+ return dpo;
+}
+
+
+template <class POLY, int dim, int spacedim>
+bool
+FE_DGVector<POLY,dim,spacedim>::has_support_on_face (
+ const unsigned int,
+ const unsigned int) const
+{
+ return true;
+}
+
+
+template <class POLY, int dim, int spacedim>
+void
+FE_DGVector<POLY,dim,spacedim>::interpolate(
+ std::vector<double>&,
+ const std::vector<double>&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template <class POLY, int dim, int spacedim>
+void
+FE_DGVector<POLY,dim,spacedim>::interpolate(
+ std::vector<double>& /*local_dofs*/,
+ const std::vector<Vector<double> >& /*values*/,
+ unsigned int /*offset*/) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+DEAL_II_NAMESPACE_CLOSE
* to use this function mostly.
*/
template <int dim, typename number, int spacedim>
- static void compute_projection_matrices(const FiniteElement<dim,spacedim> &fe,
- std::vector<std::vector<FullMatrix<number> > >& matrices);
+ static void compute_projection_matrices(
+ const FiniteElement<dim,spacedim> &fe,
+ std::vector<std::vector<FullMatrix<number> > >& matrices,
+ const bool isotropic_only = false);
/**
* Projects scalar data defined in
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <fe/fe_dg_vector.templates.h>
+#include <base/polynomials_abf.h>
+#include <base/polynomials_bdm.h>
+#include <base/polynomials_nedelec.h>
+#include <base/polynomials_raviart_thomas.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template class FE_DGVector<PolynomialsABF<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsBDM<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsNedelec<deal_II_dimension>, deal_II_dimension>;
+template class FE_DGVector<PolynomialsRaviartThomas<deal_II_dimension>, deal_II_dimension>;
+
+DEAL_II_NAMESPACE_CLOSE
+
// $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
template <int dim>
bool
-FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
- const unsigned int face_index) const
+FE_RaviartThomas<dim>::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));
template <>
void
FETools::compute_projection_matrices(const FiniteElement<1,2>&,
- std::vector<std::vector<FullMatrix<double> > >& )
+ std::vector<std::vector<FullMatrix<double> > >&, bool)
{
Assert(false, ExcNotImplemented());
}
template <>
void
FETools::compute_projection_matrices(const FiniteElement<2,3>&,
- std::vector<std::vector<FullMatrix<double> > >& )
+ std::vector<std::vector<FullMatrix<double> > >&, bool)
{
Assert(false, ExcNotImplemented());
}
template <int dim, typename number, int spacedim>
void
FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
- std::vector<std::vector<FullMatrix<number> > >& matrices)
+ std::vector<std::vector<FullMatrix<number> > >& matrices,
+ const bool isotropic_only)
{
const unsigned int n = fe.dofs_per_cell;
const unsigned int nd = fe.n_components();
mass.gauss_jordan();
}
- // loop over all possible refinement cases
- for (unsigned int ref_case = RefinementCase<dim>::cut_x;
- ref_case < RefinementCase<dim>::isotropic_refinement+1; ++ref_case)
+ // loop over all possible
+ // refinement cases
+ unsigned int ref_case = (isotropic_only)
+ ? RefinementCase<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+ for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
{
const unsigned int
nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
template
void FETools::compute_projection_matrices<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&);
+(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&, bool);
template
void FETools::interpolate<deal_II_dimension>
<h3>deal.II</h3>
<ol>
+
+ <li><p>New: FE_DGVector implements discontinuous elements based on
+ vector valued polynomial spaces.
+ <br>
+ (GK 2010/09/17)
+ </p></li>
+
<li>
<p>
Fixed: The methods VectorTools::interpolate_boundary_values and