DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
vtk_utils.h
1// Copyright (C) 2020-2022 Garth N. Wells and Jørgen S. Dokken
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 "cells.h"
10#include "vtk_utils.h"
11#include <algorithm>
12#include <array>
13#include <basix/mdspan.hpp>
14#include <concepts>
15#include <cstdint>
16#include <dolfinx/common/IndexMap.h>
17#include <dolfinx/fem/DofMap.h>
18#include <dolfinx/fem/FiniteElement.h>
19#include <dolfinx/fem/FunctionSpace.h>
20#include <dolfinx/mesh/Geometry.h>
21#include <dolfinx/mesh/Mesh.h>
22#include <dolfinx/mesh/Topology.h>
23#include <span>
24#include <tuple>
25#include <vector>
26
27namespace dolfinx
28{
29namespace fem
30{
31template <std::floating_point T>
33}
34
35namespace mesh
36{
37enum class CellType;
38template <std::floating_point T>
39class Geometry;
40} // namespace mesh
41
42namespace io
43{
44namespace impl
45{
57template <typename T>
58std::tuple<std::vector<T>, std::array<std::size_t, 2>,
59 std::vector<std::int64_t>, std::vector<std::uint8_t>>
60tabulate_lagrange_dof_coordinates(const fem::FunctionSpace<T>& V)
61{
62 auto mesh = V.mesh();
63 assert(mesh);
64 auto topology = mesh->topology();
65 assert(topology);
66 const std::size_t gdim = mesh->geometry().dim();
67 const int tdim = topology->dim();
68
69 // Get dofmap data
70 auto dofmap = V.dofmap();
71 assert(dofmap);
72 auto map_dofs = dofmap->index_map;
73 assert(map_dofs);
74 const int index_map_bs = dofmap->index_map_bs();
75 const int dofmap_bs = dofmap->bs();
76
77 // Get element data
78 auto element = V.element();
79 assert(element);
80 const int e_bs = element->block_size();
81 const std::size_t scalar_dofs = element->space_dimension() / e_bs;
82 const std::size_t num_nodes
83 = index_map_bs * (map_dofs->size_local() + map_dofs->num_ghosts())
84 / dofmap_bs;
85
86 // Get the dof coordinates on the reference element and the mesh
87 // coordinate map
88 const auto [X, Xshape] = element->interpolation_points();
89 const fem::CoordinateElement<T>& cmap = mesh->geometry().cmap();
90
91 // Prepare cell geometry
92 auto dofmap_x = mesh->geometry().dofmap();
93 std::span<const T> x_g = mesh->geometry().x();
94 const std::size_t num_dofs_g = cmap.dim();
95
96 std::span<const std::uint32_t> cell_info;
97 if (element->needs_dof_transformations())
98 {
99 mesh->topology_mutable()->create_entity_permutations();
100 cell_info = std::span(mesh->topology()->get_cell_permutation_info());
101 }
102
103 // Transformation from reference element basis function data to
104 // conforming element basis function function
105 auto apply_dof_transformation
106 = element->template dof_transformation_fn<T>(fem::doftransform::standard);
107
108 using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
109 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
110 using cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
111 T, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>>;
112
113 // Tabulate basis functions at node reference coordinates
114 const std::array<std::size_t, 4> phi_shape
115 = cmap.tabulate_shape(0, Xshape[0]);
116 std::vector<T> phi_b(
117 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
118 cmdspan4_t phi_full(phi_b.data(), phi_shape);
119 cmap.tabulate(0, X, Xshape, phi_b);
120 auto phi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
121 phi_full, 0, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent,
122 MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
123
124 // Loop over cells and tabulate dofs
125 auto map = topology->index_map(tdim);
126 assert(map);
127 const std::int32_t num_cells = map->size_local() + map->num_ghosts();
128 std::vector<T> x_b(scalar_dofs * gdim);
129 mdspan2_t x(x_b.data(), scalar_dofs, gdim);
130 std::vector<T> coordinate_dofs_b(num_dofs_g * gdim);
131 mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
132
133 std::vector<T> coords(num_nodes * 3, 0.0);
134 std::array<std::size_t, 2> cshape = {num_nodes, 3};
135 for (std::int32_t c = 0; c < num_cells; ++c)
136 {
137 // Extract cell geometry
138 for (std::size_t i = 0; i < dofmap_x.extent(1); ++i)
139 for (std::size_t j = 0; j < gdim; ++j)
140 coordinate_dofs(i, j) = x_g[3 * dofmap_x(c, i) + j];
141
142 // Tabulate dof coordinates on cell
143 cmap.push_forward(x, coordinate_dofs, phi);
144 apply_dof_transformation(x_b, std::span(cell_info.data(), cell_info.size()),
145 c, x.extent(1));
146
147 // Copy dof coordinates into vector
148 auto dofs = dofmap->cell_dofs(c);
149 for (std::size_t i = 0; i < dofs.size(); ++i)
150 {
151 std::int32_t dof = dofs[i];
152 for (std::size_t j = 0; j < gdim; ++j)
153 coords[3 * dof + j] = x(i, j);
154 }
155 }
156
157 // Original point IDs
158 std::vector<std::int64_t> x_id(num_nodes);
159 std::array<std::int64_t, 2> range = map_dofs->local_range();
160 std::int32_t size_local = range[1] - range[0];
161 std::iota(x_id.begin(), std::next(x_id.begin(), size_local), range[0]);
162 std::span ghosts = map_dofs->ghosts();
163 std::ranges::copy(ghosts, std::next(x_id.begin(), size_local));
164
165 // Ghosts
166 std::vector<std::uint8_t> id_ghost(num_nodes, 0);
167 std::fill(std::next(id_ghost.begin(), size_local), id_ghost.end(), 1);
168
169 return {std::move(coords), cshape, std::move(x_id), std::move(id_ghost)};
170}
171
172} // namespace impl
173
188template <typename T>
189std::tuple<std::vector<T>, std::array<std::size_t, 2>,
190 std::vector<std::int64_t>, std::vector<std::uint8_t>,
191 std::vector<std::int64_t>, std::array<std::size_t, 2>>
193{
194 auto mesh = V.mesh();
195 assert(mesh);
196 auto topology = mesh->topology();
197 assert(topology);
198 const int tdim = topology->dim();
199
200 assert(V.element());
201 if (V.element()->is_mixed())
202 throw std::runtime_error("Can't create VTK mesh from a mixed element");
203
204 const auto [x, xshape, x_id, x_ghost]
205 = impl::tabulate_lagrange_dof_coordinates(V);
206 auto map = topology->index_map(tdim);
207 const std::size_t num_cells = map->size_local() + map->num_ghosts();
208
209 // Create permutation from DOLFINx dof ordering to VTK
210 auto dofmap = V.dofmap();
211 assert(dofmap);
212 const int element_block_size = V.element()->block_size();
213 const std::uint32_t num_nodes
214 = V.element()->space_dimension() / element_block_size;
215 const std::vector<std::uint16_t> vtkmap = io::cells::transpose(
216 io::cells::perm_vtk(topology->cell_type(), num_nodes));
217
218 // Extract topology for all local cells as
219 // [v0_0, ...., v0_N0, v1_0, ...., v1_N1, ....]
220 std::array<std::size_t, 2> shape = {num_cells, num_nodes};
221 std::vector<std::int64_t> vtk_topology(shape[0] * shape[1]);
222 for (std::size_t c = 0; c < shape[0]; ++c)
223 {
224 auto dofs = dofmap->cell_dofs(c);
225 for (std::size_t i = 0; i < dofs.size(); ++i)
226 vtk_topology[c * shape[1] + i] = dofs[vtkmap[i]];
227 }
228
229 return {std::move(x),
230 xshape,
231 std::move(x_id),
232 std::move(x_ghost),
233 std::move(vtk_topology),
234 shape};
235}
236
251std::pair<std::vector<std::int64_t>, std::array<std::size_t, 2>>
253 MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
254 const std::int32_t,
255 MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>
256 dofmap_x,
257 mesh::CellType cell_type);
258} // namespace io
259} // namespace dolfinx
Degree-of-freedom map representations and tools.
Definition XDMFFile.h:27
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:55
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:48
int dim() const
The dimension of the coordinate element space.
Definition CoordinateElement.cpp:206
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:188
This class represents a finite element function space defined by a mesh, a finite element,...
Definition vtk_utils.h:32
std::vector< std::uint16_t > perm_vtk(mesh::CellType type, int num_nodes)
Permutation array to map from VTK to DOLFINx node ordering.
Definition cells.cpp:577
std::vector< std::uint16_t > transpose(std::span< const std::uint16_t > map)
Compute the transpose of a re-ordering map.
Definition cells.cpp:654
std::tuple< std::vector< T >, std::array< std::size_t, 2 >, std::vector< std::int64_t >, std::vector< std::uint8_t >, std::vector< std::int64_t >, std::array< std::size_t, 2 > > vtk_mesh_from_space(const fem::FunctionSpace< T > &V)
Given a FunctionSpace, create a topology and geometry based on the dof coordinates.
Definition vtk_utils.h:192
std::pair< std::vector< std::int64_t >, std::array< std::size_t, 2 > > extract_vtk_connectivity(MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan< const std::int32_t, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents< std::size_t, 2 > > dofmap_x, mesh::CellType cell_type)
Extract the cell topology (connectivity) in VTK ordering for all cells the mesh. The 'topology' inclu...
Definition vtk_utils.cpp:23
CellType
Cell type identifier.
Definition cell_types.h:22
Top-level namespace.
Definition defines.h:12