From fc53f901898a2b8b7add421293ab4ebf8aaaa04c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 16 Jan 2017 12:06:18 -0700 Subject: [PATCH] Better track which UpdateFlags are required for which operation. --- source/fe/mapping_q_generic.cc | 166 ++++++++++++++++++--------------- 1 file changed, 91 insertions(+), 75 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 333dae05e7..1df72d9163 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2000 - 2016 by the deal.II authors +// Copyright (C) 2000 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -848,7 +848,11 @@ initialize_face (const UpdateFlags update_flags, if (dim > 1) { - if (this->update_each & update_boundary_forms) + if (this->update_each & (update_boundary_forms | + update_normal_vectors | + update_jacobians | + update_JxW_values | + update_inverse_jacobians)) { aux.resize (dim-1, std::vector > (n_original_q_points)); @@ -1961,10 +1965,12 @@ MappingQGeneric::requires_update_flags (const UpdateFlags in) cons // contravariant or the Piola transforms, we add the JxW values // to the list of flags to be updated for each cell. if (out & update_contravariant_transformation) - out |= update_JxW_values; + out |= update_volume_elements; + // the same is true when computing normal vectors: they require + // the determinant of the Jacobian if (out & update_normal_vectors) - out |= update_JxW_values; + out |= update_volume_elements; } return out; @@ -2778,14 +2784,24 @@ namespace internal { const UpdateFlags update_flags = data.update_each; - if (update_flags & update_boundary_forms) + if (update_flags & (update_boundary_forms | + update_normal_vectors | + update_jacobians | + update_JxW_values | + update_inverse_jacobians)) { - AssertDimension (output_data.boundary_forms.size(), n_q_points); + if (update_flags & update_boundary_forms) + AssertDimension (output_data.boundary_forms.size(), n_q_points); if (update_flags & update_normal_vectors) AssertDimension (output_data.normal_vectors.size(), n_q_points); if (update_flags & update_JxW_values) AssertDimension (output_data.JxW_values.size(), n_q_points); + Assert (data.aux.size() >= dim-1, ExcInternalError()); + + // first compute some common data that is used for evaluating + // all of the flags below + // map the unit tangentials to the real cell. checking for d!=dim-1 // eliminates compiler warnings regarding unsigned int expressions < // 0. @@ -2804,91 +2820,91 @@ namespace internal make_array_view(data.aux[d])); } - // if dim==spacedim, we can use the unit tangentials to compute the - // boundary form by simply taking the cross product - if (dim == spacedim) + if (update_flags & update_boundary_forms) { - for (unsigned int i=0; i DX_t = - data.contravariant[point].transpose(); + if (dim==2) + { + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); - // then compute the face normal from the face tangent - // and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } } } } - if (update_flags & (update_normal_vectors - | update_JxW_values)) + if (update_flags & update_JxW_values) for (unsigned int i=0; i::subface_ratio( - cell->subface_case(face_no), subface_no); - output_data.JxW_values[i] *= area_ratio; - } + if (subface_no != numbers::invalid_unsigned_int) + { + const double area_ratio = GeometryInfo::subface_ratio(cell->subface_case(face_no), + subface_no); + output_data.JxW_values[i] *= area_ratio; } - - if (update_flags & update_normal_vectors) - output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / - output_data.boundary_forms[i].norm()); } + if (update_flags & update_normal_vectors) + for (unsigned int i=0; i(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); + if (update_flags & update_jacobians) for (unsigned int point=0; point