fill_inverse_JxW_values(
AlignedVector<VectorizedArray<Number>> &inverse_jxw) const
{
- constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
- Assert(inverse_jxw.size() > 0 && inverse_jxw.size() % dofs_per_cell == 0,
+ constexpr unsigned int dofs_per_component_on_cell =
+ Utilities::pow(fe_degree + 1, dim);
+ Assert(inverse_jxw.size() > 0 &&
+ inverse_jxw.size() % dofs_per_component_on_cell == 0,
ExcMessage(
"Expected diagonal to be a multiple of scalar dof per cells"));
// temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
// from fe_eval (will not reallocate any memory)
- const unsigned int previous_size = inverse_jxw.size();
- inverse_jxw.resize(dofs_per_cell);
- fe_eval.fill_JxW_values(inverse_jxw);
-
- // invert
- inverse_jxw.resize_fast(previous_size);
- for (unsigned int q = 0; q < dofs_per_cell; ++q)
- inverse_jxw[q] = 1. / inverse_jxw[q];
+ for (unsigned int q = 0; q < dofs_per_component_on_cell; ++q)
+ inverse_jxw[q] = 1. / fe_eval.JxW(q);
// copy values to rest of vector
- for (unsigned int q = dofs_per_cell; q < previous_size;)
- for (unsigned int i = 0; i < dofs_per_cell; ++i, ++q)
+ for (unsigned int q = dofs_per_component_on_cell; q < inverse_jxw.size();)
+ for (unsigned int i = 0; i < dofs_per_component_on_cell; ++i, ++q)
inverse_jxw[q] = inverse_jxw[i];
}
double error = 0, error2 = 0, abs = 0;
- QGauss<dim> quad(2);
- FEValues<dim> fe_values(fe, quad, update_JxW_values);
- FEEvaluation<dim, 1> fe_eval(mf_data);
- AlignedVector<VectorizedArray<double>> jxw_values_manual(fe_eval.n_q_points);
+ QGauss<dim> quad(2);
+ FEValues<dim> fe_values(fe, quad, update_JxW_values);
+ FEEvaluation<dim, 1> fe_eval(mf_data);
for (unsigned int cell = 0; cell < mf_data.n_macro_cells(); ++cell)
{
fe_eval.reinit(cell);
- fe_eval.fill_JxW_values(jxw_values_manual);
for (unsigned int v = 0; v < mf_data.n_components_filled(cell); ++v)
{
fe_values.reinit(mf_data.get_cell_iterator(cell, v));
{
abs += fe_values.JxW(q);
error += std::abs(fe_values.JxW(q) - fe_eval.JxW(q)[v]);
- error2 += std::abs(fe_values.JxW(q) - jxw_values_manual[q][v]);
}
}
}
- deallog << "Norm of difference: " << error / abs << " " << error2 / abs
- << std::endl
- << std::endl;
+ deallog << "Norm of difference: " << error / abs << std::endl << std::endl;
}