DOLFINx
DOLFINx C++ interface
Expression.h
1// Copyright (C) 2020 - 2021 Jack S. Hale and Michal Habera
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "Function.h"
10#include <array>
11#include <dolfinx/common/utils.h>
12#include <dolfinx/mesh/Mesh.h>
13#include <functional>
14#include <span>
15#include <utility>
16#include <vector>
17
18namespace dolfinx::fem
19{
20template <typename T>
21class Constant;
22
31
32template <typename T>
34{
35 template <typename X, typename = void>
36 struct scalar_value_type
37 {
38 typedef X value_type;
39 };
40 template <typename X>
41 struct scalar_value_type<X, std::void_t<typename X::value_type>>
42 {
43 typedef typename X::value_type value_type;
44 };
45 using scalar_value_type_t = typename scalar_value_type<T>::value_type;
46
47public:
62 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
63 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
64 std::span<const double> X, std::array<std::size_t, 2> Xshape,
65 const std::function<void(T*, const T*, const T*,
66 const scalar_value_type_t*, const int*,
67 const uint8_t*)>
68 fn,
69 const std::vector<int>& value_shape,
70 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
71 const std::shared_ptr<const FunctionSpace> argument_function_space
72 = nullptr)
73 : _coefficients(coefficients), _constants(constants), _mesh(mesh),
74 _x_ref(std::vector<double>(X.begin(), X.end()), Xshape), _fn(fn),
75 _value_shape(value_shape),
76 _argument_function_space(argument_function_space)
77 {
78 // Extract mesh from argument's function space
79 if (!_mesh and argument_function_space)
80 _mesh = argument_function_space->mesh();
81 if (argument_function_space and _mesh != argument_function_space->mesh())
82 throw std::runtime_error("Incompatible mesh");
83 if (!_mesh)
84 {
85 throw std::runtime_error(
86 "No mesh could be associated with the Expression.");
87 }
88 }
89
91 Expression(Expression&& form) = default;
92
94 virtual ~Expression() = default;
95
98 std::shared_ptr<const FunctionSpace> argument_function_space() const
99 {
100 return _argument_function_space;
101 };
102
105 const std::vector<std::shared_ptr<const Function<T>>>& coefficients() const
106 {
107 return _coefficients;
108 }
109
114 const std::vector<std::shared_ptr<const Constant<T>>>& constants() const
115 {
116 return _constants;
117 }
118
122 std::vector<int> coefficient_offsets() const
123 {
124 std::vector<int> n{0};
125 for (auto& c : _coefficients)
126 {
127 if (!c)
128 throw std::runtime_error("Not all form coefficients have been set.");
129 n.push_back(n.back() + c->function_space()->element()->space_dimension());
130 }
131 return n;
132 }
133
139 template <typename U>
140 void eval(const std::span<const std::int32_t>& cells, U& values) const
141 {
142 // Extract data from Expression
143 assert(_mesh);
144
145 // Prepare coefficients and constants
146 const auto [coeffs, cstride] = pack_coefficients(*this, cells);
147 const std::vector<T> constant_data = pack_constants(*this);
148 const auto& fn = this->get_tabulate_expression();
149
150 // Prepare cell geometry
152 = _mesh->geometry().dofmap();
153 const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
154 std::span<const double> x_g = _mesh->geometry().x();
155
156 // Create data structures used in evaluation
157 std::vector<scalar_value_type_t> coordinate_dofs(3 * num_dofs_g);
158
159 int num_argument_dofs = 1;
160 std::span<const std::uint32_t> cell_info;
161 std::function<void(const std::span<T>&,
162 const std::span<const std::uint32_t>&, std::int32_t,
163 int)>
164 dof_transform_to_transpose
165 = [](const std::span<T>&, const std::span<const std::uint32_t>&,
166 std::int32_t, int)
167 {
168 // Do nothing
169 };
170
171 if (_argument_function_space)
172 {
173 num_argument_dofs
174 = _argument_function_space->dofmap()->element_dof_layout().num_dofs();
175 auto element = _argument_function_space->element();
176
177 assert(element);
178 if (element->needs_dof_transformations())
179 {
180 _mesh->topology_mutable().create_entity_permutations();
181 cell_info = std::span(_mesh->topology().get_cell_permutation_info());
182 dof_transform_to_transpose
183 = element
184 ->template get_dof_transformation_to_transpose_function<T>();
185 }
186 }
187
188 const int size0 = _x_ref.second[0] * value_size();
189 std::vector<T> values_local(size0 * num_argument_dofs, 0);
190 const std::span<T> _values_local(values_local);
191
192 // Iterate over cells and 'assemble' into values
193 for (std::size_t c = 0; c < cells.size(); ++c)
194 {
195 const std::int32_t cell = cells[c];
196
197 auto x_dofs = x_dofmap.links(cell);
198 for (std::size_t i = 0; i < x_dofs.size(); ++i)
199 {
200 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
201 std::next(coordinate_dofs.begin(), 3 * i));
202 }
203
204 const T* coeff_cell = coeffs.data() + c * cstride;
205 std::fill(values_local.begin(), values_local.end(), 0.0);
206 _fn(values_local.data(), coeff_cell, constant_data.data(),
207 coordinate_dofs.data(), nullptr, nullptr);
208
209 dof_transform_to_transpose(_values_local, cell_info, c, size0);
210
211 for (std::size_t j = 0; j < values_local.size(); ++j)
212 values(c, j) = values_local[j];
213 }
214 }
215
218 const std::function<void(T*, const T*, const T*, const scalar_value_type_t*,
219 const int*, const uint8_t*)>&
221 {
222 return _fn;
223 }
224
227 std::shared_ptr<const mesh::Mesh> mesh() const { return _mesh; }
228
231 int value_size() const
232 {
233 return std::reduce(_value_shape.begin(), _value_shape.end(), 1,
234 std::multiplies{});
235 }
236
239 const std::vector<int>& value_shape() const { return _value_shape; }
240
243 std::pair<std::vector<double>, std::array<std::size_t, 2>> X() const
244 {
245 return _x_ref;
246 }
247
249 using scalar_type = T;
250
251private:
252 // Function space for Argument
253 std::shared_ptr<const FunctionSpace> _argument_function_space;
254
255 // Coefficients associated with the Expression
256 std::vector<std::shared_ptr<const Function<T>>> _coefficients;
257
258 // Constants associated with the Expression
259 std::vector<std::shared_ptr<const Constant<T>>> _constants;
260
261 // Function to evaluate the Expression
262 std::function<void(T*, const T*, const T*, const scalar_value_type_t*,
263 const int*, const uint8_t*)>
264 _fn;
265
266 // The mesh
267 std::shared_ptr<const mesh::Mesh> _mesh;
268
269 // Shape of the evaluated expression
270 std::vector<int> _value_shape;
271
272 // Evaluation points on reference cell. Synonymous with X in public interface.
273 std::pair<std::vector<double>, std::array<std::size_t, 2>> _x_ref;
274};
275} // namespace dolfinx::fem
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:20
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:34
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Expression.h:122
int value_size() const
Get value size.
Definition: Expression.h:231
void eval(const std::span< const std::int32_t > &cells, U &values) const
Evaluate the expression on cells.
Definition: Expression.h:140
const std::vector< std::shared_ptr< const Function< T > > > & coefficients() const
Get coefficients.
Definition: Expression.h:105
const std::vector< std::shared_ptr< const Constant< T > > > & constants() const
Get constants.
Definition: Expression.h:114
Expression(Expression &&form)=default
Move constructor.
std::pair< std::vector< double >, std::array< std::size_t, 2 > > X() const
Evaluation points on the reference cell.
Definition: Expression.h:243
const std::vector< int > & value_shape() const
Get value shape.
Definition: Expression.h:239
virtual ~Expression()=default
Destructor.
std::shared_ptr< const FunctionSpace > argument_function_space() const
Get argument function space.
Definition: Expression.h:98
std::shared_ptr< const mesh::Mesh > mesh() const
Get mesh.
Definition: Expression.h:227
T scalar_type
Scalar type (T)
Definition: Expression.h:249
Expression(const std::vector< std::shared_ptr< const Function< T > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, std::span< const double > X, std::array< std::size_t, 2 > Xshape, const std::function< void(T *, const T *, const T *, const scalar_value_type_t *, const int *, const uint8_t *)> fn, const std::vector< int > &value_shape, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr, const std::shared_ptr< const FunctionSpace > argument_function_space=nullptr)
Create an Expression.
Definition: Expression.h:61
const std::function< void(T *, const T *, const T *, const scalar_value_type_t *, const int *, const uint8_t *)> & get_tabulate_expression() const
Get function for tabulate_expression.
Definition: Expression.h:220
This class represents a function in a finite element function space , given by.
Definition: Function.h:45
This class provides a static adjacency list data structure. It is commonly used to store directed gra...
Definition: AdjacencyList.h:26
std::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:111
Finite element method functionality.
Definition: assemble_matrix_impl.h:25
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:970
void pack_coefficients(const Form< T > &form, IntegralType integral_type, int id, const std::span< T > &c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition: utils.h:738