8#ifndef GRIDFORMAT_TRAITS_DOLFINX_HPP_
9#define GRIDFORMAT_TRAITS_DOLFINX_HPP_
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>
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>
37struct Cell { std::int32_t index; };
38struct Point { std::int32_t index; };
42 CellType cell_type(dolfinx::mesh::CellType 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;
54 throw NotImplemented(
"Support for dolfinx cell type '" + dolfinx::mesh::to_string(ct) +
"'");
57 bool is_cellwise_constant(
const dolfinx::fem::FunctionSpace& space) {
58 return space.dofmap()->element_dof_layout().num_dofs() == 1;
62 bool is_cellwise_constant(
const dolfinx::fem::Function<T>& f) {
63 assert(f.function_space());
64 return is_cellwise_constant(*f.function_space());
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) {
86 static GridFormat::CellType get(
const dolfinx::mesh::Mesh& mesh,
const DolfinX::Cell&) {
87 return DolfinX::Detail::cell_type(mesh.topology().cell_type());
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())
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) {
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) {
118 static std::ranges::range
auto get(
const dolfinx::mesh::Mesh& mesh,
const DolfinX::Point& point) {
120 mesh.geometry().x()[point.index*3],
121 mesh.geometry().x()[point.index*3 + 1],
122 mesh.geometry().x()[point.index*3 + 2]
129 static std::integral
auto get(
const dolfinx::mesh::Mesh&,
const DolfinX::Point& point) {
135struct NumberOfPoints<dolfinx::mesh::Mesh> {
136 static std::integral
auto get(
const dolfinx::mesh::Mesh& mesh) {
137 return mesh.geometry().x().size()/3;
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();
152 static std::integral
auto get(
const dolfinx::mesh::Mesh& mesh,
const DolfinX::Cell& cell) {
153 return mesh.geometry().dofmap().links(cell.index).size();
170 if (!space.mesh() || !space.element())
171 throw ValueError(
"Cannot construct mesh from space without mesh or element");
173 _cell_type = space.mesh()->topology().cell_type();
174 _mesh = space.mesh();
175 _element = space.element();
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);
186 void update(
const dolfinx::fem::FunctionSpace& space) {
193 _node_coords.clear();
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]; }
207 std::int64_t id(
const Point& p)
const {
208 return _node_ids[p.index];
211 auto position(
const Point& p)
const {
212 assert(_node_coords_shape[1] == 3);
214 _node_coords[p.index*3],
215 _node_coords[p.index*3 + 1],
216 _node_coords[p.index*3 + 2]
220 std::ranges::range
auto points()
const {
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)};
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])};
235 std::ranges::range
auto cells()
const {
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)};
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);
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);
255 auto cell_type()
const {
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;
272 void _check_built()
const {
274 throw InvalidState(
"Mesh has not been built");
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));
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});
288 f.x()->array().data() + i*f_components,
289 std::min(dim, f_components),
294 throw NotImplemented(
"Tensor evaluation");
295 return std::array<std::array<T, dim>, dim>{};
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};
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;
319template<
typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
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");
330 const auto block_size = f.function_space()->element()->block_size();
331 const auto dim = f.function_space()->mesh()->geometry().dim();
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);
337 writer.
set_point_field(name, [&] (
const auto p) {
return writer.
grid().
template evaluate<2>(f, p); }, prec);
348template<
typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
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");
359 const auto block_size = f.function_space()->element()->block_size();
360 const auto dim = f.function_space()->mesh()->geometry().dim();
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);
366 writer.
set_cell_field(name, [&] (
const auto p) {
return writer.
grid().
template evaluate<2>(f, p); }, prec);
377template<
typename Writer, Concepts::Scalar T, Concepts::Scalar P = T>
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);
385 set_point_function(f, writer, name, prec);
393struct Cells<DolfinX::LagrangePolynomialGrid> {
402 return DolfinX::Detail::cell_type(mesh.cell_type());
409 return mesh.points(cell);
414struct Points<DolfinX::LagrangePolynomialGrid> {
416 return mesh.points();
423 return mesh.position(point);
430 return mesh.id(point);
435struct NumberOfPoints<DolfinX::LagrangePolynomialGrid> {
437 return mesh.number_of_points();
442struct NumberOfCells<DolfinX::LagrangePolynomialGrid> {
444 return mesh.number_of_cells();
451 return mesh.number_of_cell_points();
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