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 <dolfinx/common/utils.h>
11#include <dolfinx/mesh/Mesh.h>
12#include <functional>
13#include <utility>
14#include <vector>
15#include <xtensor/xtensor.hpp>
16#include <xtl/xspan.hpp>
17
18namespace dolfinx::fem
19{
20template <typename T>
21class Constant;
22
31
32template <typename T>
34{
35public:
49 const std::vector<std::shared_ptr<const Function<T>>>& coefficients,
50 const std::vector<std::shared_ptr<const Constant<T>>>& constants,
51 const xt::xtensor<double, 2>& X,
52 const std::function<void(T*, const T*, const T*, const double*,
53 const int*, const uint8_t*)>
54 fn,
55 const std::vector<int>& value_shape,
56 const std::shared_ptr<const mesh::Mesh>& mesh = nullptr,
57 const std::shared_ptr<const FunctionSpace> argument_function_space
58 = nullptr)
59 : _coefficients(coefficients), _constants(constants), _mesh(mesh),
60 _x_ref(X), _fn(fn), _value_shape(value_shape),
61 _argument_function_space(argument_function_space)
62 {
63 // Extract mesh from argument's function space
64 if (!_mesh and argument_function_space)
65 _mesh = argument_function_space->mesh();
66 if (argument_function_space and _mesh != argument_function_space->mesh())
67 throw std::runtime_error("Incompatible mesh");
68 if (!_mesh)
69 throw std::runtime_error(
70 "No mesh could be associated with the Expression.");
71 }
72
74 Expression(Expression&& form) = default;
75
77 virtual ~Expression() = default;
78
81 std::shared_ptr<const fem::FunctionSpace> argument_function_space() const
82 {
83 return _argument_function_space;
84 };
85
88 const std::vector<std::shared_ptr<const fem::Function<T>>>&
90 {
91 return _coefficients;
92 }
93
98 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants() const
99 {
100 return _constants;
101 }
102
106 std::vector<int> coefficient_offsets() const
107 {
108 std::vector<int> n{0};
109 for (auto& c : _coefficients)
110 {
111 if (!c)
112 throw std::runtime_error("Not all form coefficients have been set.");
113 n.push_back(n.back() + c->function_space()->element()->space_dimension());
114 }
115 return n;
116 }
117
123 template <typename U>
124 void eval(const xtl::span<const std::int32_t>& cells, U& values) const
125 {
126 // Extract data from Expression
127 assert(_mesh);
128
129 // Prepare coefficients and constants
130 const auto [coeffs, cstride] = pack_coefficients(*this, cells);
131 const std::vector<T> constant_data = pack_constants(*this);
132 const auto& fn = this->get_tabulate_expression();
133
134 // Prepare cell geometry
136 = _mesh->geometry().dofmap();
137 const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
138 xtl::span<const double> x_g = _mesh->geometry().x();
139
140 // Create data structures used in evaluation
141 std::vector<double> coordinate_dofs(3 * num_dofs_g);
142
143 int num_argument_dofs = 1;
144 xtl::span<const std::uint32_t> cell_info;
145 std::function<void(const xtl::span<T>&,
146 const xtl::span<const std::uint32_t>&, std::int32_t,
147 int)>
148 dof_transform_to_transpose
149 = [](const xtl::span<T>&, const xtl::span<const std::uint32_t>&,
150 std::int32_t, int)
151 {
152 // Do nothing
153 };
154
155 if (_argument_function_space)
156 {
157 num_argument_dofs
158 = _argument_function_space->dofmap()->element_dof_layout().num_dofs();
159 auto element = _argument_function_space->element();
160
161 assert(element);
162 if (element->needs_dof_transformations())
163 {
164 _mesh->topology_mutable().create_entity_permutations();
165 cell_info = xtl::span(_mesh->topology().get_cell_permutation_info());
166 dof_transform_to_transpose
167 = element
168 ->template get_dof_transformation_to_transpose_function<T>();
169 }
170 }
171
172 const int size0 = _x_ref.shape(0) * value_size();
173 std::vector<T> values_local(size0 * num_argument_dofs, 0);
174 const xtl::span<T> _values_local(values_local);
175
176 // Iterate over cells and 'assemble' into values
177 for (std::size_t c = 0; c < cells.size(); ++c)
178 {
179 const std::int32_t cell = cells[c];
180
181 auto x_dofs = x_dofmap.links(cell);
182 for (std::size_t i = 0; i < x_dofs.size(); ++i)
183 {
184 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
185 std::next(coordinate_dofs.begin(), 3 * i));
186 }
187
188 const T* coeff_cell = coeffs.data() + c * cstride;
189 std::fill(values_local.begin(), values_local.end(), 0.0);
190 _fn(values_local.data(), coeff_cell, constant_data.data(),
191 coordinate_dofs.data(), nullptr, nullptr);
192
193 dof_transform_to_transpose(_values_local, cell_info, c, size0);
194
195 for (std::size_t j = 0; j < values_local.size(); ++j)
196 values(c, j) = values_local[j];
197 }
198 }
199
202 const std::function<void(T*, const T*, const T*, const double*, const int*,
203 const uint8_t*)>&
205 {
206 return _fn;
207 }
208
211 std::shared_ptr<const mesh::Mesh> mesh() const { return _mesh; }
212
215 int value_size() const
216 {
217 return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
218 std::multiplies<int>());
219 }
220
223 const std::vector<int>& value_shape() const { return _value_shape; }
224
227 const xt::xtensor<double, 2>& X() const { return _x_ref; }
228
230 using scalar_type = T;
231
232private:
233 // Function space for Argument
234 std::shared_ptr<const fem::FunctionSpace> _argument_function_space;
235
236 // Coefficients associated with the Expression
237 std::vector<std::shared_ptr<const fem::Function<T>>> _coefficients;
238
239 // Constants associated with the Expression
240 std::vector<std::shared_ptr<const fem::Constant<T>>> _constants;
241
242 // Function to evaluate the Expression
243 std::function<void(T*, const T*, const T*, const double*, const int*,
244 const uint8_t*)>
245 _fn;
246
247 // The mesh
248 std::shared_ptr<const mesh::Mesh> _mesh;
249
250 // Shape of the evaluated expression
251 std::vector<int> _value_shape;
252
253 // Evaluation points on reference cell. Synonymous with X in public interface.
254 xt::xtensor<double, 2> _x_ref;
255};
256} // namespace dolfinx::fem
Constant value which can be attached to a Form. Constants may be scalar (rank 0), vector (rank 1),...
Definition: Constant.h:21
Represents a mathematical expression evaluated at a pre-defined set of points on the reference cell....
Definition: Expression.h:34
void eval(const xtl::span< const std::int32_t > &cells, U &values) const
Evaluate the expression on cells.
Definition: Expression.h:124
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:106
const std::vector< std::shared_ptr< const fem::Constant< T > > > & constants() const
Get constants.
Definition: Expression.h:98
int value_size() const
Get value size.
Definition: Expression.h:215
Expression(const std::vector< std::shared_ptr< const Function< T > > > &coefficients, const std::vector< std::shared_ptr< const Constant< T > > > &constants, const xt::xtensor< double, 2 > &X, const std::function< void(T *, const T *, const T *, const double *, 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:48
Expression(Expression &&form)=default
Move constructor.
const std::vector< int > & value_shape() const
Get value shape.
Definition: Expression.h:223
const xt::xtensor< double, 2 > & X() const
Get evaluation points on reference cell.
Definition: Expression.h:227
virtual ~Expression()=default
Destructor.
std::shared_ptr< const mesh::Mesh > mesh() const
Get mesh.
Definition: Expression.h:211
const std::function< void(T *, const T *, const T *, const double *, const int *, const uint8_t *)> & get_tabulate_expression() const
Get function for tabulate_expression.
Definition: Expression.h:204
std::shared_ptr< const fem::FunctionSpace > argument_function_space() const
Get argument function space.
Definition: Expression.h:81
const std::vector< std::shared_ptr< const fem::Function< T > > > & coefficients() const
Get coefficients.
Definition: Expression.h:89
T scalar_type
Scalar type (T)
Definition: Expression.h:230
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:46
xtl::span< T > links(int node)
Get the links (edges) for given node.
Definition: AdjacencyList.h:131
Finite element method functionality.
Definition: assemble_matrix_impl.h:24
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:891
void pack_coefficients(const Form< T > &form, IntegralType integral_type, int id, const xtl::span< T > &c, int cstride)
Pack coefficients of a Form for a given integral type and domain id.
Definition: utils.h:656