From f3377ce95257c20a405b50769991980f6d153332 Mon Sep 17 00:00:00 2001 From: Maien Hamed Date: Tue, 11 Aug 2015 14:31:26 +0200 Subject: [PATCH] Initialize scratch vectors for transformed values and gradients in get_data. Added entry to doc/news/changes.h reflecting analytical hessian computations --- doc/news/changes.h | 7 + include/deal.II/fe/fe_poly_tensor.h | 8 + source/fe/fe_poly_tensor.cc | 270 ++++++++++++++-------------- 3 files changed, 145 insertions(+), 140 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index c637d7b79a..2db0318d82 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -90,6 +90,13 @@ inconvenience this causes.
    +
  1. Improved: Some finite elements compute hessians analytically rather than + by finite differencing. Namely, these are finite elements that are subclasses + of FEPoly as well as FESystem with those as base elements. +
    + (Maien Hamed, 2015/08/01-2015/08/09) +
  2. +
  3. New: There is now a function Mapping::project_real_point_to_unit_point_on_face() that calls Mapping::transform_real_to_unit_cell() and then projects the result to a provided face. diff --git a/include/deal.II/fe/fe_poly_tensor.h b/include/deal.II/fe/fe_poly_tensor.h index 8793dbfe01..0c9c5ed9b3 100644 --- a/include/deal.II/fe/fe_poly_tensor.h +++ b/include/deal.II/fe/fe_poly_tensor.h @@ -239,6 +239,14 @@ protected: * point. */ std::vector< std::vector< DerivativeForm<1, dim, spacedim> > > shape_grads; + + /** + * Scratch arrays for intermediate computations + */ + mutable std::vector sign_change; + mutable std::vector > transformed_shape_values; + mutable std::vector > transformed_shape_grads; + mutable std::vector > untransformed_shape_grads; }; /** diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 5b50d5c30f..57bdd8b80b 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -305,6 +305,8 @@ FE_PolyTensor::get_data ( std::vector > grads(0); std::vector > grad_grads(0); + data->sign_change.resize (this->dofs_per_cell); + // initialize fields only if really // necessary. otherwise, don't // allocate memory @@ -314,6 +316,8 @@ FE_PolyTensor::get_data ( data->shape_values.resize (this->dofs_per_cell); for (unsigned int i=0; idofs_per_cell; ++i) data->shape_values[i].resize (n_q_points); + if (mapping_type != mapping_none) + data->transformed_shape_values.resize (n_q_points); } if (flags & update_gradients) @@ -322,6 +326,13 @@ FE_PolyTensor::get_data ( data->shape_grads.resize (this->dofs_per_cell); for (unsigned int i=0; idofs_per_cell; ++i) data->shape_grads[i].resize (n_q_points); + data->transformed_shape_grads.resize (n_q_points); + if ( (mapping_type == mapping_raviart_thomas) + || + (mapping_type == mapping_piola) + || + (mapping_type == mapping_nedelec)) + data->untransformed_shape_grads.resize(n_q_points); } // if second derivatives through @@ -420,12 +431,12 @@ fill_fe_values (const Mapping &mapping, // Compute eventual sign changes depending on the neighborhood // between two faces. - std::vector sign_change (this->dofs_per_cell, 1.0); + std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); + get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) @@ -444,59 +455,58 @@ fill_fe_values (const Mapping &mapping, || ((mapping_type == mapping_piola) || (mapping_type == mapping_raviart_thomas) || (mapping_type == mapping_nedelec)))) - switch (mapping_type) - { - case mapping_none: - { - for (unsigned int k=0; k > shape_values (n_q_points); - mapping.transform(fe_data.shape_values[i], - shape_values, mapping_internal, mapping_type); + case mapping_covariant: + case mapping_contravariant: + { + mapping.transform(fe_data.shape_values[i], + fe_data.transformed_shape_values, mapping_internal, mapping_type); - for (unsigned int k=0; k > shape_values (n_q_points); - mapping.transform(fe_data.shape_values[i], shape_values, - mapping_internal, mapping_piola); - for (unsigned int k=0; k > shape_values (n_q_points); - mapping.transform (fe_data.shape_values[i], shape_values, - mapping_internal, mapping_covariant); + case mapping_nedelec: + { + mapping.transform (fe_data.shape_values[i], fe_data.transformed_shape_values, + mapping_internal, mapping_covariant); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; + for (unsigned int k = 0; k < n_q_points; ++k) + for (unsigned int d = 0; d < dim; ++d) + output_data.shape_values(first+d,k) = fe_data.sign_change[i] * fe_data.transformed_shape_values[k][d]; - break; - } + break; + } - default: - Assert(false, ExcNotImplemented()); - } + default: + Assert(false, ExcNotImplemented()); + } + } // update gradients. apply the same logic as above if (flags & update_gradients @@ -507,57 +517,50 @@ fill_fe_values (const Mapping &mapping, || (mapping_type == mapping_nedelec)))) { - std::vector > shape_grads1 (n_q_points); - //std::vector > shape_grads2 (n_q_points); - switch (mapping_type) { case mapping_none: { - mapping.transform(fe_data.shape_grads[i], shape_grads1, + mapping.transform(fe_data.shape_grads[i], fe_data.transformed_shape_grads, mapping_internal, mapping_covariant); for (unsigned int k=0; k > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k &mapping, // (J^-T)Dv(J^-1) is the // value we want to have on // the real cell. - std::vector > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k &mapping, // Compute eventual sign changes depending // on the neighborhood between two faces. - std::vector sign_change (this->dofs_per_cell, 1.0); + std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); + get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { @@ -665,6 +668,7 @@ fill_fe_face_values (const Mapping &mapping, if (flags & update_values) { + VectorSlice< std::vector > > transformed_shape_values(fe_data.transformed_shape_values, offset, n_q_points); switch (mapping_type) { case mapping_none: @@ -680,38 +684,36 @@ fill_fe_face_values (const Mapping &mapping, { // Use auxiliary vector // for transformation - std::vector > shape_values (n_q_points); mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_type); + transformed_shape_values, mapping_internal, mapping_type); for (unsigned int k=0; k > shape_values (n_q_points); + mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_piola); + transformed_shape_values, mapping_internal, mapping_piola); for (unsigned int k=0; k > shape_values (n_q_points); mapping.transform (make_slice (fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_covariant); + transformed_shape_values, mapping_internal, mapping_covariant); for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; + output_data.shape_values(first+d,k) = fe_data.sign_change[i] * transformed_shape_values[k][d]; break; } @@ -723,56 +725,51 @@ fill_fe_face_values (const Mapping &mapping, if (flags & update_gradients) { - std::vector > shape_grads1 (n_q_points); - // std::vector > shape_grads2 (n_q_points); - - + VectorSlice< std::vector > > transformed_shape_grads (fe_data.transformed_shape_grads, offset, n_q_points); switch (mapping_type) { case mapping_none: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_internal, mapping_covariant); + transformed_shape_grads, mapping_internal, mapping_covariant); for (unsigned int k=0; k > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k &mapping, // (J^-T)Dv(J^-1) is the // value we want to have on // the real cell. - std::vector > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k &mapping, // Compute eventual sign changes depending // on the neighborhood between two faces. - std::vector sign_change (this->dofs_per_cell, 1.0); + std::fill( fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0 ); if (mapping_type == mapping_raviart_thomas) - get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + get_face_sign_change_rt (cell, this->dofs_per_face, fe_data.sign_change); else if (mapping_type == mapping_nedelec) - get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); + get_face_sign_change_nedelec (cell, this->dofs_per_face, fe_data.sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { @@ -876,6 +873,7 @@ fill_fe_subface_values (const Mapping &mapping, if (flags & update_values) { + VectorSlice< std::vector > > transformed_shape_values (fe_data.transformed_shape_values, offset, n_q_points); switch (mapping_type) { case mapping_none: @@ -891,37 +889,34 @@ fill_fe_subface_values (const Mapping &mapping, { // Use auxiliary vector for // transformation - std::vector > shape_values (n_q_points); mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_type); + transformed_shape_values, mapping_internal, mapping_type); for (unsigned int k=0; k > shape_values (n_q_points); mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_piola); + transformed_shape_values, mapping_internal, mapping_piola); for (unsigned int k=0; k > shape_values (n_q_points); mapping.transform (make_slice (fe_data.shape_values[i], offset, n_q_points), - shape_values, mapping_internal, mapping_covariant); + transformed_shape_values, mapping_internal, mapping_covariant); for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; + output_data.shape_values(first+d,k) = fe_data.sign_change[i] * transformed_shape_values[k][d]; break; } @@ -932,58 +927,54 @@ fill_fe_subface_values (const Mapping &mapping, if (flags & update_gradients) { - std::vector > shape_grads1 (n_q_points); - std::vector > shape_grads2 (n_q_points); - + VectorSlice< std::vector > > transformed_shape_grads (fe_data.transformed_shape_grads, offset, n_q_points); switch (mapping_type) { case mapping_none: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), - shape_grads1, mapping_internal, mapping_covariant); + transformed_shape_grads, mapping_internal, mapping_covariant); for (unsigned int k=0; k > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k &mapping, // (J^-T)Dv(J^-1) is the // value we want to have on // the real cell. - std::vector > input(fe_data.shape_grads[i].size()); - for (unsigned int k=0; k