DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
discreteoperators.h
1// Copyright (C) 2015-2025 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 <algorithm>
13#include <array>
14#include <concepts>
15#include <dolfinx/common/IndexMap.h>
16#include <dolfinx/common/math.h>
17#include <dolfinx/la/utils.h>
18#include <dolfinx/mesh/Mesh.h>
19#include <memory>
20#include <span>
21#include <vector>
22
23namespace dolfinx::fem
24{
25
85template <std::floating_point T, dolfinx::scalar U = T>
87 la::MatSet<U> auto&& mat_set)
88{
89 // Get mesh
90 auto mesh = V0.mesh();
91 assert(mesh);
92 assert(V1.mesh());
93 if (mesh != V1.mesh())
94 throw std::runtime_error("Meshes must be the same.");
95
96 if (mesh->geometry().dim() != 3)
97 throw std::runtime_error("Geometric must be equal to 3..");
98 if (mesh->geometry().dim() != mesh->topology()->dim())
99 {
100 throw std::runtime_error(
101 "Geometric and topological dimensions must be equal.");
102 }
103 constexpr int gdim = 3;
104
105 // Get elements
106 std::shared_ptr<const FiniteElement<T>> e0 = V0.element();
107 assert(e0);
108 if (e0->map_type() != basix::maps::type::covariantPiola)
109 {
110 throw std::runtime_error(
111 "Finite element for parent space must be covariant Piola.");
112 }
113
114 std::shared_ptr<const FiniteElement<T>> e1 = V1.element();
115 assert(e1);
116 if (e1->map_type() != basix::maps::type::contravariantPiola)
117 {
118 throw std::runtime_error(
119 "Finite element for target space must be contracovariant Piola.");
120 }
121
122 // Get cell orientation information
123 std::span<const std::uint32_t> cell_info;
124 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
125 {
126 mesh->topology_mutable()->create_entity_permutations();
127 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
128 }
129
130 // Get dofmaps
131 auto dofmap0 = V0.dofmap();
132 assert(dofmap0);
133 auto dofmap1 = V1.dofmap();
134 assert(dofmap1);
135
136 // Get dof transformation operators
137 auto apply_dof_transformation0
138 = e0->template dof_transformation_fn<T>(doftransform::standard, false);
139 auto apply_inverse_dof_transform1 = e1->template dof_transformation_fn<U>(
141
142 // Get sizes of elements
143 const std::size_t space_dim0 = e0->space_dimension();
144 const std::size_t space_dim1 = e1->space_dimension();
145 if (e0->reference_value_size() != 3)
146 throw std::runtime_error("Value size for parent space should be 3.");
147 if (e1->reference_value_size() != 3)
148 throw std::runtime_error("Value size for target space should be 3.");
149
150 // Get the V1 reference interpolation points
151 const auto [X, Xshape] = e1->interpolation_points();
152
153 // Get/compute geometry map and evaluate at interpolation points
154 const CoordinateElement<T>& cmap = mesh->geometry().cmap();
155 auto x_dofmap = mesh->geometry().dofmap();
156 const std::size_t num_dofs_g = cmap.dim();
157 std::span<const T> x_g = mesh->geometry().x();
158 std::array<std::size_t, 4> Phi_g_shape = cmap.tabulate_shape(1, Xshape[0]);
159 std::vector<T> Phi_g_b(std::reduce(Phi_g_shape.begin(), Phi_g_shape.end(), 1,
160 std::multiplies{}));
161 // (derivative, point index, phi, component)
162 md::mdspan<const T, md::extents<std::size_t, gdim + 1, md::dynamic_extent,
163 md::dynamic_extent, 1>>
164 Phi_g(Phi_g_b.data(), Phi_g_shape);
165 cmap.tabulate(1, X, Xshape, Phi_g_b);
166
167 // Geometry data structures
168 std::vector<T> coord_dofs_b(num_dofs_g * gdim);
169 md::mdspan<T, md::extents<std::size_t, md::dynamic_extent, 3>> coord_dofs(
170 coord_dofs_b.data(), num_dofs_g, gdim);
171 std::vector<T> J_b(Xshape[0] * gdim * gdim);
172 md::mdspan<T, md::extents<std::size_t, md::dynamic_extent, 3, 3>> J(
173 J_b.data(), Xshape[0], gdim, gdim);
174 std::vector<T> K_b(Xshape[0] * gdim * gdim);
175 md::mdspan<T, typename decltype(J)::extents_type> K(K_b.data(), J.extents());
176 std::vector<T> detJ(Xshape[0]);
177 std::vector<T> det_scratch(2 * gdim * gdim);
178
179 // Evaluate V0 basis function derivatives at reference interpolation
180 // points for V1, (deriv, pt_idx, phi (dof), comp)
181 const auto [Phi0_b, Phi0_shape] = e0->tabulate(X, Xshape, 1);
182 md::mdspan<const T, md::extents<std::size_t, 4, md::dynamic_extent,
183 md::dynamic_extent, 3>>
184 Phi0(Phi0_b.data(), Phi0_shape);
185
186 // Create working arrays, (point, phi (dof), deriv, comp)
187 md::extents<std::size_t, md::dynamic_extent, md::dynamic_extent, 3, 3>
188 dPhi0_ext(Phi0.extent(1), Phi0.extent(2), Phi0.extent(0) - 1,
189 Phi0.extent(3));
190 std::vector<T> dPhi0_b(dPhi0_ext.extent(0) * dPhi0_ext.extent(1)
191 * dPhi0_ext.extent(2) * dPhi0_ext.extent(3));
192 md::mdspan<T, decltype(dPhi0_ext)> dPhi0(dPhi0_b.data(), dPhi0_ext);
193
194 // Get the interpolation operator (matrix) Pi that maps a function
195 // evaluated at the interpolation points to the V1 element degrees of
196 // freedom, i.e. dofs = Pi f_x
197 const auto [Pi1_b, pi_shape] = e1->interpolation_operator();
198 md::mdspan<const T, md::dextents<std::size_t, 2>> Pi_1(Pi1_b.data(),
199 pi_shape);
200
201 // curl data structure, (pt_idx, phi (dof), comp)
202 std::vector<T> curl_b(dPhi0.extent(0) * dPhi0.extent(1) * dPhi0.extent(3));
203 md::mdspan<
204 T, md::extents<std::size_t, md::dynamic_extent, md::dynamic_extent, 3>>
205 curl(curl_b.data(), dPhi0.extent(0), dPhi0.extent(1), dPhi0.extent(3));
206
207 std::vector<U> Ab(space_dim0 * space_dim1);
208 std::vector<U> local1(space_dim1);
209
210 // Iterate over mesh and interpolate on each cell
211 assert(mesh->topology()->index_map(gdim));
212 for (std::int32_t c = 0; c < mesh->topology()->index_map(gdim)->size_local();
213 ++c)
214 {
215 // Get cell geometry (coordinate dofs)
216 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
217 for (std::size_t i = 0; i < x_dofs.size(); ++i)
218 for (std::size_t j = 0; j < gdim; ++j)
219 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
220
221 // Compute Jacobians at reference points for current cell
222 std::ranges::fill(J_b, 0);
223 for (std::size_t p = 0; p < Xshape[0]; ++p)
224 {
225 auto dPhi_g
226 = md::submdspan(Phi_g, std::pair(1, gdim + 1), p, md::full_extent, 0);
227 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
228 cmap.compute_jacobian(dPhi_g, coord_dofs, _J);
229 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
230 cmap.compute_jacobian_inverse(_J, _K);
231 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
232 }
233
234 // TODO: re-order loops and/or re-pack Phi0 to allow a simple flat
235 // copy?
236
237 // Copy (d)Phi0 (on reference) and apply DOF transformation
238 // Phi0: (deriv, pt_idx, phi (dof), comp)
239 // dPhi0: (pt_idx, phi (dov), deriv, comp)
240 for (std::size_t p = 0; p < Phi0.extent(1); ++p) // point
241 for (std::size_t phi = 0; phi < Phi0.extent(2); ++phi) // phi_i
242 for (std::size_t d = 0; d < Phi0.extent(3); ++d) // Comp. of phi
243 for (std::size_t dx = 0; dx < dPhi0.extent(2); ++dx) // dx
244 dPhi0(p, phi, dx, d) = Phi0(dx + 1, p, phi, d);
245
246 for (std::size_t p = 0; p < dPhi0.extent(0); ++p) // point
247 {
248 // Size: num_phi * num_derivs * num_components
249 std::size_t size = dPhi0.extent(1) * dPhi0.extent(2) * dPhi0.extent(3);
250 std::size_t offset = p * size; // Offset for point p
251
252 // Shape: (num_phi , (value_size * num_derivs))
253 apply_dof_transformation0(std::span(dPhi0.data_handle() + offset, size),
254 cell_info, c,
255 dPhi0.extent(2) * dPhi0.extent(3));
256 }
257
258 // Compute curl
259 // dPhi0: (pt_idx, phi_idx, deriv, comp)
260 // curl: (pt_idx, phi_idx, comp)
261 for (std::size_t p = 0; p < curl.extent(0); ++p) // point
262 {
263 for (std::size_t i = 0; i < curl.extent(1); ++i) // phi_i
264 {
265 curl(p, i, 0) = dPhi0(p, i, 1, 2) - dPhi0(p, i, 2, 1);
266 curl(p, i, 1) = dPhi0(p, i, 2, 0) - dPhi0(p, i, 0, 2);
267 curl(p, i, 2) = dPhi0(p, i, 0, 1) - dPhi0(p, i, 1, 0);
268 }
269 }
270
271 // Apply interpolation matrix to basis derivatives values of V0 at
272 // the interpolation points of V1
273 for (std::size_t i = 0; i < curl.extent(1); ++i)
274 {
275 auto values = md::submdspan(curl, md::full_extent, i, md::full_extent);
276 impl::interpolation_apply(Pi_1, values, std::span(local1), 1);
277 for (std::size_t j = 0; j < local1.size(); ++j)
278 Ab[space_dim0 * j + i] = local1[j];
279 }
280
281 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
282 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);
283 }
284}
285
312template <dolfinx::scalar T, std::floating_point U = dolfinx::scalar_value_t<T>>
314 std::pair<std::reference_wrapper<const FiniteElement<U>>,
315 std::reference_wrapper<const DofMap>>
316 V0,
317 std::pair<std::reference_wrapper<const FiniteElement<U>>,
318 std::reference_wrapper<const DofMap>>
319 V1,
320 auto&& mat_set)
321{
322 auto& e0 = V0.first.get();
323 const DofMap& dofmap0 = V0.second.get();
324 auto& e1 = V1.first.get();
325 const DofMap& dofmap1 = V1.second.get();
326
327 using cmdspan2_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
328 using cmdspan4_t = md::mdspan<const U, md::dextents<std::size_t, 4>>;
329
330 // Check elements
331 if (e0.map_type() != basix::maps::type::identity)
332 throw std::runtime_error("Wrong finite element space for V0.");
333 if (e0.block_size() != 1)
334 throw std::runtime_error("Block size is greater than 1 for V0.");
335 if (e0.reference_value_size() != 1)
336 throw std::runtime_error("Wrong value size for V0.");
337
338 if (e1.map_type() != basix::maps::type::covariantPiola)
339 throw std::runtime_error("Wrong finite element space for V1.");
340 if (e1.block_size() != 1)
341 throw std::runtime_error("Block size is greater than 1 for V1.");
342
343 // Get V1 (H(curl)) space interpolation points
344 const auto [X, Xshape] = e1.interpolation_points();
345
346 // Tabulate first derivatives of Lagrange space at V1 interpolation
347 // points
348 const int ndofs0 = e0.space_dimension();
349 const int tdim = topology.dim();
350 std::vector<U> phi0_b((tdim + 1) * Xshape[0] * ndofs0 * 1);
351 cmdspan4_t phi0(phi0_b.data(), tdim + 1, Xshape[0], ndofs0, 1);
352 e0.tabulate(phi0_b, X, Xshape, 1);
353
354 // Reshape Lagrange basis derivatives as a matrix of shape (tdim *
355 // num_points, num_dofs_per_cell)
356 cmdspan2_t dphi_reshaped(
357 phi0_b.data() + phi0.extent(3) * phi0.extent(2) * phi0.extent(1),
358 tdim * phi0.extent(1), phi0.extent(2));
359
360 // Get inverse DOF transform function
361 auto apply_inverse_dof_transform = e1.template dof_transformation_fn<T>(
363
364 // Generate cell permutations
366 const std::vector<std::uint32_t>& cell_info
367 = topology.get_cell_permutation_info();
368
369 // Create element kernel function
370
371 // Build the element interpolation matrix
372 std::vector<T> Ab(e1.space_dimension() * ndofs0);
373 {
374 md::mdspan<T, md::dextents<std::size_t, 2>> A(Ab.data(),
375 e1.space_dimension(), ndofs0);
376 const auto [Pi, shape] = e1.interpolation_operator();
377 cmdspan2_t _Pi(Pi.data(), shape);
378 math::dot(_Pi, dphi_reshaped, A);
379 }
380
381 // Insert local interpolation matrix for each cell
382 auto cell_map = topology.index_map(tdim);
383 assert(cell_map);
384 std::int32_t num_cells = cell_map->size_local();
385 std::vector<T> Ae(Ab.size());
386 for (std::int32_t c = 0; c < num_cells; ++c)
387 {
388 std::ranges::copy(Ab, Ae.begin());
389 apply_inverse_dof_transform(Ae, cell_info, c, ndofs0);
390 mat_set(dofmap1.cell_dofs(c), dofmap0.cell_dofs(c), Ae);
391 }
392}
393
409template <dolfinx::scalar T, std::floating_point U>
411 const FunctionSpace<U>& V1, auto&& mat_set)
412{
413 // Get mesh
414 auto mesh = V0.mesh();
415 assert(mesh);
416
417 // Mesh dims
418 const int tdim = mesh->topology()->dim();
419 const int gdim = mesh->geometry().dim();
420
421 // Get elements
422 std::shared_ptr<const FiniteElement<U>> e0 = V0.element();
423 assert(e0);
424 std::shared_ptr<const FiniteElement<U>> e1 = V1.element();
425 assert(e1);
426
427 std::span<const std::uint32_t> cell_info;
428 if (e1->needs_dof_transformations() or e0->needs_dof_transformations())
429 {
430 mesh->topology_mutable()->create_entity_permutations();
431 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
432 }
433
434 // Get dofmaps
435 auto dofmap0 = V0.dofmap();
436 assert(dofmap0);
437 auto dofmap1 = V1.dofmap();
438 assert(dofmap1);
439
440 // Get block sizes and dof transformation operators
441 const int bs0 = e0->block_size();
442 const int bs1 = e1->block_size();
443 auto apply_dof_transformation0
444 = e0->template dof_transformation_fn<U>(doftransform::standard, false);
445 auto apply_inverse_dof_transform1 = e1->template dof_transformation_fn<T>(
447
448 // Get sizes of elements
449 const std::size_t space_dim0 = e0->space_dimension();
450 const std::size_t space_dim1 = e1->space_dimension();
451 const std::size_t dim0 = space_dim0 / bs0;
452 const std::size_t value_size_ref0 = e0->reference_value_size();
453 const std::size_t value_size0 = V0.element()->reference_value_size();
454
455 // Get geometry data
456 const CoordinateElement<U>& cmap = mesh->geometry().cmap();
457 auto x_dofmap = mesh->geometry().dofmap();
458 const std::size_t num_dofs_g = cmap.dim();
459 std::span<const U> x_g = mesh->geometry().x();
460
461 using mdspan2_t = md::mdspan<U, md::dextents<std::size_t, 2>>;
462 using cmdspan2_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
463 using cmdspan4_t = md::mdspan<const U, md::dextents<std::size_t, 4>>;
464 using mdspan3_t = md::mdspan<U, md::dextents<std::size_t, 3>>;
465
466 // Evaluate coordinate map basis at reference interpolation points
467 const auto [X, Xshape] = e1->interpolation_points();
468 std::array<std::size_t, 4> phi_shape = cmap.tabulate_shape(1, Xshape[0]);
469 std::vector<U> phi_b(
470 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
471 cmdspan4_t phi(phi_b.data(), phi_shape);
472 cmap.tabulate(1, X, Xshape, phi_b);
473
474 // Evaluate V0 basis functions at reference interpolation points for V1
475 std::vector<U> basis_derivatives_reference0_b(Xshape[0] * dim0
476 * value_size_ref0);
477 cmdspan4_t basis_derivatives_reference0(basis_derivatives_reference0_b.data(),
478 1, Xshape[0], dim0, value_size_ref0);
479 e0->tabulate(basis_derivatives_reference0_b, X, Xshape, 0);
480
481 // Clamp values
482 std::ranges::transform(
483 basis_derivatives_reference0_b, basis_derivatives_reference0_b.begin(),
484 [atol = 1e-14](auto x) { return std::abs(x) < atol ? 0.0 : x; });
485
486 // Create working arrays
487 std::vector<U> basis_reference0_b(Xshape[0] * dim0 * value_size_ref0);
488 mdspan3_t basis_reference0(basis_reference0_b.data(), Xshape[0], dim0,
489 value_size_ref0);
490 std::vector<U> J_b(Xshape[0] * gdim * tdim);
491 mdspan3_t J(J_b.data(), Xshape[0], gdim, tdim);
492 std::vector<U> K_b(Xshape[0] * tdim * gdim);
493 mdspan3_t K(K_b.data(), Xshape[0], tdim, gdim);
494 std::vector<U> detJ(Xshape[0]);
495 std::vector<U> det_scratch(2 * tdim * gdim);
496
497 // Get the interpolation operator (matrix) `Pi` that maps a function
498 // evaluated at the interpolation points to the element degrees of
499 // freedom, i.e. dofs = Pi f_x
500 const auto [_Pi_1, pi_shape] = e1->interpolation_operator();
501 cmdspan2_t Pi_1(_Pi_1.data(), pi_shape);
502
503 bool interpolation_ident = e1->interpolation_ident();
504
505 using u_t = md::mdspan<U, md::dextents<std::size_t, 2>>;
506 using U_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
507 using J_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
508 using K_t = md::mdspan<const U, md::dextents<std::size_t, 2>>;
509 auto push_forward_fn0
510 = e0->basix_element().template map_fn<u_t, U_t, J_t, K_t>();
511
512 // Basis values of Lagrange space unrolled for block size
513 // (num_quadrature_points, Lagrange dof, value_size)
514 std::vector<U> basis_values_b(Xshape[0] * bs0 * dim0
515 * V1.element()->value_size());
516 mdspan3_t basis_values(basis_values_b.data(), Xshape[0], bs0 * dim0,
517 V1.element()->value_size());
518 std::vector<U> mapped_values_b(Xshape[0] * bs0 * dim0
519 * V1.element()->value_size());
520 mdspan3_t mapped_values(mapped_values_b.data(), Xshape[0], bs0 * dim0,
521 V1.element()->value_size());
522
523 auto pull_back_fn1
524 = e1->basix_element().template map_fn<u_t, U_t, K_t, J_t>();
525
526 std::vector<U> coord_dofs_b(num_dofs_g * gdim);
527 mdspan2_t coord_dofs(coord_dofs_b.data(), num_dofs_g, gdim);
528 std::vector<U> basis0_b(Xshape[0] * dim0 * value_size0);
529 mdspan3_t basis0(basis0_b.data(), Xshape[0], dim0, value_size0);
530
531 // Buffers
532 std::vector<T> Ab(space_dim0 * space_dim1);
533 std::vector<T> local1(space_dim1);
534
535 // Iterate over mesh and interpolate on each cell
536 auto cell_map = mesh->topology()->index_map(tdim);
537 assert(cell_map);
538 std::int32_t num_cells = cell_map->size_local();
539 for (std::int32_t c = 0; c < num_cells; ++c)
540 {
541 // Get cell geometry (coordinate dofs)
542 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
543 for (std::size_t i = 0; i < x_dofs.size(); ++i)
544 {
545 for (int j = 0; j < gdim; ++j)
546 coord_dofs(i, j) = x_g[3 * x_dofs[i] + j];
547 }
548
549 // Compute Jacobians and reference points for current cell
550 std::ranges::fill(J_b, 0);
551 for (std::size_t p = 0; p < Xshape[0]; ++p)
552 {
553 auto dphi
554 = md::submdspan(phi, std::pair(1, tdim + 1), p, md::full_extent, 0);
555 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
556 cmap.compute_jacobian(dphi, coord_dofs, _J);
557 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
558 cmap.compute_jacobian_inverse(_J, _K);
559 detJ[p] = cmap.compute_jacobian_determinant(_J, det_scratch);
560 }
561
562 // Copy evaluated basis on reference, apply DOF transformations, and
563 // push forward to physical element
564 for (std::size_t k0 = 0; k0 < basis_reference0.extent(0); ++k0)
565 for (std::size_t k1 = 0; k1 < basis_reference0.extent(1); ++k1)
566 for (std::size_t k2 = 0; k2 < basis_reference0.extent(2); ++k2)
567 basis_reference0(k0, k1, k2)
568 = basis_derivatives_reference0(0, k0, k1, k2);
569 for (std::size_t p = 0; p < Xshape[0]; ++p)
570 {
571 apply_dof_transformation0(
572 std::span(basis_reference0.data_handle() + p * dim0 * value_size_ref0,
573 dim0 * value_size_ref0),
574 cell_info, c, value_size_ref0);
575 }
576
577 for (std::size_t p = 0; p < basis0.extent(0); ++p)
578 {
579 auto _u = md::submdspan(basis0, p, md::full_extent, md::full_extent);
580 auto _U = md::submdspan(basis_reference0, p, md::full_extent,
581 md::full_extent);
582 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
583 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
584 push_forward_fn0(_u, _U, _J, detJ[p], _K);
585 }
586
587 // Unroll basis function for input space for block size
588 for (std::size_t p = 0; p < Xshape[0]; ++p)
589 for (std::size_t i = 0; i < dim0; ++i)
590 for (std::size_t j = 0; j < value_size0; ++j)
591 for (int k = 0; k < bs0; ++k)
592 basis_values(p, i * bs0 + k, j * bs0 + k) = basis0(p, i, j);
593
594 // Pull back the physical values to the reference of output space
595 for (std::size_t p = 0; p < basis_values.extent(0); ++p)
596 {
597 auto _u
598 = md::submdspan(basis_values, p, md::full_extent, md::full_extent);
599 auto _U
600 = md::submdspan(mapped_values, p, md::full_extent, md::full_extent);
601 auto _K = md::submdspan(K, p, md::full_extent, md::full_extent);
602 auto _J = md::submdspan(J, p, md::full_extent, md::full_extent);
603 pull_back_fn1(_U, _u, _K, 1.0 / detJ[p], _J);
604 }
605
606 // Apply interpolation matrix to basis values of V0 at the
607 // interpolation points of V1
608 if (interpolation_ident)
609 {
610 md::mdspan<T, md::dextents<std::size_t, 3>> A(
611 Ab.data(), Xshape[0], V1.element()->value_size(), space_dim0);
612 for (std::size_t i = 0; i < mapped_values.extent(0); ++i)
613 for (std::size_t j = 0; j < mapped_values.extent(1); ++j)
614 for (std::size_t k = 0; k < mapped_values.extent(2); ++k)
615 A(i, k, j) = mapped_values(i, j, k);
616 }
617 else
618 {
619 for (std::size_t i = 0; i < mapped_values.extent(1); ++i)
620 {
621 auto values
622 = md::submdspan(mapped_values, md::full_extent, i, md::full_extent);
623 impl::interpolation_apply(Pi_1, values, std::span(local1), bs1);
624 for (std::size_t j = 0; j < local1.size(); j++)
625 Ab[space_dim0 * j + i] = local1[j];
626 }
627 }
628
629 apply_inverse_dof_transform1(Ab, cell_info, c, space_dim0);
630 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ab);
631 }
632}
633
634} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Definition CoordinateElement.h:38
static void compute_jacobian(const U &dphi, const V &cell_geometry, W &&J)
Definition CoordinateElement.h:133
void tabulate(int nd, std::span< const T > X, std::array< std::size_t, 2 > shape, std::span< T > basis) const
Evaluate basis values and derivatives at set of points.
Definition CoordinateElement.cpp:55
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition CoordinateElement.h:142
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 tabulate.
Definition CoordinateElement.cpp:48
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:205
static double compute_jacobian_determinant(const U &J, std::span< typename U::value_type > w)
Compute the determinant of the Jacobian.
Definition CoordinateElement.h:159
Degree-of-freedom map.
Definition DofMap.h:73
std::span< const std::int32_t > cell_dofs(std::int32_t c) const
Local-to-global mapping of dofs on a cell.
Definition DofMap.h:127
Model of a finite element.
Definition FiniteElement.h:57
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:361
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:336
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:342
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:46
std::shared_ptr< const common::IndexMap > index_map(int dim) const
Get the IndexMap that described the parallel distribution of the mesh entities.
Definition Topology.cpp:823
void create_entity_permutations()
Compute entity permutations and reflections.
Definition Topology.cpp:990
const std::vector< std::uint32_t > & get_cell_permutation_info() const
Returns the permutation information.
Definition Topology.cpp:852
int dim() const noexcept
Topological dimension of the mesh.
Definition Topology.cpp:785
Matrix accumulate/set concept for functions that can be used in assemblers to accumulate or set value...
Definition utils.h:28
Finite element method functionality.
Definition assemble_expression_impl.h:23
void discrete_curl(const FunctionSpace< T > &V0, const FunctionSpace< T > &V1, la::MatSet< U > auto &&mat_set)
Assemble a discrete curl operator.
Definition discreteoperators.h:86
void discrete_gradient(mesh::Topology &topology, std::pair< std::reference_wrapper< const FiniteElement< U > >, std::reference_wrapper< const DofMap > > V0, std::pair< std::reference_wrapper< const FiniteElement< U > >, std::reference_wrapper< const DofMap > > V1, auto &&mat_set)
Assemble a discrete gradient operator.
Definition discreteoperators.h:313
@ inverse_transpose
Transpose inverse.
Definition FiniteElement.h:30
@ standard
Standard.
Definition FiniteElement.h:27
void interpolation_matrix(const FunctionSpace< U > &V0, const FunctionSpace< U > &V1, auto &&mat_set)
Assemble an interpolation operator matrix.
Definition discreteoperators.h:410
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32