*/
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_dgp_[123]d.cc}
+ * in the source directory.
+ */
+ struct Matrices
+ {
+ /**
+ * Pointers to the embedding
+ * matrices, one for each
+ * polynomial degree starting
+ * from constant elements
+ */
+ static const double * const embedding[][GeometryInfo<dim>::children_per_cell];
+
+ /**
+ * Number of elements (first
+ * index) the above field
+ * has. Equals the highest
+ * polynomial degree plus one
+ * for which the embedding
+ * matrices have been
+ * computed.
+ */
+ static const unsigned int n_embedding_matrices;
+
+ /**
+ * As @p{embedding} but for
+ * projection matrices.
+ */
+ static const double * const projection_matrices[][GeometryInfo<dim>::children_per_cell];
+
+ /**
+ * As
+ * @p{n_embedding_matrices}
+ * but for projection
+ * matrices.
+ */
+ static const unsigned int n_projection_matrices;
+ };
+
protected:
/**
FEValuesData<dim>& data) const ;
private:
-
- /**
- * 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_dgp_[123]d.cc}
- * in the source directory.
- */
- struct Matrices
- {
- /**
- * Pointers to the embedding
- * matrices, one for each
- * polynomial degree starting
- * from constant elements
- */
- static const double * const embedding[][GeometryInfo<dim>::children_per_cell];
-
- /**
- * Number of elements (first
- * index) the above field
- * has. Equals the highest
- * polynomial degree plus one
- * for which the embedding
- * matrices have been
- * computed.
- */
- static const unsigned int n_embedding_matrices;
-
- /**
- * As @p{embedding} but for
- * projection matrices.
- */
- static const double * const projection_matrices[][GeometryInfo<dim>::children_per_cell];
-
- /**
- * As
- * @p{n_embedding_matrices}
- * but for projection
- * matrices.
- */
- static const unsigned int n_projection_matrices;
- };
-
/**
* Only for internal use. Its
*/
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_dgq_[123]d.cc}
+ * in the source directory.
+ */
+ struct Matrices
+ {
+ /**
+ * Pointers to the embedding
+ * matrices, one for each
+ * polynomial degree starting
+ * from constant elements
+ */
+ static const double * const embedding[];
+
+ /**
+ * Number of elements (first
+ * index) the above field
+ * has. Equals the highest
+ * polynomial degree plus one
+ * for which the embedding
+ * matrices have been
+ * computed.
+ */
+ static const unsigned int n_embedding_matrices;
+
+ /**
+ * As @p{embedding} but for
+ * projection matrices.
+ */
+ static const double * const projection_matrices[];
+
+ /**
+ * As
+ * @p{n_embedding_matrices}
+ * but for projection
+ * matrices.
+ */
+ static const unsigned int n_projection_matrices;
+ };
+
+
protected:
/**
FEValuesData<dim>& data) const ;
private:
-
- /**
- * 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_dgq_[123]d.cc}
- * in the source directory.
- */
- struct Matrices
- {
- /**
- * Pointers to the embedding
- * matrices, one for each
- * polynomial degree starting
- * from constant elements
- */
- static const double * const embedding[];
-
- /**
- * Number of elements (first
- * index) the above field
- * has. Equals the highest
- * polynomial degree plus one
- * for which the embedding
- * matrices have been
- * computed.
- */
- static const unsigned int n_embedding_matrices;
-
- /**
- * As @p{embedding} but for
- * projection matrices.
- */
- static const double * const projection_matrices[];
-
- /**
- * As
- * @p{n_embedding_matrices}
- * but for projection
- * matrices.
- */
- static const unsigned int n_projection_matrices;
- };
-
-
/**
* Only for internal use. Its
* full name is
*/
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_nedelec_[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<dim>::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
*/
FEValuesData<dim>& data) const ;
private:
-
- /**
- * 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_nedelec_[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<dim>::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;
- };
/**
* Only for internal use. Its
*/
virtual unsigned int memory_consumption () const;
- protected:
- /**
- * @p{clone} function instead of
- * a copy constructor.
- *
- * This function is needed by the
- * constructors of @p{FESystem}.
- */
- virtual FiniteElement<dim> * clone() const;
-
- /**
- * Prepare internal data
- * structures and fill in values
- * independent of the cell.
- */
- virtual
- typename Mapping<dim>::InternalDataBase *
- get_data (const UpdateFlags,
- const Mapping<dim>& mapping,
- const Quadrature<dim>& quadrature) const ;
-
- /**
- * Implementation of the same
- * function in
- * @ref{FiniteElement}.
- */
- virtual void
- fill_fe_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_internal,
- typename Mapping<dim>::InternalDataBase &fe_internal,
- FEValuesData<dim>& data) const;
-
- /**
- * Implementation of the same
- * function in
- * @ref{FiniteElement}.
- */
- virtual void
- fill_fe_face_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_internal,
- typename Mapping<dim>::InternalDataBase &fe_internal,
- FEValuesData<dim>& data) const ;
-
- /**
- * Implementation of the same
- * function in
- * @ref{FiniteElement}.
- */
- virtual void
- fill_fe_subface_values (const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int sub_no,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_internal,
- typename Mapping<dim>::InternalDataBase &fe_internal,
- FEValuesData<dim>& data) const ;
-
- private:
-
/**
* Declare a nested class which
* will hold static definitions of
*/
static const unsigned int n_constraint_matrices;
};
+
+ protected:
+ /**
+ * @p{clone} function instead of
+ * a copy constructor.
+ *
+ * This function is needed by the
+ * constructors of @p{FESystem}.
+ */
+ virtual FiniteElement<dim> * clone() const;
+
+ /**
+ * Prepare internal data
+ * structures and fill in values
+ * independent of the cell.
+ */
+ virtual
+ typename Mapping<dim>::InternalDataBase *
+ get_data (const UpdateFlags,
+ const Mapping<dim>& mapping,
+ const Quadrature<dim>& quadrature) const ;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_face_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const ;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_subface_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const ;
+
+ private:
/**
* Only for internal use. Its
*/
unsigned int get_degree () const;
- protected:
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values) const ;
-
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim-1>& quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &exterior_form,
- typename std::vector<Point<dim> > &normal_vectors) const ;
-
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int sub_no,
- const Quadrature<dim-1>& quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &exterior_form,
- typename std::vector<Point<dim> > &normal_vectors) const ;
-
/**
* Storage for internal data of
* Q_degree transformation.
typename MappingQ1<dim>::InternalData mapping_q1_data;
};
+ protected:
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values) const ;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim-1>& quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ typename std::vector<Tensor<1,dim> > &exterior_form,
+ typename std::vector<Point<dim> > &normal_vectors) const ;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1>& quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ typename std::vector<Tensor<1,dim> > &exterior_form,
+ typename std::vector<Point<dim> > &normal_vectors) const ;
+
/**
* For @p{dim=2,3}. Append the
* support points of all shape
*/
virtual UpdateFlags update_each (const UpdateFlags) const;
- protected:
-
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim>& quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values) const ;
-
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim-1>& quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &boundary_form,
- typename std::vector<Point<dim> > &normal_vectors) const ;
-
- /**
- * Implementation of the interface in
- * @ref{Mapping}.
- */
- virtual void
- fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int sub_no,
- const Quadrature<dim-1>& quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &boundary_form,
- typename std::vector<Point<dim> > &normal_vectors) const ;
-
/**
* Storage for internal data of
* d-linear transformation.
unsigned int n_shape_functions;
};
+ protected:
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim>& quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values) const ;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim-1>& quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ typename std::vector<Tensor<1,dim> > &boundary_form,
+ typename std::vector<Point<dim> > &normal_vectors) const ;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1>& quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ typename std::vector<Tensor<1,dim> > &boundary_form,
+ typename std::vector<Point<dim> > &normal_vectors) const ;
+
/**
* Compute shape values and/or
* derivatives.
const std::vector<std::vector<bool> > &nonzero_components)
:
FiniteElementData<dim> (fe_data),
- system_to_component_table(dofs_per_cell),
- face_system_to_component_table(dofs_per_face),
- system_to_base_table(dofs_per_cell),
- face_system_to_base_table(dofs_per_face),
- component_to_system_table(components, std::vector<unsigned>(dofs_per_cell)),
- face_component_to_system_table(components, std::vector<unsigned>(dofs_per_face)),
- component_to_base_table (components, std::make_pair(0U, 0U)),
- restriction_is_additive_flags(restriction_is_additive_flags),
- nonzero_components (nonzero_components),
- n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
- cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
- n_nonzero_components_table.end(),
- std::bind2nd(std::not_equal_to<unsigned int>(),
- 1U))
- ==
- n_nonzero_components_table.end())
+ system_to_component_table(this->dofs_per_cell),
+ face_system_to_component_table(this->dofs_per_face),
+ system_to_base_table(this->dofs_per_cell),
+ face_system_to_base_table(this->dofs_per_face),
+ component_to_system_table(this->components,
+ std::vector<unsigned>(this->dofs_per_cell)),
+ face_component_to_system_table(this->components,
+ std::vector<unsigned>(this->dofs_per_face)),
+ component_to_base_table (this->components,
+ std::make_pair(0U, 0U)),
+ restriction_is_additive_flags(restriction_is_additive_flags),
+ nonzero_components (nonzero_components),
+ n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
+ cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
+ n_nonzero_components_table.end(),
+ std::bind2nd(std::not_equal_to<unsigned int>(),
+ 1U))
+ ==
+ n_nonzero_components_table.end())
{
Assert (restriction_is_additive_flags.size()==fe_data.components,
ExcDimensionMismatch(restriction_is_additive_flags.size(),
fe_data.components));
- Assert (nonzero_components.size() == dofs_per_cell,
+ Assert (nonzero_components.size() == this->dofs_per_cell,
ExcInternalError());
for (unsigned int i=0; i<nonzero_components.size(); ++i)
{
- Assert (nonzero_components[i].size() == n_components(),
+ Assert (nonzero_components[i].size() == this->n_components(),
ExcInternalError());
Assert (std::count (nonzero_components[i].begin(),
nonzero_components[i].end(),
ExcInternalError());
Assert (n_nonzero_components_table[i] >= 1,
ExcInternalError());
- Assert (n_nonzero_components_table[i] <= n_components(),
+ Assert (n_nonzero_components_table[i] <= this->n_components(),
ExcInternalError());
};
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
{
- restriction[i].reinit (dofs_per_cell, dofs_per_cell);
- prolongation[i].reinit (dofs_per_cell, dofs_per_cell);
+ restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell);
+ prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell);
};
// first set sizes of some
switch (dim)
{
case 1:
- Assert ((interface_constraints.m() == 0) &&
- (interface_constraints.n() == 0),
- ExcInternalError());
- break;
+ Assert ((interface_constraints.m() == 0) &&
+ (interface_constraints.n() == 0),
+ ExcInternalError());
+ break;
case 2:
- interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line,
- dofs_per_face);
- break;
+ interface_constraints.reinit (this->dofs_per_vertex
+ +2*this->dofs_per_line,
+ this->dofs_per_face);
+ break;
case 3:
- interface_constraints.reinit (5*dofs_per_vertex +
- 12*dofs_per_line +
- 4*dofs_per_quad,
- dofs_per_face);
- break;
+ interface_constraints.reinit (5*this->dofs_per_vertex +
+ 12*this->dofs_per_line +
+ 4*this->dofs_per_quad,
+ this->dofs_per_face);
+ break;
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
// this is the default way, if there is only
// one component; if there are several, then
// the constructor of the derived class needs
// to fill these arrays
- for (unsigned int j=0 ; j<dofs_per_cell ; ++j)
+ for (unsigned int j=0 ; j<this->dofs_per_cell ; ++j)
{
system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
component_to_system_table[0][j] = j;
}
- for (unsigned int j=0 ; j<dofs_per_face ; ++j)
+ for (unsigned int j=0 ; j<this->dofs_per_face ; ++j)
{
face_system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
const FullMatrix<double> &
FiniteElementBase<dim>::constraints () const
{
- Assert ((dofs_per_face == 0) || (interface_constraints.m() != 0),
+ Assert ((this->dofs_per_face == 0) || (interface_constraints.m() != 0),
ExcConstraintsVoid());
if (dim==1)
// there are as many as there are
// degrees of freedom
Assert ((unit_support_points.size() == 0) ||
- (unit_support_points.size() == dofs_per_cell),
+ (unit_support_points.size() == this->dofs_per_cell),
ExcInternalError());
return unit_support_points;
};
// there are as many as there are
// degrees of freedom on a face
Assert ((unit_face_support_points.size() == 0) ||
- (unit_face_support_points.size() == dofs_per_face),
+ (unit_face_support_points.size() == this->dofs_per_face),
ExcInternalError());
return unit_face_support_points;
};
Assert ((fe_internal.update_each | fe_internal.update_once)
& update_second_derivatives,
ExcInternalError());
- Assert (data.shape_2nd_derivatives.size() == dofs_per_cell,
+ Assert (data.shape_2nd_derivatives.size() == this->dofs_per_cell,
ExcInternalError());
// Number of quadrature points
const unsigned int n_q_points = data.shape_2nd_derivatives[0].size();
// for all shape functions at all
// quadrature points and difference
// quotients in all directions:
- for (unsigned int shape=0; shape<dofs_per_cell; ++shape)
+ for (unsigned int shape=0; shape<this->dofs_per_cell; ++shape)
{
for (unsigned int d1=0; d1<dim; ++d1)
for (unsigned int q=0; q<n_q_points; ++q)
std::vector<bool>(1,true),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
- degree(degree),
- polynomial_space (Legendre<double>::generate_complete_basis(degree))
+ degree(degree),
+ polynomial_space (Legendre<double>::generate_complete_basis(degree))
{
// if defined, copy over matrices
// from precomputed arrays
if ((degree < Matrices::n_embedding_matrices) &&
(Matrices::embedding[degree][0] != 0))
for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
- prolongation[c].fill (Matrices::embedding[degree][c]);
+ this->prolongation[c].fill (Matrices::embedding[degree][c]);
else
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- prolongation[i].reinit(0,0);
+ this->prolongation[i].reinit(0,0);
// // same as above: copy over matrix
// // from predefined values and
FE_DGP<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_value(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_value(i, p);
}
FE_DGP<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_grad(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad(i, p);
}
FE_DGP<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_grad_grad(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad_grad(i, p);
}
// allocate memory
if (flags & update_values)
{
- values.resize (dofs_per_cell);
- data->shape_values.resize(dofs_per_cell,
+ values.resize (this->dofs_per_cell);
+ data->shape_values.resize(this->dofs_per_cell,
std::vector<double>(n_q_points));
}
if (flags & update_gradients)
{
- grads.resize (dofs_per_cell);
- data->shape_gradients.resize(dofs_per_cell,
+ grads.resize (this->dofs_per_cell);
+ data->shape_gradients.resize(this->dofs_per_cell,
std::vector<Tensor<1,dim> >(n_q_points));
}
{
polynomial_space.compute(quadrature.point(i),
values, grads, grad_grads);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
if (flags & update_values)
data->shape_values[k][i] = values[k];
const UpdateFlags flags(fe_data.current_update_flags());
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
if (flags & update_values)
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
std::vector<bool>(1,true),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
- degree(degree),
- polynomial_space (LagrangeEquidistant::generate_complete_basis(degree))
+ degree(degree),
+ polynomial_space (LagrangeEquidistant::generate_complete_basis(degree))
{
// generate permutation/rotation
// index sets to generate some
if ((degree < Matrices::n_embedding_matrices) &&
(Matrices::embedding[degree] != 0))
{
- prolongation[0].fill (Matrices::embedding[degree]);
+ this->prolongation[0].fill (Matrices::embedding[degree]);
switch (dim)
{
case 1:
- prolongation[1].fill_permutation (prolongation[0],
- right, right);
- break;
+ this->prolongation[1].fill_permutation (this->prolongation[0],
+ right, right);
+ break;
case 2:
- prolongation[1].fill_permutation (prolongation[0],
- right, right);
- prolongation[2].fill_permutation (prolongation[1],
- right, right);
- prolongation[3].fill_permutation (prolongation[2],
- right, right);
- break;
+ this->prolongation[1].fill_permutation (this->prolongation[0],
+ right, right);
+ this->prolongation[2].fill_permutation (this->prolongation[1],
+ right, right);
+ this->prolongation[3].fill_permutation (this->prolongation[2],
+ right, right);
+ break;
case 3:
- prolongation[1].fill_permutation (prolongation[0],
- right, right);
- prolongation[5].fill_permutation (prolongation[1],
- right, right);
- prolongation[4].fill_permutation (prolongation[5],
- right, right);
- prolongation[7].fill_permutation (prolongation[4],
- top, top);
- prolongation[3].fill_permutation (prolongation[7],
- top, top);
- prolongation[6].fill_permutation (prolongation[5],
- top, top);
- prolongation[2].fill_permutation (prolongation[6],
- top, top);
- break;
+ this->prolongation[1].fill_permutation (this->prolongation[0],
+ right, right);
+ this->prolongation[5].fill_permutation (this->prolongation[1],
+ right, right);
+ this->prolongation[4].fill_permutation (this->prolongation[5],
+ right, right);
+ this->prolongation[7].fill_permutation (this->prolongation[4],
+ top, top);
+ this->prolongation[3].fill_permutation (this->prolongation[7],
+ top, top);
+ this->prolongation[6].fill_permutation (this->prolongation[5],
+ top, top);
+ this->prolongation[2].fill_permutation (this->prolongation[6],
+ top, top);
+ break;
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
}
else
// matrix undefined, set size to zero
for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
- prolongation[i].reinit(0);
+ this->prolongation[i].reinit(0);
// same as above: copy over matrix
// from predefined values and
if ((degree < Matrices::n_projection_matrices) &&
(Matrices::projection_matrices[degree] != 0))
{
- restriction[0].fill (Matrices::projection_matrices[degree]);
+ this->restriction[0].fill (Matrices::projection_matrices[degree]);
switch (dim)
{
case 1:
- restriction[1].fill_permutation (restriction[0],
- right, right);
- break;
+ this->restriction[1].fill_permutation (this->restriction[0],
+ right, right);
+ break;
case 2:
- restriction[1].fill_permutation (restriction[0],
- right, right);
- restriction[2].fill_permutation (restriction[1],
- right, right);
- restriction[3].fill_permutation (restriction[2],
- right, right);
- break;
+ this->restriction[1].fill_permutation (this->restriction[0],
+ right, right);
+ this->restriction[2].fill_permutation (this->restriction[1],
+ right, right);
+ this->restriction[3].fill_permutation (this->restriction[2],
+ right, right);
+ break;
case 3:
- restriction[1].fill_permutation (restriction[0],
- right, right);
- restriction[5].fill_permutation (restriction[1],
- right, right);
- restriction[4].fill_permutation (restriction[5],
- right, right);
- restriction[7].fill_permutation (restriction[4],
- top, top);
- restriction[3].fill_permutation (restriction[7],
- top, top);
- restriction[6].fill_permutation (restriction[5],
- top, top);
- restriction[2].fill_permutation (restriction[6],
- top, top);
- break;
+ this->restriction[1].fill_permutation (this->restriction[0],
+ right, right);
+ this->restriction[5].fill_permutation (this->restriction[1],
+ right, right);
+ this->restriction[4].fill_permutation (this->restriction[5],
+ right, right);
+ this->restriction[7].fill_permutation (this->restriction[4],
+ top, top);
+ this->restriction[3].fill_permutation (this->restriction[7],
+ top, top);
+ this->restriction[6].fill_permutation (this->restriction[5],
+ top, top);
+ this->restriction[2].fill_permutation (this->restriction[6],
+ top, top);
+ break;
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
}
else
// matrix undefined, set size to zero
for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
- restriction[i].reinit(0);
+ this->restriction[i].reinit(0);
// finally fill in support points
{
// constant elements, take
// midpoint
- unit_support_points.resize(1);
+ this->unit_support_points.resize(1);
for (unsigned int i=0; i<dim; ++i)
- unit_support_points[0](i) = 0.5;
+ this->unit_support_points[0](i) = 0.5;
}
else
{
for (unsigned int i=1; i<dim; ++i)
n *= degree+1;
- unit_support_points.resize(n);
+ this->unit_support_points.resize(n);
const double step = 1./degree;
Point<dim> p;
if (dim>2)
p(2) = iz * step;
- unit_support_points[k++] = p;
+ this->unit_support_points[k++] = p;
};
};
FE_DGQ<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_value(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_value(i, p);
}
FE_DGQ<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_grad(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad(i, p);
}
FE_DGQ<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
return polynomial_space.compute_grad_grad(i, p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad_grad(i, p);
}
// Rotate xy-plane
// counter-clockwise
case 'z':
- for (unsigned int iz=0;iz<((dim>2) ? n:1);++iz)
- for (unsigned int j=0;j<n;++j)
- for (unsigned int i=0;i<n;++i)
- {
- unsigned int k = n*i-j+n-1 + n*n*iz;
- numbers[l++] = k;
- }
- break;
- // Rotate xy-plane
- // clockwise
+ for (unsigned int iz=0;iz<((dim>2) ? n:1);++iz)
+ for (unsigned int j=0;j<n;++j)
+ for (unsigned int i=0;i<n;++i)
+ {
+ unsigned int k = n*i-j+n-1 + n*n*iz;
+ numbers[l++] = k;
+ }
+ break;
+ // Rotate xy-plane
+ // clockwise
case 'Z':
- for (unsigned int iz=0;iz<((dim>2) ? n:1);++iz)
- for (unsigned int iy=0;iy<n;++iy)
- for (unsigned int ix=0;ix<n;++ix)
- {
- unsigned int k = n*ix-iy+n-1 + n*n*iz;
- numbers[k] = l++;
- }
- break;
- // Rotate yz-plane
- // counter-clockwise
+ for (unsigned int iz=0;iz<((dim>2) ? n:1);++iz)
+ for (unsigned int iy=0;iy<n;++iy)
+ for (unsigned int ix=0;ix<n;++ix)
+ {
+ unsigned int k = n*ix-iy+n-1 + n*n*iz;
+ numbers[k] = l++;
+ }
+ break;
+ // Rotate yz-plane
+ // counter-clockwise
case 'x':
- Assert (dim>2,
- typename FiniteElementData<dim>::ExcSpaceDimensionMismatch (dim,3));
- for (unsigned int iz=0;iz<n;++iz)
- for (unsigned int iy=0;iy<n;++iy)
- for (unsigned int ix=0;ix<n;++ix)
- {
- unsigned int k = n*(n*iy-iz+n-1) + ix;
- numbers[l++] = k;
- }
- break;
- // Rotate yz-plane
- // clockwise
+ Assert (dim>2,
+ typename FiniteElementData<dim>::ExcSpaceDimensionMismatch (dim,3));
+ for (unsigned int iz=0;iz<n;++iz)
+ for (unsigned int iy=0;iy<n;++iy)
+ for (unsigned int ix=0;ix<n;++ix)
+ {
+ unsigned int k = n*(n*iy-iz+n-1) + ix;
+ numbers[l++] = k;
+ }
+ break;
+ // Rotate yz-plane
+ // clockwise
case 'X':
- Assert (dim>2,
- typename FiniteElementData<dim>::ExcSpaceDimensionMismatch (dim,3));
- for (unsigned int iz=0;iz<n;++iz)
- for (unsigned int iy=0;iy<n;++iy)
- for (unsigned int ix=0;ix<n;++ix)
- {
- unsigned int k = n*(n*iy-iz+n-1) + ix;
- numbers[k] = l++;
- }
- break;
+ Assert (dim>2,
+ typename FiniteElementData<dim>::ExcSpaceDimensionMismatch (dim,3));
+ for (unsigned int iz=0;iz<n;++iz)
+ for (unsigned int iy=0;iy<n;++iy)
+ for (unsigned int ix=0;ix<n;++ix)
+ {
+ unsigned int k = n*(n*iy-iz+n-1) + ix;
+ numbers[k] = l++;
+ }
+ break;
default:
Assert (false, ExcNotImplemented ());
}
// allocate memory
if (flags & update_values)
{
- values.resize (dofs_per_cell);
- data->shape_values.resize(dofs_per_cell,
+ values.resize (this->dofs_per_cell);
+ data->shape_values.resize(this->dofs_per_cell,
std::vector<double>(n_q_points));
}
if (flags & update_gradients)
{
- grads.resize (dofs_per_cell);
- data->shape_gradients.resize(dofs_per_cell,
+ grads.resize (this->dofs_per_cell);
+ data->shape_gradients.resize(this->dofs_per_cell,
std::vector<Tensor<1,dim> >(n_q_points));
}
{
polynomial_space.compute(quadrature.point(i),
values, grads, grad_grads);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
if (flags & update_values)
data->shape_values[k][i] = values[k];
const UpdateFlags flags(fe_data.current_update_flags());
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
if (flags & update_values)
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
FE_DGQ<dim>::has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const
{
- Assert (shape_index < dofs_per_cell,
- ExcIndexRange (shape_index, 0, dofs_per_cell));
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
Assert (face_index < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
switch (dim)
{
case 1:
- {
- // in 1d, things are simple. since
- // there is only one degree of
- // freedom per vertex in this
- // class, the first is on vertex 0
- // (==face 0 in some sense), the
- // second on face 1:
- return (((shape_index == 0) && (face_index == 0)) ||
- ((shape_index == 1) && (face_index == 1)));
- };
+ {
+ // in 1d, things are simple. since
+ // there is only one degree of
+ // freedom per vertex in this
+ // class, the first is on vertex 0
+ // (==face 0 in some sense), the
+ // second on face 1:
+ return (((shape_index == 0) && (face_index == 0)) ||
+ ((shape_index == 1) && (face_index == 1)));
+ };
case 2:
- {
- if (face_index==0 && shape_index < n)
- return true;
- if (face_index==1 && (shape_index % n) == degree)
- return true;
- if (face_index==2 && shape_index >= dofs_per_cell-n)
- return true;
- if (face_index==3 && (shape_index % n) == 0)
- return true;
- return false;
- };
+ {
+ if (face_index==0 && shape_index < n)
+ return true;
+ if (face_index==1 && (shape_index % n) == degree)
+ return true;
+ if (face_index==2 && shape_index >= this->dofs_per_cell-n)
+ return true;
+ if (face_index==3 && (shape_index % n) == 0)
+ return true;
+ return false;
+ };
case 3:
- {
- const unsigned int in2 = shape_index % n2;
+ {
+ const unsigned int in2 = shape_index % n2;
- // y=0
- if (face_index==0 && in2 < n )
- return true;
- // y=1
- if (face_index==1 && in2 >= n2-n)
- return true;
- // z=0
- if (face_index==2 && shape_index < n2)
- return true;
- // x=1
- if (face_index==3 && (shape_index % n) == n-1)
- return true;
- // z=1
- if (face_index==4 && shape_index >= dofs_per_cell - n2)
- return true;
- // x=0
- if (face_index==5 && (shape_index % n) == 0)
- return true;
- return false;
- };
+ // y=0
+ if (face_index==0 && in2 < n )
+ return true;
+ // y=1
+ if (face_index==1 && in2 >= n2-n)
+ return true;
+ // z=0
+ if (face_index==2 && shape_index < n2)
+ return true;
+ // x=1
+ if (face_index==3 && (shape_index % n) == n-1)
+ return true;
+ // z=1
+ if (face_index==4 && shape_index >= this->dofs_per_cell - n2)
+ return true;
+ // x=0
+ if (face_index==5 && (shape_index % n) == 0)
+ return true;
+ return false;
+ };
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
return true;
}
std::vector<bool> (dim,false),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,
std::vector<bool>(dim,true))),
- degree(degree)
+ degree(degree)
{
Assert (dim >= 2, ExcNotUsefulInThisDimension());
// are defined. otherwise set them
// to invalid size
if (degree<Matrices::n_constraint_matrices+1)
- interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+ this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
else
- interface_constraints.reinit(0,0);
+ this->interface_constraints.reinit(0,0);
// 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<GeometryInfo<dim>::children_per_cell; ++c)
- prolongation[c].fill (Matrices::embedding[degree-1][c]);
+ this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
else
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- prolongation[i].reinit(0,0);
+ this->prolongation[i].reinit(0,0);
// then fill restriction
// matrices. they are hardcoded for
switch (dim)
{
case 2: // 2d
- {
- switch (degree)
- {
- case 1:
- {
- // DoF on bottom line
- // of coarse cell will
- // be mean value of
- // bottom DoFs on the
- // two adjacent child
- // cells
- restriction[0](0,0) = 0.5;
- restriction[1](0,0) = 0.5;
- // same for other DoFs
- restriction[1](1,1) = 0.5;
- restriction[2](1,1) = 0.5;
-
- restriction[2](2,2) = 0.5;
- restriction[3](2,2) = 0.5;
-
- restriction[3](3,3) = 0.5;
- restriction[0](3,3) = 0.5;
-
- break;
- };
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // DoF on bottom line
+ // of coarse cell will
+ // be mean value of
+ // bottom DoFs on the
+ // two adjacent child
+ // cells
+ this->restriction[0](0,0) = 0.5;
+ this->restriction[1](0,0) = 0.5;
+ // same for other DoFs
+ this->restriction[1](1,1) = 0.5;
+ this->restriction[2](1,1) = 0.5;
+
+ this->restriction[2](2,2) = 0.5;
+ this->restriction[3](2,2) = 0.5;
+
+ this->restriction[3](3,3) = 0.5;
+ this->restriction[0](3,3) = 0.5;
+
+ break;
+ };
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- restriction[i].reinit(0,0);
- };
- };
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ };
- break;
- };
+ break;
+ };
case 3: // 3d
- {
- switch (degree)
- {
- case 1:
- {
- // same principle as in
- // 2d
- restriction[0](0,0) = 0.5;
- restriction[1](0,0) = 0.5;
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // same principle as in
+ // 2d
+ this->restriction[0](0,0) = 0.5;
+ this->restriction[1](0,0) = 0.5;
- restriction[1](1,1) = 0.5;
- restriction[2](1,1) = 0.5;
+ this->restriction[1](1,1) = 0.5;
+ this->restriction[2](1,1) = 0.5;
- restriction[2](2,2) = 0.5;
- restriction[3](2,2) = 0.5;
+ this->restriction[2](2,2) = 0.5;
+ this->restriction[3](2,2) = 0.5;
- restriction[3](3,3) = 0.5;
- restriction[0](3,3) = 0.5;
+ this->restriction[3](3,3) = 0.5;
+ this->restriction[0](3,3) = 0.5;
- restriction[4](4,4) = 0.5;
- restriction[5](4,4) = 0.5;
+ this->restriction[4](4,4) = 0.5;
+ this->restriction[5](4,4) = 0.5;
- restriction[5](5,5) = 0.5;
- restriction[6](5,5) = 0.5;
+ this->restriction[5](5,5) = 0.5;
+ this->restriction[6](5,5) = 0.5;
- restriction[6](6,6) = 0.5;
- restriction[7](6,6) = 0.5;
+ this->restriction[6](6,6) = 0.5;
+ this->restriction[7](6,6) = 0.5;
- restriction[7](7,7) = 0.5;
- restriction[4](7,7) = 0.5;
+ this->restriction[7](7,7) = 0.5;
+ this->restriction[4](7,7) = 0.5;
- restriction[1](8,8) = 0.5;
- restriction[5](8,8) = 0.5;
+ this->restriction[1](8,8) = 0.5;
+ this->restriction[5](8,8) = 0.5;
- restriction[2](9,9) = 0.5;
- restriction[6](9,9) = 0.5;
+ this->restriction[2](9,9) = 0.5;
+ this->restriction[6](9,9) = 0.5;
- restriction[3](10,10) = 0.5;
- restriction[7](10,10) = 0.5;
+ this->restriction[3](10,10) = 0.5;
+ this->restriction[7](10,10) = 0.5;
- restriction[0](11,11) = 0.5;
- restriction[5](11,11) = 0.5;
+ this->restriction[0](11,11) = 0.5;
+ this->restriction[5](11,11) = 0.5;
- break;
- };
+ break;
+ };
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- restriction[i].reinit(0,0);
- };
- };
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ };
- break;
- };
+ break;
+ };
default:
- Assert (false,ExcNotImplemented());
+ Assert (false,ExcNotImplemented());
}
// finally fill in support points
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component < dim, ExcIndexRange (component, 0, dim));
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements
- case 1:
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ switch (i)
{
- switch (i)
- {
- // (1-y, 0)
- case 0: return (component == 0 ? 1-p(1) : 0);
- // (0,x)
- case 1: return (component == 0 ? 0 : p(0));
- // (y, 0)
- case 2: return (component == 0 ? p(1) : 0);
- // (0, 1-x)
- case 3: return (component == 0 ? 0 : 1-p(0));
-
- // there are
- // only four
- // shape
- // functions!?
- default:
- Assert (false, ExcInternalError());
- return 0;
- };
+ // (1-y, 0)
+ case 0: return (component == 0 ? 1-p(1) : 0);
+ // (0,x)
+ case 1: return (component == 0 ? 0 : p(0));
+ // (y, 0)
+ case 2: return (component == 0 ? p(1) : 0);
+ // (0, 1-x)
+ case 3: return (component == 0 ? 0 : 1-p(0));
+
+ // there are
+ // only four
+ // shape
+ // functions!?
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
};
+ };
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
case 3: // 3D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements
- case 1:
- {
- // note that the
- // degrees of freedom
- // on opposite faces
- // have a common vector
- // direction, so simply
- // that a little. these
- // directions are:
- //
- // for lines 0, 2, 4, 6:
- // (1,0,0)
- // for lines 1, 3, 5, 7:
- // (0,0,1)
- // for lines 8, 9, 10, 11:
- // (0,1,0)
- //
- // thus, sort out all
- // those cases where
- // the component is
- // zero anyway, and
- // only otherwise
- // compute the
- // spatially dependent
- // part which is then
- // also the return
- // value
- if (((i<8) && (((i%2==0) && (component!=0)) ||
- ((i%2==1) && (component!=2)))) ||
- ((i>=8) && (component != 1)))
- return 0;
-
- // now we know that the
- // only non-zero
- // component is
- // requested:
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ // note that the
+ // degrees of freedom
+ // on opposite faces
+ // have a common vector
+ // direction, so simply
+ // that a little. these
+ // directions are:
+ //
+ // for lines 0, 2, 4, 6:
+ // (1,0,0)
+ // for lines 1, 3, 5, 7:
+ // (0,0,1)
+ // for lines 8, 9, 10, 11:
+ // (0,1,0)
+ //
+ // thus, sort out all
+ // those cases where
+ // the component is
+ // zero anyway, and
+ // only otherwise
+ // compute the
+ // spatially dependent
+ // part which is then
+ // also the return
+ // value
+ if (((i<8) && (((i%2==0) && (component!=0)) ||
+ ((i%2==1) && (component!=2)))) ||
+ ((i>=8) && (component != 1)))
+ return 0;
+
+ // now we know that the
+ // only non-zero
+ // component is
+ // requested:
//[Anna]:checked: OK
- const double x = p(0),
- y = p(1),
- z = p(2);
- switch (i)
- {
- case 0: return (1-y)*(1-z);
- case 2: return (1-y)*z;
- case 1: return x*(1-y);
- case 3: return (1-x)*(1-y);
-
- case 4: return y*(1-z);
- case 6: return y*z;
- case 5: return x*y;
- case 7: return (1-x)*y;
+ const double x = p(0),
+ y = p(1),
+ z = p(2);
+ switch (i)
+ {
+ case 0: return (1-y)*(1-z);
+ case 2: return (1-y)*z;
+ case 1: return x*(1-y);
+ case 3: return (1-x)*(1-y);
+
+ case 4: return y*(1-z);
+ case 6: return y*z;
+ case 5: return x*y;
+ case 7: return (1-x)*y;
- case 8: return (1-x)*(1-z);
- case 9: return x*(1-z);
- case 10: return x*z;
- case 11: return (1-x)*z;
- default:
- Assert (false, ExcInternalError());
- return 0;
- };
+ case 8: return (1-x)*(1-z);
+ case 9: return x*(1-z);
+ case 10: return x*z;
+ case 11: return (1-x)*z;
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
};
+ };
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return 0;
const Point<dim> &/*p*/,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component < dim, ExcIndexRange (component, 0, dim));
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // 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, thrid
- // index=component
- // within gradient
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // 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, thrid
+ // index=component
+ // within gradient
//[Anna]: checked:OK
- static const double unit_gradients[4][2][2]
- = { { {0.,-1.}, {0.,0.} },
- { {0.,0.}, {1.,0.} },
- { {0.,+1.}, {0.,0.} },
- { {0.,0.}, {-1.,0.} } };
- return Tensor<1,dim>(unit_gradients[i][component]);
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ static const double unit_gradients[4][2][2]
+ = { { {0.,-1.}, {0.,0.} },
+ { {0.,0.}, {1.,0.} },
+ { {0.,+1.}, {0.,0.} },
+ { {0.,0.}, {-1.,0.} } };
+ return Tensor<1,dim>(unit_gradients[i][component]);
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return Tensor<1,dim>();
const Point<dim> &/*p*/,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component < dim, ExcIndexRange (component, 0, dim));
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements. their second
- // derivatives on the
- // unit cell are zero
- case 1:
- {
- return Tensor<2,dim>();
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
case 3: // 3D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements. their second
- // derivatives on the
- // unit cell are zero
- case 1:
- {
- return Tensor<2,dim>();
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return Tensor<2,dim>();
// all degrees of freedom are on
// edges, and their order is the
// same as the edges themselves
- unit_support_points.resize(GeometryInfo<dim>::lines_per_cell * degree);
+ this->unit_support_points.resize(GeometryInfo<dim>::lines_per_cell * degree);
unsigned int index = 0;
for (unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
{
// between the vertices of each
// line
for (unsigned int d=0; d<degree; ++d, ++index)
- unit_support_points[index]
+ this->unit_support_points[index]
= (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1);
};
};
// is this correct? all DoFs on lines, none on faces or bubbles?
// do this the same as above, but
// for one dimension less
- unit_face_support_points.resize(GeometryInfo<dim-1>::lines_per_cell * degree);
+ this->unit_face_support_points.resize(GeometryInfo<dim-1>::lines_per_cell * degree);
unsigned int index = 0;
for (unsigned int line=0; line<GeometryInfo<dim-1>::lines_per_cell; ++line)
{
// between the vertices of each
// line
for (unsigned int d=0; d<degree; ++d, ++index)
- unit_face_support_points[index]
+ this->unit_face_support_points[index]
= (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1);
};
};
FE_Nedelec<dim>::has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const
{
- Assert (shape_index < dofs_per_cell,
- ExcIndexRange (shape_index, 0, dofs_per_cell));
+ Assert (shape_index < this->dofs_per_cell,
+ ExcIndexRange (shape_index, 0, this->dofs_per_cell));
Assert (face_index < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
std::vector<bool> (1,false),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
- degree(degree),
- renumber(dofs_per_cell, 0),
- renumber_inverse(dofs_per_cell, 0),
- face_renumber(dofs_per_face, 0),
- polynomial_space(LagrangeEquidistant::generate_complete_basis(degree))
+ degree(degree),
+ renumber(this->dofs_per_cell, 0),
+ renumber_inverse(this->dofs_per_cell, 0),
+ face_renumber(this->dofs_per_face, 0),
+ polynomial_space(LagrangeEquidistant::generate_complete_basis(degree))
{
// do some internal book-keeping on
// cells and faces. if in 1d, the
// build inverse of renumbering
// vector
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
renumber_inverse[renumber[i]]=i;
// copy constraint matrices if they
// are defined. otherwise set them
// to invalid size
if ((dim>1) && (degree<Matrices::n_constraint_matrices+1))
- interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+ this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
else
- interface_constraints.reinit(0,0);
+ this->interface_constraints.reinit(0,0);
// 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<GeometryInfo<dim>::children_per_cell; ++c)
- prolongation[c].fill (Matrices::embedding[degree-1][c]);
+ this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
else
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- prolongation[i].reinit(0,0);
+ this->prolongation[i].reinit(0,0);
// then fill restriction
// matrices. they are hardcoded for
switch (dim)
{
case 1: // 1d
- {
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](2,1) = 1;
- restriction[1](1,1) = 1;
- restriction[1](1,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](2,3) = 1;
- restriction[1](1,1) = 1;
- restriction[1](3,2) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](2,3) = 1;
- restriction[0](3,1) = 1;
- restriction[1](1,1) = 1;
- restriction[1](3,0) = 1;
- restriction[1](4,3) = 1;
- break;
+ {
+ switch (degree)
+ {
+ case 1:
+ this->restriction[0](0,0) = 1;
+ this->restriction[1](1,1) = 1;
+ break;
+ case 2:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](2,1) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](1,1) = 1;
+ break;
+ case 3:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](2,3) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](3,2) = 1;
+ break;
+ case 4:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](2,3) = 1;
+ this->restriction[0](3,1) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](3,0) = 1;
+ this->restriction[1](4,3) = 1;
+ break;
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- restriction[i].reinit(0,0);
- };
- }
- break;
- };
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
case 2: // 2d
- {
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- restriction[2](2,2) = 1;
- restriction[3](3,3) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](4,1) = 1;
- restriction[0](7,3) = 1;
- restriction[0](8,2) = 1;
- restriction[1](1,1) = 1;
- restriction[1](4,0) = 1;
- restriction[1](5,2) = 1;
- restriction[1](8,3) = 1;
- restriction[2](2,2) = 1;
- restriction[2](5,1) = 1;
- restriction[2](6,3) = 1;
- restriction[2](8,0) = 1;
- restriction[3](3,3) = 1;
- restriction[3](6,2) = 1;
- restriction[3](7,0) = 1;
- restriction[3](8,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](4,5) = 1;
- restriction[0](10,11) = 1;
- restriction[0](12,15) = 1;
- restriction[1](1,1) = 1;
- restriction[1](5,4) = 1;
- restriction[1](6,7) = 1;
- restriction[1](13,14) = 1;
- restriction[2](2,2) = 1;
- restriction[2](7,6) = 1;
- restriction[2](9,8) = 1;
- restriction[2](15,12) = 1;
- restriction[3](3,3) = 1;
- restriction[3](8,9) = 1;
- restriction[3](11,10) = 1;
- restriction[3](14,13) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](4,5) = 1;
- restriction[0](5,1) = 1;
- restriction[0](13,14) = 1;
- restriction[0](14,3) = 1;
- restriction[0](16,20) = 1;
- restriction[0](17,8) = 1;
- restriction[0](19,11) = 1;
- restriction[0](20,2) = 1;
- restriction[1](1,1) = 1;
- restriction[1](5,0) = 1;
- restriction[1](6,5) = 1;
- restriction[1](7,8) = 1;
- restriction[1](8,2) = 1;
- restriction[1](17,14) = 1;
- restriction[1](18,20) = 1;
- restriction[1](20,3) = 1;
- restriction[1](21,11) = 1;
- restriction[2](2,2) = 1;
- restriction[2](8,1) = 1;
- restriction[2](9,8) = 1;
- restriction[2](11,3) = 1;
- restriction[2](12,11) = 1;
- restriction[2](20,0) = 1;
- restriction[2](21,5) = 1;
- restriction[2](23,14) = 1;
- restriction[2](24,20) = 1;
- restriction[3](3,3) = 1;
- restriction[3](10,11) = 1;
- restriction[3](11,2) = 1;
- restriction[3](14,0) = 1;
- restriction[3](15,14) = 1;
- restriction[3](19,5) = 1;
- restriction[3](20,1) = 1;
- restriction[3](22,20) = 1;
- restriction[3](23,8) = 1;
- break;
-
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- restriction[i].reinit(0,0);
- };
- }
- break;
- };
+ {
+ switch (degree)
+ {
+ case 1:
+ this->restriction[0](0,0) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[3](3,3) = 1;
+ break;
+ case 2:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](4,1) = 1;
+ this->restriction[0](7,3) = 1;
+ this->restriction[0](8,2) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](4,0) = 1;
+ this->restriction[1](5,2) = 1;
+ this->restriction[1](8,3) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](5,1) = 1;
+ this->restriction[2](6,3) = 1;
+ this->restriction[2](8,0) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](6,2) = 1;
+ this->restriction[3](7,0) = 1;
+ this->restriction[3](8,1) = 1;
+ break;
+ case 3:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](4,5) = 1;
+ this->restriction[0](10,11) = 1;
+ this->restriction[0](12,15) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](5,4) = 1;
+ this->restriction[1](6,7) = 1;
+ this->restriction[1](13,14) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](7,6) = 1;
+ this->restriction[2](9,8) = 1;
+ this->restriction[2](15,12) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](8,9) = 1;
+ this->restriction[3](11,10) = 1;
+ this->restriction[3](14,13) = 1;
+ break;
+ case 4:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](4,5) = 1;
+ this->restriction[0](5,1) = 1;
+ this->restriction[0](13,14) = 1;
+ this->restriction[0](14,3) = 1;
+ this->restriction[0](16,20) = 1;
+ this->restriction[0](17,8) = 1;
+ this->restriction[0](19,11) = 1;
+ this->restriction[0](20,2) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](5,0) = 1;
+ this->restriction[1](6,5) = 1;
+ this->restriction[1](7,8) = 1;
+ this->restriction[1](8,2) = 1;
+ this->restriction[1](17,14) = 1;
+ this->restriction[1](18,20) = 1;
+ this->restriction[1](20,3) = 1;
+ this->restriction[1](21,11) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](8,1) = 1;
+ this->restriction[2](9,8) = 1;
+ this->restriction[2](11,3) = 1;
+ this->restriction[2](12,11) = 1;
+ this->restriction[2](20,0) = 1;
+ this->restriction[2](21,5) = 1;
+ this->restriction[2](23,14) = 1;
+ this->restriction[2](24,20) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](10,11) = 1;
+ this->restriction[3](11,2) = 1;
+ this->restriction[3](14,0) = 1;
+ this->restriction[3](15,14) = 1;
+ this->restriction[3](19,5) = 1;
+ this->restriction[3](20,1) = 1;
+ this->restriction[3](22,20) = 1;
+ this->restriction[3](23,8) = 1;
+ break;
+
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
case 3: // 3d
- {
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- restriction[2](2,2) = 1;
- restriction[3](3,3) = 1;
- restriction[4](4,4) = 1;
- restriction[5](5,5) = 1;
- restriction[6](6,6) = 1;
- restriction[7](7,7) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](8,1) = 1;
- restriction[0](11,3) = 1;
- restriction[0](16,4) = 1;
- restriction[0](20,2) = 1;
- restriction[0](22,5) = 1;
- restriction[0](25,7) = 1;
- restriction[0](26,6) = 1;
- restriction[1](1,1) = 1;
- restriction[1](8,0) = 1;
- restriction[1](9,2) = 1;
- restriction[1](17,5) = 1;
- restriction[1](20,3) = 1;
- restriction[1](22,4) = 1;
- restriction[1](23,6) = 1;
- restriction[1](26,7) = 1;
- restriction[2](2,2) = 1;
- restriction[2](9,1) = 1;
- restriction[2](10,3) = 1;
- restriction[2](18,6) = 1;
- restriction[2](20,0) = 1;
- restriction[2](23,5) = 1;
- restriction[2](24,7) = 1;
- restriction[2](26,4) = 1;
- restriction[3](3,3) = 1;
- restriction[3](10,2) = 1;
- restriction[3](11,0) = 1;
- restriction[3](19,7) = 1;
- restriction[3](20,1) = 1;
- restriction[3](24,6) = 1;
- restriction[3](25,4) = 1;
- restriction[3](26,5) = 1;
- restriction[4](4,4) = 1;
- restriction[4](12,5) = 1;
- restriction[4](15,7) = 1;
- restriction[4](16,0) = 1;
- restriction[4](21,6) = 1;
- restriction[4](22,1) = 1;
- restriction[4](25,3) = 1;
- restriction[4](26,2) = 1;
- restriction[5](5,5) = 1;
- restriction[5](12,4) = 1;
- restriction[5](13,6) = 1;
- restriction[5](17,1) = 1;
- restriction[5](21,7) = 1;
- restriction[5](22,0) = 1;
- restriction[5](23,2) = 1;
- restriction[5](26,3) = 1;
- restriction[6](6,6) = 1;
- restriction[6](13,5) = 1;
- restriction[6](14,7) = 1;
- restriction[6](18,2) = 1;
- restriction[6](21,4) = 1;
- restriction[6](23,1) = 1;
- restriction[6](24,3) = 1;
- restriction[6](26,0) = 1;
- restriction[7](7,7) = 1;
- restriction[7](14,6) = 1;
- restriction[7](15,4) = 1;
- restriction[7](19,3) = 1;
- restriction[7](21,5) = 1;
- restriction[7](24,2) = 1;
- restriction[7](25,0) = 1;
- restriction[7](26,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](8,9) = 1;
- restriction[0](14,15) = 1;
- restriction[0](24,25) = 1;
- restriction[0](32,35) = 1;
- restriction[0](40,43) = 1;
- restriction[0](52,55) = 1;
- restriction[0](56,63) = 1;
- restriction[1](1,1) = 1;
- restriction[1](9,8) = 1;
- restriction[1](10,11) = 1;
- restriction[1](26,27) = 1;
- restriction[1](33,34) = 1;
- restriction[1](41,42) = 1;
- restriction[1](44,47) = 1;
- restriction[1](57,62) = 1;
- restriction[2](2,2) = 1;
- restriction[2](11,10) = 1;
- restriction[2](13,12) = 1;
- restriction[2](28,29) = 1;
- restriction[2](35,32) = 1;
- restriction[2](46,45) = 1;
- restriction[2](49,50) = 1;
- restriction[2](61,58) = 1;
- restriction[3](3,3) = 1;
- restriction[3](12,13) = 1;
- restriction[3](15,14) = 1;
- restriction[3](30,31) = 1;
- restriction[3](34,33) = 1;
- restriction[3](48,51) = 1;
- restriction[3](54,53) = 1;
- restriction[3](60,59) = 1;
- restriction[4](4,4) = 1;
- restriction[4](16,17) = 1;
- restriction[4](22,23) = 1;
- restriction[4](25,24) = 1;
- restriction[4](36,39) = 1;
- restriction[4](42,41) = 1;
- restriction[4](53,54) = 1;
- restriction[4](58,61) = 1;
- restriction[5](5,5) = 1;
- restriction[5](17,16) = 1;
- restriction[5](18,19) = 1;
- restriction[5](27,26) = 1;
- restriction[5](37,38) = 1;
- restriction[5](43,40) = 1;
- restriction[5](45,46) = 1;
- restriction[5](59,60) = 1;
- restriction[6](6,6) = 1;
- restriction[6](19,18) = 1;
- restriction[6](21,20) = 1;
- restriction[6](29,28) = 1;
- restriction[6](39,36) = 1;
- restriction[6](47,44) = 1;
- restriction[6](51,48) = 1;
- restriction[6](63,56) = 1;
- restriction[7](7,7) = 1;
- restriction[7](20,21) = 1;
- restriction[7](23,22) = 1;
- restriction[7](31,30) = 1;
- restriction[7](38,37) = 1;
- restriction[7](50,49) = 1;
- restriction[7](55,52) = 1;
- restriction[7](62,57) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](8,9) = 1;
- restriction[0](9,1) = 1;
- restriction[0](17,18) = 1;
- restriction[0](18,3) = 1;
- restriction[0](32,33) = 1;
- restriction[0](33,4) = 1;
- restriction[0](44,48) = 1;
- restriction[0](45,12) = 1;
- restriction[0](47,15) = 1;
- restriction[0](48,2) = 1;
- restriction[0](62,66) = 1;
- restriction[0](63,36) = 1;
- restriction[0](65,21) = 1;
- restriction[0](66,5) = 1;
- restriction[0](89,93) = 1;
- restriction[0](90,30) = 1;
- restriction[0](92,42) = 1;
- restriction[0](93,7) = 1;
- restriction[0](98,111) = 1;
- restriction[0](99,75) = 1;
- restriction[0](101,57) = 1;
- restriction[0](102,24) = 1;
- restriction[0](107,84) = 1;
- restriction[0](108,39) = 1;
- restriction[0](110,27) = 1;
- restriction[0](111,6) = 1;
- restriction[1](1,1) = 1;
- restriction[1](9,0) = 1;
- restriction[1](10,9) = 1;
- restriction[1](11,12) = 1;
- restriction[1](12,2) = 1;
- restriction[1](35,36) = 1;
- restriction[1](36,5) = 1;
- restriction[1](45,18) = 1;
- restriction[1](46,48) = 1;
- restriction[1](48,3) = 1;
- restriction[1](49,15) = 1;
- restriction[1](63,33) = 1;
- restriction[1](64,66) = 1;
- restriction[1](66,4) = 1;
- restriction[1](67,21) = 1;
- restriction[1](71,75) = 1;
- restriction[1](72,24) = 1;
- restriction[1](74,39) = 1;
- restriction[1](75,6) = 1;
- restriction[1](99,93) = 1;
- restriction[1](100,111) = 1;
- restriction[1](102,30) = 1;
- restriction[1](103,57) = 1;
- restriction[1](108,42) = 1;
- restriction[1](109,84) = 1;
- restriction[1](111,7) = 1;
- restriction[1](112,27) = 1;
- restriction[2](2,2) = 1;
- restriction[2](12,1) = 1;
- restriction[2](13,12) = 1;
- restriction[2](15,3) = 1;
- restriction[2](16,15) = 1;
- restriction[2](38,39) = 1;
- restriction[2](39,6) = 1;
- restriction[2](48,0) = 1;
- restriction[2](49,9) = 1;
- restriction[2](51,18) = 1;
- restriction[2](52,48) = 1;
- restriction[2](74,36) = 1;
- restriction[2](75,5) = 1;
- restriction[2](77,75) = 1;
- restriction[2](78,24) = 1;
- restriction[2](81,42) = 1;
- restriction[2](82,84) = 1;
- restriction[2](84,7) = 1;
- restriction[2](85,27) = 1;
- restriction[2](108,33) = 1;
- restriction[2](109,66) = 1;
- restriction[2](111,4) = 1;
- restriction[2](112,21) = 1;
- restriction[2](117,93) = 1;
- restriction[2](118,111) = 1;
- restriction[2](120,30) = 1;
- restriction[2](121,57) = 1;
- restriction[3](3,3) = 1;
- restriction[3](14,15) = 1;
- restriction[3](15,2) = 1;
- restriction[3](18,0) = 1;
- restriction[3](19,18) = 1;
- restriction[3](41,42) = 1;
- restriction[3](42,7) = 1;
- restriction[3](47,9) = 1;
- restriction[3](48,1) = 1;
- restriction[3](50,48) = 1;
- restriction[3](51,12) = 1;
- restriction[3](80,84) = 1;
- restriction[3](81,39) = 1;
- restriction[3](83,27) = 1;
- restriction[3](84,6) = 1;
- restriction[3](92,33) = 1;
- restriction[3](93,4) = 1;
- restriction[3](95,93) = 1;
- restriction[3](96,30) = 1;
- restriction[3](107,66) = 1;
- restriction[3](108,36) = 1;
- restriction[3](110,21) = 1;
- restriction[3](111,5) = 1;
- restriction[3](116,111) = 1;
- restriction[3](117,75) = 1;
- restriction[3](119,57) = 1;
- restriction[3](120,24) = 1;
- restriction[4](4,4) = 1;
- restriction[4](20,21) = 1;
- restriction[4](21,5) = 1;
- restriction[4](29,30) = 1;
- restriction[4](30,7) = 1;
- restriction[4](33,0) = 1;
- restriction[4](34,33) = 1;
- restriction[4](53,57) = 1;
- restriction[4](54,24) = 1;
- restriction[4](56,27) = 1;
- restriction[4](57,6) = 1;
- restriction[4](65,9) = 1;
- restriction[4](66,1) = 1;
- restriction[4](68,66) = 1;
- restriction[4](69,36) = 1;
- restriction[4](90,18) = 1;
- restriction[4](91,93) = 1;
- restriction[4](93,3) = 1;
- restriction[4](94,42) = 1;
- restriction[4](101,48) = 1;
- restriction[4](102,12) = 1;
- restriction[4](104,111) = 1;
- restriction[4](105,75) = 1;
- restriction[4](110,15) = 1;
- restriction[4](111,2) = 1;
- restriction[4](113,84) = 1;
- restriction[4](114,39) = 1;
- restriction[5](5,5) = 1;
- restriction[5](21,4) = 1;
- restriction[5](22,21) = 1;
- restriction[5](23,24) = 1;
- restriction[5](24,6) = 1;
- restriction[5](36,1) = 1;
- restriction[5](37,36) = 1;
- restriction[5](54,30) = 1;
- restriction[5](55,57) = 1;
- restriction[5](57,7) = 1;
- restriction[5](58,27) = 1;
- restriction[5](66,0) = 1;
- restriction[5](67,9) = 1;
- restriction[5](69,33) = 1;
- restriction[5](70,66) = 1;
- restriction[5](72,12) = 1;
- restriction[5](73,75) = 1;
- restriction[5](75,2) = 1;
- restriction[5](76,39) = 1;
- restriction[5](102,18) = 1;
- restriction[5](103,48) = 1;
- restriction[5](105,93) = 1;
- restriction[5](106,111) = 1;
- restriction[5](111,3) = 1;
- restriction[5](112,15) = 1;
- restriction[5](114,42) = 1;
- restriction[5](115,84) = 1;
- restriction[6](6,6) = 1;
- restriction[6](24,5) = 1;
- restriction[6](25,24) = 1;
- restriction[6](27,7) = 1;
- restriction[6](28,27) = 1;
- restriction[6](39,2) = 1;
- restriction[6](40,39) = 1;
- restriction[6](57,4) = 1;
- restriction[6](58,21) = 1;
- restriction[6](60,30) = 1;
- restriction[6](61,57) = 1;
- restriction[6](75,1) = 1;
- restriction[6](76,36) = 1;
- restriction[6](78,12) = 1;
- restriction[6](79,75) = 1;
- restriction[6](84,3) = 1;
- restriction[6](85,15) = 1;
- restriction[6](87,42) = 1;
- restriction[6](88,84) = 1;
- restriction[6](111,0) = 1;
- restriction[6](112,9) = 1;
- restriction[6](114,33) = 1;
- restriction[6](115,66) = 1;
- restriction[6](120,18) = 1;
- restriction[6](121,48) = 1;
- restriction[6](123,93) = 1;
- restriction[6](124,111) = 1;
- restriction[7](7,7) = 1;
- restriction[7](26,27) = 1;
- restriction[7](27,6) = 1;
- restriction[7](30,4) = 1;
- restriction[7](31,30) = 1;
- restriction[7](42,3) = 1;
- restriction[7](43,42) = 1;
- restriction[7](56,21) = 1;
- restriction[7](57,5) = 1;
- restriction[7](59,57) = 1;
- restriction[7](60,24) = 1;
- restriction[7](83,15) = 1;
- restriction[7](84,2) = 1;
- restriction[7](86,84) = 1;
- restriction[7](87,39) = 1;
- restriction[7](93,0) = 1;
- restriction[7](94,33) = 1;
- restriction[7](96,18) = 1;
- restriction[7](97,93) = 1;
- restriction[7](110,9) = 1;
- restriction[7](111,1) = 1;
- restriction[7](113,66) = 1;
- restriction[7](114,36) = 1;
- restriction[7](119,48) = 1;
- restriction[7](120,12) = 1;
- restriction[7](122,111) = 1;
- restriction[7](123,75) = 1;
- break;
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- restriction[i].reinit(0,0);
- };
- }
- break;
- };
+ {
+ switch (degree)
+ {
+ case 1:
+ this->restriction[0](0,0) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[4](4,4) = 1;
+ this->restriction[5](5,5) = 1;
+ this->restriction[6](6,6) = 1;
+ this->restriction[7](7,7) = 1;
+ break;
+ case 2:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](8,1) = 1;
+ this->restriction[0](11,3) = 1;
+ this->restriction[0](16,4) = 1;
+ this->restriction[0](20,2) = 1;
+ this->restriction[0](22,5) = 1;
+ this->restriction[0](25,7) = 1;
+ this->restriction[0](26,6) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](8,0) = 1;
+ this->restriction[1](9,2) = 1;
+ this->restriction[1](17,5) = 1;
+ this->restriction[1](20,3) = 1;
+ this->restriction[1](22,4) = 1;
+ this->restriction[1](23,6) = 1;
+ this->restriction[1](26,7) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](9,1) = 1;
+ this->restriction[2](10,3) = 1;
+ this->restriction[2](18,6) = 1;
+ this->restriction[2](20,0) = 1;
+ this->restriction[2](23,5) = 1;
+ this->restriction[2](24,7) = 1;
+ this->restriction[2](26,4) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](10,2) = 1;
+ this->restriction[3](11,0) = 1;
+ this->restriction[3](19,7) = 1;
+ this->restriction[3](20,1) = 1;
+ this->restriction[3](24,6) = 1;
+ this->restriction[3](25,4) = 1;
+ this->restriction[3](26,5) = 1;
+ this->restriction[4](4,4) = 1;
+ this->restriction[4](12,5) = 1;
+ this->restriction[4](15,7) = 1;
+ this->restriction[4](16,0) = 1;
+ this->restriction[4](21,6) = 1;
+ this->restriction[4](22,1) = 1;
+ this->restriction[4](25,3) = 1;
+ this->restriction[4](26,2) = 1;
+ this->restriction[5](5,5) = 1;
+ this->restriction[5](12,4) = 1;
+ this->restriction[5](13,6) = 1;
+ this->restriction[5](17,1) = 1;
+ this->restriction[5](21,7) = 1;
+ this->restriction[5](22,0) = 1;
+ this->restriction[5](23,2) = 1;
+ this->restriction[5](26,3) = 1;
+ this->restriction[6](6,6) = 1;
+ this->restriction[6](13,5) = 1;
+ this->restriction[6](14,7) = 1;
+ this->restriction[6](18,2) = 1;
+ this->restriction[6](21,4) = 1;
+ this->restriction[6](23,1) = 1;
+ this->restriction[6](24,3) = 1;
+ this->restriction[6](26,0) = 1;
+ this->restriction[7](7,7) = 1;
+ this->restriction[7](14,6) = 1;
+ this->restriction[7](15,4) = 1;
+ this->restriction[7](19,3) = 1;
+ this->restriction[7](21,5) = 1;
+ this->restriction[7](24,2) = 1;
+ this->restriction[7](25,0) = 1;
+ this->restriction[7](26,1) = 1;
+ break;
+ case 3:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](8,9) = 1;
+ this->restriction[0](14,15) = 1;
+ this->restriction[0](24,25) = 1;
+ this->restriction[0](32,35) = 1;
+ this->restriction[0](40,43) = 1;
+ this->restriction[0](52,55) = 1;
+ this->restriction[0](56,63) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](9,8) = 1;
+ this->restriction[1](10,11) = 1;
+ this->restriction[1](26,27) = 1;
+ this->restriction[1](33,34) = 1;
+ this->restriction[1](41,42) = 1;
+ this->restriction[1](44,47) = 1;
+ this->restriction[1](57,62) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](11,10) = 1;
+ this->restriction[2](13,12) = 1;
+ this->restriction[2](28,29) = 1;
+ this->restriction[2](35,32) = 1;
+ this->restriction[2](46,45) = 1;
+ this->restriction[2](49,50) = 1;
+ this->restriction[2](61,58) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](12,13) = 1;
+ this->restriction[3](15,14) = 1;
+ this->restriction[3](30,31) = 1;
+ this->restriction[3](34,33) = 1;
+ this->restriction[3](48,51) = 1;
+ this->restriction[3](54,53) = 1;
+ this->restriction[3](60,59) = 1;
+ this->restriction[4](4,4) = 1;
+ this->restriction[4](16,17) = 1;
+ this->restriction[4](22,23) = 1;
+ this->restriction[4](25,24) = 1;
+ this->restriction[4](36,39) = 1;
+ this->restriction[4](42,41) = 1;
+ this->restriction[4](53,54) = 1;
+ this->restriction[4](58,61) = 1;
+ this->restriction[5](5,5) = 1;
+ this->restriction[5](17,16) = 1;
+ this->restriction[5](18,19) = 1;
+ this->restriction[5](27,26) = 1;
+ this->restriction[5](37,38) = 1;
+ this->restriction[5](43,40) = 1;
+ this->restriction[5](45,46) = 1;
+ this->restriction[5](59,60) = 1;
+ this->restriction[6](6,6) = 1;
+ this->restriction[6](19,18) = 1;
+ this->restriction[6](21,20) = 1;
+ this->restriction[6](29,28) = 1;
+ this->restriction[6](39,36) = 1;
+ this->restriction[6](47,44) = 1;
+ this->restriction[6](51,48) = 1;
+ this->restriction[6](63,56) = 1;
+ this->restriction[7](7,7) = 1;
+ this->restriction[7](20,21) = 1;
+ this->restriction[7](23,22) = 1;
+ this->restriction[7](31,30) = 1;
+ this->restriction[7](38,37) = 1;
+ this->restriction[7](50,49) = 1;
+ this->restriction[7](55,52) = 1;
+ this->restriction[7](62,57) = 1;
+ break;
+ case 4:
+ this->restriction[0](0,0) = 1;
+ this->restriction[0](8,9) = 1;
+ this->restriction[0](9,1) = 1;
+ this->restriction[0](17,18) = 1;
+ this->restriction[0](18,3) = 1;
+ this->restriction[0](32,33) = 1;
+ this->restriction[0](33,4) = 1;
+ this->restriction[0](44,48) = 1;
+ this->restriction[0](45,12) = 1;
+ this->restriction[0](47,15) = 1;
+ this->restriction[0](48,2) = 1;
+ this->restriction[0](62,66) = 1;
+ this->restriction[0](63,36) = 1;
+ this->restriction[0](65,21) = 1;
+ this->restriction[0](66,5) = 1;
+ this->restriction[0](89,93) = 1;
+ this->restriction[0](90,30) = 1;
+ this->restriction[0](92,42) = 1;
+ this->restriction[0](93,7) = 1;
+ this->restriction[0](98,111) = 1;
+ this->restriction[0](99,75) = 1;
+ this->restriction[0](101,57) = 1;
+ this->restriction[0](102,24) = 1;
+ this->restriction[0](107,84) = 1;
+ this->restriction[0](108,39) = 1;
+ this->restriction[0](110,27) = 1;
+ this->restriction[0](111,6) = 1;
+ this->restriction[1](1,1) = 1;
+ this->restriction[1](9,0) = 1;
+ this->restriction[1](10,9) = 1;
+ this->restriction[1](11,12) = 1;
+ this->restriction[1](12,2) = 1;
+ this->restriction[1](35,36) = 1;
+ this->restriction[1](36,5) = 1;
+ this->restriction[1](45,18) = 1;
+ this->restriction[1](46,48) = 1;
+ this->restriction[1](48,3) = 1;
+ this->restriction[1](49,15) = 1;
+ this->restriction[1](63,33) = 1;
+ this->restriction[1](64,66) = 1;
+ this->restriction[1](66,4) = 1;
+ this->restriction[1](67,21) = 1;
+ this->restriction[1](71,75) = 1;
+ this->restriction[1](72,24) = 1;
+ this->restriction[1](74,39) = 1;
+ this->restriction[1](75,6) = 1;
+ this->restriction[1](99,93) = 1;
+ this->restriction[1](100,111) = 1;
+ this->restriction[1](102,30) = 1;
+ this->restriction[1](103,57) = 1;
+ this->restriction[1](108,42) = 1;
+ this->restriction[1](109,84) = 1;
+ this->restriction[1](111,7) = 1;
+ this->restriction[1](112,27) = 1;
+ this->restriction[2](2,2) = 1;
+ this->restriction[2](12,1) = 1;
+ this->restriction[2](13,12) = 1;
+ this->restriction[2](15,3) = 1;
+ this->restriction[2](16,15) = 1;
+ this->restriction[2](38,39) = 1;
+ this->restriction[2](39,6) = 1;
+ this->restriction[2](48,0) = 1;
+ this->restriction[2](49,9) = 1;
+ this->restriction[2](51,18) = 1;
+ this->restriction[2](52,48) = 1;
+ this->restriction[2](74,36) = 1;
+ this->restriction[2](75,5) = 1;
+ this->restriction[2](77,75) = 1;
+ this->restriction[2](78,24) = 1;
+ this->restriction[2](81,42) = 1;
+ this->restriction[2](82,84) = 1;
+ this->restriction[2](84,7) = 1;
+ this->restriction[2](85,27) = 1;
+ this->restriction[2](108,33) = 1;
+ this->restriction[2](109,66) = 1;
+ this->restriction[2](111,4) = 1;
+ this->restriction[2](112,21) = 1;
+ this->restriction[2](117,93) = 1;
+ this->restriction[2](118,111) = 1;
+ this->restriction[2](120,30) = 1;
+ this->restriction[2](121,57) = 1;
+ this->restriction[3](3,3) = 1;
+ this->restriction[3](14,15) = 1;
+ this->restriction[3](15,2) = 1;
+ this->restriction[3](18,0) = 1;
+ this->restriction[3](19,18) = 1;
+ this->restriction[3](41,42) = 1;
+ this->restriction[3](42,7) = 1;
+ this->restriction[3](47,9) = 1;
+ this->restriction[3](48,1) = 1;
+ this->restriction[3](50,48) = 1;
+ this->restriction[3](51,12) = 1;
+ this->restriction[3](80,84) = 1;
+ this->restriction[3](81,39) = 1;
+ this->restriction[3](83,27) = 1;
+ this->restriction[3](84,6) = 1;
+ this->restriction[3](92,33) = 1;
+ this->restriction[3](93,4) = 1;
+ this->restriction[3](95,93) = 1;
+ this->restriction[3](96,30) = 1;
+ this->restriction[3](107,66) = 1;
+ this->restriction[3](108,36) = 1;
+ this->restriction[3](110,21) = 1;
+ this->restriction[3](111,5) = 1;
+ this->restriction[3](116,111) = 1;
+ this->restriction[3](117,75) = 1;
+ this->restriction[3](119,57) = 1;
+ this->restriction[3](120,24) = 1;
+ this->restriction[4](4,4) = 1;
+ this->restriction[4](20,21) = 1;
+ this->restriction[4](21,5) = 1;
+ this->restriction[4](29,30) = 1;
+ this->restriction[4](30,7) = 1;
+ this->restriction[4](33,0) = 1;
+ this->restriction[4](34,33) = 1;
+ this->restriction[4](53,57) = 1;
+ this->restriction[4](54,24) = 1;
+ this->restriction[4](56,27) = 1;
+ this->restriction[4](57,6) = 1;
+ this->restriction[4](65,9) = 1;
+ this->restriction[4](66,1) = 1;
+ this->restriction[4](68,66) = 1;
+ this->restriction[4](69,36) = 1;
+ this->restriction[4](90,18) = 1;
+ this->restriction[4](91,93) = 1;
+ this->restriction[4](93,3) = 1;
+ this->restriction[4](94,42) = 1;
+ this->restriction[4](101,48) = 1;
+ this->restriction[4](102,12) = 1;
+ this->restriction[4](104,111) = 1;
+ this->restriction[4](105,75) = 1;
+ this->restriction[4](110,15) = 1;
+ this->restriction[4](111,2) = 1;
+ this->restriction[4](113,84) = 1;
+ this->restriction[4](114,39) = 1;
+ this->restriction[5](5,5) = 1;
+ this->restriction[5](21,4) = 1;
+ this->restriction[5](22,21) = 1;
+ this->restriction[5](23,24) = 1;
+ this->restriction[5](24,6) = 1;
+ this->restriction[5](36,1) = 1;
+ this->restriction[5](37,36) = 1;
+ this->restriction[5](54,30) = 1;
+ this->restriction[5](55,57) = 1;
+ this->restriction[5](57,7) = 1;
+ this->restriction[5](58,27) = 1;
+ this->restriction[5](66,0) = 1;
+ this->restriction[5](67,9) = 1;
+ this->restriction[5](69,33) = 1;
+ this->restriction[5](70,66) = 1;
+ this->restriction[5](72,12) = 1;
+ this->restriction[5](73,75) = 1;
+ this->restriction[5](75,2) = 1;
+ this->restriction[5](76,39) = 1;
+ this->restriction[5](102,18) = 1;
+ this->restriction[5](103,48) = 1;
+ this->restriction[5](105,93) = 1;
+ this->restriction[5](106,111) = 1;
+ this->restriction[5](111,3) = 1;
+ this->restriction[5](112,15) = 1;
+ this->restriction[5](114,42) = 1;
+ this->restriction[5](115,84) = 1;
+ this->restriction[6](6,6) = 1;
+ this->restriction[6](24,5) = 1;
+ this->restriction[6](25,24) = 1;
+ this->restriction[6](27,7) = 1;
+ this->restriction[6](28,27) = 1;
+ this->restriction[6](39,2) = 1;
+ this->restriction[6](40,39) = 1;
+ this->restriction[6](57,4) = 1;
+ this->restriction[6](58,21) = 1;
+ this->restriction[6](60,30) = 1;
+ this->restriction[6](61,57) = 1;
+ this->restriction[6](75,1) = 1;
+ this->restriction[6](76,36) = 1;
+ this->restriction[6](78,12) = 1;
+ this->restriction[6](79,75) = 1;
+ this->restriction[6](84,3) = 1;
+ this->restriction[6](85,15) = 1;
+ this->restriction[6](87,42) = 1;
+ this->restriction[6](88,84) = 1;
+ this->restriction[6](111,0) = 1;
+ this->restriction[6](112,9) = 1;
+ this->restriction[6](114,33) = 1;
+ this->restriction[6](115,66) = 1;
+ this->restriction[6](120,18) = 1;
+ this->restriction[6](121,48) = 1;
+ this->restriction[6](123,93) = 1;
+ this->restriction[6](124,111) = 1;
+ this->restriction[7](7,7) = 1;
+ this->restriction[7](26,27) = 1;
+ this->restriction[7](27,6) = 1;
+ this->restriction[7](30,4) = 1;
+ this->restriction[7](31,30) = 1;
+ this->restriction[7](42,3) = 1;
+ this->restriction[7](43,42) = 1;
+ this->restriction[7](56,21) = 1;
+ this->restriction[7](57,5) = 1;
+ this->restriction[7](59,57) = 1;
+ this->restriction[7](60,24) = 1;
+ this->restriction[7](83,15) = 1;
+ this->restriction[7](84,2) = 1;
+ this->restriction[7](86,84) = 1;
+ this->restriction[7](87,39) = 1;
+ this->restriction[7](93,0) = 1;
+ this->restriction[7](94,33) = 1;
+ this->restriction[7](96,18) = 1;
+ this->restriction[7](97,93) = 1;
+ this->restriction[7](110,9) = 1;
+ this->restriction[7](111,1) = 1;
+ this->restriction[7](113,66) = 1;
+ this->restriction[7](114,36) = 1;
+ this->restriction[7](119,48) = 1;
+ this->restriction[7](120,12) = 1;
+ this->restriction[7](122,111) = 1;
+ this->restriction[7](123,75) = 1;
+ break;
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
default:
- Assert (false,ExcNotImplemented());
+ Assert (false,ExcNotImplemented());
}
// finally fill in support points
FE_Q<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
return polynomial_space.compute_value(renumber_inverse[i], p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_value(renumber_inverse[i], p);
}
FE_Q<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
return polynomial_space.compute_grad(renumber_inverse[i], p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad(renumber_inverse[i], p);
}
FE_Q<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
}
for (unsigned int i=1; i<dim; ++i)
n *= degree+1;
- unit_support_points.resize(n);
+ this->unit_support_points.resize(n);
const double step = 1./degree;
Point<dim> p;
if (dim>2)
p(2) = iz * step;
- unit_support_points[renumber[k++]] = p;
+ this->unit_support_points[renumber[k++]] = p;
};
};
for (unsigned int i=1; i<codim; ++i)
n *= degree+1;
- unit_face_support_points.resize(n);
+ this->unit_face_support_points.resize(n);
const double step = 1./degree;
Point<codim> p;
if (codim>2)
p(2) = iz * step;
- unit_face_support_points[face_renumber[k++]] = p;
+ this->unit_face_support_points[face_renumber[k++]] = p;
};
};
switch (dim)
{
case 1:
- {
- const unsigned int values[GeometryInfo<1>::vertices_per_cell]
- = { 0, degree };
- index = values[i];
- break;
- };
+ {
+ const unsigned int values[GeometryInfo<1>::vertices_per_cell]
+ = { 0, degree };
+ index = values[i];
+ break;
+ };
case 2:
- {
- const unsigned int values[GeometryInfo<2>::vertices_per_cell]
- = { 0, degree, n*degree+degree, n*degree };
- index = values[i];
- break;
- };
+ {
+ const unsigned int values[GeometryInfo<2>::vertices_per_cell]
+ = { 0, degree, n*degree+degree, n*degree };
+ index = values[i];
+ break;
+ };
case 3:
- {
- const unsigned int values[GeometryInfo<3>::vertices_per_cell]
- = { 0, degree,
- n*n*degree + degree, n*n*degree,
- n*degree, n*degree+degree,
- n*n*degree + n*degree+degree, n*n*degree + n*degree};
- index = values[i];
- break;
- };
+ {
+ const unsigned int values[GeometryInfo<3>::vertices_per_cell]
+ = { 0, degree,
+ n*n*degree + degree, n*n*degree,
+ n*degree, n*degree+degree,
+ n*n*degree + n*degree+degree, n*n*degree + n*degree};
+ index = values[i];
+ break;
+ };
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
renumber[index] = i;
case 100:
case 200: case 202:
case 300: case 302: case 304: case 306:
- incr = 1;
- break;
- // lines in y-direction
+ incr = 1;
+ break;
+ // lines in y-direction
case 201: case 203:
case 308: case 309: case 310: case 311:
- incr = n;
- break;
- // lines in z-direction
+ incr = n;
+ break;
+ // lines in z-direction
case 301: case 303: case 305: case 307:
- incr = n*n;
- break;
+ incr = n*n;
+ break;
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
switch (i+100*dim)
{
case 100:
case 200: case 203:
case 300: case 303: case 308:
- tensorstart = 0;
- break;
- // x=1 y=z=0
+ tensorstart = 0;
+ break;
+ // x=1 y=z=0
case 201:
case 301: case 309:
- tensorstart = degree;
- break;
- // y=1 x=z=0
+ tensorstart = degree;
+ break;
+ // y=1 x=z=0
case 202:
case 304: case 307:
- tensorstart = n*degree;
- break;
- // x=z=1 y=0
+ tensorstart = n*degree;
+ break;
+ // x=z=1 y=0
case 310:
- tensorstart = n*n*degree+degree;
- break;
- // z=1 x=y=0
+ tensorstart = n*n*degree+degree;
+ break;
+ // z=1 x=y=0
case 302: case 311:
- tensorstart = n*n*degree;
- break;
- // x=y=1 z=0
+ tensorstart = n*n*degree;
+ break;
+ // x=y=1 z=0
case 305:
- tensorstart = n*degree+degree;
- break;
- // y=z=1 x=0
+ tensorstart = n*degree+degree;
+ break;
+ // y=z=1 x=0
case 306:
- tensorstart = n*n*n-n;
- break;
+ tensorstart = n*n*n-n;
+ break;
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
for (unsigned int jx = 1; jx<degree ;++jx)
switch (i)
{
case 0:
- tensorstart = 0; incx = 1;
- if (dim==2)
- incy = n;
- else
- incy = n*n;
- break;
+ tensorstart = 0; incx = 1;
+ if (dim==2)
+ incy = n;
+ else
+ incy = n*n;
+ break;
case 1:
- tensorstart = n*degree; incx = 1; incy = n*n;
- break;
+ tensorstart = n*degree; incx = 1; incy = n*n;
+ break;
case 2:
- tensorstart = 0; incx = 1; incy = n;
- break;
+ tensorstart = 0; incx = 1; incy = n;
+ break;
case 3:
- tensorstart = degree; incx = n; incy = n*n;
- break;
+ tensorstart = degree; incx = n; incy = n*n;
+ break;
case 4:
- tensorstart = n*n*degree; incx = 1; incy = n;
- break;
+ tensorstart = n*n*degree; incx = 1; incy = n;
+ break;
case 5:
- tensorstart = 0; incx = n; incy = n*n;
- break;
+ tensorstart = 0; incx = n; incy = n*n;
+ break;
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
for (unsigned int jy = 1; jy<degree; jy++)
// allocate memory
if (flags & update_values)
{
- values.resize (dofs_per_cell);
- data->shape_values.resize(dofs_per_cell,
+ values.resize (this->dofs_per_cell);
+ data->shape_values.resize(this->dofs_per_cell,
std::vector<double>(n_q_points));
}
if (flags & update_gradients)
{
- grads.resize (dofs_per_cell);
- data->shape_gradients.resize(dofs_per_cell,
+ grads.resize (this->dofs_per_cell);
+ data->shape_gradients.resize(this->dofs_per_cell,
std::vector<Tensor<1,dim> >(n_q_points));
}
values, grads, grad_grads);
if (flags & update_values)
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
data->shape_values[renumber[k]][i] = values[k];
if (flags & update_gradients)
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
data->shape_gradients[renumber[k]][i] = grads[k];
}
return data;
const UpdateFlags flags(fe_data.current_update_flags());
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
if (flags & update_values)
for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
+ for (unsigned int k=0; k<this->dofs_per_cell; ++k)
{
for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
if (flags & update_values)
FE_Q<dim>::has_support_on_face (const unsigned int shape_index_,
const unsigned int face_index) const
{
- Assert (shape_index_ < dofs_per_cell,
- ExcIndexRange (shape_index_, 0, dofs_per_cell));
+ Assert (shape_index_ < this->dofs_per_cell,
+ ExcIndexRange (shape_index_, 0, this->dofs_per_cell));
Assert (face_index < GeometryInfo<dim>::faces_per_cell,
ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
Assert (dim<=3, ExcNotImplemented());
const unsigned int cell_start = (dim==2)
- ? first_quad_index
- : first_hex_index;
+ ? this->first_quad_index
+ : this->first_hex_index;
const unsigned int face_start = (dim==2)
- ? first_line_index
- : first_quad_index;
+ ? this->first_line_index
+ : this->first_quad_index;
// Interior degrees of
// freedom correspond to
// equal to the face index.
if (shape_index >= face_start)
{
- shape_index -= first_line_index;
- shape_index /= dofs_per_face;
+ shape_index -= this->first_line_index;
+ shape_index /= this->dofs_per_face;
return (shape_index == face_index);
}
// Only degrees of freedom on
// a vertex are left.
- shape_index /= dofs_per_vertex;
+ shape_index /= this->dofs_per_vertex;
// Use a table to figure out
// which face is neighbor to
// which vertex.
case 314: case 324: case 354: case 315:
case 325: case 335: case 316: case 336:
case 346: case 317: case 347: case 357:
- return true;
+ return true;
default:
- return false;
+ return false;
}
return true;
};
FiniteElement<dim> (multiply_dof_numbers(fe, n_elements),
compute_restriction_is_additive_flags (fe, n_elements),
compute_nonzero_components(fe, n_elements)),
- base_elements(1)
+ base_elements(1)
{
base_elements[0] = ElementPair(fe.clone(), n_elements);
base_elements[0].first->subscribe ();
fe2, n2),
compute_nonzero_components(fe1, n1,
fe2, n2)),
- base_elements(2)
+ base_elements(2)
{
base_elements[0] = ElementPair(fe1.clone(), n1);
base_elements[0].first->subscribe ();
compute_nonzero_components(fe1, n1,
fe2, n2,
fe3, n3)),
- base_elements(3)
+ base_elements(3)
{
base_elements[0] = ElementPair(fe1.clone(), n1);
base_elements[0].first->subscribe ();
switch (n_base_elements())
{
case 1:
- return new FESystem(base_element(0),
- element_multiplicity(0));
+ return new FESystem(base_element(0),
+ element_multiplicity(0));
case 2:
- return new FESystem(base_element(0),
- element_multiplicity(0),
- base_element(1),
- element_multiplicity(1));
+ return new FESystem(base_element(0),
+ element_multiplicity(0),
+ base_element(1),
+ element_multiplicity(1));
case 3:
- return new FESystem(base_element(0),
- element_multiplicity(0),
- base_element(1),
- element_multiplicity(1),
- base_element(2),
- element_multiplicity(2));
+ return new FESystem(base_element(0),
+ element_multiplicity(0),
+ base_element(1),
+ element_multiplicity(1),
+ base_element(2),
+ element_multiplicity(2));
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
return 0;
}
FESystem<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (is_primitive(i),
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (this->is_primitive(i),
typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
- return (base_element(system_to_base_table[i].first.first)
- .shape_value(system_to_base_table[i].second, p));
+ return (base_element(this->system_to_base_table[i].first.first)
+ .shape_value(this->system_to_base_table[i].second, p));
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (component < n_components(),
- ExcIndexRange (component, 0, n_components()));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (component < this->n_components(),
+ ExcIndexRange (component, 0, this->n_components()));
// first find out to which of the
// base elements this desired
// component belongs, and which
// component within this base
// element it is
- const unsigned int base = component_to_base(i).first;
- const unsigned int component_in_base = component_to_base(i).second;
+ const unsigned int base = this->component_to_base(i).first;
+ const unsigned int component_in_base = this->component_to_base(i).second;
// then get value from base
// element. note that that will
// primitive; thus, there is no
// need to check this here
return (base_element(base).
- shape_value_component(system_to_base_table[i].second,
+ shape_value_component(this->system_to_base_table[i].second,
p,
component_in_base));
}
FESystem<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (is_primitive(i),
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (this->is_primitive(i),
typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
- return (base_element(system_to_base_table[i].first.first)
- .shape_grad(system_to_base_table[i].second, p));
+ return (base_element(this->system_to_base_table[i].first.first)
+ .shape_grad(this->system_to_base_table[i].second, p));
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (component < n_components(),
- ExcIndexRange (component, 0, n_components()));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (component < this->n_components(),
+ ExcIndexRange (component, 0, this->n_components()));
// first find out to which of the
// base elements this desired
// component belongs, and which
// component within this base
// element it is
- const unsigned int base = component_to_base(i).first;
- const unsigned int component_in_base = component_to_base(i).second;
+ const unsigned int base = this->component_to_base(i).first;
+ const unsigned int component_in_base = this->component_to_base(i).second;
// then get value from base
// element. note that that will
// primitive; thus, there is no
// need to check this here
return (base_element(base).
- shape_grad_component(system_to_base_table[i].second,
+ shape_grad_component(this->system_to_base_table[i].second,
p,
component_in_base));
}
FESystem<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (is_primitive(i),
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (this->is_primitive(i),
typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
- return (base_element(system_to_base_table[i].first.first)
- .shape_grad_grad(system_to_base_table[i].second, p));
+ return (base_element(this->system_to_base_table[i].first.first)
+ .shape_grad_grad(this->system_to_base_table[i].second, p));
}
const Point<dim> &p,
const unsigned int component) const
{
- Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
- Assert (component < n_components(),
- ExcIndexRange (component, 0, n_components()));
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell));
+ Assert (component < this->n_components(),
+ ExcIndexRange (component, 0, this->n_components()));
// first find out to which of the
// base elements this desired
// component belongs, and which
// component within this base
// element it is
- const unsigned int base = component_to_base(i).first;
- const unsigned int component_in_base = component_to_base(i).second;
+ const unsigned int base = this->component_to_base(i).first;
+ const unsigned int component_in_base = this->component_to_base(i).second;
// then get value from base
// element. note that that will
// primitive; thus, there is no
// need to check this here
return (base_element(base).
- shape_grad_grad_component(system_to_base_table[i].second,
+ shape_grad_grad_component(this->system_to_base_table[i].second,
p,
component_in_base));
}
for (unsigned int point=0; point<quadrature.n_quadrature_points; ++point)
for (unsigned int k=0; k<base_fe.dofs_per_cell; ++k)
{
- const unsigned int system_index=component_to_system_index(comp,k);
+ const unsigned int system_index=this->component_to_system_index(comp,k);
if (base_flags & update_values)
data.shape_values(system_index, point)=
base_data.shape_values(k,point);
for (unsigned int base=0; base < n_base_elements(); ++base)
for (unsigned int m = 0; m < element_multiplicity(base); ++m)
for (unsigned int k=0; k<base_element(base).n_components(); ++k)
- component_to_base_table[total_index++] = std::make_pair(base,k);
- Assert (total_index == component_to_base_table.size(),
+ this->component_to_base_table[total_index++] = std::make_pair(base,k);
+ Assert (total_index == this->component_to_base_table.size(),
ExcInternalError());
// Initialize index tables.
= (base_element(base).dofs_per_vertex*vertex_number +
local_index);
- system_to_base_table[total_index]
+ this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base, m), index_in_base);
if (base_element(base).is_primitive(index_in_base))
= comp_start + comp_in_base;
const unsigned int index_in_comp
= base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
+ this->system_to_component_table[total_index]
= std::make_pair (comp, index_in_comp);
}
else
- system_to_component_table[total_index] = non_primitive_index;
+ this->system_to_component_table[total_index] = non_primitive_index;
}
}
local_index +
base_element(base).first_line_index);
- system_to_base_table[total_index]
+ this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
if (base_element(base).is_primitive(index_in_base))
= base_element(base).system_to_component_index(index_in_base).first;
const unsigned int comp
= comp_start + comp_in_base;
- const unsigned int index_in_comp
- = base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ this->system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
}
else
- system_to_component_table[total_index] = non_primitive_index;
+ this->system_to_component_table[total_index] = non_primitive_index;
}
}
local_index +
base_element(base).first_quad_index);
- system_to_base_table[total_index]
+ this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
if (base_element(base).is_primitive(index_in_base))
= base_element(base).system_to_component_index(index_in_base).first;
const unsigned int comp
= comp_start + comp_in_base;
- const unsigned int index_in_comp
- = base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ this->system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
}
else
- system_to_component_table[total_index] = non_primitive_index;
+ this->system_to_component_table[total_index] = non_primitive_index;
}
}
local_index +
base_element(base).first_hex_index);
- system_to_base_table[total_index]
+ this->system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), index_in_base);
if (base_element(base).is_primitive(index_in_base))
= base_element(base).system_to_component_index(index_in_base).first;
const unsigned int comp
= comp_start + comp_in_base;
- const unsigned int index_in_comp
- = base_element(base).system_to_component_index(index_in_base).second;
- system_to_component_table[total_index]
- = std::make_pair (comp, index_in_comp);
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ this->system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
}
else
- system_to_component_table[total_index] = non_primitive_index;
+ this->system_to_component_table[total_index] = non_primitive_index;
}
}
// components to linear
// index. Fortunately, this is the
// inverse of what we just did.
- std::vector<unsigned int> dofs_per_component (n_components(), 0);
- for (unsigned int sys=0; sys<dofs_per_cell; ++sys)
- ++dofs_per_component[system_to_component_index(sys).first];
- for (unsigned int component=0; component<n_components(); ++component)
- component_to_system_table[component].resize(dofs_per_component[component]);
+ std::vector<unsigned int> dofs_per_component (this->n_components(), 0);
+ for (unsigned int sys=0; sys<this->dofs_per_cell; ++sys)
+ ++dofs_per_component[this->system_to_component_index(sys).first];
+ for (unsigned int component=0; component<this->n_components(); ++component)
+ this->component_to_system_table[component].resize(dofs_per_component[component]);
// then go the reverse way to fill the array
- for (unsigned int sys=0; sys<dofs_per_cell; ++sys)
+ for (unsigned int sys=0; sys<this->dofs_per_cell; ++sys)
{
const unsigned int
- comp = system_to_component_index(sys).first,
- index_in_comp = system_to_component_index(sys).second;
+ comp = this->system_to_component_index(sys).first,
+ index_in_comp = this->system_to_component_index(sys).second;
- Assert (comp < component_to_system_table.size(),
+ Assert (comp < this->component_to_system_table.size(),
ExcInternalError());
- Assert (index_in_comp < component_to_system_table[comp].size(),
+ Assert (index_in_comp < this->component_to_system_table[comp].size(),
ExcInternalError());
- component_to_system_table[comp][index_in_comp] = sys;
+ this->component_to_system_table[comp][index_in_comp] = sys;
};
};
= (base_element(base).dofs_per_vertex*vertex_number +
local_index);
- face_system_to_base_table[total_index]
+ this->face_system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), face_index_in_base);
if (base_element(base).is_primitive(index_in_base))
= comp_start + comp_in_base;
const unsigned int face_index_in_comp
= base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
+ this->face_system_to_component_table[total_index]
= std::make_pair (comp, face_index_in_comp);
}
else
- face_system_to_component_table[total_index] = non_primitive_index;
+ this->face_system_to_component_table[total_index] = non_primitive_index;
}
}
base_element(base).dofs_per_line * line_number +
local_index);
- face_system_to_base_table[total_index]
+ this->face_system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), face_index_in_base);
if (base_element(base).is_primitive(index_in_base))
= comp_start + comp_in_base;
const unsigned int face_index_in_comp
= base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
+ this->face_system_to_component_table[total_index]
= std::make_pair (comp, face_index_in_comp);
}
else
- face_system_to_component_table[total_index] = non_primitive_index;
+ this->face_system_to_component_table[total_index] = non_primitive_index;
}
}
base_element(base).dofs_per_quad * quad_number +
local_index);
- face_system_to_base_table[total_index]
+ this->face_system_to_base_table[total_index]
= std::make_pair (std::make_pair(base,m), face_index_in_base);
if (base_element(base).is_primitive(index_in_base))
= comp_start + comp_in_base;
const unsigned int face_index_in_comp
= base_element(base).face_system_to_component_index(face_index_in_base).second;
- face_system_to_component_table[total_index]
+ this->face_system_to_component_table[total_index]
= std::make_pair (comp, face_index_in_comp);
}
else
- face_system_to_component_table[total_index] = non_primitive_index;
+ this->face_system_to_component_table[total_index] = non_primitive_index;
}
}
- Assert (total_index == dofs_per_face, ExcInternalError());
- Assert (total_index == face_system_to_component_table.size(),
+ Assert (total_index == this->dofs_per_face, ExcInternalError());
+ Assert (total_index == this->face_system_to_component_table.size(),
ExcInternalError());
- Assert (total_index == face_system_to_base_table.size(),
+ Assert (total_index == this->face_system_to_base_table.size(),
ExcInternalError());
// finally, initialize reverse mapping
- std::vector<unsigned int> dofs_per_component (n_components(), 0);
- for (unsigned int sys=0; sys<dofs_per_face; ++sys)
- ++dofs_per_component[face_system_to_component_index(sys).first];
- for (unsigned int component=0; component<n_components(); ++component)
- face_component_to_system_table[component].resize(dofs_per_component[component]);
+ std::vector<unsigned int> dofs_per_component (this->n_components(), 0);
+ for (unsigned int sys=0; sys<this->dofs_per_face; ++sys)
+ ++dofs_per_component[this->face_system_to_component_index(sys).first];
+ for (unsigned int component=0; component<this->n_components(); ++component)
+ this->face_component_to_system_table[component].resize(dofs_per_component[component]);
// then go the reverse way to fill
// the array
- for (unsigned int sys=0; sys<dofs_per_face; ++sys)
+ for (unsigned int sys=0; sys<this->dofs_per_face; ++sys)
{
const unsigned int
- comp = face_system_to_component_index(sys).first,
- index_in_comp = face_system_to_component_index(sys).second;
+ comp = this->face_system_to_component_index(sys).first,
+ index_in_comp = this->face_system_to_component_index(sys).second;
- Assert (comp < face_component_to_system_table.size(),
+ Assert (comp < this->face_component_to_system_table.size(),
ExcInternalError());
- Assert (index_in_comp < face_component_to_system_table[comp].size(),
+ Assert (index_in_comp < this->face_component_to_system_table[comp].size(),
ExcInternalError());
- face_component_to_system_table[comp][index_in_comp] = sys;
+ this->face_component_to_system_table[comp][index_in_comp] = sys;
};
};
// function in conjunction with the
// numbers
// first_{line,quad,...}_index
- for (unsigned int n=0; n<interface_constraints.n(); ++n)
- for (unsigned int m=0; m<interface_constraints.m(); ++m)
+ for (unsigned int n=0; n<this->interface_constraints.n(); ++n)
+ for (unsigned int m=0; m<this->interface_constraints.m(); ++m)
{
// for the pair (n,m) find
// out which base element
// base element), second is
// index within this instance
const std::pair<std::pair<unsigned int,unsigned int>, unsigned int> n_index
- = face_system_to_base_table[n];
+ = this->face_system_to_base_table[n];
// likewise for the m
// index. this is more
switch (dim)
{
case 1:
- {
- // we should never get here!
- // (in 1d, the constraints matrix
- // should be of size zero)
- Assert (false, ExcInternalError());
- break;
- };
+ {
+ // we should never get here!
+ // (in 1d, the constraints matrix
+ // should be of size zero)
+ Assert (false, ExcInternalError());
+ break;
+ };
case 2:
- {
- // the indices m=0..d_v-1 are
- // from the center vertex.
- // their order is the same
- // as for the first vertex
- // of the whole cell, so we
- // can use the
- // system_to_base_table
- // variable (using the
- // face_s_t_base_t function would
- // yield the same)
- if (m < dofs_per_vertex)
- m_index = system_to_base_table[m];
- else
- // then come the two sets of
- // line indices
+ {
+ // the indices m=0..d_v-1 are
+ // from the center vertex.
+ // their order is the same
+ // as for the first vertex
+ // of the whole cell, so we
+ // can use the
+ // system_to_base_table
+ // variable (using the
+ // face_s_t_base_t function would
+ // yield the same)
+ if (m < this->dofs_per_vertex)
+ m_index = this->system_to_base_table[m];
+ else
+ // then come the two sets of
+ // line indices
+ {
+ const unsigned int index_in_line
+ = (m-this->dofs_per_vertex) % this->dofs_per_line;
+ const unsigned int sub_line
+ = (m-this->dofs_per_vertex) / this->dofs_per_line;
+ Assert (sub_line < 2, ExcInternalError());
+
+ // from this
+ // information, try
+ // to get base
+ // element and
+ // instance of base
+ // element. we do
+ // so by
+ // constructing the
+ // corresponding
+ // face index of m
+ // in the present
+ // element, then
+ // use
+ // face_system_to_base_table
+ const unsigned int tmp1 = 2*this->dofs_per_vertex+index_in_line;
+ m_index.first = this->face_system_to_base_table[tmp1].first;
+
+ // what we are
+ // still missing is
+ // the index of m
+ // within the base
+ // elements
+ // interface_constraints
+ // table
+ //
+ // here, the second
+ // value of
+ // face_system_to_base_table
+ // can help: it
+ // denotes the face
+ // index of that
+ // shape function
+ // within the base
+ // element. since
+ // we know that it
+ // is a line dof,
+ // we can construct
+ // the rest: tmp2
+ // will denote the
+ // index of this
+ // shape function
+ // among the line
+ // shape functions:
+ Assert (this->face_system_to_base_table[tmp1].second >=
+ 2*base_element(m_index.first.first).dofs_per_vertex,
+ ExcInternalError());
+ const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
+ 2*base_element(m_index.first.first).dofs_per_vertex;
+ Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
+ ExcInternalError());
+ m_index.second = base_element(m_index.first.first).dofs_per_vertex +
+ base_element(m_index.first.first).dofs_per_line*sub_line +
+ tmp2;
+ };
+ break;
+ };
+
+ case 3:
+ {
+ // same way as above,
+ // although a little
+ // more complicated...
+
+ // the indices
+ // m=0..5*d_v-1 are
+ // from the center and
+ // the four subline
+ // vertices. their
+ // order is the same as
+ // for the first vertex
+ // of the whole cell,
+ // so we can use the
+ // simple arithmetic
+ if (m < 5*this->dofs_per_vertex)
+ m_index = this->system_to_base_table[m];
+ else
+ // then come the 12 sets of
+ // line indices
+ if (m < 5*this->dofs_per_vertex + 12*this->dofs_per_line)
{
+ // for the
+ // meaning of all
+ // this, see the
+ // 2d part
const unsigned int index_in_line
- = (m-dofs_per_vertex) % dofs_per_line;
+ = (m-5*this->dofs_per_vertex) % this->dofs_per_line;
const unsigned int sub_line
- = (m-dofs_per_vertex) / dofs_per_line;
- Assert (sub_line < 2, ExcInternalError());
-
- // from this
- // information, try
- // to get base
- // element and
- // instance of base
- // element. we do
- // so by
- // constructing the
- // corresponding
- // face index of m
- // in the present
- // element, then
- // use
- // face_system_to_base_table
- const unsigned int tmp1 = 2*dofs_per_vertex+index_in_line;
- m_index.first = face_system_to_base_table[tmp1].first;
-
- // what we are
- // still missing is
- // the index of m
- // within the base
- // elements
- // interface_constraints
- // table
- //
- // here, the second
- // value of
- // face_system_to_base_table
- // can help: it
- // denotes the face
- // index of that
- // shape function
- // within the base
- // element. since
- // we know that it
- // is a line dof,
- // we can construct
- // the rest: tmp2
- // will denote the
- // index of this
- // shape function
- // among the line
- // shape functions:
- Assert (face_system_to_base_table[tmp1].second >=
- 2*base_element(m_index.first.first).dofs_per_vertex,
+ = (m-5*this->dofs_per_vertex) / this->dofs_per_line;
+ Assert (sub_line < 12, ExcInternalError());
+
+ const unsigned int tmp1 = 4*this->dofs_per_vertex+index_in_line;
+ m_index.first = this->face_system_to_base_table[tmp1].first;
+
+ Assert (this->face_system_to_base_table[tmp1].second >=
+ 4*base_element(m_index.first.first).dofs_per_vertex,
ExcInternalError());
- const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
- 2*base_element(m_index.first.first).dofs_per_vertex;
+ const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
+ 4*base_element(m_index.first.first).dofs_per_vertex;
Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
ExcInternalError());
- m_index.second = base_element(m_index.first.first).dofs_per_vertex +
+ m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
base_element(m_index.first.first).dofs_per_line*sub_line +
tmp2;
- };
- break;
- };
-
- case 3:
- {
- // same way as above,
- // although a little
- // more complicated...
-
- // the indices
- // m=0..5*d_v-1 are
- // from the center and
- // the four subline
- // vertices. their
- // order is the same as
- // for the first vertex
- // of the whole cell,
- // so we can use the
- // simple arithmetic
- if (m < 5*dofs_per_vertex)
- m_index = system_to_base_table[m];
+ }
else
- // then come the 12 sets of
- // line indices
- if (m < 5*dofs_per_vertex + 12*dofs_per_line)
- {
- // for the
- // meaning of all
- // this, see the
- // 2d part
- const unsigned int index_in_line
- = (m-5*dofs_per_vertex) % dofs_per_line;
- const unsigned int sub_line
- = (m-5*dofs_per_vertex) / dofs_per_line;
- Assert (sub_line < 12, ExcInternalError());
-
- const unsigned int tmp1 = 4*dofs_per_vertex+index_in_line;
- m_index.first = face_system_to_base_table[tmp1].first;
-
- Assert (face_system_to_base_table[tmp1].second >=
- 4*base_element(m_index.first.first).dofs_per_vertex,
- ExcInternalError());
- const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
- 4*base_element(m_index.first.first).dofs_per_vertex;
- Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
- ExcInternalError());
- m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
- base_element(m_index.first.first).dofs_per_line*sub_line +
- tmp2;
- }
- else
- // on one of the four sub-quads
- {
- // for the
- // meaning of all
- // this, see the
- // 2d part
- const unsigned int index_in_quad
- = (m-5*dofs_per_vertex-12*dofs_per_line) % dofs_per_line;
- const unsigned int sub_quad
- = (m-5*dofs_per_vertex-12*dofs_per_line) / dofs_per_line;
- Assert (sub_quad < 4, ExcInternalError());
-
- const unsigned int tmp1 = 4*dofs_per_vertex+4*dofs_per_line+index_in_quad;
- m_index.first = face_system_to_base_table[tmp1].first;
-
- Assert (face_system_to_base_table[tmp1].second >=
- 4*base_element(m_index.first.first).dofs_per_vertex +
- 4*base_element(m_index.first.first).dofs_per_line,
- ExcInternalError());
- const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
- 4*base_element(m_index.first.first).dofs_per_vertex -
- 4*base_element(m_index.first.first).dofs_per_line;
- Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad,
- ExcInternalError());
- m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
- 12*base_element(m_index.first.first).dofs_per_line +
- base_element(m_index.first.first).dofs_per_quad*sub_quad +
- tmp2;
- };
+ // on one of the four sub-quads
+ {
+ // for the
+ // meaning of all
+ // this, see the
+ // 2d part
+ const unsigned int index_in_quad
+ = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) % this->dofs_per_line;
+ const unsigned int sub_quad
+ = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) / this->dofs_per_line;
+ Assert (sub_quad < 4, ExcInternalError());
+
+ const unsigned int tmp1 = 4*this->dofs_per_vertex+4*this->dofs_per_line+index_in_quad;
+ m_index.first = this->face_system_to_base_table[tmp1].first;
+
+ Assert (this->face_system_to_base_table[tmp1].second >=
+ 4*base_element(m_index.first.first).dofs_per_vertex +
+ 4*base_element(m_index.first.first).dofs_per_line,
+ ExcInternalError());
+ const unsigned int tmp2 = this->face_system_to_base_table[tmp1].second -
+ 4*base_element(m_index.first.first).dofs_per_vertex -
+ 4*base_element(m_index.first.first).dofs_per_line;
+ Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad,
+ ExcInternalError());
+ m_index.second = 5*base_element(m_index.first.first).dofs_per_vertex +
+ 12*base_element(m_index.first.first).dofs_per_line +
+ base_element(m_index.first.first).dofs_per_quad*sub_quad +
+ tmp2;
+ };
- break;
- };
+ break;
+ };
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
// now that we gathered all
// definitely will be no
// coupling
if (n_index.first == m_index.first)
- interface_constraints(m,n)
+ this->interface_constraints(m,n)
= (base_element(n_index.first.first).constraints()(m_index.second,
n_index.second));
};
// element as well
if (!do_restriction)
for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
- restriction[i].reinit(0,0);
+ this->restriction[i].reinit(0,0);
if (!do_prolongation)
for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
- prolongation[i].reinit(0,0);
+ this->prolongation[i].reinit(0,0);
// distribute the matrices of the
// base finite elements to the
// couple. only DoFs that belong to
// the same instance of a base
// element may couple
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
{
// first find out to which
// base element indices i and
// element, then they cannot
// couple, so go on with the
// next index
- if (system_to_base_table[i].first !=
- system_to_base_table[j].first)
+ if (this->system_to_base_table[i].first !=
+ this->system_to_base_table[j].first)
continue;
// so get the common base
// element and the indices
// therein:
const unsigned int
- base = system_to_base_table[i].first.first;
+ base = this->system_to_base_table[i].first.first;
const unsigned int
- base_index_i = system_to_base_table[i].second,
- base_index_j = system_to_base_table[j].second;
+ base_index_i = this->system_to_base_table[i].second,
+ base_index_j = this->system_to_base_table[j].second;
// if we are sure that DoFs i
// and j may couple, then
for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
{
if (do_restriction)
- restriction[child] (i,j)
+ this->restriction[child] (i,j)
= base_element(base).restrict(child)(base_index_i,base_index_j);
if (do_prolongation)
- prolongation[child] (i,j)
+ this->prolongation[child] (i,j)
= base_element(base).prolongate(child)(base_index_i,base_index_j);
};
};
FESystem<dim>::
initialize_unit_support_points ()
{
- // if one of the base elements
- // has no support points, then
- // it makes no sense to define
- // support points for the
- // composed element, so return
- // an empty array to
- // demonstrate that
- // fact
+ // if one of the base elements
+ // has no support points, then
+ // it makes no sense to define
+ // support points for the
+ // composed element, so return
+ // an empty array to
+ // demonstrate that
+ // fact
for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
if (!base_element(base_el).has_support_points())
{
- unit_support_points.resize(0);
+ this->unit_support_points.resize(0);
return;
};
// generate unit support points
// from unit support points of sub
// elements
- unit_support_points.resize(dofs_per_cell);
+ this->unit_support_points.resize(this->dofs_per_cell);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
const unsigned int
- base = system_to_base_table[i].first.first,
- base_index = system_to_base_table[i].second;
+ base = this->system_to_base_table[i].first.first,
+ base_index = this->system_to_base_table[i].second;
Assert (base<n_base_elements(), ExcInternalError());
Assert (base_index<base_element(base).unit_support_points.size(),
ExcInternalError());
- unit_support_points[i] = base_element(base).unit_support_points[base_index];
+ this->unit_support_points[i] = base_element(base).unit_support_points[base_index];
};
}
for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
if (!base_element(base_el).has_support_points())
{
- unit_face_support_points.resize(0);
+ this->unit_face_support_points.resize(0);
return;
};
// generate unit face support points
// from unit support points of sub
// elements
- unit_face_support_points.resize(dofs_per_face);
+ this->unit_face_support_points.resize(this->dofs_per_face);
- for (unsigned int i=0; i<dofs_per_face; ++i)
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
{
- const unsigned int base_i = face_system_to_base_table[i].first.first;
- const unsigned int index_in_base = face_system_to_base_table[i].second;
+ const unsigned int base_i = this->face_system_to_base_table[i].first.first;
+ const unsigned int index_in_base = this->face_system_to_base_table[i].second;
Assert (index_in_base < base_element(base_i).unit_face_support_points.size(),
ExcInternalError());
- unit_face_support_points[i]
+ this->unit_face_support_points[i]
= base_element(base_i).unit_face_support_points[index_in_base];
};
}
const unsigned int face_index) const
{
const std::pair<unsigned int, unsigned int> component
- = system_to_component_index(shape_index);
- const unsigned int base = component_to_base(component.first).first;
+ = this->system_to_component_index(shape_index);
+ const unsigned int base = this->component_to_base(component.first).first;
return base_element(base).has_support_on_face(component.second,
face_index);
}
unsigned int
FESystem<dim>::memory_consumption () const
{
- // neglect size of data stored in
- // @p{base_elements} due to some
- // problems with teh
- // compiler. should be neglectable
- // after all, considering the size
- // of the data of the subelements
+ // neglect size of data stored in
+ // @p{base_elements} due to some
+ // problems with teh
+ // compiler. should be neglectable
+ // after all, considering the size
+ // of the data of the subelements
unsigned int mem = (FiniteElement<dim>::memory_consumption () +
sizeof (base_elements));
for (unsigned int i=0; i<base_elements.size(); ++i)
const FiniteElement<dim> &fe,
const UpdateFlags flags)
{
- update_flags = flags;
+ this->update_flags = flags;
// initialize the table mapping
// from shape function number to
// the total number of non-zero
// components accumulated over all
// shape functions
- shape_function_to_row_table.resize (fe.dofs_per_cell);
+ this->shape_function_to_row_table.resize (fe.dofs_per_cell);
unsigned int row = 0;
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
{
- shape_function_to_row_table[i] = row;
+ this->shape_function_to_row_table[i] = row;
row += fe.n_nonzero_components (i);
};
// that we will need to their
// correct size
if (flags & update_values)
- shape_values.reinit(n_nonzero_shape_components, n_quadrature_points);
+ this->shape_values.reinit(n_nonzero_shape_components, n_quadrature_points);
if (flags & update_gradients)
{
- shape_gradients.resize(n_nonzero_shape_components);
+ this->shape_gradients.resize(n_nonzero_shape_components);
for (unsigned int i=0; i<n_nonzero_shape_components; ++i)
- shape_gradients[i].resize(n_quadrature_points);
+ this->shape_gradients[i].resize(n_quadrature_points);
}
if (flags & update_second_derivatives)
{
- shape_2nd_derivatives.resize(n_nonzero_shape_components);
+ this->shape_2nd_derivatives.resize(n_nonzero_shape_components);
for (unsigned int i=0; i<n_nonzero_shape_components; ++i)
- shape_2nd_derivatives[i].resize(n_quadrature_points);
+ this->shape_2nd_derivatives[i].resize(n_quadrature_points);
}
if (flags & update_q_points)
- quadrature_points.resize(n_quadrature_points);
+ this->quadrature_points.resize(n_quadrature_points);
if (flags & update_JxW_values)
- JxW_values.resize(n_quadrature_points);
+ this->JxW_values.resize(n_quadrature_points);
if (flags & update_boundary_forms)
- boundary_forms.resize(n_quadrature_points);
+ this->boundary_forms.resize(n_quadrature_points);
if (flags & update_normal_vectors)
- normal_vectors.resize(n_quadrature_points);
+ this->normal_vectors.resize(n_quadrature_points);
}
mapping_data(0),
fe_data(0)
{
- update_flags = flags;
+ this->update_flags = flags;
}
void FEValuesBase<dim>::get_function_values (const InputVector &fe_function,
typename std::vector<number> &values) const
{
- Assert (update_flags & update_values, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
Assert (fe->n_components() == 1,
ExcWrongNoOfComponents());
Assert (values.size() == n_quadrature_points,
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
values[point] += (dof_values(shape_func) *
- shape_value(shape_func, point));
+ this->shape_value(shape_func, point));
};
for (unsigned i=0;i<values.size();++i)
Assert (values[i].size() == n_components, ExcWrongNoOfComponents());
- Assert (update_flags & update_values, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_values, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
const typename std::vector<Point<dim> > &
FEValuesBase<dim>::get_quadrature_points () const
{
- Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
- return quadrature_points;
+ Assert (this->update_flags & update_q_points, ExcAccessToUninitializedField());
+ return this->quadrature_points;
};
const std::vector<double> &
FEValuesBase<dim>::get_JxW_values () const
{
- Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
- return JxW_values;
+ Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField());
+ return this->JxW_values;
}
get_function_grads (const InputVector &fe_function,
typename std::vector<Tensor<1,dim> > &gradients) const
{
- Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField());
Assert (fe->n_components() == 1,
ExcWrongNoOfComponents());
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
- Tensor<1,dim> tmp = shape_grad(shape_func,point);
+ Tensor<1,dim> tmp = this->shape_grad(shape_func,point);
tmp *= dof_values(shape_func);
gradients[point] += tmp;
};
for (unsigned i=0; i<gradients.size(); ++i)
Assert (gradients[i].size() == n_components, ExcWrongNoOfComponents());
- Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
if (fe->is_primitive (shape_func))
{
- Tensor<1,dim> tmp = shape_grad(shape_func,point);
+ Tensor<1,dim> tmp = this->shape_grad(shape_func,point);
tmp *= dof_values(shape_func);
gradients[point][fe->system_to_component_index(shape_func).first]
+= tmp;
else
for (unsigned int c=0; c<n_components; ++c)
{
- Tensor<1,dim> tmp = shape_grad_component(shape_func,point,c);
+ Tensor<1,dim> tmp = this->shape_grad_component(shape_func,point,c);
tmp *= dof_values(shape_func);
gradients[point][c] += tmp;
};
ExcWrongNoOfComponents());
Assert (second_derivatives.size() == n_quadrature_points,
ExcWrongVectorSize(second_derivatives.size(), n_quadrature_points));
- Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
- Tensor<2,dim> tmp = shape_2nd_derivative(shape_func,point);
+ Tensor<2,dim> tmp = this->shape_2nd_derivative(shape_func,point);
tmp *= dof_values(shape_func);
second_derivatives[point] += tmp;
};
for (unsigned i=0;i<second_derivs.size();++i)
Assert (second_derivs[i].size() == n_components, ExcWrongNoOfComponents());
- Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
else
for (unsigned int c=0; c<n_components; ++c)
{
- Tensor<2,dim> tmp = shape_2nd_derivative_component(shape_func,point,c);
+ Tensor<2,dim> tmp = this->shape_2nd_derivative_component(shape_func,point,c);
tmp *= dof_values(shape_func);
second_derivs[point][c] += tmp;
};
const Point<dim> &
FEValuesBase<dim>::quadrature_point (const unsigned int i) const
{
- Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
- Assert (i<quadrature_points.size(), ExcIndexRange(i, 0, quadrature_points.size()));
+ Assert (this->update_flags & update_q_points, ExcAccessToUninitializedField());
+ Assert (i<this->quadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size()));
- return quadrature_points[i];
+ return this->quadrature_points[i];
};
template <int dim>
double FEValuesBase<dim>::JxW (const unsigned int i) const
{
- Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
- Assert (i<JxW_values.size(), ExcIndexRange(i, 0, JxW_values.size()));
+ Assert (this->update_flags & update_JxW_values, ExcAccessToUninitializedField());
+ Assert (i<this->JxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size()));
- return JxW_values[i];
+ return this->JxW_values[i];
};
unsigned int
FEValuesBase<dim>::memory_consumption () const
{
- return (MemoryConsumption::memory_consumption (shape_values) +
- MemoryConsumption::memory_consumption (shape_gradients) +
- MemoryConsumption::memory_consumption (shape_2nd_derivatives) +
- MemoryConsumption::memory_consumption (JxW_values) +
- MemoryConsumption::memory_consumption (quadrature_points) +
- MemoryConsumption::memory_consumption (normal_vectors) +
- MemoryConsumption::memory_consumption (boundary_forms) +
- sizeof(update_flags) +
+ return (MemoryConsumption::memory_consumption (this->shape_values) +
+ MemoryConsumption::memory_consumption (this->shape_gradients) +
+ MemoryConsumption::memory_consumption (this->shape_2nd_derivatives) +
+ MemoryConsumption::memory_consumption (this->JxW_values) +
+ MemoryConsumption::memory_consumption (this->quadrature_points) +
+ MemoryConsumption::memory_consumption (this->normal_vectors) +
+ MemoryConsumption::memory_consumption (this->boundary_forms) +
+ sizeof(this->update_flags) +
MemoryConsumption::memory_consumption (present_cell) +
MemoryConsumption::memory_consumption (n_quadrature_points) +
MemoryConsumption::memory_consumption (dofs_per_cell) +
MemoryConsumption::memory_consumption (*mapping_data) +
MemoryConsumption::memory_consumption (fe_data) +
MemoryConsumption::memory_consumption (*fe_data) +
- MemoryConsumption::memory_consumption (shape_function_to_row_table));
+ MemoryConsumption::memory_consumption (this->shape_function_to_row_table));
};
update_default,
mapping,
fe),
- quadrature (q)
+ quadrature (q)
{
initialize (update_flags);
};
FEValuesBase<dim> (q.n_quadrature_points,
fe.dofs_per_cell,
update_default,
- get_default_mapping(),
+ this->get_default_mapping(),
fe),
- quadrature (q)
+ quadrature (q)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
initialize (update_flags);
Assert ((update_flags & update_normal_vectors) == false,
typename FEValuesBase<dim>::ExcInvalidUpdateFlag());
- const UpdateFlags flags = compute_update_flags (update_flags);
+ const UpdateFlags flags = this->compute_update_flags (update_flags);
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
// calls to reinit
- mapping_data = mapping->get_data(flags, quadrature);
- fe_data = fe->get_data(flags, *mapping, quadrature);
+ this->mapping_data = this->mapping->get_data(flags, quadrature);
+ this->fe_data = this->fe->get_data(flags, *this->mapping, quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize (n_quadrature_points, *fe, flags);
+ FEValuesData<dim>::initialize (this->n_quadrature_points, *this->fe, flags);
};
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*fe) ==
+ Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
typename FEValuesBase<dim>::ExcFEDontMatch());
- present_cell = cell;
+ this->present_cell = cell;
- get_mapping().fill_fe_values(cell,
- quadrature,
- *mapping_data,
- quadrature_points,
- JxW_values);
+ this->get_mapping().fill_fe_values(cell,
+ quadrature,
+ *this->mapping_data,
+ this->quadrature_points,
+ this->JxW_values);
- get_fe().fill_fe_values(get_mapping(),
- cell,
- quadrature,
- *mapping_data,
- *fe_data,
- *this);
+ this->get_fe().fill_fe_values(this->get_mapping(),
+ cell,
+ quadrature,
+ *this->mapping_data,
+ *this->fe_data,
+ *this);
}
update_default,
mapping,
fe),
- quadrature(quadrature)
+ quadrature(quadrature)
{};
const typename std::vector<Point<dim> > &
FEFaceValuesBase<dim>::get_normal_vectors () const
{
- Assert (update_flags & update_normal_vectors,
+ Assert (this->update_flags & update_normal_vectors,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
- return normal_vectors;
+ return this->normal_vectors;
};
const typename std::vector<Tensor<1,dim> > &
FEFaceValuesBase<dim>::get_boundary_forms () const
{
- Assert (update_flags & update_boundary_forms,
+ Assert (this->update_flags & update_boundary_forms,
FEValuesBase<dim>::ExcAccessToUninitializedField());
- return boundary_forms;
+ return this->boundary_forms;
};
FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
fe.dofs_per_cell,
update_flags,
- get_default_mapping(),
+ this->get_default_mapping(),
fe, quadrature)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
void
FEFaceValues<dim>::initialize (const UpdateFlags update_flags)
{
- const UpdateFlags flags = compute_update_flags (update_flags);
+ const UpdateFlags flags = this->compute_update_flags (update_flags);
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
// calls to reinit
- mapping_data = mapping->get_face_data(flags, quadrature);
- fe_data = fe->get_face_data(flags, *mapping, quadrature);
+ this->mapping_data = this->mapping->get_face_data(flags, this->quadrature);
+ this->fe_data = this->fe->get_face_data(flags, *this->mapping, this->quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize(n_quadrature_points, *fe, flags);
+ FEValuesData<dim>::initialize(this->n_quadrature_points, *this->fe, flags);
};
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*fe) ==
+ Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
typename FEValuesBase<dim>::ExcFEDontMatch());
- present_cell = cell;
- present_face = cell->face(face_no);
+ this->present_cell = cell;
+ this->present_face = cell->face(face_no);
- get_mapping().fill_fe_face_values(cell, face_no,
- quadrature,
- *mapping_data,
- quadrature_points,
- JxW_values,
- boundary_forms,
- normal_vectors);
+ this->get_mapping().fill_fe_face_values(cell, face_no,
+ this->quadrature,
+ *this->mapping_data,
+ this->quadrature_points,
+ this->JxW_values,
+ this->boundary_forms,
+ this->normal_vectors);
- get_fe().fill_fe_face_values(get_mapping(),
- cell, face_no,
- quadrature,
- *mapping_data,
- *fe_data,
- *this);
+ this->get_fe().fill_fe_face_values(this->get_mapping(),
+ cell, face_no,
+ this->quadrature,
+ *this->mapping_data,
+ *this->fe_data,
+ *this);
};
FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
fe.dofs_per_cell,
update_flags,
- get_default_mapping(),
+ this->get_default_mapping(),
fe, quadrature)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
void
FESubfaceValues<dim>::initialize (const UpdateFlags update_flags)
{
- const UpdateFlags flags = compute_update_flags (update_flags);
+ const UpdateFlags flags = this->compute_update_flags (update_flags);
// then get objects into which the
// FE and the Mapping can store
// intermediate data used across
// calls to reinit
- mapping_data = mapping->get_subface_data(flags, quadrature);
- fe_data = fe->get_subface_data(flags, *mapping, quadrature);
+ this->mapping_data = this->mapping->get_subface_data(flags, this->quadrature);
+ this->fe_data = this->fe->get_subface_data(flags,
+ *this->mapping,
+ this->quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize(n_quadrature_points, *fe, flags);
+ FEValuesData<dim>::initialize(this->n_quadrature_points, *this->fe, flags);
};
// passed to the constructor and
// used by the DoFHandler used by
// this cell, are the same
- Assert (static_cast<const FiniteElementData<dim>&>(*fe) ==
+ Assert (static_cast<const FiniteElementData<dim>&>(*this->fe) ==
static_cast<const FiniteElementData<dim>&>(cell->get_dof_handler().get_fe()),
typename FEValuesBase<dim>::ExcFEDontMatch());
Assert (face_no < GeometryInfo<dim>::faces_per_cell,
Assert (subface_no < GeometryInfo<dim>::subfaces_per_face,
ExcIndexRange (subface_no, 0, GeometryInfo<dim>::subfaces_per_face));
- present_cell = cell;
- present_face = cell->face(face_no);
-
- get_mapping().fill_fe_subface_values(cell, face_no, subface_no,
- quadrature,
- *mapping_data,
- quadrature_points,
- JxW_values,
- boundary_forms,
- normal_vectors);
+ this->present_cell = cell;
+ this->present_face = cell->face(face_no);
+
+ this->get_mapping().fill_fe_subface_values(cell, face_no, subface_no,
+ this->quadrature,
+ *this->mapping_data,
+ this->quadrature_points,
+ this->JxW_values,
+ this->boundary_forms,
+ this->normal_vectors);
- get_fe().fill_fe_subface_values(get_mapping(),
- cell, face_no, subface_no,
- quadrature,
- *mapping_data,
- *fe_data,
- *this);
+ this->get_fe().fill_fe_subface_values(this->get_mapping(),
+ cell, face_no, subface_no,
+ this->quadrature,
+ *this->mapping_data,
+ *this->fe_data,
+ *this);
};
MappingCartesian<dim>::get_data (const UpdateFlags update_flags,
const Quadrature<dim> &q) const
{
- // Assert (flags & update_normal_vectors == 0, ExcNotImplemented());
+ // Assert (flags & update_normal_vectors == 0, ExcNotImplemented());
InternalData* data = new InternalData (q);
{
// 2D
case 200:
- n (1) = -1.;
- break;
+ n (1) = -1.;
+ break;
case 201:
- n (0) = 1.;
- break;
+ n (0) = 1.;
+ break;
case 202:
- n (1) = 1.;
- break;
+ n (1) = 1.;
+ break;
case 203:
- n (0) = -1.;
- break;
- // 3D
+ n (0) = -1.;
+ break;
+ // 3D
case 300:
- n (1) = -1.;
- break;
+ n (1) = -1.;
+ break;
case 301:
- n (1) = 1.;
- break;
+ n (1) = 1.;
+ break;
case 302:
- n (2) = -1.;
- break;
+ n (2) = -1.;
+ break;
case 303:
- n (0) = 1.;
- break;
+ n (0) = 1.;
+ break;
case 304:
- n (2) = 1.;
- break;
+ n (2) = 1.;
+ break;
case 305:
- n (0) = -1.;
- break;
+ n (0) = -1.;
+ break;
default:
- Assert (false, ExcInternalError());
+ Assert (false, ExcInternalError());
}
// furthermore, all normal
// vectors on a face are equal
// since the jacobian is diagonal
double J = 1.;
for (unsigned int d=0;d<dim;++d)
- if (d != (normal_directions[face_no]/2))
+ if (d != (this->normal_directions[face_no]/2))
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
template <int dim>
void
MappingCartesian<dim>::fill_fe_subface_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int sub_no,
- const Quadrature<dim-1> &q,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &boundary_forms,
- typename std::vector<Point<dim> > &normal_vectors) const
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1> &q,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ typename std::vector<Tensor<1,dim> > &boundary_forms,
+ typename std::vector<Point<dim> > &normal_vectors) const
{
// convert data object to internal
// data for this class. fails with
// since the jacobian is diagonal
double J = 1.;
for (unsigned int d=0;d<dim;++d)
- if (d != (normal_directions[face_no]/2))
+ if (d != (this->normal_directions[face_no]/2))
J *= data.length[d];
if (data.current_update_flags() & update_JxW_values)
MappingQ<dim>::InternalData::InternalData (const unsigned int n_shape_functions)
:
MappingQ1<dim>::InternalData(n_shape_functions),
- use_mapping_q1_on_current_cell(false),
- mapping_q1_data(1 << dim)
+ use_mapping_q1_on_current_cell(false),
+ mapping_q1_data(1 << dim)
{
- is_mapping_q1_data=false;
+ this->is_mapping_q1_data=false;
}
}
}
else
- // not precomputed, then do so now
+ // not precomputed, then do so now
compute_laplace_vector(lohvs);
// the sum of weights of the points
a.resize(GeometryInfo<dim>::vertices_per_cell);
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(vertex_mapping[i]);
+ a[i] = cell->vertex(this->vertex_mapping[i]);
}
}
switch (dim)
{
case 2:
- // in 2d, add the
- // points on the four
- // bounding lines to
- // the exterior (outer)
- // points
- add_line_support_points (cell, a);
- apply_laplace_vector (laplace_on_quad_vector,a);
- break;
+ // in 2d, add the
+ // points on the four
+ // bounding lines to
+ // the exterior (outer)
+ // points
+ add_line_support_points (cell, a);
+ apply_laplace_vector (laplace_on_quad_vector,a);
+ break;
case 3:
- // in 3d also add the
- // points located on
- // the boundary faces
- add_line_support_points (cell, a);
- add_quad_support_points (cell, a);
- apply_laplace_vector (laplace_on_hex_vector, a);
- break;
+ // in 3d also add the
+ // points located on
+ // the boundary faces
+ add_line_support_points (cell, a);
+ add_quad_support_points (cell, a);
+ apply_laplace_vector (laplace_on_hex_vector, a);
+ break;
default:
- Assert(false, ExcNotImplemented());
- break;
+ Assert(false, ExcNotImplemented());
+ break;
};
}
switch (degree)
{
case 2:
- {
- a.push_back(middle);
- break;
- };
+ {
+ a.push_back(middle);
+ break;
+ };
case 3:
- {
- // The four points in the
- // cell are located at
- // the midpoint between
- // the middle point and
- // the 4 vertices
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.);
- break;
- };
+ {
+ // The four points in the
+ // cell are located at
+ // the midpoint between
+ // the middle point and
+ // the 4 vertices
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ a.push_back(middle*2./3.+cell->vertex(this->vertex_mapping[i])/3.);
+ break;
+ };
case 4:
- {
- Assert(a.size()==16, ExcInternalError());
- a.insert(a.end(), 9, Point<dim>());
+ {
+ Assert(a.size()==16, ExcInternalError());
+ a.insert(a.end(), 9, Point<dim>());
- const unsigned int inner_map[8]=
- { 0, 1, 2, 5, 8, 7, 6, 3 };
+ const unsigned int inner_map[8]=
+ { 0, 1, 2, 5, 8, 7, 6, 3 };
- // The nine points in the
- // cell are located at the
- // midpoint between the
- // middle point and (the 4
- // vertices and the face
- // midpoints)
+ // The nine points in the
+ // cell are located at the
+ // midpoint between the
+ // middle point and (the 4
+ // vertices and the face
+ // midpoints)
- a[16+4]=middle;
- for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- {
- a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
- a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
- }
- break;
- };
+ a[16+4]=middle;
+ for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
+ a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
+ }
+ break;
+ };
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
};
};
else
tensor = data->covariant.begin();
}
- while (begin!=end)
+ while (begin!=end)
{
contract (*(begin++), *(src++), *(tensor++));
}
else
tensor = data->contravariant.begin();
}
- while (begin!=end)
+ while (begin!=end)
{
contract (*(begin++), *(tensor++), *(src++));
}
Assert(mdata!=0, ExcInternalError());
mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
- || cell->has_boundary_lines());
+ || cell->has_boundary_lines());
typename MappingQ1<dim>::InternalData *p_data=0;
if (mdata->use_mapping_q1_on_current_cell)
template<>
const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[2] =
{ 0, 1
-};
+ };
template<>
MappingQ1<2>::vertex_mapping[4] =
{
0, 1, 3, 2
-};
+ };
template<>
const unsigned int MappingQ1<deal_II_dimension>::vertex_mapping[8] =
{
0, 1, 4, 5, 3, 2, 7, 6
-};
+ };
template<>
// ignored for the interior of a
// cell.
if (out & (update_JxW_values
- |update_normal_vectors))
+ |update_normal_vectors))
out |= update_boundary_forms;
if (out & (update_covariant_transformation
- | update_JxW_values
- | update_boundary_forms
- | update_normal_vectors))
+ | update_JxW_values
+ | update_boundary_forms
+ | update_normal_vectors))
out |= update_contravariant_transformation;
return out;
// dim so we can subtract 1
// without getting negative
// values.
- unsigned int nindex = normal_directions[i]/2 + dim;
+ unsigned int nindex = this->normal_directions[i]/2 + dim;
// First tangential has a
// non-zero in component
// (i+1)%dim, if normal is
// non-zero in i.
Tensor<1,dim> tangential;
- tangential[(nindex+1)%dim] = (normal_directions[i]%2) ? -1 : 1;
+ tangential[(nindex+1)%dim] = (this->normal_directions[i]%2) ? -1 : 1;
data.unit_tangentials[i].resize(n_original_q_points);
std::fill (data.unit_tangentials[i].begin(),
data.unit_tangentials[i].end(),
Assert (JxW_values.size() == npts,
ExcDimensionMismatch(JxW_values.size(), npts));
for (unsigned int point=0; point<npts; ++point)
- JxW_values[point]
- = determinant(data.contravariant[point])*weights[point];
+ JxW_values[point]
+ = determinant(data.contravariant[point])*weights[point];
}
}
switch (dim)
{
case 2:
- {
- for (; result != end; ++result, ++tang1)
- cross_product (*result, *tang1);
- break;
- };
+ {
+ for (; result != end; ++result, ++tang1)
+ cross_product (*result, *tang1);
+ break;
+ };
case 3:
- {
- transform_contravariant(data.aux[1].begin(),
- data.aux[1].end(),
- data.unit_tangentials[
- face_no+GeometryInfo<dim>::faces_per_cell].begin(),
- data);
- typename std::vector<Tensor<1,dim> >::const_iterator
- tang2 = data.aux[1].begin();
- for (;result != end; ++result, ++tang1, ++tang2)
- cross_product (*result, *tang1, *tang2);
- break;
- };
+ {
+ transform_contravariant(data.aux[1].begin(),
+ data.aux[1].end(),
+ data.unit_tangentials[
+ face_no+GeometryInfo<dim>::faces_per_cell].begin(),
+ data);
+ typename std::vector<Tensor<1,dim> >::const_iterator
+ tang2 = data.aux[1].begin();
+ for (;result != end; ++result, ++tang1, ++tang2)
+ cross_product (*result, *tang1, *tang2);
+ break;
+ };
default:
- Assert(false, ExcNotImplemented());
+ Assert(false, ExcNotImplemented());
}
if (update_flags & (update_normal_vectors
// are always numbered first,
// we can access them easily
for (unsigned int j=0; j<dim; ++j)
- shift_vector[j] = mapping_values(vertex_mapping[i]*dim+j);
+ shift_vector[j] = mapping_values(this->vertex_mapping[i]*dim+j);
// compute new support point by
// old (reference) value and
// added shift
- a[i] = cell->vertex(vertex_mapping[i]) + shift_vector;
+ a[i] = cell->vertex(this->vertex_mapping[i]) + shift_vector;
}
}