From 874d89dbaf62321c6b6a21f468a8f5cd561b41fd Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 14 Jan 2021 18:08:32 +0100 Subject: [PATCH] Implemented 2nd order correction. --- source/fe/mapping_fe.cc | 32 +++++++++++++------------------- 1 file changed, 13 insertions(+), 19 deletions(-) diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 7d605bde60..7a3822b5e0 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -917,14 +917,12 @@ MappingFE::transform_real_to_unit_cell( Point p_unit; - Tensor<1, dim> grad_FT_residual; - unsigned int loop = 0; // This loop solves the following problem: - // grad_F^T residual + (grad_F^T grad_F + grad_F^T hesss_F dp) dp = 0 + // grad_F^T residual + (grad_F^T grad_F + grad_F^T hess_F^T dp) dp = 0 // where the term - // (grad_F^T hess_F dp) is approximated by (hess_F * residual) + // (grad_F^T hess_F dp) is approximated by (-hess_F * residual) // This is basically a second order approximation of Newton method, where the // Jacobian is corrected with a higher order term coming from the hessian. do @@ -933,24 +931,20 @@ MappingFE::transform_real_to_unit_cell( // Transpose of the gradient map DerivativeForm<1, spacedim, dim> grad_FT; + Tensor<1, dim> grad_FT_residual; Tensor<2, dim> corrected_metric_tensor; + DerivativeForm<2, spacedim, dim> hess_FT; - // [TODO] - // When 2nd derivatives are implemented, we'll uncomment this - // DerivativeForm<2, spacedim, dim> hess_FT; for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i) { mapped_point += support_points[i] * this->fe->shape_value(i, p_unit); - const auto grad_F_i = this->fe->shape_grad(i, p_unit); - // [TODO] - // When 2nd derivatives are implemented, we'll uncomment this - // auto hessian_k = this->fe->shape_grad_grad(i, p_unit); + const auto grad_F_i = this->fe->shape_grad(i, p_unit); + const auto hessian_F_i = this->fe->shape_grad_grad(i, p_unit); for (unsigned int j = 0; j < dim; ++j) { grad_FT[j] += grad_F_i[j] * support_points[i]; - // [TODO] - // When 2nd derivatives are implemented, we'll uncomment this - // hess_FT[j] += hessian_k[j] * support_points[i]; + for (unsigned int l = 0; l < dim; ++l) + hess_FT[j][l] += hessian_F_i[j][l] * support_points[i]; } } @@ -966,12 +960,12 @@ MappingFE::transform_real_to_unit_cell( break; // Now compute the (corrected) metric tensor - corrected_metric_tensor = -apply_transformation(grad_FT, grad_FT); - - // [TODO] - // When 2nd derivatives are implemented, we'll uncomment this - // corrected_metric_tensor += apply_transformation(hess_FT, residual); + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + corrected_metric_tensor[j][l] = + -grad_FT[j] * grad_FT[l] + hess_FT[j][l] * residual; + // And compute the update auto g_inverse = invert(corrected_metric_tensor); p_unit -= Point(g_inverse * grad_FT_residual); -- 2.39.5