#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
+#include <deal.II/base/thread_local_storage.h>
#include <deal.II/lac/vector.h>
* 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
std::vector<std::vector<unsigned int> > &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.
*/
const Mapping<dim> & 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
dh(&mydh, "FEFieldFunction"),
data_vector(myv),
mapping(mymapping),
+ cell_hint(dh->end()),
n_components(mydh.get_fe().n_components())
{
- cell = dh->begin_active();
}
FEFieldFunction<dim, DH, VECTOR>::
set_active_cell(typename DH::active_cell_iterator &newcell)
{
- cell = newcell;
+ cell_hint.get() = newcell;
}
-
+
template <int dim, typename DH, typename VECTOR>
void FEFieldFunction<dim, DH, VECTOR>::vector_value (const Point<dim> &p,
Vector<double> &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<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
// Check if we already have all we need
cell = my_pair.first;
qp = my_pair.second;
}
+
+ cell_hint.get() = cell;
// Now we can find out about the point
Quadrature<dim> quad(qp);
{
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<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
// Check if we already have all we need
cell = my_pair.first;
qp = my_pair.second;
}
-
+
+ cell_hint.get() = cell;
+
// Now we can find out about the point
Quadrature<dim> quad(qp);
FEValues<dim> fe_v(mapping, dh->get_fe(), quad,
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<dim> qp = mapping.transform_real_to_unit_cell(cell, points[0]);
// Check if we already have a