14#include <dolfinx/common/types.h>
15#include <dolfinx/mesh/Mesh.h>
23template <dolfinx::scalar T>
39 std::floating_point U = dolfinx::scalar_value_type_t<T>>
65 const std::vector<std::shared_ptr<
69 std::span<const geometry_type>
X, std::array<std::size_t, 2> Xshape,
78 _x_ref(std::vector<
geometry_type>(
X.begin(),
X.end()), Xshape), _fn(fn),
82 for (
auto& c : _coefficients)
85 if (c->function_space()->mesh()
86 != _coefficients.front()->function_space()->mesh())
88 throw std::runtime_error(
"Coefficients not all defined on same mesh.");
102 std::shared_ptr<const FunctionSpace<geometry_type>>
105 return _argument_function_space;
111 std::shared_ptr<const Function<scalar_type, geometry_type>>>&
114 return _coefficients;
121 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
134 std::vector<int> n{0};
135 for (
auto& c : _coefficients)
138 throw std::runtime_error(
"Not all form coefficients have been set.");
139 n.push_back(n.back() + c->function_space()->element()->space_dimension());
154 std::span<const std::int32_t> entities,
155 std::span<scalar_type> values,
156 std::array<std::size_t, 2> vshape)
const
159 if (mesh.topology()->dim() == _x_ref.second[1])
163 else if (mesh.topology()->dim() == _x_ref.second[1] + 1)
169 throw std::runtime_error(
"Invalid dimension of evaluation points.");
173 const std::vector<scalar_type> constant_data =
pack_constants(*
this);
177 auto x_dofmap = mesh.geometry().dofmap();
180 auto& cmap = mesh.geometry().cmap();
182 const std::size_t num_dofs_g = cmap.dim();
183 auto x_g = mesh.geometry().x();
186 std::vector<geometry_type> coord_dofs(3 * num_dofs_g);
188 int num_argument_dofs = 1;
189 std::span<const std::uint32_t> cell_info;
190 std::function<void(std::span<scalar_type>, std::span<const std::uint32_t>,
193 = [](std::span<scalar_type>, std::span<const std::uint32_t>,
199 if (_argument_function_space)
202 = _argument_function_space->dofmap()->element_dof_layout().num_dofs();
203 auto element = _argument_function_space->element();
206 if (element->needs_dof_transformations())
208 mesh.topology_mutable()->create_entity_permutations();
209 cell_info = std::span(mesh.topology()->get_cell_permutation_info());
212 ->template get_post_dof_transformation_function<scalar_type>(
218 std::function<
const std::int32_t*(std::span<const std::int32_t>,
221 = []([[maybe_unused]] std::span<const std::int32_t> entities,
222 [[maybe_unused]] std::size_t idx) {
return nullptr; };
226 = [](std::span<const std::int32_t> entities, std::size_t idx)
227 {
return entities.data() + 2 * idx + 1; };
231 const int size0 = _x_ref.second[0] *
value_size();
232 std::vector<scalar_type> values_local(size0 * num_argument_dofs, 0);
233 for (std::size_t e = 0; e < entities.size() / estride; ++e)
235 const std::int32_t entity = entities[e * estride];
236 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
237 x_dofmap, entity, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
238 for (std::size_t i = 0; i < x_dofs.size(); ++i)
240 std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3,
241 std::next(coord_dofs.begin(), 3 * i));
244 const scalar_type* coeff_cell = coeffs.data() + e * cstride;
245 const int* entity_index = get_entity_index(entities, e);
247 std::fill(values_local.begin(), values_local.end(), 0);
248 _fn(values_local.data(), coeff_cell, constant_data.data(),
249 coord_dofs.data(), entity_index,
nullptr);
251 post_dof_transform(values_local, cell_info, e, size0);
252 for (std::size_t j = 0; j < values_local.size(); ++j)
253 values[e * vshape[1] + j] = values_local[j];
270 return std::reduce(_value_shape.begin(), _value_shape.end(), 1,
276 const std::vector<int>&
value_shape()
const {
return _value_shape; }
280 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>>
X()
const
287 std::shared_ptr<const FunctionSpace<geometry_type>> _argument_function_space;
290 std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>>
294 std::vector<std::shared_ptr<const Constant<scalar_type>>> _constants;
302 std::vector<int> _value_shape;
306 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _x_ref;
Constant value which can be attached to a Form.
Definition Constant.h:23
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell.
Definition Expression.h:41
U geometry_type
Geometry type of the points.
Definition Expression.h:50
std::pair< std::vector< geometry_type >, std::array< std::size_t, 2 > > X() const
Evaluation points on the reference cell.
Definition Expression.h:280
const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > & coefficients() const
Get coefficients.
Definition Expression.h:112
Expression(const std::vector< std::shared_ptr< const Function< scalar_type, geometry_type > > > &coefficients, const std::vector< std::shared_ptr< const Constant< scalar_type > > > &constants, std::span< const geometry_type > X, std::array< std::size_t, 2 > Xshape, std::function< void(scalar_type *, const scalar_type *, const scalar_type *, const geometry_type *, const int *, const uint8_t *)> fn, const std::vector< int > &value_shape, std::shared_ptr< const FunctionSpace< geometry_type > > argument_function_space=nullptr)
Create an Expression.
Definition Expression.h:64
void eval(const mesh::Mesh< geometry_type > &mesh, std::span< const std::int32_t > entities, std::span< scalar_type > values, std::array< std::size_t, 2 > vshape) const
Evaluate Expression on cells or facets.
Definition Expression.h:153
const std::vector< std::shared_ptr< const Constant< scalar_type > > > & constants() const
Get constants.
Definition Expression.h:122
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell.
Definition Expression.h:132
int value_size() const
Get value size.
Definition Expression.h:268
Expression(Expression &&e)=default
Move constructor.
const std::vector< int > & value_shape() const
Get value shape.
Definition Expression.h:276
virtual ~Expression()=default
Destructor.
T scalar_type
Scalar type.
Definition Expression.h:47
const std::function< void(scalar_type *, const scalar_type *, const scalar_type *, const geometry_type *, const int *, const uint8_t *)> & get_tabulate_expression() const
Get function for tabulate_expression.
Definition Expression.h:261
std::shared_ptr< const FunctionSpace< geometry_type > > argument_function_space() const
Get argument function space.
Definition Expression.h:103
Finite Element, containing the dof layout on a reference element, and various methods for evaluating ...
Definition FiniteElement.h:29
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
This class represents a function in a finite element function space , given by.
Definition Function.h:45
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
This concept is used to constrain the a template type to floating point real or complex types....
Definition types.h:20
Finite element method functionality.
Definition assemble_matrix_impl.h:26
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u into a single array ready for assembly.
Definition utils.h:1249
void pack_coefficients(const Form< T, U > &form, IntegralType integral_type, int id, std::span< T > c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition utils.h:1019