Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Base/Data.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Generic data storage with different memory layouts
9 : : \details Generic data storage with different memory layouts. See also the
10 : : rationale discussed in the [design](layout.html) document.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef Data_h
14 : : #define Data_h
15 : :
16 : : #include <array>
17 : : #include <string>
18 : : #include <cstdint>
19 : : #include <vector>
20 : : #include <set>
21 : : #include <algorithm>
22 : :
23 : : #include "Types.hpp"
24 : : #include "Keywords.hpp"
25 : : #include "Exception.hpp"
26 : :
27 : : #include "NoWarning/pup_stl.hpp"
28 : :
29 : : namespace tk {
30 : :
31 : : //! Tags for selecting data layout policies
32 : : const uint8_t UnkEqComp = 0;
33 : : const uint8_t EqCompUnk = 1;
34 : :
35 : : //! Zero-runtime-cost data-layout wrappers with type-based compile-time dispatch
36 : : template< uint8_t Layout >
37 : : class Data {
38 : :
39 : : private:
40 : : //! \brief Inherit type of number of components from keyword 'ncomp', used
41 : : //! also for type of offset
42 : : using ncomp_t = kw::ncomp::info::expect::type;
43 : :
44 : : public:
45 : : //! Default constructor (required for Charm++ migration)
46 : 102442 : explicit Data() : m_vec(), m_nunk(), m_nprop() {}
47 : :
48 : : //! Constructor
49 : : //! \param[in] nu Number of unknowns to allocate memory for
50 : : //! \param[in] np Total number of properties, i.e., scalar variables or
51 : : //! components, per unknown
52 : 5662130 : explicit Data( ncomp_t nu, ncomp_t np ) :
53 : : m_vec( nu*np ),
54 : : m_nunk( nu ),
55 [ + - ]: 5662130 : m_nprop( np ) {}
56 : :
57 : : //! Const data access dispatch
58 : : //! \details Public interface to const-ref data access to a single real
59 : : //! value. Use it as Data(p,c,o), where p is the unknown index, c is
60 : : //! the component index specifying the scalar equation within a system of
61 : : //! equations, and o is the offset specifying the position at which the
62 : : //! system resides among other systems. Requirement: offset + component <
63 : : //! nprop, unknown < nunk, enforced with an assert in DEBUG mode, see also
64 : : //! the constructor.
65 : : //! \param[in] unknown Unknown index
66 : : //! \param[in] component Component index, i.e., position of a scalar within
67 : : //! a system
68 : : //! \param[in] offset System offset specifying the position of the system of
69 : : //! equations among other systems
70 : : //! \return Const reference to data of type tk::real
71 : : const tk::real&
72 :35148481364 : operator()( ncomp_t unknown, ncomp_t component, ncomp_t offset ) const
73 :35148481364 : { return access( unknown, component, offset, int2type< Layout >() ); }
74 : :
75 : : //! Non-const data access dispatch
76 : : //! \details Public interface to non-const-ref data access to a single real
77 : : //! value. Use it as Data(p,c,o), where p is the unknown index, c is
78 : : //! the component index specifying the scalar equation within a system of
79 : : //! equations, and o is the offset specifying the position at which the
80 : : //! system resides among other systems. Requirement: offset + component <
81 : : //! nprop, unknown < nunk, enforced with an assert in DEBUG mode, see also
82 : : //! the constructor.
83 : : //! \param[in] unknown Unknown index
84 : : //! \param[in] component Component index, i.e., position of a scalar within
85 : : //! a system
86 : : //! \param[in] offset System offset specifying the position of the system of
87 : : //! equations among other systems
88 : : //! \return Non-const reference to data of type tk::real
89 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
90 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
91 : : tk::real&
92 :21224106536 : operator()( ncomp_t unknown, ncomp_t component, ncomp_t offset ) {
93 :21224106536 : return const_cast< tk::real& >(
94 : : static_cast< const Data& >( *this ).
95 :21224106536 : operator()( unknown, component, offset ) );
96 : : }
97 : :
98 : : //! Const ptr to physical variable access dispatch
99 : : //! \details Public interface to the first half of a physical variable
100 : : //! access. cptr() and var() are two member functions intended to be used
101 : : //! together in case when component and offset would be expensive to
102 : : //! compute for data access via the function call operator, i.e., cptr()
103 : : //! can be used to pre-compute part of the address, which returns a
104 : : //! pointer and var() can be used to finish the data access using the
105 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
106 : : //! address known based on component and offset and intended to be used in
107 : : //! a setup phase. Then var() takes this partial address and finishes the
108 : : //! address calculation given the unknown id. Thus the following two data
109 : : //! accesses are equivalent (modulo constness):
110 : : //! * real& value = operator()( unk, comp, offs ); and
111 : : //! * const real* p = cptr( comp, offs ); and
112 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
113 : : //! Requirement: offset + component < nprop, enforced with an assert in
114 : : //! DEBUG mode, see also the constructor.
115 : : //! \param[in] component Component index, i.e., position of a scalar within
116 : : //! a system
117 : : //! \param[in] offset System offset specifying the position of the system of
118 : : //! equations among other systems
119 : : //! \return Pointer to data of type tk::real for use with var()
120 : : //! \see Example client code in Statistics::setupOrdinary() and
121 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
122 : : const tk::real*
123 : 151367071 : cptr( ncomp_t component, ncomp_t offset ) const
124 : 151367071 : { return cptr( component, offset, int2type< Layout >() ); }
125 : :
126 : : //! Const-ref data-access dispatch
127 : : //! \details Public interface to the second half of a physical variable
128 : : //! access. cptr() and var() are two member functions intended to be used
129 : : //! together in case when component and offset would be expensive to
130 : : //! compute for data access via the function call operator, i.e., cptr()
131 : : //! can be used to pre-compute part of the address, which returns a
132 : : //! pointer and var() can be used to finish the data access using the
133 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
134 : : //! address known based on component and offset and intended to be used in
135 : : //! a setup phase. Then var() takes this partial address and finishes the
136 : : //! address calculation given the unknown id. Thus the following two data
137 : : //! accesses are equivalent (modulo constness):
138 : : //! * real& value = operator()( unk, comp, offs ); and
139 : : //! * const real* p = cptr( comp, offs ); and
140 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
141 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
142 : : //! see also the constructor.
143 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
144 : : //! \param[in] unknown Unknown index
145 : : //! \return Const reference to data of type tk::real
146 : : //! \see Example client code in Statistics::setupOrdinary() and
147 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
148 : : const tk::real&
149 :24931069007 : var( const tk::real* pt, ncomp_t unknown ) const
150 :24931069007 : { return var( pt, unknown, int2type< Layout >() ); }
151 : :
152 : : //! Non-const-ref data-access dispatch
153 : : //! \details Public interface to the second half of a physical variable
154 : : //! access. cptr() and var() are two member functions intended to be used
155 : : //! together in case when component and offset would be expensive to
156 : : //! compute for data access via the function call operator, i.e., cptr()
157 : : //! can be used to pre-compute part of the address, which returns a
158 : : //! pointer and var() can be used to finish the data access using the
159 : : //! pointer returned by cptr(). In other words, cptr() returns part of the
160 : : //! address known based on component and offset and intended to be used in
161 : : //! a setup phase. Then var() takes this partial address and finishes the
162 : : //! address calculation given the unknown id. Thus the following two data
163 : : //! accesses are equivalent (modulo constness):
164 : : //! * real& value = operator()( unk, comp, offs ); and
165 : : //! * const real* p = cptr( comp, offs ); and
166 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
167 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
168 : : //! see also the constructor.
169 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
170 : : //! \param[in] unknown Unknown index
171 : : //! \return Non-const reference to data of type tk::real
172 : : //! \see Example client code in Statistics::setupOrdinary() and
173 : : //! Statistics::accumulateOrd() in Statistics/Statistics.C.
174 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
175 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
176 : : tk::real&
177 : 5403902689 : var( const tk::real* pt, ncomp_t unknown ) {
178 : 5403902689 : return const_cast< tk::real& >(
179 : 5403902689 : static_cast< const Data& >( *this ).var( pt, unknown ) );
180 : : }
181 : :
182 : : //! Access to number of unknowns
183 : : //! \return Number of unknowns
184 : 213617289 : ncomp_t nunk() const noexcept { return m_nunk; }
185 : :
186 : : //! Access to number of properties
187 : : //! \details This is the total number of scalar components per unknown
188 : : //! \return Number of propertes/unknown
189 : 1483687307 : ncomp_t nprop() const noexcept { return m_nprop; }
190 : :
191 : : //! Extract flat vector of all unknowns
192 : : //! \return Flat vector of reals
193 : : std::vector< tk::real >
194 : 677 : flat() const {
195 [ + - ]: 677 : std::vector< tk::real > w( m_nunk * m_nprop );
196 [ + + ]: 2708 : for (std::size_t j=0; j<m_nprop; ++j)
197 [ + + ]: 3754737 : for (std::size_t i=0; i<m_nunk; ++i)
198 [ + - ]: 3752706 : w[i*m_nprop+j] = operator()( i, j, 0 );
199 : 677 : return w;
200 : : }
201 : :
202 : : //! Assign from flat vector
203 : : //! \param[in] rhs Flat vector to assign from
204 : : //! \note This is the opposite of flat().
205 : : void operator=( const std::vector< tk::real >& rhs ) {
206 : : Assert( rhs.size() == m_nunk * m_nprop, "Size mismatch" );
207 : : for (std::size_t j=0; j<m_nprop; ++j)
208 : : for (std::size_t i=0; i<m_nunk; ++i)
209 : : operator()( i, j, 0 ) = rhs[i*m_nprop+j];
210 : : }
211 : :
212 : : //! Assign from array(3) of vectors
213 : : //! \param[in] rhs Array(3) of vectors to assign from
214 : : void operator=( const std::array< std::vector< tk::real >, 3 >& rhs ) {
215 : : Assert( m_nprop == 3, "Size mismatch" );
216 : : Assert( rhs[0].size() == m_nunk, "Size mismatch" );
217 : : Assert( rhs[1].size() == m_nunk, "Size mismatch" );
218 : : Assert( rhs[2].size() == m_nunk, "Size mismatch" );
219 : : for (std::size_t j=0; j<3; ++j)
220 : : for (std::size_t i=0; i<m_nunk; ++i)
221 : : operator()( i, j, 0 ) = rhs[j][i];
222 : : }
223 : :
224 : : //! Extract vector of unknowns given component and offset
225 : : //! \details Requirement: offset + component < nprop, enforced with an
226 : : //! assert in DEBUG mode, see also the constructor.
227 : : //! \param[in] component Component index, i.e., position of a scalar within
228 : : //! a system
229 : : //! \param[in] offset System offset specifying the position of the system of
230 : : //! equations among other systems
231 : : //! \return A vector of unknowns given by component at offset (length:
232 : : //! nunk(), i.e., the first constructor argument)
233 : : std::vector< tk::real >
234 : 134932 : extract( ncomp_t component, ncomp_t offset ) const {
235 [ + - ]: 134932 : std::vector< tk::real > w( m_nunk );
236 [ + + ]: 29647749 : for (ncomp_t i=0; i<m_nunk; ++i)
237 [ + - ]: 29512817 : w[i] = operator()( i, component, offset );
238 : 134932 : return w;
239 : : }
240 : :
241 : : //! Extract (a copy of) all components for an unknown
242 : : //! \details Requirement: unknown < nunk, enforced with an assert in DEBUG
243 : : //! mode, see also the constructor.
244 : : //! \param[in] unknown Index of unknown
245 : : //! \return A vector of components for a single unknown (length: nprop,
246 : : //! i.e., the second constructor argument)
247 : : std::vector< tk::real >
248 : 137515387 : extract( ncomp_t unknown ) const {
249 [ + - ]: 137515387 : std::vector< tk::real > w( m_nprop );
250 [ + + ][ + - ]: 745522118 : for (ncomp_t i=0; i<m_nprop; ++i) w[i] = operator()( unknown, i, 0 );
251 : 137515387 : return w;
252 : : }
253 : :
254 : : //! Extract all components for unknown
255 : : //! \details Requirement: unknown < nunk, enforced with an assert in DEBUG
256 : : //! mode, see also the constructor.
257 : : //! \param[in] unknown Index of unknown
258 : : //! \return A vector of components for a single unknown (length: nprop,
259 : : //! i.e., the second constructor argument)
260 : : //! \note This is simply an alias for extract( unknown )
261 : : std::vector< tk::real >
262 : 137507259 : operator[]( ncomp_t unknown ) const { return extract( unknown ); }
263 : :
264 : : //! Extract (a copy of) four values of unknowns
265 : : //! \details Requirement: offset + component < nprop, [A,B,C,D] < nunk,
266 : : //! enforced with an assert in DEBUG mode, see also the constructor.
267 : : //! \param[in] component Component index, i.e., position of a scalar within
268 : : //! a system
269 : : //! \param[in] offset System offset specifying the position of the system of
270 : : //! equations among other systems
271 : : //! \param[in] A Index of 1st unknown
272 : : //! \param[in] B Index of 2nd unknown
273 : : //! \param[in] C Index of 3rd unknown
274 : : //! \param[in] D Index of 4th unknown
275 : : //! \return Array of the four values of component at offset
276 : : std::array< tk::real, 4 >
277 : 89763342 : extract( ncomp_t component, ncomp_t offset,
278 : : ncomp_t A, ncomp_t B, ncomp_t C, ncomp_t D ) const
279 : : {
280 : 89763342 : auto p = cptr( component, offset );
281 : 89763342 : return {{ var(p,A), var(p,B), var(p,C), var(p,D) }};
282 : : }
283 : :
284 : : //! Extract (a copy of) four values of unknowns
285 : : //! \details Requirement: offset + component < nprop, for all N[i] < nunk,
286 : : //! enforced with an assert in DEBUG mode, see also the constructor.
287 : : //! \param[in] component Component index, i.e., position of a scalar within
288 : : //! a system
289 : : //! \param[in] offset System offset specifying the position of the system of
290 : : //! equations among other systems
291 : : //! \param[in] N Indices of the 4 unknowns
292 : : //! \return Array of the four values of component at offset
293 : : std::array< tk::real, 4 >
294 : 89763342 : extract( ncomp_t component, ncomp_t offset,
295 : : const std::array< ncomp_t, 4 >& N ) const
296 : : {
297 : 89763342 : return extract( component, offset, N[0], N[1], N[2], N[3] );
298 : : }
299 : :
300 : : //! Extract (a copy of) three values of unknowns
301 : : //! \details Requirement: offset + component < nprop, [A,B,C] < nunk,
302 : : //! enforced with an assert in DEBUG mode, see also the constructor.
303 : : //! \param[in] component Component index, i.e., position of a scalar within
304 : : //! a system
305 : : //! \param[in] offset System offset specifying the position of the system of
306 : : //! equations among other systems
307 : : //! \param[in] A Index of 1st unknown
308 : : //! \param[in] B Index of 2nd unknown
309 : : //! \param[in] C Index of 3rd unknown
310 : : //! \return Array of the four values of component at offset
311 : : std::array< tk::real, 3 >
312 : 91500 : extract( ncomp_t component, ncomp_t offset,
313 : : ncomp_t A, ncomp_t B, ncomp_t C ) const
314 : : {
315 : 91500 : auto p = cptr( component, offset );
316 : 91500 : return {{ var(p,A), var(p,B), var(p,C) }};
317 : : }
318 : :
319 : : //! Extract (a copy of) three values of unknowns
320 : : //! \details Requirement: offset + component < nprop, for all N[i] < nunk,
321 : : //! enforced with an assert in DEBUG mode, see also the constructor.
322 : : //! \param[in] component Component index, i.e., position of a scalar within
323 : : //! a system
324 : : //! \param[in] offset System offset specifying the position of the system of
325 : : //! equations among other systems
326 : : //! \param[in] N Indices of the 3 unknowns
327 : : //! \return Array of the three values of component at offset
328 : : std::array< tk::real, 3 >
329 : 91500 : extract( ncomp_t component, ncomp_t offset,
330 : : const std::array< ncomp_t, 3 >& N ) const
331 : : {
332 : 91500 : return extract( component, offset, N[0], N[1], N[2] );
333 : : }
334 : :
335 : : //! Const-ref accessor to underlying raw data as a std::vector
336 : : //! \return Constant reference to underlying raw data
337 : 1165633 : const std::vector< tk::real >& vec() const { return m_vec; }
338 : :
339 : : //! Non-const-ref accessor to underlying raw data as a std::vector
340 : : //! \return Non-constant reference to underlying raw data
341 : 0 : std::vector< tk::real >& vec() { return m_vec; }
342 : :
343 : : //! Compound operator-=
344 : : //! \param[in] rhs Data object to subtract
345 : : //! \return Reference to ourselves after subtraction
346 : 499 : Data< Layout >& operator-= ( const Data< Layout >& rhs ) {
347 [ - + ][ - - ]: 499 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
348 [ - + ][ - - ]: 499 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
349 : 499 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
350 : : m_vec.cbegin(), m_vec.begin(),
351 : 884765 : []( tk::real s, tk::real d ){ return d-s; } );
352 : 499 : return *this;
353 : : }
354 : : //! Operator -
355 : : //! \param[in] rhs Data object to subtract
356 : : //! \return Copy of Data object after rhs has been subtracted
357 : : //! \details Implemented in terms of compound operator-=
358 : 499 : Data< Layout > operator- ( const Data< Layout >& rhs )
359 [ + - ][ + - ]: 998 : const { return Data< Layout >( *this ) -= rhs; }
360 : :
361 : : //! Compound operator+=
362 : : //! \param[in] rhs Data object to add
363 : : //! \return Reference to ourselves after addition
364 : 101055 : Data< Layout >& operator+= ( const Data< Layout >& rhs ) {
365 [ - + ][ - - ]: 101055 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
366 [ - + ][ - - ]: 101055 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
367 : 101055 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
368 : : m_vec.cbegin(), m_vec.begin(),
369 : 40242306 : []( tk::real s, tk::real d ){ return d+s; } );
370 : 101055 : return *this;
371 : : }
372 : : //! Operator +
373 : : //! \param[in] rhs Data object to add
374 : : //! \return Copy of Data object after rhs has been multiplied with
375 : : //! \details Implemented in terms of compound operator+=
376 : 101055 : Data< Layout > operator+ ( const Data< Layout >& rhs )
377 [ + - ][ + - ]: 202110 : const { return Data< Layout >( *this ) += rhs; }
378 : :
379 : : //! Compound operator*= multiplying by another Data object item by item
380 : : //! \param[in] rhs Data object to multiply with
381 : : //! \return Reference to ourselves after multiplication
382 : : Data< Layout >& operator*= ( const Data< Layout >& rhs ) {
383 : : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
384 : : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
385 : : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
386 : : m_vec.cbegin(), m_vec.begin(),
387 : : []( tk::real s, tk::real d ){ return d*s; } );
388 : : return *this;
389 : : }
390 : : //! Operator * multiplying by another Data object item by item
391 : : //! \param[in] rhs Data object to multiply with
392 : : //! \return Copy of Data object after rhs has been multiplied with
393 : : //! \details Implemented in terms of compound operator*=
394 : : Data< Layout > operator* ( const Data< Layout >& rhs )
395 : : const { return Data< Layout >( *this ) *= rhs; }
396 : :
397 : : //! Compound operator*= multiplying all items by a scalar
398 : : //! \param[in] rhs Scalar to multiply with
399 : : //! \return Reference to ourselves after multiplication
400 : 46206 : Data< Layout >& operator*= ( tk::real rhs ) {
401 : : // cppcheck-suppress useStlAlgorithm
402 [ + + ]: 25709661 : for (auto& v : m_vec) v *= rhs;
403 : 46206 : return *this;
404 : : }
405 : : //! Operator * multiplying all items by a scalar
406 : : //! \param[in] rhs Scalar to multiply with
407 : : //! \return Copy of Data object after rhs has been multiplied with
408 : : //! \details Implemented in terms of compound operator*=
409 : : Data< Layout > operator* ( tk::real rhs )
410 : : const { return Data< Layout >( *this ) *= rhs; }
411 : :
412 : : //! Compound operator/=
413 : : //! \param[in] rhs Data object to divide by
414 : : //! \return Reference to ourselves after division
415 : 36566 : Data< Layout >& operator/= ( const Data< Layout >& rhs ) {
416 [ - + ][ - - ]: 36566 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
417 [ - + ][ - - ]: 36566 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
418 : 36566 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
419 : : m_vec.cbegin(), m_vec.begin(),
420 : 9719234 : []( tk::real s, tk::real d ){ return d/s; } );
421 : 36566 : return *this;
422 : : }
423 : : //! Operator /
424 : : //! \param[in] rhs Data object to divide by
425 : : //! \return Copy of Data object after rhs has been divided by
426 : : //! \details Implemented in terms of compound operator/=
427 : 36566 : Data< Layout > operator/ ( const Data< Layout >& rhs )
428 [ + - ][ + - ]: 73132 : const { return Data< Layout >( *this ) /= rhs; }
429 : :
430 : : //! Compound operator/= dividing all items by a scalar
431 : : //! \param[in] rhs Scalar to divide with
432 : : //! \return Reference to ourselves after division
433 : : Data< Layout >& operator/= ( tk::real rhs ) {
434 : : // cppcheck-suppress useStlAlgorithm
435 : : for (auto& v : m_vec) v /= rhs;
436 : : return *this;
437 : : }
438 : : //! Operator / dividing all items by a scalar
439 : : //! \param[in] rhs Scalar to divide with
440 : : //! \return Copy of Data object after rhs has been divided by
441 : : //! \details Implemented in terms of compound operator/=
442 : : Data< Layout > operator/ ( tk::real rhs )
443 : : const { return Data< Layout >( *this ) /= rhs; }
444 : :
445 : : //! Add new unknown at the end of the container
446 : : //! \param[in] prop Vector of properties to initialize the new unknown with
447 : 479971 : void push_back( const std::vector< tk::real >& prop )
448 : 479971 : { return push_back( prop, int2type< Layout >() ); }
449 : :
450 : : //! Resize data store to contain 'count' elements
451 : : //! \param[in] count Resize store to contain count * nprop elements
452 : : //! \param[in] value Value to initialize new data with (default: 0.0)
453 : : //! \note This works for both shrinking and enlarging, as this simply
454 : : //! translates to std::vector::resize(). Note that count changes, nprop
455 : : //! does not, see the private overload resize().
456 : 15016 : void resize( std::size_t count, tk::real value = 0.0 )
457 : 15016 : { resize( count, value, int2type< Layout >() ); }
458 : :
459 : : //! Remove a number of unknowns
460 : : //! \param[in] unknown Set of indices of unknowns to remove
461 : 380 : void rm( const std::set< ncomp_t >& unknown ) {
462 : 107580 : auto remove = [ &unknown ]( std::size_t i ) -> bool {
463 [ + - ][ - + ]: 53600 : if (unknown.find(i) != end(unknown)) return true;
464 : 53600 : return false;
465 : : };
466 : 380 : std::size_t last = 0;
467 [ + + ]: 53980 : for(std::size_t i=0; i<m_nunk; ++i, ++last) {
468 [ + - ][ - + ]: 53600 : while( remove(i) ) ++i;
469 [ - + ]: 53600 : if (i >= m_nunk) break;
470 [ + + ]: 266080 : for (ncomp_t p = 0; p<m_nprop; ++p)
471 : 212480 : m_vec[ last*m_nprop+p ] = m_vec[ i*m_nprop+p ];
472 : : }
473 [ + - ]: 380 : m_vec.resize( last*m_nprop );
474 : 380 : m_nunk -= unknown.size();
475 : 380 : }
476 : :
477 : : //! Fill vector of unknowns with the same value
478 : : //! \details Requirement: offset + component < nprop, enforced with an
479 : : //! assert in DEBUG mode, see also the constructor.
480 : : //! \param[in] component Component index, i.e., position of a scalar within
481 : : //! a system
482 : : //! \param[in] offset System offset specifying the position of the system of
483 : : //! equations among other systems
484 : : //! \param[in] value Value to fill vector of unknowns with
485 : 146550 : inline void fill( ncomp_t component, ncomp_t offset, tk::real value ) {
486 : 146550 : auto p = cptr( component, offset );
487 [ + + ]: 26144805 : for (ncomp_t i=0; i<m_nunk; ++i) var(p,i) = value;
488 : 146550 : }
489 : :
490 : : //! Fill full data storage with value
491 : : //! \param[in] value Value to fill data with
492 : 286172 : void fill( tk::real value )
493 : 286172 : { std::fill( begin(m_vec), end(m_vec), value ); }
494 : :
495 : : //! Check if vector of unknowns is empty
496 : 611300633 : bool empty() const noexcept { return m_vec.empty(); }
497 : :
498 : : //! Layout name dispatch
499 : : //! \return The name of the data layout used
500 : 307 : static std::string layout() { return layout( int2type< Layout >() ); }
501 : :
502 : : /** @name Pack/Unpack: Serialize Data object for Charm++ */
503 : : ///@{
504 : : //! \brief Pack/Unpack serialize member function
505 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
506 : 234177 : void pup( PUP::er &p ) {
507 : 234177 : p | m_vec;
508 : 234177 : p | m_nunk;
509 : 234177 : p | m_nprop;
510 : 234177 : }
511 : : //! \brief Pack/Unpack serialize operator|
512 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
513 : : //! \param[in,out] d DataLyaout object reference
514 : 234177 : friend void operator|( PUP::er& p, Data& d ) { d.pup(p); }
515 : : //@}
516 : :
517 : : private:
518 : : //! Transform a compile-time uint8_t into a type, used for dispatch
519 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
520 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
521 : : template< uint8_t m > struct int2type { enum { value = m }; };
522 : :
523 : : //! Overloads for the various const data accesses
524 : : //! \details Requirement: offset + component < nprop, unknown < nunk,
525 : : //! enforced with an assert in DEBUG mode, see also the constructor.
526 : : //! \param[in] unknown Unknown index
527 : : //! \param[in] component Component index, i.e., position of a scalar within
528 : : //! a system
529 : : //! \param[in] offset System offset specifying the position of the system of
530 : : //! equations among other systems
531 : : //! \return Const reference to data of type tk::real
532 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
533 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
534 : : const tk::real&
535 :35148481364 : access( ncomp_t unknown, ncomp_t component, ncomp_t offset,
536 : : int2type< UnkEqComp > ) const
537 : : {
538 [ - + ][ - - ]:35148481364 : Assert( offset + component < m_nprop, "Out-of-bounds access: offset + "
[ - - ][ - - ]
539 : : "component < number of properties" );
540 [ - + ][ - - ]:35148481364 : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
[ - - ][ - - ]
541 : : "unknowns" );
542 :35148481364 : return m_vec[ unknown*m_nprop + offset + component ];
543 : : }
544 : : const tk::real&
545 : : access( ncomp_t unknown, ncomp_t component, ncomp_t offset,
546 : : int2type< EqCompUnk > ) const
547 : : {
548 : : Assert( offset + component < m_nprop, "Out-of-bounds access: offset + "
549 : : "component < number of properties" );
550 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
551 : : "unknowns" );
552 : : return m_vec[ (offset+component)*m_nunk + unknown ];
553 : : }
554 : :
555 : : // Overloads for the various const ptr to physical variable accesses
556 : : //! \details Requirement: offset + component < nprop, unknown < nunk,
557 : : //! enforced with an assert in DEBUG mode, see also the constructor.
558 : : //! \param[in] component Component index, i.e., position of a scalar within
559 : : //! a system
560 : : //! \param[in] offset System offset specifying the position of the system of
561 : : //! equations among other systems
562 : : //! \return Pointer to data of type tk::real for use with var()
563 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
564 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
565 : : const tk::real*
566 : 151367071 : cptr( ncomp_t component, ncomp_t offset, int2type< UnkEqComp > ) const {
567 [ - + ][ - - ]: 151367071 : Assert( offset + component < m_nprop, "Out-of-bounds access: offset + "
[ - - ][ - - ]
568 : : "component < number of properties" );
569 : 151367071 : return m_vec.data() + component + offset;
570 : : }
571 : : const tk::real*
572 : : cptr( ncomp_t component, ncomp_t offset, int2type< EqCompUnk > ) const {
573 : : Assert( offset + component < m_nprop, "Out-of-bounds access: offset + "
574 : : "component < number of properties" );
575 : : return m_vec.data() + (offset+component)*m_nunk;
576 : : }
577 : :
578 : : // Overloads for the various const physical variable accesses
579 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
580 : : //! see also the constructor.
581 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
582 : : //! \param[in] unknown Unknown index
583 : : //! \return Const reference to data of type tk::real
584 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
585 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
586 : : inline const tk::real&
587 :24931069007 : var( const tk::real* const pt, ncomp_t unknown, int2type< UnkEqComp > )
588 : : const {
589 [ - + ][ - - ]:24931069007 : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
[ - - ][ - - ]
590 : : "unknowns" );
591 :24931069007 : return *(pt + unknown*m_nprop);
592 : : }
593 : : inline const tk::real&
594 : : var( const tk::real* const pt, ncomp_t unknown, int2type< EqCompUnk > )
595 : : const {
596 : : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
597 : : "unknowns" );
598 : : return *(pt + unknown);
599 : : }
600 : :
601 : : //! Add new unknown
602 : : //! \param[in] prop Vector of properties to initialize the new unknown with
603 : : //! \note Only the UnkEqComp overload is provided as this operation would be
604 : : //! too inefficient with the EqCompUnk data layout.
605 : 479971 : void push_back( const std::vector< tk::real >& prop, int2type< UnkEqComp > )
606 : : {
607 [ - + ][ - - ]: 479971 : Assert( prop.size() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
608 : 479971 : m_vec.resize( (m_nunk+1) * m_nprop );
609 : 479971 : ncomp_t u = m_nunk;
610 : 479971 : ++m_nunk;
611 [ + + ]: 2879826 : for (ncomp_t i=0; i<m_nprop; ++i) operator()( u, i, 0 ) = prop[i];
612 : 479971 : }
613 : :
614 : : void push_back( const std::vector< tk::real >&, int2type< EqCompUnk > )
615 : : { Throw( "Not implented. It would be inefficient" ); }
616 : :
617 : : //! Resize data store to contain 'count' elements
618 : : //! \param[in] count Resize store to contain 'count' elements
619 : : //! \param[in] value Value to initialize new data with
620 : : //! \note Only the UnkEqComp overload is provided as this operation would be
621 : : //! too inefficient with the EqCompUnk data layout.
622 : : //! \note This works for both shrinking and enlarging, as this simply
623 : : //! translates to std::vector::resize().
624 : 15016 : void resize( std::size_t count, tk::real value, int2type< UnkEqComp > ) {
625 : 15016 : m_vec.resize( count * m_nprop, value );
626 : 15016 : m_nunk = count;
627 : 15016 : }
628 : :
629 : : void resize( std::size_t, tk::real, int2type< EqCompUnk > ) {
630 : : Throw( "Not implemented. It would be inefficient" );
631 : : }
632 : :
633 : : // Overloads for the name-queries of data lauouts
634 : : //! \return The name of the data layout used
635 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
636 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
637 : 307 : static std::string layout( int2type< UnkEqComp > )
638 [ + - ]: 307 : { return "unknown-major"; }
639 : : static std::string layout( int2type< EqCompUnk > )
640 : : { return "equation-major"; }
641 : :
642 : : std::vector< tk::real > m_vec; //!< Data pointer
643 : : ncomp_t m_nunk; //!< Number of unknowns
644 : : ncomp_t m_nprop; //!< Number of properties/unknown
645 : : };
646 : :
647 : : //! Operator * multiplying all items by a scalar from the left
648 : : //! \param[in] lhs Scalar to multiply with
649 : : //! \param[in] rhs Date object to multiply
650 : : //! \return New Data object with all items multipled with lhs
651 : : template< uint8_t Layout >
652 : 46206 : Data< Layout > operator* ( tk::real lhs, const Data< Layout >& rhs ) {
653 [ + - ]: 92412 : return Data< Layout >( rhs ) *= lhs;
654 : : }
655 : :
656 : : //! Operator min between two Data objects
657 : : //! \param[in] a 1st Data object
658 : : //! \param[in] b 2nd Data object
659 : : //! \return New Data object containing the minimum of all values for each
660 : : //! value in _a_ and _b_
661 : : //! \note The Data objects _a_ and _b_ must have the same number of
662 : : //! unknowns and properties.
663 : : //! \note As opposed to std::min, this function creates and returns a new object
664 : : //! instead of returning a reference to one of the operands.
665 : : template< uint8_t Layout >
666 : : Data< Layout > min( const Data< Layout >& a, const Data< Layout >& b ) {
667 : : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
668 : : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
669 : : Data< Layout > r( a.nunk(), a.nprop() );
670 : : std::transform( a.vec().cbegin(), a.vec().cend(),
671 : : b.vec().cbegin(), r.vec().begin(),
672 : : []( tk::real s, tk::real d ){ return std::min(s,d); } );
673 : :
674 : : return r;
675 : : }
676 : :
677 : : //! Operator max between two Data objects
678 : : //! \param[in] a 1st Data object
679 : : //! \param[in] b 2nd Data object
680 : : //! \return New Data object containing the maximum of all values for each
681 : : //! value in _a_ and _b_
682 : : //! \note The Data objects _a_ and _b_ must have the same number of
683 : : //! unknowns and properties.
684 : : //! \note As opposed to std::max, this function creates and returns a new object
685 : : //! instead of returning a reference to one of the operands.
686 : : template< uint8_t Layout >
687 : : Data< Layout > max( const Data< Layout >& a, const Data< Layout >& b ) {
688 : : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
689 : : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
690 : : Data< Layout > r( a.nunk(), a.nprop() );
691 : : std::transform( a.vec().cbegin(), a.vec().cend(),
692 : : b.vec().cbegin(), r.vec().begin(),
693 : : []( tk::real s, tk::real d ){ return std::max(s,d); } );
694 : : return r;
695 : : }
696 : :
697 : : //! Operator == between two Data objects
698 : : //! \param[in] lhs Data object to compare
699 : : //! \param[in] rhs Data object to compare
700 : : //! \return True if all entries are equal up to epsilon
701 : : template< uint8_t Layout >
702 : : bool operator== ( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
703 : : Assert( rhs.nunk() == lhs.nunk(), "Incorrect number of unknowns" );
704 : : Assert( rhs.nprop() == lhs.nprop(), "Incorrect number of properties" );
705 : : auto l = lhs.vec().cbegin();
706 : : auto r = rhs.vec().cbegin();
707 : : while (l != lhs.vec().cend()) {
708 : : if (std::abs(*l - *r) > std::numeric_limits< tk::real >::epsilon())
709 : : return false;
710 : : ++l; ++r;
711 : : }
712 : : return true;
713 : : }
714 : :
715 : : //! Operator != between two Data objects
716 : : //! \param[in] lhs Data object to compare
717 : : //! \param[in] rhs Data object to compare
718 : : //! \return True if all entries are unequal up to epsilon
719 : : template< uint8_t Layout >
720 : : bool operator!= ( const Data< Layout >& lhs, const Data< Layout >& rhs )
721 : : { return !(lhs == rhs); }
722 : :
723 : : //! Compute the maximum difference between the elements of two Data objects
724 : : //! \param[in] lhs 1st Data object
725 : : //! \param[in] rhs 2nd Data object
726 : : //! \return The index, i.e., the raw position, of and the largest absolute value
727 : : //! of the difference between all corresponding elements of _lhs_ and _rhs_.
728 : : //! \details The position returned is the position in the underlying raw data
729 : : //! structure, independent of components, offsets, etc. If lhs == rhs with
730 : : //! precision std::numeric_limits< tk::real >::epsilon(), a pair of (0,0.0)
731 : : //! is returned.
732 : : //! \note The Data objects _lhs_ and _rhs_ must have the same number of
733 : : //! unknowns and properties.
734 : : template< uint8_t Layout >
735 : : std::pair< std::size_t, tk::real >
736 : 499 : maxdiff( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
737 [ - + ][ - - ]: 499 : Assert( lhs.nunk() == rhs.nunk(), "Number of unknowns unequal" );
[ - - ][ - - ]
738 [ - + ][ - - ]: 499 : Assert( lhs.nprop() == rhs.nprop(), "Number of properties unequal" );
[ - - ][ - - ]
739 : 499 : auto l = lhs.vec().cbegin();
740 : 499 : auto r = rhs.vec().cbegin();
741 : 499 : std::pair< std::size_t, tk::real > m( 0, std::abs(*l - *r) );
742 : 499 : ++l; ++r;
743 [ + + ]: 884765 : while (l != lhs.vec().cend()) {
744 : 884266 : const auto d = std::abs(*l - *r);
745 [ + + ][ + - ]: 884266 : if (d > m.second) m = { std::distance(lhs.vec().cbegin(),l), d };
746 : 884266 : ++l; ++r;
747 : : }
748 : 499 : return m;
749 : : }
750 : :
751 : : } // tk::
752 : :
753 : : #endif // Data_h
|
