DOLFINx
DOLFINx C++ interface
discreteoperators.h
1// Copyright (C) 2015-2022 Garth N. Wells, Jørgen S. Dokken
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 "DofMap.h"
10#include "FiniteElement.h"
11#include "FunctionSpace.h"
12#include <array>
13#include <dolfinx/common/IndexMap.h>
14#include <dolfinx/common/math.h>
15#include <dolfinx/common/utils.h>
16#include <dolfinx/mesh/Mesh.h>
17#include <memory>
18#include <span>
19#include <vector>
20
21namespace dolfinx::fem
22{
23
47template <typename T, typename U>
49 U&& mat_set)
50{
51 namespace stdex = std::experimental;
52 using cmdspan2_t
53 = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
54 using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
55 using cmdspan4_t
56 = stdex::mdspan<const double, stdex::dextents<std::size_t, 4>>;
57
58 // Get mesh
59 std::shared_ptr<const mesh::Mesh> mesh = V1.mesh();
60 assert(mesh);
61
62 // Check spaces
63 std::shared_ptr<const FiniteElement> e0 = V0.element();
64 assert(e0);
65 if (e0->map_type() != basix::maps::type::identity)
66 throw std::runtime_error("Wrong finite element space for V0.");
67 if (e0->block_size() != 1)
68 throw std::runtime_error("Block size is greather than 1 for V0.");
69 if (e0->reference_value_size() != 1)
70 throw std::runtime_error("Wrong value size for V0.");
71
72 std::shared_ptr<const FiniteElement> e1 = V1.element();
73 assert(e1);
74 if (e1->map_type() != basix::maps::type::covariantPiola)
75 throw std::runtime_error("Wrong finite element space for V1.");
76 if (e1->block_size() != 1)
77 throw std::runtime_error("Block size is greather than 1 for V1.");
78
79 // Get V0 (H(curl)) space interpolation points
80 const auto [X, Xshape] = e1->interpolation_points();
81
82 // Tabulate first order derivatives of Lagrange space at H(curl)
83 // interpolation points
84 const int ndofs0 = e0->space_dimension();
85 const int tdim = mesh->topology().dim();
86 std::vector<double> phi0_b((tdim + 1) * Xshape[0] * ndofs0 * 1);
87 cmdspan4_t phi0(phi0_b.data(), tdim + 1, Xshape[0], ndofs0, 1);
88 e0->tabulate(phi0_b, X, Xshape, 1);
89
90 // Reshape lagrange basis derivatives as a matrix of shape (tdim *
91 // num_points, num_dofs_per_cell)
92 cmdspan2_t dphi_reshaped(
93 phi0_b.data() + phi0.extent(3) * phi0.extent(2) * phi0.extent(1),
94 tdim * phi0.extent(1), phi0.extent(2));
95
96 // Get inverse DOF transform function
97 auto apply_inverse_dof_transform
98 = e1->get_dof_transformation_function<T>(true, true, false);
99
100 // Generate cell permutations
101 mesh->topology_mutable().create_entity_permutations();
102 const std::vector<std::uint32_t>& cell_info
103 = mesh->topology().get_cell_permutation_info();
104
105 // Create element kernel function
106 std::shared_ptr<const DofMap> dofmap0 = V0.dofmap();
107 assert(dofmap0);
108 std::shared_ptr<const DofMap> dofmap1 = V1.dofmap();
109 assert(dofmap1);
110
111 // Build the element interpolation matrix
112 std::vector<T> Ab(e1->space_dimension() * ndofs0);
113 {
114 stdex::mdspan<T, stdex::dextents<std::size_t, 2>> A(
115 Ab.data(), e1->space_dimension(), ndofs0);
116 const auto [Pi, shape] = e1->interpolation_operator();
117 cmdspan2_t _Pi(Pi.data(), shape);
118 math::dot(_Pi, dphi_reshaped, A);
119 }
120
121 // Insert local interpolation matrix for each cell
122 auto cell_map = mesh->topology().index_map(tdim);
123 assert(cell_map);
124 std::int32_t num_cells = cell_map->size_local();
125 std::vector<T> Ae(Ab.size());
126 for (std::int32_t c = 0; c < num_cells; ++c)
127 {
128 std::copy(Ab.cbegin(), Ab.cend(), Ae.begin());
129 apply_inverse_dof_transform(Ae, cell_info, c, ndofs0);
130 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ae);
131 }
132}
133
149template <typename T, typename U>
151 U&& mat_set)
152{
153 // Get mesh
154 auto mesh = V0.mesh();
155 assert(mesh);
156
157 // Mesh dims
158 const int tdim = mesh->topology().dim();
159 const int gdim = mesh->geometry().dim();
160
161 // Get elements
162 std::shared_ptr<const FiniteElement> e0 = V0.element();
163 assert(e0);
164 std::shared_ptr<const FiniteElement> e1 = V1.element();
165 assert(e1);
166
167 std::span<const std::uint32_t> cell_info;
168 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
169 {
170 mesh->topology_mutable().create_entity_permutations();
171 cell_info = std::span(mesh->topology().get_cell_permutation_info());
172 }
173
174 // Get dofmaps
175 auto dofmap0 = V0.dofmap();
176 assert(dofmap0);
177 auto dofmap1 = V1.dofmap();
178 assert(dofmap1);
179
180 // Get block sizes and dof transformation operators
181 const int bs0 = e0->block_size();
182 const int bs1 = e1->block_size();
183 const auto apply_dof_transformation0
184 = e0->get_dof_transformation_function<double>(false, false, false);
185 const auto apply_inverse_dof_transform1
186 = e1->get_dof_transformation_function<T>(true, true, false);
187
188 // Get sizes of elements
189 const std::size_t space_dim0 = e0->space_dimension();
190 const std::size_t space_dim1 = e1->space_dimension();
191 const std::size_t dim0 = space_dim0 / bs0;
192 const std::size_t value_size_ref0 = e0->reference_value_size() / bs0;
193 const std::size_t value_size0 = e0->value_size() / bs0;
194 const std::size_t value_size1 = e1->value_size() / bs1;
195
196 // Get geometry data
197 const CoordinateElement& cmap = mesh->geometry().cmap();
199 = mesh->geometry().dofmap();
200 const std::size_t num_dofs_g = cmap.dim();
201 std::span<const double> x_g = mesh->geometry().x();
202
203 namespace stdex = std::experimental;
204 using mdspan2_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
205 using cmdspan2_t
206 = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
207 using cmdspan3_t
208 = stdex::mdspan<const double, stdex::dextents<std::size_t, 3>>;
209 using cmdspan4_t
210 = stdex::mdspan<const double, stdex::dextents<std::size_t, 4>>;
211 using mdspan3_t = stdex::mdspan<double, stdex::dextents<std::size_t, 3>>;
212
213 // Evaluate coordinate map basis at reference interpolation points
214 const auto [X, Xshape] = e1->interpolation_points();
215 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, Xshape[0]);
216 std::vector<double> phi_b(
217 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
218 cmdspan4_t phi(phi_b.data(), phi_shape);
219 cmap.tabulate(1, X, Xshape, phi_b);
220
221 // Evaluate V0 basis functions at reference interpolation points for V1
222 std::vector<double> basis_derivatives_reference0_b(Xshape[0] * dim0
223 * value_size_ref0);
224 cmdspan4_t basis_derivatives_reference0(basis_derivatives_reference0_b.data(),
225 1, Xshape[0], dim0, value_size_ref0);
226 e0->tabulate(basis_derivatives_reference0_b, X, Xshape, 0);
227
228 // Clamp values
229 std::transform(basis_derivatives_reference0_b.begin(),
230 basis_derivatives_reference0_b.end(),
231 basis_derivatives_reference0_b.begin(),
232 [atol = 1e-14](auto x)
233 { return std::abs(x) < atol ? 0.0 : x; });
234
235 // Create working arrays
236 std::vector<double> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
237 mdspan3_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0,
238 value_size_ref0);
239 std::vector<double> J_b(Xshape[0] * gdim * tdim);
240 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
241 std::vector<double> K_b(Xshape[0] * tdim * gdim);
242 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
243 std::vector<double> detJ(Xshape[0]);
244 std::vector<double> det_scratch(2 * tdim * gdim);
245
246 // Get the interpolation operator (matrix) `Pi` that maps a function
247 // evaluated at the interpolation points to the element degrees of
248 // freedom, i.e. dofs = Pi f_x
249 const auto [_Pi_1, pi_shape] = e1->interpolation_operator();
250 cmdspan2_t Pi_1(_Pi_1.data(), pi_shape);
251
252 bool interpolation_ident = e1->interpolation_ident();
253
254 namespace stdex = std::experimental;
255 using u_t = stdex::mdspan<double, stdex::dextents<std::size_t, 2>>;
256 using U_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
257 using J_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
258 using K_t = stdex::mdspan<const double, stdex::dextents<std::size_t, 2>>;
259 auto push_forward_fn0 = e0->basix_element().map_fn<u_t, U_t, J_t, K_t>();
260
261 // Basis values of Lagrange space unrolled for block size
262 // (num_quadrature_points, Lagrange dof, value_size)
263 std::vector<double> basis_values_b(Xshape[0] * bs0 * dim0 * e1->value_size());
264 mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0,
265 e1->value_size());
266 std::vector<double> mapped_values_b(Xshape[0] * bs0 * dim0
267 * e1->value_size());
268 mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0,
269 e1->value_size());
270
271 auto pull_back_fn1 = e1->basix_element().map_fn<u_t, U_t, K_t, J_t>();
272
273 std::vector<double> coord_dofs_b(num_dofs_g * gdim);
274 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
275 std::vector<double> basis0_b(Xshape[0] * dim0 * value_size0);
276 mdspan3_t basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
277
278 // Buffers
279 std::vector<T> Ab(space_dim0 * space_dim1);
280 std::vector<T> local1(space_dim1);
281
282 // Iterate over mesh and interpolate on each cell
283 auto cell_map = mesh->topology().index_map(tdim);
284 assert(cell_map);
285 std::int32_t num_cells = cell_map->size_local();
286 for (std::int32_t c = 0; c < num_cells; ++c)
287 {
288 // Get cell geometry (coordinate dofs)
289 auto x_dofs = x_dofmap.links(c);
290 for (std::size_t i = 0; i < x_dofs.size(); ++i)
291 {
292 for (std::size_t j = 0; j < gdim; ++j)
293 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
294 }
295
296 // Compute Jacobians and reference points for current cell
297 std::fill(J_b.begin(), J_b.end(), 0);
298 for (std::size_t p = 0; p < Xshape[0]; ++p)
299 {
300 auto dphi = stdex::submdspan(phi, std::pair(1, tdim + 1), p,
301 stdex::full_extent, 0);
302 auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
303 cmap.compute_jacobian(dphi, coord_dofs, _J);
304 auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
305 cmap.compute_jacobian_inverse(_J, _K);
306 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
307 }
308
309 // Copy evaluated basis on reference, apply DOF transformations, and
310 // push forward to physical element
311 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
312 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
313 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
314 basis_reference0(k0, k1, k2)
315 = basis_derivatives_reference0(0, k0, k1, k2);
316 for (std::size_t p = 0; p < Xshape[0]; ++p)
317 {
318 apply_dof_transformation0(
319 std::span(basis_reference0.data_handle() + p * dim0 * value_size_ref0,
320 dim0 * value_size_ref0),
321 cell_info, c, value_size_ref0);
322 }
323
324 for (std::size_t p = 0; p < basis0.extent(0); ++p)
325 {
326 auto _u
327 = stdex::submdspan(basis0, p, stdex::full_extent, stdex::full_extent);
328 auto _U = stdex::submdspan(basis_reference0, p, stdex::full_extent,
329 stdex::full_extent);
330 auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
331 auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
332 push_forward_fn0(_u, _U, _J, detJ[p], _K);
333 }
334
335 // Unroll basis function for input space for block size
336 for (std::size_t p = 0; p < Xshape[0]; ++p)
337 for (std::size_t i = 0; i < dim0; ++i)
338 for (std::size_t j = 0; j < value_size0; ++j)
339 for (int k = 0; k < bs0; ++k)
340 basis_values(p, i * bs0 + k, j * bs0 + k) = basis0(p, i, j);
341
342 // Pull back the physical values to the reference of output space
343 for (std::size_t p = 0; p < basis_values.extent(0); ++p)
344 {
345 auto _u = stdex::submdspan(basis_values, p, stdex::full_extent,
346 stdex::full_extent);
347 auto _U = stdex::submdspan(mapped_values, p, stdex::full_extent,
348 stdex::full_extent);
349 auto _K = stdex::submdspan(K, p, stdex::full_extent, stdex::full_extent);
350 auto _J = stdex::submdspan(J, p, stdex::full_extent, stdex::full_extent);
351 pull_back_fn1(_U, _u, _K, 1.0 / detJ[p], _J);
352 }
353
354 // Apply interpolation matrix to basis values of V0 at the
355 // interpolation points of V1
356 if (interpolation_ident)
357 {
358 stdex::mdspan<T, stdex::dextents<std::size_t, 3>> A(
359 Ab.data(), Xshape[0], e1->value_size(), space_dim0);
360 for (std::size_t i = 0; i < mapped_values.extent(0); ++i)
361 for (std::size_t j = 0; j < mapped_values.extent(1); ++j)
362 for (std::size_t k = 0; k < mapped_values.extent(2); ++k)
363 A(i, k, j) = mapped_values(i, j, k);
364 }
365 else
366 {
367 for (std::size_t i = 0; i < mapped_values.extent(1); ++i)
368 {
369 auto values = stdex::submdspan(mapped_values, stdex::full_extent, i,
370 stdex::full_extent);
371 impl::interpolation_apply(Pi_1, values, std::span(local1), bs1);
372 for (std::size_t j = 0; j < local1.size(); j++)
373 Ab[space_dim0 * j + i] = local1[j];
374 }
375 }
376
377 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
378 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);
379 }
380}
381
382} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:32
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Compute Jacobian for a cell with given geometry using the basis functions and first order derivatives...
Definition: CoordinateElement.h:100
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:109
std::array< std::size_t, 4 > tabulate_shape(std::size_t nd, std::size_t num_points) const
Shape of array to fill when calling FiniteElement::tabulate
Definition: CoordinateElement.cpp:42
void tabulate(int nd, std::span< const double > X, std::array< std::size_t, 2 > shape, std::span< double > basis) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:47
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:183
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:126
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:31
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition: FunctionSpace.cpp:253
std::shared_ptr< const mesh::Mesh > mesh() const
The mesh.
Definition: FunctionSpace.cpp:246
std::shared_ptr< const FiniteElement > element() const
The finite element.
Definition: FunctionSpace.cpp:248
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
void interpolation_matrix(const FunctionSpace &V0, const FunctionSpace &V1, U &&mat_set)
Assemble an interpolation operator matrix.
Definition: discreteoperators.h:150
void discrete_gradient(const FunctionSpace &V0, const FunctionSpace &V1, U &&mat_set)
Assemble a discrete gradient operator.
Definition: discreteoperators.h:48