From 4d4a12fecb40ed5080fc35569b5e2ccd34a2f221 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 24 Nov 2015 11:18:05 -0600 Subject: [PATCH] Compute higher order derivative mappings in MappingCartesian. These fields were previously only computed for cells, but not for faces or subfaces. This happened to work by accident because these fields are initialized by zeros. But I'm changing this so that they are initialized by signaling NaNs, and so we need to be explicit about this. --- source/fe/mapping_cartesian.cc | 126 ++++++++++++++++++++------------- 1 file changed, 77 insertions(+), 49 deletions(-) diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 72fe571824..f8f16e989b 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -339,10 +339,8 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, const typename Mapping::InternalDataBase &internal_data, internal::FEValues::MappingRelatedData &output_data) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible + // convert data object to internal data for this class. fails with + // an exception if that is not possible Assert (dynamic_cast (&internal_data) != 0, ExcInternalError()); const InternalData &data = static_cast (internal_data); @@ -353,11 +351,8 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, output_data.quadrature_points, dummy); - // compute Jacobian - // determinant. all values are - // equal and are the product of the - // local lengths in each coordinate - // direction + // compute Jacobian determinant. all values are equal and are the + // product of the local lengths in each coordinate direction if (data.update_each & (update_JxW_values | update_volume_elements)) if (cell_similarity != CellSimilarity::translation) { @@ -369,8 +364,8 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, for (unsigned int i=0; i::cell_iterator &cell, for (unsigned int j=0; j(); + output_data.jacobian_grads[i] = DerivativeForm<2,dim,spacedim>(); if (data.update_each & update_jacobian_pushed_forward_grads) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.jacobian_pushed_forward_grads[i] = Tensor<3,spacedim>(); - // "compute" the hessian of the Jacobian - // at the quadrature points, which are all - // also zero of course + // "compute" the hessian of the Jacobian at the quadrature points, + // which are all also zero of course if (data.update_each & update_jacobian_2nd_derivatives) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.jacobian_2nd_derivatives[i] = DerivativeForm<3,dim,spacedim>(); if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.jacobian_pushed_forward_2nd_derivatives[i] = Tensor<4,spacedim>(); if (data.update_each & update_jacobian_3rd_derivatives) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.jacobian_3rd_derivatives[i] = DerivativeForm<4,dim,spacedim>(); if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.jacobian_pushed_forward_3rd_derivatives[i] = Tensor<5,spacedim>(); - // "compute" inverse Jacobian at the - // quadrature points, which are all - // the same + // "compute" inverse Jacobian at the quadrature points, which are + // all the same if (data.update_each & update_inverse_jacobians) if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); + output_data.inverse_jacobians[i] = Tensor<2,dim>(); for (unsigned int j=0; j::cell_iterator & output_data.quadrature_points, output_data.normal_vectors); - // first compute Jacobian - // determinant, which is simply the - // product of the local lengths - // since the jacobian is diagonal + // first compute Jacobian determinant, which is simply the product + // of the local lengths since the jacobian is diagonal double J = 1.; for (unsigned int d=0; d::unit_normal_direction[face_no]) @@ -488,6 +478,30 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & output_data.jacobians[i][d][d] = data.cell_extents[d]; } + if (data.update_each & update_jacobian_grads) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_grads) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_2nd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_3rd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) + for (unsigned int i=0; i(); + if (data.update_each & update_inverse_jacobians) for (unsigned int i=0; i::cell_iterato const typename Mapping::InternalDataBase &internal_data, internal::FEValues::MappingRelatedData &output_data) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible + // convert data object to internal data for this class. fails with + // an exception if that is not possible Assert (dynamic_cast (&internal_data) != 0, ExcInternalError()); const InternalData &data = static_cast (internal_data); @@ -521,10 +533,8 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato output_data.quadrature_points, output_data.normal_vectors); - // first compute Jacobian - // determinant, which is simply the - // product of the local lengths - // since the jacobian is diagonal + // first compute Jacobian determinant, which is simply the product + // of the local lengths since the jacobian is diagonal double J = 1.; for (unsigned int d=0; d::unit_normal_direction[face_no]) @@ -532,16 +542,10 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato if (data.update_each & update_JxW_values) { - // Here, - // cell->face(face_no)->n_children() - // would be the right choice, - // but unfortunately the - // current function is also - // called for faces without - // children (see - // tests/fe/mapping.cc). Add - // following switch to avoid - // diffs in tests/fe/mapping.OK + // Here, cell->face(face_no)->n_children() would be the right + // choice, but unfortunately the current function is also called + // for faces without children (see tests/fe/mapping.cc). Add + // following switch to avoid diffs in tests/fe/mapping.OK const unsigned int n_subfaces= cell->face(face_no)->has_children() ? cell->face(face_no)->n_children() : @@ -570,6 +574,30 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato output_data.jacobians[i][d][d] = data.cell_extents[d]; } + if (data.update_each & update_jacobian_grads) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_grads) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_2nd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_3rd_derivatives) + for (unsigned int i=0; i(); + + if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) + for (unsigned int i=0; i(); + if (data.update_each & update_inverse_jacobians) for (unsigned int i=0; i