DOLFINx
DOLFINx C++ interface
Loading...
Searching...
No Matches
FunctionSpace.h
1// Copyright (C) 2008-2023 Anders Logg and Garth N. Wells
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 "CoordinateElement.h"
10#include "DofMap.h"
11#include "FiniteElement.h"
12#include <boost/uuid/uuid.hpp>
13#include <boost/uuid/uuid_generators.hpp>
14#include <concepts>
15#include <cstddef>
16#include <cstdint>
17#include <dolfinx/common/IndexMap.h>
18#include <dolfinx/mesh/Geometry.h>
19#include <dolfinx/mesh/Mesh.h>
20#include <dolfinx/mesh/Topology.h>
21#include <map>
22#include <memory>
23#include <vector>
24
25namespace dolfinx::fem
26{
32template <std::floating_point T>
34{
35public:
37 using geometry_type = T;
38
45 std::shared_ptr<const FiniteElement<geometry_type>> element,
46 std::shared_ptr<const DofMap> dofmap,
47 std::vector<std::size_t> value_shape)
48 : _mesh(mesh), _element(element), _dofmap(dofmap),
49 _id(boost::uuids::random_generator()()), _root_space_id(_id),
50 _value_shape(value_shape)
51 {
52 // Do nothing
53 }
54
55 // Copy constructor (deleted)
56 FunctionSpace(const FunctionSpace& V) = delete;
57
60
62 virtual ~FunctionSpace() = default;
63
64 // Assignment operator (delete)
65 FunctionSpace& operator=(const FunctionSpace& V) = delete;
66
69
78 FunctionSpace sub(const std::vector<int>& component) const
79 {
80 assert(_mesh);
81 assert(_element);
82 assert(_dofmap);
83
84 // Check that component is valid
85 if (component.empty())
86 throw std::runtime_error("Component must be non-empty");
87
88 // Extract sub-element
89 auto element = this->_element->extract_sub_element(component);
90
91 // Extract sub dofmap
92 auto dofmap
93 = std::make_shared<DofMap>(_dofmap->extract_sub_dofmap(component));
94
95 // Create new sub space
96 FunctionSpace sub_space(_mesh, element, dofmap,
98 _mesh->topology()->dim(),
99 _mesh->geometry().dim()));
100
101 // Set root space id and component w.r.t. root
102 sub_space._root_space_id = _root_space_id;
103 sub_space._component = _component;
104 sub_space._component.insert(sub_space._component.end(), component.begin(),
105 component.end());
106 return sub_space;
107 }
108
113 bool contains(const FunctionSpace& V) const
114 {
115 if (this == std::addressof(V))
116 {
117 // Spaces are the same (same memory address)
118 return true;
119 }
120 else if (_root_space_id != V._root_space_id)
121 {
122 // Different root spaces
123 return false;
124 }
125 else if (_component.size() > V._component.size())
126 {
127 // V is a superspace of *this
128 return false;
129 }
130 else if (!std::equal(_component.begin(), _component.end(),
131 V._component.begin()))
132 {
133 // Components of 'this' are not the same as the leading components
134 // of V
135 return false;
136 }
137 else
138 {
139 // Ok, V is really our subspace
140 return true;
141 }
142 }
143
147 std::pair<FunctionSpace, std::vector<std::int32_t>> collapse() const
148 {
149 if (_component.empty())
150 throw std::runtime_error("Function space is not a subspace");
151
152 // Create collapsed DofMap
153 auto [_collapsed_dofmap, collapsed_dofs]
154 = _dofmap->collapse(_mesh->comm(), *_mesh->topology());
155 auto collapsed_dofmap
156 = std::make_shared<DofMap>(std::move(_collapsed_dofmap));
157
158 return {
159 FunctionSpace(_mesh, _element, collapsed_dofmap,
160 compute_value_shape(_element, _mesh->topology()->dim(),
161 _mesh->geometry().dim())),
162 std::move(collapsed_dofs)};
163 }
164
168 std::vector<int> component() const { return _component; }
169
181 std::vector<geometry_type> tabulate_dof_coordinates(bool transpose) const
182 {
183 if (!_component.empty())
184 {
185 throw std::runtime_error("Cannot tabulate coordinates for a "
186 "FunctionSpace that is a subspace.");
187 }
188
189 assert(_element);
190 if (_element->is_mixed())
191 {
192 throw std::runtime_error(
193 "Cannot tabulate coordinates for a mixed FunctionSpace.");
194 }
195
196 // Geometric dimension
197 assert(_mesh);
198 assert(_element);
199 const std::size_t gdim = _mesh->geometry().dim();
200 const int tdim = _mesh->topology()->dim();
201
202 // Get dofmap local size
203 assert(_dofmap);
204 std::shared_ptr<const common::IndexMap> index_map = _dofmap->index_map;
205 assert(index_map);
206 const int index_map_bs = _dofmap->index_map_bs();
207 const int dofmap_bs = _dofmap->bs();
208
209 const int element_block_size = _element->block_size();
210 const std::size_t scalar_dofs
211 = _element->space_dimension() / element_block_size;
212 const std::int32_t num_dofs
213 = index_map_bs * (index_map->size_local() + index_map->num_ghosts())
214 / dofmap_bs;
215
216 // Get the dof coordinates on the reference element
217 if (!_element->interpolation_ident())
218 {
219 throw std::runtime_error("Cannot evaluate dof coordinates - this element "
220 "does not have pointwise evaluation.");
221 }
222 const auto [X, Xshape] = _element->interpolation_points();
223
224 // Get coordinate map
225 const CoordinateElement<geometry_type>& cmap = _mesh->geometry().cmap();
226
227 // Prepare cell geometry
228 auto x_dofmap = _mesh->geometry().dofmap();
229 const std::size_t num_dofs_g = cmap.dim();
230 std::span<const geometry_type> x_g = _mesh->geometry().x();
231
232 // Array to hold coordinates to return
233 const std::size_t shape_c0 = transpose ? 3 : num_dofs;
234 const std::size_t shape_c1 = transpose ? num_dofs : 3;
235 std::vector<geometry_type> coords(shape_c0 * shape_c1, 0);
236
237 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
239 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
240 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
241 const geometry_type,
242 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
243
244 // Loop over cells and tabulate dofs
245 std::vector<geometry_type> x_b(scalar_dofs * gdim);
246 mdspan2_t x(x_b.data(), scalar_dofs, gdim);
247
248 std::vector<geometry_type> coordinate_dofs_b(num_dofs_g * gdim);
249 mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
250
251 auto map = _mesh->topology()->index_map(tdim);
252 assert(map);
253 const int num_cells = map->size_local() + map->num_ghosts();
254
255 std::span<const std::uint32_t> cell_info;
256 if (_element->needs_dof_transformations())
257 {
258 _mesh->topology_mutable()->create_entity_permutations();
259 cell_info = std::span(_mesh->topology()->get_cell_permutation_info());
260 }
261
262 auto apply_dof_transformation
263 = _element
264 ->template get_pre_dof_transformation_function<geometry_type>();
265
266 const std::array<std::size_t, 4> phi_shape
267 = cmap.tabulate_shape(0, Xshape[0]);
268 std::vector<geometry_type> phi_b(
269 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
270 cmdspan4_t phi_full(phi_b.data(), phi_shape);
271 cmap.tabulate(0, X, Xshape, phi_b);
272 auto phi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
273 phi_full, 0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
274 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
275
276 for (int c = 0; c < num_cells; ++c)
277 {
278 // Extract cell geometry
279 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
280 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
281 for (std::size_t i = 0; i < x_dofs.size(); ++i)
282 for (std::size_t j = 0; j < gdim; ++j)
283 coordinate_dofs(i, j) = x_g[3 * x_dofs[i] + j];
284
285 // Tabulate dof coordinates on cell
286 cmap.push_forward(x, coordinate_dofs, phi);
287 apply_dof_transformation(
288 x_b, std::span(cell_info.data(), cell_info.size()), c, x.extent(1));
289
290 // Get cell dofmap
291 auto dofs = _dofmap->cell_dofs(c);
292
293 // Copy dof coordinates into vector
294 if (!transpose)
295 {
296 for (std::size_t i = 0; i < dofs.size(); ++i)
297 for (std::size_t j = 0; j < gdim; ++j)
298 coords[dofs[i] * 3 + j] = x(i, j);
299 }
300 else
301 {
302 for (std::size_t i = 0; i < dofs.size(); ++i)
303 for (std::size_t j = 0; j < gdim; ++j)
304 coords[j * num_dofs + dofs[i]] = x(i, j);
305 }
306 }
307
308 return coords;
309 }
310
312 std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const
313 {
314 return _mesh;
315 }
316
318 std::shared_ptr<const FiniteElement<geometry_type>> element() const
319 {
320 return _element;
321 }
322
324 std::shared_ptr<const DofMap> dofmap() const { return _dofmap; }
325
327 std::span<const std::size_t> value_shape() const noexcept
328 {
329 return _value_shape;
330 }
331
337 int value_size() const
338 {
339 return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
340 std::multiplies{});
341 }
342
343private:
344 // The mesh
345 std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
346
347 // The finite element
348 std::shared_ptr<const FiniteElement<geometry_type>> _element;
349
350 // The dofmap
351 std::shared_ptr<const DofMap> _dofmap;
352
353 // The component w.r.t. to root space
354 std::vector<int> _component;
355
356 // Unique identifier for the space and for its root space
357 boost::uuids::uuid _id;
358 boost::uuids::uuid _root_space_id;
359
360 std::vector<std::size_t> _value_shape;
361};
362
372template <dolfinx::scalar T>
373std::array<std::vector<std::shared_ptr<const FunctionSpace<T>>>, 2>
375 const std::vector<
376 std::vector<std::array<std::shared_ptr<const FunctionSpace<T>>, 2>>>& V)
377{
378 assert(!V.empty());
379 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces0(V.size(),
380 nullptr);
381 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces1(V.front().size(),
382 nullptr);
383
384 // Loop over rows
385 for (std::size_t i = 0; i < V.size(); ++i)
386 {
387 // Loop over columns
388 for (std::size_t j = 0; j < V[i].size(); ++j)
389 {
390 auto& V0 = V[i][j][0];
391 auto& V1 = V[i][j][1];
392 if (V0 and V1)
393 {
394 if (!spaces0[i])
395 spaces0[i] = V0;
396 else
397 {
398 if (spaces0[i] != V0)
399 throw std::runtime_error("Mismatched test space for row.");
400 }
401
402 if (!spaces1[j])
403 spaces1[j] = V1;
404 else
405 {
406 if (spaces1[j] != V1)
407 throw std::runtime_error("Mismatched trial space for column.");
408 }
409 }
410 }
411 }
412
413 // Check that there are no null entries
414 if (std::find(spaces0.begin(), spaces0.end(), nullptr) != spaces0.end())
415 throw std::runtime_error("Could not deduce all block test spaces.");
416 if (std::find(spaces1.begin(), spaces1.end(), nullptr) != spaces1.end())
417 throw std::runtime_error("Could not deduce all block trial spaces.");
418
419 return {spaces0, spaces1};
420}
421
427template <std::floating_point T>
428std::vector<std::size_t> compute_value_shape(
429 std::shared_ptr<const dolfinx::fem::FiniteElement<T>> element,
430 std::size_t tdim, std::size_t gdim)
431{
432 auto rvs = element->reference_value_shape();
433 std::vector<std::size_t> value_shape(rvs.size());
434 if (element->block_size() > 1)
435 {
436 for (std::size_t i = 0; i < rvs.size(); ++i)
437 {
438 value_shape[i] = rvs[i];
439 }
440 }
441 else
442 {
443 for (std::size_t i = 0; i < rvs.size(); ++i)
444 {
445 if (rvs[i] == tdim)
446 value_shape[i] = gdim;
447 else
448 value_shape[i] = rvs[i];
449 }
450 }
451 return value_shape;
452}
453
455template <typename U, typename V, typename W, typename X>
456FunctionSpace(U mesh, V element, W dofmap, X value_shape)
457 -> FunctionSpace<typename std::remove_cvref<
458 typename U::element_type>::type::geometry_type::value_type>;
459
460} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition CoordinateElement.h:38
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:53
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:46
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:193
static void push_forward(U &&x, const V &cell_geometry, const W &phi)
Compute physical coordinates x for points X in the reference configuration.
Definition CoordinateElement.h:167
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
std::vector< int > component() const
Get the component with respect to the root superspace.
Definition FunctionSpace.h:168
FunctionSpace & operator=(FunctionSpace &&V)=default
Move assignment operator.
std::pair< FunctionSpace, std::vector< std::int32_t > > collapse() const
Collapse a subspace and return a new function space and a map from new to old dofs.
Definition FunctionSpace.h:147
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:324
FunctionSpace(std::shared_ptr< const mesh::Mesh< geometry_type > > mesh, std::shared_ptr< const FiniteElement< geometry_type > > element, std::shared_ptr< const DofMap > dofmap, std::vector< std::size_t > value_shape)
Create function space for given mesh, element and dofmap.
Definition FunctionSpace.h:44
FunctionSpace sub(const std::vector< int > &component) const
Create a subspace (view) for a specific component.
Definition FunctionSpace.h:78
bool contains(const FunctionSpace &V) const
Check whether V is subspace of this, or this itself.
Definition FunctionSpace.h:113
int value_size() const
The value size, e.g. 1 for a scalar-valued function, 2 for a 2D vector, 9 for a second-order tensor i...
Definition FunctionSpace.h:337
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:312
virtual ~FunctionSpace()=default
Destructor.
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:318
T geometry_type
Geometry type of the Mesh that the FunctionSpace is defined on.
Definition FunctionSpace.h:37
std::span< const std::size_t > value_shape() const noexcept
The shape of the value space.
Definition FunctionSpace.h:327
std::vector< geometry_type > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:181
FunctionSpace(FunctionSpace &&V)=default
Move constructor.
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Finite element method functionality.
Definition assemble_matrix_impl.h:26
std::vector< std::size_t > compute_value_shape(std::shared_ptr< const dolfinx::fem::FiniteElement< T > > element, std::size_t tdim, std::size_t gdim)
Compute the physical value shape of an element for a mesh.
Definition FunctionSpace.h:428
std::array< std::vector< std::shared_ptr< const FunctionSpace< T > > >, 2 > common_function_spaces(const std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< T > >, 2 > > > &V)
Extract FunctionSpaces for (0) rows blocks and (1) columns blocks from a rectangular array of (test,...
Definition FunctionSpace.h:374