Note: this is documentation for an old release. View the latest documentation at docs.fenicsproject.org/dolfinx/v0.9.0/cpp/doxygen/d6/df3/FunctionSpace_8h_source.html
DOLFINx 0.7.3
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{
27
33template <std::floating_point T>
35{
36public:
41 FunctionSpace(std::shared_ptr<const mesh::Mesh<T>> mesh,
42 std::shared_ptr<const FiniteElement<T>> element,
43 std::shared_ptr<const DofMap> dofmap)
44 : _mesh(mesh), _element(element), _dofmap(dofmap),
45 _id(boost::uuids::random_generator()()), _root_space_id(_id)
46 {
47 // Do nothing
48 }
49
50 // Copy constructor (deleted)
51 FunctionSpace(const FunctionSpace& V) = delete;
52
55
57 virtual ~FunctionSpace() = default;
58
59 // Assignment operator (delete)
60 FunctionSpace& operator=(const FunctionSpace& V) = delete;
61
64
68 std::shared_ptr<FunctionSpace<T>> sub(const std::vector<int>& component) const
69 {
70 assert(_mesh);
71 assert(_element);
72 assert(_dofmap);
73
74 // Check if sub space is already in the cache and not expired
75 if (auto it = _subspaces.find(component); it != _subspaces.end())
76 {
77 if (auto s = it->second.lock())
78 return s;
79 }
80
81 // Extract sub-element
82 auto element = this->_element->extract_sub_element(component);
83
84 // Extract sub dofmap
85 auto dofmap
86 = std::make_shared<DofMap>(_dofmap->extract_sub_dofmap(component));
87
88 // Create new sub space
89 auto sub_space = std::make_shared<FunctionSpace<T>>(_mesh, element, dofmap);
90
91 // Set root space id and component w.r.t. root
92 sub_space->_root_space_id = _root_space_id;
93 sub_space->_component = _component;
94 sub_space->_component.insert(sub_space->_component.end(), component.begin(),
95 component.end());
96
97 // Insert new subspace into cache
98 _subspaces.emplace(sub_space->_component, sub_space);
99
100 return sub_space;
101 }
102
107 bool contains(const FunctionSpace& V) const
108 {
109 if (this == std::addressof(V))
110 {
111 // Spaces are the same (same memory address)
112 return true;
113 }
114 else if (_root_space_id != V._root_space_id)
115 {
116 // Different root spaces
117 return false;
118 }
119 else if (_component.size() > V._component.size())
120 {
121 // V is a superspace of *this
122 return false;
123 }
124 else if (!std::equal(_component.begin(), _component.end(),
125 V._component.begin()))
126 {
127 // Components of 'this' are not the same as the leading components
128 // of V
129 return false;
130 }
131 else
132 {
133 // Ok, V is really our subspace
134 return true;
135 }
136 }
137
141 std::pair<FunctionSpace<T>, std::vector<std::int32_t>> collapse() const
142 {
143 if (_component.empty())
144 throw std::runtime_error("Function space is not a subspace");
145
146 // Create collapsed DofMap
147 auto [_collapsed_dofmap, collapsed_dofs]
148 = _dofmap->collapse(_mesh->comm(), *_mesh->topology());
149 auto collapsed_dofmap
150 = std::make_shared<DofMap>(std::move(_collapsed_dofmap));
151
152 return {FunctionSpace(_mesh, _element, collapsed_dofmap),
153 std::move(collapsed_dofs)};
154 }
155
159 std::vector<int> component() const { return _component; }
160
172 std::vector<T> tabulate_dof_coordinates(bool transpose) const
173 {
174 if (!_component.empty())
175 {
176 throw std::runtime_error("Cannot tabulate coordinates for a "
177 "FunctionSpace that is a subspace.");
178 }
179
180 assert(_element);
181 if (_element->is_mixed())
182 {
183 throw std::runtime_error(
184 "Cannot tabulate coordinates for a mixed FunctionSpace.");
185 }
186
187 // Geometric dimension
188 assert(_mesh);
189 assert(_element);
190 const std::size_t gdim = _mesh->geometry().dim();
191 const int tdim = _mesh->topology()->dim();
192
193 // Get dofmap local size
194 assert(_dofmap);
195 std::shared_ptr<const common::IndexMap> index_map = _dofmap->index_map;
196 assert(index_map);
197 const int index_map_bs = _dofmap->index_map_bs();
198 const int dofmap_bs = _dofmap->bs();
199
200 const int element_block_size = _element->block_size();
201 const std::size_t scalar_dofs
202 = _element->space_dimension() / element_block_size;
203 const std::int32_t num_dofs
204 = index_map_bs * (index_map->size_local() + index_map->num_ghosts())
205 / dofmap_bs;
206
207 // Get the dof coordinates on the reference element
208 if (!_element->interpolation_ident())
209 {
210 throw std::runtime_error("Cannot evaluate dof coordinates - this element "
211 "does not have pointwise evaluation.");
212 }
213 const auto [X, Xshape] = _element->interpolation_points();
214
215 // Get coordinate map
216 if (_mesh->geometry().cmaps().size() > 1)
217 throw std::runtime_error("Mixed topology not supported");
218 const CoordinateElement<T>& cmap = _mesh->geometry().cmaps()[0];
219
220 // Prepare cell geometry
221 auto x_dofmap = _mesh->geometry().dofmap();
222 const std::size_t num_dofs_g = cmap.dim();
223 std::span<const T> x_g = _mesh->geometry().x();
224
225 // Array to hold coordinates to return
226 const std::size_t shape_c0 = transpose ? 3 : num_dofs;
227 const std::size_t shape_c1 = transpose ? num_dofs : 3;
228 std::vector<T> coords(shape_c0 * shape_c1, 0);
229
230 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
231 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
232 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
233 const T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
234
235 // Loop over cells and tabulate dofs
236 std::vector<T> x_b(scalar_dofs * gdim);
237 mdspan2_t x(x_b.data(), scalar_dofs, gdim);
238
239 std::vector<T> coordinate_dofs_b(num_dofs_g * gdim);
240 mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
241
242 auto map = _mesh->topology()->index_map(tdim);
243 assert(map);
244 const int num_cells = map->size_local() + map->num_ghosts();
245
246 std::span<const std::uint32_t> cell_info;
247 if (_element->needs_dof_transformations())
248 {
249 _mesh->topology_mutable()->create_entity_permutations();
250 cell_info = std::span(_mesh->topology()->get_cell_permutation_info());
251 }
252
253 auto apply_dof_transformation
254 = _element->template get_dof_transformation_function<T>();
255
256 const std::array<std::size_t, 4> phi_shape
257 = cmap.tabulate_shape(0, Xshape[0]);
258 std::vector<T> phi_b(
259 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
260 cmdspan4_t phi_full(phi_b.data(), phi_shape);
261 cmap.tabulate(0, X, Xshape, phi_b);
262 auto phi = MDSPAN_IMPL_STANDARD_NAMESPACE::MDSPAN_IMPL_PROPOSED_NAMESPACE::
263 submdspan(phi_full, 0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
264 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
265
266 for (int c = 0; c < num_cells; ++c)
267 {
268 // Extract cell geometry
269 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::
270 MDSPAN_IMPL_PROPOSED_NAMESPACE::submdspan(
271 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
272 for (std::size_t i = 0; i < x_dofs.size(); ++i)
273 for (std::size_t j = 0; j < gdim; ++j)
274 coordinate_dofs(i, j) = x_g[3 * x_dofs[i] + j];
275
276 // Tabulate dof coordinates on cell
277 cmap.push_forward(x, coordinate_dofs, phi);
278 apply_dof_transformation(
279 x_b, std::span(cell_info.data(), cell_info.size()), c, x.extent(1));
280
281 // Get cell dofmap
282 auto dofs = _dofmap->cell_dofs(c);
283
284 // Copy dof coordinates into vector
285 if (!transpose)
286 {
287 for (std::size_t i = 0; i < dofs.size(); ++i)
288 for (std::size_t j = 0; j < gdim; ++j)
289 coords[dofs[i] * 3 + j] = x(i, j);
290 }
291 else
292 {
293 for (std::size_t i = 0; i < dofs.size(); ++i)
294 for (std::size_t j = 0; j < gdim; ++j)
295 coords[j * num_dofs + dofs[i]] = x(i, j);
296 }
297 }
298
299 return coords;
300 }
301
303 std::shared_ptr<const mesh::Mesh<T>> mesh() const { return _mesh; }
304
306 std::shared_ptr<const FiniteElement<T>> element() const { return _element; }
307
309 std::shared_ptr<const DofMap> dofmap() const { return _dofmap; }
310
311private:
312 // The mesh
313 std::shared_ptr<const mesh::Mesh<T>> _mesh;
314
315 // The finite element
316 std::shared_ptr<const FiniteElement<T>> _element;
317
318 // The dofmap
319 std::shared_ptr<const DofMap> _dofmap;
320
321 // The component w.r.t. to root space
322 std::vector<int> _component;
323
324 // Unique identifier for the space and for its root space
325 boost::uuids::uuid _id;
326 boost::uuids::uuid _root_space_id;
327
328 // Cache of subspaces
329 mutable std::map<std::vector<int>, std::weak_ptr<FunctionSpace>> _subspaces;
330};
331
341template <dolfinx::scalar T>
342std::array<std::vector<std::shared_ptr<const FunctionSpace<T>>>, 2>
344 const std::vector<
345 std::vector<std::array<std::shared_ptr<const FunctionSpace<T>>, 2>>>& V)
346{
347 assert(!V.empty());
348 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces0(V.size(),
349 nullptr);
350 std::vector<std::shared_ptr<const FunctionSpace<T>>> spaces1(V.front().size(),
351 nullptr);
352
353 // Loop over rows
354 for (std::size_t i = 0; i < V.size(); ++i)
355 {
356 // Loop over columns
357 for (std::size_t j = 0; j < V[i].size(); ++j)
358 {
359 auto& V0 = V[i][j][0];
360 auto& V1 = V[i][j][1];
361 if (V0 and V1)
362 {
363 if (!spaces0[i])
364 spaces0[i] = V0;
365 else
366 {
367 if (spaces0[i] != V0)
368 throw std::runtime_error("Mismatched test space for row.");
369 }
370
371 if (!spaces1[j])
372 spaces1[j] = V1;
373 else
374 {
375 if (spaces1[j] != V1)
376 throw std::runtime_error("Mismatched trial space for column.");
377 }
378 }
379 }
380 }
381
382 // Check that there are no null entries
383 if (std::find(spaces0.begin(), spaces0.end(), nullptr) != spaces0.end())
384 throw std::runtime_error("Could not deduce all block test spaces.");
385 if (std::find(spaces1.begin(), spaces1.end(), nullptr) != spaces1.end())
386 throw std::runtime_error("Could not deduce all block trial spaces.");
387
388 return {spaces0, spaces1};
389}
390
392template <typename U, typename V, typename W>
393FunctionSpace(U mesh, V element, W dofmap)
394 -> FunctionSpace<typename std::remove_cvref<
395 typename U::element_type>::type::geometry_type::value_type>;
396
397} // namespace dolfinx::fem
Degree-of-freedom map representations and tools.
A CoordinateElement manages coordinate mappings for isoparametric cells.
Definition CoordinateElement.h:39
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:168
Finite Element, containing the dof layout on a reference element, and various methods for evaluating ...
Definition FiniteElement.h:27
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:35
std::vector< int > component() const
Get the component with respect to the root superspace.
Definition FunctionSpace.h:159
FunctionSpace & operator=(FunctionSpace &&V)=default
Move assignment operator.
std::shared_ptr< const DofMap > dofmap() const
The dofmap.
Definition FunctionSpace.h:309
std::shared_ptr< const FiniteElement< T > > element() const
The finite element.
Definition FunctionSpace.h:306
bool contains(const FunctionSpace &V) const
Check whether V is subspace of this, or this itself.
Definition FunctionSpace.h:107
virtual ~FunctionSpace()=default
Destructor.
std::vector< T > tabulate_dof_coordinates(bool transpose) const
Tabulate the physical coordinates of all dofs on this process.
Definition FunctionSpace.h:172
std::pair< FunctionSpace< T >, 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:141
FunctionSpace(std::shared_ptr< const mesh::Mesh< T > > mesh, std::shared_ptr< const FiniteElement< T > > element, std::shared_ptr< const DofMap > dofmap)
Create function space for given mesh, element and dofmap.
Definition FunctionSpace.h:41
std::shared_ptr< const mesh::Mesh< T > > mesh() const
The mesh.
Definition FunctionSpace.h:303
std::shared_ptr< FunctionSpace< T > > sub(const std::vector< int > &component) const
Extract subspace for component.
Definition FunctionSpace.h:68
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:25
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:343