-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int n_q_points,
- const DataSetDescriptor data_set,
- const CellSimilarity::Similarity cell_similarity,
- const InternalData &data,
- std::vector<Point<spacedim> > &quadrature_points) const
+namespace internal
{
- const UpdateFlags update_flags(data.current_update_flags());
-
- // if necessary, recompute the
- // support points of the
- // transformation of this cell
- // (note that we need to first
- // check the triangulation pointer,
- // since otherwise the second test
- // might trigger an exception if
- // the triangulations are not the
- // same)
- if ((data.mapping_support_points.size() == 0)
- ||
- (&cell->get_triangulation() !=
- &data.cell_of_current_support_points->get_triangulation())
- ||
- (cell != data.cell_of_current_support_points))
- {
- compute_mapping_support_points(cell, data.mapping_support_points);
- data.cell_of_current_support_points = cell;
- }
-
- // first compute quadrature points
- if (update_flags & update_quadrature_points)
- {
- AssertDimension (quadrature_points.size(), n_q_points);
-
- for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double *shape = &data.shape(point+data_set,0);
- Point<spacedim> result = (shape[0] *
- data.mapping_support_points[0]);
- for (unsigned int k=1; k<data.n_shape_functions; ++k)
- for (unsigned int i=0; i<spacedim; ++i)
- result[i] += shape[k] * data.mapping_support_points[k][i];
- quadrature_points[point] = result;
- }
- }
-
- // then Jacobians
- if (update_flags & update_contravariant_transformation)
+ namespace
+ {
+ /**
+ * Compute the locations of quadrature points on the object described by
+ * the first argument (and the cell for which the mapping support points
+ * have already been set), but only if the update_flags of the @p data
+ * argument indicate so.
+ */
+ template <int dim, int spacedim>
+ void
+ maybe_compute_q_points (const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor data_set,
+ const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
+ std::vector<Point<spacedim> > &quadrature_points)
{
- AssertDimension (data.contravariant.size(), n_q_points);
+ const UpdateFlags update_flags(data.current_update_flags());
- // if the current cell is just a
- // translation of the previous one, no
- // need to recompute jacobians...
- if (cell_similarity != CellSimilarity::translation)
+ if (update_flags & update_quadrature_points)
{
- std::fill(data.contravariant.begin(), data.contravariant.end(),
- DerivativeForm<1,dim,spacedim>());
-
- Assert (data.n_shape_functions > 0, ExcInternalError());
- const Tensor<1,spacedim> *supp_pts =
- &data.mapping_support_points[0];
-
- for (unsigned int point=0; point<n_q_points; ++point)
+ for (unsigned int point=0; point<quadrature_points.size(); ++point)
{
- const Tensor<1,dim> *data_derv =
- &data.derivative(point+data_set, 0);
-
- double result [spacedim][dim];
-
- // peel away part of sum to avoid zeroing the
- // entries and adding for the first time
- for (unsigned int i=0; i<spacedim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- result[i][j] = data_derv[0][j] * supp_pts[0][i];
+ const double *shape = &data.shape(point+data_set,0);
+ Point<spacedim> result = (shape[0] *
+ data.mapping_support_points[0]);
for (unsigned int k=1; k<data.n_shape_functions; ++k)
for (unsigned int i=0; i<spacedim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- result[i][j] += data_derv[k][j] * supp_pts[k][i];
-
- // write result into contravariant data. for
- // j=dim in the case dim<spacedim, there will
- // never be any nonzero data that arrives in
- // here, so it is ok anyway because it was
- // initialized to zero at the initialization
- for (unsigned int i=0; i<spacedim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- data.contravariant[point][i][j] = result[i][j];
+ result[i] += shape[k] * data.mapping_support_points[k][i];
+ quadrature_points[point] = result;
}
}
}
- if (update_flags & update_covariant_transformation)
- {
- AssertDimension (data.covariant.size(), n_q_points);
- if (cell_similarity != CellSimilarity::translation)
- for (unsigned int point=0; point<n_q_points; ++point)
- {
- data.covariant[point] = (data.contravariant[point]).covariant_form();
- }
- }
- if (update_flags & update_volume_elements)
- if (cell_similarity != CellSimilarity::translation)
- for (unsigned int point=0; point<n_q_points; ++point)
- data.volume_elements[point] = data.contravariant[point].determinant();
+ /**
+ * Update the co- and contravariant matrices as well as their determinant, for the cell
+ * described stored in the data object, but only if the update_flags of the @p data
+ * argument indicate so.
+ *
+ * Skip the computation if possible as indicated by the first argument.
+ */
+ template <int dim, int spacedim>
+ void
+ maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity,
+ const typename dealii::MappingQ1<dim,spacedim>::DataSetDescriptor data_set,
+ const typename dealii::MappingQ1<dim,spacedim>::InternalData &data)
+ {
+ const UpdateFlags update_flags(data.current_update_flags());
+ if (update_flags & update_contravariant_transformation)
+ // if the current cell is just a
+ // translation of the previous one, no
+ // need to recompute jacobians...
+ if (cell_similarity != CellSimilarity::translation)
+ {
+ const unsigned int n_q_points = data.contravariant.size();
+
+ std::fill(data.contravariant.begin(), data.contravariant.end(),
+ DerivativeForm<1,dim,spacedim>());
+
+ Assert (data.n_shape_functions > 0, ExcInternalError());
+ const Tensor<1,spacedim> *supp_pts =
+ &data.mapping_support_points[0];
+
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const Tensor<1,dim> *data_derv =
+ &data.derivative(point+data_set, 0);
+
+ double result [spacedim][dim];
+
+ // peel away part of sum to avoid zeroing the
+ // entries and adding for the first time
+ for (unsigned int i=0; i<spacedim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ result[i][j] = data_derv[0][j] * supp_pts[0][i];
+ for (unsigned int k=1; k<data.n_shape_functions; ++k)
+ for (unsigned int i=0; i<spacedim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ result[i][j] += data_derv[k][j] * supp_pts[k][i];
+
+ // write result into contravariant data. for
+ // j=dim in the case dim<spacedim, there will
+ // never be any nonzero data that arrives in
+ // here, so it is ok anyway because it was
+ // initialized to zero at the initialization
+ for (unsigned int i=0; i<spacedim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ data.contravariant[point][i][j] = result[i][j];
+ }
+ }
+
+ if (update_flags & update_covariant_transformation)
+ if (cell_similarity != CellSimilarity::translation)
+ {
+ const unsigned int n_q_points = data.contravariant.size();
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ data.covariant[point] = (data.contravariant[point]).covariant_form();
+ }
+ }
+ if (update_flags & update_volume_elements)
+ if (cell_similarity != CellSimilarity::translation)
+ {
+ const unsigned int n_q_points = data.contravariant.size();
+ for (unsigned int point=0; point<n_q_points; ++point)
+ data.volume_elements[point] = data.contravariant[point].determinant();
+ }
+
+ }
+ }
}
-
template<int dim, int spacedim>
void
MappingQ1<dim,spacedim>::compute_mapping_support_points(
const unsigned int n_q_points=quadrature.size();
- compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity,
- data,
- output_data.quadrature_points);
+ // if necessary, recompute the support points of the transformation of this cell
+ // (note that we need to first check the triangulation pointer, since otherwise
+ // the second test might trigger an exception if the triangulations are not the
+ // same)
+ if ((data.mapping_support_points.size() == 0)
+ ||
+ (&cell->get_triangulation() !=
+ &data.cell_of_current_support_points->get_triangulation())
+ ||
+ (cell != data.cell_of_current_support_points))
+ {
+ compute_mapping_support_points(cell, data.mapping_support_points);
+ data.cell_of_current_support_points = cell;
+ }
+ internal::maybe_compute_q_points<dim,spacedim> (DataSetDescriptor::cell (),
+ data,
+ output_data.quadrature_points);
+ internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
+ DataSetDescriptor::cell (),
+ data);
const UpdateFlags update_flags(data.current_update_flags());
const std::vector<double> &weights=quadrature.get_weights();
{
namespace
{
+ /**
+ * Depending on what information is called for in the update flags of the
+ * @p data object, compute the various pieces of information that is required
+ * by the fill_fe_face_values() and fill_fe_subface_values() functions.
+ * This function simply unifies the work that would be done by
+ * those two functions.
+ *
+ * The resulting data is put into the @p output_data argument.
+ */
template <int dim, int spacedim>
void
- compute_fill_face (const dealii::MappingQ1<dim,spacedim> &mapping,
- const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int subface_no,
- const unsigned int n_q_points,
- const std::vector<double> &weights,
- const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
- FEValuesData<dim,spacedim> &output_data)
+ maybe_compute_face_data (const dealii::MappingQ1<dim,spacedim> &mapping,
+ const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const unsigned int n_q_points,
+ const std::vector<double> &weights,
+ const typename dealii::MappingQ1<dim,spacedim>::InternalData &data,
+ FEValuesData<dim,spacedim> &output_data)
{
const UpdateFlags update_flags(data.current_update_flags());
void
MappingQ1<dim,spacedim>::
compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int subface_no,
- const unsigned int n_q_points,
- const DataSetDescriptor data_set,
- const std::vector<double> &weights,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const DataSetDescriptor data_set,
+ const Quadrature<dim-1> &quadrature,
const InternalData &internal_data,
FEValuesData<dim,spacedim> &output_data) const
{
- compute_fill (cell, n_q_points, data_set, CellSimilarity::none,
- internal_data,
- output_data.quadrature_points);
- internal::compute_fill_face (*this,
- cell, face_no, subface_no, n_q_points,
- weights, internal_data,
- output_data);
+ // if necessary, recompute the support points of the transformation of this cell
+ // (note that we need to first check the triangulation pointer, since otherwise
+ // the second test might trigger an exception if the triangulations are not the
+ // same)
+ if ((internal_data.mapping_support_points.size() == 0)
+ ||
+ (&cell->get_triangulation() !=
+ &internal_data.cell_of_current_support_points->get_triangulation())
+ ||
+ (cell != internal_data.cell_of_current_support_points))
+ {
+ compute_mapping_support_points(cell, internal_data.mapping_support_points);
+ internal_data.cell_of_current_support_points = cell;
+ }
+
+ internal::maybe_compute_q_points<dim,spacedim> (data_set,
+ internal_data,
+ output_data.quadrature_points);
+ internal::maybe_update_Jacobians<dim,spacedim> (CellSimilarity::none,
+ data_set,
+ internal_data);
+ internal::maybe_compute_face_data (*this,
+ cell, face_no, subface_no, quadrature.size(),
+ quadrature.get_weights(), internal_data,
+ output_data);
}
const unsigned int n_q_points = quadrature.size();
compute_fill_face (cell, face_no, numbers::invalid_unsigned_int,
- n_q_points,
DataSetDescriptor::face (face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
n_q_points),
- quadrature.get_weights(),
+ quadrature,
data,
output_data);
}
const unsigned int n_q_points = quadrature.size();
compute_fill_face (cell, face_no, subface_no,
- n_q_points,
DataSetDescriptor::subface (face_no, subface_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
n_q_points,
cell->subface_case(face_no)),
- quadrature.get_weights(),
+ quadrature,
data,
output_data);
}