DOLFINx 0.9.0
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>
33class FunctionSpace
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,
97 compute_value_shape(element,
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
172 bool symmetric() const
173 {
174 if (_element)
175 return _element->symmetric();
176 return false;
177 }
178
190 std::vector<geometry_type> tabulate_dof_coordinates(bool transpose) const
191 {
192 if (!_component.empty())
193 {
194 throw std::runtime_error("Cannot tabulate coordinates for a "
195 "FunctionSpace that is a subspace.");
196 }
197
198 assert(_element);
199 if (_element->is_mixed())
200 {
201 throw std::runtime_error(
202 "Cannot tabulate coordinates for a mixed FunctionSpace.");
203 }
204
205 // Geometric dimension
206 assert(_mesh);
207 assert(_element);
208 const std::size_t gdim = _mesh->geometry().dim();
209 const int tdim = _mesh->topology()->dim();
210
211 // Get dofmap local size
212 assert(_dofmap);
213 std::shared_ptr<const common::IndexMap> index_map = _dofmap->index_map;
214 assert(index_map);
215 const int index_map_bs = _dofmap->index_map_bs();
216 const int dofmap_bs = _dofmap->bs();
217
218 const int element_block_size = _element->block_size();
219 const std::size_t scalar_dofs
220 = _element->space_dimension() / element_block_size;
221 const std::int32_t num_dofs
222 = index_map_bs * (index_map->size_local() + index_map->num_ghosts())
223 / dofmap_bs;
224
225 // Get the dof coordinates on the reference element
226 if (!_element->interpolation_ident())
227 {
228 throw std::runtime_error("Cannot evaluate dof coordinates - this element "
229 "does not have pointwise evaluation.");
230 }
231 const auto [X, Xshape] = _element->interpolation_points();
232
233 // Get coordinate map
234 const CoordinateElement<geometry_type>& cmap = _mesh->geometry().cmap();
235
236 // Prepare cell geometry
237 auto x_dofmap = _mesh->geometry().dofmap();
238 const std::size_t num_dofs_g = cmap.dim();
239 std::span<const geometry_type> x_g = _mesh->geometry().x();
240
241 // Array to hold coordinates to return
242 const std::size_t shape_c0 = transpose ? 3 : num_dofs;
243 const std::size_t shape_c1 = transpose ? num_dofs : 3;
244 std::vector<geometry_type> coords(shape_c0 * shape_c1, 0);
245
246 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
248 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
249 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
250 const geometry_type,
251 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
252
253 // Loop over cells and tabulate dofs
254 std::vector<geometry_type> x_b(scalar_dofs * gdim);
255 mdspan2_t x(x_b.data(), scalar_dofs, gdim);
256
257 std::vector<geometry_type> coordinate_dofs_b(num_dofs_g * gdim);
258 mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
259
260 auto map = _mesh->topology()->index_map(tdim);
261 assert(map);
262 const int num_cells = map->size_local() + map->num_ghosts();
263
264 std::span<const std::uint32_t> cell_info;
265 if (_element->needs_dof_transformations())
266 {
267 _mesh->topology_mutable()->create_entity_permutations();
268 cell_info = std::span(_mesh->topology()->get_cell_permutation_info());
269 }
270
271 const std::array<std::size_t, 4> phi_shape
272 = cmap.tabulate_shape(0, Xshape[0]);
273 std::vector<geometry_type> phi_b(
274 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
275 cmdspan4_t phi_full(phi_b.data(), phi_shape);
276 cmap.tabulate(0, X, Xshape, phi_b);
277 auto phi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
278 phi_full, 0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
279 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
280
281 // TODO: Check transform
282 // Basis function reference-to-conforming transformation function
283 auto apply_dof_transformation
284 = _element->template dof_transformation_fn<geometry_type>(
286
287 for (int c = 0; c < num_cells; ++c)
288 {
289 // Extract cell geometry 'dofs'
290 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
291 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
292 for (std::size_t i = 0; i < x_dofs.size(); ++i)
293 for (std::size_t j = 0; j < gdim; ++j)
294 coordinate_dofs(i, j) = x_g[3 * x_dofs[i] + j];
295
296 // Tabulate dof coordinates on cell
297 cmap.push_forward(x, coordinate_dofs, phi);
298 apply_dof_transformation(
299 x_b, std::span(cell_info.data(), cell_info.size()), c, x.extent(1));
300
301 // Get cell dofmap
302 auto dofs = _dofmap->cell_dofs(c);
303
304 // Copy dof coordinates into vector
305 if (!transpose)
306 {
307 for (std::size_t i = 0; i < dofs.size(); ++i)
308 for (std::size_t j = 0; j < gdim; ++j)
309 coords[dofs[i] * 3 + j] = x(i, j);
310 }
311 else
312 {
313 for (std::size_t i = 0; i < dofs.size(); ++i)
314 for (std::size_t j = 0; j < gdim; ++j)
315 coords[j * num_dofs + dofs[i]] = x(i, j);
316 }
317 }
318
319 return coords;
320 }
321
323 std::shared_ptr<const mesh::Mesh<geometry_type>> mesh() const
324 {
325 return _mesh;
326 }
327
329 std::shared_ptr<const FiniteElement<geometry_type>> element() const
330 {
331 return _element;
332 }
333
335 std::shared_ptr<const DofMap> dofmap() const { return _dofmap; }
336
338 std::span<const std::size_t> value_shape() const noexcept
339 {
340 return _value_shape;
341 }
342
348 int value_size() const
349 {
350 return std::accumulate(_value_shape.begin(), _value_shape.end(), 1,
351 std::multiplies{});
352 }
353
354private:
355 // The mesh
356 std::shared_ptr<const mesh::Mesh<geometry_type>> _mesh;
357
358 // The finite element
359 std::shared_ptr<const FiniteElement<geometry_type>> _element;
360
361 // The dofmap
362 std::shared_ptr<const DofMap> _dofmap;
363
364 // The component w.r.t. to root space
365 std::vector<int> _component;
366
367 // Unique identifier for the space and for its root space
368 boost::uuids::uuid _id;
369 boost::uuids::uuid _root_space_id;
370
371 std::vector<std::size_t> _value_shape;
372};
373
383template <dolfinx::scalar T>
384std::array<std::vector<std::shared_ptr<const FunctionSpace<T>>>, 2>
386 const std::vector<
387 std::vector<std::array<std::shared_ptr<const FunctionSpace<T>>, 2>>>& V)
388{
389 assert(!V.empty());
390 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces0(V.size(),
391 nullptr);
392 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces1(V.front().size(),
393 nullptr);
394
395 // Loop over rows
396 for (std::size_t i = 0; i < V.size(); ++i)
397 {
398 // Loop over columns
399 for (std::size_t j = 0; j < V[i].size(); ++j)
400 {
401 auto& V0 = V[i][j][0];
402 auto& V1 = V[i][j][1];
403 if (V0 and V1)
404 {
405 if (!spaces0[i])
406 spaces0[i] = V0;
407 else
408 {
409 if (spaces0[i] != V0)
410 throw std::runtime_error("Mismatched test space for row.");
411 }
412
413 if (!spaces1[j])
414 spaces1[j] = V1;
415 else
416 {
417 if (spaces1[j] != V1)
418 throw std::runtime_error("Mismatched trial space for column.");
419 }
420 }
421 }
422 }
423
424 // Check that there are no null entries
425 if (std::find(spaces0.begin(), spaces0.end(), nullptr) != spaces0.end())
426 throw std::runtime_error("Could not deduce all block test spaces.");
427 if (std::find(spaces1.begin(), spaces1.end(), nullptr) != spaces1.end())
428 throw std::runtime_error("Could not deduce all block trial spaces.");
429
430 return {spaces0, spaces1};
431}
432
438template <std::floating_point T>
439std::vector<std::size_t> compute_value_shape(
440 std::shared_ptr<const dolfinx::fem::FiniteElement<T>> element,
441 std::size_t tdim, std::size_t gdim)
442{
443 auto rvs = element->reference_value_shape();
444 std::vector<std::size_t> value_shape(rvs.size());
445 if (element->block_size() > 1)
446 {
447 for (std::size_t i = 0; i < rvs.size(); ++i)
448 {
449 value_shape[i] = rvs[i];
450 }
451 }
452 else
453 {
454 for (std::size_t i = 0; i < rvs.size(); ++i)
455 {
456 if (rvs[i] == tdim)
457 value_shape[i] = gdim;
458 else
459 value_shape[i] = rvs[i];
460 }
461 }
462 return value_shape;
463}
464
466template <typename U, typename V, typename W, typename X>
467FunctionSpace(U mesh, V element, W dofmap, X value_shape)
468 -> FunctionSpace<typename std::remove_cvref<
469 typename U::element_type>::type::geometry_type::value_type>;
470
471} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
Definition XDMFFile.h:27
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
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:206
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:188
Model of a finite element.
Definition FiniteElement.h:38
This class represents a finite element function space defined by a mesh, a finite element,...
Definition vtk_utils.h:32
std::vector< int > component() const
Get the component with respect to the root superspace.
Definition FunctionSpace.h:168
bool symmetric() const
Indicate whether this function space represents a symmetric 2-tensor.
Definition FunctionSpace.h:172
FunctionSpace & operator=(FunctionSpace &&V)=default
Move assignment operator.
std::pair< FunctionSpace, std::vector< std::int32_t > > collapse() const
Definition FunctionSpace.h:147
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:335
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
Definition FunctionSpace.h:348
std::shared_ptr< const mesh::Mesh< geometry_type > > mesh() const
The mesh.
Definition FunctionSpace.h:323
virtual ~FunctionSpace()=default
Destructor.
std::shared_ptr< const FiniteElement< geometry_type > > element() const
The finite element.
Definition FunctionSpace.h:329
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:338
std::vector< geometry_type > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:190
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
FunctionSpace(U mesh, V element, W dofmap, X value_shape) -> FunctionSpace< typename std::remove_cvref< typename U::element_type >::type::geometry_type::value_type >
Type deduction.
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:439
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)
Definition FunctionSpace.h:385