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 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Generic data storage with different memory layouts
10 : : \details Generic data storage with different memory layouts. See also the
11 : : rationale discussed in the [design](layout.html) document.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef Data_h
15 : : #define Data_h
16 : :
17 : : #include <array>
18 : : #include <string>
19 : : #include <cstdint>
20 : : #include <vector>
21 : : #include <set>
22 : : #include <algorithm>
23 : :
24 : : #include "Types.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 : : using ncomp_t = uint64_t;
41 : :
42 : : public:
43 : : //! Default constructor (required for Charm++ migration)
44 : 40622 : explicit Data() : m_vec(), m_nunk(), m_nprop() {}
45 : :
46 : : //! Constructor
47 : : //! \param[in] nu Number of unknowns to allocate memory for
48 : : //! \param[in] np Total number of properties, i.e., scalar variables or
49 : : //! components, per unknown
50 : 14251 : explicit Data( ncomp_t nu, ncomp_t np ) :
51 [ + - ]: 14251 : m_vec( nu*np ),
52 : 14251 : m_nunk( nu ),
53 : 14251 : m_nprop( np ) {}
54 : :
55 : : //! Const data access dispatch
56 : : //! \details Public interface to const-ref data access to a single real
57 : : //! value. Use it as Data(p,c), where p is the unknown index, and c is
58 : : //! the component index specifying the scalar equation within a system of
59 : : //! equations. Requirement: component < nprop, unknown < nunk, enforced
60 : : //! with an assert in DEBUG mode, see also the constructor.
61 : : //! \param[in] unknown Unknown index
62 : : //! \param[in] component Component index, i.e., position of a scalar within
63 : : //! a system
64 : : //! \return Const reference to data of type tk::real
65 : : const tk::real&
66 : 6008401291 : operator()( ncomp_t unknown, ncomp_t component ) const
67 : 6008401291 : { return access( unknown, component, int2type< Layout >() ); }
68 : :
69 : : //! Non-const data access dispatch
70 : : //! \details Public interface to non-const-ref data access to a single real
71 : : //! value. Use it as Data(p,c,o), where p is the unknown index, and c is
72 : : //! the component index specifying the scalar equation within a system of
73 : : //! equations. Requirement: component < nprop, unknown < nunk, enforced
74 : : //! with an assert in DEBUG mode, see also the constructor.
75 : : //! \param[in] unknown Unknown index
76 : : //! \param[in] component Component index, i.e., position of a scalar within
77 : : //! a system
78 : : //! \return Non-const reference to data of type tk::real
79 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
80 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
81 : : tk::real&
82 : 3825186790 : operator()( ncomp_t unknown, ncomp_t component ) {
83 : : return const_cast< tk::real& >(
84 : : static_cast< const Data& >( *this ).
85 : 3825186790 : operator()( unknown, component ) );
86 : : }
87 : :
88 : : //! Access to number of unknowns
89 : : //! \return Number of unknowns
90 : 35210451 : ncomp_t nunk() const noexcept { return m_nunk; }
91 : :
92 : : //! Access to number of properties
93 : : //! \details This is the total number of scalar components per unknown
94 : : //! \return Number of propertes/unknown
95 : 192343370 : ncomp_t nprop() const noexcept { return m_nprop; }
96 : :
97 : : //! Extract vector of unknowns given component
98 : : //! \details Requirement: component < nprop, enforced with an assert in
99 : : //! DEBUG mode, see also the constructor.
100 : : //! \param[in] component Component index, i.e., position of a scalar within
101 : : //! a system
102 : : //! \return A vector of unknowns given by component (length: nunk(), i.e.,
103 : : //! the first constructor argument)
104 : : std::vector< tk::real >
105 : 26816 : extract( ncomp_t component ) const {
106 [ + - ]: 26816 : std::vector< tk::real > w( m_nunk );
107 [ + - ][ + + ]: 5819369 : for (ncomp_t i=0; i<m_nunk; ++i) w[i] = operator()( i, component );
108 : 26816 : return w;
109 : 0 : }
110 : :
111 : : //! Extract all components for unknown
112 : : //! \details Requirement: unknown < nunk, enforced with an assert in DEBUG
113 : : //! mode, see also the constructor.
114 : : //! \param[in] unknown Index of unknown
115 : : //! \return A vector of components for a single unknown (length: nprop,
116 : : //! i.e., the second constructor argument)
117 : : std::vector< tk::real >
118 : 11069924 : operator[]( ncomp_t unknown ) const {
119 [ + - ]: 11069924 : std::vector< tk::real > w( m_nprop );
120 [ + - ][ + + ]: 102054971 : for (ncomp_t i=0; i<m_nprop; ++i) w[i] = operator()( unknown, i );
121 : 11069924 : return w;
122 : 0 : }
123 : :
124 : : //! Extract (a copy of) four values of unknowns
125 : : //! \details Requirement: component < nprop, for all N[i] < nunk, enforced
126 : : //! with an assert in DEBUG mode, see also the constructor.
127 : : //! \param[in] component Component index, i.e., position of a scalar within
128 : : //! a system
129 : : //! \param[in] N Indices of the 4 unknowns
130 : : //! \return Array of the four values of component
131 : : std::array< tk::real, 4 >
132 : 6058952 : extract( ncomp_t component, const std::array< ncomp_t, 4 >& N ) const {
133 : 6058952 : auto p = cptr( component );
134 : 6058952 : return {{ var(p,N[0]), var(p,N[1]), var(p,N[2]), var(p,N[3]) }};
135 : : }
136 : :
137 : : //! Const-ref accessor to underlying raw data as a std::vector
138 : : //! \return Constant reference to underlying raw data
139 : 126 : const std::vector< tk::real >& vec() const { return m_vec; }
140 : :
141 : : //! Non-const-ref accessor to underlying raw data as a std::vector
142 : : //! \return Non-constant reference to underlying raw data
143 : 44 : std::vector< tk::real >& vec() { return m_vec; }
144 : :
145 : : //! Compound operator-=
146 : : //! \param[in] rhs Data object to subtract
147 : : //! \return Reference to ourselves after subtraction
148 : 4 : Data< Layout >& operator-= ( const Data< Layout >& rhs ) {
149 [ - + ][ - - ]: 4 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
150 [ - + ][ - - ]: 4 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
151 : 4 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
152 : : m_vec.cbegin(), m_vec.begin(),
153 : 24 : []( tk::real s, tk::real d ){ return d-s; } );
154 : 4 : return *this;
155 : : }
156 : : //! Operator -
157 : : //! \param[in] rhs Data object to subtract
158 : : //! \return Copy of Data object after rhs has been subtracted
159 : : //! \details Implemented in terms of compound operator-=
160 : 2 : Data< Layout > operator- ( const Data< Layout >& rhs )
161 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) -= rhs; }
[ + - ]
162 : :
163 : : //! Compound operator+=
164 : : //! \param[in] rhs Data object to add
165 : : //! \return Reference to ourselves after addition
166 : 4 : Data< Layout >& operator+= ( const Data< Layout >& rhs ) {
167 [ - + ][ - - ]: 4 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
168 [ - + ][ - - ]: 4 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
169 : 4 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
170 : : m_vec.cbegin(), m_vec.begin(),
171 : 24 : []( tk::real s, tk::real d ){ return d+s; } );
172 : 4 : return *this;
173 : : }
174 : : //! Operator +
175 : : //! \param[in] rhs Data object to add
176 : : //! \return Copy of Data object after rhs has been multiplied with
177 : : //! \details Implemented in terms of compound operator+=
178 : 2 : Data< Layout > operator+ ( const Data< Layout >& rhs )
179 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) += rhs; }
[ + - ]
180 : :
181 : : //! Compound operator*= multiplying by another Data object item by item
182 : : //! \param[in] rhs Data object to multiply with
183 : : //! \return Reference to ourselves after multiplication
184 : 4 : Data< Layout >& operator*= ( const Data< Layout >& rhs ) {
185 [ - + ][ - - ]: 4 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
186 [ - + ][ - - ]: 4 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
187 : 4 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
188 : : m_vec.cbegin(), m_vec.begin(),
189 : 24 : []( tk::real s, tk::real d ){ return d*s; } );
190 : 4 : return *this;
191 : : }
192 : : //! Operator * multiplying by another Data object item by item
193 : : //! \param[in] rhs Data object to multiply with
194 : : //! \return Copy of Data object after rhs has been multiplied with
195 : : //! \details Implemented in terms of compound operator*=
196 : 2 : Data< Layout > operator* ( const Data< Layout >& rhs )
197 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) *= rhs; }
[ + - ]
198 : :
199 : : //! Compound operator*= multiplying all items by a scalar
200 : : //! \param[in] rhs Scalar to multiply with
201 : : //! \return Reference to ourselves after multiplication
202 : 6 : Data< Layout >& operator*= ( tk::real rhs ) {
203 : : // cppcheck-suppress useStlAlgorithm
204 [ + + ]: 42 : for (auto& v : m_vec) v *= rhs;
205 : 6 : return *this;
206 : : }
207 : : //! Operator * multiplying all items by a scalar
208 : : //! \param[in] rhs Scalar to multiply with
209 : : //! \return Copy of Data object after rhs has been multiplied with
210 : : //! \details Implemented in terms of compound operator*=
211 : 2 : Data< Layout > operator* ( tk::real rhs )
212 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) *= rhs; }
213 : :
214 : : //! Compound operator/=
215 : : //! \param[in] rhs Data object to divide by
216 : : //! \return Reference to ourselves after division
217 : 4 : Data< Layout >& operator/= ( const Data< Layout >& rhs ) {
218 [ - + ][ - - ]: 4 : Assert( rhs.nunk() == m_nunk, "Incorrect number of unknowns" );
[ - - ][ - - ]
219 [ - + ][ - - ]: 4 : Assert( rhs.nprop() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
220 : 4 : std::transform( rhs.vec().cbegin(), rhs.vec().cend(),
221 : : m_vec.cbegin(), m_vec.begin(),
222 : 24 : []( tk::real s, tk::real d ){ return d/s; } );
223 : 4 : return *this;
224 : : }
225 : : //! Operator /
226 : : //! \param[in] rhs Data object to divide by
227 : : //! \return Copy of Data object after rhs has been divided by
228 : : //! \details Implemented in terms of compound operator/=
229 : 2 : Data< Layout > operator/ ( const Data< Layout >& rhs )
230 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) /= rhs; }
[ + - ]
231 : :
232 : : //! Compound operator/= dividing all items by a scalar
233 : : //! \param[in] rhs Scalar to divide with
234 : : //! \return Reference to ourselves after division
235 : 4 : Data< Layout >& operator/= ( tk::real rhs ) {
236 : : // cppcheck-suppress useStlAlgorithm
237 [ + + ]: 28 : for (auto& v : m_vec) v /= rhs;
238 : 4 : return *this;
239 : : }
240 : : //! Operator / dividing all items by a scalar
241 : : //! \param[in] rhs Scalar to divide with
242 : : //! \return Copy of Data object after rhs has been divided by
243 : : //! \details Implemented in terms of compound operator/=
244 : 2 : Data< Layout > operator/ ( tk::real rhs )
245 [ + - ][ + - ]: 4 : const { return Data< Layout >( *this ) /= rhs; }
246 : :
247 : : //! Add new unknown at the end of the container
248 : : //! \param[in] prop Vector of properties to initialize the new unknown with
249 : 1 : void push_back( const std::vector< tk::real >& prop )
250 : 1 : { return push_back( prop, int2type< Layout >() ); }
251 : :
252 : : //! Resize data store to contain 'count' elements
253 : : //! \param[in] count Resize store to contain count * nprop elements
254 : : //! \param[in] value Value to initialize new data with (default: 0.0)
255 : : //! \note This works for both shrinking and enlarging, as this simply
256 : : //! translates to std::vector::resize(). Note that count changes, nprop
257 : : //! does not, see the private overload resize().
258 : 2 : void resize( std::size_t count, tk::real value = 0.0 )
259 : 2 : { resize( count, value, int2type< Layout >() ); }
260 : :
261 : : //! Resize data given number of unknowns and number of properties
262 : : //! \param[in] nu Number of unknowns to allocate memory for
263 : : //! \param[in] np Total number of properties, i.e., scalar variables or
264 : : //! components, per unknown
265 : 6744 : void resize( std::size_t nu, std::size_t np ) {
266 : 6744 : m_vec.resize( nu * np );
267 : 6744 : m_nunk = nu;
268 : 6744 : m_nprop = np;
269 : 6744 : }
270 : :
271 : : //! Remove a number of unknowns
272 : : //! \param[in] unknown Set of indices of unknowns to remove
273 : 3 : void rm( const std::set< ncomp_t >& unknown ) {
274 : 25 : auto remove = [ &unknown ]( std::size_t i ) -> bool {
275 [ + - ][ + + ]: 11 : if (unknown.find(i) != end(unknown)) return true;
276 : 6 : return false;
277 : : };
278 : 3 : std::size_t last = 0;
279 [ + + ]: 7 : for(std::size_t i=0; i<m_nunk; ++i, ++last) {
280 [ + - ][ + + ]: 11 : while( remove(i) ) ++i;
281 [ + + ]: 6 : if (i >= m_nunk) break;
282 [ + + ]: 11 : for (ncomp_t p = 0; p<m_nprop; ++p)
283 : 7 : m_vec[ last*m_nprop+p ] = m_vec[ i*m_nprop+p ];
284 : : }
285 [ + - ]: 3 : m_vec.resize( last*m_nprop );
286 : 3 : m_nunk -= unknown.size();
287 : 3 : }
288 : :
289 : : //! Fill vector of unknowns with the same value
290 : : //! \details Requirement: component < nprop, enforced with an assert in
291 : : //! DEBUG mode, see also the constructor.
292 : : //! \param[in] component Component index, i.e., position of a scalar within
293 : : //! a system
294 : : //! \param[in] value Value to fill vector of unknowns with
295 : 4 : inline void fill( ncomp_t component, tk::real value ) {
296 : 4 : auto p = cptr( component );
297 [ + + ]: 16 : for (ncomp_t i=0; i<m_nunk; ++i) var(p,i) = value;
298 : 4 : }
299 : :
300 : : //! Fill full data storage with value
301 : : //! \param[in] value Value to fill data with
302 : 176009 : void fill( tk::real value )
303 : 176009 : { std::fill( begin(m_vec), end(m_vec), value ); }
304 : :
305 : : //! Check if vector of unknowns is empty
306 : : bool empty() const noexcept { return m_vec.empty(); }
307 : :
308 : : //! Layout name dispatch
309 : : //! \return The name of the data layout used
310 : 2 : static std::string layout() { return layout( int2type< Layout >() ); }
311 : :
312 : : /** @name Pack/Unpack: Serialize Data object for Charm++ */
313 : : ///@{
314 : : //! \brief Pack/Unpack serialize member function
315 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
316 : : // cppcheck-suppress constParameter
317 : 60541 : void pup( PUP::er &p ) {
318 : 60541 : p | m_vec;
319 : 60541 : p | m_nunk;
320 : 60541 : p | m_nprop;
321 : 60541 : }
322 : : //! \brief Pack/Unpack serialize operator|
323 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
324 : : //! \param[in,out] d DataLyaout object reference
325 : 60541 : friend void operator|( PUP::er& p, Data& d ) { d.pup(p); }
326 : : //@}
327 : :
328 : : private:
329 : : //! Transform a compile-time uint8_t into a type, used for dispatch
330 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
331 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
332 : : template< uint8_t m > struct int2type { enum { value = m }; };
333 : :
334 : : //! Const ptr to physical variable access dispatch
335 : : //! \details Public interface to the first half of a physical variable
336 : : //! access. cptr() and var() are two member functions intended to be used
337 : : //! together in case when component would be expensive to compute for data
338 : : //! access via the function call operator, i.e., cptr(), can be used to
339 : : //! pre-compute part of the address, which returns a pointer and var() can
340 : : //! be used to finish the data access using the pointer returned by
341 : : //! cptr(). In other words, cptr() returns part of the address known based
342 : : //! on component and intended to be used in a setup phase. Then var()
343 : : //! takes this partial address and finishes the address calculation given
344 : : //! the unknown id. Thus the following two data accesses are equivalent
345 : : //! (modulo constness):
346 : : //! * real& value = operator()( unk, comp ); and
347 : : //! * const real* p = cptr( comp ); and
348 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
349 : : //! Requirement: component < nprop, enforced with an assert in DEBUG mode,
350 : : //! see also the constructor.
351 : : //! \param[in] component Component index, i.e., position of a scalar within
352 : : //! a system
353 : : //! \return Pointer to data of type tk::real for use with var()
354 : : const tk::real*
355 : 6058956 : cptr( ncomp_t component ) const
356 : 6058956 : { return cptr( component, int2type< Layout >() ); }
357 : :
358 : : //! Const-ref data-access dispatch
359 : : //! \details Public interface to the second half of a physical variable
360 : : //! access. cptr() and var() are two member functions intended to be used
361 : : //! together in case when component would be expensive to compute for data
362 : : //! access via the function call operator, i.e., cptr(), can be used to
363 : : //! pre-compute part of the address, which returns a pointer and var() can
364 : : //! be used to finish the data access using the pointer returned by
365 : : //! cptr(). In other words, cptr() returns part of the address known based
366 : : //! on component and intended to be used in a setup phase. Then var()
367 : : //! takes this partial address and finishes the address calculation given
368 : : //! the unknown id. Thus the following two data accesses are equivalent
369 : : //! (modulo constness):
370 : : //! * real& value = operator()( unk, comp ); and
371 : : //! * const real* p = cptr( comp ); and
372 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
373 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
374 : : //! see also the constructor.
375 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
376 : : //! \param[in] unknown Unknown index
377 : : //! \return Const reference to data of type tk::real
378 : : const tk::real&
379 : 24235820 : var( const tk::real* pt, ncomp_t unknown ) const
380 : 24235820 : { return var( pt, unknown, int2type< Layout >() ); }
381 : :
382 : : //! Non-const-ref data-access dispatch
383 : : //! \details Public interface to the second half of a physical variable
384 : : //! access. cptr() and var() are two member functions intended to be used
385 : : //! together in case when component would be expensive to compute for data
386 : : //! access via the function call operator, i.e., cptr(), can be used to
387 : : //! pre-compute part of the address, which returns a pointer and var() can
388 : : //! be used to finish the data access using the pointer returned by
389 : : //! cptr(). In other words, cptr() returns part of the address known based
390 : : //! on component and intended to be used in a setup phase. Then var()
391 : : //! takes this partial address and finishes the address calculation given
392 : : //! the unknown id. Thus the following two data accesses are equivalent
393 : : //! (modulo constness):
394 : : //! * real& value = operator()( unk, comp ); and
395 : : //! * const real* p = cptr( comp ); and
396 : : //! const real& value = var( p, unk ); or real& value = var( p, unk );
397 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
398 : : //! see also the constructor.
399 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
400 : : //! \param[in] unknown Unknown index
401 : : //! \return Non-const reference to data of type tk::real
402 : : //! \see "Avoid Duplication in const and Non-const Member Function," and
403 : : //! "Use const whenever possible," Scott Meyers, Effective C++, 3d ed.
404 : : tk::real&
405 : 12 : var( const tk::real* pt, ncomp_t unknown ) {
406 : : return const_cast< tk::real& >(
407 : 12 : static_cast< const Data& >( *this ).var( pt, unknown ) );
408 : : }
409 : :
410 : : //! Overloads for the various const data accesses
411 : : //! \details Requirement: component < nprop, unknown < nunk, enforced with
412 : : //! an assert in DEBUG mode, see also the constructor.
413 : : //! \param[in] unknown Unknown index
414 : : //! \param[in] component Component index, i.e., position of a scalar within
415 : : //! a system
416 : : //! \return Const reference to data of type tk::real
417 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
418 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
419 : : const tk::real&
420 : 6008400916 : access( ncomp_t unknown, ncomp_t component, int2type< UnkEqComp > ) const {
421 [ + + ][ + - ]: 6008400916 : Assert( component < m_nprop,
[ + - ][ + - ]
422 : : "Out-of-bounds access: component < number of properties" );
423 [ + + ][ + - ]: 6008400910 : Assert( unknown < m_nunk,
[ + - ][ + - ]
424 : : "Out-of-bounds access: unknown < number of unknowns" );
425 : 6008400902 : return m_vec[ unknown*m_nprop + component ];
426 : : }
427 : : const tk::real&
428 : 375 : access( ncomp_t unknown, ncomp_t component, int2type< EqCompUnk > ) const {
429 [ + + ][ + - ]: 375 : Assert( component < m_nprop,
[ + - ][ + - ]
430 : : "Out-of-bounds access: component < number of properties" );
431 [ + + ][ + - ]: 369 : Assert( unknown < m_nunk,
[ + - ][ + - ]
432 : : "Out-of-bounds access: unknown < number of unknowns" );
433 : 362 : return m_vec[ component*m_nunk + unknown ];
434 : : }
435 : :
436 : : // Overloads for the various const ptr to physical variable accesses
437 : : //! \details Requirement: component < nprop, unknown < nunk, enforced with
438 : : //! an assert in DEBUG mode, see also the constructor.
439 : : //! \param[in] component Component index, i.e., position of a scalar within
440 : : //! a system
441 : : //! \return Pointer to data of type tk::real for use with var()
442 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
443 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
444 : : const tk::real*
445 : 6058952 : cptr( ncomp_t component, int2type< UnkEqComp > ) const {
446 [ - + ][ - - ]: 6058952 : Assert( component < m_nprop,
[ - - ][ - - ]
447 : : "Out-of-bounds access: component < number of properties" );
448 : 6058952 : return m_vec.data() + component;
449 : : }
450 : : const tk::real*
451 : 4 : cptr( ncomp_t component, int2type< EqCompUnk > ) const {
452 [ - + ][ - - ]: 4 : Assert( component < m_nprop,
[ - - ][ - - ]
453 : : "Out-of-bounds access: component < number of properties" );
454 : 4 : return m_vec.data() + component*m_nunk;
455 : : }
456 : :
457 : : // Overloads for the various const physical variable accesses
458 : : //! Requirement: unknown < nunk, enforced with an assert in DEBUG mode,
459 : : //! see also the constructor.
460 : : //! \param[in] pt Pointer to data of type tk::real as returned from cptr()
461 : : //! \param[in] unknown Unknown index
462 : : //! \return Const reference to data of type tk::real
463 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
464 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
465 : : inline const tk::real&
466 : 24235806 : var( const tk::real* const pt, ncomp_t unknown, int2type< UnkEqComp > )
467 : : const {
468 [ - + ][ - - ]: 24235806 : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
[ - - ][ - - ]
469 : : "unknowns" );
470 : 24235806 : return *(pt + unknown*m_nprop);
471 : : }
472 : : inline const tk::real&
473 : 14 : var( const tk::real* const pt, ncomp_t unknown, int2type< EqCompUnk > )
474 : : const {
475 [ - + ][ - - ]: 14 : Assert( unknown < m_nunk, "Out-of-bounds access: unknown < number of "
[ - - ][ - - ]
476 : : "unknowns" );
477 : 14 : return *(pt + unknown);
478 : : }
479 : :
480 : : //! Add new unknown
481 : : //! \param[in] prop Vector of properties to initialize the new unknown with
482 : : //! \note Only the UnkEqComp overload is provided as this operation would be
483 : : //! too inefficient with the EqCompUnk data layout.
484 : 1 : void push_back( const std::vector< tk::real >& prop, int2type< UnkEqComp > )
485 : : {
486 [ - + ][ - - ]: 1 : Assert( prop.size() == m_nprop, "Incorrect number of properties" );
[ - - ][ - - ]
487 : 1 : m_vec.resize( (m_nunk+1) * m_nprop );
488 : 1 : ncomp_t u = m_nunk;
489 : 1 : ++m_nunk;
490 [ + + ]: 3 : for (ncomp_t i=0; i<m_nprop; ++i) operator()( u, i ) = prop[i];
491 : 1 : }
492 : :
493 : : void push_back( const std::vector< tk::real >&, int2type< EqCompUnk > )
494 : : { Throw( "Not implented. It would be inefficient" ); }
495 : :
496 : : //! Resize data store to contain 'count' elements
497 : : //! \param[in] count Resize store to contain 'count' elements
498 : : //! \param[in] value Value to initialize new data with
499 : : //! \note Only the UnkEqComp overload is provided as this operation would be
500 : : //! too inefficient with the EqCompUnk data layout.
501 : : //! \note This works for both shrinking and enlarging, as this simply
502 : : //! translates to std::vector::resize().
503 : 2 : void resize( std::size_t count, tk::real value, int2type< UnkEqComp > ) {
504 : 2 : m_vec.resize( count * m_nprop, value );
505 : 2 : m_nunk = count;
506 : 2 : }
507 : :
508 : : void resize( std::size_t, tk::real, int2type< EqCompUnk > ) {
509 : : Throw( "Not implemented. It would be inefficient" );
510 : : }
511 : :
512 : : // Overloads for the name-queries of data lauouts
513 : : //! \return The name of the data layout used
514 : : //! \see A. Alexandrescu, Modern C++ Design: Generic Programming and Design
515 : : //! Patterns Applied, Addison-Wesley Professional, 2001.
516 : 1 : static std::string layout( int2type< UnkEqComp > )
517 [ + - ]: 1 : { return "unknown-major"; }
518 : 1 : static std::string layout( int2type< EqCompUnk > )
519 [ + - ]: 1 : { return "equation-major"; }
520 : :
521 : : std::vector< tk::real > m_vec; //!< Data pointer
522 : : ncomp_t m_nunk; //!< Number of unknowns
523 : : ncomp_t m_nprop; //!< Number of properties/unknown
524 : : };
525 : :
526 : : //! Operator * multiplying all items by a scalar from the left
527 : : //! \param[in] lhs Scalar to multiply with
528 : : //! \param[in] rhs Date object to multiply
529 : : //! \return New Data object with all items multipled with lhs
530 : : template< uint8_t Layout >
531 : 2 : Data< Layout > operator* ( tk::real lhs, const Data< Layout >& rhs ) {
532 [ + - ][ + - ]: 4 : return Data< Layout >( rhs ) *= lhs;
533 : : }
534 : :
535 : : //! Operator min between two Data objects
536 : : //! \param[in] a 1st Data object
537 : : //! \param[in] b 2nd Data object
538 : : //! \return New Data object containing the minimum of all values for each
539 : : //! value in _a_ and _b_
540 : : //! \note The Data objects _a_ and _b_ must have the same number of
541 : : //! unknowns and properties.
542 : : //! \note As opposed to std::min, this function creates and returns a new object
543 : : //! instead of returning a reference to one of the operands.
544 : : template< uint8_t Layout >
545 : 2 : Data< Layout > min( const Data< Layout >& a, const Data< Layout >& b ) {
546 [ - + ][ - - ]: 2 : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
[ - - ][ - - ]
547 [ - + ][ - - ]: 2 : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
[ - - ][ - - ]
548 : 2 : Data< Layout > r( a.nunk(), a.nprop() );
549 [ + - ]: 4 : std::transform( a.vec().cbegin(), a.vec().cend(),
550 : 4 : b.vec().cbegin(), r.vec().begin(),
551 : 12 : []( tk::real s, tk::real d ){ return std::min(s,d); } );
552 : :
553 : 2 : return r;
554 : 0 : }
555 : :
556 : : //! Operator max between two Data objects
557 : : //! \param[in] a 1st Data object
558 : : //! \param[in] b 2nd Data object
559 : : //! \return New Data object containing the maximum of all values for each
560 : : //! value in _a_ and _b_
561 : : //! \note The Data objects _a_ and _b_ must have the same number of
562 : : //! unknowns and properties.
563 : : //! \note As opposed to std::max, this function creates and returns a new object
564 : : //! instead of returning a reference to one of the operands.
565 : : template< uint8_t Layout >
566 : 2 : Data< Layout > max( const Data< Layout >& a, const Data< Layout >& b ) {
567 [ - + ][ - - ]: 2 : Assert( a.nunk() == b.nunk(), "Number of unknowns unequal" );
[ - - ][ - - ]
568 [ - + ][ - - ]: 2 : Assert( a.nprop() == b.nprop(), "Number of properties unequal" );
[ - - ][ - - ]
569 : 2 : Data< Layout > r( a.nunk(), a.nprop() );
570 [ + - ]: 4 : std::transform( a.vec().cbegin(), a.vec().cend(),
571 : 4 : b.vec().cbegin(), r.vec().begin(),
572 : 12 : []( tk::real s, tk::real d ){ return std::max(s,d); } );
573 : 2 : return r;
574 : 0 : }
575 : :
576 : : //! Operator == between two Data objects
577 : : //! \param[in] lhs Data object to compare
578 : : //! \param[in] rhs Data object to compare
579 : : //! \return True if all entries are equal up to epsilon
580 : : template< uint8_t Layout >
581 : 8 : bool operator== ( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
582 [ - + ][ - - ]: 8 : Assert( rhs.nunk() == lhs.nunk(), "Incorrect number of unknowns" );
[ - - ][ - - ]
583 [ - + ][ - - ]: 8 : Assert( rhs.nprop() == lhs.nprop(), "Incorrect number of properties" );
[ - - ][ - - ]
584 : 8 : auto l = lhs.vec().cbegin();
585 : 8 : auto r = rhs.vec().cbegin();
586 [ + + ]: 32 : while (l != lhs.vec().cend()) {
587 [ + + ]: 28 : if (std::abs(*l - *r) > std::numeric_limits< tk::real >::epsilon())
588 : 4 : return false;
589 : 24 : ++l; ++r;
590 : : }
591 : 4 : return true;
592 : : }
593 : :
594 : : //! Operator != between two Data objects
595 : : //! \param[in] lhs Data object to compare
596 : : //! \param[in] rhs Data object to compare
597 : : //! \return True if all entries are unequal up to epsilon
598 : : template< uint8_t Layout >
599 : 4 : bool operator!= ( const Data< Layout >& lhs, const Data< Layout >& rhs )
600 : 4 : { return !(lhs == rhs); }
601 : :
602 : : //! Compute the maximum difference between the elements of two Data objects
603 : : //! \param[in] lhs 1st Data object
604 : : //! \param[in] rhs 2nd Data object
605 : : //! \return The index, i.e., the raw position, of and the largest absolute value
606 : : //! of the difference between all corresponding elements of _lhs_ and _rhs_.
607 : : //! \details The position returned is the position in the underlying raw data
608 : : //! structure, independent of components. If lhs == rhs with precision
609 : : //! std::numeric_limits< tk::real >::epsilon(), a pair of (0,0.0) is returned.
610 : : //! \note The Data objects _lhs_ and _rhs_ must have the same number of
611 : : //! unknowns and properties.
612 : : template< uint8_t Layout >
613 : : std::pair< std::size_t, tk::real >
614 : 4 : maxdiff( const Data< Layout >& lhs, const Data< Layout >& rhs ) {
615 [ - + ][ - - ]: 4 : Assert( lhs.nunk() == rhs.nunk(), "Number of unknowns unequal" );
[ - - ][ - - ]
616 [ - + ][ - - ]: 4 : Assert( lhs.nprop() == rhs.nprop(), "Number of properties unequal" );
[ - - ][ - - ]
617 : 4 : auto l = lhs.vec().cbegin();
618 : 4 : auto r = rhs.vec().cbegin();
619 : 4 : std::pair< std::size_t, tk::real > m( 0, std::abs(*l - *r) );
620 : 4 : ++l; ++r;
621 [ + + ]: 24 : while (l != lhs.vec().cend()) {
622 : 20 : const auto d = std::abs(*l - *r);
623 [ + + ][ + - ]: 20 : if (d > m.second) m = { std::distance(lhs.vec().cbegin(),l), d };
624 : 20 : ++l; ++r;
625 : : }
626 : 4 : return m;
627 : : }
628 : :
629 : : } // tk::
630 : :
631 : : #endif // Data_h
|
