BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_Database.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California,
4 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
12#ifndef BELFEM_CL_DATABASE_HPP
13#define BELFEM_CL_DATABASE_HPP
14
15#include "typedefs.hpp"
16#include "cl_Vector.hpp"
17#include "cl_TensorMeshFactory.hpp"
18#include "cl_Mesh.hpp"
19#include "hdf5_tools.hpp"
20
21namespace belfem
22{
29 class Database
30 {
31 const proc_t mCommRank ;
32 const proc_t mCommSize ;
33 const string mLabel ;
34
35 Vector< real > mValues ;
36 Matrix< index_t > mTopology ;
37 TensorMeshConfig * mConfig = nullptr ;
38
39 void ( Database::*mFunction2D ) ( const real xi, const real eta, real * N ) const = nullptr ;
40 void ( Database::*mFunction3D ) ( const real xi, const real eta, const real zeta, real * N ) const = nullptr ;
41
42 void ( Database::*mdFunction2Ddxi ) ( const real xi, const real eta, real * N ) const = nullptr ;
43 void ( Database::*mdFunction3Ddxi ) ( const real xi, const real eta, const real zeta, real * N ) const = nullptr ;
44
45 void ( Database::*mdFunction2Ddeta ) ( const real xi, const real eta, real * N ) const = nullptr ;
46 void ( Database::*mdFunction3Ddeta ) ( const real xi, const real eta, const real zeta, real * N ) const = nullptr ;
47
48 void ( Database::*mdFunction3Ddzeta ) ( const real xi, const real eta, const real zeta, real * N ) const = nullptr ;
49
50
51 real * mN = nullptr ;
52 uint mNumNodesPerElement = 0 ;
53
54 public:
55
56 Database( const string & aFilePath, const string & aMaterial );
57
58 Database( const hid_t aHDF5, const string & aLabel );
59
60 Database( Mesh * aMesh, const string & aField, const bool aProject = true, const string aMaterial="" );
61
62 ~Database();
63
64 real
65 evaluate( const real aX, const real aY ) const ;
66
67 real
68 evaluate( const real aX, const real aY, const real aZ ) const ;
69
70 real
71 evaluate_derivx( const real aX, const real aY ) const ;
72
73 real
74 evaluate_derivy( const real aX, const real aY ) const ;
75
76 real
77 evaluate_derivx( const real aX, const real aY, const real aZ ) const ;
78
79 real
80 evaluate_derivy( const real aX, const real aY, const real aZ ) const ;
81
82 real
83 evaluate_derivz( const real aX, const real aY, const real aZ ) const ;
84
85 real
86 min( const uint aDimension ) const ;
87
88 real
89 max( const uint aDimension ) const ;
90
93 const Vector< real > &
94 values() const
95 {
96 return mValues ;
97 }
98
99 void
100 save( const hid_t aHDF5 );
101
102 private:
103
104 void
105 load( const hid_t aHDF5 );
106
107 void
108 set_element_type();
109
110 index_t
111 nidx( const index_t i, const index_t j ) const;
112
113 index_t
114 nidx( const index_t i, const index_t j, const index_t k ) const;
115
116 void
117 eval_quad4( const real xi, const real eta, real * N ) const ;
118
119 void
120 eval_quad9( const real xi, const real eta, real * N ) const ;
121
122 void
123 eval_quad16( const real xi, const real eta, real * N ) const ;
124
125 void
126 eval_hex8( const real xi, const real eta, const real zeta, real * N ) const ;
127
128 void
129 eval_hex27( const real xi, const real eta, const real zeta, real * N ) const ;
130
131 void
132 eval_hex64( const real xi, const real eta, const real zeta, real * N ) const ;
133
134 void
135 deval_quad4dxi( const real xi, const real eta, real * N ) const ;
136
137 void
138 deval_quad4deta( const real xi, const real eta, real * N ) const ;
139
140 void
141 deval_quad9dxi( const real xi, const real eta, real * N ) const ;
142
143 void
144 deval_quad9deta( const real xi, const real eta, real * N ) const ;
145
146
147 void
148 deval_quad16dxi( const real xi, const real eta, real * N ) const ;
149
150 void
151 deval_quad16deta( const real xi, const real eta, real * N ) const ;
152
153 void
154 deval_hex8dxi( const real xi, const real eta, const real zeta, real * N ) const ;
155
156 void
157 deval_hex8deta( const real xi, const real eta, const real zeta, real * N ) const ;
158
159 void
160 deval_hex8dzeta( const real xi, const real eta, const real zeta, real * N ) const ;
161
162 void
163 deval_hex27dxi( const real xi, const real eta, const real zeta, real * N ) const ;
164
165 void
166 deval_hex64dxi( const real xi, const real eta, const real zeta, real * N ) const ;
167
168
169 void
170 deval_hex27deta( const real xi, const real eta, const real zeta, real * N ) const ;
171
172 void
173 deval_hex64deta( const real xi, const real eta, const real zeta, real * N ) const ;
174
175 void
176 deval_hex27dzeta( const real xi, const real eta, const real zeta, real * N ) const ;
177
178 void
179 deval_hex64dzeta( const real xi, const real eta, const real zeta, real * N ) const ;
180 };
181
182 inline index_t
183 Database::nidx( const index_t i, const index_t j ) const
184 {
185 return i * mConfig->num_nodes( 1 ) + j ;
186 }
187
188 index_t
189 inline Database::nidx( const index_t i, const index_t j, const index_t k ) const
190 {
191 return mConfig->num_nodes( 2 ) * ( mConfig->num_nodes( 1 ) * i + j ) + k ;
192 }
193
194 inline
195 real Database::evaluate( const real aX, const real aY ) const
196 {
197 // find element index
198 index_t i = mConfig->element_ijk( 0, aX );
199 index_t j = mConfig->element_ijk( 1, aY );
200 index_t e = mConfig->element_index( i, j );
201
202 // get center point of element
203 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
204 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
205
206 // compute parameter coordinates
207 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
208 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
209
210 // evaluate shape function
211 ( this->*mFunction2D )( xi, eta, mN );
212
213 // interpolate value
214 real aValue = 0.0 ;
215 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
216 {
217 aValue += mN[n] * mValues( mTopology( n, e ) );
218 }
219 return aValue;
220 }
221
222 inline
223 real Database::evaluate_derivx( const real aX, const real aY ) const
224 {
225 // find element index
226 index_t i = mConfig->element_ijk( 0, aX );
227 index_t j = mConfig->element_ijk( 1, aY );
228 index_t e = mConfig->element_index( i, j );
229
230 // get center point of element
231 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
232 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
233
234 // compute parameter coordinates
235 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
236 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
237
238 // evaluate shape function
239 ( this->*mdFunction2Ddxi )( xi, eta, mN );
240
241 // interpolate value
242 real aValue = 0.0 ;
243 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
244 {
245 aValue += mN[n] * mValues( mTopology( n, e ) );
246 }
247 aValue *= 2. * mConfig->inv_element_step( 0 ) ;
248 return aValue;
249 }
250
251 inline
252 real Database::evaluate_derivy( const real aX, const real aY ) const
253 {
254 // find element index
255 index_t i = mConfig->element_ijk( 0, aX );
256 index_t j = mConfig->element_ijk( 1, aY );
257 index_t e = mConfig->element_index( i, j );
258
259 // get center point of element
260 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
261 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
262
263 // compute parameter coordinates
264 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
265 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
266
267 // evaluate shape function
268 ( this->*mdFunction2Ddeta )( xi, eta, mN );
269
270 // interpolate value
271 real aValue = 0.0 ;
272 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
273 {
274 aValue += mN[n] * mValues( mTopology( n, e ) );
275 }
276 aValue *= 2. * mConfig->inv_element_step( 1 ) ;
277 return aValue;
278 }
279
280 inline
281 real Database::evaluate( const real aX, const real aY, const real aZ ) const
282 {
283 // find element index
284 index_t i = mConfig->element_ijk( 0, aX );
285 index_t j = mConfig->element_ijk( 1, aY );
286 index_t k = mConfig->element_ijk( 2, aZ );
287 index_t e = mConfig->element_index( i, j, k );
288
289 // get center point of element
290 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
291 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
292 real zc = mConfig->min( 2 ) + ( k + 0.5 ) * mConfig->element_step( 2 );
293
294 // compute parameter coordinates
295 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
296 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
297 real zeta = 2. * ( aZ - zc ) * mConfig->inv_element_step( 2 ) ;
298
299 // evaluate the shape function
300 ( this->*mFunction3D )( xi, eta, zeta, mN );
301
302 // interpolate value
303 real aValue = 0.0 ;
304 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
305 {
306 aValue += mN[n] * mValues( mTopology( n, e ) );
307 }
308 return aValue;
309 }
310
311 inline
312 real Database::evaluate_derivx( const real aX, const real aY, const real aZ ) const
313 {
314 // find element index
315 index_t i = mConfig->element_ijk( 0, aX );
316 index_t j = mConfig->element_ijk( 1, aY );
317 index_t k = mConfig->element_ijk( 2, aZ );
318 index_t e = mConfig->element_index( i, j, k );
319
320 // get center point of element
321 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
322 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
323 real zc = mConfig->min( 2 ) + ( k + 0.5 ) * mConfig->element_step( 2 );
324
325 // compute parameter coordinates
326 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
327 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
328 real zeta = 2. * ( aZ - zc ) * mConfig->inv_element_step( 2 ) ;
329
330 // evaluate the shape function
331 ( this->*mdFunction3Ddxi )( xi, eta, zeta, mN );
332
333 // interpolate value
334 real aValue = 0.0 ;
335 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
336 {
337 aValue += mN[n] * mValues( mTopology( n, e ) );
338 }
339 aValue *= 2. * mConfig->inv_element_step( 0 ) ;
340 return aValue;
341 }
342
343 inline
344 real Database::evaluate_derivy( const real aX, const real aY, const real aZ ) const
345 {
346 // find element index
347 index_t i = mConfig->element_ijk( 0, aX );
348 index_t j = mConfig->element_ijk( 1, aY );
349 index_t k = mConfig->element_ijk( 2, aZ );
350 index_t e = mConfig->element_index( i, j, k );
351
352 // get center point of element
353 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
354 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
355 real zc = mConfig->min( 2 ) + ( k + 0.5 ) * mConfig->element_step( 2 );
356
357 // compute parameter coordinates
358 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
359 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
360 real zeta = 2. * ( aZ - zc ) * mConfig->inv_element_step( 2 ) ;
361
362 // evaluate the shape function
363 ( this->*mdFunction3Ddeta )( xi, eta, zeta, mN );
364
365 // interpolate value
366 real aValue = 0.0 ;
367 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
368 {
369 aValue += mN[n] * mValues( mTopology( n, e ) );
370 }
371 aValue *= 2. * mConfig->inv_element_step( 1 ) ;
372 return aValue;
373 }
374
375 inline
376 real Database::evaluate_derivz( const real aX, const real aY, const real aZ ) const
377 {
378 // find element index
379 index_t i = mConfig->element_ijk( 0, aX );
380 index_t j = mConfig->element_ijk( 1, aY );
381 index_t k = mConfig->element_ijk( 2, aZ );
382 index_t e = mConfig->element_index( i, j, k );
383
384 // get center point of element
385 real xc = mConfig->min( 0 ) + ( i + 0.5 ) * mConfig->element_step( 0 );
386 real yc = mConfig->min( 1 ) + ( j + 0.5 ) * mConfig->element_step( 1 );
387 real zc = mConfig->min( 2 ) + ( k + 0.5 ) * mConfig->element_step( 2 );
388
389 // compute parameter coordinates
390 real xi = 2. * ( aX - xc ) * mConfig->inv_element_step( 0 ) ;
391 real eta = 2. * ( aY - yc ) * mConfig->inv_element_step( 1 ) ;
392 real zeta = 2. * ( aZ - zc ) * mConfig->inv_element_step( 2 ) ;
393
394 // evaluate the shape function
395 ( this->*mdFunction3Ddzeta )( xi, eta, zeta, mN );
396
397 // interpolate value
398 real aValue = 0.0 ;
399 for ( index_t n = 0; n < mNumNodesPerElement ; ++n )
400 {
401 aValue += mN[n] * mValues( mTopology( n, e ) );
402 }
403 aValue *= 2. * mConfig->inv_element_step( 2 ) ;
404 return aValue;
405 }
406
407
408 inline void
409 Database::eval_quad4( const real xi, const real eta, real * N ) const
410 {
411 N[0] = ( ( 1.0 - xi ) * ( 1.0 - eta ) ) * 0.25;
412 N[1] = ( ( 1.0 + xi ) * ( 1.0 - eta ) ) * 0.25;
413 N[2] = ( ( 1.0 + xi ) * ( 1.0 + eta ) ) * 0.25;
414 N[3] = ( ( 1.0 - xi ) * ( 1.0 + eta ) ) * 0.25;
415 }
416
417 inline void
418 Database::eval_quad9( const real xi, const real eta, real * N ) const
419 {
420 const real c = xi * eta * 0.25;
421 const real xi2 = xi*xi;
422 const real eta2 = eta*eta;
423
424 N[0] = ( c * ( eta - 1.0 ) * (xi - 1.0) );
425 N[1] = ( c * ( eta - 1.0 ) * (xi + 1.0) );
426 N[2] = ( c * ( eta + 1.0 ) * (xi + 1.0) );
427 N[3] = ( c * ( eta + 1.0 ) * (xi - 1.0) );
428 N[4] = ( eta * ( 1.0 - xi2 ) * ( eta - 1.0 ) ) * 0.5;
429 N[5] = ( xi * ( 1.0 - eta2)*( xi + 1.0 ) )*0.5;
430 N[6] = ( eta * (1.0 - xi2)*( eta + 1.0 ) )*0.5;
431 N[7] = ( xi*( 1.0 - eta2 )*( xi - 1.0 ) )*0.5;
432 N[8] = ( eta2 - 1.0 )*( xi2 - 1.0 );
433 }
434
435 inline void
436 Database::eval_quad16( const real xi, const real eta, real * N ) const
437 {
438 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 )*0.0625;
439 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) )*0.0625;
440 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) )*0.0625;
441 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 )*0.0625;
442
443 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 )*0.0625;
444 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) )*0.0625;
445 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) )*0.0625;
446 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 )*0.0625;
447
448 N[ 0 ] = a0*b0;
449 N[ 1 ] = a3*b0;
450 N[ 2 ] = a3*b3;
451 N[ 3 ] = a0*b3;
452 N[ 4 ] = a1*b0;
453 N[ 5 ] = a2*b0;
454 N[ 6 ] = a3*b1;
455 N[ 7 ] = a3*b2;
456 N[ 8 ] = a2*b3;
457 N[ 9 ] = a1*b3;
458 N[ 10 ] = a0*b2;
459 N[ 11 ] = a0*b1;
460 N[ 12 ] = a1*b1;
461 N[ 13 ] = a2*b1;
462 N[ 14 ] = a2*b2;
463 N[ 15 ] = a1*b2;
464 }
465
466 inline void
467 Database::eval_hex8( const real xi, const real eta, const real zeta, real * N ) const
468 {
469 N[0] = - ( eta - 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 ) * 0.125;
470 N[1] = ( eta - 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 ) * 0.125;
471 N[2] = - ( eta + 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 ) * 0.125;
472 N[3] = ( eta + 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 ) * 0.125;
473 N[4] = ( eta - 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 ) * 0.125;
474 N[5] = - ( eta - 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 ) * 0.125;
475 N[6] = ( eta + 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 ) * 0.125;
476 N[7] = - ( eta + 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 ) * 0.125;
477 }
478
479 inline void
480 Database::eval_hex27( const real xi, const real eta, const real zeta, real * N ) const
481 {
482 const real xi2 = xi*xi;
483 const real eta2 = eta*eta;
484 const real zeta2 = zeta*zeta;
485
486 const real a = -0.25 * eta * zeta;
487 const real b = -0.25 * xi * zeta;
488 const real c = -0.25 * xi * eta;
489 const real d = 0.125 * xi * eta * zeta;
490
491 N[ 0 ] = d * ( eta - 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 );
492 N[ 1 ] = d * ( eta - 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 );
493 N[ 2 ] = d * ( eta + 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 );
494 N[ 3 ] = d * ( eta + 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 );
495 N[ 4 ] = d * ( eta - 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 );
496 N[ 5 ] = d * ( eta - 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 );
497 N[ 6 ] = d * ( eta + 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 );
498 N[ 7 ] = d * ( eta + 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 );
499 N[ 8 ] = a * ( xi2 - 1.0 ) * ( eta - 1.0 ) * ( zeta - 1.0 );
500 N[ 9 ] = b * ( eta2 - 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 );
501 N[ 10 ] = a * ( xi2 - 1.0 ) * ( eta + 1.0 ) * ( zeta - 1.0 );
502 N[ 11 ] = b * ( eta2 - 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 );
503 N[ 12 ] = c * ( zeta2 - 1.0 ) * ( eta - 1.0 ) * ( xi - 1.0 );
504 N[ 13 ] = c * ( zeta2 - 1.0 ) * ( eta - 1.0 ) * ( xi + 1.0 );
505 N[ 14 ] = c * ( zeta2 - 1.0 ) * ( eta + 1.0 ) * ( xi + 1.0 );
506 N[ 15 ] = c * ( zeta2 - 1.0 ) * ( eta + 1.0 ) * ( xi - 1.0 );
507 N[ 16 ] = a * ( xi2 - 1.0 ) * ( eta - 1.0 ) * ( zeta + 1.0 );
508 N[ 17 ] = b * ( eta2 - 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 );
509 N[ 18 ] = a * ( xi2 - 1.0 ) * ( eta + 1.0 ) * ( zeta + 1.0 );
510 N[ 19 ] = b * ( eta2 - 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 );
511 N[ 20 ] = -( eta2 - 1.0 ) * ( xi2 - 1.0 ) * ( zeta2 - 1.0 );
512 N[ 21 ] = ( zeta * ( eta2 - 1.0 ) * ( xi2 - 1.0 ) * ( zeta - 1.0 ) ) * 0.5;
513 N[ 22 ] = ( zeta * ( eta2 - 1.0 ) * ( xi2 - 1.0 ) * ( zeta + 1.0 ) ) * 0.5;
514 N[ 23 ] = ( xi * ( eta2 - 1.0 ) * ( zeta2 - 1.0 ) * ( xi - 1.0 ) ) * 0.5;
515 N[ 24 ] = ( xi * ( eta2 - 1.0 ) * ( zeta2 - 1.0 ) * ( xi + 1.0 ) ) * 0.5;
516 N[ 25 ] = ( eta * ( xi2 - 1.0 ) * ( zeta2 - 1.0 ) * ( eta - 1.0 ) ) * 0.5;
517 N[ 26 ] = ( eta * ( xi2 - 1.0 ) * ( zeta2 - 1.0 ) * ( eta + 1.0 ) ) * 0.5;
518 }
519
520 inline void
521 Database::eval_hex64( const real xi, const real eta, const real zeta, real * N ) const
522 {
523 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 )*0.0625;
524 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) )*0.0625;
525 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) )*0.0625;
526 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 )*0.0625;
527
528 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 )*0.0625;
529 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) )*0.0625;
530 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) )*0.0625;
531 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 )*0.0625;
532
533 const real c0 = ( zeta*( 1.0 + 9.0 * zeta * ( 1.0 - zeta ) ) - 1.0 )*0.0625;
534 const real c1 = ( 9.0 - zeta * ( 27.0 + zeta*( 9.0 - 27.0*zeta ) ) )*0.0625;
535 const real c2 = ( 9.0 + zeta * ( 27.0 - zeta*( 9.0 + 27.0*zeta ) ) )*0.0625;
536 const real c3 = ( -zeta*( 1.0 - 9.0 * zeta * ( 1.0 + zeta ) ) - 1.0 )*0.0625;
537
538 N[ 0 ] = a0 * b0 * c0;
539 N[ 1 ] = a3 * b0 * c0;
540 N[ 2 ] = a3 * b3 * c0;
541 N[ 3 ] = a0 * b3 * c0;
542 N[ 4 ] = a0 * b0 * c3;
543 N[ 5 ] = a3 * b0 * c3;
544 N[ 6 ] = a3 * b3 * c3;
545 N[ 7 ] = a0 * b3 * c3;
546 N[ 8 ] = a1 * b0 * c0;
547 N[ 9 ] = a2 * b0 * c0;
548 N[ 10 ] = a0 * b1 * c0;
549 N[ 11 ] = a0 * b2 * c0;
550 N[ 12 ] = a0 * b0 * c1;
551 N[ 13 ] = a0 * b0 * c2;
552 N[ 14 ] = a3 * b1 * c0;
553 N[ 15 ] = a3 * b2 * c0;
554 N[ 16 ] = a3 * b0 * c1;
555 N[ 17 ] = a3 * b0 * c2;
556 N[ 18 ] = a2 * b3 * c0;
557 N[ 19 ] = a1 * b3 * c0;
558 N[ 20 ] = a3 * b3 * c1;
559 N[ 21 ] = a3 * b3 * c2;
560 N[ 22 ] = a0 * b3 * c1;
561 N[ 23 ] = a0 * b3 * c2;
562 N[ 24 ] = a1 * b0 * c3;
563 N[ 25 ] = a2 * b0 * c3;
564 N[ 26 ] = a0 * b1 * c3;
565 N[ 27 ] = a0 * b2 * c3;
566 N[ 28 ] = a3 * b1 * c3;
567 N[ 29 ] = a3 * b2 * c3;
568 N[ 30 ] = a2 * b3 * c3;
569 N[ 31 ] = a1 * b3 * c3;
570 N[ 32 ] = a1 * b1 * c0;
571 N[ 33 ] = a1 * b2 * c0;
572 N[ 34 ] = a2 * b2 * c0;
573 N[ 35 ] = a2 * b1 * c0;
574 N[ 36 ] = a1 * b0 * c1;
575 N[ 37 ] = a2 * b0 * c1;
576 N[ 38 ] = a2 * b0 * c2;
577 N[ 39 ] = a1 * b0 * c2;
578 N[ 40 ] = a0 * b1 * c1;
579 N[ 41 ] = a0 * b1 * c2;
580 N[ 42 ] = a0 * b2 * c2;
581 N[ 43 ] = a0 * b2 * c1;
582 N[ 44 ] = a3 * b1 * c1;
583 N[ 45 ] = a3 * b2 * c1;
584 N[ 46 ] = a3 * b2 * c2;
585 N[ 47 ] = a3 * b1 * c2;
586 N[ 48 ] = a2 * b3 * c1;
587 N[ 49 ] = a1 * b3 * c1;
588 N[ 50 ] = a1 * b3 * c2;
589 N[ 51 ] = a2 * b3 * c2;
590 N[ 52 ] = a1 * b1 * c3;
591 N[ 53 ] = a2 * b1 * c3;
592 N[ 54 ] = a2 * b2 * c3;
593 N[ 55 ] = a1 * b2 * c3;
594 N[ 56 ] = a1 * b1 * c1;
595 N[ 57 ] = a2 * b1 * c1;
596 N[ 58 ] = a2 * b2 * c1;
597 N[ 59 ] = a1 * b2 * c1;
598 N[ 60 ] = a1 * b1 * c2;
599 N[ 61 ] = a2 * b1 * c2;
600 N[ 62 ] = a2 * b2 * c2;
601 N[ 63 ] = a1 * b2 * c2;
602 }
603
604 inline void
605 Database::deval_quad4dxi( const real xi, const real eta, real * N ) const
606 {
607 N[0] = - 0.25 * ( 1.0 - eta ) ;
608 N[1] = 0.25 * ( 1.0 - eta ) ;
609 N[2] = 0.25 * ( 1.0 + eta ) ;
610 N[3] = - 0.25 * ( 1.0 + eta ) ;
611 }
612
613 inline void
614 Database::deval_quad4deta( const real xi, const real eta, real * N ) const
615 {
616 N[ 0 ] = 0.25 * ( xi - 1.0 );
617 N[ 1 ] = -0.25 * ( xi + 1.0 );
618 N[ 2 ] = 0.25 * ( xi + 1.0 );
619 N[ 3 ] = -0.25 * ( xi - 1.0 );
620 }
621
622 inline void
623 Database::deval_quad9dxi( const real xi, const real eta, real * N ) const
624 {
625 const real c = xi * eta ;
626 //const real xi2 = xi*xi;
627 const real eta2 = eta*eta;
628
629 N[0] = ( eta * ( 2.0 * xi - 1.0 ) * ( eta - 1.0 ) ) * 0.25;
630 N[1] = ( eta * ( 2.0 * xi + 1.0 ) * ( eta - 1.0 ) ) * 0.25;
631 N[2] = ( eta * ( 2.0 * xi + 1.0 ) * ( eta + 1.0 ) ) * 0.25;
632 N[3] = ( eta * ( 2.0 * xi - 1.0 ) * ( eta + 1.0 ) ) * 0.25;
633 N[4] = - c * ( eta - 1.0 );
634 N[5] = -( ( eta2 - 1.0 ) * ( 2.0 * xi + 1.0 ) ) * 0.5;
635 N[6] = - c * ( eta + 1.0 );
636 N[7] = -( ( eta2 - 1.0 ) * ( 2.0 * xi - 1.0 ) ) * 0.5;
637 N[8] = 2.0 * xi * ( eta2 - 1.0 );
638 }
639
640 inline void
641 Database::deval_quad9deta( const real xi, const real eta, real * N ) const
642 {
643 const real c = xi * eta ;
644 const real xi2 = xi*xi;
645 // const real eta2 = eta*eta;
646
647 N[ 0 ] = ( xi * ( 2.0 * eta - 1.0 ) * ( xi - 1.0 ) ) * 0.25;
648 N[ 1 ] = ( xi * ( 2.0 * eta - 1.0 ) * ( xi + 1.0 ) ) * 0.25;
649 N[ 2 ] = ( xi * ( 2.0 * eta + 1.0 ) * ( xi + 1.0 ) ) * 0.25;
650 N[ 3 ] = ( xi * ( 2.0 * eta + 1.0 ) * ( xi - 1.0 ) ) * 0.25;
651 N[ 4 ] = -( ( 2.0 * eta - 1.0 ) * ( xi2 - 1.0 ) ) * 0.5;
652 N[ 5 ] = - c * ( xi + 1.0 );
653 N[ 6 ] = -( ( 2.0 * eta + 1.0 ) * ( xi2 - 1.0 ) ) * 0.5;
654 N[ 7 ] = - c * ( xi - 1.0 );
655 N[ 8 ] = 2.0 * eta * ( xi2 - 1.0 );
656 }
657
658 inline void
659 Database::deval_quad16dxi( const real xi, const real eta, real * N ) const
660 {
661 const real da0 = ( 1.0 + xi*( 18.0 - 27.0*xi )) * 0.0625;
662 const real da1 = ( -27.0 - xi*( 18.0 - 81.0*xi )) * 0.0625;
663 const real da2 = ( 27.0 - xi*( 18.0 + 81.0*xi )) * 0.0625;
664 const real da3 = ( -1.0 + xi*( 18.0 + 27.0*xi )) * 0.0625;
665
666 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 )*0.0625;
667 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) )*0.0625;
668 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) )*0.0625;
669 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 )*0.0625;
670
671 N[ 0 ] = da0*b0;
672 N[ 1 ] = da3*b0;
673 N[ 2 ] = da3*b3;
674 N[ 3 ] = da0*b3;
675 N[ 4 ] = da1*b0;
676 N[ 5 ] = da2*b0;
677 N[ 6 ] = da3*b1;
678 N[ 7 ] = da3*b2;
679 N[ 8 ] = da2*b3;
680 N[ 9 ] = da1*b3;
681 N[ 10 ] = da0*b2;
682 N[ 11 ] = da0*b1;
683 N[ 12 ] = da1*b1;
684 N[ 13 ] = da2*b1;
685 N[ 14 ] = da2*b2;
686 N[ 15 ] = da1*b2;
687 }
688
689
690
691
692 inline void
693 Database::deval_quad16deta( const real xi, const real eta, real * N ) const
694 {
695 // often used parameters
696 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 ) * 0.0625;
697 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) ) * 0.0625;
698 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) ) * 0.0625;
699 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 ) * 0.0625;
700
701 const real db0 = ( 1.0 + eta*( 18.0 - 27.0*eta )) * 0.0625;
702 const real db1 = ( -27.0 - eta*( 18.0 - 81.0*eta )) * 0.0625;
703 const real db2 = ( 27.0 - eta*( 18.0 + 81.0*eta )) * 0.0625;
704 const real db3 = ( -1.0 + eta*( 18.0 + 27.0*eta )) * 0.0625;
705
706 // populate output matrix
707 N[ 0 ] = a0*db0;
708 N[ 1 ] = a3*db0;
709 N[ 2 ] = a3*db3;
710 N[ 3 ] = a0*db3;
711 N[ 4 ] = a1*db0;
712 N[ 5 ] = a2*db0;
713 N[ 6 ] = a3*db1;
714 N[ 7 ] = a3*db2;
715 N[ 8 ] = a2*db3;
716 N[ 9 ] = a1*db3;
717 N[ 10 ] = a0*db2;
718 N[ 11 ] = a0*db1;
719 N[ 12 ] = a1*db1;
720 N[ 13 ] = a2*db1;
721 N[ 14 ] = a2*db2;
722 N[ 15 ] = a1*db2;
723
724 }
725
726
727 inline void
728 Database::deval_hex8dxi( const real xi, const real eta, const real zeta, real * N ) const
729 {
730 N[0] = -( eta - 1 ) * ( zeta - 1 ) * 0.125;
731 N[1] = ( eta - 1 ) * ( zeta - 1 ) * 0.125;
732 N[2] = -( eta + 1 ) * ( zeta - 1 ) * 0.125;
733 N[3] = ( eta + 1 ) * ( zeta - 1 ) * 0.125;
734 N[4] = ( eta - 1 ) * ( zeta + 1 ) * 0.125;
735 N[5] = -( eta - 1 ) * ( zeta + 1 ) * 0.125;
736 N[6] = ( eta + 1 ) * ( zeta + 1 ) * 0.125;
737 N[7] = -( eta + 1 ) * ( zeta + 1 ) * 0.125;
738 }
739
740 inline void
741 Database::deval_hex8deta( const real xi, const real eta, const real zeta, real * N ) const
742 {
743 N[ 0 ] = -( xi - 1 ) * ( zeta - 1 ) * 0.125;
744 N[ 1 ] = ( xi + 1 ) * ( zeta - 1 ) * 0.125;
745 N[ 2 ] = -( xi + 1 ) * ( zeta - 1 ) * 0.125;
746 N[ 3 ] = ( xi - 1 ) * ( zeta - 1 ) * 0.125;
747 N[ 4 ] = ( xi - 1 ) * ( zeta + 1 ) * 0.125;
748 N[ 5 ] = -( xi + 1 ) * ( zeta + 1 ) * 0.125;
749 N[ 6 ] = ( xi + 1 ) * ( zeta + 1 ) * 0.125;
750 N[ 7 ] = -( xi - 1 ) * ( zeta + 1 ) * 0.125;
751 }
752
753 inline void
754 Database::deval_hex8dzeta( const real xi, const real eta, const real zeta, real * N ) const
755 {
756 N[ 0 ] = -( eta - 1 ) * ( xi - 1 ) * 0.125;
757 N[ 1 ] = ( eta - 1 ) * ( xi + 1 ) * 0.125;
758 N[ 2 ] = -( eta + 1 ) * ( xi + 1 ) * 0.125;
759 N[ 3 ] = ( eta + 1 ) * ( xi - 1 ) * 0.125;
760 N[ 4 ] = ( eta - 1 ) * ( xi - 1 ) * 0.125;
761 N[ 5 ] = -( eta - 1 ) * ( xi + 1 ) * 0.125;
762 N[ 6 ] = ( eta + 1 ) * ( xi + 1 ) * 0.125;
763 N[ 7 ] = -( eta + 1 ) * ( xi - 1 ) * 0.125;
764 }
765
766 inline void
767 Database::deval_hex27dxi( const real xi, const real eta, const real zeta, real * N ) const
768 {
769 const real eta2 = eta*eta;
770 const real zeta2 = zeta*zeta;
771
772 const real a = 0.125*eta*zeta;
773 const real b = 0.125*xi*zeta;
774 const real c = 0.125*xi*eta;
775 const real d = -0.5*xi*eta*zeta;
776
777 N[ 0 ] = a * ( 2.0 * xi - 1.0 ) * ( eta - 1.0 ) * ( zeta - 1.0 );
778 N[ 1 ] = a * ( 2.0 * xi + 1.0 ) * ( eta - 1.0 ) * ( zeta - 1.0 );
779 N[ 2 ] = a * ( 2.0 * xi + 1.0 ) * ( eta + 1.0 ) * ( zeta - 1.0 );
780 N[ 3 ] = a * ( 2.0 * xi - 1.0 ) * ( eta + 1.0 ) * ( zeta - 1.0 );
781 N[ 4 ] = a * ( 2.0 * xi - 1.0 ) * ( eta - 1.0 ) * ( zeta + 1.0 );
782 N[ 5 ] = a * ( 2.0 * xi + 1.0 ) * ( eta - 1.0 ) * ( zeta + 1.0 );
783 N[ 6 ] = a * ( 2.0 * xi + 1.0 ) * ( eta + 1.0 ) * ( zeta + 1.0 );
784 N[ 7 ] = a * ( 2.0 * xi - 1.0 ) * ( eta + 1.0 ) * ( zeta + 1.0 );
785 N[ 8 ] = d * ( eta - 1.0 ) * ( zeta - 1.0 );
786 N[ 9 ] = - ( zeta * ( eta2 - 1.0 ) * ( 2.0 * xi + 1.0 ) * ( zeta - 1.0 ) ) * 0.25;
787 N[ 10 ] = d * ( eta + 1.0 ) * ( zeta - 1.0 );
788 N[ 11 ] = - ( zeta * ( eta2 - 1.0 ) * ( 2.0 * xi - 1.0 ) * ( zeta - 1.0 ) ) * 0.25;
789 N[ 12 ] = - ( eta * ( 2.0 * xi - 1.0 ) * ( zeta2 - 1.0 ) * ( eta - 1.0 ) ) * 0.25;
790 N[ 13 ] = - ( eta * ( 2.0 * xi + 1.0 ) * ( zeta2 - 1.0 ) * ( eta - 1.0 ) ) * 0.25;
791 N[ 14 ] = - ( eta * ( 2.0 * xi + 1.0 ) * ( zeta2 - 1.0 ) * ( eta + 1.0 ) ) * 0.25;
792 N[ 15 ] = - ( eta * ( 2.0 * xi - 1.0 ) * ( zeta2 - 1.0 ) * ( eta + 1.0 ) ) * 0.25;
793 N[ 16 ] = d * ( eta - 1.0 ) * ( zeta + 1.0 );
794 N[ 17 ] = - ( zeta * ( eta2 - 1.0 ) * ( 2.0 * xi + 1.0 ) * ( zeta + 1.0 ) ) * 0.25;
795 N[ 18 ] = d * ( eta + 1.0 ) * ( zeta + 1.0 );
796 N[ 19 ] = - ( zeta * ( eta2 - 1.0 ) * ( 2.0 * xi - 1.0 ) * ( zeta + 1.0 ) ) * 0.25;
797 N[ 20 ] = - 2.0 * xi * ( eta2 - 1.0 ) * ( zeta2 - 1.0 );
798 N[ 21 ] = 8.0 * b * ( eta2 - 1.0 ) * ( zeta - 1.0 );
799 N[ 22 ] = 8.0 * b * ( eta2 - 1.0 ) * ( zeta + 1.0 );
800 N[ 23 ] = ( ( eta2 - 1.0 ) * ( 2.0 * xi - 1.0 ) * ( zeta2 - 1.0 ) ) * 0.5;
801 N[ 24 ] = ( ( eta2 - 1.0 ) * ( 2.0 * xi + 1.0 ) * ( zeta2 - 1.0 ) ) * 0.5;
802 N[ 25 ] = 8.0 * c * ( zeta2 - 1.0 ) * ( eta - 1.0 );
803 N[ 26 ] = 8.0 * c * ( zeta2 - 1.0 ) * ( eta + 1.0 );
804 }
805
806 inline void
807 Database::deval_hex27deta( const real xi, const real eta, const real zeta, real * N ) const
808 {
809 const real xi2 = xi * xi ;
810 const real zeta2 = zeta*zeta;
811
812 const real a = 0.125*eta*zeta;
813 const real b = 0.125*xi*zeta;
814 const real c = 0.125*xi*eta;
815 const real d = -0.5*xi*eta*zeta;
816
817 N[ 0 ] = b * ( 2.0 * eta - 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 );
818 N[ 1 ] = b * ( 2.0 * eta - 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 );
819 N[ 2 ] = b * ( 2.0 * eta + 1.0 ) * ( xi + 1.0 ) * ( zeta - 1.0 );
820 N[ 3 ] = b * ( 2.0 * eta + 1.0 ) * ( xi - 1.0 ) * ( zeta - 1.0 );
821 N[ 4 ] = b * ( 2.0 * eta - 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 );
822 N[ 5 ] = b * ( 2.0 * eta - 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 );
823 N[ 6 ] = b * ( 2.0 * eta + 1.0 ) * ( xi + 1.0 ) * ( zeta + 1.0 );
824 N[ 7 ] = b * ( 2.0 * eta + 1.0 ) * ( xi - 1.0 ) * ( zeta + 1.0 );
825 N[ 8 ] = - ( zeta * ( 2.0 * eta - 1.0 ) * ( xi2 - 1.0 ) * ( zeta - 1.0 ) ) * 0.25;
826 N[ 9 ] = d * ( xi + 1.0 ) * ( zeta - 1.0 );
827 N[ 10 ] = - ( zeta * ( 2.0 * eta + 1.0 ) * ( xi2 - 1.0 ) * ( zeta - 1.0 ) ) * 0.25;
828 N[ 11 ] = d * ( xi - 1.0 ) * ( zeta - 1.0 );
829 N[ 12 ] = - ( xi * ( 2.0 * eta - 1.0 ) * ( zeta2 - 1.0 ) * ( xi - 1.0 ) ) * 0.25;
830 N[ 13 ] = - ( xi * ( 2.0 * eta - 1.0 ) * ( zeta2 - 1.0 ) * ( xi + 1.0 ) ) * 0.25;
831 N[ 14 ] = - ( xi * ( 2.0 * eta + 1.0 ) * ( zeta2 - 1.0 ) * ( xi + 1.0 ) ) * 0.25;
832 N[ 15 ] = - ( xi * ( 2.0 * eta + 1.0 ) * ( zeta2 - 1.0 ) * ( xi - 1.0 ) ) * 0.25;
833 N[ 16 ] = - ( zeta * ( 2.0 * eta - 1.0 ) * ( xi2 - 1.0 ) * ( zeta + 1.0 ) ) * 0.25;
834 N[ 17 ] = d * ( xi + 1.0 ) * ( zeta + 1.0 );
835 N[ 18 ] = - ( zeta * ( 2.0 * eta + 1.0 ) * ( xi2 - 1.0 ) * ( zeta + 1.0 ) ) * 0.25;
836 N[ 19 ] = d * ( xi - 1.0 ) * ( zeta + 1.0 );
837 N[ 20 ] = - 2.0 * eta * ( xi2 - 1.0 ) * ( zeta2 - 1.0 );
838 N[ 21 ] = 8.0 * a * ( xi2 - 1.0 ) * ( zeta - 1.0 );
839 N[ 22 ] = 8.0 * a * ( xi2 - 1.0 ) * ( zeta + 1.0 );
840 N[ 23 ] = 8.0 * c * ( zeta2 - 1.0 ) * ( xi - 1.0 );
841 N[ 24 ] = 8.0 * c * ( zeta2 - 1.0 ) * ( xi + 1.0 );
842 N[ 25 ] = ( ( 2.0 * eta - 1.0 ) * ( xi2 - 1.0 ) * ( zeta2 - 1.0 ) ) * 0.5;
843 N[ 26 ] = ( ( 2.0 * eta + 1.0 ) * ( xi2 - 1.0 ) * ( zeta2 - 1.0 ) ) * 0.5;
844 }
845
846 inline void
847 Database::deval_hex27dzeta( const real xi, const real eta, const real zeta, real * N ) const
848 {
849 const real xi2 = xi * xi ;
850 const real eta2 = eta*eta;
851
852 const real a = 0.125*eta*zeta;
853 const real b = 0.125*xi*zeta;
854 const real c = 0.125*xi*eta;
855 const real d = -0.5*xi*eta*zeta;
856
857 N[ 0 ] = c * ( 2.0 * zeta - 1.0 ) * ( eta - 1.0 ) * ( xi - 1.0 );
858 N[ 1 ] = c * ( 2.0 * zeta - 1.0 ) * ( eta - 1.0 ) * ( xi + 1.0 );
859 N[ 2 ] = c * ( 2.0 * zeta - 1.0 ) * ( eta + 1.0 ) * ( xi + 1.0 );
860 N[ 3 ] = c * ( 2.0 * zeta - 1.0 ) * ( eta + 1.0 ) * ( xi - 1.0 );
861 N[ 4 ] = c * ( 2.0 * zeta + 1.0 ) * ( eta - 1.0 ) * ( xi - 1.0 );
862 N[ 5 ] = c * ( 2.0 * zeta + 1.0 ) * ( eta - 1.0 ) * ( xi + 1.0 );
863 N[ 6 ] = c * ( 2.0 * zeta + 1.0 ) * ( eta + 1.0 ) * ( xi + 1.0 );
864 N[ 7 ] = c * ( 2.0 * zeta + 1.0 ) * ( eta + 1.0 ) * ( xi - 1.0 );
865 N[ 8 ] = - ( eta * ( xi2 - 1.0 ) * ( 2.0 * zeta - 1.0 ) * ( eta - 1.0 ) ) * 0.25;
866 N[ 9 ] = - ( xi * ( eta2 - 1.0 ) * ( 2.0 * zeta - 1.0 ) * ( xi + 1.0 ) ) * 0.25;
867 N[ 10 ] = - ( eta * ( xi2 - 1.0 ) * ( 2.0 * zeta - 1.0 ) * ( eta + 1.0 ) ) * 0.25;
868 N[ 11 ] = - ( xi * ( eta2 - 1.0 ) * ( 2.0 * zeta - 1.0 ) * ( xi - 1.0 ) ) * 0.25;
869 N[ 12 ] = d * ( eta - 1.0 ) * ( xi - 1.0 );
870 N[ 13 ] = d * ( eta - 1.0 ) * ( xi + 1.0 );
871 N[ 14 ] = d * ( eta + 1.0 ) * ( xi + 1.0 );
872 N[ 15 ] = d * ( eta + 1.0 ) * ( xi - 1.0 );
873 N[ 16 ] = - ( eta * ( xi2 - 1.0 ) * ( 2.0 * zeta + 1.0 ) * ( eta - 1.0 ) ) * 0.25;
874 N[ 17 ] = - ( xi * ( eta2 - 1.0 ) * ( 2.0 * zeta + 1.0 ) * ( xi + 1.0 ) ) * 0.25;
875 N[ 18 ] = - ( eta * ( xi2 - 1.0 ) * ( 2.0 * zeta + 1.0 ) * ( eta + 1.0 ) ) * 0.25;
876 N[ 19 ] = - ( xi * ( eta2 - 1.0 ) * ( 2.0 * zeta + 1.0 ) * ( xi - 1.0 ) ) * 0.25;
877 N[ 20 ] = - 2.0 * zeta * ( eta2 - 1.0 ) * ( xi2 - 1.0 );
878 N[ 21 ] = ( ( eta2 - 1.0 ) * ( xi2 - 1.0 ) * ( 2.0 * zeta - 1.0 ) ) * 0.5;
879 N[ 22 ] = ( ( eta2 - 1.0 ) * ( xi2 - 1.0 ) * ( 2.0 * zeta + 1.0 ) ) * 0.5;
880 N[ 23 ] = 8.0 * b * ( eta2 - 1.0 ) * ( xi - 1.0 );
881 N[ 24 ] = 8.0 * b * ( eta2 - 1.0 ) * ( xi + 1.0 );
882 N[ 25 ] = 8.0 * a * ( xi2 - 1.0 ) * ( eta - 1.0 );
883 N[ 26 ] = 8.0 * a * ( xi2 - 1.0 ) * ( eta + 1.0 );
884
885 }
886
887 inline void
888 Database::deval_hex64dxi( const real xi, const real eta, const real zeta, real * N ) const
889 {
890 const real da0 = ( 1.0 + xi*( 18.0 - 27.0*xi )) * 0.0625;
891 const real da1 = ( -27.0 - xi*( 18.0 - 81.0*xi )) * 0.0625;
892 const real da2 = ( 27.0 - xi*( 18.0 + 81.0*xi )) * 0.0625;
893 const real da3 = ( -1.0 + xi*( 18.0 + 27.0*xi )) * 0.0625;
894
895 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 )*0.0625;
896 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) )*0.0625;
897 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) )*0.0625;
898 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 )*0.0625;
899
900 const real c0 = ( zeta*( 1.0 + 9.0 * zeta * ( 1.0 - zeta ) ) - 1.0 )*0.0625;
901 const real c1 = ( 9.0 - zeta * ( 27.0 + zeta*( 9.0 - 27.0*zeta ) ) )*0.0625;
902 const real c2 = ( 9.0 + zeta * ( 27.0 - zeta*( 9.0 + 27.0*zeta ) ) )*0.0625;
903 const real c3 = ( -zeta*( 1.0 - 9.0 * zeta * ( 1.0 + zeta ) ) - 1.0 )*0.0625;
904
905 N[ 0 ] = da0 * b0 * c0;
906 N[ 1 ] = da3 * b0 * c0;
907 N[ 2 ] = da3 * b3 * c0;
908 N[ 3 ] = da0 * b3 * c0;
909 N[ 4 ] = da0 * b0 * c3;
910 N[ 5 ] = da3 * b0 * c3;
911 N[ 6 ] = da3 * b3 * c3;
912 N[ 7 ] = da0 * b3 * c3;
913 N[ 8 ] = da1 * b0 * c0;
914 N[ 9 ] = da2 * b0 * c0;
915 N[ 10 ] = da0 * b1 * c0;
916 N[ 11 ] = da0 * b2 * c0;
917 N[ 12 ] = da0 * b0 * c1;
918 N[ 13 ] = da0 * b0 * c2;
919 N[ 14 ] = da3 * b1 * c0;
920 N[ 15 ] = da3 * b2 * c0;
921 N[ 16 ] = da3 * b0 * c1;
922 N[ 17 ] = da3 * b0 * c2;
923 N[ 18 ] = da2 * b3 * c0;
924 N[ 19 ] = da1 * b3 * c0;
925 N[ 20 ] = da3 * b3 * c1;
926 N[ 21 ] = da3 * b3 * c2;
927 N[ 22 ] = da0 * b3 * c1;
928 N[ 23 ] = da0 * b3 * c2;
929 N[ 24 ] = da1 * b0 * c3;
930 N[ 25 ] = da2 * b0 * c3;
931 N[ 26 ] = da0 * b1 * c3;
932 N[ 27 ] = da0 * b2 * c3;
933 N[ 28 ] = da3 * b1 * c3;
934 N[ 29 ] = da3 * b2 * c3;
935 N[ 30 ] = da2 * b3 * c3;
936 N[ 31 ] = da1 * b3 * c3;
937 N[ 32 ] = da1 * b1 * c0;
938 N[ 33 ] = da1 * b2 * c0;
939 N[ 34 ] = da2 * b2 * c0;
940 N[ 35 ] = da2 * b1 * c0;
941 N[ 36 ] = da1 * b0 * c1;
942 N[ 37 ] = da2 * b0 * c1;
943 N[ 38 ] = da2 * b0 * c2;
944 N[ 39 ] = da1 * b0 * c2;
945 N[ 40 ] = da0 * b1 * c1;
946 N[ 41 ] = da0 * b1 * c2;
947 N[ 42 ] = da0 * b2 * c2;
948 N[ 43 ] = da0 * b2 * c1;
949 N[ 44 ] = da3 * b1 * c1;
950 N[ 45 ] = da3 * b2 * c1;
951 N[ 46 ] = da3 * b2 * c2;
952 N[ 47 ] = da3 * b1 * c2;
953 N[ 48 ] = da2 * b3 * c1;
954 N[ 49 ] = da1 * b3 * c1;
955 N[ 50 ] = da1 * b3 * c2;
956 N[ 51 ] = da2 * b3 * c2;
957 N[ 52 ] = da1 * b1 * c3;
958 N[ 53 ] = da2 * b1 * c3;
959 N[ 54 ] = da2 * b2 * c3;
960 N[ 55 ] = da1 * b2 * c3;
961 N[ 56 ] = da1 * b1 * c1;
962 N[ 57 ] = da2 * b1 * c1;
963 N[ 58 ] = da2 * b2 * c1;
964 N[ 59 ] = da1 * b2 * c1;
965 N[ 60 ] = da1 * b1 * c2;
966 N[ 61 ] = da2 * b1 * c2;
967 N[ 62 ] = da2 * b2 * c2;
968 N[ 63 ] = da1 * b2 * c2;
969 }
970
971 inline void
972 Database::deval_hex64deta( const real xi, const real eta, const real zeta, real * N ) const
973 {
974 // often used parameters
975 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 ) * 0.0625;
976 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) ) * 0.0625;
977 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) ) * 0.0625;
978 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 ) * 0.0625;
979
980 const real c0 = ( zeta*( 1.0 + 9.0 * zeta * ( 1.0 - zeta ) ) - 1.0 )*0.0625;
981 const real c1 = ( 9.0 - zeta * ( 27.0 + zeta*( 9.0 - 27.0*zeta ) ) )*0.0625;
982 const real c2 = ( 9.0 + zeta * ( 27.0 - zeta*( 9.0 + 27.0*zeta ) ) )*0.0625;
983 const real c3 = ( -zeta*( 1.0 - 9.0 * zeta * ( 1.0 + zeta ) ) - 1.0 )*0.0625;
984
985 const real db0 = ( 1.0 + eta*( 18.0 - 27.0*eta )) * 0.0625;
986 const real db1 = ( -27.0 - eta*( 18.0 - 81.0*eta )) * 0.0625;
987 const real db2 = ( 27.0 - eta*( 18.0 + 81.0*eta )) * 0.0625;
988 const real db3 = ( -1.0 + eta*( 18.0 + 27.0*eta )) * 0.0625;
989
990 N[ 0 ] = a0*c0*db0;
991 N[ 1 ] = a3*c0*db0;
992 N[ 2 ] = a3*c0*db3;
993 N[ 3 ] = a0*c0*db3;
994 N[ 4 ] = a0*c3*db0;
995 N[ 5 ] = a3*c3*db0;
996 N[ 6 ] = a3*c3*db3;
997 N[ 7 ] = a0*c3*db3;
998 N[ 8 ] = a1*c0*db0;
999 N[ 9 ] = a2*c0*db0;
1000 N[ 10 ] = a0*c0*db1;
1001 N[ 11 ] = a0*c0*db2;
1002 N[ 12 ] = a0*c1*db0;
1003 N[ 13 ] = a0*c2*db0;
1004 N[ 14 ] = a3*c0*db1;
1005 N[ 15 ] = a3*c0*db2;
1006 N[ 16 ] = a3*c1*db0;
1007 N[ 17 ] = a3*c2*db0;
1008 N[ 18 ] = a2*c0*db3;
1009 N[ 19 ] = a1*c0*db3;
1010 N[ 20 ] = a3*c1*db3;
1011 N[ 21 ] = a3*c2*db3;
1012 N[ 22 ] = a0*c1*db3;
1013 N[ 23 ] = a0*c2*db3;
1014 N[ 24 ] = a1*c3*db0;
1015 N[ 25 ] = a2*c3*db0;
1016 N[ 26 ] = a0*c3*db1;
1017 N[ 27 ] = a0*c3*db2;
1018 N[ 28 ] = a3*c3*db1;
1019 N[ 29 ] = a3*c3*db2;
1020 N[ 30 ] = a2*c3*db3;
1021 N[ 31 ] = a1*c3*db3;
1022 N[ 32 ] = a1*c0*db1;
1023 N[ 33 ] = a1*c0*db2;
1024 N[ 34 ] = a2*c0*db2;
1025 N[ 35 ] = a2*c0*db1;
1026 N[ 36 ] = a1*c1*db0;
1027 N[ 37 ] = a2*c1*db0;
1028 N[ 38 ] = a2*c2*db0;
1029 N[ 39 ] = a1*c2*db0;
1030 N[ 40 ] = a0*c1*db1;
1031 N[ 41 ] = a0*c2*db1;
1032 N[ 42 ] = a0*c2*db2;
1033 N[ 43 ] = a0*c1*db2;
1034 N[ 44 ] = a3*c1*db1;
1035 N[ 45 ] = a3*c1*db2;
1036 N[ 46 ] = a3*c2*db2;
1037 N[ 47 ] = a3*c2*db1;
1038 N[ 48 ] = a2*c1*db3;
1039 N[ 49 ] = a1*c1*db3;
1040 N[ 50 ] = a1*c2*db3;
1041 N[ 51 ] = a2*c2*db3;
1042 N[ 52 ] = a1*c3*db1;
1043 N[ 53 ] = a2*c3*db1;
1044 N[ 54 ] = a2*c3*db2;
1045 N[ 55 ] = a1*c3*db2;
1046 N[ 56 ] = a1*c1*db1;
1047 N[ 57 ] = a2*c1*db1;
1048 N[ 58 ] = a2*c1*db2;
1049 N[ 59 ] = a1*c1*db2;
1050 N[ 60 ] = a1*c2*db1;
1051 N[ 61 ] = a2*c2*db1;
1052 N[ 62 ] = a2*c2*db2;
1053 N[ 63 ] = a1*c2*db2;
1054 }
1055
1056 inline void
1057 Database::deval_hex64dzeta( const real xi, const real eta, const real zeta, real * N ) const
1058 {
1059 // often used parameters
1060 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 ) * 0.0625;
1061 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) ) * 0.0625;
1062 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) ) * 0.0625;
1063 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 ) * 0.0625;
1064
1065 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 ) * 0.0625;
1066 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) ) * 0.0625;
1067 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) ) * 0.0625;
1068 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 ) * 0.0625;
1069
1070 const real dc0 = ( 1.0 + zeta*( 18.0 - 27.0*zeta )) * 0.0625;
1071 const real dc1 = ( -27.0 - zeta*( 18.0 - 81.0*zeta )) * 0.0625;
1072 const real dc2 = ( 27.0 - zeta*( 18.0 + 81.0*zeta )) * 0.0625;
1073 const real dc3 = ( -1.0 + zeta*( 18.0 + 27.0*zeta )) * 0.0625;
1074
1075 N[ 0 ] = a0*b0*dc0;
1076 N[ 1 ] = a3*b0*dc0;
1077 N[ 2 ] = a3*b3*dc0;
1078 N[ 3 ] = a0*b3*dc0;
1079 N[ 4 ] = a0*b0*dc3;
1080 N[ 5 ] = a3*b0*dc3;
1081 N[ 6 ] = a3*b3*dc3;
1082 N[ 7 ] = a0*b3*dc3;
1083 N[ 8 ] = a1*b0*dc0;
1084 N[ 9 ] = a2*b0*dc0;
1085 N[ 10 ] = a0*b1*dc0;
1086 N[ 11 ] = a0*b2*dc0;
1087 N[ 12 ] = a0*b0*dc1;
1088 N[ 13 ] = a0*b0*dc2;
1089 N[ 14 ] = a3*b1*dc0;
1090 N[ 15 ] = a3*b2*dc0;
1091 N[ 16 ] = a3*b0*dc1;
1092 N[ 17 ] = a3*b0*dc2;
1093 N[ 18 ] = a2*b3*dc0;
1094 N[ 19 ] = a1*b3*dc0;
1095 N[ 20 ] = a3*b3*dc1;
1096 N[ 21 ] = a3*b3*dc2;
1097 N[ 22 ] = a0*b3*dc1;
1098 N[ 23 ] = a0*b3*dc2;
1099 N[ 24 ] = a1*b0*dc3;
1100 N[ 25 ] = a2*b0*dc3;
1101 N[ 26 ] = a0*b1*dc3;
1102 N[ 27 ] = a0*b2*dc3;
1103 N[ 28 ] = a3*b1*dc3;
1104 N[ 29 ] = a3*b2*dc3;
1105 N[ 30 ] = a2*b3*dc3;
1106 N[ 31 ] = a1*b3*dc3;
1107 N[ 32 ] = a1*b1*dc0;
1108 N[ 33 ] = a1*b2*dc0;
1109 N[ 34 ] = a2*b2*dc0;
1110 N[ 35 ] = a2*b1*dc0;
1111 N[ 36 ] = a1*b0*dc1;
1112 N[ 37 ] = a2*b0*dc1;
1113 N[ 38 ] = a2*b0*dc2;
1114 N[ 39 ] = a1*b0*dc2;
1115 N[ 40 ] = a0*b1*dc1;
1116 N[ 41 ] = a0*b1*dc2;
1117 N[ 42 ] = a0*b2*dc2;
1118 N[ 43 ] = a0*b2*dc1;
1119 N[ 44 ] = a3*b1*dc1;
1120 N[ 45 ] = a3*b2*dc1;
1121 N[ 46 ] = a3*b2*dc2;
1122 N[ 47 ] = a3*b1*dc2;
1123 N[ 48 ] = a2*b3*dc1;
1124 N[ 49 ] = a1*b3*dc1;
1125 N[ 50 ] = a1*b3*dc2;
1126 N[ 51 ] = a2*b3*dc2;
1127 N[ 52 ] = a1*b1*dc3;
1128 N[ 53 ] = a2*b1*dc3;
1129 N[ 54 ] = a2*b2*dc3;
1130 N[ 55 ] = a1*b2*dc3;
1131 N[ 56 ] = a1*b1*dc1;
1132 N[ 57 ] = a2*b1*dc1;
1133 N[ 58 ] = a2*b2*dc1;
1134 N[ 59 ] = a1*b2*dc1;
1135 N[ 60 ] = a1*b1*dc2;
1136 N[ 61 ] = a2*b1*dc2;
1137 N[ 62 ] = a2*b2*dc2;
1138 N[ 63 ] = a1*b2*dc2;
1139 }
1140
1141
1142}
1143#endif //BELFEM_CL_DATABASE_HPP
cl_Mesh.hpp
cl_TensorMeshFactory.hpp
cl_Vector.hpp
belfem::Database::evaluate_derivy
real evaluate_derivy(const real aX, const real aY) const
Definition cl_Database.hpp:252
belfem::Database::evaluate_derivz
real evaluate_derivz(const real aX, const real aY, const real aZ) const
Definition cl_Database.hpp:376
belfem::Database::save
void save(const hid_t aHDF5)
Definition cl_Database.cpp:130
belfem::Database::evaluate
real evaluate(const real aX, const real aY) const
Definition cl_Database.hpp:195
belfem::Database::min
real min(const uint aDimension) const
Definition cl_Database.cpp:264
belfem::Database::max
real max(const uint aDimension) const
Definition cl_Database.cpp:270
belfem::Database::evaluate_derivx
real evaluate_derivx(const real aX, const real aY) const
Definition cl_Database.hpp:223
belfem::Database::values
const Vector< real > & values() const
Definition cl_Database.hpp:94
belfem::Database::Database
Database(const string &aFilePath, const string &aMaterial)
Definition cl_Database.cpp:21
belfem::Mesh
Top-level container for all mesh entities.
Definition cl_Mesh.hpp:60
belfem::TensorMeshConfig
Definition cl_TensorMeshConfig.hpp:22
belfem::TensorMeshConfig::num_nodes
uint num_nodes(const uint aDimension=BELFEM_UINT_MAX) const
Definition cl_TensorMeshConfig.hpp:102
hdf5_tools.hpp
belfem::constant::c
const real c
speed of light in m/s ( exact ) http://physics.nist.gov/cgi-bin/cuu/Value?c
Definition constants.hpp:65
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::hid_t
int hid_t
Definition hdf5_types.hpp:20
belfem::uint
unsigned int uint
Definition typedefs.hpp:30
belfem::proc_t
int proc_t
Definition commtypes.hpp:29
belfem::index_t
uint32_t index_t
Definition typedefs.hpp:52
belfem::real
double real
Definition typedefs.hpp:36
corc.tests.test_curve_frame.a
a
Definition test_curve_frame.py:117
corc.tests.test_twist_crosscheck.d
d
Definition test_twist_crosscheck.py:91