DOLFINx 0.10.0.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>
32class FunctionSpace;
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 = md::mdspan<T, md::dextents<std::size_t, 2>>;
109 using cmdspan4_t = md::mdspan<T, md::dextents<std::size_t, 4>>;
110
111 // Tabulate basis functions at node reference coordinates
112 const std::array<std::size_t, 4> phi_shape
113 = cmap.tabulate_shape(0, Xshape[0]);
114 std::vector<T> phi_b(
115 std::reduce(phi_shape.begin(), phi_shape.end(), 1, std::multiplies{}));
116 cmdspan4_t phi_full(phi_b.data(), phi_shape);
117 cmap.tabulate(0, X, Xshape, phi_b);
118 auto phi = md::submdspan(phi_full, 0, md::full_extent, md::full_extent, 0);
119
120 // Loop over cells and tabulate dofs
121 auto map = topology->index_map(tdim);
122 assert(map);
123 const std::int32_t num_cells = map->size_local() + map->num_ghosts();
124 std::vector<T> x_b(scalar_dofs * gdim);
125 mdspan2_t x(x_b.data(), scalar_dofs, gdim);
126 std::vector<T> coordinate_dofs_b(num_dofs_g * gdim);
127 mdspan2_t coordinate_dofs(coordinate_dofs_b.data(), num_dofs_g, gdim);
128
129 std::vector<T> coords(num_nodes * 3, 0.0);
130 std::array<std::size_t, 2> cshape = {num_nodes, 3};
131 for (std::int32_t c = 0; c < num_cells; ++c)
132 {
133 // Extract cell geometry
134 for (std::size_t i = 0; i < dofmap_x.extent(1); ++i)
135 for (std::size_t j = 0; j < gdim; ++j)
136 coordinate_dofs(i, j) = x_g[3 * dofmap_x(c, i) + j];
137
138 // Tabulate dof coordinates on cell
139 cmap.push_forward(x, coordinate_dofs, phi);
140 apply_dof_transformation(x_b, std::span(cell_info.data(), cell_info.size()),
141 c, x.extent(1));
142
143 // Copy dof coordinates into vector
144 auto dofs = dofmap->cell_dofs(c);
145 for (std::size_t i = 0; i < dofs.size(); ++i)
146 {
147 std::int32_t dof = dofs[i];
148 for (std::size_t j = 0; j < gdim; ++j)
149 coords[3 * dof + j] = x(i, j);
150 }
151 }
152
153 // Original point IDs
154 std::vector<std::int64_t> x_id(num_nodes);
155 std::array<std::int64_t, 2> range = map_dofs->local_range();
156 std::int32_t size_local = range[1] - range[0];
157 std::iota(x_id.begin(), std::next(x_id.begin(), size_local), range[0]);
158 std::span ghosts = map_dofs->ghosts();
159 std::ranges::copy(ghosts, std::next(x_id.begin(), size_local));
160
161 // Ghosts
162 std::vector<std::uint8_t> id_ghost(num_nodes, 0);
163 std::fill(std::next(id_ghost.begin(), size_local), id_ghost.end(), 1);
164
165 return {std::move(coords), cshape, std::move(x_id), std::move(id_ghost)};
166}
167
168} // namespace impl
169
184template <typename T>
185std::tuple<std::vector<T>, std::array<std::size_t, 2>,
186 std::vector<std::int64_t>, std::vector<std::uint8_t>,
187 std::vector<std::int64_t>, std::array<std::size_t, 2>>
189{
190 auto mesh = V.mesh();
191 assert(mesh);
192 auto topology = mesh->topology();
193 assert(topology);
194 const int tdim = topology->dim();
195
196 assert(V.element());
197 if (V.element()->is_mixed())
198 throw std::runtime_error("Can't create VTK mesh from a mixed element");
199
200 const auto [x, xshape, x_id, x_ghost]
201 = impl::tabulate_lagrange_dof_coordinates(V);
202 auto map = topology->index_map(tdim);
203 const std::size_t num_cells = map->size_local() + map->num_ghosts();
204
205 // Create permutation from DOLFINx dof ordering to VTK
206 auto dofmap = V.dofmap();
207 assert(dofmap);
208 const int element_block_size = V.element()->block_size();
209 const std::uint32_t num_nodes
210 = V.element()->space_dimension() / element_block_size;
211 const std::vector<std::uint16_t> vtkmap = io::cells::transpose(
212 io::cells::perm_vtk(topology->cell_type(), num_nodes));
213
214 // Extract topology for all local cells as
215 // [v0_0, ...., v0_N0, v1_0, ...., v1_N1, ....]
216 std::array<std::size_t, 2> shape = {num_cells, num_nodes};
217 std::vector<std::int64_t> vtk_topology(shape[0] * shape[1]);
218 for (std::size_t c = 0; c < shape[0]; ++c)
219 {
220 auto dofs = dofmap->cell_dofs(c);
221 for (std::size_t i = 0; i < dofs.size(); ++i)
222 vtk_topology[c * shape[1] + i] = dofs[vtkmap[i]];
223 }
224
225 return {std::move(x),
226 xshape,
227 std::move(x_id),
228 std::move(x_ghost),
229 std::move(vtk_topology),
230 shape};
231}
232
247std::pair<std::vector<std::int64_t>, std::array<std::size_t, 2>>
249 md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> dofmap_x,
250 mesh::CellType cell_type);
251} // namespace io
252} // namespace dolfinx
Degree-of-freedom map representations and tools.
This class represents a finite element function space defined by a mesh, a finite element,...
Definition FunctionSpace.h:34
Geometry stores the geometry imposed on a mesh.
Definition Geometry.h:34
Finite element method functionality.
Definition assemble_expression_impl.h:23
@ standard
Standard.
Definition FiniteElement.h:27
FunctionSpace(U mesh, V element, W dofmap) -> FunctionSpace< typename std::remove_cvref< typename U::element_type >::type::geometry_type::value_type >
Type deduction.
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:516
std::vector< std::uint16_t > transpose(std::span< const std::uint16_t > map)
Compute the transpose of a re-ordering map.
Definition cells.cpp:670
Support for file IO.
Definition cells.h:119
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:188
std::pair< std::vector< std::int64_t >, std::array< std::size_t, 2 > > extract_vtk_connectivity(md::mdspan< const std::int32_t, md::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
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
CellType
Cell type identifier.
Definition cell_types.h:22
Top-level namespace.
Definition defines.h:12