From fea705afd4de89870c6e7e9abefec463693be1a6 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 7 Apr 2003 15:47:33 +0000 Subject: [PATCH] Remove again from main branch until element is actually tested a little. git-svn-id: https://svn.dealii.org/trunk@7371 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 deletions(-) delete mode 100644 deal.II/deal.II/include/fe/fe_raviart_thomas.h delete mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas.cc delete mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc delete mode 100644 deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc delete 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 deleted file mode 100644 index 73ee24cc7a..0000000000 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ /dev/null @@ -1,581 +0,0 @@ -//--------------------------------------------------------------- -// $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 deleted file mode 100644 index 64e44f5f5e..0000000000 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ /dev/null @@ -1,1237 +0,0 @@ -//---------------------------------------------------------------- -// $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 deleted file mode 100644 index a47a21a4a6..0000000000 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc +++ /dev/null @@ -1,53 +0,0 @@ -//---------------------------------------------------------------- -// $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 deleted file mode 100644 index 70e9484b54..0000000000 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc +++ /dev/null @@ -1,137 +0,0 @@ -//---------------------------------------------------------------- -// $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 deleted file mode 100644 index b5d60f4c43..0000000000 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc +++ /dev/null @@ -1,242 +0,0 @@ -//---------------------------------------------------------------- -// $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