#include <deal.II/base/tensor_polynomials_base.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <mutex>
#include <vector>
+
DEAL_II_NAMESPACE_OPEN
/**
*/
PolynomialsRaviartThomas(const unsigned int k);
+ /**
+ * Copy constructor.
+ */
+ PolynomialsRaviartThomas(const PolynomialsRaviartThomas &other);
+
/**
* Compute the value and the first and second derivatives of each Raviart-
* Thomas polynomial at @p unit_point.
*/
static std::vector<std::vector<Polynomials::Polynomial<double>>>
create_polynomials(const unsigned int k);
+
+ /**
+ * A mutex with which to guard access to the following `mutable`
+ * variables.
+ */
+ mutable std::mutex scratch_arrays_mutex;
+
+ /**
+ * The following arrays are used as scratch data in the evaluate() function.
+ * Since that function is `const`, they have to be marked as `mutable`, and
+ * since we want that function to be thread-safe, access to these scratch
+ * arrays is guarded by the mutex above.
+ */
+ mutable std::vector<double> scratch_values;
+ mutable std::vector<Tensor<1, dim>> scratch_grads;
+ mutable std::vector<Tensor<2, dim>> scratch_grad_grads;
+ mutable std::vector<Tensor<3, dim>> scratch_third_derivatives;
+ mutable std::vector<Tensor<4, dim>> scratch_fourth_derivatives;
};
+
template <int dim>
inline std::string
PolynomialsRaviartThomas<dim>::name() const
+template <int dim>
+PolynomialsRaviartThomas<dim>::PolynomialsRaviartThomas(
+ const PolynomialsRaviartThomas &other)
+ : TensorPolynomialsBase<dim>(other)
+ , polynomial_space(other.polynomial_space)
+// no need to copy the scratch data, or the mutex
+{}
+
+
+
template <int dim>
std::vector<std::vector<Polynomials::Polynomial<double>>>
PolynomialsRaviartThomas<dim>::create_polynomials(const unsigned int k)
}
+
template <int dim>
void
PolynomialsRaviartThomas<dim>::evaluate(
fourth_derivatives.size() == 0,
ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
- // have a few scratch
- // arrays. because we don't want to
- // re-allocate them every time this
- // function is called, we make them
- // static. however, in return we
- // have to ensure that the calls to
- // the use of these variables is
- // locked with a mutex. if the
- // mutex is removed, several tests
- // (notably
- // deal.II/create_mass_matrix_05)
- // will start to produce random
- // results in multithread mode
- static std::mutex mutex;
- std::lock_guard<std::mutex> lock(mutex);
-
- static std::vector<double> p_values;
- static std::vector<Tensor<1, dim>> p_grads;
- static std::vector<Tensor<2, dim>> p_grad_grads;
- static std::vector<Tensor<3, dim>> p_third_derivatives;
- static std::vector<Tensor<4, dim>> p_fourth_derivatives;
+ // In the following, we will access the scratch arrays declared as 'mutable'
+ // in the class. In order to do so from this function, we have to make sure
+ // that we guard this access by a mutex so that the invokation of this 'const'
+ // function is thread-safe.
+ std::lock_guard<std::mutex> lock(scratch_arrays_mutex);
const unsigned int n_sub = polynomial_space.n();
- p_values.resize((values.size() == 0) ? 0 : n_sub);
- p_grads.resize((grads.size() == 0) ? 0 : n_sub);
- p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
- p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
- p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
+ scratch_values.resize((values.size() == 0) ? 0 : n_sub);
+ scratch_grads.resize((grads.size() == 0) ? 0 : n_sub);
+ scratch_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
+ scratch_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
+ scratch_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 :
+ n_sub);
for (unsigned int d = 0; d < dim; ++d)
{
p(c) = unit_point((c + d) % dim);
polynomial_space.evaluate(p,
- p_values,
- p_grads,
- p_grad_grads,
- p_third_derivatives,
- p_fourth_derivatives);
+ scratch_values,
+ scratch_grads,
+ scratch_grad_grads,
+ scratch_third_derivatives,
+ scratch_fourth_derivatives);
- for (unsigned int i = 0; i < p_values.size(); ++i)
- values[i + d * n_sub][d] = p_values[i];
+ for (unsigned int i = 0; i < scratch_values.size(); ++i)
+ values[i + d * n_sub][d] = scratch_values[i];
- for (unsigned int i = 0; i < p_grads.size(); ++i)
+ for (unsigned int i = 0; i < scratch_grads.size(); ++i)
for (unsigned int d1 = 0; d1 < dim; ++d1)
- grads[i + d * n_sub][d][(d1 + d) % dim] = p_grads[i][d1];
+ grads[i + d * n_sub][d][(d1 + d) % dim] = scratch_grads[i][d1];
- for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
+ for (unsigned int i = 0; i < scratch_grad_grads.size(); ++i)
for (unsigned int d1 = 0; d1 < dim; ++d1)
for (unsigned int d2 = 0; d2 < dim; ++d2)
grad_grads[i + d * n_sub][d][(d1 + d) % dim][(d2 + d) % dim] =
- p_grad_grads[i][d1][d2];
+ scratch_grad_grads[i][d1][d2];
- for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
+ for (unsigned int i = 0; i < scratch_third_derivatives.size(); ++i)
for (unsigned int d1 = 0; d1 < dim; ++d1)
for (unsigned int d2 = 0; d2 < dim; ++d2)
for (unsigned int d3 = 0; d3 < dim; ++d3)
third_derivatives[i + d * n_sub][d][(d1 + d) % dim]
[(d2 + d) % dim][(d3 + d) % dim] =
- p_third_derivatives[i][d1][d2][d3];
+ scratch_third_derivatives[i][d1][d2][d3];
- for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
+ for (unsigned int i = 0; i < scratch_fourth_derivatives.size(); ++i)
for (unsigned int d1 = 0; d1 < dim; ++d1)
for (unsigned int d2 = 0; d2 < dim; ++d2)
for (unsigned int d3 = 0; d3 < dim; ++d3)
fourth_derivatives[i + d * n_sub][d][(d1 + d) % dim]
[(d2 + d) % dim][(d3 + d) % dim]
[(d4 + d) % dim] =
- p_fourth_derivatives[i][d1][d2][d3][d4];
+ scratch_fourth_derivatives[i][d1][d2][d3]
+ [d4];
}
}