*/
const unsigned int total_dofs;
+ /**
+ * Number of basis functions used for the
+ * transformation from unit cell to real
+ * cell. For a linear mapping, this number
+ * equals the number of vertices.
+ */
+ const unsigned int n_transform_functions;
+
/**
* Default constructor. Constructs an element
* which is not so useful. Checking
* #total_dofs# is therefore a good way to
* check if something went wrong.
*/
- FiniteElementData () :
- dofs_per_vertex(0),
- dofs_per_line(0),
- dofs_per_face(0),
- total_dofs(0) {};
+ FiniteElementData ();
/**
* A more useful version to construct
* an object of this type.
*/
FiniteElementData (const unsigned int dofs_per_vertex,
- const unsigned int dofs_per_line) :
+ const unsigned int dofs_per_line,
+ const unsigned int n_transform_functions) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_face(dofs_per_vertex),
- total_dofs (2*dofs_per_vertex+dofs_per_line) {};
+ total_dofs (2*dofs_per_vertex+dofs_per_line),
+ n_transform_functions(n_transform_functions) {};
/**
* Comparison operator. It is not clear to
* this one explicitely.
*/
bool operator == (const FiniteElementData<1> &) const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInternalError);
};
*/
const unsigned int total_dofs;
+ /**
+ * Number of basis functions used for the
+ * transformation from unit cell to real
+ * cell. For a linear mapping, this number
+ * equals the number of vertices.
+ */
+ const unsigned int n_transform_functions;
+
/**
* Default constructor. Constructs an element
* which is not so useful. Checking
* #total_dofs# is therefore a good way to
* check if something went wrong.
*/
- FiniteElementData () :
- dofs_per_vertex(0),
- dofs_per_line(0),
- dofs_per_quad(0),
- dofs_per_face(0),
- total_dofs(0) {};
+ FiniteElementData ();
/**
* A more useful version to construct
*/
FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad) :
+ const unsigned int dofs_per_quad,
+ const unsigned int n_transform_functions) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_quad(dofs_per_quad),
dofs_per_line),
total_dofs (4*dofs_per_vertex+
4*dofs_per_line+
- dofs_per_quad) {};
+ dofs_per_quad),
+ n_transform_functions (n_transform_functions) {};
/**
* Comparison operator. It is not clear to
* this one explicitely.
*/
bool operator == (const FiniteElementData<2> &) const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInternalError);
};
*/
FiniteElementBase (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad=0);
+ const unsigned int dofs_per_quad,
+ const unsigned int n_transform_functions);
/**
* Return a readonly reference to the
<< "The interface matrix has a size of " << arg1
<< "x" << arg2
<< ", which is not reasonable in the present dimension.");
- /**
- * Exception
- */
- DeclException0 (ExcInternalError);
protected:
/**
*/
FiniteElement (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad=0) :
+ const unsigned int dofs_per_quad,
+ const unsigned int n_transform_functions) :
FiniteElementBase<dim> (dofs_per_vertex,
dofs_per_line,
- dofs_per_quad) {};
+ dofs_per_quad,
+ n_transform_functions) {};
/**
* Destructor. Only declared to have a
* #p# is a point on the reference element.
*/
virtual double shape_value (const unsigned int i,
- const Point<dim>& p) const = 0;
+ const Point<dim> &p) const = 0;
/**
* Return the gradient of the #i#th shape
* #p# is a point on the reference element,
*/
virtual Point<dim> shape_grad (const unsigned int i,
- const Point<dim>& p) const = 0;
+ const Point<dim> &p) const = 0;
+ /**
+ * Return the value of the #i#th shape
+ * function of the transformation mapping
+ * from unit cell to real cell. For
+ * isoparametric elements, this function
+ * is the same as the ansatz functions,
+ * but for sublinear or other mappings,
+ * they differ.
+ */
+ virtual double shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const = 0;
+
+ /**
+ * Same as above: return gradient of the
+ * #i#th shape function for the mapping
+ * from unit to real cell.
+ */
+ virtual Point<dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const = 0;
+
/**
* Compute the Jacobian matrix and the
* quadrature points as well as the ansatz
* unit cell to the real cell, which
* makes things much more complicated.
*
+ * The #shape_values/grads_transform#
+ * arrays store the values and gradients
+ * of the transformation basis functions.
+ * While this information is not necessary
+ * for the computation of the other fields,
+ * it allows for significant speedups, since
+ * the values and gradients of the transform
+ * functions at the quadrature points
+ * need not be recomputed each time this
+ * function is called.
+ *
* The function assumes that the fields
* already have the right number of
* elements.
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
const bool compute_face_jacobians,
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
const bool compute_face_jacobians,
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grads_transform,
const Boundary<dim> &boundary) const;
/**
const unsigned int dofs_per_quad=0) :
FiniteElement<dim> (dofs_per_vertex,
dofs_per_line,
- dofs_per_quad) {};
+ dofs_per_quad,
+ GeometryInfo<dim>::vertices_per_cell) {};
+
+ /**
+ * Return the value of the #i#th shape
+ * function at point #p# on the unit cell.
+ * Here, the (bi-)linear basis functions
+ * are meant, which are used for the
+ * computation of the transformation from
+ * unit cell to real space cell.
+ */
+ virtual double shape_value_transform (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Return the gradient of the #i#th shape
+ * function at point #p# on the unit cell.
+ * Here, the (bi-)linear basis functions
+ * are meant, which are used for the
+ * computation of the transformation from
+ * unit cell to real space cell.
+ */
+ virtual Point<dim> shape_grad_transform (const unsigned int i,
+ const Point<dim> &p) const;
/**
* Refer to the base class for detailed
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const;
-
- protected:
- /**
- * Return the value of the #i#th shape
- * function at point #p# on the unit cell.
- * Here, the (bi-)linear basis functions
- * are meant, which are used for the
- * computation of the transformation from
- * unit cell to real space cell.
- */
- double linear_shape_value(const unsigned int i,
- const Point<dim>& p) const;
-
- /**
- * Return the gradient of the #i#th shape
- * function at point #p# on the unit cell.
- * Here, the (bi-)linear basis functions
- * are meant, which are used for the
- * computation of the transformation from
- * unit cell to real space cell.
- */
- Point<dim> linear_shape_grad(const unsigned int i,
- const Point<dim>& p) const;
};
* the derived classes and the #get_values# function for the exact usage of
* this variable.
*
+ * For many of the actual computations done by the #fill_fe_*# functions of
+ * the #FiniteElement# class and its decendants, the values and gradients of
+ * the transformation functions are needed. For example, for the computation
+ * of the real location of a quadrature point from the location on the unit
+ * cell, the values are needed, while for the computation of the Jacobian
+ * matrix the gradient is needed. While for linear elements the transformation
+ * functions coincide with the ansatz functions, this does not hold for higher
+ * order elements with subparametric mappings and for other types of elements
+ * such as non-conforming ones, etc, such that the precomputed values and
+ * gradients of the ansatz functions (#unit_shape_values# and
+ * #unit_shape_grads#) cannot be used for the present purpose.
+ * In principle, these values could be computed each time the #fill_fe_*#
+ * function is called; however, this computation is highly redundant, since
+ * only the values on the unit cell and only at the quadrature points are
+ * needed, i.e. they are the same for each cell that #fill_fe_*# is called.
+ * Therefore, two additional arrays, #unit_shape_values_transform# and
+ * #unit_shape_grads_transform# are provided, which are filled upon construction
+ * of an object of this type, which the actual finite element may or may not
+ * use. Later on, the #fill_fe_*# functions are passed pointers to these
+ * arrays, which they may use to extract the values and gradients of the
+ * transform functions. If a concrete finite element choses not to use this
+ * field, it shall set its field #n_transform_functions# to zero.
+ *
+ * The #unit_shape_grads_transform# array is provided by the derived classes
+ * to allow for the inclusion of multiple faces, etc.
+ *
*
* \subsection{Definitions}
*
* not assume that someone who wants to get the #JxW_values# must set the
* #update_jacobians# flag besides the #update_JxW_values# flag.
*
- *
* It is also forbidden that the constructor of a class set the
* #update_jacobians# flag if the user specifies #update_JxW_values#. This is
* since derived classes may be able to compute the #JxW_values# field without
*/
const unsigned int total_dofs;
+ /**
+ * Number of basis functions for the
+ * transformation from the unit cell
+ * to the real cell. See the docs for
+ * more information on this field.
+ */
+ const unsigned int n_transform_functions;
+
/**
* Constructor. Set up the array sizes
* with #n_q_points# quadrature points,
FEValuesBase (const unsigned int n_q_points,
const unsigned int n_ansatz_points,
const unsigned int n_dofs,
+ const unsigned int n_transform_functions,
const unsigned int n_values_array,
const UpdateFlags update_flags);
*/
vector<dFMatrix> jacobi_matrices;
+ /**
+ * Store the values of the basis functions
+ * of the transformation from unit cell
+ * to real cell at the quadrature points.
+ *
+ * This field stores some data which is not
+ * really needed for the assemblage of
+ * matrices and vectors but makes that
+ * operation much faster. Each time the
+ * #FEValues::reinit# function calls
+ * the #FiniteElemenet::fill_fe_values#
+ * function, this and the next array are
+ * passed. The #fill_fe_values# function
+ * may or may not use this field.
+ *
+ * The element #(i,j)# denotes the value
+ * of the #i#th transfer basis function
+ * at the #j#th quadrature point.
+ */
+ vector<dFMatrix> shape_values_transform;
+
/**
* Store which of the data sets in the
* #shape_values# array is presently
* Store which fields are to be updated by
* the reinit function.
*/
- UpdateFlags update_flags;
+ UpdateFlags update_flags;
/**
* Store the cell selected last time
*/
vector<vector<Point<dim> > > unit_shape_gradients;
-
+ /**
+ * Gradients of the basis
+ * functions of the transformation.
+ * Analogous to the #shape_values_transform#
+ * array of the base class.
+ */
+ vector<vector<Point<dim> > > unit_shape_gradients_transform;
+
/**
* Array of quadrature points in the unit
* cell. This array is set up upon
FEFaceValuesBase (const unsigned int n_q_points,
const unsigned int n_ansatz_points,
const unsigned int n_dofs,
+ const unsigned int n_transform_functions,
const unsigned int n_faces_or_subfaces,
const UpdateFlags update_flags);
* #unit_shape_gradients[face][dof][q_point]#
*/
vector<vector<vector<Point<dim> > > > unit_shape_gradients;
+
+ /**
+ * Gradients of the basis
+ * functions of the transformation.
+ * Analogous to the #shape_values_transform#
+ * array of the base class.
+ */
+ vector<vector<vector<Point<dim> > > > unit_shape_gradients_transform;
/**
* Array of quadrature points on the
#if deal_II_dimension == 1
+FiniteElementData<1>::FiniteElementData () :
+ dofs_per_vertex(0),
+ dofs_per_line(0),
+ dofs_per_face(0),
+ total_dofs(0),
+ n_transform_functions(0) {
+ Assert (false, ExcInternalError());
+};
+
+
+
bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const {
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
#if deal_II_dimension == 2
+
+FiniteElementData<2>::FiniteElementData () :
+ dofs_per_vertex(0),
+ dofs_per_line(0),
+ dofs_per_quad(0),
+ dofs_per_face(0),
+ total_dofs(0),
+ n_transform_functions(0) {
+ Assert (false, ExcInternalError());
+};
+
+
+
bool FiniteElementData<2>::operator== (const FiniteElementData<2> &f) const {
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
template <>
FiniteElementBase<1>::FiniteElementBase (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad) :
+ const unsigned int dofs_per_quad,
+ const unsigned int n_transform_funcs) :
FiniteElementData<1> (dofs_per_vertex,
- dofs_per_line)
+ dofs_per_line,
+ n_transform_functs)
{
Assert (dofs_per_quad==0, ExcInternalError());
template <>
FiniteElementBase<2>::FiniteElementBase (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad) :
+ const unsigned int dofs_per_quad,
+ const unsigned int n_transform_funcs) :
FiniteElementData<2> (dofs_per_vertex,
dofs_per_line,
- dofs_per_quad)
+ dofs_per_quad,
+ n_transform_funcs)
{
const unsigned int dim=2;
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grad_transform,
const Boundary<1> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
const bool ,
vector<Point<1> > &,
const bool,
+ const dFMatrix &,
+ const vector<vector<Point<1> > > &,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
};
const bool ,
vector<Point<1> > &,
const bool,
+ const dFMatrix &,
+ const vector<vector<Point<1> > > &,
const Boundary<1> &) const {
Assert (false, ExcNotImplemented());
};
const bool,
vector<Point<dim> > &,
const bool,
+ const dFMatrix &,
+ const vector<vector<Point<dim> > > &,
const Boundary<dim> &) const {
Assert (false, ExcPureFunctionCalled());
};
const bool compute_face_jacobians,
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_gradients_transform,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
jacobians, compute_jacobians,
dummy, false,
q_points, compute_q_points,
+ shape_values_transform, shape_gradients_transform,
boundary);
if (compute_ansatz_points)
const bool compute_face_jacobians,
vector<Point<dim> > &normal_vectors,
const bool compute_normal_vectors,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_gradients_transform,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
jacobians, compute_jacobians,
dummy, false,
q_points, compute_q_points,
+ shape_values_transform, shape_gradients_transform,
boundary);
if (compute_face_jacobians)
+
+/*---------------------------- FELinearMapping ----------------------------------*/
+
#if deal_II_dimension == 1
template <>
inline
double
-FELinearMapping<1>::linear_shape_value(const unsigned int i,
- const Point<1> &p) const
+FELinearMapping<1>::shape_value_transform (const unsigned int i,
+ const Point<1> &p) const
{
Assert((i<2), ExcInvalidIndex(i));
const double xi = p(0);
template <>
inline
Point<1>
-FELinearMapping<1>::linear_shape_grad(const unsigned int i,
- const Point<1>&) const
+FELinearMapping<1>::shape_grad_transform(const unsigned int i,
+ const Point<1>&) const
{
Assert((i<2), ExcInvalidIndex(i));
switch (i)
template <>
-void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
+void FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
const Boundary<1> &,
const vector<Point<0> > &,
vector<double> &) const {
template <>
-void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
+void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
const unsigned int ,
const vector<Point<0> > &,
vector<double> &) const {
template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
const unsigned int,
const Boundary<1> &,
const vector<Point<0> > &,
template <>
-void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
+void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
const unsigned int,
const unsigned int,
const vector<Point<0> > &,
const bool compute_ansatz_points,
vector<Point<1> > &q_points,
const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_gradients_transform,
const Boundary<1> &boundary) const {
// simply pass down
FiniteElement<1>::fill_fe_values (cell, unit_points,
jacobians, compute_jacobians,
ansatz_points, compute_ansatz_points,
- q_points, compute_q_points, boundary);
+ q_points, compute_q_points,
+ shape_values_transform, shape_gradients_transform,
+ boundary);
};
template <>
inline
double
-FELinearMapping<2>::linear_shape_value (const unsigned int i,
- const Point<2>& p) const
+FELinearMapping<2>::shape_value_transform (const unsigned int i,
+ const Point<2>& p) const
{
Assert((i<4), ExcInvalidIndex(i));
switch (i)
template <>
inline
Point<2>
-FELinearMapping<2>::linear_shape_grad (const unsigned int i,
- const Point<2>& p) const
+FELinearMapping<2>::shape_grad_transform (const unsigned int i,
+ const Point<2>& p) const
{
Assert((i<4), ExcInvalidIndex(i));
switch (i)
const bool compute_ansatz_points,
vector<Point<dim> > &q_points,
const bool compute_q_points,
+ const dFMatrix &shape_values_transform,
+ const vector<vector<Point<dim> > > &shape_grad_transform,
const Boundary<dim> &boundary) const {
Assert (jacobians.size() == unit_points.size(),
ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
// use a subparametric mapping
for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
for (unsigned int l=0; l<n_points; ++l)
- q_points[l] += vertices[j] * linear_shape_value(j, unit_points[l]);
+ q_points[l] += vertices[j] * shape_values_transform(j, l);
};
{
// we want the linear transform,
// so use that function
- const Point<dim> gradient = linear_shape_grad (s, unit_points[l]);
+ const Point<dim> gradient = shape_grad_transform[s][l];
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
M(i,j) += vertices[s](i) * gradient(j);
/*------------------------------- Explicit Instantiations -------------*/
-template class FiniteElementData<deal_II_dimension>;
template class FiniteElementBase<deal_II_dimension>;
template class FiniteElement<deal_II_dimension>;
template class FELinearMapping<deal_II_dimension>;
FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
const unsigned int n_ansatz_points,
const unsigned int n_dofs,
+ const unsigned int n_transform_functions,
const unsigned int n_values_arrays,
const UpdateFlags update_flags) :
n_quadrature_points (n_q_points),
total_dofs (n_dofs),
+ n_transform_functions (n_transform_functions),
shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
shape_gradients (n_dofs, vector<Point<dim> >(n_q_points)),
weights (n_q_points, 0),
quadrature_points (n_q_points, Point<dim>()),
ansatz_points (n_ansatz_points, Point<dim>()),
jacobi_matrices (n_q_points, dFMatrix(dim,dim)),
+ shape_values_transform (n_values_arrays,
+ dFMatrix(n_transform_functions,
+ n_quadrature_points)),
selected_dataset (0),
update_flags (update_flags) {};
FEValuesBase<dim> (quadrature.n_quadrature_points,
fe.total_dofs,
fe.total_dofs,
+ fe.n_transform_functions,
1,
update_flags),
unit_shape_gradients(fe.total_dofs,
vector<Point<dim> >(quadrature.n_quadrature_points)),
+ unit_shape_gradients_transform(fe.n_transform_functions,
+ vector<Point<dim> >(quadrature.n_quadrature_points)),
unit_quadrature_points(quadrature.get_quad_points())
{
Assert ((update_flags & update_normal_vectors) == false,
= fe.shape_grad(i, unit_quadrature_points[j]);
};
+ for (unsigned int i=0; i<n_transform_functions; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ shape_values_transform[0] (i,j)
+ = fe.shape_value_transform (i, unit_quadrature_points[j]);
+ unit_shape_gradients_transform[i][j]
+ = fe.shape_grad_transform(i, unit_quadrature_points[j]);
+ };
+
weights = quadrature.get_weights ();
};
update_flags & update_ansatz_points,
quadrature_points,
update_flags & update_q_points,
+ shape_values_transform[0], unit_shape_gradients_transform,
boundary);
// compute gradients on real element if
FEFaceValuesBase<dim>::FEFaceValuesBase (const unsigned int n_q_points,
const unsigned int n_ansatz_points,
const unsigned int n_dofs,
+ const unsigned int n_transform_functions,
const unsigned int n_faces_or_subfaces,
const UpdateFlags update_flags) :
FEValuesBase<dim> (n_q_points,
n_ansatz_points,
n_dofs,
+ n_transform_functions,
n_faces_or_subfaces,
update_flags),
unit_shape_gradients (n_faces_or_subfaces,
vector<vector<Point<dim> > >(n_dofs,
vector<Point<dim> >(n_q_points))),
+ unit_shape_gradients_transform (n_faces_or_subfaces,
+ vector<vector<Point<dim> > >(n_transform_functions,
+ vector<Point<dim> >(n_q_points))),
unit_face_quadrature_points (n_q_points, Point<dim-1>()),
unit_quadrature_points (n_faces_or_subfaces,
vector<Point<dim> >(n_q_points, Point<dim>())),
FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
fe.dofs_per_face,
fe.total_dofs,
- 2*dim,
+ fe.n_transform_functions,
+ GeometryInfo<dim>::faces_per_cell,
update_flags)
{
unit_face_quadrature_points = quadrature.get_quad_points();
// of the space with #dim# dimensions.
// the points are still on the unit
// cell, not on the real cell.
- for (unsigned int face=0; face<2*dim; ++face)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
QProjector<dim>::project_to_face (quadrature, face, unit_quadrature_points[face]);
- for (unsigned int face=0; face<2*dim; ++face)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
{
unit_shape_gradients[face][i][j]
= fe.shape_grad(i, unit_quadrature_points[face][j]);
};
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<n_transform_functions; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ shape_values_transform[face] (i,j)
+ = fe.shape_value_transform (i, unit_quadrature_points[face][j]);
+ unit_shape_gradients_transform[face][i][j]
+ = fe.shape_grad_transform(i, unit_quadrature_points[face][j]);
+ };
};
update_flags & update_JxW_values,
normal_vectors,
update_flags & update_normal_vectors,
+ shape_values_transform[face_no],
+ unit_shape_gradients_transform[face_no],
boundary);
// compute gradients on real element if
FEFaceValuesBase<dim> (quadrature.n_quadrature_points,
0,
fe.total_dofs,
- 2*dim*(1<<(dim-1)),
+ fe.n_transform_functions,
+ GeometryInfo<dim>::faces_per_cell * GeometryInfo<dim>::subfaces_per_face,
update_flags)
{
Assert ((update_flags & update_ansatz_points) == false,
// of the space with #dim# dimensions.
// the points are still on the unit
// cell, not on the real cell.
- for (unsigned int face=0; face<2*dim; ++face)
- for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
QProjector<dim>::project_to_subface (quadrature,
face, subface,
unit_quadrature_points[face*(1<<(dim-1))+subface]);
- for (unsigned int face=0; face<2*dim; ++face)
- for (unsigned int subface=0; subface<(1<<(dim-1)); ++subface)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
{
- shape_values[face*(1<<(dim-1))+subface](i,j)
- = fe.shape_value(i, unit_quadrature_points[face*(1<<(dim-1))+subface][j]);
- unit_shape_gradients[face*(1<<(dim-1))+subface][i][j]
- = fe.shape_grad(i, unit_quadrature_points[face*(1<<(dim-1))+subface][j]);
+ shape_values[face*GeometryInfo<dim>::subfaces_per_face+subface](i,j)
+ = fe.shape_value(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+ unit_shape_gradients[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+ = fe.shape_grad(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+ };
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
+ for (unsigned int i=0; i<n_transform_functions; ++i)
+ for (unsigned int j=0; j<n_quadrature_points; ++j)
+ {
+ shape_values_transform[face*GeometryInfo<dim>::subfaces_per_face+subface] (i,j)
+ = fe.shape_value_transform (i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
+ unit_shape_gradients_transform[face*GeometryInfo<dim>::subfaces_per_face+subface][i][j]
+ = fe.shape_grad_transform(i, unit_quadrature_points[face*GeometryInfo<dim>::subfaces_per_face+subface][j]);
};
};
update_flags & update_JxW_values,
normal_vectors,
update_flags & update_normal_vectors,
+ shape_values_transform[selected_dataset],
+ unit_shape_gradients_transform[selected_dataset],
boundary);
// compute gradients on real element if