From: Niklas Wik Date: Fri, 3 Jun 2022 14:22:20 +0000 (+0200) Subject: Non-affine matrix-free Piola transform X-Git-Tag: v9.5.0-rc1~1081^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d8537e4fac2531d1f30d5583f2468989c03cae44;p=dealii.git Non-affine matrix-free Piola transform --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 5e93016525..730eb0bc18 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3022,6 +3022,10 @@ inline FEEvaluationBasemapped_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(); } @@ -3086,6 +3090,10 @@ operator=(const FEEvaluationBasemapped_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(); } @@ -5865,9 +5873,99 @@ FEEvaluationAccess:: 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> &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; } @@ -5914,14 +6012,17 @@ FEEvaluationAccess:: 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; @@ -5929,11 +6030,6 @@ FEEvaluationAccess:: divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det; } - else - { - // General cell - Assert(false, ExcNotImplemented()); - } } else { @@ -6192,15 +6288,109 @@ FEEvaluationAccess:: 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> &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 @@ -6258,37 +6448,36 @@ FEEvaluationAccess:: 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 { @@ -7193,6 +7382,8 @@ FEEvaluationJ_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); @@ -7265,7 +7456,9 @@ FEEvaluationcell_type <= internal::MatrixFreeFunctions::GeometryType::affine) { @@ -7278,6 +7471,9 @@ FEEvaluationmapping_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); } @@ -7292,6 +7488,10 @@ FEEvaluationmapping_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); } @@ -7300,7 +7500,9 @@ FEEvaluationjacobian = 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) @@ -7340,6 +7542,14 @@ FEEvaluationmapping_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] = @@ -7381,6 +7591,15 @@ FEEvaluationmapping_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 <= @@ -7611,6 +7830,12 @@ FEEvaluationdata->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:: @@ -7900,6 +8125,26 @@ FEEvaluationdata->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:: @@ -8137,6 +8382,11 @@ FEFaceEvaluationjacobian_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) { @@ -8294,6 +8544,11 @@ FEFaceEvaluationjacobian_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] @@ -8428,6 +8683,12 @@ FEFaceEvaluationdata->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:: template run(n_components, @@ -8565,6 +8826,26 @@ FEFaceEvaluationdata->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:: template run(n_components, diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index abd3f0e70c..230231f44a 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -756,6 +756,13 @@ protected: 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 @@ -810,6 +817,20 @@ protected: */ 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 @@ -980,6 +1001,12 @@ protected: internal::MatrixFreeFunctions::MappingDataOnTheFly> 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 @@ -1039,6 +1066,7 @@ inline FEEvaluationData::FEEvaluationData( , quadrature_points(nullptr) , jacobian(nullptr) , jacobian_gradients(nullptr) + , jacobian_gradients_non_inverse(nullptr) , J_value(nullptr) , normal_vectors(nullptr) , normal_x_jacobian(nullptr) @@ -1063,6 +1091,7 @@ inline FEEvaluationData::FEEvaluationData( internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) , subface_index(0) , cell_type(internal::MatrixFreeFunctions::general) + , divergence_is_requested(false) {} @@ -1088,6 +1117,7 @@ inline FEEvaluationData::FEEvaluationData( , quadrature_points(nullptr) , jacobian(nullptr) , jacobian_gradients(nullptr) + , jacobian_gradients_non_inverse(nullptr) , J_value(nullptr) , normal_vectors(nullptr) , normal_x_jacobian(nullptr) @@ -1098,12 +1128,16 @@ inline FEEvaluationData::FEEvaluationData( , 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(); } @@ -1121,17 +1155,18 @@ FEEvaluationData::operator=(const FEEvaluationData &other) 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; @@ -1151,10 +1186,11 @@ FEEvaluationData::operator=(const FEEvaluationData &other) 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; } @@ -1179,7 +1215,7 @@ FEEvaluationData::set_data_pointers( 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; @@ -1190,14 +1226,18 @@ FEEvaluationData::set_data_pointers( // 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); } diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index c63f563e41..81b1849743 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -796,10 +796,31 @@ namespace internal 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); + } } } @@ -943,6 +964,12 @@ namespace internal 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() + @@ -1263,6 +1290,28 @@ namespace internal 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]); } } } @@ -1368,8 +1417,12 @@ namespace internal 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()); @@ -2152,6 +2205,29 @@ namespace internal 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]); } }; @@ -2440,6 +2516,10 @@ namespace internal 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()); @@ -2678,7 +2758,10 @@ namespace internal 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) { @@ -2810,6 +2893,8 @@ namespace internal { 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); @@ -3008,8 +3093,14 @@ namespace internal face_data_by_cells[my_q].normals_times_jacobians[1].resize_fast( storage_length * GeometryInfo::faces_per_cell); if (update_flags & update_jacobian_grads) - face_data_by_cells[my_q].jacobian_gradients[0].resize_fast( - storage_length * GeometryInfo::faces_per_cell); + { + face_data_by_cells[my_q].jacobian_gradients[0].resize_fast( + storage_length * GeometryInfo::faces_per_cell); + face_data_by_cells[my_q] + .jacobian_gradients_non_inverse[0] + .resize_fast(storage_length * + GeometryInfo::faces_per_cell); + } if (update_flags & update_quadrature_points) face_data_by_cells[my_q].quadrature_points.resize_fast( diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h index 67fd4ba04c..4cf01ba24b 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.h +++ b/include/deal.II/matrix_free/mapping_info_storage.h @@ -259,6 +259,27 @@ namespace internal 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 diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h index e34d09d4e4..fd468200c1 100644 --- a/include/deal.II/matrix_free/mapping_info_storage.templates.h +++ b/include/deal.II/matrix_free/mapping_info_storage.templates.h @@ -112,6 +112,7 @@ namespace internal { jacobians[i].clear(); jacobian_gradients[i].clear(); + jacobian_gradients_non_inverse[i].clear(); normals_times_jacobians[i].clear(); } quadrature_point_offsets.clear(); @@ -183,6 +184,10 @@ namespace internal 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) + @@ -218,7 +223,11 @@ namespace internal 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); diff --git a/tests/matrix_free/matrix_vector_rt_01.output b/tests/matrix_free/matrix_vector_rt_01.output index 4c1b0e9e9a..1264f4e3bc 100644 --- a/tests/matrix_free/matrix_vector_rt_01.output +++ b/tests/matrix_free/matrix_vector_rt_01.output @@ -4,50 +4,50 @@ 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: 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:: diff --git a/tests/matrix_free/matrix_vector_rt_02.output b/tests/matrix_free/matrix_vector_rt_02.output index ecae4aa430..fa646f4404 100644 --- a/tests/matrix_free/matrix_vector_rt_02.output +++ b/tests/matrix_free/matrix_vector_rt_02.output @@ -4,38 +4,38 @@ DEAL:2d::Number of cells: 4 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:: diff --git a/tests/matrix_free/matrix_vector_rt_04.cc b/tests/matrix_free/matrix_vector_rt_04.cc new file mode 100644 index 0000000000..567f829ce9 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_04.cc @@ -0,0 +1,172 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#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 +class DeformedCubeManifold : public dealii::ChartManifold +{ +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 + push_forward(dealii::Point 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 space_point; + for (unsigned int d = 0; d < dim; ++d) + space_point(d) = chart_point(d) + sinval; + return space_point; + } + + dealii::Point + pull_back(dealii::Point const &space_point) const override + { + dealii::Point x = space_point; + dealii::Point 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> + clone() const override + { + return std::make_unique>(left, + right, + deformation, + frequency); + } + +private: + double const left; + double const right; + double const deformation; + unsigned int const frequency; +}; + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + unsigned int const frequency = 2; + double const deformation = 0.05; + static DeformedCubeManifold manifold(0.0, 1.0, deformation, frequency); + tria.set_all_manifold_ids(1); + tria.set_manifold(1, manifold); + + std::vector 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 &vertex = cell->vertex(v); + Point new_point = manifold.push_forward(vertex); + vertex = new_point; + vertex_touched[cell->vertex_index(v)] = true; + } + } + } + + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints 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(dof, constraints, TestType::values); + do_test(dof, constraints, TestType::gradients); + do_test(dof, constraints, TestType::divergence); +} diff --git a/tests/matrix_free/matrix_vector_rt_04.output b/tests/matrix_free/matrix_vector_rt_04.output new file mode 100644 index 0000000000..eec9900828 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_04.output @@ -0,0 +1,53 @@ + +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:: diff --git a/tests/matrix_free/matrix_vector_rt_face_01.output b/tests/matrix_free/matrix_vector_rt_face_01.output index 2d1aea3333..2046705869 100644 --- a/tests/matrix_free/matrix_vector_rt_face_01.output +++ b/tests/matrix_free/matrix_vector_rt_face_01.output @@ -4,28 +4,28 @@ DEAL:2d::Number of cells: 4 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:: diff --git a/tests/matrix_free/matrix_vector_rt_face_02.output b/tests/matrix_free/matrix_vector_rt_face_02.output index d107fc460f..2504d6d8eb 100644 --- a/tests/matrix_free/matrix_vector_rt_face_02.output +++ b/tests/matrix_free/matrix_vector_rt_face_02.output @@ -4,28 +4,28 @@ DEAL:2d::Number of cells: 4 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:: diff --git a/tests/matrix_free/matrix_vector_rt_face_04.cc b/tests/matrix_free/matrix_vector_rt_face_04.cc new file mode 100644 index 0000000000..cc4ee5ada1 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_04.cc @@ -0,0 +1,173 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#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 +class DeformedCubeManifold : public dealii::ChartManifold +{ +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 + push_forward(dealii::Point 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 space_point; + for (unsigned int d = 0; d < dim; ++d) + space_point(d) = chart_point(d) + sinval; + return space_point; + } + + dealii::Point + pull_back(dealii::Point const &space_point) const override + { + dealii::Point x = space_point; + dealii::Point 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> + clone() const override + { + return std::make_unique>(left, + right, + deformation, + frequency); + } + +private: + double const left; + double const right; + double const deformation; + unsigned int const frequency; +}; + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + unsigned int const frequency = 2; + double const deformation = 0.05; + static DeformedCubeManifold manifold(0.0, 1.0, deformation, frequency); + tria.set_all_manifold_ids(1); + tria.set_manifold(1, manifold); + + std::vector 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 &vertex = cell->vertex(v); + Point new_point = manifold.push_forward(vertex); + vertex = new_point; + vertex_touched[cell->vertex_index(v)] = true; + } + } + } + + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints 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(dof, constraints, TestType::values); + do_test(dof, constraints, TestType::gradients); + do_test(dof, constraints, TestType::divergence); +} diff --git a/tests/matrix_free/matrix_vector_rt_face_04.output b/tests/matrix_free/matrix_vector_rt_face_04.output new file mode 100644 index 0000000000..63b6507137 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_04.output @@ -0,0 +1,40 @@ + +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::