From c417172df106750daeea12a79c8b6c92b51bee33 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 7 Apr 2003 15:43:52 +0000 Subject: [PATCH] Add initial versions of these files. git-svn-id: https://svn.dealii.org/trunk@7370 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/fe/fe_raviart_thomas.h | 581 ++++++++ .../deal.II/source/fe/fe_raviart_thomas.cc | 1237 +++++++++++++++++ .../deal.II/source/fe/fe_raviart_thomas_1d.cc | 53 + .../deal.II/source/fe/fe_raviart_thomas_2d.cc | 137 ++ .../deal.II/source/fe/fe_raviart_thomas_3d.cc | 242 ++++ 5 files changed, 2250 insertions(+) create mode 100644 deal.II/deal.II/include/fe/fe_raviart_thomas.h create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas.cc create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc create mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h new file mode 100644 index 0000000000..73ee24cc7a --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -0,0 +1,581 @@ +//--------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002, 2003 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_raviart_thomas_h +#define __deal2__fe_raviart_thomas_h + +#include +#include +#include +#include +#include + +template class TensorProductPolynomials; +template class MappingQ; + + + +/** + * Implementation of continuous Raviart-Thomas elements for the space + * H_div. Note, however, that continuity only concerns the normal + * component of the vector field. + * + * The constructor of this class takes the degree @p{p} of this finite + * element. However, presently, only lowest order elements + * (i.e. @p{p==1}) are implemented. + * + * + * @sect3{Interpolation to finer and coarser meshes} + * + * Each finite element class in deal.II provides matrices that are + * used to interpolate from coarser to finer meshes and the other way + * round. Interpolation from a mother cell to its children is usually + * trivial, since finite element spaces are normally nested and this + * kind of interpolation is therefore exact. On the other hand, when + * we interpolate from child cells to the mother cell, we usually have + * to throw away some information. + * + * For continuous elements, this transfer usually happens by + * interpolating the values on the child cells at the support points + * of the shape functions of the mother cell. However, for + * discontinuous elements, we often use a projection from the child + * cells to the mother cell. The projection approach is only possible + * for discontinuous elements, since it cannot be guaranteed that the + * values of the projected functions on one cell and its neighbor + * match. In this case, only an interpolation can be + * used. (Internally, whether the values of a shape function are + * interpolated or projected, or better: whether the matrices the + * finite element provides are to be treated with the properties of a + * projection or of an interpolation, is controlled by the + * @p{restriction_is_additive} flag. See there for more information.) + * + * Here, things are not so simple: since the element has some + * continuity requirements across faces, we can only resort to some + * kind of interpolation. On the other hand, for the lowest order + * elements, the values of generating functionals are the (constant) + * tangential values of the shape functions. We would therefore really + * like to take the mean value of the tangential values of the child + * faces, and make this the value of the mother face. Then, however, + * taking a mean value of two piecewise constant function is not an + * interpolation, but a restriction. Since this is not possible, we + * cannot use this. + * + * To make a long story somewhat shorter, when interpolating from + * refined edges to a coarse one, we do not take the mean value, but + * pick only one (the one from the first child edge). While this is + * not optimal, it is certainly a valid choice (using an interpolation + * point that is not in the middle of the cell, but shifted to one + * side), and it also preserves the order of the interpolation. + * + * + * @sect3{Numbering of the degrees of freedom (DoFs)} + * + * Nedelec elements have their degrees of freedom on edges, with shape + * functions being vector valued and pointing in tangential + * direction. We use the standard enumeration and direction of edges + * in deal.II, yielding the following shape functions in 2d: + * + * @begin{verbatim} + * 2 + * *--->---* + * | | + * 3^ ^1 + * | | + * *--->---* + * 0 + * @end{verbatim} + * + * For the 3d case, the ordering follows the same scheme: the lines + * are numbered as described in the documentation of the + * @ref{Triangulation} class, i.e. + * @begin{verbatim} + * *---6---* *---6---* + * /| | / /| + * 11 | 5 11 10 5 + * / 7 | / / | + * * | | *---2---* | + * | *---4---* | | * + * | / / | 1 / + * 3 8 9 3 | 9 + * |/ / | |/ + * *---0---* *---0---* + * @end{verbatim} + * and their directions are as follows: + * @begin{verbatim} + * *--->---* *--->---* + * /| | / /| + * ^ | ^ ^ ^ ^ + * / ^ | / / | + * * | | *--->---* | + * | *--->---* | | * + * | / / | ^ / + * ^ ^ ^ ^ | ^ + * |/ / | |/ + * *--->---* *--->---* + * @end{verbatim} + * + * The element does not make much sense in 1d, so it is not + * implemented there. + * + * + * @author Wolfgang Bangerth, 2003 + */ +template +class FE_RaviartThomas : public FiniteElement +{ + public: + /** + * Constructor for the Nedelec + * element of degree @p{p}. + */ + FE_RaviartThomas (const unsigned int p); + + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the second derivative + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the polynomial degree + * of this finite element, + * i.e. the value passed to the + * constructor. + */ + unsigned int get_degree () 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, @p{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, + * @p{element_multiplicity(0)} + * returns one, and all other + * indices will throw an error. + */ + virtual unsigned int element_multiplicity (const unsigned int index) const; + + /** + * This function returns + * @p{true}, if the shape + * function @p{shape_index} has + * non-zero values on the face + * @p{face_index}. For the lowest + * order Nedelec elements, this + * is actually the case for the + * one on which the shape + * function is defined and all + * neighboring ones. + * + * Implementation of the + * interface in + * @ref{FiniteElement} + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + * + * This function is made virtual, + * since finite element objects + * are usually accessed through + * pointers to their base class, + * rather than the class itself. + */ + virtual unsigned int memory_consumption () const; + + + /** + * Declare a nested class which + * will hold static definitions + * of various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files + * @p{fe_raviart_thomas_[23]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Embedding matrices. For + * each element type (the + * first index) there are as + * many embedding matrices as + * there are children per + * cell. The first index + * starts with linear + * elements and goes up in + * polynomial degree. The + * array may grow in the + * future with the number of + * elements for which these + * matrices have been + * computed. If for some + * element, the matrices have + * not been computed then you + * may use the element + * nevertheless but can not + * access the respective + * fields. + */ + static const double * const + embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree for + * which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As the + * @p{embedding_matrices} + * field, but for the + * interface constraints. One + * for each element for which + * it has been computed. + */ + static const double * const constraint_matrices[]; + + /** + * Like + * @p{n_embedding_matrices}, + * but for the number of + * interface constraint + * matrices. + */ + static const unsigned int n_constraint_matrices; + }; + /** + * Exception + */ + DeclException0 (ExcNotUsefulInThisDimension); + + protected: + /** + * @p{clone} function instead of + * a copy constructor. + * + * This function is needed by the + * constructors of @p{FESystem}. + */ + virtual FiniteElement * clone() const; + + /** + * Prepare internal data + * structures and fill in values + * independent of the cell. + */ + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) 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{unit_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_support_points (); + + /** + * Initialize the + * @p{unit_face_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_face_support_points (); + + /** + * 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; + + /** + * Degree of the polynomials. + */ + const unsigned int degree; + + /** + * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. + */ + class InternalData : public FiniteElementBase::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. + */ + Table<2,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. + */ + Table<2,Tensor<2,dim> > shape_gradients; + }; + + /** + * Allow access from other + * dimensions. + */ + template friend class FE_RaviartThomas; +}; + + +/* -------------- declaration of explicit specializations ------------- */ + +template <> void FE_RaviartThomas<1>::initialize_unit_face_support_points (); + +// declaration of explicit specializations of member variables, if the +// compiler allows us to do that (the standard says we must) +#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG +template <> +const double * const +FE_RaviartThomas<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell]; + +template <> +const unsigned int FE_RaviartThomas<1>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_RaviartThomas<1>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_RaviartThomas<1>::Matrices::n_constraint_matrices; + +template <> +const double * const +FE_RaviartThomas<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell]; + +template <> +const unsigned int FE_RaviartThomas<2>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_RaviartThomas<2>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_RaviartThomas<2>::Matrices::n_constraint_matrices; + +template <> +const double * const +FE_RaviartThomas<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell]; + +template <> +const unsigned int FE_RaviartThomas<3>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_RaviartThomas<3>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_RaviartThomas<3>::Matrices::n_constraint_matrices; + +#endif + +#endif diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc new file mode 100644 index 0000000000..64e44f5f5e --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -0,0 +1,1237 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 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 + + +template +FE_RaviartThomas::FE_RaviartThomas (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree), + dim), + std::vector (FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell,false), + std::vector >(FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell, + std::vector(dim,true))), + degree(degree) +{ + Assert (dim >= 2, ExcNotUsefulInThisDimension()); + + // copy constraint matrices if they + // are defined. otherwise leave + // them at zero size + if (degreeinterface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + }; + + // next copy over embedding + // matrices if they are defined + if ((degree < Matrices::n_embedding_matrices+1) && + (Matrices::embedding[degree-1][0] != 0)) + for (unsigned int c=0; c::children_per_cell; ++c) + { + // copy + this->prolongation[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + this->prolongation[c].fill (Matrices::embedding[degree-1][c]); + // and make sure that the row + // sum is 0.5 (for usual + // elements, the row sum must + // be 1, but here the shape + // function is multiplied by + // the inverse of the + // Jacobian, which introduces + // a factor of 1/2 when going + // from mother to child) + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[c](row,col); + Assert (std::fabs(sum-.5) < 1e-14, + ExcInternalError()); + }; + }; + + // then fill restriction + // matrices. they are hardcoded for + // the first few elements + switch (dim) + { + case 2: // 2d + { + switch (degree) + { + case 1: + { + // this is a strange + // element, since it is + // both additive and + // then it is also + // not. ideally, we + // would like to have + // the value of the + // shape function on + // the coarse line to + // be the mean value of + // that on the two + // child ones. thus, + // one should make it + // additive. however, + // additivity only + // works if an element + // does not have any + // continuity + // requirements, since + // otherwise degrees of + // freedom are shared + // between adjacent + // elements, and when + // we make the element + // additive, that would + // mean that we end up + // adding up + // contributions not + // only from the child + // cells of this cell, + // but also from the + // child cells of the + // neighbor, and since + // we cannot know + // whether there even + // exists a neighbor we + // cannot simply make + // the element + // additive. + // + // so, until someone + // comes along with a + // better alternative, + // we do the following: + // make the element + // non-additive, and + // simply pick the + // value of one of the + // child lines for the + // value of the mother + // line (note that we + // have to multiply by + // two, since the shape + // functions scale with + // the inverse + // Jacobian). we thus + // throw away the + // information of one + // of the child lines, + // but there seems to + // be no other way than + // that... + // + // note: to make things + // consistent, and + // restriction + // independent of the + // order in which we + // travel across the + // cells of the coarse + // grid, we have to + // make sure that we + // take the same small + // line when visiting + // its two neighbors, + // to get the value for + // the mother line. we + // take the first line + // always, in the + // canonical direction + // of lines + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + + this->restriction[0](0,0) = 2.; + this->restriction[1](1,1) = 2.; + this->restriction[3](2,2) = 2.; + this->restriction[0](3,3) = 2.; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), leave them + // empty. this does not + // prevent the use of + // this FE, but will + // prevent the use of + // these matrices + break; + }; + }; + + break; + }; + + + case 3: // 3d + { + switch (degree) + { + case 1: + { + // same principle as in + // 2d, take one child + // cell to get at the + // values of each of + // the 12 lines + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + this->restriction[0](0,0) = 2.; + this->restriction[0](3,3) = 2.; + this->restriction[1](1,1) = 2.; + this->restriction[3](2,2) = 2.; + + this->restriction[4](4,4) = 2.; + this->restriction[4](7,7) = 2.; + this->restriction[5](5,5) = 2.; + this->restriction[7](6,6) = 2.; + + this->restriction[0](8,8) = 2.; + this->restriction[1](9,9) = 2.; + this->restriction[2](10,10) = 2.; + this->restriction[3](11,11) = 2.; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), leave them + // empty. this does not + // prevent the use of + // this FE, but will + // prevent the use of + // these matrices + break; + }; + }; + + break; + }; + + default: + Assert (false,ExcNotImplemented()); + } + + // finally fill in support points + // on cell and face + initialize_unit_support_points (); + initialize_unit_face_support_points (); + + // then make + // system_to_component_table + // invalid, since this has no + // meaning for the present element + std::vector > tmp1, tmp2; + this->system_to_component_table.swap (tmp1); + this->face_system_to_component_table.swap (tmp2); +} + + + +template +FiniteElement * +FE_RaviartThomas::clone() const +{ + return new FE_RaviartThomas(degree); +} + + +#if deal_II_dimension == 1 + +template <> +double +FE_RaviartThomas<1>::shape_value_component (const unsigned int , + const Point<1> &, + const unsigned int ) const +{ + Assert (false, ExcNotImplemented()); + return 0.; +} + +#endif + +#if deal_II_dimension == 2 + +template <> +double +FE_RaviartThomas<2>::shape_value_component (const unsigned int i, + const Point<2> &p, + const unsigned int component) const +{ + const unsigned int dim = 2; + + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // first order Raviart-Thomas elements + case 1: + { + switch (i) + { + // (0, 1-y) + case 0: return (component == 0 ? 0: 1-p(1)); + // (x,0) + case 1: return (component == 0 ? p(0) : 0); + // (0, y) + case 2: return (component == 0 ? 0: p(1)); + // (1-x, 0) + case 3: return (component == 0 ? 1-p(0) : 0); + + // there are only + // four shape + // functions!? + default: + Assert (false, ExcInternalError()); + return 0; + }; + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return 0; +} + +#endif + +#if deal_II_dimension == 3 + +template <> +double +FE_RaviartThomas<3>::shape_value_component (const unsigned int i, + const Point<3> &/*p*/, + const unsigned int component) const +{ + const unsigned int dim = 3; + + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return 0; +} + +#endif + +#if deal_II_dimension == 1 + +template <> +Tensor<1,1> +FE_RaviartThomas<1>::shape_grad_component (const unsigned int , + const Point<1> &, + const unsigned int ) const +{ + Assert (false, ExcNotImplemented()); + return Tensor<1,1>(); +} + +#endif + +#if deal_II_dimension == 2 + +template <> +Tensor<1,2> +FE_RaviartThomas<2>::shape_grad_component (const unsigned int i, + const Point<2> &, + const unsigned int component) const +{ + const unsigned int dim = 2; + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // first order Raviart-Thomas elements + case 1: + { + // on the unit cell, the + // gradients of these shape + // functions are constant, so + // we pack them into a table + // for simpler lookup + // + // the format is: first + // index=shape function + // number; second + // index=vector component, + // third index=component + // within gradient + static const double unit_gradients[4][2][2] + = { { {0.,0.} , {0.,-1.} }, + { {1.,0.} , {0.,0.} }, + { {0.,0.} , {0.,+1.} }, + { {-1.,0.}, {0.,0.} } }; + return Tensor<1,dim>(unit_gradients[i][component]); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<1,dim>(); +} + +#endif + +#if deal_II_dimension == 3 + +template <> +Tensor<1,3> +FE_RaviartThomas<3>::shape_grad_component (const unsigned int i, + const Point<3> &/*p*/, + const unsigned int component) const +{ + const unsigned int dim = 3; + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<1,dim>(); +} + +#endif + + +#if deal_II_dimension == 1 + +template <> +Tensor<2,1> +FE_RaviartThomas<1>::shape_grad_grad_component (const unsigned int , + const Point<1> &, + const unsigned int ) const +{ + Assert (false, ExcNotImplemented()); + return Tensor<2,1>(); +} + +#endif + + +#if deal_II_dimension == 2 + +template <> +Tensor<2,2> +FE_RaviartThomas<2>::shape_grad_grad_component (const unsigned int i, + const Point<2> &/*p*/, + const unsigned int component) const +{ + const unsigned int dim = 2; + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // first order Raviart-Thomas + // elements. their second + // derivatives on the unit cell + // are zero + case 1: + { + return Tensor<2,dim>(); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<2,dim>(); +} + +#endif + +#if deal_II_dimension == 3 + +template <> +Tensor<2,3> +FE_RaviartThomas<3>::shape_grad_grad_component (const unsigned int i, + const Point<3> &/*p*/, + const unsigned int component) const +{ + const unsigned int dim = 3; + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component < dim, ExcIndexRange (component, 0, dim)); + + switch (degree) + { + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<2,dim>(); +} + +#endif + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + + +template +void FE_RaviartThomas::initialize_unit_support_points () +{ + switch (degree) + { + case 1: + { + // all degrees of freedom are + // on edges, and their order + // is the same as the edges + // themselves + this->unit_support_points.resize(GeometryInfo::lines_per_cell); + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs right + // between the vertices + // of each line + this->unit_support_points[line] = (vertex_0 + vertex_1) / 2; + }; + + break; + }; + + default: + // no higher order + // elements implemented + // right now + Assert (false, ExcNotImplemented()); + }; +} + + +#if deal_II_dimension == 1 + +template <> +void FE_RaviartThomas<1>::initialize_unit_face_support_points () +{ + // no faces in 1d, so nothing to do +} + +#endif + + +template +void FE_RaviartThomas::initialize_unit_face_support_points () +{ + switch (degree) + { + case 1: + { + // do this the same as above, but + // for one dimension less + this->unit_face_support_points.resize(GeometryInfo::lines_per_cell); + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs right + // between the vertices of each + // line + this->unit_face_support_points[line] = (vertex_0 + vertex_1) / 2; + }; + break; + }; + + default: + // no higher order + // elements implemented + // right now + Assert (false, ExcNotImplemented()); + }; +} + + + +template +std::vector +FE_RaviartThomas::get_dpo_vector(const unsigned int degree) +{ + Assert (degree == 1, ExcNotImplemented()); + + // for degree==1, put all degrees + // of freedom on the lines, and in + // particular @p{degree} DoFs per + // line: + std::vector dpo(dim+1, 0U); + dpo[1] = degree; + + return dpo; +} + + + +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_covariant_transformation; + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_second_derivatives) + out |= update_second_derivatives | update_covariant_transformation; + + return out; +} + + + +//---------------------------------------------------------------------- +// Data field initialization +//---------------------------------------------------------------------- + +template +typename Mapping::InternalDataBase * +FE_RaviartThomas::get_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature) const +{ + // generate a new data object and + // initialize some fields + InternalData* data = new InternalData; + + // check what needs to be + // initialized only once and what + // on every cell/face/subface we + // visit + data->update_once = update_once(update_flags); + data->update_each = update_each(update_flags); + data->update_flags = data->update_once | data->update_each; + + const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory + if (flags & update_values) + data->shape_values.reinit (this->dofs_per_cell, n_q_points); + + if (flags & update_gradients) + data->shape_gradients.reinit (this->dofs_per_cell, n_q_points); + + // if second derivatives through + // finite differencing is required, + // then initialize some objects for + // that + if (flags & update_second_derivatives) + data->initialize_2nd (this, mapping, quadrature); + + // next already fill those fields + // of which we have information by + // now. note that the shape values + // and gradients are only those on + // the unit cell, and need to be + // transformed when visiting an + // actual cell + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int q=0; qshape_values[i][q][c] + = shape_value_component(i,quadrature.point(q),c); + + if (flags & update_gradients) + for (unsigned int c=0; cshape_gradients[i][q][c] + = shape_grad_component(i,quadrature.point(q),c); + } + + return data; +} + + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_RaviartThomas::fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &fe_data = dynamic_cast (fedata); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + Assert (fe_data.shape_values[k].size() == n_q_points, + ExcInternalError()); + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin(), + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + Assert (fe_data.shape_gradients[k].size() == n_q_points, + ExcInternalError()); + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin(), + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, 0, mapping_data, fe_data, data); +} + + + +template +void +FE_RaviartThomas::fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &fe_data = dynamic_cast (fedata); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = face * quadrature.n_quadrature_points; + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin()+offset, + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + // + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin()+offset, + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +template +void +FE_RaviartThomas::fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const +{ + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible + InternalData &fe_data = dynamic_cast (fedata); + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = ((face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points); + + // get the flags indicating the + // fields that have to be filled + const UpdateFlags flags(fe_data.current_update_flags()); + + const unsigned int n_q_points = quadrature.n_quadrature_points; + + // fill shape function + // values. these are vector-valued, + // so we have to transform + // them. since the output format + // (in data.shape_values) is a + // sequence of doubles (one for + // each non-zero shape function + // value, and for each quadrature + // point, rather than a sequence of + // small vectors, we have to use a + // number of conversions + if (flags & update_values) + { + Assert (fe_data.shape_values.n_cols() == + GeometryInfo::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_values (n_q_points); + + Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_values.n_cols() == n_q_points, + ExcInternalError()); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // first transform shape + // values... + mapping.transform_covariant(&*shape_values.begin(), + &*shape_values.end(), + fe_data.shape_values[k].begin()+offset, + mapping_data); + + // then copy over to target: + for (unsigned int q=0; q::faces_per_cell * n_q_points, + ExcInternalError()); + + std::vector > shape_grads1 (n_q_points); + std::vector > shape_grads2 (n_q_points); + + Assert (data.shape_gradients.n_rows() == this->dofs_per_cell * dim, + ExcInternalError()); + Assert (data.shape_gradients.n_cols() == n_q_points, + ExcInternalError()); + + // loop over all shape + // functions, and treat the + // gradients of each shape + // function at all quadrature + // points + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // treat the gradients of + // this particular shape + // function at all + // q-points. if Dv is the + // gradient of the shape + // function on the unit + // cell, then + // (J^-T)Dv(J^-1) is the + // value we want to have on + // the real cell. so, we + // will have to apply a + // covariant transformation + // to Dv twice. since the + // interface only allows + // multiplication with + // (J^-1) from the right, + // we have to trick a + // little in between + // + // do first transformation + mapping.transform_covariant(&*shape_grads1.begin(), + &*shape_grads1.end(), + fe_data.shape_gradients[k].begin()+offset, + mapping_data); + // transpose matrix + for (unsigned int q=0; qcompute_2nd (mapping, cell, offset, mapping_data, fe_data, data); +} + + + +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, + const unsigned int face_index) const +{ + Assert (shape_index < this->dofs_per_cell, + ExcIndexRange (shape_index, 0, this->dofs_per_cell)); + Assert (face_index < GeometryInfo::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); + + switch (degree) + { + case 1: + { + switch (dim) + { + case 2: + { + // only on the one + // non-adjacent face + // are the values + // actually zero. list + // these in a table + const unsigned int + opposite_faces[GeometryInfo<2>::faces_per_cell] + = { 2, 3, 0, 1}; + + return (face_index != opposite_faces[shape_index]); + }; + + default: Assert (false, ExcNotImplemented()); + }; + }; + + default: // other degree + Assert (false, ExcNotImplemented()); + }; + + return true; +} + + + +template +unsigned int +FE_RaviartThomas::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_RaviartThomas::get_degree () const +{ + return degree; +} + + + +template class FE_RaviartThomas; diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc new file mode 100644 index 0000000000..a47a21a4a6 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc @@ -0,0 +1,53 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 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. +// +//---------------------------------------------------------------- + + + +// only compile this file if in 1d. note that Raviart-Thomas elements +// do not make much sense in 1d anyway, so this file only contains +// dummy implementations to avoid linker errors due to missing symbols +#if deal_II_dimension == 1 + + +#include + + +template <> +const double * const +FE_RaviartThomas<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] = +{}; + + +template <> +const unsigned int +FE_RaviartThomas<1>::Matrices::n_embedding_matrices = 0; + + + +// No constraints in 1d +template <> +const unsigned int +FE_RaviartThomas<1>::Matrices::n_constraint_matrices = 0; + + +template <> +const double * const +FE_RaviartThomas<1>::Matrices::constraint_matrices[] = {}; + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {} } +#endif // #if deal_II_dimension == 1 + diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc new file mode 100644 index 0000000000..70e9484b54 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc @@ -0,0 +1,137 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 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. +// +//---------------------------------------------------------------- + + +// only compile this file if in 2d +#if deal_II_dimension == 2 + + +#include + +// Transfer matrices for finite elements: have one matrix for each of +// the four child cells which tells us how the degrees of freedom on +// the child cell are obtained from the degrees of freedom on the +// mother cell +// +// note the following: since the shape functions themselves and not +// only the gradients are transformed using the mapping object from +// the unit cell to the real cell, the actual values of the function +// on the real cell is degree of freedom times value of the shape +// function on the unit cell times inverse Jacobian. Thus, what has +// the DoF value 1 on the mother cell must have the DoF value 1/2 on +// the child cell since the latter is smaller by a (linear scaling) +// factor of two. +namespace FE_RaviartThomas_2d +{ + static const double q1_into_q1_refined_0[] = + { + .5, 0, 0 , 0, + 0, 0.25,0, 0.25, + 0.25, 0, 0.25,0, + 0, 0, 0, .5 + }; + + static const double q1_into_q1_refined_1[] = + { + .5, 0., 0., 0., + 0., .5, 0., 0., + 0.25, 0., 0.25, 0., + 0., 0.25, 0., 0.25, + }; + + static const double q1_into_q1_refined_2[] = + { + 0.25, 0., 0.25, 0., + 0., .5, 0., 0., + 0., 0., .5, 0., + 0., 0.25, 0., 0.25, + }; + + static const double q1_into_q1_refined_3[] = + { + 0.25, 0., 0.25, 0., + 0., 0.25, 0., 0.25, + 0., 0., .5, 0., + 0., 0., 0., .5, + }; +} // namespace FE_RaviartThomas_2d + + +// embedding matrices + +template <> +const double * const +FE_RaviartThomas<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] = +{ + { FE_RaviartThomas_2d::q1_into_q1_refined_0, FE_RaviartThomas_2d::q1_into_q1_refined_1, + FE_RaviartThomas_2d::q1_into_q1_refined_2, FE_RaviartThomas_2d::q1_into_q1_refined_3 } +}; + + +template <> +const unsigned int +FE_RaviartThomas<2>::Matrices::n_embedding_matrices += sizeof(FE_RaviartThomas<2>::Matrices::embedding) / +sizeof(FE_RaviartThomas<2>::Matrices::embedding[0]); + + +// Constraint matrices: how do the new value on child faces depend on +// the values on the mother face if that face has a hanging node +// +// Here, the same applies as for the embedding matrices: since the DoF +// values are not only multiplied by the values of the shape function +// on the unit cell, but also by the transformation, we have to +// multiply the value on the large face by 1/2 to get the same value +// back on the small face. in other words, if a DoF has weight 1 on +// the big cell, then it has to have weight 1/2 on the small ones, in +// order to give the same value of the shape function in real space +namespace FE_RaviartThomas_2d +{ + static const double constraint_q1[] = + { + // the function is constant + // along each edge, so each + // degree of freedom on the + // refined edge has the same + // value as that on the + // coarse edge, modulo the + // issue with the + // transformation described + // above + 1./2., 1./2. + }; + +} + + +template <> +const double * const +FE_RaviartThomas<2>::Matrices::constraint_matrices[] = +{ + FE_RaviartThomas_2d::constraint_q1 +}; + + +template <> +const unsigned int +FE_RaviartThomas<2>::Matrices::n_constraint_matrices += sizeof(FE_RaviartThomas<2>::Matrices::constraint_matrices) / +sizeof(FE_RaviartThomas<2>::Matrices::constraint_matrices[0]); + + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {} } +#endif // #if deal_II_dimension == 2 diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc new file mode 100644 index 0000000000..b5d60f4c43 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc @@ -0,0 +1,242 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003 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. +// +//---------------------------------------------------------------- + +// Transfer matrices for finite elements + + +// only compile this file if in 3d +#if deal_II_dimension == 3 + +#include + +// Transfer matrices for finite elements: have one matrix for each of +// the four child cells which tells us how the degrees of freedom on +// the child cell are obtained from the degrees of freedom on the +// mother cell +// +// note the following: since the shape functions themselves and not +// only the gradients are transformed using the mapping object from +// the unit cell to the real cell, the actual values of the function +// on the real cell is degree of freedom times value of the shape +// function on the unit cell times Jacobian. Thus, what has the DoF +// value 1 on the mother cell must have the DoF value 2 on the child +// cell since the latter is smaller by a (linear scaling) factor of +// two. +namespace FE_RaviartThomas_3d +{ + static const double q1_into_q1_refined_0[] = + { + .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., + 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., + 0., 0., 0., .5, 0., 0., 0., 0.,0.,0.,0.,0., + 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, .0, 0., 0.25, + }; + + static const double q1_into_q1_refined_1[] = + { + .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., + 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., + 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., + 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + }; + + static const double q1_into_q1_refined_2[] = + { + 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., + 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., + 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, + }; + + static const double q1_into_q1_refined_3[] = + { + 0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0., + 0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0., + 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25, + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, + }; + + static const double q1_into_q1_refined_4[] = + { + 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0., 0., 0.25, 0., 0., 0., 0.25, 0., 0., 0., 0., + 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25, + }; + + static const double q1_into_q1_refined_5[] = + { + 0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0., + 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + + }; + + static const double q1_into_q1_refined_6[] = + { + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0., + 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0.25, 0., 0.25, 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, + + + }; + + static const double q1_into_q1_refined_7[] = + { + 0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0., + 0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0., + 0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0., + 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0., 0., 0.25, + 0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, + }; + +} // namespace FE_RaviartThomas_3d + + +// embedding matrices + +template <> +const double * const +FE_RaviartThomas<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] = +{ + { FE_RaviartThomas_3d::q1_into_q1_refined_0, FE_RaviartThomas_3d::q1_into_q1_refined_1, + FE_RaviartThomas_3d::q1_into_q1_refined_2, FE_RaviartThomas_3d::q1_into_q1_refined_3, + FE_RaviartThomas_3d::q1_into_q1_refined_4, FE_RaviartThomas_3d::q1_into_q1_refined_5, + FE_RaviartThomas_3d::q1_into_q1_refined_6, FE_RaviartThomas_3d::q1_into_q1_refined_7 } +}; + + +template <> +const unsigned int +FE_RaviartThomas<3>::Matrices::n_embedding_matrices += sizeof(FE_RaviartThomas<3>::Matrices::embedding) / +sizeof(FE_RaviartThomas<3>::Matrices::embedding[0]); + + + +// Constraint matrices: how do the new value on child faces depend on +// the values on the mother face if that face has a hanging node +// +// Here, the same applies as for the embedding matrices: since the DoF +// values are not only multiplied by the values of the shape function +// on the unit cell, but also by the transformation, we have to +// multiply the value on the large face by 1/2 to get the same value +// back on the small face +namespace FE_RaviartThomas_3d +{ + static const double constraint_q1[] = + { + 0, .25, 0, .25, // first the four interior lines + .25, 0, .25, 0, + 0, .25, 0, .25, + .25, 0, .25, 0, + .5, 0, 0, 0, // then the two child lines of each of the four outer + .5, 0, 0, 0, // ones. since the shape functions are constant on each + 0, .5, 0, 0, // line, the two child lines get the same weights, modulo + 0, .5, 0, 0, // the issue with the division by length scaling + 0, 0, .5, 0, + 0, 0, .5, 0, + 0, 0, 0, .5, + 0, 0, 0, .5 + }; +} + + + +template <> +const double * const +FE_RaviartThomas<3>::Matrices::constraint_matrices[] = +{ + FE_RaviartThomas_3d::constraint_q1 +}; + + + +template <> +const unsigned int +FE_RaviartThomas<3>::Matrices::n_constraint_matrices += sizeof(FE_RaviartThomas<3>::Matrices::constraint_matrices) / +sizeof(FE_RaviartThomas<3>::Matrices::constraint_matrices[0]); + + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {} } +#endif // #if deal_II_dimension == 3 -- 2.39.5