DOLFINx
DOLFINx C++ interface
Loading...
Searching...
No Matches
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 "Constant.h"
10#include "Function.h"
11#include <algorithm>
12#include <array>
13#include <concepts>
14#include <dolfinx/common/types.h>
15#include <dolfinx/mesh/Mesh.h>
16#include <functional>
17#include <span>
18#include <utility>
19#include <vector>
20
21namespace dolfinx::fem
22{
23template <dolfinx::scalar T>
24class Constant;
25
38template <dolfinx::scalar T,
39 std::floating_point U = dolfinx::scalar_value_type_t<T>>
41{
42public:
47 using scalar_type = T;
48
50 using geometry_type = U;
51
65 const std::vector<std::shared_ptr<
67 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
69 std::span<const geometry_type> X, std::array<std::size_t, 2> Xshape,
70 std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
71 const geometry_type*, const int*, const uint8_t*)>
72 fn,
73 const std::vector<int>& value_shape,
74 std::shared_ptr<const FunctionSpace<geometry_type>>
76 = nullptr)
77 : _coefficients(coefficients), _constants(constants),
78 _x_ref(std::vector<geometry_type>(X.begin(), X.end()), Xshape), _fn(fn),
79 _value_shape(value_shape),
80 _argument_function_space(argument_function_space)
81 {
82 for (auto& c : _coefficients)
83 {
84 assert(c);
85 if (c->function_space()->mesh()
86 != _coefficients.front()->function_space()->mesh())
87 {
88 throw std::runtime_error("Coefficients not all defined on same mesh.");
89 }
90 }
91 }
92
94 Expression(Expression&& e) = default;
95
97 virtual ~Expression() = default;
98
102 std::shared_ptr<const FunctionSpace<geometry_type>>
104 {
105 return _argument_function_space;
106 };
107
110 const std::vector<
111 std::shared_ptr<const Function<scalar_type, geometry_type>>>&
113 {
114 return _coefficients;
115 }
116
121 const std::vector<std::shared_ptr<const Constant<scalar_type>>>&
122 constants() const
123 {
124 return _constants;
125 }
126
132 std::vector<int> coefficient_offsets() const
133 {
134 std::vector<int> n{0};
135 for (auto& c : _coefficients)
136 {
137 if (!c)
138 throw std::runtime_error("Not all form coefficients have been set.");
139 n.push_back(n.back() + c->function_space()->element()->space_dimension());
140 }
141 return n;
142 }
143
154 std::span<const std::int32_t> entities,
155 std::span<scalar_type> values,
156 std::array<std::size_t, 2> vshape) const
157 {
158 std::size_t estride;
159 if (mesh.topology()->dim() == _x_ref.second[1])
160 {
161 estride = 1;
162 }
163 else if (mesh.topology()->dim() == _x_ref.second[1] + 1)
164 {
165 estride = 2;
166 }
167 else
168 {
169 throw std::runtime_error("Invalid dimension of evaluation points.");
170 }
171 // Prepare coefficients and constants
172 const auto [coeffs, cstride] = pack_coefficients(*this, entities, estride);
173 const std::vector<scalar_type> constant_data = pack_constants(*this);
174 auto fn = this->get_tabulate_expression();
175
176 // Prepare cell geometry
177 auto x_dofmap = mesh.geometry().dofmap();
178
179 // Get geometry data
180 auto& cmap = mesh.geometry().cmap();
181
182 const std::size_t num_dofs_g = cmap.dim();
183 auto x_g = mesh.geometry().x();
184
185 // Create data structures used in evaluation
186 std::vector<geometry_type> coord_dofs(3 * num_dofs_g);
187
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>,
191 std::int32_t, int)>
192 post_dof_transform
193 = [](std::span<scalar_type>, std::span<const std::uint32_t>,
194 std::int32_t, int)
195 {
196 // Do nothing
197 };
198
199 if (_argument_function_space)
200 {
201 num_argument_dofs
202 = _argument_function_space->dofmap()->element_dof_layout().num_dofs();
203 auto element = _argument_function_space->element();
204
205 assert(element);
206 if (element->needs_dof_transformations())
207 {
208 mesh.topology_mutable()->create_entity_permutations();
209 cell_info = std::span(mesh.topology()->get_cell_permutation_info());
210 post_dof_transform
211 = element
212 ->template get_post_dof_transformation_function<scalar_type>(
214 }
215 }
216
217 // Create get entity index function
218 std::function<const std::int32_t*(std::span<const std::int32_t>,
219 std::size_t)>
220 get_entity_index
221 = []([[maybe_unused]] std::span<const std::int32_t> entities,
222 [[maybe_unused]] std::size_t idx) { return nullptr; };
223 if (estride == 2)
224 {
225 get_entity_index
226 = [](std::span<const std::int32_t> entities, std::size_t idx)
227 { return entities.data() + 2 * idx + 1; };
228 }
229
230 // Iterate over cells and 'assemble' into values
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)
234 {
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)
239 {
240 std::copy_n(std::next(x_g.begin(), 3 * x_dofs[i]), 3,
241 std::next(coord_dofs.begin(), 3 * i));
242 }
243
244 const scalar_type* coeff_cell = coeffs.data() + e * cstride;
245 const int* entity_index = get_entity_index(entities, e);
246
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);
250
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];
254 }
255 }
256
259 const std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
260 const geometry_type*, const int*, const uint8_t*)>&
262 {
263 return _fn;
264 }
265
268 int value_size() const
269 {
270 return std::reduce(_value_shape.begin(), _value_shape.end(), 1,
271 std::multiplies{});
272 }
273
276 const std::vector<int>& value_shape() const { return _value_shape; }
277
280 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> X() const
281 {
282 return _x_ref;
283 }
284
285private:
286 // Function space for Argument
287 std::shared_ptr<const FunctionSpace<geometry_type>> _argument_function_space;
288
289 // Coefficients associated with the Expression
290 std::vector<std::shared_ptr<const Function<scalar_type, geometry_type>>>
291 _coefficients;
292
293 // Constants associated with the Expression
294 std::vector<std::shared_ptr<const Constant<scalar_type>>> _constants;
295
296 // Function to evaluate the Expression
297 std::function<void(scalar_type*, const scalar_type*, const scalar_type*,
298 const geometry_type*, const int*, const uint8_t*)>
299 _fn;
300
301 // Shape of the evaluated expression
302 std::vector<int> _value_shape;
303
304 // Evaluation points on reference cell. Synonymous with X in public
305 // interface.
306 std::pair<std::vector<geometry_type>, std::array<std::size_t, 2>> _x_ref;
307};
308} // namespace dolfinx::fem
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