GCC Code Coverage Report

Directory: gridformat/
File: gridformat/traits/dolfinx.hpp
Date: 2026-10-12 17:50:57
Exec Total Coverage
Lines: 193 201 96.0%
Functions: 56 59 94.9%
Branches: 111 219 50.7%
Line Branch Exec Source
1 // SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2 // SPDX-License-Identifier: MIT
3 /*!
4 * \file
5 * \ingroup PredefinedTraits
6 * \brief Traits specializations for <a href="https://docs.fenicsproject.org/dolfinx/v0.11.0/cpp/">dolfinx meshes</a>
7 * \note These traits are written for and tested with dolfinx 0.11.
8 */
9 #ifndef GRIDFORMAT_TRAITS_DOLFINX_HPP_
10 #define GRIDFORMAT_TRAITS_DOLFINX_HPP_
11
12 #include <array>
13 #include <ranges>
14 #include <cstdint>
15 #include <utility>
16 #include <memory>
17 #include <algorithm>
18 #include <concepts>
19 #include <span>
20 #include <vector>
21
22 #include <dolfinx/io/cells.h>
23 #include <dolfinx/io/vtk_utils.h>
24 #include <dolfinx/mesh/cell_types.h>
25 #include <dolfinx/mesh/Mesh.h>
26 #include <dolfinx/fem/Function.h>
27 #include <dolfinx/fem/FunctionSpace.h>
28
29 #include <gridformat/common/ranges.hpp>
30 #include <gridformat/common/exceptions.hpp>
31 #include <gridformat/common/precision.hpp>
32 #include <gridformat/grid/cell_type.hpp>
33 #include <gridformat/grid/traits.hpp>
34 #include <gridformat/grid/writer.hpp>
35
36
37 namespace GridFormat {
38
39 namespace DolfinX {
40
41 struct Cell { std::int32_t index; };
42 struct Point { std::int32_t index; };
43
44 #ifndef DOXYGEN
45 namespace Detail {
46 3139 CellType cell_type(dolfinx::mesh::CellType ct) {
47 // since dolfinx supports higher-order cells, let's simply return the lagrange variants always
48
5/9
✗ Branch 0 not taken.
✓ Branch 1 taken 31 times.
✓ Branch 2 taken 193 times.
✓ Branch 3 taken 97 times.
✓ Branch 4 taken 2305 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 513 times.
✗ Branch 8 not taken.
 3139 switch (ct) {
49 ✗ case dolfinx::mesh::CellType::point: return GridFormat::CellType::vertex;
50 31 case dolfinx::mesh::CellType::interval: return GridFormat::CellType::lagrange_segment;
51 193 case dolfinx::mesh::CellType::triangle: return GridFormat::CellType::lagrange_triangle;
52 97 case dolfinx::mesh::CellType::quadrilateral: return GridFormat::CellType::lagrange_quadrilateral;
53 2305 case dolfinx::mesh::CellType::tetrahedron: return GridFormat::CellType::lagrange_tetrahedron;
54 ✗ case dolfinx::mesh::CellType::pyramid: break;
55 ✗ case dolfinx::mesh::CellType::prism: break;
56 513 case dolfinx::mesh::CellType::hexahedron: return GridFormat::CellType::lagrange_hexahedron;
57 }
58 ✗ throw NotImplemented("Support for dolfinx cell type '" + dolfinx::mesh::to_string(ct) + "'");
59 }
60
61 template<std::floating_point T>
62 11923 bool is_cellwise_constant(const dolfinx::fem::FunctionSpace<T>& space) {
63
2/4
✓ Branch 1 taken 11923 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 11923 times.
✗ Branch 7 not taken.
 11923 return space.dofmap()->element_dof_layout().num_dofs() == 1;
64 }
65
66 template<typename T, std::floating_point U>
67 11914 bool is_cellwise_constant(const dolfinx::fem::Function<T, U>& f) {
68
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 11914 times.
 11914 assert(f.function_space());
69
1/2
✓ Branch 3 taken 11914 times.
✗ Branch 4 not taken.
 11914 return is_cellwise_constant(*f.function_space());
70 }
71
72 //! Return the geometry nodes of the given cell (for meshes with a single cell type)
73 template<std::floating_point T>
74 9033 std::span<const std::int32_t> cell_nodes(const dolfinx::mesh::Mesh<T>& mesh, std::int32_t cell) {
75
1/2
✓ Branch 2 taken 9033 times.
✗ Branch 3 not taken.
 9033 const auto dofmap = mesh.geometry().dofmaps().front();
76 9033 const auto num_nodes = dofmap.extent(1);
77 9033 return {dofmap.data_handle() + static_cast<std::size_t>(cell)*num_nodes, num_nodes};
78 }
79 } // namespace Detail
80 #endif // DOXYGEN
81
82 } // namespace DolfinX
83
84 namespace Traits {
85
86 template<std::floating_point T>
87 struct Cells<dolfinx::mesh::Mesh<T>> {
88 360 static std::ranges::range auto get(const dolfinx::mesh::Mesh<T>& mesh) {
89
1/2
✓ Branch 6 taken 360 times.
✗ Branch 7 not taken.
 360 auto map = mesh.topology()->index_map(mesh.topology()->dim());
90
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 360 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 360 if (!map) throw GridFormat::ValueError("Cell index map not available");
91
3/6
✓ Branch 1 taken 360 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 360 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 360 times.
✗ Branch 10 not taken.
 360 return std::views::iota(0, map->size_local()) | std::views::transform([] (std::int32_t i) {
92 18066 return DolfinX::Cell{i};
93 720 });
94 360 }
95 };
96
97 template<std::floating_point T>
98 struct CellType<dolfinx::mesh::Mesh<T>, DolfinX::Cell> {
99 3011 static GridFormat::CellType get(const dolfinx::mesh::Mesh<T>& mesh, const DolfinX::Cell&) {
100
2/4
✓ Branch 3 taken 3011 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3011 times.
✗ Branch 7 not taken.
 3011 return DolfinX::Detail::cell_type(mesh.topology()->cell_type());
101 }
102 };
103
104 template<std::floating_point T>
105 struct CellPoints<dolfinx::mesh::Mesh<T>, DolfinX::Cell> {
106 9033 static std::ranges::range auto get(const dolfinx::mesh::Mesh<T>& mesh, const DolfinX::Cell& cell) {
107
1/2
✓ Branch 1 taken 9033 times.
✗ Branch 2 not taken.
 9033 const auto nodes = DolfinX::Detail::cell_nodes(mesh, cell.index);
108
1/2
✓ Branch 2 taken 9033 times.
✗ Branch 3 not taken.
 9033 auto permutation = dolfinx::io::cells::transpose(
109
2/4
✓ Branch 4 taken 9033 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9033 times.
✗ Branch 8 not taken.
 18066 dolfinx::io::cells::perm_vtk(mesh.topology()->cell_type(), nodes.size())
110 );
111
1/2
✓ Branch 2 taken 9033 times.
✗ Branch 3 not taken.
 9033 return std::views::iota(std::size_t{0}, nodes.size())
112
1/4
✓ Branch 3 taken 9033 times.
✗ Branch 4 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 18066 | std::views::transform([perm = std::move(permutation), nodes=nodes] (std::size_t i) {
113 26726 return DolfinX::Point{nodes[perm[i]]};
114 }
115
1/2
✓ Branch 1 taken 9033 times.
✗ Branch 2 not taken.
 18066 );
116 9033 }
117 };
118
119 template<std::floating_point T>
120 struct Points<dolfinx::mesh::Mesh<T>> {
121 180 static std::ranges::range auto get(const dolfinx::mesh::Mesh<T>& mesh) {
122 180 const auto num_points = mesh.geometry().x().size()/3;
123
3/6
✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 180 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 180 times.
✗ Branch 8 not taken.
 360 return std::views::iota(std::size_t{0}, num_points) | std::views::transform([] (std::size_t i) {
124 6321 return DolfinX::Point{static_cast<std::int32_t>(i)};
125 360 });
126 }
127 };
128
129 template<std::floating_point T>
130 struct PointCoordinates<dolfinx::mesh::Mesh<T>, DolfinX::Point> {
131 17577 static std::ranges::range auto get(const dolfinx::mesh::Mesh<T>& mesh, const DolfinX::Point& point) {
132 17577 const auto x = mesh.geometry().x();
133 17577 return std::array{x[point.index*3], x[point.index*3 + 1], x[point.index*3 + 2]};
134 }
135 };
136
137 template<std::floating_point T>
138 struct PointId<dolfinx::mesh::Mesh<T>, DolfinX::Point> {
139 17577 static std::integral auto get(const dolfinx::mesh::Mesh<T>&, const DolfinX::Point& point) {
140 17577 return point.index;
141 }
142 };
143
144 template<std::floating_point T>
145 struct NumberOfPoints<dolfinx::mesh::Mesh<T>> {
146 880 static std::integral auto get(const dolfinx::mesh::Mesh<T>& mesh) {
147 880 return mesh.geometry().x().size()/3;
148 }
149 };
150
151 template<std::floating_point T>
152 struct NumberOfCells<dolfinx::mesh::Mesh<T>> {
153 880 static std::integral auto get(const dolfinx::mesh::Mesh<T>& mesh) {
154
1/2
✓ Branch 6 taken 880 times.
✗ Branch 7 not taken.
 880 auto map = mesh.topology()->index_map(mesh.topology()->dim());
155
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 880 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 880 if (!map) throw GridFormat::ValueError("Cell index map not available");
156 1760 return map->size_local();
157 880 }
158 };
159
160 template<std::floating_point T>
161 struct NumberOfCellPoints<dolfinx::mesh::Mesh<T>, DolfinX::Cell> {
162 9033 static std::integral auto get(const dolfinx::mesh::Mesh<T>& mesh, const DolfinX::Cell&) {
163
1/2
✓ Branch 2 taken 9033 times.
✗ Branch 3 not taken.
 9033 return mesh.geometry().dofmaps().front().extent(1);
164 }
165 };
166
167 } // namespace Traits
168
169 namespace DolfinX {
170
171 /*!
172 * \ingroup PredefinedTraits
173 * \brief Wrapper around a nodal dolfinx::FunctionSpace, exposing it as a mesh
174 * composed of lagrange elements with the order of the given function space.
175 * \tparam T The geometry (coordinate) type of the underlying dolfinx mesh.
176 */
177 template<std::floating_point T>
178 class LagrangePolynomialGrid {
179 public:
180 LagrangePolynomialGrid() = default;
181 9 LagrangePolynomialGrid(const dolfinx::fem::FunctionSpace<T>& space) {
182
6/16
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 9 times.
✓ Branch 10 taken 9 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 9 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 9 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
 9 if (!space.mesh() || !space.element())
183 ✗ throw ValueError("Cannot construct mesh from space without mesh or element");
184 // dolfinx does not check this precondition of vtk_mesh_from_space and fails with undefined behaviour
185
3/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 6 times.
 9 if (Detail::is_cellwise_constant(space))
186
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 3 throw ValueError("Cannot construct mesh from a space with cell-wise constant (P0) functions");
187
188
1/2
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 6 _cell_type = space.mesh()->topology()->cell_type();
189 6 _mesh = space.mesh();
190
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 _element = space.element();
191
192
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 auto [x, xshape, xids, _ghosts, cells, cells_shape] = dolfinx::io::vtk_mesh_from_space(space);
193 6 _node_coords = std::move(x);
194 6 _node_coords_shape = std::move(xshape);
195 6 _node_ids = std::move(xids);
196 6 _cells = std::move(cells);
197 6 _cells_shape = std::move(cells_shape);
198 6 _set = true;
199 21 }
200
201 3 void update(const dolfinx::fem::FunctionSpace<T>& space) {
202
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 *this = LagrangePolynomialGrid{space};
203 3 }
204
205 3 void clear() {
206 3 _mesh = nullptr;
207 3 _element = nullptr;
208 3 _node_coords.clear();
209 3 _node_ids.clear();
210 3 _cells.clear();
211 3 _set = false;
212 3 }
213
214 3 static LagrangePolynomialGrid from(const dolfinx::fem::FunctionSpace<T>& space) {
215 3 return {space};
216 }
217
218 277 std::integral auto number_of_points() const { return _node_coords_shape[0]; }
219 313 std::integral auto number_of_cells() const { return _cells_shape[0]; }
220 256 std::integral auto number_of_cell_points() const { return _cells_shape[1]; }
221
222 6534 std::int64_t id(const Point& p) const {
223 6534 return _node_ids[p.index];
224 }
225
226 1539 auto position(const Point& p) const {
227
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1539 times.
 1539 assert(_node_coords_shape[1] == 3);
228 return std::array{
229 1539 _node_coords[p.index*3],
230 1539 _node_coords[p.index*3 + 1],
231 1539 _node_coords[p.index*3 + 2]
232 1539 };
233 }
234
235 33 std::ranges::range auto points() const {
236 33 _check_built();
237
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 return std::views::iota(std::size_t{0}, _node_coords_shape[0])
238
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
 90 | std::views::transform([] (std::size_t i) {
239 13851 return Point{static_cast<std::int32_t>(i)};
240 60 });
241 }
242
243 128 std::ranges::range auto points(const Cell& cell) const {
244
1/2
✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
 128 return std::views::iota(std::size_t{0}, _cells_shape[1])
245
2/4
✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 128 times.
✗ Branch 6 not taken.
 384 | std::views::transform([&, i=cell.index, ncorners=_cells_shape[1]] (std::size_t local_point_index) {
246 3456 return Point{static_cast<std::int32_t>(_cells[i*ncorners + local_point_index])};
247 256 });
248 }
249
250 45 std::ranges::range auto cells() const {
251 45 _check_built();
252
1/2
✓ Branch 2 taken 42 times.
✗ Branch 3 not taken.
 42 return std::views::iota(std::size_t{0}, _cells_shape[0])
253
2/4
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 42 times.
✗ Branch 5 not taken.
 126 | std::views::transform([] (std::size_t i) {
254 1536 return Cell{static_cast<std::int32_t>(i)};
255 84 });
256 }
257
258 template<int rank = 0, int dim = 3, Concepts::Scalar V>
259 2048 auto evaluate(const dolfinx::fem::Function<V, T>& f, const Cell& c) const {
260
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1024 times.
 2048 assert(is_compatible(f));
261 2048 return _evaluate<rank, dim>(f, c.index);
262 }
263
264 template<int rank = 0, int dim = 3, Concepts::Scalar V>
265 21546 auto evaluate(const dolfinx::fem::Function<V, T>& f, const Point& p) const {
266
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 10773 times.
 21546 assert(is_compatible(f));
267 21546 return _evaluate<rank, dim>(f, p.index);
268 }
269
270 128 auto cell_type() const {
271 128 return _cell_type;
272 }
273
274 template<Concepts::Scalar V>
275 11854 bool is_compatible(const dolfinx::fem::Function<V, T>& f) const {
276
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11854 times.
 11854 if (!_set) return false;
277
1/2
✗ Branch 6 not taken.
✓ Branch 7 taken 11854 times.
 11854 if (!f.function_space()->mesh()) return false;
278
2/2
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 11845 times.
 11854 if (f.function_space()->mesh() != _mesh) return false;
279
2/2
✓ Branch 1 taken 10797 times.
✓ Branch 2 taken 1048 times.
 11845 if (!Detail::is_cellwise_constant(f)) {
280
2/4
✓ Branch 3 taken 10797 times.
✗ Branch 4 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 10797 times.
 10797 if (!f.function_space()->element()) return false;
281
3/6
✓ Branch 4 taken 10797 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 10797 times.
✗ Branch 9 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 10797 times.
 10797 if (*f.function_space()->element() != *_element) return false;
282 }
283 11845 return true;
284 }
285
286 private:
287 78 void _check_built() const {
288
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 72 times.
 78 if (!_set)
289
1/2
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 6 throw InvalidState("Mesh has not been built");
290 72 }
291
292 template<int rank = 0, int dim = 3, Concepts::Scalar V>
293 23594 auto _evaluate(const dolfinx::fem::Function<V, T>& f, const std::integral auto i) const {
294
1/2
✓ Branch 3 taken 11797 times.
✗ Branch 4 not taken.
 23594 const auto f_components = f.function_space()->element()->block_size();
295 23594 const auto& values = f.x()->array();
296
3/6
✓ Branch 3 taken 11797 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 11797 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 11797 times.
 23594 assert(f.function_space()->element()->value_shape().size() == rank);
297
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 11797 times.
 23594 assert(values.size() >= static_cast<std::size_t>(i*f_components + f_components));
298
299 if constexpr (rank == 0)
300 13592 return values[i*f_components];
301 else if constexpr (rank == 1) {
302 10002 auto result = Ranges::filled_array<dim>(V{0});
303 30006 std::copy_n(
304
1/2
✓ Branch 2 taken 5001 times.
✗ Branch 3 not taken.
 10002 values.data() + i*f_components,
305 10002 std::min(dim, f_components),
306 result.begin()
307 );
308 20004 return result;
309 } else {
310 ✗ throw NotImplemented("Tensor evaluation");
311 return std::array<std::array<V, dim>, dim>{}; // for return type deduction
312 }
313 }
314
315 dolfinx::mesh::CellType _cell_type;
316 std::shared_ptr<const dolfinx::mesh::Mesh<T>> _mesh{nullptr};
317 std::shared_ptr<const dolfinx::fem::FiniteElement<T>> _element{nullptr};
318
319 std::vector<T> _node_coords;
320 std::array<std::size_t, 2> _node_coords_shape;
321 std::vector<std::int64_t> _node_ids;
322 std::vector<std::int64_t> _cells;
323 std::array<std::size_t, 2> _cells_shape;
324 bool _set = false;
325 };
326
327 template<std::floating_point T>
328 LagrangePolynomialGrid(const dolfinx::fem::FunctionSpace<T>&) -> LagrangePolynomialGrid<T>;
329
330 /*!
331 * \ingroup PredefinedTraits
332 * \brief Create a LagrangePolynomialGrid from the given function space.
333 */
334 template<std::floating_point T>
335 3 LagrangePolynomialGrid<T> make_lagrange_grid(const dolfinx::fem::FunctionSpace<T>& space) {
336 3 return LagrangePolynomialGrid<T>{space};
337 }
338
339 /*!
340 * \ingroup PredefinedTraits
341 * \brief Insert the given function into the writer as point field.
342 * \param f The function to be inserted
343 * \param writer The writer in which to insert it
344 * \param name The name of the field (defaults to `f.name`)
345 * \param prec The precision with which to write the field (defaults to the function's scalar type)
346 */
347 template<typename Writer, Concepts::Scalar T, std::floating_point U, Concepts::Scalar P = T>
348 54 void set_point_function(const dolfinx::fem::Function<T, U>& f,
349 Writer& writer,
350 std::string name = "",
351 const Precision<P>& prec = {}) {
352
2/2
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 24 times.
 54 if (!writer.grid().is_compatible(f))
353
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 6 throw ValueError("Grid passed to writer is incompatible with the given function");
354
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 21 times.
 48 if (Detail::is_cellwise_constant(f))
355
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 6 throw ValueError("Given function is not node-based");
356
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 15 times.
 42 if (name.empty())
357 12 name = f.name;
358
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 42 const auto block_size = f.function_space()->element()->block_size();
359 42 const auto dim = f.function_space()->mesh()->geometry().dim();
360
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 9 times.
 42 if (block_size == 1)
361
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 4641 writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<0>(f, p); }, prec);
362
1/2
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
 18 else if (dim >= block_size)
363
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 3096 writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<1>(f, p); }, prec);
364 else
365 ✗ writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec);
366 42 }
367
368 /*!
369 * \ingroup PredefinedTraits
370 * \brief Insert the given function into the writer as cell field.
371 * \param f The function to be inserted
372 * \param writer The writer in which to insert it
373 * \param name The name of the field (defaults to `f.name`)
374 * \param prec The precision with which to write the field (defaults to the function's scalar type)
375 */
376 template<typename Writer, Concepts::Scalar T, std::floating_point U, Concepts::Scalar P = T>
377 60 void set_cell_function(const dolfinx::fem::Function<T, U>& f,
378 Writer& writer,
379 std::string name = "",
380 const Precision<P>& prec = {}) {
381
2/2
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 24 times.
 60 if (!writer.grid().is_compatible(f))
382
1/2
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
 12 throw ValueError("Grid passed to writer is incompatible with the given function");
383
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 21 times.
 48 if (!Detail::is_cellwise_constant(f))
384
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 6 throw ValueError("Given function is not constant per grid cell");
385
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
 42 if (name.empty())
386 6 name = f.name;
387
1/2
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
 42 const auto block_size = f.function_space()->element()->block_size();
388 42 const auto dim = f.function_space()->mesh()->geometry().dim();
389
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
 42 if (block_size == 1)
390
1/2
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
 542 writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<0>(f, p); }, prec);
391
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 12 else if (dim >= block_size)
392
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 268 writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<1>(f, p); }, prec);
393 else
394 ✗ writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec);
395 42 }
396
397 /*!
398 * \ingroup PredefinedTraits
399 * \brief Insert the given function into the writer as field.
400 * \param f The function to be inserted
401 * \param writer The writer in which to insert it
402 * \param name The name of the field (defaults to `f.name`)
403 * \param prec The precision with which to write the field (defaults to the function's scalar type)
404 */
405 template<typename Writer, Concepts::Scalar T, std::floating_point U, Concepts::Scalar P = T>
406 42 void set_function(const dolfinx::fem::Function<T, U>& f,
407 Writer& writer,
408 const std::string& name = "",
409 const Precision<P>& prec = {}) {
410
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 9 times.
 42 if (Detail::is_cellwise_constant(f))
411
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✓ Branch 5 taken 3 times.
 30 set_cell_function(f, writer, name, prec);
412 else
413
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 18 set_point_function(f, writer, name, prec);
414 36 }
415
416 } // namespace DolfinX
417
418 namespace Traits {
419
420 template<std::floating_point T>
421 struct Cells<DolfinX::LagrangePolynomialGrid<T>> {
422 39 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh) {
423 39 return mesh.cells();
424 }
425 };
426
427 template<std::floating_point T>
428 struct CellType<DolfinX::LagrangePolynomialGrid<T>, DolfinX::Cell> {
429 128 static GridFormat::CellType get(const DolfinX::LagrangePolynomialGrid<T>& mesh, const DolfinX::Cell&) {
430 128 return DolfinX::Detail::cell_type(mesh.cell_type());
431 }
432 };
433
434 template<std::floating_point T>
435 struct CellPoints<DolfinX::LagrangePolynomialGrid<T>, DolfinX::Cell> {
436 128 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh, const DolfinX::Cell& cell) {
437 128 return mesh.points(cell);
438 }
439 };
440
441 template<std::floating_point T>
442 struct Points<DolfinX::LagrangePolynomialGrid<T>> {
443 27 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh) {
444 27 return mesh.points();
445 }
446 };
447
448 template<std::floating_point T>
449 struct PointCoordinates<DolfinX::LagrangePolynomialGrid<T>, DolfinX::Point> {
450 1539 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh, const DolfinX::Point& point) {
451 1539 return mesh.position(point);
452 }
453 };
454
455 template<std::floating_point T>
456 struct PointId<DolfinX::LagrangePolynomialGrid<T>, DolfinX::Point> {
457 6534 static std::integral auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh, const DolfinX::Point& point) {
458 6534 return mesh.id(point);
459 }
460 };
461
462 template<std::floating_point T>
463 struct NumberOfPoints<DolfinX::LagrangePolynomialGrid<T>> {
464 277 static std::integral auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh) {
465 277 return mesh.number_of_points();
466 }
467 };
468
469 template<std::floating_point T>
470 struct NumberOfCells<DolfinX::LagrangePolynomialGrid<T>> {
471 313 static std::integral auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh) {
472 313 return mesh.number_of_cells();
473 }
474 };
475
476 template<std::floating_point T>
477 struct NumberOfCellPoints<DolfinX::LagrangePolynomialGrid<T>, DolfinX::Cell> {
478 256 static std::integral auto get(const DolfinX::LagrangePolynomialGrid<T>& mesh, const DolfinX::Cell&) {
479 256 return mesh.number_of_cell_points();
480 }
481 };
482
483 } // namespace Traits
484 } // namespace GridFormat
485
486 #endif // GRIDFORMAT_TRAITS_DOLFINX_HPP_
487