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 <vector>
19#include <xtensor/xadapt.hpp>
20#include <xtensor/xbuilder.hpp>
21#include <xtensor/xindex_view.hpp>
22#include <xtensor/xtensor.hpp>
23#include <xtl/xspan.hpp>
24
25namespace dolfinx::fem
26{
27
51template <typename T, typename U>
53 const fem::FunctionSpace& V1, U&& mat_set)
54{
55 // Get mesh
56 std::shared_ptr<const mesh::Mesh> mesh = V1.mesh();
57 assert(mesh);
58
59 // Check spaces
60 std::shared_ptr<const FiniteElement> e0 = V0.element();
61 assert(e0);
62 if (e0->map_type() != basix::maps::type::identity)
63 throw std::runtime_error("Wrong finite element space for V0.");
64 if (e0->block_size() != 1)
65 throw std::runtime_error("Block size is greather than 1 for V0.");
66 if (e0->reference_value_size() != 1)
67 throw std::runtime_error("Wrong value size for V0.");
68
69 std::shared_ptr<const FiniteElement> e1 = V1.element();
70 assert(e1);
71 if (e1->map_type() != basix::maps::type::covariantPiola)
72 throw std::runtime_error("Wrong finite element space for V1.");
73 if (e1->block_size() != 1)
74 throw std::runtime_error("Block size is greather than 1 for V1.");
75
76 // Get V0 (H(curl)) space interpolation points
77 const xt::xtensor<double, 2> X = e1->interpolation_points();
78
79 // Tabulate first order derivatives of Lagrange space at H(curl)
80 // interpolation points
81 const int ndofs0 = e0->space_dimension();
82 const int tdim = mesh->topology().dim();
83 xt::xtensor<double, 4> phi0
84 = xt::empty<double>({tdim + 1, int(X.shape(0)), ndofs0, 1});
85 e0->tabulate(phi0, X, 1);
86
87 // Reshape lagrange basis derivatives as a matrix of shape (tdim *
88 // num_points, num_dofs_per_cell)
89 auto dphi0 = xt::view(phi0, xt::xrange(std::size_t(1), phi0.shape(0)),
90 xt::all(), xt::all(), 0);
91 auto dphi_reshaped
92 = xt::reshape_view(dphi0, {tdim * phi0.shape(1), phi0.shape(2)});
93
94 // Get inverse DOF transform function
95 auto apply_inverse_dof_transform
96 = e1->get_dof_transformation_function<T>(true, true, false);
97
98 // Generate cell permutations
99 mesh->topology_mutable().create_entity_permutations();
100 const std::vector<std::uint32_t>& cell_info
101 = mesh->topology().get_cell_permutation_info();
102
103 // Create element kernel function
104 std::shared_ptr<const DofMap> dofmap0 = V0.dofmap();
105 assert(dofmap0);
106 std::shared_ptr<const DofMap> dofmap1 = V1.dofmap();
107 assert(dofmap1);
108
109 // Build the element interpolation matrix
110 std::vector<T> A(e1->space_dimension() * ndofs0);
111 {
112 auto _A = xt::adapt(A, std::vector<int>{e1->space_dimension(), ndofs0});
113 const xt::xtensor<double, 2> Pi = e1->interpolation_operator();
114 math::dot(Pi, dphi_reshaped, _A);
115 }
116
117 // Insert local interpolation matrix for each cell
118 auto cell_map = mesh->topology().index_map(tdim);
119 assert(cell_map);
120 std::int32_t num_cells = cell_map->size_local();
121 std::vector<T> Ae(A.size());
122 for (std::int32_t c = 0; c < num_cells; ++c)
123 {
124 std::copy(A.cbegin(), A.cend(), Ae.begin());
125 apply_inverse_dof_transform(Ae, cell_info, c, ndofs0);
126 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), Ae);
127 }
128}
129
145template <typename T, typename U>
147 const fem::FunctionSpace& V1, U&& mat_set)
148{
149 // Get mesh
150 auto mesh = V0.mesh();
151 assert(mesh);
152
153 // Mesh dims
154 const int tdim = mesh->topology().dim();
155 const int gdim = mesh->geometry().dim();
156
157 // Get elements
158 std::shared_ptr<const FiniteElement> element0 = V0.element();
159 assert(element0);
160 std::shared_ptr<const FiniteElement> element1 = V1.element();
161 assert(element1);
162
163 xtl::span<const std::uint32_t> cell_info;
164 if (element1->needs_dof_transformations()
165 or element0->needs_dof_transformations())
166 {
167 mesh->topology_mutable().create_entity_permutations();
168 cell_info = xtl::span(mesh->topology().get_cell_permutation_info());
169 }
170
171 // Get dofmaps
172 auto dofmap0 = V0.dofmap();
173 auto dofmap1 = V1.dofmap();
174
175 // Get block sizes and dof transformation operators
176 const int bs0 = element0->block_size();
177 const int bs1 = element1->block_size();
178 const auto apply_dof_transformation0
179 = element0->get_dof_transformation_function<double>(false, false, false);
180 const auto apply_inverse_dof_transform1
181 = element1->get_dof_transformation_function<T>(true, true, false);
182
183 // Get sizes of elements
184 const std::size_t space_dim0 = element0->space_dimension();
185 const std::size_t space_dim1 = element1->space_dimension();
186 const std::size_t dim0 = space_dim0 / bs0;
187 const std::size_t value_size_ref0 = element0->reference_value_size() / bs0;
188 const std::size_t value_size0 = element0->value_size() / bs0;
189 const std::size_t value_size1 = element1->value_size() / bs1;
190
191 // Get geometry data
192 const fem::CoordinateElement& cmap = mesh->geometry().cmap();
194 = mesh->geometry().dofmap();
195 const std::size_t num_dofs_g = cmap.dim();
196 xtl::span<const double> x_g = mesh->geometry().x();
197
198 // Evaluate coordinate map basis at reference interpolation points
199 const xt::xtensor<double, 2> X = element1->interpolation_points();
200 xt::xtensor<double, 4> phi(cmap.tabulate_shape(1, X.shape(0)));
201 cmap.tabulate(1, X, phi);
202 xt::xtensor<double, 2> dphi
203 = xt::view(phi, xt::range(1, tdim + 1), 0, xt::all(), 0);
204
205 // Evaluate V0 basis functions at reference interpolation points for V1
206 xt::xtensor<double, 4> basis_derivatives_reference0(
207 {1, X.shape(0), dim0, value_size_ref0});
208 element0->tabulate(basis_derivatives_reference0, X, 0);
209
210 double rtol = 1e-14;
211 double atol = 1e-14;
212 auto inds = xt::isclose(basis_derivatives_reference0, 0.0, rtol, atol);
213 xt::filtration(basis_derivatives_reference0, inds) = 0.0;
214
215 // Create working arrays
216 xt::xtensor<double, 3> basis_reference0({X.shape(0), dim0, value_size_ref0});
217 xt::xtensor<double, 3> J({X.shape(0), gdim, tdim});
218 xt::xtensor<double, 3> K({X.shape(0), tdim, gdim});
219 std::vector<double> detJ(X.shape(0));
220
221 // Get the interpolation operator (matrix) `Pi` that maps a function
222 // evaluated at the interpolation points to the element degrees of
223 // freedom, i.e. dofs = Pi f_x
224 const xt::xtensor<double, 2>& Pi_1 = element1->interpolation_operator();
225 bool interpolation_ident = element1->interpolation_ident();
226
227 using u_t = xt::xview<decltype(basis_reference0)&, std::size_t,
228 xt::xall<std::size_t>, xt::xall<std::size_t>>;
229 using U_t = xt::xview<decltype(basis_reference0)&, std::size_t,
230 xt::xall<std::size_t>, xt::xall<std::size_t>>;
231 using J_t = xt::xview<decltype(J)&, std::size_t, xt::xall<std::size_t>,
232 xt::xall<std::size_t>>;
233 using K_t = xt::xview<decltype(K)&, std::size_t, xt::xall<std::size_t>,
234 xt::xall<std::size_t>>;
235 auto push_forward_fn0 = element0->map_fn<u_t, U_t, J_t, K_t>();
236
237 // Basis values of Lagrange space unrolled for block size
238 // (num_quadrature_points, Lagrange dof, value_size)
239 xt::xtensor<double, 3> basis_values = xt::zeros<double>(
240 {X.shape(0), bs0 * dim0, (std::size_t)element1->value_size()});
241 xt::xtensor<double, 3> mapped_values(
242 {X.shape(0), bs0 * dim0, (std::size_t)element1->value_size()});
243
244 using u1_t = xt::xview<decltype(basis_values)&, std::size_t,
245 xt::xall<std::size_t>, xt::xall<std::size_t>>;
246 using U1_t = xt::xview<decltype(mapped_values)&, std::size_t,
247 xt::xall<std::size_t>, xt::xall<std::size_t>>;
248 auto pull_back_fn1 = element1->map_fn<U1_t, u1_t, K_t, J_t>();
249
250 xt::xtensor<double, 2> coordinate_dofs({num_dofs_g, 3});
251 xt::xtensor<double, 3> basis0({X.shape(0), dim0, value_size0});
252 std::vector<T> A(space_dim0 * space_dim1);
253 std::vector<T> local1(space_dim1);
254
255 std::vector<std::size_t> shape
256 = {X.shape(0), (std::size_t)element1->value_size(), space_dim0};
257 auto _A = xt::adapt(A, shape);
258
259 // Iterate over mesh and interpolate on each cell
260 auto cell_map = mesh->topology().index_map(tdim);
261 assert(cell_map);
262 std::int32_t num_cells = cell_map->size_local();
263 auto _coordinate_dofs
264 = xt::view(coordinate_dofs, xt::all(), xt::xrange(0, gdim));
265 for (std::int32_t c = 0; c < num_cells; ++c)
266 {
267 // Get cell geometry (coordinate dofs)
268 auto x_dofs = x_dofmap.links(c);
269 for (std::size_t i = 0; i < x_dofs.size(); ++i)
270 {
271 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
272 std::next(coordinate_dofs.begin(), 3 * i));
273 }
274
275 // Compute Jacobians and reference points for current cell
276 J.fill(0);
277 for (std::size_t p = 0; p < X.shape(0); ++p)
278 {
279 auto _J = xt::view(J, p, xt::all(), xt::all());
280 cmap.compute_jacobian(dphi, _coordinate_dofs, _J);
281 cmap.compute_jacobian_inverse(_J, xt::view(K, p, xt::all(), xt::all()));
282 detJ[p] = cmap.compute_jacobian_determinant(_J);
283 }
284
285 // Get evaluated basis on reference, apply DOF transformations, and
286 // push forward to physical element
287 basis_reference0 = xt::view(basis_derivatives_reference0, 0, xt::all(),
288 xt::all(), xt::all());
289 for (std::size_t p = 0; p < X.shape(0); ++p)
290 {
291 apply_dof_transformation0(
292 xtl::span(basis_reference0.data() + p * dim0 * value_size_ref0,
293 dim0 * value_size_ref0),
294 cell_info, c, value_size_ref0);
295 }
296
297 for (std::size_t i = 0; i < basis0.shape(0); ++i)
298 {
299 auto _K = xt::view(K, i, xt::all(), xt::all());
300 auto _J = xt::view(J, i, xt::all(), xt::all());
301 auto _u = xt::view(basis0, i, xt::all(), xt::all());
302 auto _U = xt::view(basis_reference0, i, xt::all(), xt::all());
303 push_forward_fn0(_u, _U, _J, detJ[i], _K);
304 }
305
306 // Unroll basis function for input space for block size
307 for (std::size_t p = 0; p < X.shape(0); ++p)
308 for (std::size_t i = 0; i < dim0; ++i)
309 for (std::size_t j = 0; j < value_size0; ++j)
310 for (int k = 0; k < bs0; ++k)
311 basis_values(p, i * bs0 + k, j * bs0 + k) = basis0(p, i, j);
312
313 // Pull back the physical values to the reference of output space
314 for (std::size_t p = 0; p < basis_values.shape(0); ++p)
315 {
316 auto _K = xt::view(K, p, xt::all(), xt::all());
317 auto _J = xt::view(J, p, xt::all(), xt::all());
318 auto _u = xt::view(basis_values, p, xt::all(), xt::all());
319 auto _U = xt::view(mapped_values, p, xt::all(), xt::all());
320 pull_back_fn1(_U, _u, _K, 1.0 / detJ[p], _J);
321 }
322
323 // Apply interpolation matrix to basis values of V0 at the interpolation
324 // points of V1
325 if (interpolation_ident)
326 _A.assign(xt::transpose(mapped_values, {0, 2, 1}));
327 else
328 {
329 for (std::size_t i = 0; i < mapped_values.shape(1); ++i)
330 {
331 auto _mapped_values = xt::view(mapped_values, xt::all(), i, xt::all());
332 impl::interpolation_apply(Pi_1, _mapped_values, local1, bs1);
333 for (std::size_t j = 0; j < local1.size(); j++)
334 A[space_dim0 * j + i] = local1[j];
335 }
336 }
337 apply_inverse_dof_transform1(A, cell_info, c, space_dim0);
338 mat_set(dofmap1->cell_dofs(c), dofmap0->cell_dofs(c), A);
339 }
340}
341
342} // namespace dolfinx::fem
Degree-of-freedeom map representations ans tools.
Definition: CoordinateElement.h:31
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:107
static double compute_jacobian_determinant(const U &J)
Compute the determinant of the Jacobian.
Definition: CoordinateElement.h:130
static void compute_jacobian_inverse(const U &J, V &&K)
Compute the inverse of the Jacobian.
Definition: CoordinateElement.h:116
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:44
int dim() const
The dimension of the geometry element space.
Definition: CoordinateElement.cpp:207
xt::xtensor< double, 4 > tabulate(int nd, const xt::xtensor< double, 2 > &X) const
Evaluate basis values and derivatives at set of points.
Definition: CoordinateElement.cpp:50
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:32
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition: FunctionSpace.cpp:248
std::shared_ptr< const mesh::Mesh > mesh() const
The mesh.
Definition: FunctionSpace.cpp:241
std::shared_ptr< const FiniteElement > element() const
The finite element.
Definition: FunctionSpace.cpp:243
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
void interpolation_matrix(const fem::FunctionSpace &V0, const fem::FunctionSpace &V1, U &&mat_set)
Assemble an interpolation operator matrix.
Definition: discreteoperators.h:146
void discrete_gradient(const fem::FunctionSpace &V0, const fem::FunctionSpace &V1, U &&mat_set)
Assemble a discrete gradient operator.
Definition: discreteoperators.h:52