GCC Code Coverage Report

Line	Branch	Exec	Source
1			// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2			// vi: set et ts=4 sw=4 sts=4:
3			//
4			// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5			// SPDX-License-Identifier: GPL-3.0-or-later
6			//
7			/*!
8			* \file
9			* \ingroup InputOutput
10			* \brief A data handle for communicating grid data for VTK grids
11			*/
12			#ifndef DUMUX_VTK_GRID_DATA_HANDLE_HH
13			#define DUMUX_VTK_GRID_DATA_HANDLE_HH
14
15			#include <memory>
16			#include <algorithm>
17			#include <map>
18			#include <numeric>
19
20			#include <dune/common/parallel/communication.hh>
21			#include <dune/geometry/dimension.hh>
22			#include <dune/grid/common/partitionset.hh>
23			#include <dune/grid/common/datahandleif.hh>
24
25			namespace Dumux {
26
27			/*!
28			* \ingroup InputOutput
29			* \brief A data handle for communicating grid data for VTK grids
30			*/
31			template<class Grid, class GridFactory, class Data>
32			struct VtkGridDataHandle
33			: public Dune::CommDataHandleIF<VtkGridDataHandle<Grid, GridFactory, Data>, typename Data::value_type>
34			{
35			using GridView = typename Grid::LevelGridView;
36
37		✗	VtkGridDataHandle(const Grid& grid, const GridFactory& gridFactory, VTKReader::Data& cellData, VTKReader::Data& pointData)
38			: gridView_(grid.levelGridView(0))
39			, idSet_(grid.localIdSet())
40			, userCellData_(cellData)
41		✗	, userPointData_(pointData)
42			{
43			// For the following to work we assume a sorted map of keys to values in the user data.
44			// This is not guaranteed by the VTKReader, so we need to sort the data first.
45		✗	for (const auto& [key, data] : userCellData_)
46		✗	cellData_[key] = std::move(userCellData_[key]);
47		✗	for (const auto& [key, data] : userPointData_)
48		✗	pointData_[key] = std::move(userPointData_[key]);
49
50			// assume all data is on rank 0 (see grid manager)
51			// First broadcast how many keys we have
52		✗	std::array<std::size_t, 2> numKeys{{ cellData_.size(), pointData_.size() }};
53		✗	Dune::MPIHelper::getCommunication().broadcast(numKeys.data(), 2, 0);
54
55			// then broadcast the length of the individual key strings
56			// and the number of component associated with each key (e.g. vector/tensor fields)
57		✗	std::vector<std::size_t> keyLengthAndComponents(2*(numKeys[0] + numKeys[1]), 0);
58			{
59		✗	int n = 0;
60
61			// key length
62		✗	for (const auto& [key, data] : cellData_)
63		✗	keyLengthAndComponents[n++] = key.size();
64		✗	for (const auto& [key, data] : pointData_)
65		✗	keyLengthAndComponents[n++] = key.size();
66
67			// number of components
68		✗	for (const auto& [key, data] : cellData_)
69		✗	keyLengthAndComponents[n++] = gridView_.size(0) > 0 ? data.size()/gridView_.size(0) : 0;
70		✗	for (const auto& [key, data] : pointData_)
71		✗	keyLengthAndComponents[n++] = gridView_.size(Grid::dimension) > 0 ? data.size()/gridView_.size(Grid::dimension) : 0;
72
73			// entries only exist on rank 0 and the data containers are empty on other ranks
74		✗	assert((Dune::MPIHelper::instance().rank() == 0) == (n == keyLengthAndComponents.size()));
75
76		✗	Dune::MPIHelper::getCommunication().broadcast(keyLengthAndComponents.data(), keyLengthAndComponents.size(), 0);
77			}
78
79			// save the number of components for each cell and point data array
80		✗	const auto begin = keyLengthAndComponents.begin() + numKeys[0] + numKeys[1];
81		✗	cellDataComponents_.assign(begin, begin + numKeys[0]);
82		✗	pointDataComponents_.assign(begin + numKeys[0], keyLengthAndComponents.end());
83		✗	numCellDataPerElement_ = std::accumulate(cellDataComponents_.begin(), cellDataComponents_.end(), 0UL);
84		✗	numPointDataPerVertex_ = std::accumulate(pointDataComponents_.begin(), pointDataComponents_.end(), 0UL);
85
86			// then broadcast the actual keys
87		✗	std::string keys; keys.resize(std::accumulate(keyLengthAndComponents.begin(), begin, 0UL));
88			{
89		✗	int n = 0;
90		✗	for (const auto& [key, data] : cellData_)
91		✗	for (const auto& c : key)
92		✗	keys[n++] = c;
93		✗	for (const auto& [key, data] : pointData_)
94		✗	for (const auto& c : key)
95		✗	keys[n++] = c;
96
97			// entries only exist on rank 0 and the data containers are empty on other ranks
98		✗	assert((Dune::MPIHelper::instance().rank() == 0) == (n == keys.size()));
99
100		✗	Dune::MPIHelper::getCommunication().broadcast(keys.data(), keys.size(), 0);
101			}
102
103			// create the entries in the cellData and pointData maps on all processes
104			std::size_t offset = 0;
105		✗	for (int keyIdx = 0; keyIdx < numKeys[0]; ++keyIdx)
106			{
107		✗	if (std::string key{ keys, offset, keyLengthAndComponents[keyIdx] }; cellData_.count(key) == 0)
108		✗	cellData_[key] = Data{};
109
110		✗	offset += keyLengthAndComponents[keyIdx];
111			}
112		✗	for (int keyIdx = numKeys[0]; keyIdx < numKeys[0] + numKeys[1]; ++keyIdx)
113			{
114		✗	if (std::string key{ keys, offset, keyLengthAndComponents[keyIdx] }; pointData_.count(key) == 0)
115		✗	pointData_[key] = Data{};
116
117		✗	offset += keyLengthAndComponents[keyIdx];
118			}
119
120			// write data into an id map
121		✗	for (const auto& element : elements(gridView_, Dune::Partitions::interior))
122			{
123		✗	data_[idSet_.id(element)].resize(numCellDataPerElement_);
124
125		✗	int n = 0, l = 0;
126		✗	for (const auto& [key, data] : cellData_)
127			{
128		✗	const auto nComp = cellDataComponents_[l++];
129		✗	for (int k = 0; k < nComp; ++k)
130		✗	std::swap(cellData_[key][k + nComp*gridFactory.insertionIndex(element)], data_[idSet_.id(element)][n++]);
131			}
132
133		✗	assert(n == numCellDataPerElement_);
134			}
135
136		✗	for (const auto& vertex : vertices(gridView_))
137			{
138		✗	data_[idSet_.id(vertex)].resize(numPointDataPerVertex_);
139
140		✗	int n = 0, l = 0;
141		✗	for (const auto& [key, data] : pointData_)
142			{
143		✗	const auto nComp = pointDataComponents_[l++];
144		✗	for (int k = 0; k < nComp; ++k)
145		✗	std::swap(pointData_[key][k + nComp*gridFactory.insertionIndex(vertex)], data_[idSet_.id(vertex)][n++]);
146			}
147
148		✗	assert(n == numPointDataPerVertex_);
149			}
150		✗	}
151
152		✗	~VtkGridDataHandle()
153			{
154			// resize arrays and unpack communicated data
155		✗	const auto& indexSet = gridView_.indexSet();
156
157			{
158		✗	int n = 0;
159		✗	for (const auto& [key, data] : cellData_)
160		✗	cellData_[key].resize(indexSet.size(0)*cellDataComponents_[n++]);
161			}
162			{
163		✗	int n = 0;
164		✗	for (const auto& [key, data] : pointData_)
165		✗	pointData_[key].resize(indexSet.size(GridView::dimension)*pointDataComponents_[n++]);
166			}
167
168		✗	for (const auto& element : elements(gridView_))
169			{
170		✗	int n = 0, l = 0;
171		✗	for (const auto& [key, data] : cellData_)
172			{
173		✗	const auto nComp = cellDataComponents_[l++];
174		✗	for (int k = 0; k < nComp; ++k)
175		✗	std::swap(cellData_[key][k + nComp*indexSet.index(element)], data_[idSet_.id(element)][n++]);
176			}
177			}
178
179		✗	for (const auto& vertex : vertices(gridView_))
180			{
181		✗	int n = 0, l = 0;
182		✗	for (const auto& [key, data] : pointData_)
183			{
184		✗	const auto nComp = pointDataComponents_[l++];
185		✗	for (int k = 0; k < nComp; ++k)
186		✗	std::swap(pointData_[key][k + nComp*indexSet.index(vertex)], data_[idSet_.id(vertex)][n++]);
187			}
188			}
189
190			// move data back from sorted internal storage to user data
191		✗	for (const auto& [key, data] : cellData_)
192		✗	userCellData_[key] = std::move(cellData_[key]);
193		✗	for (const auto& [key, data] : pointData_)
194		✗	userPointData_[key] = std::move(pointData_[key]);
195		✗	}
196
197			Dune::CommDataHandleIF<VtkGridDataHandle<Grid, GridFactory, Data>, typename Data::value_type>& interface()
198			{ return *this; }
199
200		✗	bool contains (int dim, int codim) const
201		✗	{ return codim == 0 || codim == dim; }
202
203			//! returns true if size per entity of given dim and codim is a constant
204		✗	bool fixedSize(int dim, int codim) const
205		✗	{ return true; }
206
207			template<class Entity>
208		✗	std::size_t size (const Entity&) const
209			{
210			if constexpr (Entity::codimension == 0)
211		✗	return numCellDataPerElement_;
212			else if constexpr (Entity::codimension == GridView::dimension)
213		✗	return numPointDataPerVertex_;
214			else
215		✗	return 0;
216			}
217
218			template<class MessageBufferImp, class Entity>
219		✗	void gather (MessageBufferImp& buff, const Entity& e) const
220			{
221			if constexpr (Entity::codimension == 0)
222			{
223		✗	const auto& data = data_[idSet_.id(e)];
224		✗	for (int n = 0; n < numCellDataPerElement_; ++n)
225		✗	buff.write(data[n]);
226			}
227
228			if constexpr (Entity::codimension == GridView::dimension)
229			{
230		✗	const auto& data = data_[idSet_.id(e)];
231		✗	for (int n = 0; n < numPointDataPerVertex_; ++n)
232		✗	buff.write(data[n]);
233			}
234		✗	}
235
236			template<class MessageBufferImp, class Entity>
237		✗	void scatter (MessageBufferImp& buff, const Entity& e, std::size_t n)
238			{
239		✗	auto& data = data_[idSet_.id(e)];
240		✗	data.resize(n);
241
242			if constexpr (Entity::codimension == 0)
243		✗	for (int k = 0; k < numCellDataPerElement_; ++k)
244		✗	buff.read(data[k]);
245
246			if constexpr (Entity::codimension == GridView::dimension)
247		✗	for (int k = 0; k < numPointDataPerVertex_; ++k)
248		✗	buff.read(data[k]);
249		✗	}
250
251			private:
252			using IdSet = typename Grid::LocalIdSet;
253
254			const GridView gridView_;
255			const IdSet &idSet_;
256
257			VTKReader::Data& userCellData_;
258			VTKReader::Data& userPointData_;
259
260			std::map<std::string, VTKReader::Data::mapped_type> cellData_;
261			std::map<std::string, VTKReader::Data::mapped_type> pointData_;
262
263			std::vector<std::size_t> cellDataComponents_;
264			std::vector<std::size_t> pointDataComponents_;
265
266			std::size_t numCellDataPerElement_;
267			std::size_t numPointDataPerVertex_;
268
269			mutable std::map< typename IdSet::IdType, std::vector<typename Data::value_type> > data_;
270			};
271
272			} // namespace Dumux
273
274			#endif
275