From: bangerth Date: Wed, 7 Sep 2011 23:21:05 +0000 (+0000) Subject: Fix by Patrick Sodre: FEFieldFunction was not thread-safe due to the cell_hint cache... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a40aeff25ecf24b5bc016bcb9f8cfba9dd5eb8e5;p=dealii-svn.git Fix by Patrick Sodre: FEFieldFunction was not thread-safe due to the cell_hint cache (a mutable variable) that different threads kept stomping on. This is now fixed using a thread-local variable. git-svn-id: https://svn.dealii.org/trunk@24281 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c234605355..6adb15686d 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -283,6 +283,13 @@ and DoF handlers embedded in higher dimensional space have been added.

Specific improvements

    +
  1. Fixed: FEFieldFunction class was not thread-safe because it keeps a +cache on the side that was invalidated when different +threads kept pouncing on it. This is now fixed. +
    +(Patrick Sodré, 2011/09/07) + +
  2. New: There is now a function GridTools::volume computing the volume of a triangulation.
    diff --git a/deal.II/include/deal.II/numerics/fe_field_function.h b/deal.II/include/deal.II/numerics/fe_field_function.h index c57fca493c..151c5255e9 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.h @@ -19,6 +19,7 @@ #include #include #include +#include #include @@ -150,10 +151,10 @@ namespace Functions * points lie, you can tell * this object by calling this * function. This will speed - * things up a little. + * things up a little. */ void set_active_cell(typename DH::active_cell_iterator &newcell); - + /** * Get ONE vector value at the * given point. It is @@ -329,6 +330,12 @@ namespace Functions std::vector > &maps) const; private: + /** + * Typedef holding the local cell_hint. + */ + typedef Threads::ThreadLocalStorage < + typename DH::active_cell_iterator > cell_hint_t; + /** * Pointer to the dof handler. */ @@ -347,10 +354,9 @@ namespace Functions const Mapping & mapping; /** - * The current cell in which we - * are evaluating. + * The latest cell hint. */ - mutable typename DH::active_cell_iterator cell; + mutable cell_hint_t cell_hint; /** * Store the number of diff --git a/deal.II/include/deal.II/numerics/fe_field_function.templates.h b/deal.II/include/deal.II/numerics/fe_field_function.templates.h index 2a220a7e78..70dce9ed44 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.templates.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.templates.h @@ -31,9 +31,9 @@ namespace Functions dh(&mydh, "FEFieldFunction"), data_vector(myv), mapping(mymapping), + cell_hint(dh->end()), n_components(mydh.get_fe().n_components()) { - cell = dh->begin_active(); } @@ -43,17 +43,20 @@ namespace Functions FEFieldFunction:: set_active_cell(typename DH::active_cell_iterator &newcell) { - cell = newcell; + cell_hint.get() = newcell; } - + template void FEFieldFunction::vector_value (const Point &p, Vector &values) const { Assert (values.size() == n_components, ExcDimensionMismatch(values.size(), n_components)); + typename DH::active_cell_iterator cell = cell_hint.get(); + if (cell == dh->end()) + cell = dh->begin_active(); Point qp = mapping.transform_real_to_unit_cell(cell, p); // Check if we already have all we need @@ -64,6 +67,8 @@ namespace Functions cell = my_pair.first; qp = my_pair.second; } + + cell_hint.get() = cell; // Now we can find out about the point Quadrature quad(qp); @@ -96,6 +101,9 @@ namespace Functions { Assert (gradients.size() == n_components, ExcDimensionMismatch(gradients.size(), n_components)); + typename DH::active_cell_iterator cell = cell_hint.get(); + if(cell == dh->end()) + cell = dh->begin_active(); Point qp = mapping.transform_real_to_unit_cell(cell, p); // Check if we already have all we need @@ -106,7 +114,9 @@ namespace Functions cell = my_pair.first; qp = my_pair.second; } - + + cell_hint.get() = cell; + // Now we can find out about the point Quadrature quad(qp); FEValues fe_v(mapping, dh->get_fe(), quad, @@ -268,6 +278,9 @@ namespace Functions bool left_over = true; // Current quadrature point + typename DH::active_cell_iterator cell = cell_hint.get(); + if (cell == dh->end()) + cell = dh->begin_active(); Point qp = mapping.transform_real_to_unit_cell(cell, points[0]); // Check if we already have a