| 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.6.0/cpp/">dolfinx meshes</a> | ||
| 7 | */ | ||
| 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> | ||
| 30 | #include <gridformat/grid/writer.hpp> | ||
| 31 | |||
| 32 | |||
| 33 | namespace GridFormat { | ||
| 34 | |||
| 35 | namespace DolfinX { | ||
| 36 | |||
| 37 | struct Cell { std::int32_t index; }; | ||
| 38 | struct Point { std::int32_t index; }; | ||
| 39 | |||
| 40 | #ifndef DOXYGEN | ||
| 41 | namespace Detail { | ||
| 42 | 4642 | CellType cell_type(dolfinx::mesh::CellType ct) { | |
| 43 | // since dolfinx supports higher-order cells, let's simply return the lagrange variants always | ||
| 44 |
5/9✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 289 times.
✓ Branch 3 taken 145 times.
✓ Branch 4 taken 3457 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 705 times.
✗ Branch 8 not taken.
|
4642 | switch (ct) { |
| 45 | ✗ | case dolfinx::mesh::CellType::point: return GridFormat::CellType::vertex; | |
| 46 | 46 | case dolfinx::mesh::CellType::interval: return GridFormat::CellType::lagrange_segment; | |
| 47 | 289 | case dolfinx::mesh::CellType::triangle: return GridFormat::CellType::lagrange_triangle; | |
| 48 | 145 | case dolfinx::mesh::CellType::quadrilateral: return GridFormat::CellType::lagrange_quadrilateral; | |
| 49 | 3457 | 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 | 705 | 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 | 11914 | bool is_cellwise_constant(const dolfinx::fem::FunctionSpace& space) { | |
| 58 |
2/4✓ Branch 1 taken 11914 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 11914 times.
✗ Branch 7 not taken.
|
11914 | return space.dofmap()->element_dof_layout().num_dofs() == 1; |
| 59 | } | ||
| 60 | |||
| 61 | template<typename T> | ||
| 62 | 11914 | bool is_cellwise_constant(const dolfinx::fem::Function<T>& f) { | |
| 63 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 11914 times.
|
11914 | assert(f.function_space()); |
| 64 |
1/2✓ Branch 3 taken 11914 times.
✗ Branch 4 not taken.
|
11914 | return is_cellwise_constant(*f.function_space()); |
| 65 | } | ||
| 66 | } // namespace Detail | ||
| 67 | #endif // DOXYGEN | ||
| 68 | |||
| 69 | } // namespace DolfinX | ||
| 70 | |||
| 71 | namespace Traits { | ||
| 72 | |||
| 73 | template<> | ||
| 74 | struct Cells<dolfinx::mesh::Mesh> { | ||
| 75 | 520 | static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh) { | |
| 76 |
3/6✓ Branch 1 taken 520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 520 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 520 times.
✗ Branch 9 not taken.
|
520 | auto map = mesh.topology().index_map(mesh.topology().dim()); |
| 77 |
1/4✗ Branch 1 not taken.
✓ Branch 2 taken 520 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
520 | if (!map) throw GridFormat::ValueError("Cell index map not available"); |
| 78 |
2/4✓ Branch 1 taken 520 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 520 times.
✗ Branch 7 not taken.
|
520 | return std::views::iota(0, map->size_local()) | std::views::transform([] (std::int32_t i) { |
| 79 | 27084 | return DolfinX::Cell{i}; | |
| 80 |
1/2✓ Branch 1 taken 520 times.
✗ Branch 2 not taken.
|
1040 | }); |
| 81 | 520 | } | |
| 82 | }; | ||
| 83 | |||
| 84 | template<> | ||
| 85 | struct CellType<dolfinx::mesh::Mesh, DolfinX::Cell> { | ||
| 86 | 4514 | static GridFormat::CellType get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell&) { | |
| 87 | 4514 | return DolfinX::Detail::cell_type(mesh.topology().cell_type()); | |
| 88 | } | ||
| 89 | }; | ||
| 90 | |||
| 91 | template<> | ||
| 92 | struct CellPoints<dolfinx::mesh::Mesh, DolfinX::Cell> { | ||
| 93 | 13542 | static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell& cell) { | |
| 94 |
2/4✓ Branch 1 taken 13542 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13542 times.
✗ Branch 5 not taken.
|
13542 | std::span links = mesh.geometry().dofmap().links(cell.index); |
| 95 | auto permutation = dolfinx::io::cells::transpose( | ||
| 96 |
2/4✓ Branch 2 taken 13542 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 13542 times.
✗ Branch 7 not taken.
|
27084 | dolfinx::io::cells::perm_vtk(mesh.topology().cell_type(), links.size()) |
| 97 |
1/2✓ Branch 1 taken 13542 times.
✗ Branch 2 not taken.
|
13542 | ); |
| 98 |
1/2✓ Branch 2 taken 13542 times.
✗ Branch 3 not taken.
|
13542 | return std::views::iota(std::size_t{0}, links.size()) |
| 99 |
1/4✓ Branch 3 taken 13542 times.
✗ Branch 4 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
27084 | | std::views::transform([perm = std::move(permutation), links=links] (std::size_t i) { |
| 100 | 40034 | return DolfinX::Point{links[perm[i]]}; | |
| 101 | } | ||
| 102 |
1/2✓ Branch 1 taken 13542 times.
✗ Branch 2 not taken.
|
27084 | ); |
| 103 | 13542 | } | |
| 104 | }; | ||
| 105 | |||
| 106 | template<> | ||
| 107 | struct Points<dolfinx::mesh::Mesh> { | ||
| 108 | 260 | static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh) { | |
| 109 |
2/4✓ Branch 1 taken 260 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 260 times.
✗ Branch 5 not taken.
|
260 | const auto num_points = mesh.geometry().x().size()/3; |
| 110 |
2/4✓ Branch 1 taken 260 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 260 times.
✗ Branch 5 not taken.
|
520 | return std::views::iota(std::size_t{0}, num_points) | std::views::transform([] (std::size_t i) { |
| 111 | 9558 | return DolfinX::Point{static_cast<std::int32_t>(i)}; | |
| 112 |
1/2✓ Branch 1 taken 260 times.
✗ Branch 2 not taken.
|
520 | }); |
| 113 | } | ||
| 114 | }; | ||
| 115 | |||
| 116 | template<> | ||
| 117 | struct PointCoordinates<dolfinx::mesh::Mesh, DolfinX::Point> { | ||
| 118 | 26389 | static std::ranges::range auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Point& point) { | |
| 119 | return std::array{ | ||
| 120 |
2/4✓ Branch 1 taken 26389 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26389 times.
✗ Branch 5 not taken.
|
26389 | mesh.geometry().x()[point.index*3], |
| 121 |
1/2✓ Branch 1 taken 26389 times.
✗ Branch 2 not taken.
|
26389 | mesh.geometry().x()[point.index*3 + 1], |
| 122 |
1/2✓ Branch 1 taken 26389 times.
✗ Branch 2 not taken.
|
26389 | mesh.geometry().x()[point.index*3 + 2] |
| 123 |
2/4✓ Branch 1 taken 26389 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26389 times.
✗ Branch 5 not taken.
|
79167 | }; |
| 124 | } | ||
| 125 | }; | ||
| 126 | |||
| 127 | template<> | ||
| 128 | struct PointId<dolfinx::mesh::Mesh, DolfinX::Point> { | ||
| 129 | 26389 | static std::integral auto get(const dolfinx::mesh::Mesh&, const DolfinX::Point& point) { | |
| 130 | 26389 | return point.index; | |
| 131 | } | ||
| 132 | }; | ||
| 133 | |||
| 134 | template<> | ||
| 135 | struct NumberOfPoints<dolfinx::mesh::Mesh> { | ||
| 136 | 1280 | static std::integral auto get(const dolfinx::mesh::Mesh& mesh) { | |
| 137 |
2/4✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
|
1280 | return mesh.geometry().x().size()/3; |
| 138 | } | ||
| 139 | }; | ||
| 140 | |||
| 141 | template<> | ||
| 142 | struct NumberOfCells<dolfinx::mesh::Mesh> { | ||
| 143 | 1280 | static std::integral auto get(const dolfinx::mesh::Mesh& mesh) { | |
| 144 |
3/6✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1280 times.
✗ Branch 9 not taken.
|
1280 | auto map = mesh.topology().index_map(mesh.topology().dim()); |
| 145 |
1/4✗ Branch 1 not taken.
✓ Branch 2 taken 1280 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1280 | if (!map) throw GridFormat::ValueError("Cell index map not available"); |
| 146 | 2560 | return map->size_local(); | |
| 147 | 1280 | } | |
| 148 | }; | ||
| 149 | |||
| 150 | template<> | ||
| 151 | struct NumberOfCellPoints<dolfinx::mesh::Mesh, DolfinX::Cell> { | ||
| 152 | 13542 | static std::integral auto get(const dolfinx::mesh::Mesh& mesh, const DolfinX::Cell& cell) { | |
| 153 |
2/4✓ Branch 1 taken 13542 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13542 times.
✗ Branch 5 not taken.
|
13542 | return mesh.geometry().dofmap().links(cell.index).size(); |
| 154 | } | ||
| 155 | }; | ||
| 156 | |||
| 157 | } // namespace Traits | ||
| 158 | |||
| 159 | namespace DolfinX { | ||
| 160 | |||
| 161 | /*! | ||
| 162 | * \ingroup PredefinedTraits | ||
| 163 | * \brief Wrapper around a nodal dolfinx::FunctionSpace, exposing it as a mesh | ||
| 164 | * composed of lagrange elements with the order of the given function space. | ||
| 165 | */ | ||
| 166 | class LagrangePolynomialGrid { | ||
| 167 | public: | ||
| 168 | LagrangePolynomialGrid() = default; | ||
| 169 | 9 | LagrangePolynomialGrid(const dolfinx::fem::FunctionSpace& space) { | |
| 170 |
7/18✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 9 times.
✓ Branch 12 taken 9 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 9 times.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 9 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
9 | if (!space.mesh() || !space.element()) |
| 171 | ✗ | throw ValueError("Cannot construct mesh from space without mesh or element"); | |
| 172 | |||
| 173 |
2/4✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
|
9 | _cell_type = space.mesh()->topology().cell_type(); |
| 174 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | _mesh = space.mesh(); |
| 175 |
1/2✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
|
9 | _element = space.element(); |
| 176 | |||
| 177 |
2/2✓ Branch 1 taken 6 times.
✓ Branch 2 taken 3 times.
|
9 | auto [x, xshape, xids, _ghosts, cells, cells_shape] = dolfinx::io::vtk_mesh_from_space(space); |
| 178 | 6 | _node_coords = std::move(x); | |
| 179 | 6 | _node_coords_shape = std::move(xshape); | |
| 180 | 6 | _node_ids = std::move(xids); | |
| 181 | 6 | _cells = std::move(cells); | |
| 182 | 6 | _cells_shape = std::move(cells_shape); | |
| 183 | 6 | _set = true; | |
| 184 | 21 | } | |
| 185 | |||
| 186 | 3 | void update(const dolfinx::fem::FunctionSpace& space) { | |
| 187 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | *this = LagrangePolynomialGrid{space}; |
| 188 | 3 | } | |
| 189 | |||
| 190 | 3 | void clear() { | |
| 191 | 3 | _mesh = nullptr; | |
| 192 | 3 | _element = nullptr; | |
| 193 | 3 | _node_coords.clear(); | |
| 194 | 3 | _node_ids.clear(); | |
| 195 | 3 | _cells.clear(); | |
| 196 | 3 | _set = false; | |
| 197 | 3 | } | |
| 198 | |||
| 199 | 6 | static LagrangePolynomialGrid from(const dolfinx::fem::FunctionSpace& space) { | |
| 200 | 6 | return {space}; | |
| 201 | } | ||
| 202 | |||
| 203 | 277 | std::integral auto number_of_points() const { return _node_coords_shape[0]; } | |
| 204 | 313 | std::integral auto number_of_cells() const { return _cells_shape[0]; } | |
| 205 | 256 | std::integral auto number_of_cell_points() const { return _cells_shape[1]; } | |
| 206 | |||
| 207 | 6534 | std::int64_t id(const Point& p) const { | |
| 208 | 6534 | return _node_ids[p.index]; | |
| 209 | } | ||
| 210 | |||
| 211 | 1539 | auto position(const Point& p) const { | |
| 212 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1539 times.
|
1539 | assert(_node_coords_shape[1] == 3); |
| 213 | return std::array{ | ||
| 214 | 1539 | _node_coords[p.index*3], | |
| 215 | 1539 | _node_coords[p.index*3 + 1], | |
| 216 | 1539 | _node_coords[p.index*3 + 2] | |
| 217 | 1539 | }; | |
| 218 | } | ||
| 219 | |||
| 220 | 33 | std::ranges::range auto points() const { | |
| 221 | 33 | _check_built(); | |
| 222 |
1/2✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
|
30 | return std::views::iota(std::size_t{0}, _node_coords_shape[0]) |
| 223 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
60 | | std::views::transform([] (std::size_t i) { |
| 224 | 13851 | return Point{static_cast<std::int32_t>(i)}; | |
| 225 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
60 | }); |
| 226 | } | ||
| 227 | |||
| 228 | 128 | std::ranges::range auto points(const Cell& cell) const { | |
| 229 |
1/2✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
|
128 | return std::views::iota(std::size_t{0}, _cells_shape[1]) |
| 230 |
1/2✓ Branch 2 taken 128 times.
✗ Branch 3 not taken.
|
256 | | std::views::transform([&, i=cell.index, ncorners=_cells_shape[1]] (std::size_t local_point_index) { |
| 231 | 3456 | return Point{static_cast<std::int32_t>(_cells[i*ncorners + local_point_index])}; | |
| 232 |
1/2✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
|
256 | }); |
| 233 | } | ||
| 234 | |||
| 235 | 45 | std::ranges::range auto cells() const { | |
| 236 | 45 | _check_built(); | |
| 237 |
1/2✓ Branch 2 taken 42 times.
✗ Branch 3 not taken.
|
42 | return std::views::iota(std::size_t{0}, _cells_shape[0]) |
| 238 |
1/2✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
|
84 | | std::views::transform([] (std::size_t i) { |
| 239 | 1536 | return Cell{static_cast<std::int32_t>(i)}; | |
| 240 |
1/2✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
|
84 | }); |
| 241 | } | ||
| 242 | |||
| 243 | template<int rank = 0, int dim = 3, Concepts::Scalar T> | ||
| 244 | 2048 | auto evaluate(const dolfinx::fem::Function<T>& f, const Cell& c) const { | |
| 245 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 1024 times.
|
2048 | assert(is_compatible(f)); |
| 246 | 2048 | return _evaluate<rank, dim>(f, c.index); | |
| 247 | } | ||
| 248 | |||
| 249 | template<int rank = 0, int dim = 3, Concepts::Scalar T> | ||
| 250 | 21546 | auto evaluate(const dolfinx::fem::Function<T>& f, const Point& p) const { | |
| 251 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 10773 times.
|
21546 | assert(is_compatible(f)); |
| 252 | 21546 | return _evaluate<rank, dim>(f, p.index); | |
| 253 | } | ||
| 254 | |||
| 255 | 128 | auto cell_type() const { | |
| 256 | 128 | return _cell_type; | |
| 257 | } | ||
| 258 | |||
| 259 | template<Concepts::Scalar T> | ||
| 260 | 11854 | bool is_compatible(const dolfinx::fem::Function<T>& f) const { | |
| 261 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11854 times.
|
11854 | if (!_set) return false; |
| 262 |
2/4✓ Branch 3 taken 11854 times.
✗ Branch 4 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11854 times.
|
11854 | if (!f.function_space()->mesh()) return false; |
| 263 |
3/4✓ Branch 3 taken 11854 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 9 times.
✓ Branch 9 taken 11845 times.
|
11854 | if (f.function_space()->mesh() != _mesh) return false; |
| 264 |
2/2✓ Branch 1 taken 10797 times.
✓ Branch 2 taken 1048 times.
|
11845 | if (!Detail::is_cellwise_constant(f)) { |
| 265 |
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; |
| 266 |
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; |
| 267 | } | ||
| 268 | 11845 | return true; | |
| 269 | } | ||
| 270 | |||
| 271 | private: | ||
| 272 | 78 | void _check_built() const { | |
| 273 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 72 times.
|
78 | if (!_set) |
| 274 |
1/2✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
|
6 | throw InvalidState("Mesh has not been built"); |
| 275 | 72 | } | |
| 276 | |||
| 277 | template<int rank = 0, int dim = 3, Concepts::Scalar T> | ||
| 278 | 23594 | auto _evaluate(const dolfinx::fem::Function<T>& f, const std::integral auto i) const { | |
| 279 |
1/2✓ Branch 3 taken 11797 times.
✗ Branch 4 not taken.
|
23594 | const auto f_components = f.function_space()->element()->block_size(); |
| 280 |
2/4✓ Branch 3 taken 11797 times.
✗ Branch 4 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 11797 times.
|
23594 | assert(f.function_space()->element()->value_shape().size() == rank); |
| 281 |
1/2✗ Branch 4 not taken.
✓ Branch 5 taken 11797 times.
|
23594 | assert(f.x()->array().size() >= static_cast<std::size_t>(i*f_components + f_components)); |
| 282 | |||
| 283 | if constexpr (rank == 0) | ||
| 284 | 13592 | return f.x()->array()[i*f_components]; | |
| 285 | else if constexpr (rank == 1) { | ||
| 286 | 10002 | auto result = Ranges::filled_array<dim>(T{0}); | |
| 287 | 30006 | std::copy_n( | |
| 288 |
1/2✓ Branch 5 taken 5001 times.
✗ Branch 6 not taken.
|
10002 | f.x()->array().data() + i*f_components, |
| 289 | 10002 | std::min(dim, f_components), | |
| 290 | result.begin() | ||
| 291 | ); | ||
| 292 | 20004 | 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 | |||
| 311 | /*! | ||
| 312 | * \ingroup PredefinedTraits | ||
| 313 | * \brief Insert the given function into the writer as point field. | ||
| 314 | * \param f The function to be inserted | ||
| 315 | * \param writer The writer in which to insert it | ||
| 316 | * \param name The name of the field (defaults to `f.name`) | ||
| 317 | * \param prec The precision with which to write the field (defaults to the function's scalar type) | ||
| 318 | */ | ||
| 319 | template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T> | ||
| 320 | 54 | void set_point_function(const dolfinx::fem::Function<T>& f, | |
| 321 | Writer& writer, | ||
| 322 | std::string name = "", | ||
| 323 | const Precision<P>& prec = {}) { | ||
| 324 |
2/2✓ Branch 2 taken 3 times.
✓ Branch 3 taken 24 times.
|
54 | if (!writer.grid().is_compatible(f)) |
| 325 |
1/2✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
6 | throw ValueError("Grid passed to writer is incompatible with the given function"); |
| 326 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 21 times.
|
48 | if (Detail::is_cellwise_constant(f)) |
| 327 |
1/2✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
6 | throw ValueError("Given function is not node-based"); |
| 328 |
2/2✓ Branch 1 taken 6 times.
✓ Branch 2 taken 15 times.
|
42 | if (name.empty()) |
| 329 | 12 | name = f.name; | |
| 330 |
1/2✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
42 | const auto block_size = f.function_space()->element()->block_size(); |
| 331 |
3/6✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 21 times.
✗ Branch 11 not taken.
|
42 | const auto dim = f.function_space()->mesh()->geometry().dim(); |
| 332 |
2/2✓ Branch 0 taken 12 times.
✓ Branch 1 taken 9 times.
|
42 | if (block_size == 1) |
| 333 |
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); |
| 334 |
1/2✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
|
18 | else if (dim >= block_size) |
| 335 |
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); |
| 336 | else | ||
| 337 | ✗ | writer.set_point_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec); | |
| 338 | 42 | } | |
| 339 | |||
| 340 | /*! | ||
| 341 | * \ingroup PredefinedTraits | ||
| 342 | * \brief Insert the given function into the writer as cell field. | ||
| 343 | * \param f The function to be inserted | ||
| 344 | * \param writer The writer in which to insert it | ||
| 345 | * \param name The name of the field (defaults to `f.name`) | ||
| 346 | * \param prec The precision with which to write the field (defaults to the function's scalar type) | ||
| 347 | */ | ||
| 348 | template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T> | ||
| 349 | 60 | void set_cell_function(const dolfinx::fem::Function<T>& f, | |
| 350 | Writer& writer, | ||
| 351 | std::string name = "", | ||
| 352 | const Precision<P>& prec = {}) { | ||
| 353 |
2/2✓ Branch 2 taken 6 times.
✓ Branch 3 taken 24 times.
|
60 | if (!writer.grid().is_compatible(f)) |
| 354 |
1/2✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
|
12 | throw ValueError("Grid passed to writer is incompatible with the given function"); |
| 355 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 21 times.
|
48 | if (!Detail::is_cellwise_constant(f)) |
| 356 |
1/2✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
6 | throw ValueError("Given function is not constant per grid cell"); |
| 357 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
|
42 | if (name.empty()) |
| 358 | 6 | name = f.name; | |
| 359 |
1/2✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
42 | const auto block_size = f.function_space()->element()->block_size(); |
| 360 |
3/6✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 21 times.
✗ Branch 11 not taken.
|
42 | const auto dim = f.function_space()->mesh()->geometry().dim(); |
| 361 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 6 times.
|
42 | if (block_size == 1) |
| 362 |
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); |
| 363 |
1/2✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
|
12 | else if (dim >= block_size) |
| 364 |
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); |
| 365 | else | ||
| 366 | ✗ | writer.set_cell_field(name, [&] (const auto p) { return writer.grid().template evaluate<2>(f, p); }, prec); | |
| 367 | 42 | } | |
| 368 | |||
| 369 | /*! | ||
| 370 | * \ingroup PredefinedTraits | ||
| 371 | * \brief Insert the given function into the writer as field. | ||
| 372 | * \param f The function to be inserted | ||
| 373 | * \param writer The writer in which to insert it | ||
| 374 | * \param name The name of the field (defaults to `f.name`) | ||
| 375 | * \param prec The precision with which to write the field (defaults to the function's scalar type) | ||
| 376 | */ | ||
| 377 | template<typename Writer, Concepts::Scalar T, Concepts::Scalar P = T> | ||
| 378 | 42 | void set_function(const dolfinx::fem::Function<T>& f, | |
| 379 | Writer& writer, | ||
| 380 | const std::string& name = "", | ||
| 381 | const Precision<P>& prec = {}) { | ||
| 382 |
2/2✓ Branch 1 taken 12 times.
✓ Branch 2 taken 9 times.
|
42 | if (Detail::is_cellwise_constant(f)) |
| 383 |
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); |
| 384 | else | ||
| 385 |
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); |
| 386 | 36 | } | |
| 387 | |||
| 388 | } // namespace DolfinX | ||
| 389 | |||
| 390 | namespace Traits { | ||
| 391 | |||
| 392 | template<> | ||
| 393 | struct Cells<DolfinX::LagrangePolynomialGrid> { | ||
| 394 | 39 | static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh) { | |
| 395 | 39 | return mesh.cells(); | |
| 396 | } | ||
| 397 | }; | ||
| 398 | |||
| 399 | template<> | ||
| 400 | struct CellType<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> { | ||
| 401 | 128 | static GridFormat::CellType get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell&) { | |
| 402 | 128 | return DolfinX::Detail::cell_type(mesh.cell_type()); | |
| 403 | } | ||
| 404 | }; | ||
| 405 | |||
| 406 | template<> | ||
| 407 | struct CellPoints<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> { | ||
| 408 | 128 | static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell& cell) { | |
| 409 | 128 | return mesh.points(cell); | |
| 410 | } | ||
| 411 | }; | ||
| 412 | |||
| 413 | template<> | ||
| 414 | struct Points<DolfinX::LagrangePolynomialGrid> { | ||
| 415 | 27 | static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh) { | |
| 416 | 27 | return mesh.points(); | |
| 417 | } | ||
| 418 | }; | ||
| 419 | |||
| 420 | template<> | ||
| 421 | struct PointCoordinates<DolfinX::LagrangePolynomialGrid, DolfinX::Point> { | ||
| 422 | 1539 | static std::ranges::range auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Point& point) { | |
| 423 | 1539 | return mesh.position(point); | |
| 424 | } | ||
| 425 | }; | ||
| 426 | |||
| 427 | template<> | ||
| 428 | struct PointId<DolfinX::LagrangePolynomialGrid, DolfinX::Point> { | ||
| 429 | 6534 | static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Point& point) { | |
| 430 | 6534 | return mesh.id(point); | |
| 431 | } | ||
| 432 | }; | ||
| 433 | |||
| 434 | template<> | ||
| 435 | struct NumberOfPoints<DolfinX::LagrangePolynomialGrid> { | ||
| 436 | 277 | static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh) { | |
| 437 | 277 | return mesh.number_of_points(); | |
| 438 | } | ||
| 439 | }; | ||
| 440 | |||
| 441 | template<> | ||
| 442 | struct NumberOfCells<DolfinX::LagrangePolynomialGrid> { | ||
| 443 | 313 | static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh) { | |
| 444 | 313 | return mesh.number_of_cells(); | |
| 445 | } | ||
| 446 | }; | ||
| 447 | |||
| 448 | template<> | ||
| 449 | struct NumberOfCellPoints<DolfinX::LagrangePolynomialGrid, DolfinX::Cell> { | ||
| 450 | 256 | static std::integral auto get(const DolfinX::LagrangePolynomialGrid& mesh, const DolfinX::Cell&) { | |
| 451 | 256 | return mesh.number_of_cell_points(); | |
| 452 | } | ||
| 453 | }; | ||
| 454 | |||
| 455 | } // namespace Traits | ||
| 456 | } // namespace GridFormat | ||
| 457 | |||
| 458 | #endif // GRIDFORMAT_TRAITS_DOLFINX_HPP_ | ||
| 459 |