template <int dim>
inline Point<dim>
do_transform_real_to_unit_cell_internal_codim1(
- const typename dealii::Triangulation<dim, dim + 1>::cell_iterator &cell,
- const Point<dim + 1> & p,
- const Point<dim> & initial_p_unit,
- typename dealii::MappingQ<dim, dim + 1>::InternalData &mdata)
+ const Point<dim + 1> & p,
+ const Point<dim> & initial_p_unit,
+ const std::vector<Point<dim + 1>> & points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> & renumber)
{
const unsigned int spacedim = dim + 1;
- const unsigned int n_shapes = mdata.shape_values.size();
- (void)n_shapes;
- Assert(n_shapes != 0, ExcInternalError());
- Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
- Assert(mdata.shape_second_derivatives.size() == n_shapes,
- ExcInternalError());
-
- std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
- Assert(points.size() == n_shapes, ExcInternalError());
-
- Point<spacedim> p_minus_F;
+ AssertDimension(points.size(),
+ Utilities::pow(polynomials_1d.size(), dim));
- Tensor<1, spacedim> DF[dim];
- Tensor<1, spacedim> D2F[dim][dim];
+ Point<dim> p_unit = initial_p_unit;
- Point<dim> p_unit = initial_p_unit;
- Point<dim> f;
- Tensor<2, dim> df;
+ const double eps = 1.e-12;
+ const unsigned int loop_limit = 10;
- // Evaluate first and second derivatives
- mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
+ unsigned int loop = 0;
+ double f_weighted_norm_square = 1.;
- for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
+ while (f_weighted_norm_square > eps * eps && loop++ < loop_limit)
{
- const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
- const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
- const Point<spacedim> &point_k = points[k];
-
+ const auto p_real =
+ internal::evaluate_tensor_product_value_and_gradient(
+ polynomials_1d,
+ points,
+ p_unit,
+ polynomials_1d.size() == 2,
+ renumber);
+ Tensor<1, spacedim> p_minus_F = p - p_real.first;
+ const DerivativeForm<1, spacedim, dim> DF = p_real.second;
+
+ const auto hessian = internal::evaluate_tensor_product_hessian(
+ polynomials_1d, points, p_unit, renumber);
+ Point<dim> f;
+ Tensor<2, dim> df;
for (unsigned int j = 0; j < dim; ++j)
{
- DF[j] += grad_phi_k[j] * point_k;
+ f[j] = DF[j] * p_minus_F;
for (unsigned int l = 0; l < dim; ++l)
- D2F[j][l] += hessian_k[j][l] * point_k;
+ df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F;
}
- }
-
- p_minus_F = p;
- p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
-
-
- for (unsigned int j = 0; j < dim; ++j)
- f[j] = DF[j] * p_minus_F;
-
- for (unsigned int j = 0; j < dim; ++j)
- {
- f[j] = DF[j] * p_minus_F;
- for (unsigned int l = 0; l < dim; ++l)
- df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
- }
-
- const double eps = 1.e-12 * cell->diameter();
- const unsigned int loop_limit = 10;
-
- unsigned int loop = 0;
-
- while (f.norm() > eps && loop++ < loop_limit)
- {
// Solve [df(x)]d=f(x)
const Tensor<1, dim> d =
invert(df) * static_cast<const Tensor<1, dim> &>(f);
+ f_weighted_norm_square = d.norm_square();
p_unit -= d;
-
- for (unsigned int j = 0; j < dim; ++j)
- {
- DF[j].clear();
- for (unsigned int l = 0; l < dim; ++l)
- D2F[j][l].clear();
- }
-
- mdata.compute_shape_function_values(
- std::vector<Point<dim>>(1, p_unit));
-
- for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
- {
- const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
- const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k);
- const Point<spacedim> &point_k = points[k];
-
- for (unsigned int j = 0; j < dim; ++j)
- {
- DF[j] += grad_phi_k[j] * point_k;
- for (unsigned int l = 0; l < dim; ++l)
- D2F[j][l] += hessian_k[j][l] * point_k;
- }
- }
-
- // TODO: implement a line search here in much the same way as for
- // the corresponding function above that does so for dim==spacedim
- p_minus_F = p;
- p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
-
- for (unsigned int j = 0; j < dim; ++j)
- {
- f[j] = DF[j] * p_minus_F;
- for (unsigned int l = 0; l < dim; ++l)
- df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
- }
}
constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes;
constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2;
+ // Since MappingQ::InternalData does not have separate arrays for the
+ // covariant and contravariant transformations, but uses the arrays in
+ // the `MappingRelatedData`, it can happen that vectors do not have the
+ // right size
+ if (update_flags & update_contravariant_transformation)
+ jacobians.resize(n_q_points);
+ if (update_flags & update_covariant_transformation)
+ inverse_jacobians.resize(n_q_points);
+
EvaluationFlags::EvaluationFlags evaluation_flag =
(update_flags & update_quadrature_points ? EvaluationFlags::values :
EvaluationFlags::nothing) |
Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
data.n_shape_functions > 0,
ExcInternalError());
- Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
- n_q_points == data.contravariant.size(),
- ExcDimensionMismatch(n_q_points, data.contravariant.size()));
Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
n_q_points == jacobian_grads.size(),
ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
if (evaluation_flag & EvaluationFlags::gradients)
{
- std::fill(data.contravariant.begin(),
- data.contravariant.end(),
- DerivativeForm<1, dim, spacedim>());
// We need to reinterpret the data after evaluate has been applied.
for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
for (unsigned int point = 0; point < n_q_points; ++point)
const unsigned int total_number = point * dim + j;
const unsigned int new_comp = total_number / n_q_points;
const unsigned int new_point = total_number % n_q_points;
- data.contravariant[new_point][out_comp * n_lanes + in_comp]
- [new_comp] =
+ jacobians[new_point][out_comp * n_lanes +
+ in_comp][new_comp] =
eval.begin_gradients()[(out_comp * n_q_points + point) *
dim +
j][in_comp];
}
}
- if (update_flags & update_covariant_transformation)
- 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();
+ data.volume_elements[point] = jacobians[point].determinant();
// copy values from InternalData to vector given by reference
- if (update_flags & update_jacobians)
+ if (update_flags & update_covariant_transformation)
{
- 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();
+ inverse_jacobians[point] =
+ jacobians[point].covariant_form().transpose();
}
if (evaluation_flag & EvaluationFlags::hessians)
+ template <int dim, int spacedim>
+ inline DEAL_II_ALWAYS_INLINE void
+ store_vectorized_tensor(
+ const unsigned int n_points,
+ const unsigned int cur_index,
+ const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+ & derivative,
+ std::vector<DerivativeForm<1, dim, spacedim>> &result_array)
+ {
+ AssertDimension(result_array.size(), n_points);
+ constexpr unsigned int n_lanes = VectorizedArray<double>::size();
+ if (cur_index + n_lanes <= n_points)
+ {
+ std::array<unsigned int, n_lanes> indices;
+ for (unsigned int j = 0; j < n_lanes; ++j)
+ indices[j] = j * dim * spacedim;
+ const unsigned int even = (dim * spacedim) / 4 * 4;
+ double * result_ptr = &result_array[cur_index][0][0];
+ const VectorizedArray<double> *derivative_ptr = &derivative[0][0];
+ vectorized_transpose_and_store(
+ false, even, derivative_ptr, indices.data(), result_ptr);
+ for (unsigned int d = even; d < dim * spacedim; ++d)
+ for (unsigned int j = 0; j < n_lanes; ++j)
+ result_ptr[j * dim * spacedim + d] = derivative_ptr[d][j];
+ }
+ else
+ for (unsigned int j = 0; j < n_lanes && cur_index + j < n_points; ++j)
+ for (unsigned int d = 0; d < spacedim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ result_array[cur_index + j][d][e] = derivative[d][e][j];
+ }
+
+
+
template <int dim, int spacedim>
inline void
maybe_update_q_points_Jacobians_generic(
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
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,
const unsigned int n_points = unit_points.size();
const unsigned int n_lanes = VectorizedArray<double>::size();
+ // Since MappingQ::InternalData does not have separate arrays for the
+ // covariant and contravariant transformations, but uses the arrays in
+ // the `MappingRelatedData`, it can happen that vectors do not have the
+ // right size
+ if (update_flags & update_contravariant_transformation)
+ jacobians.resize(n_points);
+ if (update_flags & update_covariant_transformation)
+ inverse_jacobians.resize(n_points);
+
// 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)
polynomials_1d,
support_points,
p_vec,
- polynomial_degree == 1,
+ polynomials_1d.size() == 2,
renumber_lexicographic_to_hierarchic);
+ DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+ derivative;
+ for (unsigned int d = 0; d < spacedim; ++d)
+ for (unsigned int e = 0; e < dim; ++e)
+ derivative[d][e] = result.second[e][d];
if (update_flags & update_quadrature_points)
for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
continue;
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];
+ store_vectorized_tensor(n_points, i, derivative, jacobians);
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];
+ {
+ const auto determinant = derivative.determinant();
+ for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+ data.volume_elements[i + j] = determinant[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();
+ {
+ const auto covariant = derivative.covariant_form();
+ store_vectorized_tensor(n_points,
+ i,
+ covariant.transpose(),
+ inverse_jacobians);
+ }
}
else
{
polynomials_1d,
support_points,
unit_points[i],
- polynomial_degree == 1,
+ polynomials_1d.size() == 2,
renumber_lexicographic_to_hierarchic);
if (update_flags & update_quadrature_points)
if (cell_similarity == CellSimilarity::translation)
continue;
- if (update_flags & update_contravariant_transformation)
- {
- DerivativeForm<1, spacedim, dim> jac_transposed = result.second;
- data.contravariant[i] = jac_transposed.transpose();
- }
-
- if (update_flags & update_volume_elements)
- data.volume_elements[i] = data.contravariant[i].determinant();
+ if (dim == spacedim && update_flags & update_volume_elements)
+ data.volume_elements[i] =
+ DerivativeForm<1, spacedim, dim>(result.second).determinant();
- if (update_flags & update_jacobians)
- jacobians[i] = data.contravariant[i];
+ if (update_flags & update_contravariant_transformation)
+ jacobians[i] =
+ DerivativeForm<1, spacedim, dim>(result.second).transpose();
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();
+ inverse_jacobians[i] = jacobians[i].covariant_form().transpose();
}
}
template <int dim, int spacedim>
inline void
maybe_update_jacobian_grads(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_grads)
+ if (data.update_each & update_jacobian_grads)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points = jacobian_grads.size();
if (cell_similarity != CellSimilarity::translation)
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<2, dim> *second =
- &data.second_derivative(point + data_set, 0);
- double result[spacedim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- result[i][j][l] =
- (second[0][j][l] * 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)
- for (unsigned int l = 0; l < dim; ++l)
- result[i][j][l] +=
- (second[k][j][l] * data.mapping_support_points[k][i]);
+ const SymmetricTensor<2, dim, Tensor<1, spacedim>> second =
+ internal::evaluate_tensor_product_hessian(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic);
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim; ++j)
for (unsigned int l = 0; l < dim; ++l)
- jacobian_grads[point][i][j][l] = result[i][j][l];
+ jacobian_grads[point][i][j][l] = second[j][l][i];
}
}
}
template <int dim, int spacedim>
inline void
maybe_update_jacobian_pushed_forward_grads(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_pushed_forward_grads)
+ if (data.update_each & update_jacobian_pushed_forward_grads)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
if (cell_similarity != CellSimilarity::translation)
double tmp[spacedim][spacedim][spacedim];
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<2, dim> *second =
- &data.second_derivative(point + data_set, 0);
- double result[spacedim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- result[i][j][l] =
- (second[0][j][l] * 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)
- for (unsigned int l = 0; l < dim; ++l)
- result[i][j][l] +=
- (second[k][j][l] *
- data.mapping_support_points[k][i]);
+ const SymmetricTensor<2, dim, Tensor<1, spacedim>> second =
+ internal::evaluate_tensor_product_hessian(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic);
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[point].transpose();
// first push forward the j-components
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < spacedim; ++j)
for (unsigned int l = 0; l < dim; ++l)
{
- tmp[i][j][l] =
- result[i][0][l] * data.covariant[point][j][0];
+ tmp[i][j][l] = second[0][l][i] * covariant[j][0];
for (unsigned int jr = 1; jr < dim; ++jr)
{
tmp[i][j][l] +=
- result[i][jr][l] * data.covariant[point][j][jr];
+ second[jr][l][i] * covariant[j][jr];
}
}
for (unsigned int l = 0; l < spacedim; ++l)
{
jacobian_pushed_forward_grads[point][i][j][l] =
- tmp[i][j][0] * data.covariant[point][l][0];
+ tmp[i][j][0] * covariant[l][0];
for (unsigned int lr = 1; lr < dim; ++lr)
{
jacobian_pushed_forward_grads[point][i][j][l] +=
- tmp[i][j][lr] * data.covariant[point][l][lr];
+ tmp[i][j][lr] * covariant[l][lr];
}
}
}
+ template <int dim, int spacedim, int length_tensor>
+ inline DerivativeForm<3, dim, spacedim>
+ expand_3rd_derivatives(
+ const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
+ {
+ Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
+ DerivativeForm<3, dim, spacedim> result;
+ for (unsigned int i = 0; i < spacedim; ++i)
+ {
+ if (dim == 1)
+ result[i][0][0][0] = compressed[0][i];
+ else if (dim >= 2)
+ {
+ for (unsigned int d = 0; d < 2; ++d)
+ for (unsigned int e = 0; e < 2; ++e)
+ for (unsigned int f = 0; f < 2; ++f)
+ result[i][d][e][f] = compressed[d + e + f][i];
+ if (dim == 3)
+ {
+ AssertDimension(length_tensor, 10);
+
+ // the derivatives are ordered in rows by increasing z
+ // derivative, and in each row we have x^{(n-j)} y^{(j)} as
+ // j walks through the columns
+ for (unsigned int d = 0; d < 2; ++d)
+ for (unsigned int e = 0; e < 2; ++e)
+ {
+ result[i][d][e][2] = compressed[4 + d + e][i];
+ result[i][d][2][e] = compressed[4 + d + e][i];
+ result[i][2][d][e] = compressed[4 + d + e][i];
+ }
+ for (unsigned int d = 0; d < 2; ++d)
+ {
+ result[i][d][2][2] = compressed[7 + d][i];
+ result[i][2][d][2] = compressed[7 + d][i];
+ result[i][2][2][d] = compressed[7 + d][i];
+ }
+ result[i][2][2][2] = compressed[9][i];
+ }
+ }
+ }
+ return result;
+ }
+
+
+
/**
* Update the third derivatives of the transformation from unit to real
* cell, the Jacobian hessians.
template <int dim, int spacedim>
inline void
maybe_update_jacobian_2nd_derivatives(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_2nd_derivatives)
+ if (data.update_each & update_jacobian_2nd_derivatives)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points = jacobian_2nd_derivatives.size();
if (cell_similarity != CellSimilarity::translation)
{
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<3, dim> *third =
- &data.third_derivative(point + data_set, 0);
- double result[spacedim][dim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- result[i][j][l][m] =
- (third[0][j][l][m] *
- 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)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- result[i][j][l][m] +=
- (third[k][j][l][m] *
- data.mapping_support_points[k][i]);
-
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- jacobian_2nd_derivatives[point][i][j][l][m] =
- result[i][j][l][m];
+ jacobian_2nd_derivatives[point] = expand_3rd_derivatives<dim>(
+ internal::evaluate_tensor_product_higher_derivatives<3>(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic));
}
}
}
template <int dim, int spacedim>
inline void
maybe_update_jacobian_pushed_forward_2nd_derivatives(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
+ if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points =
jacobian_pushed_forward_2nd_derivatives.size();
if (cell_similarity != CellSimilarity::translation)
{
- double tmp[spacedim][spacedim][spacedim][spacedim];
+ dealii::ndarray<double, spacedim, spacedim, dim, dim> tmp;
+ dealii::ndarray<double, spacedim, spacedim, spacedim, dim> tmp2;
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<3, dim> *third =
- &data.third_derivative(point + data_set, 0);
- double result[spacedim][dim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- result[i][j][l][m] =
- (third[0][j][l][m] *
- 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)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- result[i][j][l][m] +=
- (third[k][j][l][m] *
- data.mapping_support_points[k][i]);
+ const DerivativeForm<3, dim, spacedim> third =
+ expand_3rd_derivatives<dim>(
+ internal::evaluate_tensor_product_higher_derivatives<3>(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic));
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[point].transpose();
// push forward the j-coordinate
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int l = 0; l < dim; ++l)
for (unsigned int m = 0; m < dim; ++m)
{
- jacobian_pushed_forward_2nd_derivatives
- [point][i][j][l][m] = result[i][0][l][m] *
- data.covariant[point][j][0];
+ tmp[i][j][l][m] =
+ third[i][0][l][m] * covariant[j][0];
for (unsigned int jr = 1; jr < dim; ++jr)
- jacobian_pushed_forward_2nd_derivatives[point][i]
- [j][l]
- [m] +=
- result[i][jr][l][m] *
- data.covariant[point][j][jr];
+ tmp[i][j][l][m] +=
+ third[i][jr][l][m] * covariant[j][jr];
}
// push forward the l-coordinate
for (unsigned int l = 0; l < spacedim; ++l)
for (unsigned int m = 0; m < dim; ++m)
{
- tmp[i][j][l][m] =
- jacobian_pushed_forward_2nd_derivatives[point][i]
- [j][0][m] *
- data.covariant[point][l][0];
+ tmp2[i][j][l][m] =
+ tmp[i][j][0][m] * covariant[l][0];
for (unsigned int lr = 1; lr < dim; ++lr)
- tmp[i][j][l][m] +=
- jacobian_pushed_forward_2nd_derivatives[point]
- [i][j]
- [lr][m] *
- data.covariant[point][l][lr];
+ tmp2[i][j][l][m] +=
+ tmp[i][j][lr][m] * covariant[l][lr];
}
// push forward the m-coordinate
{
jacobian_pushed_forward_2nd_derivatives
[point][i][j][l][m] =
- tmp[i][j][l][0] * data.covariant[point][m][0];
+ tmp2[i][j][l][0] * covariant[m][0];
for (unsigned int mr = 1; mr < dim; ++mr)
jacobian_pushed_forward_2nd_derivatives[point][i]
[j][l]
[m] +=
- tmp[i][j][l][mr] * data.covariant[point][m][mr];
+ tmp2[i][j][l][mr] * covariant[m][mr];
}
}
}
+ template <int dim, int spacedim, int length_tensor>
+ inline DerivativeForm<4, dim, spacedim>
+ expand_4th_derivatives(
+ const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
+ {
+ Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
+ DerivativeForm<4, dim, spacedim> result;
+ for (unsigned int i = 0; i < spacedim; ++i)
+ {
+ if (dim == 1)
+ result[i][0][0][0][0] = compressed[0][i];
+ else if (dim >= 2)
+ {
+ for (unsigned int d = 0; d < 2; ++d)
+ for (unsigned int e = 0; e < 2; ++e)
+ for (unsigned int f = 0; f < 2; ++f)
+ for (unsigned int g = 0; g < 2; ++g)
+ result[i][d][e][f][g] = compressed[d + e + f + g][i];
+ if (dim == 3)
+ {
+ AssertDimension(length_tensor, 15);
+
+ // the derivatives are ordered in rows by increasing z
+ // derivative, and in each row we have x^{(n-j)} y^{(j)} as
+ // j walks through the columns
+ for (unsigned int d = 0; d < 2; ++d)
+ for (unsigned int e = 0; e < 2; ++e)
+ for (unsigned int f = 0; f < 2; ++f)
+ {
+ result[i][d][e][f][2] = compressed[5 + d + e + f][i];
+ result[i][d][e][2][f] = compressed[5 + d + e + f][i];
+ result[i][d][2][e][f] = compressed[5 + d + e + f][i];
+ result[i][2][d][e][f] = compressed[5 + d + e + f][i];
+ }
+ for (unsigned int d = 0; d < 2; ++d)
+ for (unsigned int e = 0; e < 2; ++e)
+ {
+ result[i][d][e][2][2] = compressed[9 + d + e][i];
+ result[i][d][2][e][2] = compressed[9 + d + e][i];
+ result[i][d][2][2][e] = compressed[9 + d + e][i];
+ result[i][2][d][e][2] = compressed[9 + d + e][i];
+ result[i][2][d][2][e] = compressed[9 + d + e][i];
+ result[i][2][2][d][e] = compressed[9 + d + e][i];
+ }
+ for (unsigned int d = 0; d < 2; ++d)
+ {
+ result[i][d][2][2][2] = compressed[12 + d][i];
+ result[i][2][d][2][2] = compressed[12 + d][i];
+ result[i][2][2][d][2] = compressed[12 + d][i];
+ result[i][2][2][2][d] = compressed[12 + d][i];
+ }
+ result[i][2][2][2][2] = compressed[14][i];
+ }
+ }
+ }
+ return result;
+ }
+
+
+
/**
* Update the fourth derivatives of the transformation from unit to real
* cell, the Jacobian hessian gradients.
template <int dim, int spacedim>
inline void
maybe_update_jacobian_3rd_derivatives(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_3rd_derivatives)
+ if (data.update_each & update_jacobian_3rd_derivatives)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points = jacobian_3rd_derivatives.size();
if (cell_similarity != CellSimilarity::translation)
{
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<4, dim> *fourth =
- &data.fourth_derivative(point + data_set, 0);
- double result[spacedim][dim][dim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- for (unsigned int n = 0; n < dim; ++n)
- result[i][j][l][m][n] =
- (fourth[0][j][l][m][n] *
- 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)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- for (unsigned int n = 0; n < dim; ++n)
- result[i][j][l][m][n] +=
- (fourth[k][j][l][m][n] *
- data.mapping_support_points[k][i]);
-
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- for (unsigned int n = 0; n < dim; ++n)
- jacobian_3rd_derivatives[point][i][j][l][m][n] =
- result[i][j][l][m][n];
+ jacobian_3rd_derivatives[point] = expand_4th_derivatives<dim>(
+ internal::evaluate_tensor_product_higher_derivatives<4>(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic));
}
}
}
template <int dim, int spacedim>
inline void
maybe_update_jacobian_pushed_forward_3rd_derivatives(
- const CellSimilarity::Similarity cell_similarity,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ const CellSimilarity::Similarity cell_similarity,
const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+ const ArrayView<const Point<dim>> & unit_points,
+ const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+ const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
{
- const UpdateFlags update_flags = data.update_each;
- if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
+ if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
{
+ const std::vector<Point<spacedim>> &support_points =
+ data.mapping_support_points;
const unsigned int n_q_points =
jacobian_pushed_forward_3rd_derivatives.size();
if (cell_similarity != CellSimilarity::translation)
{
- double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
+ dealii::
+ ndarray<double, spacedim, spacedim, spacedim, spacedim, dim>
+ tmp;
+ dealii::ndarray<double, spacedim, spacedim, spacedim, dim, dim>
+ tmp2;
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const Tensor<4, dim> *fourth =
- &data.fourth_derivative(point + data_set, 0);
- double result[spacedim][dim][dim][dim][dim];
- for (unsigned int i = 0; i < spacedim; ++i)
- for (unsigned int j = 0; j < dim; ++j)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- for (unsigned int n = 0; n < dim; ++n)
- result[i][j][l][m][n] =
- (fourth[0][j][l][m][n] *
- 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)
- for (unsigned int l = 0; l < dim; ++l)
- for (unsigned int m = 0; m < dim; ++m)
- for (unsigned int n = 0; n < dim; ++n)
- result[i][j][l][m][n] +=
- (fourth[k][j][l][m][n] *
- data.mapping_support_points[k][i]);
+ const DerivativeForm<4, dim, spacedim> fourth =
+ expand_4th_derivatives<dim>(
+ internal::evaluate_tensor_product_higher_derivatives<4>(
+ polynomials_1d,
+ support_points,
+ unit_points[point],
+ renumber_lexicographic_to_hierarchic));
+
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[point].transpose();
// push-forward the j-coordinate
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int m = 0; m < dim; ++m)
for (unsigned int n = 0; n < dim; ++n)
{
- tmp[i][j][l][m][n] = result[i][0][l][m][n] *
- data.covariant[point][j][0];
+ tmp[i][j][l][m][n] =
+ fourth[i][0][l][m][n] * covariant[j][0];
for (unsigned int jr = 1; jr < dim; ++jr)
tmp[i][j][l][m][n] +=
- result[i][jr][l][m][n] *
- data.covariant[point][j][jr];
+ fourth[i][jr][l][m][n] * covariant[j][jr];
}
// push-forward the l-coordinate
for (unsigned int m = 0; m < dim; ++m)
for (unsigned int n = 0; n < dim; ++n)
{
- jacobian_pushed_forward_3rd_derivatives
- [point][i][j][l][m][n] =
- tmp[i][j][0][m][n] *
- data.covariant[point][l][0];
+ tmp2[i][j][l][m][n] =
+ tmp[i][j][0][m][n] * covariant[l][0];
for (unsigned int lr = 1; lr < dim; ++lr)
- jacobian_pushed_forward_3rd_derivatives[point]
- [i][j][l]
- [m][n] +=
- tmp[i][j][lr][m][n] *
- data.covariant[point][l][lr];
+ tmp2[i][j][l][m][n] +=
+ tmp[i][j][lr][m][n] * covariant[l][lr];
}
// push-forward the m-coordinate
for (unsigned int n = 0; n < dim; ++n)
{
tmp[i][j][l][m][n] =
- jacobian_pushed_forward_3rd_derivatives[point]
- [i][j][l]
- [0][n] *
- data.covariant[point][m][0];
+ tmp2[i][j][l][0][n] * covariant[m][0];
for (unsigned int mr = 1; mr < dim; ++mr)
tmp[i][j][l][m][n] +=
- jacobian_pushed_forward_3rd_derivatives
- [point][i][j][l][mr][n] *
- data.covariant[point][m][mr];
+ tmp2[i][j][l][mr][n] * covariant[m][mr];
}
// push-forward the n-coordinate
{
jacobian_pushed_forward_3rd_derivatives
[point][i][j][l][m][n] =
- tmp[i][j][l][m][0] *
- data.covariant[point][n][0];
+ tmp[i][j][l][m][0] * covariant[n][0];
for (unsigned int nr = 1; nr < dim; ++nr)
jacobian_pushed_forward_3rd_derivatives[point]
[i][j][l]
[m][n] +=
- tmp[i][j][l][m][nr] *
- data.covariant[point][n][nr];
+ tmp[i][j][l][m][nr] * covariant[n][nr];
}
}
}
const UpdateFlags update_flags = data.update_each;
if (update_flags &
- (update_boundary_forms | update_normal_vectors | update_jacobians |
- update_JxW_values | update_inverse_jacobians))
+ (update_boundary_forms | update_normal_vectors | update_JxW_values))
{
if (update_flags & update_boundary_forms)
AssertDimension(output_data.boundary_forms.size(), n_q_points);
// 0.
for (unsigned int d = 0; d != dim - 1; ++d)
{
- Assert(face_no + GeometryInfo<dim>::faces_per_cell * d <
- data.unit_tangentials.size(),
+ const unsigned int vector_index =
+ face_no + GeometryInfo<dim>::faces_per_cell * d;
+ Assert(vector_index < data.unit_tangentials.size(),
ExcInternalError());
- Assert(
- data.aux[d].size() <=
- data
- .unit_tangentials[face_no +
- GeometryInfo<dim>::faces_per_cell * d]
- .size(),
- ExcInternalError());
-
- mapping.transform(
- make_array_view(
- data.unit_tangentials[face_no +
- GeometryInfo<dim>::faces_per_cell * d]),
- mapping_contravariant,
- data,
- make_array_view(data.aux[d].begin(), data.aux[d].end()));
+ Assert(data.aux[d].size() <=
+ data.unit_tangentials[vector_index].size(),
+ ExcInternalError());
+ mapping.transform(make_array_view(
+ data.unit_tangentials[vector_index]),
+ mapping_contravariant,
+ data,
+ make_array_view(data.aux[d]));
}
if (update_flags & update_boundary_forms)
//
// to compute the cell normal, use the same method used in
// fill_fe_values for cells above
- AssertDimension(data.contravariant.size(), n_q_points);
+ AssertDimension(data.output_data->jacobians.size(),
+ n_q_points);
for (unsigned int point = 0; point < n_q_points; ++point)
{
+ const DerivativeForm<1, dim, spacedim> contravariant =
+ data.output_data->jacobians[point];
if (dim == 1)
{
// J is a tangent vector
output_data.boundary_forms[point] =
- data.contravariant[point].transpose()[0];
+ contravariant.transpose()[0];
output_data.boundary_forms[point] /=
(face_no == 0 ? -1. : +1.) *
output_data.boundary_forms[point].norm();
if (dim == 2)
{
const DerivativeForm<1, spacedim, dim> DX_t =
- data.contravariant[point].transpose();
+ contravariant.transpose();
Tensor<1, spacedim> cell_normal =
cross_product_3d(DX_t[0], DX_t[1]);
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)
{
+ const ArrayView<const Point<dim>> quadrature_points(
+ &data.quadrature_points[data_set], quadrature.size());
if (dim > 1 && data.tensor_product_quadrature)
{
maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
maybe_update_q_points_Jacobians_generic(
CellSimilarity::none,
data,
- make_array_view(&data.quadrature_points[data_set],
- &data.quadrature_points[data_set] +
- quadrature.size()),
+ quadrature_points,
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);
+ CellSimilarity::none,
+ data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
+ output_data.jacobian_grads);
}
maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
CellSimilarity::none,
- data_set,
data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_grads);
maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
CellSimilarity::none,
- data_set,
data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_2nd_derivatives);
maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
CellSimilarity::none,
- data_set,
data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_2nd_derivatives);
maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
CellSimilarity::none,
- data_set,
data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_3rd_derivatives);
maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
CellSimilarity::none,
- data_set,
data,
+ quadrature_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_3rd_derivatives);
maybe_compute_face_data(mapping,
"update_contravariant_transformation"));
for (unsigned int i = 0; i < output.size(); ++i)
- output[i] =
- apply_transformation(data.contravariant[i], input[i]);
+ output[i] = apply_transformation(data.output_data->jacobians[i],
+ input[i]);
return;
}
for (unsigned int i = 0; i < output.size(); ++i)
{
output[i] =
- apply_transformation(data.contravariant[i], input[i]);
+ apply_transformation(data.output_data->jacobians[i],
+ input[i]);
output[i] /= data.volume_elements[i];
}
return;
"update_covariant_transformation"));
for (unsigned int i = 0; i < output.size(); ++i)
- output[i] = apply_transformation(data.covariant[i], input[i]);
+ output[i] = apply_transformation(
+ data.output_data->inverse_jacobians[i].transpose(), input[i]);
return;
}
for (unsigned int i = 0; i < output.size(); ++i)
{
const DerivativeForm<1, spacedim, dim> A =
- apply_transformation(data.contravariant[i],
+ apply_transformation(data.output_data->jacobians[i],
transpose(input[i]));
- output[i] =
- apply_transformation(data.covariant[i], A.transpose());
+ output[i] = apply_transformation(
+ data.output_data->inverse_jacobians[i].transpose(),
+ A.transpose());
}
return;
for (unsigned int i = 0; i < output.size(); ++i)
{
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[i].transpose();
const DerivativeForm<1, spacedim, dim> A =
- apply_transformation(data.covariant[i],
- transpose(input[i]));
- output[i] =
- apply_transformation(data.covariant[i], A.transpose());
+ apply_transformation(covariant, transpose(input[i]));
+ output[i] = apply_transformation(covariant, A.transpose());
}
return;
for (unsigned int i = 0; i < output.size(); ++i)
{
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[i].transpose();
const DerivativeForm<1, spacedim, dim> A =
- apply_transformation(data.covariant[i], input[i]);
+ apply_transformation(covariant, input[i]);
const Tensor<2, spacedim> T =
- apply_transformation(data.contravariant[i], A.transpose());
+ apply_transformation(data.output_data->jacobians[i],
+ A.transpose());
output[i] = transpose(T);
output[i] /= data.volume_elements[i];
"update_contravariant_transformation"));
for (unsigned int q = 0; q < output.size(); ++q)
- for (unsigned int i = 0; i < spacedim; ++i)
- {
- double tmp1[dim][dim];
- for (unsigned int J = 0; J < dim; ++J)
- for (unsigned int K = 0; K < dim; ++K)
- {
- tmp1[J][K] =
- data.contravariant[q][i][0] * input[q][0][J][K];
- for (unsigned int I = 1; I < dim; ++I)
- tmp1[J][K] +=
- data.contravariant[q][i][I] * input[q][I][J][K];
- }
- for (unsigned int j = 0; j < spacedim; ++j)
- {
- double tmp2[dim];
+ {
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[q].transpose();
+ const DerivativeForm<1, dim, spacedim> contravariant =
+ data.output_data->jacobians[q];
+
+ for (unsigned int i = 0; i < spacedim; ++i)
+ {
+ double tmp1[dim][dim];
+ for (unsigned int J = 0; J < dim; ++J)
for (unsigned int K = 0; K < dim; ++K)
{
- tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
- for (unsigned int J = 1; J < dim; ++J)
- tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+ tmp1[J][K] =
+ contravariant[i][0] * input[q][0][J][K];
+ for (unsigned int I = 1; I < dim; ++I)
+ tmp1[J][K] +=
+ contravariant[i][I] * input[q][I][J][K];
}
- for (unsigned int k = 0; k < spacedim; ++k)
- {
- output[q][i][j][k] =
- data.covariant[q][k][0] * tmp2[0];
- for (unsigned int K = 1; K < dim; ++K)
- output[q][i][j][k] +=
- data.covariant[q][k][K] * tmp2[K];
- }
- }
- }
+ for (unsigned int j = 0; j < spacedim; ++j)
+ {
+ double tmp2[dim];
+ for (unsigned int K = 0; K < dim; ++K)
+ {
+ tmp2[K] = covariant[j][0] * tmp1[0][K];
+ for (unsigned int J = 1; J < dim; ++J)
+ tmp2[K] += covariant[j][J] * tmp1[J][K];
+ }
+ for (unsigned int k = 0; k < spacedim; ++k)
+ {
+ output[q][i][j][k] = covariant[k][0] * tmp2[0];
+ for (unsigned int K = 1; K < dim; ++K)
+ output[q][i][j][k] += covariant[k][K] * tmp2[K];
+ }
+ }
+ }
+ }
return;
}
"update_covariant_transformation"));
for (unsigned int q = 0; q < output.size(); ++q)
- for (unsigned int i = 0; i < spacedim; ++i)
- {
- double tmp1[dim][dim];
- for (unsigned int J = 0; J < dim; ++J)
- for (unsigned int K = 0; K < dim; ++K)
- {
- tmp1[J][K] =
- data.covariant[q][i][0] * input[q][0][J][K];
- for (unsigned int I = 1; I < dim; ++I)
- tmp1[J][K] +=
- data.covariant[q][i][I] * input[q][I][J][K];
- }
- for (unsigned int j = 0; j < spacedim; ++j)
- {
- double tmp2[dim];
+ {
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[q].transpose();
+
+ for (unsigned int i = 0; i < spacedim; ++i)
+ {
+ double tmp1[dim][dim];
+ for (unsigned int J = 0; J < dim; ++J)
for (unsigned int K = 0; K < dim; ++K)
{
- tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
- for (unsigned int J = 1; J < dim; ++J)
- tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
- }
- for (unsigned int k = 0; k < spacedim; ++k)
- {
- output[q][i][j][k] =
- data.covariant[q][k][0] * tmp2[0];
- for (unsigned int K = 1; K < dim; ++K)
- output[q][i][j][k] +=
- data.covariant[q][k][K] * tmp2[K];
+ tmp1[J][K] = covariant[i][0] * input[q][0][J][K];
+ for (unsigned int I = 1; I < dim; ++I)
+ tmp1[J][K] += covariant[i][I] * input[q][I][J][K];
}
- }
- }
+ for (unsigned int j = 0; j < spacedim; ++j)
+ {
+ double tmp2[dim];
+ for (unsigned int K = 0; K < dim; ++K)
+ {
+ tmp2[K] = covariant[j][0] * tmp1[0][K];
+ for (unsigned int J = 1; J < dim; ++J)
+ tmp2[K] += covariant[j][J] * tmp1[J][K];
+ }
+ for (unsigned int k = 0; k < spacedim; ++k)
+ {
+ output[q][i][j][k] = covariant[k][0] * tmp2[0];
+ for (unsigned int K = 1; K < dim; ++K)
+ output[q][i][j][k] += covariant[k][K] * tmp2[K];
+ }
+ }
+ }
+ }
return;
}
"update_volume_elements"));
for (unsigned int q = 0; q < output.size(); ++q)
- for (unsigned int i = 0; i < spacedim; ++i)
- {
- double factor[dim];
- for (unsigned int I = 0; I < dim; ++I)
- factor[I] =
- data.contravariant[q][i][I] / data.volume_elements[q];
- double tmp1[dim][dim];
- for (unsigned int J = 0; J < dim; ++J)
- for (unsigned int K = 0; K < dim; ++K)
- {
- tmp1[J][K] = factor[0] * input[q][0][J][K];
- for (unsigned int I = 1; I < dim; ++I)
- tmp1[J][K] += factor[I] * input[q][I][J][K];
- }
- for (unsigned int j = 0; j < spacedim; ++j)
- {
- double tmp2[dim];
+ {
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.output_data->inverse_jacobians[q].transpose();
+ const DerivativeForm<1, dim, spacedim> contravariant =
+ data.output_data->jacobians[q];
+ for (unsigned int i = 0; i < spacedim; ++i)
+ {
+ double factor[dim];
+ for (unsigned int I = 0; I < dim; ++I)
+ factor[I] =
+ contravariant[i][I] * (1. / data.volume_elements[q]);
+ double tmp1[dim][dim];
+ for (unsigned int J = 0; J < dim; ++J)
for (unsigned int K = 0; K < dim; ++K)
{
- tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
- for (unsigned int J = 1; J < dim; ++J)
- tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+ tmp1[J][K] = factor[0] * input[q][0][J][K];
+ for (unsigned int I = 1; I < dim; ++I)
+ tmp1[J][K] += factor[I] * input[q][I][J][K];
}
- for (unsigned int k = 0; k < spacedim; ++k)
- {
- output[q][i][j][k] =
- data.covariant[q][k][0] * tmp2[0];
- for (unsigned int K = 1; K < dim; ++K)
- output[q][i][j][k] +=
- data.covariant[q][k][K] * tmp2[K];
- }
- }
- }
+ for (unsigned int j = 0; j < spacedim; ++j)
+ {
+ double tmp2[dim];
+ for (unsigned int K = 0; K < dim; ++K)
+ {
+ tmp2[K] = covariant[j][0] * tmp1[0][K];
+ for (unsigned int J = 1; J < dim; ++J)
+ tmp2[K] += covariant[j][J] * tmp1[J][K];
+ }
+ for (unsigned int k = 0; k < spacedim; ++k)
+ {
+ output[q][i][j][k] = covariant[k][0] * tmp2[0];
+ for (unsigned int K = 1; K < dim; ++K)
+ output[q][i][j][k] += covariant[k][K] * tmp2[K];
+ }
+ }
+ }
+ }
return;
}
"update_covariant_transformation"));
for (unsigned int i = 0; i < output.size(); ++i)
- output[i] = apply_transformation(data.covariant[i], input[i]);
+ output[i] = apply_transformation(
+ data.output_data->inverse_jacobians[i].transpose(), input[i]);
return;
}
, n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
, line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
, tensor_product_quadrature(false)
+ , output_data(nullptr)
{}
{
return (
Mapping<dim, spacedim>::InternalDataBase::memory_consumption() +
- MemoryConsumption::memory_consumption(shape_values) +
- MemoryConsumption::memory_consumption(shape_derivatives) +
- MemoryConsumption::memory_consumption(covariant) +
- MemoryConsumption::memory_consumption(contravariant) +
+ MemoryConsumption::memory_consumption(quadrature_points) +
MemoryConsumption::memory_consumption(unit_tangentials) +
MemoryConsumption::memory_consumption(aux) +
MemoryConsumption::memory_consumption(mapping_support_points) +
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_contravariant_transformation)
- contravariant.resize(n_original_q_points);
-
if (this->update_each & update_volume_elements)
volume_elements.resize(n_original_q_points);
}
}
}
-
- 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 || needs_higher_order_terms_generic || needs_higher_order_terms)
- {
- // 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_second_derivatives.resize(n_shape_functions * n_q_points);
-
- if (this->update_each & (update_jacobian_2nd_derivatives |
- update_jacobian_pushed_forward_2nd_derivatives))
- shape_third_derivatives.resize(n_shape_functions * n_q_points);
-
- if (this->update_each & (update_jacobian_3rd_derivatives |
- update_jacobian_pushed_forward_3rd_derivatives))
- shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
-
- // now also fill the various fields with their correct values
- compute_shape_function_values(quadrature.get_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];
+ quadrature_points = quadrature.get_points();
if (dim > 1 && tensor_product_quadrature)
{
if (dim > 1)
{
if (this->update_each &
- (update_boundary_forms | update_normal_vectors | update_jacobians |
- update_JxW_values | update_inverse_jacobians))
+ (update_boundary_forms | update_normal_vectors | update_JxW_values))
{
aux.resize(dim - 1,
AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
-template <int dim, int spacedim>
-void
-MappingQ<dim, spacedim>::InternalData::compute_shape_function_values(
- const std::vector<Point<dim>> &unit_points)
-{
- const unsigned int n_points = unit_points.size();
-
- // Construct the tensor product polynomials used as shape functions for
- // the Qp mapping of cells at the boundary.
- const TensorProductPolynomials<dim> tensor_pols(
- Polynomials::generate_complete_Lagrange_basis(
- line_support_points.get_points()));
- Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
-
- // then also construct the mapping from lexicographic to the Qp shape
- // function numbering
- const std::vector<unsigned int> renumber =
- FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
-
- std::vector<double> values;
- std::vector<Tensor<1, dim>> grads;
- if (shape_values.size() != 0)
- {
- Assert(shape_values.size() == n_shape_functions * n_points,
- ExcInternalError());
- values.resize(n_shape_functions);
- }
- if (shape_derivatives.size() != 0)
- {
- Assert(shape_derivatives.size() == n_shape_functions * n_points,
- ExcInternalError());
- grads.resize(n_shape_functions);
- }
-
- std::vector<Tensor<2, dim>> grad2;
- if (shape_second_derivatives.size() != 0)
- {
- Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
- ExcInternalError());
- grad2.resize(n_shape_functions);
- }
-
- std::vector<Tensor<3, dim>> grad3;
- if (shape_third_derivatives.size() != 0)
- {
- Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
- ExcInternalError());
- grad3.resize(n_shape_functions);
- }
-
- std::vector<Tensor<4, dim>> grad4;
- if (shape_fourth_derivatives.size() != 0)
- {
- Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
- ExcInternalError());
- grad4.resize(n_shape_functions);
- }
-
-
- if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
- shape_second_derivatives.size() != 0 ||
- shape_third_derivatives.size() != 0 ||
- shape_fourth_derivatives.size() != 0)
- for (unsigned int point = 0; point < n_points; ++point)
- {
- tensor_pols.evaluate(
- unit_points[point], values, grads, grad2, grad3, grad4);
-
- if (shape_values.size() != 0)
- for (unsigned int i = 0; i < n_shape_functions; ++i)
- shape(point, i) = values[renumber[i]];
-
- if (shape_derivatives.size() != 0)
- for (unsigned int i = 0; i < n_shape_functions; ++i)
- derivative(point, i) = grads[renumber[i]];
-
- if (shape_second_derivatives.size() != 0)
- for (unsigned int i = 0; i < n_shape_functions; ++i)
- second_derivative(point, i) = grad2[renumber[i]];
-
- if (shape_third_derivatives.size() != 0)
- for (unsigned int i = 0; i < n_shape_functions; ++i)
- third_derivative(point, i) = grad3[renumber[i]];
-
- if (shape_fourth_derivatives.size() != 0)
- for (unsigned int i = 0; i < n_shape_functions; ++i)
- fourth_derivative(point, i) = grad4[renumber[i]];
- }
-}
-
-
-
template <int dim, int spacedim>
MappingQ<dim, spacedim>::MappingQ(const unsigned int p)
: polynomial_degree(p)
template <int dim, int spacedim>
MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
- : polynomial_degree(p)
- , line_support_points(
- QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
- , polynomials_1d(
- Polynomials::generate_complete_Lagrange_basis(line_support_points))
- , renumber_lexicographic_to_hierarchic(
- FETools::lexicographic_to_hierarchic_numbering<dim>(p))
- , unit_cell_support_points(
- internal::MappingQImplementation::unit_support_points<dim>(
- line_support_points,
- renumber_lexicographic_to_hierarchic))
- , support_point_weights_perimeter_to_interior(
- internal::MappingQImplementation::
- compute_support_point_weights_perimeter_to_interior(
- this->polynomial_degree,
- dim))
- , support_point_weights_cell(
- internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
- this->polynomial_degree))
-{
- Assert(p >= 1,
- ExcMessage("It only makes sense to create polynomial mappings "
- "with a polynomial degree greater or equal to one."));
-}
+ : MappingQ<dim, spacedim>(p)
+{}
// dispatch to the various specializations for spacedim=dim,
// spacedim=dim+1, etc
return internal::MappingQImplementation::
- do_transform_real_to_unit_cell_internal_codim1<1>(cell,
- p,
- initial_p_unit,
- *mdata);
+ do_transform_real_to_unit_cell_internal_codim1<1>(
+ p,
+ initial_p_unit,
+ mdata->mapping_support_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic);
}
// dispatch to the various specializations for spacedim=dim,
// spacedim=dim+1, etc
return internal::MappingQImplementation::
- do_transform_real_to_unit_cell_internal_codim1<2>(cell,
- p,
- initial_p_unit,
- *mdata);
+ do_transform_real_to_unit_cell_internal_codim1<2>(
+ p,
+ initial_p_unit,
+ mdata->mapping_support_points,
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic);
}
+
+
template <>
Point<1>
MappingQ<1, 3>::transform_real_to_unit_cell_internal(
if (out & (update_JxW_values | update_normal_vectors))
out |= update_boundary_forms;
- if (out & (update_covariant_transformation | update_JxW_values |
- update_jacobians | update_jacobian_grads |
- update_boundary_forms | update_normal_vectors))
+ if (out &
+ (update_covariant_transformation | update_jacobian_grads |
+ update_jacobians | update_boundary_forms | update_normal_vectors))
out |= update_contravariant_transformation;
if (out &
Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(internal_data);
+ data.output_data = &output_data;
const unsigned int n_q_points = quadrature.size();
data,
make_array_view(quadrature.get_points()),
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data.quadrature_points,
output_data.jacobians,
internal::MappingQImplementation::maybe_update_jacobian_grads<dim,
spacedim>(
computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_grads);
}
internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads<
dim,
spacedim>(computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_grads);
internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives<
dim,
spacedim>(computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_2nd_derivatives);
internal::MappingQImplementation::
maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_2nd_derivatives);
internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives<
dim,
spacedim>(computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_3rd_derivatives);
internal::MappingQImplementation::
maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
computed_cell_similarity,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_3rd_derivatives);
const UpdateFlags update_flags = data.update_each;
{
if (dim == spacedim)
{
- const double det = data.contravariant[point].determinant();
+ const double det = data.volume_elements[point];
// check for distorted cells.
Tensor<1, spacedim> DX_t[dim];
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < dim; ++j)
- DX_t[j][i] = data.contravariant[point][i][j];
+ DX_t[j][i] = output_data.jacobians[point][i][j];
Tensor<2, dim> G; // First fundamental form
for (unsigned int i = 0; i < dim; ++i)
Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(internal_data);
+ data.output_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)
+ // 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()) ||
quadrature[0],
data,
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data);
}
Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(internal_data);
+ data.output_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)
+ // 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()) ||
quadrature,
data,
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data);
}
Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(internal_data);
+ data.output_data = &output_data;
const unsigned int n_q_points = quadrature.size();
data,
make_array_view(quadrature.get_points()),
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data.quadrature_points,
output_data.jacobians,
internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_grads);
internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads<
dim,
spacedim>(CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_grads);
internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives<
dim,
spacedim>(CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_2nd_derivatives);
internal::MappingQImplementation::
maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_2nd_derivatives);
internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives<
dim,
spacedim>(CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_3rd_derivatives);
internal::MappingQImplementation::
maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
CellSimilarity::none,
- QProjector<dim>::DataSetDescriptor::cell(),
data,
+ make_array_view(quadrature.get_points()),
+ polynomials_1d,
+ renumber_lexicographic_to_hierarchic,
output_data.jacobian_pushed_forward_3rd_derivatives);
const UpdateFlags update_flags = data.update_each;
for (unsigned int point = 0; point < n_q_points; ++point)
{
- const double det = data.contravariant[point].determinant();
+ const double det = data.volume_elements[point];
// check for distorted cells.
// The normals are n = J^{-T} * \hat{n} before normalizing.
Tensor<1, spacedim> normal;
for (unsigned int d = 0; d < spacedim; d++)
- normal[d] =
- data.covariant[point][d] * quadrature.normal_vector(point);
+ normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
+ quadrature.normal_vector(point);
output_data.JxW_values[point] = weights[point] * det * normal.norm();
Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
unit_points.end())));
const InternalData &data = static_cast<const InternalData &>(*internal_data);
+ data.output_data = &output_data;
data.mapping_support_points = this->compute_mapping_support_points(cell);
internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
data,
unit_points,
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data.quadrature_points,
output_data.jacobians,
const InternalData &data = static_cast<const InternalData &>(internal_data);
data.mapping_support_points = this->compute_mapping_support_points(cell);
+ data.output_data = &output_data;
internal::MappingQImplementation::do_fill_fe_face_values(
*this,
face_quadrature,
data,
polynomials_1d,
- polynomial_degree,
renumber_lexicographic_to_hierarchic,
output_data);
}
AssertDimension(input.size(), output.size());
Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
ExcInternalError());
- const InternalData &data = static_cast<const InternalData &>(mapping_data);
+ const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &data = *static_cast<const InternalData &>(mapping_data).output_data;
switch (mapping_kind)
{
case mapping_covariant_gradient:
{
- Assert(data.update_each & update_contravariant_transformation,
+ Assert(!data.inverse_jacobians.empty(),
typename FEValuesBase<dim>::ExcAccessToUninitializedField(
"update_covariant_transformation"));
for (unsigned int i = 0; i < spacedim; ++i)
for (unsigned int j = 0; j < spacedim; ++j)
{
- double tmp[dim];
+ double tmp[dim];
+ const DerivativeForm<1, dim, spacedim> covariant =
+ data.inverse_jacobians[q].transpose();
for (unsigned int K = 0; K < dim; ++K)
{
- tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
+ tmp[K] = covariant[j][0] * input[q][i][0][K];
for (unsigned int J = 1; J < dim; ++J)
- tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
+ tmp[K] += covariant[j][J] * input[q][i][J][K];
}
for (unsigned int k = 0; k < spacedim; ++k)
{
- output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
+ output[q][i][j][k] = covariant[k][0] * tmp[0];
for (unsigned int K = 1; K < dim; ++K)
- output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
+ output[q][i][j][k] += covariant[k][K] * tmp[K];
}
}
return;
projects / dealii.git / commitdiff