]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix a bug in MappingQGeneric::reinit().
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 11 Jul 2016 12:46:20 +0000 (07:46 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 11 Jul 2016 12:46:20 +0000 (07:46 -0500)
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.

source/fe/mapping_q_generic.cc

index 42a57423970bc895ebba2a697e2caedf2677f481..1651939d7cbbe2782492a0ec62f9327b615eeca9 100644 (file)
@@ -2542,20 +2542,18 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::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<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
                                                   data,

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.