this->mapped_geometry->get_data_storage().JxW_values.begin();
this->jacobian_gradients =
this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+ this->jacobian_gradients_non_inverse =
+ this->mapped_geometry->get_data_storage()
+ .jacobian_gradients_non_inverse[0]
+ .begin();
this->quadrature_points =
this->mapped_geometry->get_data_storage().quadrature_points.begin();
}
this->mapped_geometry->get_data_storage().JxW_values.begin();
this->jacobian_gradients =
this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+ this->jacobian_gradients_non_inverse =
+ this->mapped_geometry->get_data_storage()
+ .jacobian_gradients_non_inverse[0]
+ .begin();
this->quadrature_points =
this->mapped_geometry->get_data_storage().quadrature_points.begin();
}
else
{
// General cell
- // Here we need the jacobian gradient and not the inverse which is
- // stored in this->jacobian_gradients
- AssertThrow(false, ExcNotImplemented());
+
+ // This assert could be removed if we make sure that this is updated
+ // even though update_hessians or update_jacobian_grads is not passed,
+ // i.e make the necessary changes in
+ // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
+ Assert(this->jacobian_gradients_non_inverse != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_hessians"));
+
+ const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ this->jacobian[q_point];
+ const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
+
+ // Derivatives are reordered for faces. Need to take this into account
+ const VectorizedArrayType inv_det =
+ (is_face && dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac);
+
+ VectorizedArrayType tmp;
+ // J * grad_quad * J^-1 * det(J^-1)
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int e = 0; e < dim; ++e)
+ tmp += t_jac[f][comp] * inv_t_jac[d][e] *
+ this->gradients_quad[(f * dim + e) * nqp + q_point];
+
+ grad_out[comp][d] = tmp * inv_det;
+ }
+
+ // Contribution from values
+ {
+ // Diagonal part of jac_grad
+
+ // Add jac_grad * J^{-1} * values * det(J^{-1})
+ // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ {
+ tmp = jac_grad[0][i] * inv_t_jac[j][0] *
+ this->values_quad[q_point];
+ for (unsigned int f = 1; f < dim; ++f)
+ tmp += jac_grad[f][i] * inv_t_jac[j][f] *
+ this->values_quad[f * nqp + q_point];
+
+ grad_out[i][j] += tmp * inv_det;
+ }
+
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ {
+ tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int n = 0; n < dim; ++n)
+ for (unsigned int m = 0; m < dim; ++m)
+ tmp += inv_t_jac[m][f] * jac_grad[f][m] *
+ inv_t_jac[j][f] * t_jac[n][i] *
+ this->values_quad[n * nqp + q_point];
+ grad_out[i][j] -= tmp * inv_det;
+ }
+ }
+
+ {
+ // Off-diagonal part of jac_grad
+
+ // Add jac_grad * J^{-1} * values * det(J^{-1})
+ // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ {
+ tmp = 0;
+ for (unsigned int r = 0, f = dim; r < dim; ++r)
+ for (unsigned int k = r + 1; k < dim; ++k, ++f)
+ {
+ tmp += jac_grad[f][i] *
+ (inv_t_jac[j][k] *
+ this->values_quad[r * nqp + q_point] +
+ inv_t_jac[j][r] *
+ this->values_quad[k * nqp + q_point]);
+ for (unsigned int n = 0; n < dim; ++n)
+ for (unsigned int m = 0; m < dim; ++m)
+ tmp -= jac_grad[f][m] * t_jac[n][i] *
+ this->values_quad[n * nqp + q_point] *
+ (inv_t_jac[m][k] * inv_t_jac[j][r] +
+ inv_t_jac[m][r] * inv_t_jac[j][k]);
+ }
+ grad_out[i][j] += tmp * inv_det;
+ }
+ }
}
return grad_out;
}
divergence +=
this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
}
- else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ else
{
- // Affine cell
+ // General cell
// Derivatives are reordered for faces. Need to take this into account
const VectorizedArrayType inv_det =
- (is_face && dim == 2 && this->get_face_no() < 2) ?
- -determinant(this->jacobian[0]) :
- determinant(this->jacobian[0]);
+ determinant(
+ this->jacobian[this->cell_type >
+ internal::MatrixFreeFunctions::affine ?
+ q_point :
+ 0]) *
+ Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
// div * det(J^-1)
divergence = this->gradients_quad[q_point] * inv_det;
divergence +=
this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
}
- else
- {
- // General cell
- Assert(false, ExcNotImplemented());
- }
}
else
{
VectorizedArrayType tmp = 0;
for (unsigned int f = 0; f < dim; ++f)
for (unsigned int e = 0; e < dim; ++e)
- tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
+ tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e];
- this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+ tmp * fac;
}
}
else
{
// General cell
- AssertThrow(false, ExcNotImplemented());
+
+ const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
+ const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
+ this->jacobian[q_point];
+ const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
+
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ (!is_face) ?
+ this->quadrature_weights[q_point] :
+ this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
+ -determinant(inv_t_jac) :
+ determinant(inv_t_jac));
+
+ VectorizedArrayType tmp;
+ // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ tmp = 0;
+ for (unsigned int f = 0; f < dim; ++f)
+ for (unsigned int e = 0; e < dim; ++e)
+ tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e];
+
+ this->gradients_quad[(comp * dim + d) * nqp + q_point] =
+ tmp * fac;
+ }
+
+ // Contribution from values
+ {
+ // Diagonal part of jac_grad
+
+ // Add jac_grad * J^{-1} * values * factor
+ // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
+ for (unsigned int f = 0; f < dim; ++f)
+ {
+ tmp = 0;
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ {
+ tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
+ for (unsigned int m = 0; m < dim; ++m)
+ for (unsigned int k = 0; k < dim; ++k)
+ tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
+ inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
+ }
+ this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
+ }
+ }
+
+ {
+ // Off-diagonal part of jac_grad
+
+ // Add jac_grad * J^{-1} * values * factor
+ for (unsigned int r = 0, f = dim; r < dim; ++r)
+ for (unsigned int k = r + 1; k < dim; ++k, ++f)
+ {
+ tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
+ for (unsigned int j = 1; j < dim; ++j)
+ tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
+ for (unsigned int i = 1; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
+ this->values_from_gradients_quad[r * nqp + q_point] +=
+ tmp * fac;
+
+ tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
+ for (unsigned int j = 1; j < dim; ++j)
+ tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
+ for (unsigned int i = 1; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
+ this->values_from_gradients_quad[k * nqp + q_point] +=
+ tmp * fac;
+ }
+
+ // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
+ for (unsigned int n = 0; n < dim; ++n)
+ {
+ tmp = 0;
+ for (unsigned int r = 0, f = dim; r < dim; ++r)
+ for (unsigned int k = r + 1; k < dim; ++k, ++f)
+ for (unsigned int i = 0; i < dim; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ for (unsigned int m = 0; m < dim; ++m)
+ tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
+ (inv_t_jac[m][k] * inv_t_jac[j][r] +
+ inv_t_jac[m][r] * inv_t_jac[j][k]);
+
+ this->values_from_gradients_quad[n * nqp + q_point] -=
+ tmp * fac;
+ }
+ }
}
}
else
if (this->data->element_type ==
internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
{
- if (this->cell_type <= internal::MatrixFreeFunctions::affine)
- {
- // Affine cell
+ // General cell
+
+ // Derivatives are reordered for faces. Need to take this into account
+ // and 1/inv_det != J_value for faces
+ const VectorizedArrayType fac =
+ (!is_face) ?
+ this->quadrature_weights[q_point] * div_in :
+ (this->cell_type > internal::MatrixFreeFunctions::affine ?
+ this->J_value[q_point] :
+ this->J_value[0] * this->quadrature_weights[q_point]) *
+ div_in *
+ determinant(
+ this->jacobian[this->cell_type >
+ internal::MatrixFreeFunctions::affine ?
+ q_point :
+ 0]) *
+ Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
- // Derivatives are reordered for faces. Need to take this into account
- // and 1/inv_det != J_value for faces
- const VectorizedArrayType fac =
- ((!is_face) ?
- 1 :
- this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
- -determinant(this->jacobian[0]) :
- determinant(this->jacobian[0]))) *
- this->quadrature_weights[q_point] * div_in;
-
- for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
+ for (unsigned int e = d + 1; e < dim; ++e)
{
- this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
- for (unsigned int e = d + 1; e < dim; ++e)
- {
- this->gradients_quad[(dim * d + e) * nqp + q_point] =
- VectorizedArrayType();
- this->gradients_quad[(dim * e + d) * nqp + q_point] =
- VectorizedArrayType();
- }
+ this->gradients_quad[(dim * d + e) * nqp + q_point] =
+ VectorizedArrayType();
+ this->gradients_quad[(dim * e + d) * nqp + q_point] =
+ VectorizedArrayType();
}
}
- else
- {
- // General cell
- AssertThrow(false, ExcNotImplemented());
- }
+ this->divergence_is_requested = true;
}
else
{
this->J_value = &this->mapping_data->JxW_values[offsets];
this->jacobian_gradients =
this->mapping_data->jacobian_gradients[0].data() + offsets;
+ this->jacobian_gradients_non_inverse =
+ this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
unsigned int i = 0;
for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
auto &this_jacobian_data = mapping_storage.jacobians[0];
auto &this_J_value_data = mapping_storage.JxW_values;
auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
- auto &this_quadrature_points_data = mapping_storage.quadrature_points;
+ auto &this_jacobian_gradients_non_inverse_data =
+ mapping_storage.jacobian_gradients_non_inverse[0];
+ auto &this_quadrature_points_data = mapping_storage.quadrature_points;
if (this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine)
{
if (this->mapping_data->jacobian_gradients[0].size() > 0)
this_jacobian_gradients_data.resize_fast(1);
+ if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
+ this_jacobian_gradients_non_inverse_data.resize_fast(1);
+
if (this->mapping_data->quadrature_points.size() > 0)
this_quadrature_points_data.resize_fast(1);
}
if (this->mapping_data->jacobian_gradients[0].size() > 0)
this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
+ if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
+ this_jacobian_gradients_non_inverse_data.resize_fast(
+ this->n_quadrature_points);
+
if (this->mapping_data->quadrature_points.size() > 0)
this_quadrature_points_data.resize_fast(this->n_quadrature_points);
}
this->jacobian = this_jacobian_data.data();
this->J_value = this_J_value_data.data();
this->jacobian_gradients = this_jacobian_gradients_data.data();
- this->quadrature_points = this_quadrature_points_data.data();
+ this->jacobian_gradients_non_inverse =
+ this_jacobian_gradients_non_inverse_data.data();
+ this->quadrature_points = this_quadrature_points_data.data();
// fill internal data storage lane by lane
for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
this->mapping_data
->jacobian_gradients[0][offsets + q][i][j][lane];
+ if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
+ for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ this_jacobian_gradients_non_inverse_data[q][i][j][v] =
+ this->mapping_data
+ ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
+ [lane];
+
if (this->mapping_data->quadrature_points.size() > 0)
for (unsigned int i = 0; i < dim; ++i)
this_quadrature_points_data[q][i][v] =
this->mapping_data
->jacobian_gradients[0][offsets + q_src][i][j][lane];
+ if (this->mapping_data->jacobian_gradients_non_inverse[0].size() >
+ 0)
+ for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
+ for (unsigned int j = 0; j < dim; ++j)
+ this_jacobian_gradients_non_inverse_data[q][i][j][v] =
+ this->mapping_data
+ ->jacobian_gradients_non_inverse[0][offsets + q_src][i]
+ [j][lane];
+
if (this->mapping_data->quadrature_points.size() > 0)
{
if (cell_type <=
if (hessians_on_general_cells)
evaluation_flag_actual |= EvaluationFlags::gradients;
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
+ evaluation_flag & EvaluationFlags::gradients &&
+ (this->cell_type > internal::MatrixFreeFunctions::affine))
+ evaluation_flag_actual |= EvaluationFlags::values;
+
if (fe_degree > -1)
{
SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
}
}
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
+ integration_flag & EvaluationFlags::gradients &&
+ this->cell_type > internal::MatrixFreeFunctions::affine &&
+ this->divergence_is_requested == false)
+ {
+ unsigned int size = n_components * n_q_points;
+ if ((integration_flag & EvaluationFlags::values) != 0u)
+ {
+ for (unsigned int i = 0; i < size; ++i)
+ this->values_quad[i] += this->values_from_gradients_quad[i];
+ }
+ else
+ {
+ for (unsigned int i = 0; i < size; ++i)
+ this->values_quad[i] = this->values_from_gradients_quad[i];
+ integration_flag_actual |= EvaluationFlags::values;
+ }
+ }
+
if (fe_degree > -1)
{
SelectEvaluator<dim, fe_degree, n_q_points_1d, VectorizedArrayType>::
this->jacobian_gradients =
this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
offsets;
+ this->jacobian_gradients_non_inverse =
+ this->mapping_data
+ ->jacobian_gradients_non_inverse[!this->is_interior_face()]
+ .data() +
+ offsets;
if (this->mapping_data->quadrature_point_offsets.empty() == false)
{
this->jacobian_gradients =
this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
offsets;
+ this->jacobian_gradients_non_inverse =
+ this->mapping_data
+ ->jacobian_gradients_non_inverse[!this->is_interior_face()]
+ .data() +
+ offsets;
if (this->matrix_free->get_mapping_info()
.face_data_by_cells[this->quad_no]
if (hessians_on_general_cells)
evaluation_flag_actual |= EvaluationFlags::gradients;
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
+ evaluation_flag & EvaluationFlags::gradients &&
+ (this->cell_type > internal::MatrixFreeFunctions::affine))
+ evaluation_flag_actual |= EvaluationFlags::values;
+
if (fe_degree > -1)
internal::FEFaceEvaluationImplEvaluateSelector<dim, VectorizedArrayType>::
template run<fe_degree, n_q_points_1d>(n_components,
}
}
+ if (this->data->element_type ==
+ internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas &&
+ integration_flag & EvaluationFlags::gradients &&
+ this->cell_type > internal::MatrixFreeFunctions::affine &&
+ this->divergence_is_requested == false)
+ {
+ unsigned int size = n_components * n_q_points;
+ if ((integration_flag & EvaluationFlags::values) != 0u)
+ {
+ for (unsigned int i = 0; i < size; ++i)
+ this->values_quad[i] += this->values_from_gradients_quad[i];
+ }
+ else
+ {
+ for (unsigned int i = 0; i < size; ++i)
+ this->values_quad[i] = this->values_from_gradients_quad[i];
+ integration_flag_actual |= EvaluationFlags::values;
+ }
+ }
+
if (fe_degree > -1)
internal::FEFaceEvaluationImplIntegrateSelector<dim, VectorizedArrayType>::
template run<fe_degree, n_q_points_1d>(n_components,
const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
*jacobian_gradients;
+ /**
+ * A pointer to the gradients of the Jacobian transformation of the
+ * present cell.
+ */
+ const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
+ *jacobian_gradients_non_inverse;
+
/**
* A pointer to the Jacobian determinant of the present cell. If on a
* Cartesian cell or on a cell with constant Jacobian, this is just the
*/
Number *values_quad;
+ /**
+ * This field stores the values of the finite element function on
+ * quadrature points after applying unit cell transformations or before
+ * integrating. This field is accessed when performing the contravariant
+ * Piola transform for gradients on general cells. This is done by the
+ * functions get_gradient() and submit_gradient() when using a H(div)-
+ * conforming finite element such as FE_RaviartThomasNodal.
+ *
+ * The values of this array are stored in the start section of
+ * @p scratch_data_array. Due to its access as a thread local memory, the
+ * memory can get reused between different calls.
+ */
+ Number *values_from_gradients_quad;
+
/**
* This field stores the gradients of the finite element function on
* quadrature points after applying unit cell transformations or before
internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>>
mapped_geometry;
+ /**
+ * Bool indicating if the divergence is requested. Used internally in the case
+ * of the Piola transform.
+ */
+ bool divergence_is_requested;
+
// Make FEEvaluation and FEEvaluationBase objects friends for access to
// protected member mapped_geometry.
template <int, int, typename, bool, typename>
, quadrature_points(nullptr)
, jacobian(nullptr)
, jacobian_gradients(nullptr)
+ , jacobian_gradients_non_inverse(nullptr)
, J_value(nullptr)
, normal_vectors(nullptr)
, normal_x_jacobian(nullptr)
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
, subface_index(0)
, cell_type(internal::MatrixFreeFunctions::general)
+ , divergence_is_requested(false)
{}
, quadrature_points(nullptr)
, jacobian(nullptr)
, jacobian_gradients(nullptr)
+ , jacobian_gradients_non_inverse(nullptr)
, J_value(nullptr)
, normal_vectors(nullptr)
, normal_x_jacobian(nullptr)
, dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
, mapped_geometry(mapped_geometry)
, is_reinitialized(false)
+ , divergence_is_requested(false)
{
mapping_data = &mapped_geometry->get_data_storage();
jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
J_value = mapped_geometry->get_data_storage().JxW_values.begin();
jacobian_gradients =
mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
+ jacobian_gradients_non_inverse = mapped_geometry->get_data_storage()
+ .jacobian_gradients_non_inverse[0]
+ .begin();
quadrature_points =
mapped_geometry->get_data_storage().quadrature_points.begin();
}
AssertDimension(active_quad_index, other.active_quad_index);
AssertDimension(n_quadrature_points, descriptor->n_q_points);
- data = other.data;
- dof_info = other.dof_info;
- mapping_data = other.mapping_data;
- descriptor = other.descriptor;
- jacobian = nullptr;
- J_value = nullptr;
- normal_vectors = nullptr;
- normal_x_jacobian = nullptr;
- jacobian_gradients = nullptr;
- quadrature_points = nullptr;
- quadrature_weights = other.quadrature_weights;
+ data = other.data;
+ dof_info = other.dof_info;
+ mapping_data = other.mapping_data;
+ descriptor = other.descriptor;
+ jacobian = nullptr;
+ J_value = nullptr;
+ normal_vectors = nullptr;
+ normal_x_jacobian = nullptr;
+ jacobian_gradients = nullptr;
+ jacobian_gradients_non_inverse = nullptr;
+ quadrature_points = nullptr;
+ quadrature_weights = other.quadrature_weights;
# ifdef DEBUG
is_reinitialized = false;
internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
- face_numbers[0] = 0;
- face_orientations[0] = 0;
- subface_index = 0;
- cell_type = internal::MatrixFreeFunctions::general;
+ face_numbers[0] = 0;
+ face_orientations[0] = 0;
+ subface_index = 0;
+ cell_type = internal::MatrixFreeFunctions::general;
+ divergence_is_requested = false;
return *this;
}
2 * n_quadrature_points;
const unsigned int size_data_arrays =
n_components * dofs_per_component +
- (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+ (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 2) *
n_quadrature_points);
const unsigned int allocated_size = size_scratch_data + size_data_arrays;
// set the pointers to the correct position in the data array
values_dofs = scratch_data_array->begin();
values_quad = scratch_data_array->begin() + n_components * dofs_per_component;
- gradients_quad = scratch_data_array->begin() +
- n_components * (dofs_per_component + n_quadrature_points);
+ values_from_gradients_quad =
+ scratch_data_array->begin() +
+ n_components * (dofs_per_component + n_quadrature_points);
+ gradients_quad =
+ scratch_data_array->begin() +
+ n_components * (dofs_per_component + 2 * n_quadrature_points);
gradients_from_hessians_quad =
scratch_data_array->begin() +
- n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
+ n_components * (dofs_per_component + (dim + 2) * n_quadrature_points);
hessians_quad =
scratch_data_array->begin() +
- n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+ n_components * (dofs_per_component + (2 * dim + 2) * n_quadrature_points);
}
data.first[my_q].jacobians[0].push_back(inv_jac);
if (update_flags & update_jacobian_grads)
- data.first[my_q].jacobian_gradients[0].push_back(
- process_jacobian_gradient(inv_jac,
- inv_jac,
- jacobian_grad));
+ {
+ data.first[my_q].jacobian_gradients[0].push_back(
+ process_jacobian_gradient(inv_jac,
+ inv_jac,
+ jacobian_grad));
+ Tensor<1,
+ dim *(dim + 1) / 2,
+ Tensor<1, dim, VectorizedArrayType>>
+ jac_grad_sym;
+ // the diagonal part of Jacobian gradient comes
+ // first
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int f = 0; f < dim; ++f)
+ jac_grad_sym[d][f] = jacobian_grad[f][d][d];
+
+ // then the upper-diagonal part
+ for (unsigned int d = 0, count = dim; d < dim; ++d)
+ for (unsigned int e = d + 1; e < dim; ++e, ++count)
+ for (unsigned int f = 0; f < dim; ++f)
+ jac_grad_sym[count][f] = jacobian_grad[f][d][e];
+
+ data.first[my_q]
+ .jacobian_gradients_non_inverse[0]
+ .push_back(jac_grad_sym);
+ }
}
}
data_cells_local.jacobian_gradients[i].end(),
data_cells.jacobian_gradients[i].begin() + data_shift[0]);
data_cells_local.jacobian_gradients[i].clear();
+ std::copy(
+ data_cells_local.jacobian_gradients_non_inverse[i].begin(),
+ data_cells_local.jacobian_gradients_non_inverse[i].end(),
+ data_cells.jacobian_gradients_non_inverse[i].begin() +
+ data_shift[0]);
+ data_cells_local.jacobian_gradients_non_inverse[i].clear();
std::copy(data_cells_local.normals_times_jacobians[i].begin(),
data_cells_local.normals_times_jacobians[i].end(),
data_cells.normals_times_jacobians[i].begin() +
inv_jac_grad[d][e],
vv,
my_data.jacobian_gradients[0][idx][d][e]);
+
+ // Also store the non-inverse jacobian gradient.
+ // the diagonal part of Jacobian gradient comes
+ // first
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int f = 0; f < dim; ++f)
+ store_vectorized_array(
+ jac_grad[f][d][d],
+ vv,
+ my_data.jacobian_gradients_non_inverse[0][idx]
+ [d][f]);
+
+ // then the upper-diagonal part
+ for (unsigned int d = 0, count = dim; d < dim; ++d)
+ for (unsigned int e = d + 1; e < dim;
+ ++e, ++count)
+ for (unsigned int f = 0; f < dim; ++f)
+ store_vectorized_array(
+ jac_grad[f][d][e],
+ vv,
+ my_data.jacobian_gradients_non_inverse
+ [0][idx][count][f]);
}
}
}
cell_data[my_q].jacobians[0].resize_fast(
cell_data[my_q].JxW_values.size());
if (update_flags_cells & update_jacobian_grads)
- cell_data[my_q].jacobian_gradients[0].resize_fast(
- cell_data[my_q].JxW_values.size());
+ {
+ cell_data[my_q].jacobian_gradients[0].resize_fast(
+ cell_data[my_q].JxW_values.size());
+ cell_data[my_q].jacobian_gradients_non_inverse[0].resize_fast(
+ cell_data[my_q].JxW_values.size());
+ }
if (update_flags_cells & update_quadrature_points)
{
cell_data[my_q].quadrature_point_offsets.resize(cell_type.size());
my_data.jacobian_gradients[is_exterior]
[offset + q][d][e]);
}
+
+ // Also store the non-inverse jacobian gradient.
+ // the diagonal part of Jacobian gradient comes first.
+ // jac_grad already has its derivatives reordered,
+ // so no need to compensate for this here
+ for (unsigned int d = 0; d < dim; ++d)
+ for (unsigned int f = 0; f < dim; ++f)
+ store_vectorized_array(
+ jac_grad[f][d][d],
+ vv,
+ my_data.jacobian_gradients_non_inverse[is_exterior]
+ [offset + q]
+ [d][f]);
+
+ // then the upper-diagonal part
+ for (unsigned int d = 0, count = dim; d < dim; ++d)
+ for (unsigned int e = d + 1; e < dim; ++e, ++count)
+ for (unsigned int f = 0; f < dim; ++f)
+ store_vectorized_array(
+ jac_grad[f][d][e],
+ vv,
+ my_data.jacobian_gradients_non_inverse
+ [is_exterior][offset + q][count][f]);
}
};
face_data[my_q].JxW_values.size());
face_data[my_q].jacobian_gradients[1].resize_fast(
face_data[my_q].JxW_values.size());
+ face_data[my_q].jacobian_gradients_non_inverse[0].resize_fast(
+ face_data[my_q].JxW_values.size());
+ face_data[my_q].jacobian_gradients_non_inverse[1].resize_fast(
+ face_data[my_q].JxW_values.size());
}
face_data[my_q].normals_times_jacobians[0].resize_fast(
face_data[my_q].JxW_values.size());
my_data.JxW_values.resize_fast(max_size);
my_data.jacobians[0].resize_fast(max_size);
if (update_flags_cells & update_jacobian_grads)
- my_data.jacobian_gradients[0].resize_fast(max_size);
+ {
+ my_data.jacobian_gradients[0].resize_fast(max_size);
+ my_data.jacobian_gradients_non_inverse[0].resize_fast(max_size);
+ }
if (update_flags_cells & update_quadrature_points)
{
{
my_data.jacobian_gradients[0].resize_fast(max_size);
my_data.jacobian_gradients[1].resize_fast(max_size);
+ my_data.jacobian_gradients_non_inverse[0].resize_fast(max_size);
+ my_data.jacobian_gradients_non_inverse[1].resize_fast(max_size);
}
my_data.normals_times_jacobians[0].resize_fast(max_size);
my_data.normals_times_jacobians[1].resize_fast(max_size);
face_data_by_cells[my_q].normals_times_jacobians[1].resize_fast(
storage_length * GeometryInfo<dim>::faces_per_cell);
if (update_flags & update_jacobian_grads)
- face_data_by_cells[my_q].jacobian_gradients[0].resize_fast(
- storage_length * GeometryInfo<dim>::faces_per_cell);
+ {
+ face_data_by_cells[my_q].jacobian_gradients[0].resize_fast(
+ storage_length * GeometryInfo<dim>::faces_per_cell);
+ face_data_by_cells[my_q]
+ .jacobian_gradients_non_inverse[0]
+ .resize_fast(storage_length *
+ GeometryInfo<dim>::faces_per_cell);
+ }
if (update_flags & update_quadrature_points)
face_data_by_cells[my_q].quadrature_points.resize_fast(
2>
jacobian_gradients;
+ /**
+ * The storage of the gradients of the Jacobian transformation. Because of
+ * symmetry, only the upper diagonal and diagonal part are needed. The
+ * first index runs through the derivatives, starting with the diagonal
+ * and then continuing row-wise, i.e., $\partial^2/\partial x_1 \partial
+ * x_2$ first, then
+ * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index
+ * is the spatial coordinate.
+ *
+ * Indexed by @p data_index_offsets.
+ *
+ * Contains two fields for access from both sides for interior faces,
+ * but the default case (cell integrals or boundary integrals) only
+ * fills the zeroth component and ignores the first one.
+ */
+ std::array<
+ AlignedVector<
+ Tensor<1, spacedim *(spacedim + 1) / 2, Tensor<1, spacedim, Number>>>,
+ 2>
+ jacobian_gradients_non_inverse;
+
/**
* Stores the Jacobian transformations times the normal vector (this
* represents a shortcut that is accessed often and can thus get higher
{
jacobians[i].clear();
jacobian_gradients[i].clear();
+ jacobian_gradients_non_inverse[i].clear();
normals_times_jacobians[i].clear();
}
quadrature_point_offsets.clear();
MemoryConsumption::memory_consumption(jacobians[1]) +
MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
MemoryConsumption::memory_consumption(jacobian_gradients[1]) +
+ MemoryConsumption::memory_consumption(
+ jacobian_gradients_non_inverse[0]) +
+ MemoryConsumption::memory_consumption(
+ jacobian_gradients_non_inverse[1]) +
MemoryConsumption::memory_consumption(normals_times_jacobians[0]) +
MemoryConsumption::memory_consumption(normals_times_jacobians[1]) +
MemoryConsumption::memory_consumption(quadrature_point_offsets) +
task_info.print_memory_statistics(
out,
MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
- MemoryConsumption::memory_consumption(jacobian_gradients[1]));
+ MemoryConsumption::memory_consumption(jacobian_gradients[1]) +
+ MemoryConsumption::memory_consumption(
+ jacobian_gradients_non_inverse[0]) +
+ MemoryConsumption::memory_consumption(
+ jacobian_gradients_non_inverse[1]));
}
const std::size_t normal_size =
Utilities::MPI::sum(normal_vectors.size(), task_info.communicator);
DEAL:2d::Number of degrees of freedom: 144
DEAL:2d::
DEAL:2d::Testing Values
-DEAL:2d::Norm of difference: 6.05894e-16
+DEAL:2d::Norm of difference: 2.61453e-15
DEAL:2d::
DEAL:2d::Testing Gradients
-DEAL:2d::Norm of difference: 6.33970e-16
+DEAL:2d::Norm of difference: 6.23772e-15
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 5.24419e-16
+DEAL:2d::Norm of difference: 3.80204e-15
DEAL:2d::
DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 16
DEAL:2d::Number of degrees of freedom: 312
DEAL:2d::
DEAL:2d::Testing Values
-DEAL:2d::Norm of difference: 6.67799e-16
+DEAL:2d::Norm of difference: 7.43226e-15
DEAL:2d::
DEAL:2d::Testing Gradients
-DEAL:2d::Norm of difference: 1.06191e-15
+DEAL:2d::Norm of difference: 1.08102e-14
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 5.38383e-16
+DEAL:2d::Norm of difference: 2.64297e-15
DEAL:2d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 64
DEAL:3d::Number of degrees of freedom: 1728
DEAL:3d::
DEAL:3d::Testing Values
-DEAL:3d::Norm of difference: 6.43004e-16
+DEAL:3d::Norm of difference: 4.39080e-15
DEAL:3d::
DEAL:3d::Testing Gradients
-DEAL:3d::Norm of difference: 8.38689e-16
+DEAL:3d::Norm of difference: 1.34515e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 8.84868e-16
+DEAL:3d::Norm of difference: 5.97286e-15
DEAL:3d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
DEAL:3d::Number of cells: 64
DEAL:3d::Number of degrees of freedom: 5616
DEAL:3d::
DEAL:3d::Testing Values
-DEAL:3d::Norm of difference: 1.10927e-15
+DEAL:3d::Norm of difference: 6.20876e-15
DEAL:3d::
DEAL:3d::Testing Gradients
-DEAL:3d::Norm of difference: 1.50806e-15
+DEAL:3d::Norm of difference: 1.81453e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 1.72044e-15
+DEAL:3d::Norm of difference: 6.56897e-15
DEAL:3d::
DEAL:2d::Number of degrees of freedom: 40
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 3.77395e-16
+DEAL:2d::Norm of difference: 4.65353e-15
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 7.00345e-16
+DEAL:2d::Norm of difference: 1.63414e-15
DEAL:2d::
DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 4
DEAL:2d::Number of degrees of freedom: 84
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 7.63491e-16
+DEAL:2d::Norm of difference: 1.31449e-14
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 7.46533e-16
+DEAL:2d::Norm of difference: 2.68752e-15
DEAL:2d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 8
DEAL:3d::Number of degrees of freedom: 240
DEAL:3d::
DEAL:3d::Testing Values and Gradients
-DEAL:3d::Norm of difference: 1.22886e-15
+DEAL:3d::Norm of difference: 1.07664e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 1.16405e-15
+DEAL:3d::Norm of difference: 6.56100e-15
DEAL:3d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
DEAL:3d::Number of cells: 8
DEAL:3d::Number of degrees of freedom: 756
DEAL:3d::
DEAL:3d::Testing Values and Gradients
-DEAL:3d::Norm of difference: 2.06049e-15
+DEAL:3d::Norm of difference: 1.12673e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 1.23490e-15
+DEAL:3d::Norm of difference: 7.10070e-15
DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test is the same as matrix_vector_rt_01.cc but with non-affine cells in
+// standard orientation.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/manifold_lib.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_common.h"
+
+// This class is taken from
+// https://github.com/exadg/exadg/blob/master/include/exadg/grid/deformed_cube_manifold.h
+template <int dim>
+class DeformedCubeManifold : public dealii::ChartManifold<dim, dim, dim>
+{
+public:
+ DeformedCubeManifold(double const left,
+ double const right,
+ double const deformation,
+ unsigned int const frequency = 1)
+ : left(left)
+ , right(right)
+ , deformation(deformation)
+ , frequency(frequency)
+ {}
+
+ dealii::Point<dim>
+ push_forward(dealii::Point<dim> const &chart_point) const override
+ {
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= std::sin(frequency * dealii::numbers::PI *
+ (chart_point(d) - left) / (right - left));
+ dealii::Point<dim> space_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ space_point(d) = chart_point(d) + sinval;
+ return space_point;
+ }
+
+ dealii::Point<dim>
+ pull_back(dealii::Point<dim> const &space_point) const override
+ {
+ dealii::Point<dim> x = space_point;
+ dealii::Point<dim> one;
+ for (unsigned int d = 0; d < dim; ++d)
+ one(d) = 1.;
+
+ // Newton iteration to solve the nonlinear equation given by the point
+ dealii::Tensor<1, dim> sinvals;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ (right - left));
+
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ dealii::Tensor<1, dim> residual = space_point - x - sinval * one;
+ unsigned int its = 0;
+ while (residual.norm() > 1e-12 && its < 100)
+ {
+ dealii::Tensor<2, dim> jacobian;
+ for (unsigned int d = 0; d < dim; ++d)
+ jacobian[d][d] = 1.;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ double sinval_der = deformation * frequency / (right - left) *
+ dealii::numbers::PI *
+ std::cos(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+ for (unsigned int e = 0; e < dim; ++e)
+ if (e != d)
+ sinval_der *= sinvals[e];
+ for (unsigned int e = 0; e < dim; ++e)
+ jacobian[e][d] += sinval_der;
+ }
+
+ x += invert(jacobian) * residual;
+
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+
+ sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ residual = space_point - x - sinval * one;
+ ++its;
+ }
+ AssertThrow(residual.norm() < 1e-12,
+ dealii::ExcMessage("Newton for point did not converge."));
+ return x;
+ }
+
+ std::unique_ptr<dealii::Manifold<dim>>
+ clone() const override
+ {
+ return std::make_unique<DeformedCubeManifold<dim>>(left,
+ right,
+ deformation,
+ frequency);
+ }
+
+private:
+ double const left;
+ double const right;
+ double const deformation;
+ unsigned int const frequency;
+};
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+ unsigned int const frequency = 2;
+ double const deformation = 0.05;
+ static DeformedCubeManifold<dim> manifold(0.0, 1.0, deformation, frequency);
+ tria.set_all_manifold_ids(1);
+ tria.set_manifold(1, manifold);
+
+ std::vector<bool> vertex_touched(tria.n_vertices(), false);
+
+ for (auto cell : tria.cell_iterators())
+ {
+ for (auto const &v : cell->vertex_indices())
+ {
+ if (vertex_touched[cell->vertex_index(v)] == false)
+ {
+ Point<dim> &vertex = cell->vertex(v);
+ Point<dim> new_point = manifold.push_forward(vertex);
+ vertex = new_point;
+ vertex_touched[cell->vertex_index(v)] = true;
+ }
+ }
+ }
+
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ constraints.close();
+
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
+}
--- /dev/null
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 144
+DEAL:2d::
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 2.32427e-15
+DEAL:2d::
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 1.00134e-14
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 1.67249e-15
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 312
+DEAL:2d::
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 5.70018e-16
+DEAL:2d::
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 1.04689e-14
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 4.48610e-16
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 64
+DEAL:3d::Number of degrees of freedom: 1728
+DEAL:3d::
+DEAL:3d::Testing Values
+DEAL:3d::Norm of difference: 3.30534e-15
+DEAL:3d::
+DEAL:3d::Testing Gradients
+DEAL:3d::Norm of difference: 1.56133e-14
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 4.45455e-15
+DEAL:3d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Number of cells: 64
+DEAL:3d::Number of degrees of freedom: 5616
+DEAL:3d::
+DEAL:3d::Testing Values
+DEAL:3d::Norm of difference: 3.08613e-15
+DEAL:3d::
+DEAL:3d::Testing Gradients
+DEAL:3d::Norm of difference: 1.38898e-14
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 4.73268e-15
+DEAL:3d::
DEAL:2d::Number of degrees of freedom: 40
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 2.36910e-16
+DEAL:2d::Norm of difference: 7.46625e-15
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 4.52039e-16
+DEAL:2d::Norm of difference: 2.58308e-15
DEAL:2d::
DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 4
DEAL:2d::Number of degrees of freedom: 84
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 7.31531e-15
+DEAL:2d::Norm of difference: 7.18468e-15
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 2.97155e-15
+DEAL:2d::Norm of difference: 8.17176e-15
DEAL:2d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 8
DEAL:3d::Number of degrees of freedom: 240
DEAL:3d::
DEAL:3d::Testing Values and Gradients
-DEAL:3d::Norm of difference: 3.73456e-15
+DEAL:3d::Norm of difference: 1.38066e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 4.51598e-15
+DEAL:3d::Norm of difference: 4.05280e-15
DEAL:3d::
DEAL:2d::Number of degrees of freedom: 40
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 6.92435e-16
+DEAL:2d::Norm of difference: 1.45411e-14
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 2.84791e-16
+DEAL:2d::Norm of difference: 2.63432e-15
DEAL:2d::
DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
DEAL:2d::Number of cells: 4
DEAL:2d::Number of degrees of freedom: 84
DEAL:2d::
DEAL:2d::Testing Values and Gradients
-DEAL:2d::Norm of difference: 4.97871e-15
+DEAL:2d::Norm of difference: 1.52206e-14
DEAL:2d::
DEAL:2d::Testing Divergence
-DEAL:2d::Norm of difference: 4.05232e-15
+DEAL:2d::Norm of difference: 6.67882e-15
DEAL:2d::
DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
DEAL:3d::Number of cells: 8
DEAL:3d::Number of degrees of freedom: 240
DEAL:3d::
DEAL:3d::Testing Values and Gradients
-DEAL:3d::Norm of difference: 2.96549e-15
+DEAL:3d::Norm of difference: 2.58421e-14
DEAL:3d::
DEAL:3d::Testing Divergence
-DEAL:3d::Norm of difference: 4.91752e-15
+DEAL:3d::Norm of difference: 7.91079e-15
DEAL:3d::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// This test it the same as matrix_vector_rt_face_01.cc but with non-affine
+// cells in standard orientation.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/grid/manifold_lib.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_face_common.h"
+
+
+// This class is taken from
+// https://github.com/exadg/exadg/blob/master/include/exadg/grid/deformed_cube_manifold.h
+template <int dim>
+class DeformedCubeManifold : public dealii::ChartManifold<dim, dim, dim>
+{
+public:
+ DeformedCubeManifold(double const left,
+ double const right,
+ double const deformation,
+ unsigned int const frequency = 1)
+ : left(left)
+ , right(right)
+ , deformation(deformation)
+ , frequency(frequency)
+ {}
+
+ dealii::Point<dim>
+ push_forward(dealii::Point<dim> const &chart_point) const override
+ {
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= std::sin(frequency * dealii::numbers::PI *
+ (chart_point(d) - left) / (right - left));
+ dealii::Point<dim> space_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ space_point(d) = chart_point(d) + sinval;
+ return space_point;
+ }
+
+ dealii::Point<dim>
+ pull_back(dealii::Point<dim> const &space_point) const override
+ {
+ dealii::Point<dim> x = space_point;
+ dealii::Point<dim> one;
+ for (unsigned int d = 0; d < dim; ++d)
+ one(d) = 1.;
+
+ // Newton iteration to solve the nonlinear equation given by the point
+ dealii::Tensor<1, dim> sinvals;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI * (x(d) - left) /
+ (right - left));
+
+ double sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ dealii::Tensor<1, dim> residual = space_point - x - sinval * one;
+ unsigned int its = 0;
+ while (residual.norm() > 1e-12 && its < 100)
+ {
+ dealii::Tensor<2, dim> jacobian;
+ for (unsigned int d = 0; d < dim; ++d)
+ jacobian[d][d] = 1.;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ double sinval_der = deformation * frequency / (right - left) *
+ dealii::numbers::PI *
+ std::cos(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+ for (unsigned int e = 0; e < dim; ++e)
+ if (e != d)
+ sinval_der *= sinvals[e];
+ for (unsigned int e = 0; e < dim; ++e)
+ jacobian[e][d] += sinval_der;
+ }
+
+ x += invert(jacobian) * residual;
+
+ for (unsigned int d = 0; d < dim; ++d)
+ sinvals[d] = std::sin(frequency * dealii::numbers::PI *
+ (x(d) - left) / (right - left));
+
+ sinval = deformation;
+ for (unsigned int d = 0; d < dim; ++d)
+ sinval *= sinvals[d];
+ residual = space_point - x - sinval * one;
+ ++its;
+ }
+ AssertThrow(residual.norm() < 1e-12,
+ dealii::ExcMessage("Newton for point did not converge."));
+ return x;
+ }
+
+ std::unique_ptr<dealii::Manifold<dim>>
+ clone() const override
+ {
+ return std::make_unique<DeformedCubeManifold<dim>>(left,
+ right,
+ deformation,
+ frequency);
+ }
+
+private:
+ double const left;
+ double const right;
+ double const deformation;
+ unsigned int const frequency;
+};
+
+template <int dim, int fe_degree>
+void
+test()
+{
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(2);
+ unsigned int const frequency = 2;
+ double const deformation = 0.05;
+ static DeformedCubeManifold<dim> manifold(0.0, 1.0, deformation, frequency);
+ tria.set_all_manifold_ids(1);
+ tria.set_manifold(1, manifold);
+
+ std::vector<bool> vertex_touched(tria.n_vertices(), false);
+
+ for (auto cell : tria.cell_iterators())
+ {
+ for (auto const &v : cell->vertex_indices())
+ {
+ if (vertex_touched[cell->vertex_index(v)] == false)
+ {
+ Point<dim> &vertex = cell->vertex(v);
+ Point<dim> new_point = manifold.push_forward(vertex);
+ vertex = new_point;
+ vertex_touched[cell->vertex_index(v)] = true;
+ }
+ }
+ }
+
+
+ FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(fe);
+
+ AffineConstraints<double> constraints;
+ DoFTools::make_hanging_node_constraints(dof, constraints);
+ constraints.close();
+
+ deallog << "Using " << dof.get_fe().get_name() << std::endl;
+ deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+ << std::endl;
+ deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+ << std::endl;
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::values);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::gradients);
+ do_test<dim, fe_degree, double>(dof, constraints, TestType::divergence);
+}
--- /dev/null
+
+DEAL:2d::Using FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 144
+DEAL:2d::
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 7.70271e-15
+DEAL:2d::
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 2.06943e-14
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 5.70616e-15
+DEAL:2d::
+DEAL:2d::Using FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 312
+DEAL:2d::
+DEAL:2d::Testing Values
+DEAL:2d::Norm of difference: 6.01369e-15
+DEAL:2d::
+DEAL:2d::Testing Gradients
+DEAL:2d::Norm of difference: 3.30476e-14
+DEAL:2d::
+DEAL:2d::Testing Divergence
+DEAL:2d::Norm of difference: 3.09657e-14
+DEAL:2d::
+DEAL:3d::Using FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 64
+DEAL:3d::Number of degrees of freedom: 1728
+DEAL:3d::
+DEAL:3d::Testing Values
+DEAL:3d::Norm of difference: 1.68332e-14
+DEAL:3d::
+DEAL:3d::Testing Gradients
+DEAL:3d::Norm of difference: 3.15188e-14
+DEAL:3d::
+DEAL:3d::Testing Divergence
+DEAL:3d::Norm of difference: 1.04375e-14
+DEAL:3d::
projects / dealii.git / commitdiff