const SubQuadrature &quadrature,
const unsigned int face_no);
+ /**
+ * Compute the cell quadrature formula corresponding to using
+ * <tt>quadrature</tt> on face <tt>face_no</tt> taking into account the
+ * orientation of the face. For further details, see the general doc for this
+ * class.
+ */
+ static Quadrature<dim>
+ project_to_oriented_face(const ReferenceCell &reference_cell,
+ const SubQuadrature &quadrature,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation);
+
/**
* Compute the quadrature points on the cell if the given quadrature formula
* is used on face <tt>face_no</tt>, subface number <tt>subface_no</tt>
const RefinementCase<dim - 1> &ref_case =
RefinementCase<dim - 1>::isotropic_refinement);
+ /**
+ * Compute the cell quadrature formula corresponding to using
+ * <tt>quadrature</tt> on subface <tt>subface_no</tt> of face
+ * <tt>face_no</tt> with SubfaceCase<dim> <tt>ref_case</tt>. The last
+ * argument is only used in 3d.
+ *
+ * @note Only the points are transformed. The quadrature weights are the
+ * same as those of the original rule.
+ */
+ static Quadrature<dim>
+ project_to_oriented_subface(const ReferenceCell & reference_cell,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const internal::SubfaceCase<dim> ref_case);
+
/**
* Take a collection of face quadrature formulas and generate a cell
* quadrature formula from it where the quadrature points of the given
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const;
+ /**
+ * As opposed to the fill_fe_face_values()
+ * function that relies on pre-computed information of InternalDataBase, this
+ * function chooses the flexible evaluation path on the cell and points
+ * passed in to the current function.
+ *
+ * @param[in] cell The cell where to evaluate the mapping.
+ *
+ * @param[in] face_number The face number where to evaluate the mapping.
+ *
+ * @param[in] face_quadrature The quadrature points where the
+ * transformation (Jacobians, positions) should be computed.
+ *
+ * @param[in] internal_data A reference to an object previously created
+ * that may be used to store information the mapping can compute once on the
+ * reference cell. See the documentation of the Mapping::InternalDataBase
+ * class for an extensive description of the purpose of these objects.
+ *
+ * @param[out] output_data A struct containing the evaluated quantities such
+ * as the Jacobian resulting from application of the mapping on the given
+ * cell with its underlying manifolds.
+ */
+ void
+ fill_mapping_data_for_face_quadrature(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_number,
+ const Quadrature<dim - 1> & face_quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const;
+
/**
* @name Interface with FEValues and friends
* @{
virtual std::size_t
memory_consumption() const override;
+ /**
+ * Location of quadrature points of faces or subfaces in 3d with all
+ * possible orientations. Can be accessed with the correct offset provided
+ * via QProjector::DataSetDescriptor. Not needed/used for cells.
+ */
+ AlignedVector<Point<dim>> quadrature_points;
+
/**
* Values of shape functions. Access by function @p shape.
*
const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
std::vector<Point<spacedim>> & quadrature_points,
+ std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
+ std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians,
std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
{
const UpdateFlags update_flags = data.update_each;
data.volume_elements[point] =
data.contravariant[point].determinant();
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_jacobians)
+ {
+ const unsigned int n_q_points = data.contravariant.size();
+ AssertDimension(jacobians.size(), n_q_points);
+ if (cell_similarity != CellSimilarity::translation)
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ jacobians[point] = data.contravariant[point];
+ }
+
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_inverse_jacobians)
+ {
+ const unsigned int n_q_points = data.contravariant.size();
+ AssertDimension(inverse_jacobians.size(), n_q_points);
+ if (cell_similarity != CellSimilarity::translation)
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ inverse_jacobians[point] = data.covariant[point].transpose();
+ }
+
if (evaluation_flag & EvaluationFlags::hessians)
{
constexpr int desymmetrize_3d[6][2] = {
}
- /**
- * 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>
inline void
- maybe_compute_q_points(
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ maybe_update_q_points_Jacobians_generic(
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
- std::vector<Point<spacedim>> &quadrature_points)
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const unsigned int polynomial_degree,
+ const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
+ std::vector<Point<spacedim>> & quadrature_points,
+ std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
+ std::vector<DerivativeForm<1, spacedim, dim>> &inverse_jacobians)
{
- const UpdateFlags update_flags = data.update_each;
+ const UpdateFlags update_flags = data.update_each;
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
- if (update_flags & update_quadrature_points)
- for (unsigned int point = 0; point < quadrature_points.size(); ++point)
+ const unsigned int n_points = unit_points.size();
+ const unsigned int n_lanes = VectorizedArray<double>::size();
+
+ // Use the more heavy VectorizedArray code path if there is more than
+ // one point left to compute
+ for (unsigned int i = 0; i < n_points; i += n_lanes)
+ if (n_points - i > 1)
{
- 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;
- }
- }
+ Point<dim, VectorizedArray<double>> p_vec;
+ for (unsigned int j = 0; j < n_lanes; ++j)
+ if (i + j < n_points)
+ for (unsigned int d = 0; d < dim; ++d)
+ p_vec[d][j] = unit_points[i + j][d];
+ else
+ for (unsigned int d = 0; d < dim; ++d)
+ p_vec[d][j] = unit_points[i][d];
+ const auto result =
+ internal::evaluate_tensor_product_value_and_gradient(
+ polynomials_1d,
+ support_points,
+ p_vec,
+ polynomial_degree == 1,
+ renumber_lexicographic_to_hierarchic);
+ if (update_flags & update_quadrature_points)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ quadrature_points[i + j][d] = result.first[d][j];
- /**
- * 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>
- inline void
- maybe_update_Jacobians(
- const CellSimilarity::Similarity cell_similarity,
- const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
- const typename dealii::MappingQ<dim, spacedim>::InternalData &data)
- {
- const UpdateFlags update_flags = data.update_each;
+ if (cell_similarity == CellSimilarity::translation)
+ continue;
- 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)
+ if (update_flags & update_contravariant_transformation)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ data.contravariant[i + j][d][e] = result.second[e][d][j];
+
+ if (update_flags & update_volume_elements)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ data.volume_elements[i + j] =
+ data.contravariant[i + j].determinant();
+
+ if (update_flags & update_jacobians)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ jacobians[i + j] = data.contravariant[i + j];
+
+ if (update_flags & update_covariant_transformation)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ data.covariant[i + j] =
+ data.contravariant[i + j].covariant_form();
+
+ if (update_flags & update_inverse_jacobians)
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ inverse_jacobians[i + j] = data.covariant[i + j].transpose();
+ }
+ else
{
- const unsigned int n_q_points = data.contravariant.size();
+ const auto result =
+ internal::evaluate_tensor_product_value_and_gradient(
+ polynomials_1d,
+ support_points,
+ unit_points[i],
+ polynomial_degree == 1,
+ renumber_lexicographic_to_hierarchic);
- std::fill(data.contravariant.begin(),
- data.contravariant.end(),
- DerivativeForm<1, dim, spacedim>());
+ if (update_flags & update_quadrature_points)
+ quadrature_points[i] = result.first;
- Assert(data.n_shape_functions > 0, ExcInternalError());
+ if (cell_similarity == CellSimilarity::translation)
+ continue;
- for (unsigned int point = 0; point < n_q_points; ++point)
+ if (update_flags & update_contravariant_transformation)
{
- 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.derivative(point + data_set, 0)[j] *
- data.mapping_support_points[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.derivative(point + data_set, k)[j] *
- data.mapping_support_points[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];
+ DerivativeForm<1, spacedim, dim> jac_transposed = result.second;
+ data.contravariant[i] = jac_transposed.transpose();
}
- }
- 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)
+ data.volume_elements[i] = data.contravariant[i].determinant();
- 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();
+ if (update_flags & update_jacobians)
+ jacobians[i] = data.contravariant[i];
+
+ if (update_flags & update_covariant_transformation)
+ data.covariant[i] = data.contravariant[i].covariant_form();
+
+ if (update_flags & update_inverse_jacobians)
+ inverse_jacobians[i] = data.covariant[i].transpose();
}
}
output_data.normal_vectors[i] =
Point<spacedim>(output_data.boundary_forms[i] /
output_data.boundary_forms[i].norm());
-
- if (update_flags & update_jacobians)
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.jacobians[point] = data.contravariant[point];
-
- if (update_flags & update_inverse_jacobians)
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.inverse_jacobians[point] =
- data.covariant[point].transpose();
}
}
+
/**
* Do the work of MappingQ::fill_fe_face_values() and
* MappingQ::fill_fe_subface_values() in a generic way,
const typename QProjector<dim>::DataSetDescriptor data_set,
const Quadrature<dim - 1> & quadrature,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const unsigned int polynomial_degree,
+ const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data)
{
CellSimilarity::none,
data,
output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians,
output_data.jacobian_grads);
}
else
{
- maybe_compute_q_points<dim, spacedim>(data_set,
- data,
- output_data.quadrature_points);
- maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
- data_set,
- data);
+ internal::MappingQImplementation::
+ maybe_update_q_points_Jacobians_generic(
+ CellSimilarity::none,
+ data,
+ make_array_view(&data.quadrature_points[data_set],
+ &data.quadrature_points[data_set] +
+ quadrature.size()),
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
+ output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians);
maybe_update_jacobian_grads<dim, spacedim>(
CellSimilarity::none, data_set, data, output_data.jacobian_grads);
}
* class.
*/
void
- compute_mapping_data_for_generic_points(
+ compute_mapping_data_for_quadrature(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> & unit_points,
- MappingData & mapping_data);
+ CellSimilarity::Similarity &cell_similarity,
+ const Quadrature<dim> & quadrature,
+ MappingData & mapping_data);
/**
* Compute the mapping related data for the given @p mapping,
*/
std::vector<unsigned int> unit_points_index;
+ /**
+ * A pointer to the internal data of the underlying mapping.
+ */
+ std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+ internal_mapping_data;
+
/**
* A pointer to the underlying mapping.
*/
// always save quadrature points for now
update_flags_mapping |= update_quadrature_points;
+
+ // construct internal_mapping_data for MappingQ to be able to reuse it in
+ // reinit() calls to avoid memory allocations
+ if (const MappingQ<dim, spacedim> *mapping_q =
+ dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+ {
+ internal_mapping_data =
+ std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
+ mapping_q->get_degree());
+ }
}
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const std::vector<Point<dim>> & unit_points_in)
{
- reinit(cell, make_array_view(unit_points_in.begin(), unit_points_in.end()));
+ reinit(cell, Quadrature<dim>(unit_points_in));
}
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const ArrayView<const Point<dim>> & unit_points_in)
{
- unit_points =
- std::vector<Point<dim>>(unit_points_in.begin(), unit_points_in.end());
-
- mapping_data.resize(1);
- compute_mapping_data_for_generic_points(cell,
- unit_points_in,
- mapping_data[0]);
-
- state = State::single_cell;
- is_reinitialized();
+ reinit(cell,
+ std::vector<Point<dim>>(unit_points_in.begin(),
+ unit_points_in.end()));
}
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Quadrature<dim> & quadrature)
{
- const auto &points = quadrature.get_points();
- const auto &weights = quadrature.get_weights();
+ unit_points = quadrature.get_points();
- reinit(cell, points);
+ mapping_data.resize(1);
+ CellSimilarity::Similarity cell_similarity =
+ CellSimilarity::Similarity::none;
+ compute_mapping_data_for_quadrature(cell,
+ cell_similarity,
+ quadrature,
+ mapping_data[0]);
- if (update_flags_mapping & update_JxW_values)
- for (unsigned int q = 0; q < points.size(); ++q)
- mapping_data[0].JxW_values[q] =
- determinant(Tensor<2, dim>(mapping_data[0].jacobians[q])) *
- weights[q];
+ state = State::single_cell;
+ is_reinitialized();
}
const IteratorRange<Iterator> & cell_iterator_range,
const std::vector<std::vector<Point<dim>>> &unit_points_vector,
const unsigned int n_unfiltered_cells)
+ {
+ const unsigned int n_cells = unit_points_vector.size();
+ AssertDimension(n_cells,
+ std::distance(cell_iterator_range.begin(),
+ cell_iterator_range.end()));
+
+ std::vector<Quadrature<dim>> quadrature_vector(n_cells);
+ for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
+ quadrature_vector[cell_index] =
+ Quadrature<dim>(quadrature_vector[cell_index].get_points());
+
+ reinit_cells(cell_iterator_range, quadrature_vector, n_unfiltered_cells);
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <typename Iterator>
+ void
+ MappingInfo<dim, spacedim>::reinit_cells(
+ const IteratorRange<Iterator> & cell_iterator_range,
+ const std::vector<Quadrature<dim>> &quadrature_vector,
+ const unsigned int n_unfiltered_cells)
{
do_cell_index_compression =
n_unfiltered_cells != numbers::invalid_unsigned_int;
- const unsigned int n_cells = unit_points_vector.size();
+ const unsigned int n_cells = quadrature_vector.size();
AssertDimension(n_cells,
std::distance(cell_iterator_range.begin(),
cell_iterator_range.end()));
// fill unit points index offset vector
unit_points_index.reserve(n_cells + 1);
unit_points_index.push_back(0);
- for (const auto &unit_points : unit_points_vector)
- unit_points_index.push_back(unit_points_index.back() +
- unit_points.size());
+ for (const auto &quadrature : quadrature_vector)
+ unit_points_index.push_back(unit_points_index.back() + quadrature.size());
const unsigned int n_unit_points = unit_points_index.back();
if (do_cell_index_compression)
cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
numbers::invalid_unsigned_int);
+ CellSimilarity::Similarity cell_similarity =
+ CellSimilarity::Similarity::none;
unsigned int cell_index = 0;
for (const auto &cell : cell_iterator_range)
{
auto it = unit_points.begin() + unit_points_index[cell_index];
- for (const auto &unit_point : unit_points_vector[cell_index])
+ for (const auto &unit_point :
+ quadrature_vector[cell_index].get_points())
{
*it = unit_point;
++it;
}
- compute_mapping_data_for_generic_points(cell,
- unit_points_vector[cell_index],
- mapping_data[cell_index]);
+ compute_mapping_data_for_quadrature(cell,
+ cell_similarity,
+ quadrature_vector[cell_index],
+ mapping_data[cell_index]);
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
- template <int dim, int spacedim>
- template <typename Iterator>
- void
- MappingInfo<dim, spacedim>::reinit_cells(
- const IteratorRange<Iterator> & cell_iterator_range,
- const std::vector<Quadrature<dim>> &quadrature_vector,
- const unsigned int n_unfiltered_cells)
- {
- const unsigned int n_cells = quadrature_vector.size();
- AssertDimension(n_cells,
- std::distance(cell_iterator_range.begin(),
- cell_iterator_range.end()));
-
- std::vector<std::vector<Point<dim>>> unit_points_vector(n_cells);
- for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
- unit_points_vector[cell_index] = std::vector<Point<dim>>(
- quadrature_vector[cell_index].get_points().begin(),
- quadrature_vector[cell_index].get_points().end());
-
- reinit_cells(cell_iterator_range, unit_points_vector, n_unfiltered_cells);
-
- if (update_flags_mapping & update_JxW_values)
- for (unsigned int cell_index = 0; cell_index < n_cells; ++cell_index)
- {
- const auto &weights = quadrature_vector[cell_index].get_weights();
- for (unsigned int q = 0; q < weights.size(); ++q)
- mapping_data[cell_index].JxW_values[q] =
- determinant(
- Tensor<2, dim>(mapping_data[cell_index].jacobians[q])) *
- weights[q];
- }
- }
-
-
-
template <int dim, int spacedim>
template <typename Iterator>
void
template <int dim, int spacedim>
void
- MappingInfo<dim, spacedim>::compute_mapping_data_for_generic_points(
+ MappingInfo<dim, spacedim>::compute_mapping_data_for_quadrature(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> & unit_points,
+ CellSimilarity::Similarity & cell_similarity,
+ const Quadrature<dim> & quadrature,
MappingData & mapping_data)
{
+ update_flags_mapping |=
+ mapping->requires_update_flags(update_flags_mapping);
+
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
if (const MappingQ<dim, spacedim> *mapping_q =
dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
{
- mapping_q->fill_mapping_data_for_generic_points(cell,
- unit_points,
- update_flags_mapping,
- mapping_data);
- }
- else if (const MappingCartesian<dim, spacedim> *mapping_cartesian =
- dynamic_cast<const MappingCartesian<dim, spacedim> *>(
- &(*mapping)))
- {
- mapping_cartesian->fill_mapping_data_for_generic_points(
- cell, unit_points, update_flags_mapping, mapping_data);
+ (void)mapping_q;
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize(update_flags_mapping, quadrature, quadrature.size());
}
else
{
- FE_DGQ<dim, spacedim> dummy_fe(1);
- dealii::FEValues<dim, spacedim> fe_values(
- *mapping,
- dummy_fe,
- Quadrature<dim>(
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
- update_flags_mapping);
- fe_values.reinit(cell);
- mapping_data.initialize(unit_points.size(), update_flags_mapping);
- if (update_flags_mapping & update_jacobians)
- for (unsigned int q = 0; q < unit_points.size(); ++q)
- mapping_data.jacobians[q] = fe_values.jacobian(q);
- if (update_flags_mapping & update_inverse_jacobians)
- for (unsigned int q = 0; q < unit_points.size(); ++q)
- mapping_data.inverse_jacobians[q] = fe_values.inverse_jacobian(q);
- if (update_flags_mapping & update_quadrature_points)
- for (unsigned int q = 0; q < unit_points.size(); ++q)
- mapping_data.quadrature_points[q] = fe_values.quadrature_point(q);
+ internal_mapping_data =
+ mapping->get_data(update_flags_mapping, quadrature);
}
+
+ cell_similarity = mapping->fill_fe_values(
+ cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data);
}
update_flags_mapping |=
mapping->requires_update_flags(update_flags_mapping);
- mapping_data.initialize(quadrature.get_points().size(),
- update_flags_mapping);
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
- auto internal_mapping_data =
- mapping->get_data(update_flags_mapping, quadrature);
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
+ if (const MappingQ<dim, spacedim> *mapping_q =
+ dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ {
+ (void)mapping_q;
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize(update_flags_mapping, quadrature, quadrature.size());
+ }
+ else
+ {
+ internal_mapping_data =
+ mapping->get_data(update_flags_mapping, quadrature);
+ }
mapping->fill_fe_immersed_surface_values(cell,
quadrature,
update_flags_mapping |=
mapping->requires_update_flags(update_flags_mapping);
- mapping_data.initialize(quadrature.get_points().size(),
- update_flags_mapping);
-
- auto internal_mapping_data =
- mapping->get_face_data(update_flags_mapping,
- hp::QCollection<dim - 1>(quadrature));
+ mapping_data.initialize(quadrature.size(), update_flags_mapping);
- mapping->fill_fe_face_values(cell,
- face_no,
- hp::QCollection<dim - 1>(quadrature),
- *internal_mapping_data,
- mapping_data);
+ // reuse internal_mapping_data for MappingQ to avoid memory allocations
+ if (const MappingQ<dim, spacedim> *mapping_q =
+ dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+ {
+ auto &data =
+ dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+ *internal_mapping_data);
+ data.initialize_face(update_flags_mapping,
+ QProjector<dim>::project_to_oriented_face(
+ ReferenceCells::get_hypercube<dim>(),
+ quadrature,
+ face_no,
+ cell->face_orientation(face_no),
+ cell->face_flip(face_no),
+ cell->face_rotation(face_no)),
+ quadrature.size());
+
+ mapping_q->fill_mapping_data_for_face_quadrature(
+ cell, face_no, quadrature, *internal_mapping_data, mapping_data);
+ }
+ else
+ {
+ auto internal_mapping_data =
+ mapping->get_face_data(update_flags_mapping,
+ hp::QCollection<dim - 1>(quadrature));
+
+ mapping->fill_fe_face_values(cell,
+ face_no,
+ hp::QCollection<dim - 1>(quadrature),
+ *internal_mapping_data,
+ mapping_data);
+ }
}
} // namespace NonMatching
q_points.push_back(cell_point);
}
}
+
+ std::vector<Point<2>>
+ mutate_points_with_offset(const std::vector<Point<2>> &points,
+ const unsigned int offset)
+ {
+ switch (offset)
+ {
+ case 0:
+ return points;
+ case 1:
+ case 2:
+ case 3:
+ return rotate(points, offset);
+ case 4:
+ return reflect(points);
+ case 5:
+ case 6:
+ case 7:
+ return rotate(reflect(points), 8 - offset);
+ default:
+ Assert(false, ExcInternalError());
+ }
+ return {};
+ }
+
+ Quadrature<2>
+ mutate_quadrature(const Quadrature<2> &quadrature,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+ {
+ static const unsigned int offset[2][2][2] = {
+ {{4, 5}, // face_orientation=false; face_flip=false;
+ // face_rotation=false and true
+ {6, 7}}, // face_orientation=false; face_flip=true;
+ // face_rotation=false and true
+ {{0, 1}, // face_orientation=true; face_flip=false;
+ // face_rotation=false and true
+ {2, 3}}}; // face_orientation=true; face_flip=true;
+ // face_rotation=false and true
+
+ return Quadrature<2>(
+ mutate_points_with_offset(
+ quadrature.get_points(),
+ offset[face_orientation][face_flip][face_rotation]),
+ quadrature.get_weights());
+ }
+
+ std::pair<unsigned int, RefinementCase<2>>
+ select_subface_no_and_refinement_case(
+ const unsigned int subface_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const internal::SubfaceCase<3> ref_case)
+ {
+ constexpr int dim = 3;
+ // for each subface of a given FaceRefineCase
+ // there is a corresponding equivalent
+ // subface number of one of the "standard"
+ // RefineCases (cut_x, cut_y, cut_xy). Map
+ // the given values to those equivalent
+ // ones.
+
+ // first, define an invalid number
+ static const unsigned int e = numbers::invalid_unsigned_int;
+
+ static const RefinementCase<dim - 1>
+ equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
+ [GeometryInfo<3>::max_children_per_face] = {
+ // case_none. there should be only
+ // invalid values here. However, as
+ // this function is also called (in
+ // tests) for cells which have no
+ // refined faces, use isotropic
+ // refinement instead
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy},
+ // case_x
+ {RefinementCase<dim - 1>::cut_x,
+ RefinementCase<dim - 1>::cut_x,
+ RefinementCase<dim - 1>::no_refinement,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_x1y
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_x,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_x2y
+ {RefinementCase<dim - 1>::cut_x,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_x1y2y
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy},
+ // case_y
+ {RefinementCase<dim - 1>::cut_y,
+ RefinementCase<dim - 1>::cut_y,
+ RefinementCase<dim - 1>::no_refinement,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_y1x
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_y,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_y2x
+ {RefinementCase<dim - 1>::cut_y,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::no_refinement},
+ // case_y1x2x
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy},
+ // case_xy (case_isotropic)
+ {RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy,
+ RefinementCase<dim - 1>::cut_xy}};
+
+ static const unsigned int
+ equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic +
+ 1][GeometryInfo<3>::max_children_per_face] =
+ {// case_none, see above
+ {0, 1, 2, 3},
+ // case_x
+ {0, 1, e, e},
+ // case_x1y
+ {0, 2, 1, e},
+ // case_x2y
+ {0, 1, 3, e},
+ // case_x1y2y
+ {0, 2, 1, 3},
+ // case_y
+ {0, 1, e, e},
+ // case_y1x
+ {0, 1, 1, e},
+ // case_y2x
+ {0, 2, 3, e},
+ // case_y1x2x
+ {0, 1, 2, 3},
+ // case_xy (case_isotropic)
+ {0, 1, 2, 3}};
+
+ // If face-orientation or face_rotation are
+ // non-standard, cut_x and cut_y have to be
+ // exchanged.
+ static const RefinementCase<dim - 1> ref_case_permutation[4] = {
+ RefinementCase<dim - 1>::no_refinement,
+ RefinementCase<dim - 1>::cut_y,
+ RefinementCase<dim - 1>::cut_x,
+ RefinementCase<dim - 1>::cut_xy};
+
+ // set a corresponding (equivalent)
+ // RefineCase and subface number
+ const RefinementCase<dim - 1> equ_ref_case =
+ equivalent_refine_case[ref_case][subface_no];
+ const unsigned int equ_subface_no =
+ equivalent_subface_number[ref_case][subface_no];
+ // make sure, that we got a valid subface and RefineCase
+ Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
+ ExcInternalError());
+ Assert(equ_subface_no != e, ExcInternalError());
+ // now, finally respect non-standard faces
+ const RefinementCase<dim - 1> final_ref_case =
+ (face_orientation == face_rotation ?
+ ref_case_permutation[equ_ref_case] :
+ equ_ref_case);
+
+ const unsigned int final_subface_no =
+ GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
+ final_ref_case),
+ 4,
+ equ_subface_no,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ equ_ref_case);
+
+ return std::make_pair(final_subface_no, final_ref_case);
+ }
} // namespace
} // namespace QProjector
} // namespace internal
}
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_oriented_face(const ReferenceCell &reference_cell,
+ const Quadrature<dim - 1> &quadrature,
+ const unsigned int face_no,
+ const bool,
+ const bool,
+ const bool)
+{
+ return QProjector<dim>::project_to_face(reference_cell, quadrature, face_no);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_oriented_face(const ReferenceCell &reference_cell,
+ const Quadrature<2> &quadrature,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
+
+ const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
+ quadrature, face_orientation, face_flip, face_rotation);
+
+ return QProjector<3>::project_to_face(reference_cell, mutation, face_no);
+}
+
+
template <>
void
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_oriented_subface(
+ const ReferenceCell & reference_cell,
+ const Quadrature<dim - 1> &quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const bool,
+ const bool,
+ const bool,
+ const internal::SubfaceCase<dim>)
+{
+ return QProjector<dim>::project_to_subface(
+ reference_cell,
+ quadrature,
+ face_no,
+ subface_no,
+ RefinementCase<dim - 1>::isotropic_refinement);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_oriented_subface(
+ const ReferenceCell & reference_cell,
+ const Quadrature<2> & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const internal::SubfaceCase<3> ref_case)
+{
+ Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
+
+ const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
+ quadrature, face_orientation, face_flip, face_rotation);
+
+ const std::pair<unsigned int, RefinementCase<2>>
+ final_subface_no_and_ref_case =
+ internal::QProjector::select_subface_no_and_refinement_case(
+ subface_no, face_orientation, face_flip, face_rotation, ref_case);
+
+ return QProjector<3>::project_to_subface(
+ reference_cell,
+ mutation,
+ face_no,
+ final_subface_no_and_ref_case.first,
+ final_subface_no_and_ref_case.second);
+}
+
+
+
template <>
Quadrature<1>
QProjector<1>::project_to_all_faces(const ReferenceCell & reference_cell,
std::vector<double> weights;
weights.reserve(n_points_total);
- Table<2, std::vector<Point<2>>> mutations(quadrature.size(), 8);
-
- // do the following for all possible mutations of a face (mutation==0
- // corresponds to a face with standard orientation, no flip and no
- // rotation)
- for (unsigned int face = 0; face < quadrature.size(); ++face)
- {
- const std::vector<Point<2>> &quad = quadrature[face].get_points();
- mutations(face, 0) = quad;
- mutations(face, 1) = internal::QProjector::rotate(quad, 1);
- mutations(face, 2) = internal::QProjector::rotate(quad, 2);
- mutations(face, 3) = internal::QProjector::rotate(quad, 3);
- mutations(face, 4) = internal::QProjector::reflect(quad);
- mutations(face, 5) =
- internal::QProjector::rotate(mutations(face, 4), 3);
- mutations(face, 6) =
- internal::QProjector::rotate(mutations(face, 4), 2);
- mutations(face, 7) =
- internal::QProjector::rotate(mutations(face, 4), 1);
- }
-
- for (unsigned int i = 0; i < 8; ++i)
+ for (unsigned int offset = 0; offset < 8; ++offset)
{
+ const auto mutation = internal::QProjector::mutate_points_with_offset(
+ quadrature[0].get_points(), offset);
// project to each face and append results
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
const unsigned int q_index = quadrature.size() == 1 ? 0 : face;
internal::QProjector::project_to_hex_face_and_append(
- mutations(q_index, i), face, q_points);
+ q_index > 0 ? internal::QProjector::mutate_points_with_offset(
+ quadrature[face].get_points(), offset) :
+ mutation,
+ face,
+ q_points);
std::copy(quadrature[q_index].get_weights().begin(),
quadrature[q_index].get_weights().end(),
}
}
-
Assert(q_points.size() == n_points_total, ExcInternalError());
Assert(weights.size() == n_points_total, ExcInternalError());
Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented());
- const unsigned int dim = 3;
- const std::vector<Point<2>> &quad = quadrature.get_points();
- std::vector<Point<2>> q_reflected = internal::QProjector::reflect(quad);
- std::array<std::vector<Point<2>>, 8> mutations{
- {quad,
- internal::QProjector::rotate(quad, 1),
- internal::QProjector::rotate(quad, 2),
- internal::QProjector::rotate(quad, 3),
- q_reflected,
- internal::QProjector::rotate(q_reflected, 3),
- internal::QProjector::rotate(q_reflected, 2),
- internal::QProjector::rotate(q_reflected, 1)}};
+ const unsigned int dim = 3;
const unsigned int n_points = quadrature.size(),
n_faces = GeometryInfo<dim>::faces_per_cell,
// do the following for all possible mutations of a face (mutation==0
// corresponds to a face with standard orientation, no flip and no rotation)
- for (const auto &mutation : mutations)
+ for (unsigned int offset = 0; offset < 8; ++offset)
{
+ const auto mutation =
+ internal::QProjector::mutate_points_with_offset(quadrature.get_points(),
+ offset);
+
// project to each face and copy results
for (unsigned int face = 0; face < n_faces; ++face)
for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
0 // cut_xy
};
-
- // for each subface of a given FaceRefineCase
- // there is a corresponding equivalent
- // subface number of one of the "standard"
- // RefineCases (cut_x, cut_y, cut_xy). Map
- // the given values to those equivalent
- // ones.
-
- // first, define an invalid number
- static const unsigned int e = numbers::invalid_unsigned_int;
-
- static const RefinementCase<dim - 1>
- equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
- [GeometryInfo<3>::max_children_per_face] = {
- // case_none. there should be only
- // invalid values here. However, as
- // this function is also called (in
- // tests) for cells which have no
- // refined faces, use isotropic
- // refinement instead
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy},
- // case_x
- {RefinementCase<dim - 1>::cut_x,
- RefinementCase<dim - 1>::cut_x,
- RefinementCase<dim - 1>::no_refinement,
- RefinementCase<dim - 1>::no_refinement},
- // case_x1y
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_x,
- RefinementCase<dim - 1>::no_refinement},
- // case_x2y
- {RefinementCase<dim - 1>::cut_x,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::no_refinement},
- // case_x1y2y
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy},
- // case_y
- {RefinementCase<dim - 1>::cut_y,
- RefinementCase<dim - 1>::cut_y,
- RefinementCase<dim - 1>::no_refinement,
- RefinementCase<dim - 1>::no_refinement},
- // case_y1x
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_y,
- RefinementCase<dim - 1>::no_refinement},
- // case_y2x
- {RefinementCase<dim - 1>::cut_y,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::no_refinement},
- // case_y1x2x
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy},
- // case_xy (case_isotropic)
- {RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy,
- RefinementCase<dim - 1>::cut_xy}};
-
- static const unsigned int
- equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
- [GeometryInfo<3>::max_children_per_face] = {
- // case_none, see above
- {0, 1, 2, 3},
- // case_x
- {0, 1, e, e},
- // case_x1y
- {0, 2, 1, e},
- // case_x2y
- {0, 1, 3, e},
- // case_x1y2y
- {0, 2, 1, 3},
- // case_y
- {0, 1, e, e},
- // case_y1x
- {0, 1, 1, e},
- // case_y2x
- {0, 2, 3, e},
- // case_y1x2x
- {0, 1, 2, 3},
- // case_xy (case_isotropic)
- {0, 1, 2, 3}};
-
- // If face-orientation or face_rotation are
- // non-standard, cut_x and cut_y have to be
- // exchanged.
- static const RefinementCase<dim - 1> ref_case_permutation[4] = {
- RefinementCase<dim - 1>::no_refinement,
- RefinementCase<dim - 1>::cut_y,
- RefinementCase<dim - 1>::cut_x,
- RefinementCase<dim - 1>::cut_xy};
-
- // set a corresponding (equivalent)
- // RefineCase and subface number
- const RefinementCase<dim - 1> equ_ref_case =
- equivalent_refine_case[ref_case][subface_no];
- const unsigned int equ_subface_no =
- equivalent_subface_number[ref_case][subface_no];
- // make sure, that we got a valid subface and RefineCase
- Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
- ExcInternalError());
- Assert(equ_subface_no != e, ExcInternalError());
- // now, finally respect non-standard faces
- const RefinementCase<dim - 1> final_ref_case =
- (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
- equ_ref_case);
-
- // what we have now is the number of
- // the subface in the natural
- // orientation of the *face*. what we
- // need to know is the number of the
- // subface concerning the standard face
- // orientation as seen from the *cell*.
-
- // this mapping is not trivial, but we
- // have done exactly this stuff in the
- // child_cell_on_face function. in
- // order to reduce the amount of code
- // as well as to make maintaining the
- // functionality easier we want to
- // reuse that information. So we note
- // that on the bottom face (face 4) of
- // a hex cell the local x and y
- // coordinates of the face and the cell
- // coincide, thus also the refinement
- // case of the face corresponds to the
- // refinement case of the cell
- // (ignoring cell refinement along the
- // z direction). Using this knowledge
- // we can (ab)use the
- // child_cell_on_face function to do
- // exactly the transformation we are in
- // need of now
- const unsigned int final_subface_no =
- GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(final_ref_case),
- 4,
- equ_subface_no,
- face_orientation,
- face_flip,
- face_rotation,
- equ_ref_case);
+ const std::pair<unsigned int, RefinementCase<2>>
+ final_subface_no_and_ref_case =
+ internal::QProjector::select_subface_no_and_refinement_case(
+ subface_no, face_orientation, face_flip, face_rotation, ref_case);
return (((face_no * total_subfaces_per_face +
- ref_case_offset[final_ref_case - 1] + final_subface_no) +
+ ref_case_offset[final_subface_no_and_ref_case.second - 1] +
+ final_subface_no_and_ref_case.first) +
orientation_offset[face_orientation][face_flip][face_rotation]) *
n_quadrature_points);
}
void
MappingQ<dim, spacedim>::InternalData::initialize(
const UpdateFlags update_flags,
- const Quadrature<dim> &q,
+ const Quadrature<dim> &quadrature,
const unsigned int n_original_q_points)
{
// store the flags in the internal data object so we can access them
// in fill_fe_*_values()
this->update_each = update_flags;
- const unsigned int n_q_points = q.size();
-
- const bool needs_higher_order_terms =
- this->update_each &
- (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
- update_jacobian_pushed_forward_2nd_derivatives |
- update_jacobian_3rd_derivatives |
- update_jacobian_pushed_forward_3rd_derivatives);
+ const unsigned int n_q_points = quadrature.size();
if (this->update_each & update_covariant_transformation)
covariant.resize(n_original_q_points);
if (this->update_each & update_volume_elements)
volume_elements.resize(n_original_q_points);
- tensor_product_quadrature = q.is_tensor_product();
+ tensor_product_quadrature = quadrature.is_tensor_product();
// use of MatrixFree only for higher order elements and with more than one
// point where tensor products do not make sense
if (tensor_product_quadrature)
{
const std::array<Quadrature<1>, dim> &quad_array =
- q.get_tensor_basis();
+ quadrature.get_tensor_basis();
for (unsigned int i = 1; i < dim && tensor_product_quadrature; ++i)
{
if (quad_array[i - 1].size() != quad_array[i].size())
// numbering manually (building an FE_Q<dim> is relatively
// expensive due to constraints)
const FE_DGQ<1> fe(polynomial_degree);
- shape_info.reinit(q.get_tensor_basis()[0], fe);
+ shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
shape_info.lexicographic_numbering =
FETools::lexicographic_to_hierarchic_numbering<dim>(
polynomial_degree);
- shape_info.n_q_points = q.size();
+ shape_info.n_q_points = n_q_points;
shape_info.dofs_per_component_on_cell =
Utilities::pow(polynomial_degree + 1, dim);
}
}
}
+ const bool needs_higher_order_terms =
+ this->update_each &
+ (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
+ update_jacobian_pushed_forward_2nd_derivatives |
+ update_jacobian_3rd_derivatives |
+ update_jacobian_pushed_forward_3rd_derivatives);
+
+ const bool needs_higher_order_terms_generic =
+ !tensor_product_quadrature &&
+ (needs_higher_order_terms || this->update_each & update_jacobian_grads);
+
// Only fill the big arrays on demand in case we cannot use the tensor
// product quadrature code path
- if (dim == 1 || !tensor_product_quadrature || needs_higher_order_terms)
+ if (dim == 1 || needs_higher_order_terms_generic || needs_higher_order_terms)
{
- // see if we need the (transformation) shape function values
- // and/or gradients and resize the necessary arrays
- if (this->update_each & update_quadrature_points)
- shape_values.resize(n_shape_functions * n_q_points);
-
- if (this->update_each &
- (update_covariant_transformation |
- update_contravariant_transformation | update_JxW_values |
- update_boundary_forms | update_normal_vectors | update_jacobians |
- update_jacobian_grads | update_inverse_jacobians |
- update_jacobian_pushed_forward_grads |
- update_jacobian_2nd_derivatives |
- update_jacobian_pushed_forward_2nd_derivatives |
- update_jacobian_3rd_derivatives |
- update_jacobian_pushed_forward_3rd_derivatives))
- shape_derivatives.resize(n_shape_functions * n_q_points);
+ // compute shapes and derivatives for codim1 (for
+ // do_transform_real_to_unit_cell_internal_codim1)
+ if (dim == (spacedim - 1))
+ {
+ // see if we need the (transformation) shape function values
+ // and/or gradients and resize the necessary arrays
+ if (this->update_each & update_quadrature_points)
+ shape_values.resize(n_shape_functions * n_q_points);
+ if (this->update_each & update_jacobians)
+ shape_derivatives.resize(n_shape_functions * n_q_points);
+ }
if (this->update_each &
(update_jacobian_grads | update_jacobian_pushed_forward_grads))
shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
// now also fill the various fields with their correct values
- compute_shape_function_values(q.get_points());
+ compute_shape_function_values(quadrature.get_points());
}
}
void
MappingQ<dim, spacedim>::InternalData::initialize_face(
const UpdateFlags update_flags,
- const Quadrature<dim> &q,
+ const Quadrature<dim> &quadrature,
const unsigned int n_original_q_points)
{
- initialize(update_flags, q, n_original_q_points);
+ initialize(update_flags, quadrature, n_original_q_points);
+
+ const unsigned int n_q_points = quadrature.size();
+ quadrature_points.resize(n_q_points);
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ quadrature_points[q] = quadrature.get_points()[q];
if (dim > 1 && tensor_product_quadrature)
{
constexpr unsigned int facedim = dim - 1;
const FE_DGQ<1> fe(polynomial_degree);
- shape_info.reinit(q.get_tensor_basis()[0], fe);
+ shape_info.reinit(quadrature.get_tensor_basis()[0], fe);
shape_info.lexicographic_numbering =
FETools::lexicographic_to_hierarchic_numbering<facedim>(
polynomial_degree);
computed_cell_similarity,
data,
output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians,
output_data.jacobian_grads);
}
else
{
- internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
- QProjector<dim>::DataSetDescriptor::cell(),
- data,
- output_data.quadrature_points);
-
- internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
+ internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
- data);
+ data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
+ output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians);
internal::MappingQImplementation::maybe_update_jacobian_grads<dim,
spacedim>(
}
}
-
-
- // copy values from InternalData to vector given by reference
- if (update_flags & update_jacobians)
- {
- AssertDimension(output_data.jacobians.size(), n_q_points);
- if (computed_cell_similarity != CellSimilarity::translation)
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.jacobians[point] = data.contravariant[point];
- }
-
- // copy values from InternalData to vector given by reference
- if (update_flags & update_inverse_jacobians)
- {
- AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
- if (computed_cell_similarity != CellSimilarity::translation)
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.inverse_jacobians[point] =
- data.covariant[point].transpose();
- }
-
return computed_cell_similarity;
}
quadrature[0].size()),
quadrature[0],
data,
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
output_data);
}
cell->subface_case(face_no)),
quadrature,
data,
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
output_data);
}
data.mapping_support_points = this->compute_mapping_support_points(cell);
data.cell_of_current_support_points = cell;
- internal::MappingQImplementation::maybe_compute_q_points<dim, spacedim>(
- QProjector<dim>::DataSetDescriptor::cell(),
+ internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
+ CellSimilarity::none,
data,
- output_data.quadrature_points);
-
- internal::MappingQImplementation::maybe_update_Jacobians<dim, spacedim>(
- CellSimilarity::none, QProjector<dim>::DataSetDescriptor::cell(), data);
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
+ output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians);
internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
CellSimilarity::none,
}
}
}
-
- // copy values from InternalData to vector given by reference
- if ((update_flags & update_jacobians) != 0u)
- {
- AssertDimension(output_data.jacobians.size(), n_q_points);
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.jacobians[point] = data.contravariant[point];
- }
-
- // copy values from InternalData to vector given by reference
- if ((update_flags & update_inverse_jacobians) != 0u)
- {
- AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
- for (unsigned int point = 0; point < n_q_points; ++point)
- output_data.inverse_jacobians[point] =
- data.covariant[point].transpose();
- }
}
ExcNotImplemented());
output_data.initialize(unit_points.size(), update_flags);
- const std::vector<Point<spacedim>> support_points =
- this->compute_mapping_support_points(cell);
- const unsigned int n_points = unit_points.size();
- const unsigned int n_lanes = VectorizedArray<double>::size();
+ auto internal_data =
+ this->get_data(update_flags,
+ Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
+ unit_points.end())));
+ const InternalData &data = static_cast<const InternalData &>(*internal_data);
+ data.mapping_support_points = this->compute_mapping_support_points(cell);
- // Use the more heavy VectorizedArray code path if there is more than
- // one point left to compute
- for (unsigned int i = 0; i < n_points; i += n_lanes)
- if (n_points - i > 1)
- {
- Point<dim, VectorizedArray<double>> p_vec;
- for (unsigned int j = 0; j < n_lanes; ++j)
- if (i + j < n_points)
- for (unsigned int d = 0; d < dim; ++d)
- p_vec[d][j] = unit_points[i + j][d];
- else
- for (unsigned int d = 0; d < dim; ++d)
- p_vec[d][j] = unit_points[i][d];
-
- const auto result =
- internal::evaluate_tensor_product_value_and_gradient(
- polynomials_1d,
- support_points,
- p_vec,
- polynomial_degree == 1,
- renumber_lexicographic_to_hierarchic);
-
- if (update_flags & update_quadrature_points)
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
- for (unsigned int d = 0; d < spacedim; ++d)
- output_data.quadrature_points[i + j][d] = result.first[d][j];
+ internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
+ CellSimilarity::none,
+ data,
+ unit_points,
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
+ output_data.quadrature_points,
+ output_data.jacobians,
+ output_data.inverse_jacobians);
+}
- if (update_flags & update_jacobians)
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
- for (unsigned int d = 0; d < spacedim; ++d)
- for (unsigned int e = 0; e < dim; ++e)
- output_data.jacobians[i + j][d][e] = result.second[e][d][j];
- if (update_flags & update_inverse_jacobians)
- {
- DerivativeForm<1, spacedim, dim, VectorizedArray<double>> jac(
- result.second);
- const DerivativeForm<1, spacedim, dim, VectorizedArray<double>>
- inv_jac = jac.covariant_form();
- for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int e = 0; e < spacedim; ++e)
- output_data.inverse_jacobians[i + j][d][e] = inv_jac[d][e][j];
- }
- }
- else
- {
- const auto result =
- internal::evaluate_tensor_product_value_and_gradient(
- polynomials_1d,
- support_points,
- unit_points[i],
- polynomial_degree == 1,
- renumber_lexicographic_to_hierarchic);
-
- if (update_flags & update_quadrature_points)
- output_data.quadrature_points[i] = result.first;
-
- if (update_flags & update_jacobians)
- {
- DerivativeForm<1, spacedim, dim> jac = result.second;
- output_data.jacobians[i] = jac.transpose();
- }
- if (update_flags & update_inverse_jacobians)
- {
- DerivativeForm<1, spacedim, dim> jac(result.second);
- DerivativeForm<1, spacedim, dim> inv_jac = jac.covariant_form();
- for (unsigned int d = 0; d < dim; ++d)
- for (unsigned int e = 0; e < spacedim; ++e)
- output_data.inverse_jacobians[i][d][e] = inv_jac[d][e];
- }
- }
+template <int dim, int spacedim>
+void
+MappingQ<dim, spacedim>::fill_mapping_data_for_face_quadrature(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & face_quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const
+{
+ if (face_quadrature.get_points().empty())
+ return;
+
+ // ensure that the following static_cast is really correct:
+ Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+ ExcInternalError());
+ const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+ data.mapping_support_points = this->compute_mapping_support_points(cell);
+
+ internal::MappingQImplementation::do_fill_fe_face_values(
+ *this,
+ cell,
+ face_no,
+ numbers::invalid_unsigned_int,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ face_quadrature,
+ data,
+ polynomials_1d,
+ polynomial_degree,
+ renumber_lexicographic_to_hierarchic,
+ output_data);
}
projects / dealii.git / commitdiff