From: Wolfgang Bangerth Date: Mon, 11 Jul 2016 12:46:20 +0000 (-0500) Subject: Fix a bug in MappingQGeneric::reinit(). X-Git-Tag: v8.5.0-rc1~886^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2a7f38f57cbfdf1aaabc080b4f9a209611c8db29;p=dealii.git Fix a bug in MappingQGeneric::reinit(). Specifically, we thought that we were on the same cell as the last call to FEValues::reinit() and so did not re-compute information -- but the vertices of the triangulation had been moved in between, and consequently the previously computed data was stale. --- diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 42a5742397..1651939d7c 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -2542,20 +2542,18 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, const unsigned int n_q_points=quadrature.size(); - // if necessary, recompute the support points of the transformation of this cell - // (note that we need to first check the triangulation pointer, since otherwise - // the second test might trigger an exception if the triangulations are not the - // same) - if ((data.mapping_support_points.size() == 0) - || - (&cell->get_triangulation() != - &data.cell_of_current_support_points->get_triangulation()) - || - (cell != data.cell_of_current_support_points)) - { - data.mapping_support_points = this->compute_mapping_support_points(cell); - data.cell_of_current_support_points = cell; - } + // recompute the support points of the transformation of this + // cell. we tried to be clever here in an earlier version of the + // library by checking whether the cell is the same as the one we + // had visited last, but it turns out to be difficult to determine + // that because a cell for the purposes of a mapping is + // characterized not just by its (triangulation, level, index) + // triple, but also by the locations of its vertices, the manifold + // object attached to the cell and all of its bounding faces/edges, + // etc. to reliably test that the "cell" we are on is, therefore, + // not easily done + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), data,