GridFormat 0.4.0
I/O-Library for grid-like data structures
Loading...
Searching...
No Matches
dolfinx.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2// SPDX-License-Identifier: MIT
8#ifndef GRIDFORMAT_TRAITS_DOLFINX_HPP_
9#define GRIDFORMAT_TRAITS_DOLFINX_HPP_
10
11#include <array>
12#include <ranges>
13#include <cstdint>
14#include <utility>
15#include <memory>
16#include <algorithm>
17
18#include <dolfinx/io/cells.h>
19#include <dolfinx/io/vtk_utils.h>
20#include <dolfinx/mesh/cell_types.h>
21#include <dolfinx/mesh/Mesh.h>
22#include <dolfinx/fem/Function.h>
23#include <dolfinx/fem/FunctionSpace.h>
24
25#include <gridformat/common/ranges.hpp>
26#include <gridformat/common/exceptions.hpp>
27#include <gridformat/common/precision.hpp>
28#include <gridformat/grid/cell_type.hpp>
29#include <gridformat/grid/traits.hpp>
31
32
33namespace GridFormat {
34
35namespace DolfinX {
36
37struct Cell { std::int32_t index; };
38struct Point { std::int32_t index; };
39
40#ifndef DOXYGEN
41namespace Detail {
42 CellType cell_type(dolfinx::mesh::CellType ct) {
43 // since dolfinx supports higher-order cells, let's simply return the lagrange variants always
44 switch (ct) {
45 case dolfinx::mesh::CellType::point: return GridFormat::CellType::vertex;
46 case dolfinx::mesh::CellType::interval: return GridFormat::CellType::lagrange_segment;
47 case dolfinx::mesh::CellType::triangle: return GridFormat::CellType::lagrange_triangle;
48 case dolfinx::mesh::CellType::quadrilateral: return GridFormat::CellType::lagrange_quadrilateral;
49 case dolfinx::mesh::CellType::tetrahedron: return GridFormat::CellType::lagrange_tetrahedron;
50 case dolfinx::mesh::CellType::pyramid: break;
51 case dolfinx::mesh::CellType::prism: break;
52 case dolfinx::mesh::CellType::hexahedron: return GridFormat::CellType::lagrange_hexahedron;
53 }
54 throw NotImplemented("Support for dolfinx cell type '" + dolfinx::mesh::to_string(ct) + "'");
55 }
56
57 bool is_cellwise_constant(const dolfinx::fem::FunctionSpace& space) {
58 return space.dofmap()->element_dof_layout().num_dofs() == 1;
59 }
60
61 template<typename T>
62 bool is_cellwise_constant(const dolfinx::fem::Function<T>& f) {
63 assert(f.function_space());
64 return is_cellwise_constant(*f.function_space());
65 }
66} // namespace Detail
67#endif // DOXYGEN
68
69} // namespace DolfinX
70
71namespace Traits {
72
73template<>
74struct Cells<dolfinx::mesh::Mesh> {
75 static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh) {
76 auto map = mesh.topology().index_map(mesh.topology().dim());
77 if (!map) throw GridFormat::ValueError("Cell index map not available");
78 return std::views::iota(0, map->size_local()) | std::views::transform([] (std::int32_t i) {
79 return DolfinX::Cell{i};
80 });
81 }
82};
83
84template<>
85struct CellType<dolfinx::mesh::Mesh, DolfinX::Cell> {
86 static GridFormat::CellType get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell&) {
87 return DolfinX::Detail::cell_type(mesh.topology().cell_type());
88 }
89};
90
91template<>
92struct CellPoints<dolfinx::mesh::Mesh, DolfinX::Cell> {
93 static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell& cell) {
94 std::span links = mesh.geometry().dofmap().links(cell.index);
95 auto permutation = dolfinx::io::cells::transpose(
96 dolfinx::io::cells::perm_vtk(mesh.topology().cell_type(), links.size())
97 );
98 return std::views::iota(std::size_t{0}, links.size())
99 | std::views::transform([perm = std::move(permutation), links=links] (std::size_t i) {
100 return DolfinX::Point{links[perm[i]]};
101 }
102 );
103 }
104};
105
106template<>
107struct Points<dolfinx::mesh::Mesh> {
108 static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh) {
109 const auto num_points = mesh.geometry().x().size()/3;
110 return std::views::iota(std::size_t{0}, num_points) | std::views::transform([] (std::size_t i) {
111 return DolfinX::Point{static_cast<std::int32_t>(i)};
112 });
113 }
114};
115
116template<>
117struct PointCoordinates<dolfinx::mesh::Mesh, DolfinX::Point> {
118 static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Point& point) {
119 return std::array{
120 mesh.geometry().x()[point.index*3],
121 mesh.geometry().x()[point.index*3 + 1],
122 mesh.geometry().x()[point.index*3 + 2]
123 };
124 }
125};
126
127template<>
128struct PointId<dolfinx::mesh::Mesh, DolfinX::Point> {
129 static std::integral auto get(const dolfinx::mesh::Mesh&, const DolfinX::Point& point) {
130 return point.index;
131 }
132};
133
134template<>
135struct NumberOfPoints<dolfinx::mesh::Mesh> {
136 static std::integral auto get(const dolfinx::mesh::Mesh& mesh) {
137 return mesh.geometry().x().size()/3;
138 }
139};
140
141template<>
142struct NumberOfCells<dolfinx::mesh::Mesh> {
143 static std::integral auto get(const dolfinx::mesh::Mesh& mesh) {
144 auto map = mesh.topology().index_map(mesh.topology().dim());
145 if (!map) throw GridFormat::ValueError("Cell index map not available");
146 return map->size_local();
147 }
148};
149
150template<>
151struct NumberOfCellPoints<dolfinx::mesh::Mesh, DolfinX::Cell> {
152 static std::integral auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell& cell) {
153 return mesh.geometry().dofmap().links(cell.index).size();
154 }
155};
156
157} // namespace Traits
158
159namespace DolfinX {
160
167 public:
168 LagrangePolynomialGrid() = default;
169 LagrangePolynomialGrid(const dolfinx::fem::FunctionSpace& space) {
170 if (!space.mesh() || !space.element())
171 throw ValueError("Cannot construct mesh from space without mesh or element");
172
173 _cell_type = space.mesh()->topology().cell_type();
174 _mesh = space.mesh();
175 _element = space.element();
176
177 auto [x, xshape, xids, _ghosts, cells, cells_shape] = dolfinx::io::vtk_mesh_from_space(space);
178 _node_coords = std::move(x);
179 _node_coords_shape = std::move(xshape);
180 _node_ids = std::move(xids);
181 _cells = std::move(cells);
182 _cells_shape = std::move(cells_shape);
183 _set = true;
184 }
185
186 void update(const dolfinx::fem::FunctionSpace& space) {
187 *this = LagrangePolynomialGrid{space};
188 }
189
190 void clear() {
191 _mesh = nullptr;
192 _element = nullptr;
193 _node_coords.clear();
194 _node_ids.clear();
195 _cells.clear();
196 _set = false;
197 }
198
199 static LagrangePolynomialGrid from(const dolfinx::fem::FunctionSpace& space) {
200 return {space};
201 }
202
203 std::integral auto number_of_points() const { return _node_coords_shape[0]; }
204 std::integral auto number_of_cells() const { return _cells_shape[0]; }
205 std::integral auto number_of_cell_points() const { return _cells_shape[1]; }
206
207 std::int64_t id(const Point& p) const {
208 return _node_ids[p.index];
209 }
210
211 auto position(const Point& p) const {
212 assert(_node_coords_shape[1] == 3);
213 return std::array{
214 _node_coords[p.index*3],
215 _node_coords[p.index*3 + 1],
216 _node_coords[p.index*3 + 2]
217 };
218 }
219
220 std::ranges::range auto points() const {
221 _check_built();
222 return std::views::iota(std::size_t{0}, _node_coords_shape[0])
223 | std::views::transform([] (std::size_t i) {
224 return Point{static_cast<std::int32_t>(i)};
225 });
226 }
227
228 std::ranges::range auto points(const Cell& cell) const {
229 return std::views::iota(std::size_t{0}, _cells_shape[1])
230 | std::views::transform([&, i=cell.index, ncorners=_cells_shape[1]] (std::size_t local_point_index) {
231 return Point{static_cast<std::int32_t>(_cells[i*ncorners + local_point_index])};
232 });
233 }
234
235 std::ranges::range auto cells() const {
236 _check_built();
237 return std::views::iota(std::size_t{0}, _cells_shape[0])
238 | std::views::transform([] (std::size_t i) {
239 return Cell{static_cast<std::int32_t>(i)};
240 });
241 }
242
243 template<int rank = 0, int dim = 3, Concepts::Scalar T>
244 auto evaluate(const dolfinx::fem::Function<T>& f, const Cell& c) const {
245 assert(is_compatible(f));
246 return _evaluate<rank, dim>(f, c.index);
247 }
248
249 template<int rank = 0, int dim = 3, Concepts::Scalar T>
250 auto evaluate(const dolfinx::fem::Function<T>& f, const Point& p) const {
251 assert(is_compatible(f));
252 return _evaluate<rank, dim>(f, p.index);
253 }
254
255 auto cell_type() const {
256 return _cell_type;
257 }
258
259 template<Concepts::Scalar T>
260 bool is_compatible(const dolfinx::fem::Function<T>& f) const {
261 if (!_set) return false;
262 if (!f.function_space()->mesh()) return false;
263 if (f.function_space()->mesh() != _mesh) return false;
264 if (!Detail::is_cellwise_constant(f)) {
265 if (!f.function_space()->element()) return false;
266 if (*f.function_space()->element() != *_element) return false;
267 }
268 return true;
269 }
270
271 private:
272 void _check_built() const {
273 if (!_set)
274 throw InvalidState("Mesh has not been built");
275 }
276
277 template<int rank = 0, int dim = 3, Concepts::Scalar T>
278 auto _evaluate(const dolfinx::fem::Function<T>& f, const std::integral auto i) const {
279 const auto f_components = f.function_space()->element()->block_size();
280 assert(f.function_space()->element()->value_shape().size() == rank);
281 assert(f.x()->array().size() >= static_cast<std::size_t>(i*f_components + f_components));
282
283 if constexpr (rank == 0)
284 return f.x()->array()[i*f_components];
285 else if constexpr (rank == 1) {
286 auto result = Ranges::filled_array<dim>(T{0});
287 std::copy_n(
288 f.x()->array().data() + i*f_components,
289 std::min(dim, f_components),
290 result.begin()
291 );
292 return result;
293 } else {
294 throw NotImplemented("Tensor evaluation");
295 return std::array<std::array<T, dim>, dim>{}; // for return type deduction
296 }
297 }
298
299 dolfinx::mesh::CellType _cell_type;
300 std::shared_ptr<const dolfinx::mesh::Mesh> _mesh{nullptr};
301 std::shared_ptr<const dolfinx::fem::FiniteElement> _element{nullptr};
302
303 std::vector<double> _node_coords;
304 std::array<std::size_t, 2> _node_coords_shape;
305 std::vector<std::int64_t> _node_ids;
306 std::vector<std::int64_t> _cells;
307 std::array<std::size_t, 2> _cells_shape;
308 bool _set = false;
309};
310
319template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
320void set_point_function(const dolfinx::fem::Function<T>& f,
321 Writer& writer,
322 std::string name = "",
323 const Precision<P>& prec = {}) {
324 if (!writer.grid().is_compatible(f))
325 throw ValueError("Grid passed to writer is incompatible with the given function");
326 if (Detail::is_cellwise_constant(f))
327 throw ValueError("Given function is not node-based");
328 if (name.empty())
329 name = f.name;
330 const auto block_size = f.function_space()->element()->block_size();
331 const auto dim = f.function_space()->mesh()->geometry().dim();
332 if (block_size == 1)
333 writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<0>(f, p); }, prec);
334 else if (dim >= block_size)
335 writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<1>(f, p); }, prec);
336 else
337 writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec);
338}
339
348template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
349void set_cell_function(const dolfinx::fem::Function<T>& f,
350 Writer& writer,
351 std::string name = "",
352 const Precision<P>& prec = {}) {
353 if (!writer.grid().is_compatible(f))
354 throw ValueError("Grid passed to writer is incompatible with the given function");
355 if (!Detail::is_cellwise_constant(f))
356 throw ValueError("Given function is not constant per grid cell");
357 if (name.empty())
358 name = f.name;
359 const auto block_size = f.function_space()->element()->block_size();
360 const auto dim = f.function_space()->mesh()->geometry().dim();
361 if (block_size == 1)
362 writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<0>(f, p); }, prec);
363 else if (dim >= block_size)
364 writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<1>(f, p); }, prec);
365 else
366 writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec);
367}
368
377template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
378void set_function(const dolfinx::fem::Function<T>& f,
379 Writer& writer,
380 const std::string& name = "",
381 const Precision<P>& prec = {}) {
382 if (Detail::is_cellwise_constant(f))
383 set_cell_function(f, writer, name, prec);
384 else
385 set_point_function(f, writer, name, prec);
386}
387
388} // namespace DolfinX
389
390namespace Traits {
391
392template<>
393struct Cells<DolfinX::LagrangePolynomialGrid> {
394 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh) {
395 return mesh.cells();
396 }
397};
398
399template<>
400struct CellType<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> {
401 static GridFormat::CellType get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell&) {
402 return DolfinX::Detail::cell_type(mesh.cell_type());
403 }
404};
405
406template<>
407struct CellPoints<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> {
408 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell& cell) {
409 return mesh.points(cell);
410 }
411};
412
413template<>
414struct Points<DolfinX::LagrangePolynomialGrid> {
415 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh) {
416 return mesh.points();
417 }
418};
419
420template<>
421struct PointCoordinates<DolfinX::LagrangePolynomialGrid, DolfinX::Point> {
422 static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Point& point) {
423 return mesh.position(point);
424 }
425};
426
427template<>
428struct PointId<DolfinX::LagrangePolynomialGrid, DolfinX::Point> {
429 static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Point& point) {
430 return mesh.id(point);
431 }
432};
433
434template<>
435struct NumberOfPoints<DolfinX::LagrangePolynomialGrid> {
436 static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh) {
437 return mesh.number_of_points();
438 }
439};
440
441template<>
442struct NumberOfCells<DolfinX::LagrangePolynomialGrid> {
443 static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh) {
444 return mesh.number_of_cells();
445 }
446};
447
448template<>
449struct NumberOfCellPoints<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> {
450 static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell&) {
451 return mesh.number_of_cell_points();
452 }
453};
454
455} // namespace Traits
456} // namespace GridFormat
457
458#endif // GRIDFORMAT_TRAITS_DOLFINX_HPP_
Wrapper around a nodal dolfinx::FunctionSpace, exposing it as a mesh composed of lagrange elements wi...
Definition: dolfinx.hpp:166
Interface to the writers for all supported file formats. Depending on the chosen format,...
Definition: writer.hpp:90
void set_point_field(const std::string &name, Field &&field)
Set a point data field to be added to the output.
Definition: writer.hpp:210
const Grid & grid() const
Return a reference to the underlying grid.
Definition: writer.hpp:326
void set_cell_field(const std::string &name, Field &&field)
Set a cell data field to be added to the output.
Definition: writer.hpp:239
Base classes for grid data writers.
void set_function(const dolfinx::fem::Function< T > &f, Writer &writer, const std::string &name="", const Precision< P > &prec={})
Insert the given function into the writer as field.
Definition: dolfinx.hpp:378
void set_cell_function(const dolfinx::fem::Function< T > &f, Writer &writer, std::string name="", const Precision< P > &prec={})
Insert the given function into the writer as cell field.
Definition: dolfinx.hpp:349
void set_point_function(const dolfinx::fem::Function< T > &f, Writer &writer, std::string name="", const Precision< P > &prec={})
Insert the given function into the writer as point field.
Definition: dolfinx.hpp:320
Definition: dolfinx.hpp:37
Definition: dolfinx.hpp:38