# include <strstream>
#endif
-
-// namespace for some functions that are used in this file. they are
-// specific to numbering conventions used for the FE_RT element, and
-// are thus not very interesting to the outside world
-namespace
-{
- // auxiliary type to allow for some
- // kind of explicit template
- // specialization of the following
- // functions
- template <int dim> struct int2type {};
-
-
- // generate the j-th out of a total
- // of N points on the unit square
- // in 2d. N needs not be a square
- // number, but must be the product
- // of two integers
- //
- // there is one complication: we
- // want to generate interpolation
- // points on the unit square for
- // the shape functions for this
- // element, but for that we need to
- // make sure that these
- // interpolation points make the
- // resulting matrix rows linearly
- // independent. this is a problem
- // since we have anisotropic
- // polynomials, so for example for
- // the lowest order elements, we
- // have as polynomials in for the
- // x-component of the shape
- // functions only "x" and "1-x",
- // i.e. no y-dependence. if we
- // select as interpolation points
- // the points (.5,0) and (.5,1),
- // we're hosed!
- //
- // thus, the third parameter gives
- // the coordinate direction in
- // which the polynomial degree is
- // highest. we use this to select
- // interpolation points primarily
- // in this direction then
- inline
- Point<2> generate_unit_point (const unsigned int j,
- const unsigned int N,
- const unsigned int d,
- const int2type<2> &)
- {
- Assert (d<2, ExcInternalError());
-
- // factorize N int N1*N2. note
- // that we always have N1<=N2,
- // since the square root is
- // rounded down
- const unsigned int N1 = static_cast<unsigned int>(std::sqrt(1.*N));
- const unsigned int N2 = N/N1;
- Assert (N1*N2 == N, ExcInternalError());
-
- const unsigned int Nx = (d==0 ? N2 : N1),
- Ny = (d==1 ? N2 : N1);
-
- return Point<2> (Nx == 1 ? .5 : 1.*(j%Nx)/(Nx-1),
- Ny == 1 ? .5 : 1.*(j/Nx)/(Ny-1));
- }
-
-
- // generate the j-th out of a total
- // of N points on the unit cube
- // in 3d. N needs not be a cube
- // number, but must be the product
- // of three integers
- //
- // the same applies as above for
- // the meaning of the parameter "d"
- inline
- Point<3> generate_unit_point (const unsigned int /*j*/,
- const unsigned int N,
- const unsigned int d,
- const int2type<3> &)
- {
- Assert (d<3, ExcInternalError());
-
- const unsigned int N1 = static_cast<unsigned int>(std::pow(1.*N, 1./3.));
- const unsigned int N2 = static_cast<unsigned int>(std::sqrt(1.*N/N1));
- const unsigned int N3 = N/(N1*N2);
- Assert (N1*N2*N3 == N, ExcInternalError());
-
- Assert (false, ExcNotImplemented());
-
- return Point<3> ();
- }
-
-}
-
+#include <iostream>
+using namespace std;
template <int dim>
-FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int rt_order)
+FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
:
- FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(rt_order),
- dim, rt_order+1, FiniteElementData<dim>::Hdiv),
- get_ria_vector (rt_order),
- std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(rt_order),dim,rt_order+1).dofs_per_cell,
- std::vector<bool>(dim,true))),
- rt_order(rt_order),
- polynomials (create_polynomials(rt_order)),
- renumber (compute_renumber(rt_order))
+ FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim> (
+ deg,
+ FiniteElementData<dim>(get_dpo_vector(deg),
+ dim, deg+1, FiniteElementData<dim>::Hdiv),
+ get_ria_vector (deg),
+ std::vector<std::vector<bool> >(
+ FiniteElementData<dim>(get_dpo_vector(deg),
+ dim,deg+1).dofs_per_cell,
+ std::vector<bool>(dim,true))),
+ rt_order(deg)
{
Assert (dim >= 2, ExcImpossibleInDim(dim));
-
- // check formula (III.3.22) in the
- // book by Brezzi & Fortin about
- // the number of degrees of freedom
- // per cell
- Assert (((dim==2) &&
- (this->dofs_per_cell == 2*(rt_order+1)*(rt_order+2)))
- ||
- ((dim==3) &&
- (this->dofs_per_cell == 3*(rt_order+1)*(rt_order+1)*(rt_order+2))),
- ExcInternalError());
- Assert (renumber.size() == this->dofs_per_cell,
- ExcInternalError());
+ const unsigned int n_dofs = this->dofs_per_cell;
+
+ // First, initialize the
+ // generalized support points and
+ // quadrature weights, since they
+ // are required for interpolation.
+ initialize_support_points(deg);
+ // Now compute the inverse node
+ //matrix, generating the correct
+ //basis functions from the raw
+ //ones.
+ FullMatrix<double> M(n_dofs, n_dofs);
+ FETools::compute_node_matrix(M, *this);
+ this->inverse_node_matrix.reinit(n_dofs, n_dofs);
+ this->inverse_node_matrix.invert(M);
+ // From now on, the shape functions
+ // will be the correct ones, not
+ // the raw shape functions anymore.
+
// initialize the various matrices
- initialize_constraints ();
-
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
- this->prolongation[i].reinit (this->dofs_per_cell,
- this->dofs_per_cell);
+ this->prolongation[i].reinit (n_dofs,
+ n_dofs);
FETools::compute_embedding_matrices (*this, &this->prolongation[0]);
- initialize_restriction ();
-
- // finally fill in support points
- // on cell and face
- initialize_support_points (rt_order);
- initialize_face_support_points ();
-
+ std::vector<FullMatrix<double> >
+ face_embeddings(1<<(dim-1), FullMatrix<double>(this->dofs_per_face,
+ this->dofs_per_face));
+ FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0);
+ this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face,
+ this->dofs_per_face);
+ unsigned int target_row=0;
+ for (unsigned int d=0;d<face_embeddings.size();++d)
+ for (unsigned int i=0;i<face_embeddings[d].m();++i)
+ {
+ for (unsigned int j=0;j<face_embeddings[d].n();++j)
+ this->interface_constraints(target_row,j) = face_embeddings[d](i,j);
+ ++target_row;
+ }
+//TODO:[WB] What is this?
// then make
// system_to_component_table
// invalid, since this has no
}
-template <int dim>
-double
-FE_RaviartThomas<dim>::shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- Assert (component < dim, ExcIndexRange (component, 0, dim));
-
- // check whether this shape
- // function has a contribution in
- // this component at all, and if so
- // delegate to the respective
- // polynomial
- if (component == renumber[i].first)
- return polynomials[component].compute_value(renumber[i].second, p);
- else
- return 0;
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_RaviartThomas<dim>::shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- Assert (component < dim, ExcIndexRange (component, 0, dim));
-
- // check whether this shape
- // function has a contribution in
- // this component at all, and if so
- // delegate to the respective
- // polynomial
- if (component == renumber[i].first)
- return polynomials[component].compute_grad(renumber[i].second, p);
- else
- return Tensor<1,dim>();
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_RaviartThomas<dim>::shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- Assert (component < dim, ExcIndexRange (component, 0, dim));
-
- // check whether this shape
- // function has a contribution in
- // this component at all, and if so
- // delegate to the respective
- // polynomial
- if (component == renumber[i].first)
- return polynomials[component].compute_grad_grad(renumber[i].second, p);
- else
- return Tensor<2,dim>();
-}
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::
-get_interpolation_matrix (const FiniteElement<1> &,
- FullMatrix<double> &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-#endif
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::
-get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
- FullMatrix<double> &interpolation_matrix) const
-{
- // this is only implemented, if the
- // source FE is also a
- // Raviart-Thomas element,
- // otherwise throw an exception, as
- // the documentation says
- AssertThrow ((x_source_fe.get_name().find ("FE_RaviartThomas<") == 0)
- ||
- (dynamic_cast<const FE_RaviartThomas<dim>*>(&x_source_fe) != 0),
- typename FiniteElement<dim>::
- ExcInterpolationNotImplemented());
-
- // ok, source is a RT element, so
- // we will be able to do the work
- const FE_RaviartThomas<dim> &source_fe
- = dynamic_cast<const FE_RaviartThomas<dim>&>(x_source_fe);
-
- Assert (interpolation_matrix.m() == this->dofs_per_cell,
- ExcDimensionMismatch (interpolation_matrix.m(),
- this->dofs_per_cell));
- Assert (interpolation_matrix.n() == source_fe.dofs_per_cell,
- ExcDimensionMismatch (interpolation_matrix.m(),
- source_fe.dofs_per_cell));
-
-
- // compute the interpolation
- // matrices in much the same way as
- // we do for the embedding matrices
- // from mother to child.
- const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim;
- Assert (dofs_per_coordinate*dim == this->dofs_per_cell,
- ExcInternalError());
- for (unsigned int d=0; d<dim; ++d)
- Assert (polynomials[d].n() == dofs_per_coordinate, ExcInternalError());
-
- const unsigned int source_dofs_per_coordinate = source_fe.dofs_per_cell/dim;
- Assert (source_dofs_per_coordinate*dim == source_fe.dofs_per_cell,
- ExcInternalError());
- for (unsigned int d=0; d<dim; ++d)
- Assert (source_fe.polynomials[d].n() == source_dofs_per_coordinate, ExcInternalError());
-
- FullMatrix<double> cell_interpolation (dofs_per_coordinate,
- dofs_per_coordinate);
- FullMatrix<double> source_interpolation (dofs_per_coordinate,
- source_dofs_per_coordinate);
- FullMatrix<double> tmp (dofs_per_coordinate,
- source_dofs_per_coordinate);
- for (unsigned int d=0; d<dim; ++d)
- {
- for (unsigned int j=0; j<dofs_per_coordinate; ++j)
- {
- // generate a point on this
- // cell and evaluate the
- // shape functions there
- //
- // see the comment for that
- // function to see why the
- // third parameter is
- // necessary
- const Point<dim> p = generate_unit_point (j, dofs_per_coordinate,
- d, int2type<dim>());
- for (unsigned int i=0; i<dofs_per_coordinate; ++i)
- cell_interpolation(j,i) = polynomials[d].compute_value (i, p);
-
- for (unsigned int i=0; i<source_dofs_per_coordinate; ++i)
- source_interpolation(j,i) = source_fe.polynomials[d].compute_value (i, p);
- }
-
- // then compute the
- // interpolation matrix matrix
- // for this coordinate
- // direction
- cell_interpolation.gauss_jordan ();
- cell_interpolation.mmult (tmp, source_interpolation);
-
- // finally transfer the
- // results for this
- // coordinate into the matrix
- // corresponding to the
- // entire space on this
- // cell. cut off very small
- // values here
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- if (renumber[i].first == d)
- for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
- if (source_fe.renumber[j].first == d)
- if (std::fabs(tmp(renumber[i].second,
- source_fe.renumber[j].second)) > 1e-15)
- interpolation_matrix(i,j) = tmp(renumber[i].second,
- source_fe.renumber[j].second);
- }
-
- // if this were a Lagrange
- // interpolation element, we could
- // make sure that the row sum of
- // each of the matrices is 1 at
- // this point. note that this won't
- // work here, since we are working
- // with hierarchical elements for
- // which the shape functions don't
- // sum up to 1
- //
- // however, we can make sure that
- // only components couple that have
- // the same vector component
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
- Assert ((interpolation_matrix(i,j) == 0.) ||
- (renumber[i].first == source_fe.renumber[j].first),
- ExcInternalError());
-}
-
-
-
//---------------------------------------------------------------------------
// Auxiliary and internal functions
//---------------------------------------------------------------------------
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_constraints ()
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_constraints ()
-{
- const unsigned int dim = 2;
-
- this->interface_constraints.
- TableBase<2,double>::reinit (this->interface_constraints_size());
-
- // this case is too easy, so
- // special case it
- if (rt_order == 0)
- {
- this->interface_constraints(0,0) = this->interface_constraints(1,0) = .5;
- return;
- }
-
- // for higher orders of the
- // Raviart-Thomas element:
-
- // restricted to each face, the
- // normal component of the shape
- // functions is an element of P_{k}
- // (in 2d), or Q_{k} (in 3d), where
- // k is the degree of the element
- //
- // from this, we interpolate
- // between mother and cell
- // face. this is slightly
- // complicated by the fact that we
- // don't use Lagrange interpolation
- // polynomials, but rather
- // hierarchical polynomials, so we
- // can't just use point
- // interpolation. what we do
- // instead is to evaluate at a
- // number of points and then invert
- // the interpolation matrix
-
- // mathematically speaking, this
- // works in the following way: on
- // each subface, we want that
- // finite element solututions from
- // both sides coincide. i.e. if a
- // and b are expansion coefficients
- // for the shape functions from
- // both sides, we seek a relation
- // between x and y such that
- // sum_i a_i phi^c_i(x)
- // == sum_j b_j phi_j(x)
- // for all points x on the
- // interface. here, phi^c_i are the
- // shape functions on the small
- // cell on one side of the face,
- // and phi_j those on the big cell
- // on the other side. To get this
- // relation, it suffices to look at
- // a sufficient number of points
- // for which this has to hold. if
- // there are n functions, then we
- // need n evaluation points, and we
- // choose them equidistantly.
- //
- // what one then gets is a matrix
- // system
- // a A == b B
- // where
- // A_ij = phi^c_i(x_j)
- // B_ij = phi_i(x_j)
- // and the relation we are looking for
- // is
- // a = (A^T)^-1 B^T b
- //
- // below, we build up these
- // matrices, but rather than
- // transposing them after the
- // fact, we do so while building
- // them. A will be
- // subface_interpolation, B will be
- // face_interpolation. note that we
- // build up these matrices for all
- // faces at once, rather than
- // considering them separately. the
- // reason is that we finally will
- // want to have them in this order
- // anyway, as this is the format we
- // need inside deal.II
- const std::vector<Polynomials::Polynomial<double> >
- face_polynomials (Polynomials::Hierarchical::
- generate_complete_basis (rt_order));
- Assert (face_polynomials.size() == this->dofs_per_face, ExcInternalError());
-
- FullMatrix<double> face_interpolation (2*this->dofs_per_face, this->dofs_per_face);
- FullMatrix<double> subface_interpolation (2*this->dofs_per_face, 2*this->dofs_per_face);
-
- // generate the matrix for the
- // evaluation points on the big
- // face, and the corresponding
- // points in the coordinate system
- // of the small face. order the
- // shape functions in the same way
- // we want to have them in the
- // final matrix. extend shape
- // functions on the small faces by
- // zero to the other face on which
- // they are not defined (we do this
- // by simply not considering these
- // entries in the matrix)
- //
- // note the agreeable fact that for
- // this element, all the shape
- // functions we presently care for
- // are face-based (i.e. not vertex
- // shape functions); thus, for this
- // element, we can skip the
- // annoying index shifting for the
- // constraints matrix due to its
- // weird format
- for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- {
- const double p_face (1.*i/rt_order/2 + (subface == 0 ? 0. : .5));
- const double p_subface (1.*i/rt_order);
-
- for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
- face_interpolation(subface*this->dofs_per_face+i,
- j)
- = face_polynomials[j].value(p_face);
- subface_interpolation(subface*this->dofs_per_face+i,
- subface*this->dofs_per_face+j)
- = face_polynomials[j].value(p_subface);
- }
- }
-
- subface_interpolation.gauss_jordan ();
- subface_interpolation.mmult (this->interface_constraints,
- face_interpolation);
-
- // there is one additional thing to
- // be considered: since the shape
- // functions on the real cell
- // contain the Jacobian (actually,
- // the determinant of the inverse),
- // there is an additional factor of
- // 2 when going from the big to the
- // small cell:
- this->interface_constraints *= 1./2;
-
- // finally: constraints become
- // really messy if the matrix in
- // question has some entries that
- // are almost zero, but not
- // quite. this will happen in the
- // above procedure due to
- // round-off. let us simply delete
- // these entries
- for (unsigned int i=0; i<this->interface_constraints.m(); ++i)
- for (unsigned int j=0; j<this->interface_constraints.n(); ++j)
- if (std::fabs(this->interface_constraints(i,j)) < 1e-14)
- this->interface_constraints(i,j) = 0.;
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_constraints ()
-{
- Assert (false, ExcNotImplemented());
-}
-
-#endif
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_restriction ()
-{}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_restriction ()
-{
- const unsigned int dim = 2;
- switch (rt_order)
- {
- case 0:
- {
- // 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<GeometryInfo<dim>::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;
- };
-
-
- case 1:
- {
- for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- this->restriction[c].reinit (this->dofs_per_cell,
- this->dofs_per_cell);
-
- // first set the corner
- // nodes. note that they are
- // non-additive
- this->restriction[0](0,0) = 2.;
- this->restriction[0](6,6) = 2.;
-
- this->restriction[1](1,1) = 2.;
- this->restriction[1](2,2) = 2.;
-
- this->restriction[2](3,3) = 2.;
- this->restriction[2](5,5) = 2.;
-
- this->restriction[3](4,4) = 2.;
- this->restriction[3](7,7) = 2.;
-
- // then also set the bubble
- // nodes. they _are_
- // additive. to understand
- // what's going on, recall
- // that the bubble shape
- // functions have value -1
- // (!) at the center point,
- // by construction of the
- // polynomials, and that the
- // corner nodes have values
- // 1/2 there since they are
- // just the linears, and not
- // some interpolating
- // polynomial
- //
- // (actually, the
- // additive/non-additive
- // business shouldn't make
- // that much of a difference:
- // node 4 on cell 0 and node
- // 0 on cell 3 must have the
- // same value, since normal
- // components are
- // continuous. so we could
- // pick either and make these
- // shape functions
- // non-additive as well. we
- // choose to take the mean
- // value, which should be the
- // same as either value, and
- // make the shape function
- // additive)
- this->restriction[0](10,0) = 1.;
- this->restriction[0](10,4) = -1.;
- this->restriction[3](10,0) = -1.;
- this->restriction[3](10,4) = 1.;
-
- this->restriction[1](11,1) = 1.;
- this->restriction[1](11,5) = -1.;
- this->restriction[2](11,1) = -1.;
- this->restriction[2](11,5) = 1.;
-
- this->restriction[0](8,6) = 1.;
- this->restriction[0](8,2) = -1.;
- this->restriction[1](8,6) = -1.;
- this->restriction[1](8,2) = 1.;
-
- this->restriction[3](9,7) = 1.;
- this->restriction[3](9,3) = -1.;
- this->restriction[2](9,7) = -1.;
- this->restriction[2](9,3) = 1.;
-
- break;
- };
-
- // 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
- };
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_restriction ()
-{
- Assert (false, ExcNotImplemented());
-}
-
-#endif
-
-
#if deal_II_dimension == 1
template <int dim>
boundary_weights.reinit(n_face_points, legendre.n());
- Assert (face_points.n_quadrature_points == this->dofs_per_face,
- ExcInternalError());
+// Assert (face_points.n_quadrature_points == this->dofs_per_face,
+// ExcInternalError());
for (unsigned int k=0;k<n_face_points;++k)
{
* legendre.compute_value(i, face_points.point(k));
}
}
-
+
Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
for (;current<GeometryInfo<dim>::faces_per_cell*n_face_points;
++current)
#endif
-#if deal_II_dimension == 1
-
-template <>
-void FE_RaviartThomas<1>::initialize_face_support_points ()
-{
- // no faces in 1d, so nothing to do
-}
-
-#endif
-
-
-template <int dim>
-void FE_RaviartThomas<dim>::initialize_face_support_points ()
-{
- this->unit_face_support_points.resize (this->dofs_per_face);
-
- // like with cell
- // unit_support_points:
- // associate all of the in
- // the face mid-point, since
- // there is no other useful
- // way
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_face_support_points[i] = (dim == 2 ?
- Point<dim-1>(.5) :
- Point<dim-1>(.5,.5));
-}
-
#if deal_II_dimension == 1
}
-#if deal_II_dimension == 1
-
-template <>
-std::vector<AnisotropicPolynomials<1> >
-FE_RaviartThomas<1>::create_polynomials (const unsigned int)
-{
- Assert (false, ExcImpossibleInDim(1));
- return std::vector<AnisotropicPolynomials<1> > ();
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-std::vector<AnisotropicPolynomials<2> >
-FE_RaviartThomas<2>::create_polynomials (const unsigned int rt_order)
-{
- const unsigned int dim = 2;
-
- // use the fact that the RT(k)
- // spaces are spanned by the
- // functions
- // P_{k+1,k} \times P_{k,k+1},
- // see the book by Brezzi and
- // Fortin
- const std::vector<Polynomials::Polynomial<double> > pols[2]
- = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1),
- Polynomials::Hierarchical::generate_complete_basis (rt_order)};
-
- // create spaces (k+1,k) and (k,k+1)
- std::vector<std::vector<Polynomials::Polynomial<double> > >
- pols_vector_1(dim), pols_vector_2(dim);
- pols_vector_1[0] = pols[0];
- pols_vector_1[1] = pols[1];
-
- pols_vector_2[0] = pols[1];
- pols_vector_2[1] = pols[0];
-
- const AnisotropicPolynomials<dim> anisotropic[dim]
- = { AnisotropicPolynomials<dim> (pols_vector_1),
- AnisotropicPolynomials<dim> (pols_vector_2) };
-
- // work around a stupid bug in
- // gcc2.95 where the compiler
- // complains about reaching the end
- // of a non-void function when we
- // simply return the following
- // object unnamed, rather than
- // first creating a named object
- // and then returning it...
- const std::vector<AnisotropicPolynomials<dim> >
- ret_val (&anisotropic[0], &anisotropic[dim]);
- return ret_val;
-}
-
-#endif
-
-
-#if deal_II_dimension == 3
-
-template <>
-std::vector<AnisotropicPolynomials<3> >
-FE_RaviartThomas<3>::create_polynomials (const unsigned int rt_order)
-{
- const unsigned int dim = 3;
-
- // use the fact that the RT(k)
- // spaces are spanned by the
- // functions
- // P_{k+1,k,k} \times P_{k,k+1,k}
- // \times P_{k,k,k+1},
- // see the book by Brezzi and
- // Fortin
- const std::vector<Polynomials::Polynomial<double> > pols[2]
- = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1),
- Polynomials::Hierarchical::generate_complete_basis (rt_order)};
-
- // create spaces (k+1,k,k),
- // (k,k+1,k) and (k,k,k+1)
- std::vector<std::vector<Polynomials::Polynomial<double> > >
- pols_vector_1(dim), pols_vector_2(dim), pols_vector_3(dim);
- pols_vector_1[0] = pols[0];
- pols_vector_1[1] = pols[1];
- pols_vector_1[2] = pols[1];
-
- pols_vector_2[0] = pols[1];
- pols_vector_2[1] = pols[0];
- pols_vector_2[2] = pols[1];
-
- pols_vector_3[0] = pols[1];
- pols_vector_3[1] = pols[1];
- pols_vector_3[2] = pols[0];
-
- const AnisotropicPolynomials<dim> anisotropic[dim]
- = { AnisotropicPolynomials<dim> (pols_vector_1),
- AnisotropicPolynomials<dim> (pols_vector_2),
- AnisotropicPolynomials<dim> (pols_vector_3) };
-
- // work around a stupid bug in
- // gcc2.95 where the compiler
- // complains about reaching the end
- // of a non-void function when we
- // simply return the following
- // object unnamed, rather than
- // first creating a named object
- // and then returning it...
- const std::vector<AnisotropicPolynomials<dim> >
- ret_val (&anisotropic[0], &anisotropic[dim]);
- return ret_val;
-}
-
-#endif
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<1>::compute_renumber (const unsigned int)
-{
- Assert (false, ExcImpossibleInDim(1));
- return std::vector<std::pair<unsigned int, unsigned int> > ();
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<2>::compute_renumber (const unsigned int rt_order)
-{
- const unsigned int dim = 2;
-
- std::vector<std::pair<unsigned int, unsigned int> > ret_val;
-
- // to explain the following: the
- // first (rt_order+1) shape functions
- // are on face 0, and point in
- // y-direction, so are for the
- // second vector component. then
- // there are (rt_order+1) shape
- // functions on face 1, which is
- // for the x vector component, and
- // so on. since the order of face
- // rt_orders of freedom is arbitrary,
- // we simply use the same order as
- // that provided by the 1d
- // polynomial class on which this
- // element is based. after
- // 4*(rt_order+1), the remaining
- // shape functions are all bubbles,
- // so we can number them in any way
- // we want. we do so by first
- // numbering the x-vectors, then
- // the y-vectors
- //
- // now, we have to find a mapping
- // from the above ordering to:
- // first which vector component
- // they belong to (easy), and
- // second the index within this
- // component as provided by the
- // AnisotropicPolynomials class
- //
- // this is mostly a counting
- // argument, tedious and error
- // prone, and so boring to explain
- // that we rather not try to do so
- // here (it's simple, but boring,
- // as said), aside from a few
- // comments below
-
- // face 0
- for (unsigned int i=0; i<rt_order+1; ++i)
- ret_val.push_back (std::make_pair (1U, i));
-
- // face 1
- for (unsigned int i=0; i<rt_order+1; ++i)
- ret_val.push_back (std::make_pair (0U, (rt_order+2)*i+1));
-
- // face 2
- for (unsigned int i=0; i<rt_order+1; ++i)
- ret_val.push_back (std::make_pair (1U, (rt_order+1)+i));
-
- // face 3
- for (unsigned int i=0; i<rt_order+1; ++i)
- ret_val.push_back (std::make_pair (0U, (rt_order+2)*i));
-
- // then go on with interior bubble
- // functions, first for the
- // x-direction, then for the
- // y-direction
- for (unsigned int x=0; x<rt_order; ++x)
- for (unsigned int y=0; y<rt_order+1; ++y)
- {
- const unsigned int index_in_component = (x+2) + y*(rt_order+2);
- Assert (index_in_component < (rt_order+1)*(rt_order+2),
- ExcInternalError());
- ret_val.push_back (std::make_pair(0U, index_in_component));
- }
- for (unsigned int x=0; x<rt_order+1; ++x)
- for (unsigned int y=0; y<rt_order; ++y)
- {
- const unsigned int index_in_component = 2*(rt_order+1) + y + x*rt_order;
- Assert (index_in_component < (rt_order+1)*(rt_order+2),
- ExcInternalError());
- ret_val.push_back (std::make_pair(1U, index_in_component));
- }
-
-#ifdef DEBUG
- // make sure we have actually used
- // up all elements of the tensor
- // product polynomial
- Assert (ret_val.size() == 2*(rt_order+1)*(rt_order+2),
- ExcInternalError());
- std::vector<bool> test[dim] = { std::vector<bool>(ret_val.size()/dim, false),
- std::vector<bool>(ret_val.size()/dim, false) };
- for (unsigned int i=0; i<ret_val.size(); ++i)
- {
- Assert (ret_val[i].first < dim, ExcInternalError());
- Assert (ret_val[i].second < test[ret_val[i].first].size(),
- ExcInternalError());
- Assert (test[ret_val[i].first][ret_val[i].second] == false,
- ExcInternalError());
-
- test[ret_val[i].first][ret_val[i].second] = true;
- }
- for (unsigned int d=0; d<dim; ++d)
- Assert (std::find (test[d].begin(), test[d].end(), false) == test[d].end(),
- ExcInternalError());
-#endif
-
- return ret_val;
-}
-
-#endif
-
-
-#if deal_II_dimension == 3
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<3>::compute_renumber (const unsigned int /*rt_order*/)
-{
- Assert (false, ExcNotImplemented());
- return std::vector<std::pair<unsigned int, unsigned int> > ();
-}
-
-#endif
-
-
-
template <int dim>
UpdateFlags
return out;
}
-
-
//---------------------------------------------------------------------------
// Data field initialization
//---------------------------------------------------------------------------
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_RaviartThomas<dim>::get_data (const UpdateFlags update_flags,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &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.resize (this->dofs_per_cell,
- std::vector<Tensor<1,dim> > (n_q_points));
-
- if (flags & update_gradients)
- data->shape_gradients.resize (this->dofs_per_cell,
- std::vector<Tensor<2,dim> > (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; i<this->dofs_per_cell; ++i)
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- if (flags & update_values)
- for (unsigned int c=0; c<dim; ++c)
- data->shape_values[i][q][c]
- = shape_value_component(i,quadrature.point(q),c);
-
- if (flags & update_gradients)
- for (unsigned int c=0; c<dim; ++c)
- data->shape_gradients[i][q][c]
- = shape_grad_component(i,quadrature.point(q),c);
- }
-
- return data;
-}
-
-
-
-
-//---------------------------------------------------------------------------
-// Fill data of FEValues
-//---------------------------------------------------------------------------
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_values (const Mapping<dim> &mapping,
- const typename Triangulation<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &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<InternalData &> (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<Tensor<1,dim> > 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; k<this->dofs_per_cell; ++k)
- {
- // first transform shape
- // values...
- Assert (fe_data.shape_values[k].size() == n_q_points,
- ExcInternalError());
- mapping.transform_covariant(fe_data.shape_values[k], 0,
- shape_values,
- mapping_data);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_values[k*dim+d][q] = shape_values[q][d];
- };
- };
-
-
- if (flags & update_gradients)
- {
- std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
- std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
- Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
- ExcInternalError());
- Assert (data.shape_gradients[0].size() == 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; k<this->dofs_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(fe_data.shape_gradients[k], 0,
- shape_grads1,
- mapping_data);
- // transpose matrix
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
- // do second transformation
- mapping.transform_covariant(shape_grads2, 0, shape_grads1,
- mapping_data);
- // transpose back
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell,
- QProjector<dim>::DataSetDescriptor::cell(),
- mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_face_values (const Mapping<dim> &mapping,
- const typename Triangulation<dim>::cell_iterator &cell,
- const unsigned int face,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &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<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- face (face, cell->face_orientation(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.size() == this->dofs_per_cell,
- ExcInternalError());
- Assert (fe_data.shape_values[0].size() ==
- GeometryInfo<dim>::faces_per_cell * n_q_points *
- (dim == 3 ? 2 : 1),
- ExcInternalError());
-
- std::vector<Tensor<1,dim> > 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; k<this->dofs_per_cell; ++k)
- {
- // first transform shape
- // values...
- mapping.transform_covariant(fe_data.shape_values[k], offset,
- shape_values,
- mapping_data);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_values[k*dim+d][q] = shape_values[q][d];
- };
- };
-
-
- if (flags & update_gradients)
- {
- Assert (fe_data.shape_values.size() == this->dofs_per_cell,
- ExcInternalError());
- Assert (fe_data.shape_gradients[0].size() ==
- GeometryInfo<dim>::faces_per_cell * n_q_points *
- (dim == 3 ? 2 : 1),
- ExcInternalError());
-
- std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
- std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
- Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
- ExcInternalError());
- Assert (data.shape_gradients[0].size() == 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; k<this->dofs_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(fe_data.shape_gradients[k], offset,
- shape_grads1,
- mapping_data);
- // transpose matrix
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
- // do second transformation
- mapping.transform_covariant(shape_grads2, 0, shape_grads1,
- mapping_data);
- // transpose back
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
- const typename Triangulation<dim>::cell_iterator &cell,
- const unsigned int face,
- const unsigned int subface,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &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<InternalData &> (fedata);
-
- // offset determines which data set
- // to take (all data sets for all
- // faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset
- = (QProjector<dim>::DataSetDescriptor::
- sub_face (face, subface, cell->face_orientation(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[0].size() ==
- GeometryInfo<dim>::faces_per_cell *
- GeometryInfo<dim>::subfaces_per_face *
- n_q_points,
- ExcInternalError());
-
- std::vector<Tensor<1,dim> > 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; k<this->dofs_per_cell; ++k)
- {
- // first transform shape
- // values...
- mapping.transform_covariant(fe_data.shape_values[k], offset,
- shape_values,
- mapping_data);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_values[k*dim+d][q] = shape_values[q][d];
- };
- };
-
-
- if (flags & update_gradients)
- {
- Assert (fe_data.shape_gradients.size() ==
- GeometryInfo<dim>::faces_per_cell *
- GeometryInfo<dim>::subfaces_per_face *
- n_q_points,
- ExcInternalError());
-
- std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
- std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
- Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
- ExcInternalError());
- Assert (data.shape_gradients[0].size() == 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; k<this->dofs_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(fe_data.shape_gradients[k], offset,
- shape_grads1,
- mapping_data);
- // transpose matrix
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
- // do second transformation
- mapping.transform_covariant(shape_grads2, 0, shape_grads1,
- mapping_data);
- // transpose back
- for (unsigned int q=0; q<n_q_points; ++q)
- shape_grads2[q] = transpose(shape_grads1[q]);
-
- // then copy over to target:
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
- };
- }
-
- if (flags & update_second_derivatives)
- this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
for (unsigned int k=0;k<n_face_points;++k)
for (unsigned int i=0;i<boundary_weights.size(1);++i)
{
- local_dofs[i] += boundary_weights(k,i)
- * values[face*n_face_points+k](offset+GeometryInfo<dim>::unit_normal_direction[face]);
+ local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
+ * values[face*n_face_points+k](GeometryInfo<dim>::unit_normal_direction[face]+offset);
}
+
+ const unsigned start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+ const unsigned start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+
+ for (unsigned int k=0;k<interior_weights.size(0);++k)
+ for (unsigned int i=0;i<interior_weights.size(1);++i)
+ for (unsigned int d=0;d<dim;++d)
+ local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[k+start_cell_points](d+offset);
}
for (unsigned int k=0;k<n_face_points;++k)
for (unsigned int i=0;i<boundary_weights.size(1);++i)
{
- local_dofs[i] += boundary_weights(k,i)
+ local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
* values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+k];
}
+
+ const unsigned start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+ const unsigned start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+
+ for (unsigned int k=0;k<interior_weights.size(0);++k)
+ for (unsigned int i=0;i<interior_weights.size(1);++i)
+ for (unsigned int d=0;d<dim;++d)
+ local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[d][k+start_cell_points];
}
}
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::get_degree () const
-{
- return rt_order;
-}
-
-
-
template class FE_RaviartThomas<deal_II_dimension>;