+ template <int dim, typename Number, typename Number2>
+ SymmetricTensor<2, dim, typename ProductTypeNoPoint<Number, Number2>::type>
+ evaluate_tensor_product_hessian(
+ const std::vector<Polynomials::Polynomial<double>> &poly,
+ const std::vector<Number> & values,
+ const Point<dim, Number2> & p,
+ const std::vector<unsigned int> & renumber = {})
+ {
+ static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented");
+
+ using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
+
+ // use `int` type for this variable and the loops below to inform the
+ // compiler that the loops below will never overflow, which allows it to
+ // generate more optimized code for the variable loop bounds in the
+ // present context
+ const int n_shapes = poly.size();
+ AssertDimension(Utilities::pow(n_shapes, dim), values.size());
+ Assert(renumber.empty() || renumber.size() == values.size(),
+ ExcDimensionMismatch(renumber.size(), values.size()));
+
+ AssertIndexRange(n_shapes, 200);
+ std::array<Number2, 3 * dim * 200> shapes;
+
+ // Evaluate 1D polynomials and their derivatives
+ for (unsigned int d = 0; d < dim; ++d)
+ for (int i = 0; i < n_shapes; ++i)
+ poly[i].value(p[d], 2, shapes.data() + 3 * (d * n_shapes + i));
+
+ // Go through the tensor product of shape functions and interpolate
+ // with optimal algorithm
+ SymmetricTensor<2, dim, Number3> result;
+ for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2)
+ {
+ Number3 value_y = {}, deriv_x = {}, deriv_y = {}, deriv_xx = {},
+ deriv_xy = {}, deriv_yy = {};
+ for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1)
+ {
+ // Interpolation + derivative x direction
+ Number3 value = {}, deriv_1 = {}, deriv_2 = {};
+
+ // Distinguish the inner loop based on whether we have a
+ // renumbering or not
+ if (renumber.empty())
+ for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+ {
+ value += shapes[3 * i0] * values[i];
+ deriv_1 += shapes[3 * i0 + 1] * values[i];
+ deriv_2 += shapes[3 * i0 + 2] * values[i];
+ }
+ else
+ for (int i0 = 0; i0 < n_shapes; ++i0, ++i)
+ {
+ value += shapes[3 * i0] * values[renumber[i]];
+ deriv_1 += shapes[3 * i0 + 1] * values[renumber[i]];
+ deriv_2 += shapes[3 * i0 + 2] * values[renumber[i]];
+ }
+
+ // Interpolation + derivative in y direction
+ if (dim > 1)
+ {
+ if (dim > 2)
+ {
+ value_y += value * shapes[3 * n_shapes + 3 * i1];
+ deriv_x += deriv_1 * shapes[3 * n_shapes + 3 * i1];
+ deriv_y += value * shapes[3 * n_shapes + 3 * i1 + 1];
+ }
+ deriv_xx += deriv_2 * shapes[3 * n_shapes + 3 * i1];
+ deriv_xy += deriv_1 * shapes[3 * n_shapes + 3 * i1 + 1];
+ deriv_yy += value * shapes[3 * n_shapes + 3 * i1 + 2];
+ }
+ else
+ {
+ result[0][0] = deriv_2;
+ }
+ }
+ if (dim == 3)
+ {
+ // Interpolation + derivative in z direction
+ result[0][0] += deriv_xx * shapes[6 * n_shapes + 3 * i2];
+ result[0][1] += deriv_xy * shapes[6 * n_shapes + 3 * i2];
+ result[0][2] += deriv_x * shapes[6 * n_shapes + 3 * i2 + 1];
+ result[1][1] += deriv_yy * shapes[6 * n_shapes + 3 * i2];
+ result[1][2] += deriv_y * shapes[6 * n_shapes + 3 * i2 + 1];
+ result[2][2] += value_y * shapes[6 * n_shapes + 3 * i2 + 2];
+ }
+ else if (dim == 2)
+ {
+ result[0][0] = deriv_xx;
+ result[1][0] = deriv_xy;
+ result[1][1] = deriv_yy;
+ }
+ }
+
+ return result;
+ }
+
+
+
/**
* Same as evaluate_tensor_product_value_and_gradient() but for integration.
*/
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
#include <deal.II/non_matching/fe_values.h>
DEAL_II_NAMESPACE_OPEN
* set_active_cell()
*/
std::vector<typename VectorType::value_type> local_dof_values;
+
+ /**
+ * Description of the 1D polynomial basis for tensor product elements
+ * used for the fast path of this class using tensor product
+ * evaluators.
+ */
+ std::vector<Polynomials::Polynomial<double>> poly;
+
+ /**
+ * Renumbering for the tensor-product evaluator in the fast path.
+ */
+ std::vector<unsigned int> renumber;
+
+ /**
+ * Check whether the shape functions are linear.
+ */
+ bool polynomials_are_hat_functions;
};
// Save the element and the local dof values, since this is what we need
// to evaluate the function.
- element = &dof_handler_cell->get_fe();
+
+ // Check if we can use the fast path. In case we have a different
+ // element from the one used before we need to set up the data
+ // structures again.
+ if (element != &dof_handler_cell->get_fe())
+ {
+ poly.clear();
+ element = &dof_handler_cell->get_fe();
+
+ if (element->n_base_elements() == 1 &&
+ dealii::internal::FEPointEvaluation::is_fast_path_supported(
+ *element, 0))
+ {
+ dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
+ shape_info;
+
+ shape_info.reinit(QMidpoint<1>(), *element, 0);
+ renumber = shape_info.lexicographic_numbering;
+ poly =
+ dealii::internal::FEPointEvaluation::get_polynomial_space(
+ element->base_element(0));
+
+ polynomials_are_hat_functions =
+ (poly.size() == 2 && poly[0].value(0.) == 1. &&
+ poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
+ poly[1].value(1.) == 1.);
+ }
+ }
+ else
+ element = &dof_handler_cell->get_fe();
local_dof_indices.resize(element->dofs_per_cell);
dof_handler_cell->get_dof_indices(local_dof_indices);
AssertIndexRange(component, this->n_components);
Assert(cell_is_set(), ExcCellNotSet());
- double value = 0;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- value += local_dof_values[i] *
- element->shape_value_component(i, point, component);
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_value_and_gradient(
+ poly,
+ local_dof_values,
+ point,
+ polynomials_are_hat_functions,
+ renumber)
+ .first;
+ }
+ else
+ {
+ double value = 0;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ value += local_dof_values[i] *
+ element->shape_value_component(i, point, component);
- return value;
+ return value;
+ }
}
AssertIndexRange(component, this->n_components);
Assert(cell_is_set(), ExcCellNotSet());
- Tensor<1, dim> gradient;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- gradient += local_dof_values[i] *
- element->shape_grad_component(i, point, component);
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_value_and_gradient(
+ poly,
+ local_dof_values,
+ point,
+ polynomials_are_hat_functions,
+ renumber)
+ .second;
+ }
+ else
+ {
+ Tensor<1, dim> gradient;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ gradient += local_dof_values[i] *
+ element->shape_grad_component(i, point, component);
- return gradient;
+ return gradient;
+ }
}
AssertIndexRange(component, this->n_components);
Assert(cell_is_set(), ExcCellNotSet());
- Tensor<2, dim> hessian;
- for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
- hessian += local_dof_values[i] *
- element->shape_grad_grad_component(i, point, component);
+ if (!poly.empty() && component == 0)
+ {
+ // TODO: this could be extended to a component that is not zero
+ return dealii::internal::evaluate_tensor_product_hessian(
+ poly, local_dof_values, point, renumber);
+ }
+ else
+ {
+ Tensor<2, dim> hessian;
+ for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
+ hessian +=
+ local_dof_values[i] *
+ element->shape_grad_grad_component(i, point, component);
- return symmetrize(hessian);
+ return symmetrize(hessian);
+ }
}
} // namespace FEValuesImplementation
} // namespace internal
Assert(fe_collection->size() > 0,
ExcMessage("Incoming hp::FECollection can not be empty."));
- Assert(
- mapping_collection->size() == fe_collection->size() ||
- mapping_collection->size() == 1,
- ExcMessage(
- "Size of hp::MappingCollection must be the same as hp::FECollection or 1."));
- Assert(
- q_collection.size() == fe_collection->size() || q_collection.size() == 1,
- ExcMessage(
- "Size of hp::QCollection<dim> must be the same as hp::FECollection or 1."));
- Assert(
- q_collection_1D.size() == fe_collection->size() ||
- q_collection_1D.size() == 1,
- ExcMessage(
- "Size of hp::QCollection<1> must be the same as hp::FECollection or 1."));
+ Assert(mapping_collection->size() == fe_collection->size() ||
+ mapping_collection->size() == 1,
+ ExcMessage("Size of hp::MappingCollection must be "
+ "the same as hp::FECollection or 1."));
+ Assert(q_collection.size() == fe_collection->size() ||
+ q_collection.size() == 1,
+ ExcMessage("Size of hp::QCollection<dim> must be the "
+ "same as hp::FECollection or 1."));
+ Assert(q_collection_1D.size() == fe_collection->size() ||
+ q_collection_1D.size() == 1,
+ ExcMessage("Size of hp::QCollection<1> must be the "
+ "same as hp::FECollection or 1."));
// For each element in fe_collection, create dealii::FEValues objects to use
// on the non-intersected cells.
quadrature_generator.get_surface_quadrature();
// Even if a cell is formally intersected the number of created
- // quadrature points can be 0. Avoid creating an FEValues object if
- // that is the case.
+ // quadrature points can be 0. Avoid creating an FEValues object
+ // if that is the case.
if (inside_quadrature.size() > 0)
{
fe_values_inside.emplace((*mapping_collection)[mapping_index],
projects / dealii.git / commitdiff