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