BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_FEM_Calculator.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#ifndef BELFEM_CL_FEM_CALCULATOR_HPP
12#define BELFEM_CL_FEM_CALCULATOR_HPP
13#include <algorithm>
14#include <cmath>
15
16#include "cl_Material.hpp"
17#include "typedefs.hpp"
18#include "constants.hpp"
19#include "cl_Vector.hpp"
20#include "cl_Matrix.hpp"
21
22#include "fn_dot.hpp"
23#include "fn_det.hpp"
24#include "fn_inv2.hpp"
25#include "fn_inv3.hpp"
26
27#include "cl_Mesh.hpp"
28#include "cl_IF_IntegrationData.hpp"
29#include "en_IWGs.hpp"
30#include "nedelec/cl_EF_EdgeFunction.hpp"
31#include "fn_norm.hpp"
32
33namespace belfem
34{
35 namespace fem
36 {
37 class Element ;
38 class EdgeFunction ;
39 class Group ;
40 class Calculator ;
41 class Kernel ;
42
43//------------------------------------------------------------------------------
44
45 namespace calculator
46 {
47
48 class VectorData
49 {
50 const string mLabel ;
51 const EntityType mType ;
52 uint mIndex = BELFEM_UINT_MAX ;
53
54 Vector< real > mVectorData ;
55 public:
56
57 VectorData( const string & aLabel, const uint aSize, const EntityType aType );
58
59 ~VectorData() = default ;
60
61 const string &
62 label() const ;
63
64 void
65 set_index( const uint aIndex );
66
67 uint
68 index() const ;
69
70 Vector< real > &
71 vector();
72
73 EntityType
74 entity_type() const ;
75
76 };
77
78 class MatrixData
79 {
80 const string mLabel ;
81 uint mIndex = BELFEM_UINT_MAX ;
82 Matrix< real > mMatrixData ;
83
84 public:
85
86 MatrixData( const string & aLabel, const uint aNumRows, const uint aNumCols );
87
88
89 ~MatrixData() = default ;
90
91 const string &
92 label() const ;
93
94 void
95 set_index( const uint aIndex );
96
97 uint
98 index() const ;
99
100 Matrix< real > &
101 matrix();
102
103 };
104
105 enum class MaxwellDataValue
106 {
107 T = 0,
108 H = 1,
109 B = 2,
110 j = 3,
111 rho = 4,
112 cp = 5,
113 lambda = 6,
114 beta = 7,
115 normH = 8,
116 normB = 9,
117 normJ = 10,
118 n = 11, // <-- normal
119 drhodj = 12,
120 drhodT = 13,
121 dcpdT = 14,
122 dlambdadT = 15,
123 x = 16,
124 mu = 17,
125 dmudh = 18,
126 drhodb = 19,
127 drhodbeta = 20,
128 UNDEFINED = 21
129 };
130
131 class MaxwellData
132 {
133 Calculator * mMaxwellCalculator = nullptr ;
134 Calculator * mThermalCalculator = nullptr ;
135
136 Material * mMaterial = nullptr ;
137
138 Matrix< real > mCoords ;
139 real mX = 0.0 ;
140 real mY = 0.0 ;
141 real mZ = 0.0 ;
142
143 const real & mTime ;
144
145 Vector< real > & mH ;
146 Vector< real > & mHn ;
147 Vector< real > & mHt ;
148
149 Vector< real > & mB ;
150 Vector< real > & mBn ;
151 Vector< real > & mBt ;
152
153 Vector< real > & mJ ;
154 Vector< real > & mN ;
155
156 real mRho = BELFEM_QUIET_NAN ;
157 real mdRhodJ = BELFEM_QUIET_NAN ;
158 real mdRhodT = BELFEM_QUIET_NAN ;
159 real mdRhodB = BELFEM_QUIET_NAN ;
160 real mdRhodBeta = BELFEM_QUIET_NAN ;
161 real mBeta = 0.5 * constant::pi ;
162
163 real mT = BELFEM_QUIET_NAN ;
164
165 real mCp = BELFEM_QUIET_NAN ;
166 real mdCpdT = BELFEM_QUIET_NAN ;
167 real mLambda = BELFEM_QUIET_NAN ;
168 real mdLambdadT = BELFEM_QUIET_NAN ;
169
170 real mNormH = BELFEM_QUIET_NAN ;
171 real mNormB = BELFEM_QUIET_NAN ;
172 real mNormJ = BELFEM_QUIET_NAN ;
173
174 real mMu = constant::mu0 ;
175 real mdMudH = 0 ;
176
180 bool mTClamped = false ;
181
184 bool mRhoClamped = false ;
185
188 real mTmax = BELFEM_REAL_MAX ;
189
195 real mDensity = BELFEM_QUIET_NAN ;
196
197 real ( MaxwellData::*mFunT )( const uint aIndex ) = nullptr ;
198
199
200 real ( MaxwellData::*mFunRho )( const uint aIndex ) = nullptr ;
201 real ( MaxwellData::*mFundRhodT )( const uint aIndex ) = nullptr ;
202 real ( MaxwellData::*mFundRhodJ )( const uint aIndex ) = nullptr ;
203 real ( MaxwellData::*mFundRhodB )( const uint aIndex ) = nullptr ;
204 real ( MaxwellData::*mFundRhodBeta )( const uint aIndex ) = nullptr ;
205
208 real ( MaxwellData::*mFunHeat )( const uint aIndex ) = nullptr ;
209
210 real ( MaxwellData::*mFunMu )( const uint aIndex ) = nullptr ;
211 real ( MaxwellData::*mFundMudH )( const uint aIndex ) = nullptr ;
212
213 real ( MaxwellData::*mFunLambda )( const uint aIndex ) = nullptr ;
214 real ( MaxwellData::*mFundLambdadT )( const uint aIndex ) = nullptr ;
215
216 void ( MaxwellData::*mFunX )( const uint aIndex ) = nullptr ;
217 const Vector< real > & ( MaxwellData::*mFunH )( const uint aIndex ) = nullptr ;
218 const Vector< real > & ( MaxwellData::*mFunB )( const uint aIndex ) = nullptr ;
219
220 Vector< uint > mLastIndex ;
221
225 Calculator * mReferenceCalc = nullptr ;
226 bool mFrameCurrent = false ;
229 bool mHaveSeamT = false ;
232 Vector< real > * mBinomialVec = nullptr ;
233 Vector< real > * mTseamVec = nullptr ;
235 Vector< real > mWork ;
236
237 public:
238
239 MaxwellData( Calculator * aCalculator,
240 Kernel * aMaxwellKernel,
241 Kernel * aThermalKernel ) ;
242
243 ~MaxwellData() = default ;
244
245 const Vector< real > &
246 compute_h( const uint aIndex ) ;
247
248 const Vector< real > &
249 compute_b( const uint aIndex ) ;
250
251 const Vector< real > &
252 compute_j( const uint aIndex ) ;
253
254 void
255 compute_x( const uint aIndex ) ;
256
257 real
258 compute_rho( const uint aIndex ) ;
259
260 real
261 compute_drhodT( const uint aIndex ) ;
262
263 real
264 compute_drhodj( const uint aIndex ) ;
265
266 real
267 compute_drhodb( const uint aIndex ) ;
268
269 real
270 compute_drhodbeta( const uint aIndex ) ;
271
272 real
273 compute_T( const uint aIndex ) ;
274
275 real
276 compute_cp( const uint aIndex ) ;
277
278 real
279 compute_dcpdT( const uint aIndex ) ;
280
281 real
282 compute_lambda( const uint aIndex ) ;
283
284 real
285 compute_dlambdadT( const uint aIndex ) ;
286
287 void
288 reset();
289
290 real
291 norm_b( const uint aIndex );
292
293 real
294 norm_j( const uint aIndex );
295
296 real
297 density() const ;
298
299 bool
300 T_clamped() const ;
301
302 bool
303 rho_clamped() const ;
304
305 Calculator *
306 maxwell() ;
307
308 Calculator *
309 thermal() ;
310
311 real
312 compute_mu( const uint aIndex ) ;
313
314 real
315 compute_dmudh( const uint aIndex ) ;
316
321 real
322 compute_volumetric_heatload( const uint aIndex ) ;
323
324 private:
325
326 Vector< real > &
327 link_vector( const string & aLabel );
328
329 bool
330 is_current( MaxwellDataValue aValue, const uint aIndex ) const ;
331
332 void
333 compute_x_2d( const uint aIndex );
334
335 void
336 compute_x_3d( const uint aIndex );
337
338 void
339 set( MaxwellDataValue aValue, const uint aIndex );
340
341 real
342 compute_T_fem( const uint aIndex ) ;
343
344 real
345 compute_T_const( const uint aIndex ) ;
346
349 real
350 compute_T_side_connector( const uint aIndex ) ;
351
352 const Vector< real > &
353 compute_h_bulk_edge( const uint aIndex ) ;
354
355 const Vector< real > &
356 compute_h_ts_edge( const uint aIndex ) ;
357
358 const Vector< real > &
359 compute_h_bulk_node( const uint aIndex ) ;
360
362 const Vector< real > &
363 compute_h_side_connector( const uint aIndex ) ;
364
367 void
368 prepare_side_connector_frame() ;
369
371 void
372 side_connector_xi_eta( const real aPsi, real & aXi, real & aEta ) const ;
373
374 real
375 compute_mu_0( const uint aIndex );
376
377 real
378 compute_mu_const( const uint aIndex );
379
380 real
381 compute_mu_h( const uint aIndex );
382
383 const Vector< real > &
384 compute_b_bulk( const uint aIndex ) ;
385
386 const Vector< real > &
387 compute_b_ts( const uint aIndex ) ;
388
389 real
390 compute_lambda_bulk( const uint aIndex );
391
392 real
393 compute_dlambdadT_bulk( const uint aIndex );
394
395 real
396 compute_lambda_metal( const uint aIndex );
397
398 real
399 compute_dlambdadT_metal( const uint aIndex );
400
401
402 real
403 compute_rho_bulk( const uint aIndex );
404
405 real
406 compute_drhodT_bulk( const uint aIndex );
407
408 real
409 compute_rho_metal( const uint aIndex );
410
411 real
412 compute_drhodT_metal( const uint aIndex );
413
414 real
415 compute_drhodb_metal( const uint aIndex );
416
417 real
418 compute_drhodbeta_metal( const uint aIndex );
419
427 real
428 compute_drhodb_powerlaw_ts( const uint aIndex );
429
430 real
431 compute_drhodb_powerlaw_ts_defect( const uint aIndex );
432
433 real
434 compute_drhodb_piecewise_ts( const uint aIndex );
435
436 real
437 compute_drhodb_piecewise_ts_defect( const uint aIndex );
438
439 real
440 compute_drhodb_powerlaw_bulk( const uint aIndex );
441
442 real
443 compute_drhodb_powerlaw_bulk_defect( const uint aIndex );
444
445 real
446 compute_drhodb_piecewise_bulk( const uint aIndex );
447
448 real
449 compute_drhodb_piecewise_bulk_defect( const uint aIndex );
450
457 real
458 compute_drhodT_powerlaw_ts( const uint aIndex );
459
460 real
461 compute_drhodT_powerlaw_ts_defect( const uint aIndex );
462
463 real
464 compute_drhodT_piecewise_ts( const uint aIndex );
465
466 real
467 compute_drhodT_piecewise_ts_defect( const uint aIndex );
468
469 real
470 compute_drhodT_powerlaw_bulk( const uint aIndex );
471
472 real
473 compute_drhodT_powerlaw_bulk_defect( const uint aIndex );
474
475 real
476 compute_drhodT_piecewise_bulk( const uint aIndex );
477
478 real
479 compute_drhodT_piecewise_bulk_defect( const uint aIndex );
480
483 real
484 compute_rho_riva_ts( const uint aIndex );
485
486 real
487 compute_rho_riva_ts_defect( const uint aIndex );
488
489 real
490 compute_rho_riva_bulk( const uint aIndex );
491
492 real
493 compute_rho_riva_bulk_defect( const uint aIndex );
494
495 real
496 compute_drhodj_riva_ts( const uint aIndex );
497
498 real
499 compute_drhodj_riva_ts_defect( const uint aIndex );
500
501 real
502 compute_drhodj_riva_bulk( const uint aIndex );
503
504 real
505 compute_drhodj_riva_bulk_defect( const uint aIndex );
506
507 real
508 compute_drhodb_riva_ts( const uint aIndex );
509
510 real
511 compute_drhodb_riva_ts_defect( const uint aIndex );
512
513 real
514 compute_drhodb_riva_bulk( const uint aIndex );
515
516 real
517 compute_drhodb_riva_bulk_defect( const uint aIndex );
518
519 real
520 compute_drhodT_riva_ts( const uint aIndex );
521
522 real
523 compute_drhodT_riva_ts_defect( const uint aIndex );
524
525 real
526 compute_drhodT_riva_bulk( const uint aIndex );
527
528 real
529 compute_drhodT_riva_bulk_defect( const uint aIndex );
530
531 real
532 compute_rho_powerlaw_bulk( const uint aIndex );
533
534 real
535 compute_drhodj_powerlaw_bulk( const uint aIndex );
536
537 real
538 compute_rho_powerlaw_ts( const uint aIndex );
539
540 real
541 compute_drhodj_powerlaw_ts( const uint aIndex );
542
543 real
544 compute_rho_piecewise_bulk( const uint aIndex );
545
546 real
547 compute_drhodj_piecewise_bulk( const uint aIndex );
548
549 real
550 compute_rho_piecewise_ts( const uint aIndex );
551
552 real
553 compute_drhodj_piecewise_ts( const uint aIndex );
554
555 real
556 compute_rho_powerlaw_bulk_defect( const uint aIndex );
557
558 real
559 compute_drhodj_powerlaw_bulk_defect( const uint aIndex );
560
561 real
562 compute_rho_powerlaw_ts_defect( const uint aIndex );
563
564 real
565 compute_drhodj_powerlaw_ts_defect( const uint aIndex );
566
567 real
568 compute_rho_piecewise_bulk_defect( const uint aIndex );
569
570 real
571 compute_drhodj_piecewise_bulk_defect( const uint aIndex );
572
573 real
574 compute_rho_piecewise_ts_defect( const uint aIndex );
575
576 real
577 compute_drhodj_piecewise_ts_defect( const uint aIndex );
578
579 real
580 return_zero( const uint aIndex );
581
582 real
583 beta_dummy() const ;
584
585 real
586 compute_dmu_zero( const uint aIndex ) ;
587
588 real
589 compute_dmu_material( const uint aIndex ) ;
590
591 real
592 compute_heatload_user( const uint aIndex ) ;
593
594 };
595 } // end namespace calculator
596
597//------------------------------------------------------------------------------
598
599 class Calculator
600 {
601 // link to group
602 Group * mGroup = nullptr ;
603
604 // link to mesh
605 Mesh * mMesh = nullptr ;
606
607 const ModelDimensionality mDimensionality ;
608
610 const real & mTimestep ;
611
612 // link to current element
613 Element * mElement = nullptr ;
614
615 EdgeFunction * mEdgeFunction = nullptr ;
616 EdgeFunction * mEdgeFunctionMaster = nullptr ;
617 EdgeFunction * mEdgeFunctionSlave = nullptr ;
618
619 Cell< EdgeFunction * > mEdgeFunctionsMaster ;
620 Cell< EdgeFunction * > mEdgeFunctionsSlave ;
621
622 // switch telling if we are allocated
623 bool mIsAllocated = false ;
624
625 uint mNumberOfNodes = BELFEM_UINT_MAX ;
626 uint mNumberOfCornerNodes = BELFEM_UINT_MAX ;
627 uint mNumberOfIntegrationPoints = 0 ;
628
629 uint mNumberOfNodesOnMaster = BELFEM_UINT_MAX ;
630 uint mNumberOfNodesOnSlave = BELFEM_UINT_MAX ;
631
632 uint mNumberOfEdgesOnMaster = BELFEM_UINT_MAX ;
633 uint mNumberOfEdgesOnSlave = BELFEM_UINT_MAX ;
634
635 uint mNumberOfFacesOnMaster = BELFEM_UINT_MAX ;
636 uint mNumberOfFacesOnSlave = BELFEM_UINT_MAX ;
637 uint mMasterIndex = BELFEM_UINT_MAX ;
638
639 IntegrationData * mDomainIntegration = nullptr ;
640 IntegrationData * mLinearIntegration = nullptr ;
641 IntegrationData * mThinShellIntegration = nullptr ;
642 const IntegrationData * mMasterIntegration = nullptr ;
643 const IntegrationData * mSlaveIntegration = nullptr ;
644
645 // for enrichment
646 const IntegrationData * mMasterVolumeIntegration = nullptr ;
647 const IntegrationData * mVolumeEnrichment = nullptr ;
648 const IntegrationData * mSideSetEnrichment = nullptr ;
649
651 bool mIsLinear = false ;
652
654 bool mIsCurved = false ;
655
657 Matrix< real > mK ;
658
660 Matrix< real > mM ;
661
663 Matrix< real > mJN ;
664
666 Vector< real > mf ;
667
669 Vector< real > mq0 ;
670
672 Vector< real > mq ;
673
675 Vector< real > mqswap ;
676
678 Vector< real > mNormal ;
679 uint mNormalIndex = BELFEM_UINT_MAX ;
680 real mSurfaceIncrement = BELFEM_QUIET_NAN ;
681
683 Cell< string > mDofLabels ;
684
686 Cell< calculator::VectorData * > mVectors ;
687 Cell< calculator::MatrixData * > mMatrices ;
688
690 Map< string , calculator::VectorData * > mVectorMap ;
691 Map< string , calculator::MatrixData * > mMatrixMap ;
692
694 Vector< real > mModelParameters ;
695
696 // pointers for faster access
697 Matrix< real > mX ; // node coordinates
698 Matrix< real > mXc ; // node coordinates at corners
699 calculator::MatrixData * mJ = nullptr ; // jacobian matrix
700 calculator::MatrixData * mInvJ = nullptr ; // inverse of the jacobian matrix
701 calculator::MatrixData * mN = nullptr ; // node interpolatoin operator
702 calculator::MatrixData * mdN = nullptr ; // derivative function
703 calculator::MatrixData * mB = nullptr ; // gradient operator
704
705 // for faces
706 Matrix< real > mXm ;
707 calculator::MatrixData * mNm = nullptr ; // node interpolation operator
708 calculator::MatrixData * mJm = nullptr ;
709 calculator::MatrixData * mInvJm = nullptr ;
710 calculator::MatrixData * mBm = nullptr ;
711 calculator::MatrixData * mntBm = nullptr ; // trans( normal ) * B
712 calculator::MatrixData * mnxBm = nullptr ; // cross( normal, B )
713
714 Matrix< real > mXs ;
715 calculator::MatrixData * mNs = nullptr ; // node interpolation operator
716 calculator::MatrixData * mJs = nullptr ;
717 calculator::MatrixData * mInvJs = nullptr ;
718 calculator::MatrixData * mBs = nullptr ;
719 calculator::MatrixData * mntBs = nullptr ; // trans( normal ) * B
720 calculator::MatrixData * mnxBs = nullptr ; // cross( normal, B )
721
722 calculator::MaxwellData * mMaxwellData = nullptr ;
723
724 // kernels for the maxwell data helper, set via link_maxwell();
725 // construction of the helper is deferred to allocate()
726 Kernel * mMaxwellKernel = nullptr ;
727 Kernel * mThermalKernel = nullptr ;
728
729 real mDetJ = BELFEM_QUIET_NAN ;
730 uint mDetJIndex = BELFEM_UINT_MAX ;
731 real mRadius = BELFEM_QUIET_NAN ; // only needed if axisymmetric
732
733 // function to compute the flattened slave integration index
734 // (cumulative facet/orientation offset)
735 uint ( Calculator::*mFunSlaveIntegrationIndex )( const mesh::Facet * aFacet ) ;
736
737 real ( * mFunInvertJ )( const Matrix< real > & aA, Matrix< real > & aB );
738
739 const Vector< real > & ( Calculator::*mFunNormal )( const uint aIndex );
740
741 // function for node interpolator
742 const Matrix< real > & ( Calculator::*mFunN )( const uint aIndex );
743
744 // function for gradient operator
745 const Matrix< real > & ( Calculator::*mFunB )( const uint aIndex );
746
747 // function for node interpolator master
748 const Matrix< real > & ( Calculator::*mFunNm )( const uint aIndex );
749
750 // function for gradient operator master
751 const Matrix< real > & ( Calculator::*mFunBm )( const uint aIndex );
752
753 // function for node interpolator slave
754 const Matrix< real > & ( Calculator::*mFunNs )( const uint aIndex );
755
756 // function for gradient operator slave
757 const Matrix< real > & ( Calculator::*mFunBs )( const uint aIndex );
758
759 // volume increment
760 real ( Calculator::*mFundV )( const uint aIndex );
761
762 // surface increment
763 real ( Calculator::*mFundS )( const uint aIndex );
764
765 // inverse the jacobian
766 const Matrix< real > & ( Calculator::*mFunInvJ )( const uint aIndex );
767
768 const Vector< real > &
769 ( Calculator::*mFunCollectNodeData )( const string & aLabel );
770
771 real
772 ( Calculator::*mFunBJAngle )(
773 const Vector< real > & b,
774 const Vector< real > & j,
775 real & norm_b,
776 real & norm_j ) const ;
777
778 // for theta method
779 real mTheta = 1.0 ;
780 real mOneMinusTheta = 0.0 ;
781
782 const Vector< real > & ( Calculator::*mFunNedelecDataH )();
783 const Vector< real > & ( Calculator::*mFunNedelecDataA )();
784
785 // defaulted so that no allocate() path can ever leave it null
786 void ( Calculator::*mFunLinkElement )( Element * aElement )
787 = & Calculator::link_element_default ;
788
789 uint mIntegrationOrder = 0 ;
790
791 Cell< Vector< real > * > mQold ;
792 index_t mMaxDofFieldIndex = 0 ;
793
794//------------------------------------------------------------------------------
795 public:
796//------------------------------------------------------------------------------
797
799 Calculator( Group * aGroup, const ModelDimensionality aDimensionality );
800
801 // special constrictor used for TET10 T-Matrices
802 Calculator( Group * aGroup, Mesh * aMesh );
803
804//------------------------------------------------------------------------------
805
806 ~Calculator() ;
807
808//------------------------------------------------------------------------------
809
810 void
811 init_qold_table();
812
813//------------------------------------------------------------------------------
814
818 void
819 allocate();
820
821//------------------------------------------------------------------------------
822
823 void
824 link( Group * aGroup );
825
826//------------------------------------------------------------------------------
827
828 void
829 link( Element * aElement );
830
831//------------------------------------------------------------------------------
832
833 void
834 link( mesh::Facet * aFacet );
835
836//------------------------------------------------------------------------------
837
838 Group *
839 group() ;
840
841//------------------------------------------------------------------------------
842
843 Element *
844 element() ;
845
846//------------------------------------------------------------------------------
847
848 const Material *
849 material() const;
850
851//------------------------------------------------------------------------------
852
856 Matrix< real > &
857 K() ;
858
859//------------------------------------------------------------------------------
860
864 Matrix< real > &
865 M() ;
866
867//------------------------------------------------------------------------------
868
872 Matrix< real > &
873 JN() ;
874
875//------------------------------------------------------------------------------
876
880 Vector< real > &
881 f() ;
882
883//------------------------------------------------------------------------------
884
888 const Vector< real > &
889 q() ;
890
891//------------------------------------------------------------------------------
892
896 const Vector< real > &
897 qold( const uint aStep=0 ) ;
898
899//------------------------------------------------------------------------------
900
904 Vector< real > &
905 qswap() ;
906
907//------------------------------------------------------------------------------
908
912 const Vector< real > &
913 node_data( const string & aNodeField );
914
915//------------------------------------------------------------------------------
916
917 bool
918 vector_exists( const string aLabel ) const ;
919
920//------------------------------------------------------------------------------
921
925 Vector< real > &
926 vector( const string & aLabel );
927
928//------------------------------------------------------------------------------
929
930 bool
931 matrix_exists( const string aLabel ) const ;
932
933//------------------------------------------------------------------------------
934
938 Matrix< real > &
939 matrix( const string & aLabel );
940
941//------------------------------------------------------------------------------
942
946 const Matrix< real > &
947 X() const ;
948
949//------------------------------------------------------------------------------
950
954 const Matrix< real > &
955 J( const uint aIndex ) ;
956
957//------------------------------------------------------------------------------
958
962 const Matrix< real > &
963 invJ( const uint aIndex ) ;
964
965//------------------------------------------------------------------------------
966
970 const Matrix< real > &
971 Jm( const uint aIndex ) ;
972
973//------------------------------------------------------------------------------
974
978 const Matrix< real > &
979 Js( const uint aIndex ) ;
980
981//------------------------------------------------------------------------------
982
986 const Matrix< real > &
987 N( const uint aIndex ) ;
988
989//------------------------------------------------------------------------------
990
994 const Matrix< real > &
995 E( const uint aIndex ) ;
996
997//------------------------------------------------------------------------------
998
1002 const Matrix< real > &
1003 C( const uint aIndex ) ;
1004
1005//------------------------------------------------------------------------------
1006
1010 const Matrix< real > &
1011 G( const uint aIndex ) ;
1012
1013//------------------------------------------------------------------------------
1014
1015 uint
1016 num_nedelec_dofs() const ;
1017
1018//------------------------------------------------------------------------------
1019
1023 const Matrix< real > &
1024 Em( const uint aIndex ) ;
1025
1026//------------------------------------------------------------------------------
1027
1031 const Matrix< real > &
1032 Cm( const uint aIndex ) ;
1033
1034//------------------------------------------------------------------------------
1035
1039 const Matrix< real > &
1040 Es( const uint aIndex ) ;
1041
1042//------------------------------------------------------------------------------
1043
1047 const Matrix< real > &
1048 Cs( const uint aIndex ) ;
1049
1050//------------------------------------------------------------------------------
1051
1055 const Matrix< real > &
1056 Nm( const uint aIndex ) ;
1057
1058//------------------------------------------------------------------------------
1059
1063 const Matrix< real > &
1064 Ns( const uint aIndex ) ;
1065
1066//------------------------------------------------------------------------------
1067
1071 const Vector< real > &
1072 Nvec( const uint aIndex ) const ;
1073
1074//------------------------------------------------------------------------------
1075
1079 real
1080 node_interp( const uint aIndex, const Vector< real > & aNodeValues ) const;
1081
1082//------------------------------------------------------------------------------
1083
1087 const Matrix< real > &
1088 B( const uint aIndex=0 ) ;
1089
1090//------------------------------------------------------------------------------
1091
1095 const Matrix< real > &
1096 Bm( const uint aIndex=0 ) ;
1097
1098//------------------------------------------------------------------------------
1099
1103 const Matrix< real > &
1104 Bs( const uint aIndex=0 ) ;
1105
1106//------------------------------------------------------------------------------
1107
1111 const Matrix< real > &
1112 Xm() const ;
1113
1114//------------------------------------------------------------------------------
1115
1119 const Matrix< real > &
1120 Xs() const ;
1121
1122//------------------------------------------------------------------------------
1123
1127 real
1128 dS ( const uint aIndex=0 ) ;
1129
1130//------------------------------------------------------------------------------
1131
1135 real
1136 dV ( const uint aIndex=0 ) ;
1137
1138//------------------------------------------------------------------------------
1139
1145 const Vector< real > &
1146 normal( const uint aIndex=0 );
1147
1148//------------------------------------------------------------------------------
1149
1150 void
1151 initialize_integration( const ElementType aElementType,
1152 const InterpolationType aInterpolationType );
1153
1154//------------------------------------------------------------------------------
1155
1156 void
1157 set_integration_order( const uint aOrder );
1158
1159 uint
1160 integration_order() const ;
1161
1162//------------------------------------------------------------------------------
1163
1164 void
1165 allocate_memory();
1166
1167//------------------------------------------------------------------------------
1168
1169 const IntegrationData *
1170 integration() const ;
1171
1172//------------------------------------------------------------------------------
1173
1174 const IntegrationData *
1175 master_integration() const ;
1176
1177//------------------------------------------------------------------------------
1178
1179 const IntegrationData *
1180 slave_integration() const ;
1181
1182//------------------------------------------------------------------------------
1183
1184 const IntegrationData *
1185 volume_integration() const ;
1186
1187//------------------------------------------------------------------------------
1188
1189 const IntegrationData *
1190 volume_enrichment() const ;
1191
1192//------------------------------------------------------------------------------
1193
1194 const IntegrationData *
1195 sideset_enrichment() const ;
1196
1197//------------------------------------------------------------------------------
1198
1199 uint
1200 num_intpoints() const ;
1201
1202//------------------------------------------------------------------------------
1203
1204 real
1205 timestep() const ;
1206
1207 void
1208 set_model_parameters( const Vector< real > & aParams );
1209
1210//------------------------------------------------------------------------------
1211
1212 const Vector< real > &
1213 model_parameters() const ;
1214
1215//------------------------------------------------------------------------------
1216
1217 calculator::VectorData *
1218 create_vector( const string & aLabel,
1219 const uint aSize,
1220 const EntityType aType = EntityType::UNDEFINED );
1221
1222//------------------------------------------------------------------------------
1223
1224 calculator::MatrixData *
1225 create_matrix( const string & aLabel,
1226 const uint aNumRows,
1227 const uint aNumCols );
1228
1229//------------------------------------------------------------------------------
1230
1231 void
1232 print_dofs();
1233
1234//------------------------------------------------------------------------------
1235
1236 void
1237 print_local_dofs();
1238
1239//------------------------------------------------------------------------------
1240
1241 const Vector< real > &
1242 nedelec_data_h();
1243
1244//------------------------------------------------------------------------------
1245
1246 const Vector< real > &
1247 nedelec_data_a();
1248
1249//------------------------------------------------------------------------------
1250
1258 Calculator *
1259 get_normal_calculator(
1260 Vector< real > & aPhiM,
1261 Vector< real > & aPhiS,
1262 bool & aMasterIsConductor,
1263 bool & aSlaveIsConductor ) ;
1264
1265//------------------------------------------------------------------------------
1266
1273 const Vector< real > &
1274 nedelec_data_master_h();
1275
1276//------------------------------------------------------------------------------
1277
1282 const Vector< real > &
1283 nedelec_data_slave_h();
1284
1285//------------------------------------------------------------------------------
1286
1292 bool
1293 volume_is_conductor( const mesh::Element * aVolume ) const ;
1294
1295//------------------------------------------------------------------------------
1296
1297 Mesh *
1298 mesh() ;
1299
1300//------------------------------------------------------------------------------
1301
1302 bool
1303 element_is_linear() const ;
1304
1305//------------------------------------------------------------------------------
1306
1314 void
1315 link_maxwell( Kernel * aMaxwellKernel,
1316 Kernel * aThermalKernel = nullptr ) ;
1317
1318//------------------------------------------------------------------------------
1319
1320 // helper function for material functions
1321 // computes the angle between magnetic field and current density
1322 real
1323 bj_angle( const Vector< real > & b, const Vector< real > & j,
1324 real & norm_b, real & norm_j ) const ;
1325
1326 // helper function for material functions
1327 // computes the unfolded angle [ 0, pi ] between magnetic field and tape normal
1328 real
1329 bn_angle( const Vector< real > & b, const Vector< real > & n,
1330 real & norm_b ) const ;
1331
1332 calculator::MaxwellData *
1333 maxwell() ;
1334
1335//------------------------------------------------------------------------------
1336
1337 // getter needed by sideset connector
1338 EdgeFunction *
1339 edge_function()
1340 {
1341 return mEdgeFunction ;
1342 }
1343
1344//------------------------------------------------------------------------------
1345 private:
1346//------------------------------------------------------------------------------
1347
1348 // picks the mFunLinkElement dispatcher based on the registered
1349 // kernels; called from allocate() and again from link_maxwell()
1350 // when the thermal kernel is attached after allocation
1351 void
1352 select_link_element_dispatcher();
1353
1354 void
1355 link_element_default( Element * aElement );
1356
1357 void
1358 link_element_maxwell( Element * aElement );
1359
1360 void
1361 link_element_maxwell_thermal( Element * aElement );
1362
1363 void
1364 link_element_thermal_maxwell( Element * aElement );
1365
1366//------------------------------------------------------------------------------
1367
1368 // Dispatchers for the flattened slave integration index
1369 // (cumulative over facet, orientation). The corresponding
1370 // IntegrationData and EdgeFunction are looked up by this index
1371 // at Calculator::link() time.
1372
1373 uint
1374 slave_integration_index_2d( const mesh::Facet * aFacet ) ;
1375
1376 uint
1377 slave_integration_index_tet( const mesh::Facet * aFacet ) ;
1378
1379 uint
1380 slave_integration_index_hex( const mesh::Facet * aFacet ) ;
1381
1382 uint
1383 slave_integration_index_penta( const mesh::Facet * aFacet ) ;
1384
1385//------------------------------------------------------------------------------
1386
1387 // node interpolator for scalar fields
1388 const Matrix< real > &
1389 Nscalar( const uint aIndex ) ;
1390
1391//------------------------------------------------------------------------------
1392
1393 // node interpolator for scalar fields, master
1394 const Matrix< real > &
1395 Nscalar_master( const uint aIndex ) ;
1396
1397//------------------------------------------------------------------------------
1398
1399 // node interpolator for scalar fields, slave
1400 const Matrix< real > &
1401 Nscalar_slave( const uint aIndex ) ;
1402
1403//------------------------------------------------------------------------------
1404
1405 // node interpolator for 2D vector fields
1406 const Matrix< real > &
1407 N2D( const uint aIndex ) ;
1408
1409//------------------------------------------------------------------------------
1410
1411 // node interpolator for 3D vector fields
1412 const Matrix< real > &
1413 N3D( const uint aIndex ) ;
1414
1415//------------------------------------------------------------------------------
1416
1417 // gradient operator for scalar fields
1418 const Matrix< real > &
1419 Bscalar( const uint aIndex ) ;
1420
1421//------------------------------------------------------------------------------
1422
1423 // gradient operator for scalar fields master
1424 const Matrix< real > &
1425 Bscalar_master( const uint aIndex ) ;
1426
1427//------------------------------------------------------------------------------
1428
1429 // gradient operator for scalar fields slave
1430 const Matrix< real > &
1431 Bscalar_slave( const uint aIndex ) ;
1432
1433//------------------------------------------------------------------------------
1434
1435 // gradient operator for plane stress
1436 const Matrix< real > &
1437 Bplanestress( const uint aIndex ) ;
1438
1439//------------------------------------------------------------------------------
1440
1441 // gradient operator for 3d mech
1442 const Matrix< real > &
1443 Bvoigt( const uint aIndex ) ;
1444
1445//------------------------------------------------------------------------------
1446
1447 const Vector< real > &
1448 normal_tri_straight( const uint aIndex );
1449
1450//------------------------------------------------------------------------------
1451
1452 const Vector< real > &
1453 normal_tri_curved( const uint aIndex );
1454
1455//------------------------------------------------------------------------------
1456
1457 const Vector< real > &
1458 normal_quad_straight( const uint aIndex );
1459
1460//------------------------------------------------------------------------------
1461
1462 const Vector< real > &
1463 normal_quad_curved( const uint aIndex );
1464
1465//------------------------------------------------------------------------------
1466
1467 const Vector< real > &
1468 normal_tet_straight( const uint aIndex );
1469
1470//------------------------------------------------------------------------------
1471
1472 const Vector< real > &
1473 normal_tet_curved( const uint aIndex );
1474
1475//------------------------------------------------------------------------------
1476
1477 const Vector< real > &
1478 normal_penta( const uint aIndex );
1479
1480//------------------------------------------------------------------------------
1481
1482 const Vector< real > &
1483 normal_hex( const uint aIndex );
1484
1485//------------------------------------------------------------------------------
1486
1487 real
1488 dV_tri6_tet10( const uint aIndex );
1489
1490//------------------------------------------------------------------------------
1491
1492 real
1493 dV_ts( const uint aIndex );
1494
1495//------------------------------------------------------------------------------
1496
1497 real
1498 dV_quad4ts( const uint aIndex );
1499
1500//------------------------------------------------------------------------------
1501
1502 real
1503 dV_hex( const uint aIndex );
1504
1505//------------------------------------------------------------------------------
1506
1507 real
1508 dV_axsymmx( const uint aIndex );
1509
1510//------------------------------------------------------------------------------
1511
1512 real
1513 dV_axsymmy( const uint aIndex );
1514
1515//------------------------------------------------------------------------------
1516
1517 real
1518 dS_cartesian( const uint aIndex );
1519
1520//------------------------------------------------------------------------------
1521
1522 real
1523 dS_axsymmx( const uint aIndex );
1524
1525//------------------------------------------------------------------------------
1526
1527 real
1528 dS_axsymmy( const uint aIndex );
1529
1530//------------------------------------------------------------------------------
1531
1532 const Vector< real > &
1533 nedelec_data_linear_h();
1534
1535 const Vector< real > &
1536 nedelec_data_quadratic_h_2d();
1537
1538 const Vector< real > &
1539 nedelec_data_quadratic_h_3d();
1540
1541 const Vector< real > &
1542 nedelec_data_linear_a();
1543
1544 const Vector< real > &
1545 nedelec_data_quadratic_a_3d();
1546
1547//------------------------------------------------------------------------------
1548
1552 const Vector< real > &
1553 nedelec_data_linear( const string & aEdgeField );
1554
1555//------------------------------------------------------------------------------
1556
1557
1558//------------------------------------------------------------------------------
1559
1563 const Vector< real > &
1564 nedelec_data_quadratic_2d( const string & aEdgeField,
1565 const string & aFaceField,
1566 const string & aVectorLabel );
1567
1568//------------------------------------------------------------------------------
1569
1573 const Vector< real > &
1574 nedelec_data_quadratic_3d( const string & aEdgeField,
1575 const string & aFaceField,
1576 const string & aVectorLabel );
1577
1578//------------------------------------------------------------------------------
1579
1580 const Matrix< real > &
1581 invJ2D3D( const uint aIndex ) ;
1582
1583//------------------------------------------------------------------------------
1584
1585
1586 const Matrix< real > &
1587 invJaxsym( const uint aIndex ) ;
1588
1589 real
1590 bj_angle_2d( const Vector< real > & b, const Vector< real > & j, real & norm_b, real & norm_j ) const ;
1591
1592 real
1593 bj_angle_3d( const Vector< real > & b, const Vector< real > & j, real & norm_b, real & norm_j ) const ;
1594
1595//------------------------------------------------------------------------------
1596 };
1597//------------------------------------------------------------------------------
1598
1599 inline const string &
1600 calculator::VectorData::label() const
1601 {
1602 return mLabel ;
1603 }
1604
1605 inline void
1606 calculator::VectorData::set_index( const belfem::uint aIndex )
1607 {
1608 mIndex = aIndex ;
1609 }
1610
1611 inline uint
1612 calculator::VectorData::index() const
1613 {
1614 return mIndex ;
1615 }
1616
1617 inline Vector< real > &
1618 calculator::VectorData::vector()
1619 {
1620 return mVectorData ;
1621 }
1622
1623 inline EntityType
1624 calculator::VectorData::entity_type() const
1625 {
1626 return mType ;
1627 }
1628
1629 inline const string &
1630 calculator::MatrixData::label() const
1631 {
1632 return mLabel ;
1633 }
1634
1635 inline void
1636 calculator::MatrixData::set_index( const uint aIndex )
1637 {
1638 mIndex = aIndex ;
1639 }
1640
1641 inline uint
1642 calculator::MatrixData::index() const
1643 {
1644 return mIndex ;
1645 }
1646
1647//------------------------------------------------------------------------------
1648
1649 inline Matrix< real > &
1650 calculator::MatrixData::matrix()
1651 {
1652 return mMatrixData ;
1653 }
1654
1655//------------------------------------------------------------------------------
1656
1657 inline bool Calculator::vector_exists( const string aLabel ) const
1658 {
1659 return mVectorMap.key_exists( aLabel );
1660 }
1661
1662//------------------------------------------------------------------------------
1663
1664 inline Vector< real > &
1665 Calculator::vector( const string & aLabel )
1666 {
1667 return mVectorMap( aLabel )->vector() ;
1668 }
1669
1670//------------------------------------------------------------------------------
1671
1672 inline bool Calculator::matrix_exists( const string aLabel ) const
1673 {
1674 return mMatrixMap.key_exists( aLabel );
1675 }
1676
1677//------------------------------------------------------------------------------
1678
1679 inline Matrix< real > &
1680 Calculator::matrix( const string & aLabel )
1681 {
1682 return mMatrixMap( aLabel )->matrix() ;
1683 }
1684
1685//------------------------------------------------------------------------------
1686
1687 inline const Matrix< real > &
1688 Calculator::X() const
1689 {
1690 return mX ;
1691 }
1692
1693//------------------------------------------------------------------------------
1694
1695 inline const Matrix< real > &
1696 Calculator::Xm() const
1697 {
1698 return mXm ;
1699 }
1700
1701//------------------------------------------------------------------------------
1702
1703 inline const Matrix< real > &
1704 Calculator::Xs() const
1705 {
1706 return mXs ;
1707 }
1708
1709//------------------------------------------------------------------------------
1710
1711 inline const Matrix< real > &
1712 Calculator::J( const uint aIndex )
1713 {
1714 if ( aIndex != mJ->index() )
1715 {
1716 mJ->set_index( aIndex );
1717 mJ->matrix().matrix_data() = mIsCurved ?
1718 mDomainIntegration->dNdXi( aIndex ) * mX :
1719 mLinearIntegration->dNdXi( aIndex ) * mXc ;
1720 }
1721
1722 return mJ->matrix() ;
1723 }
1724
1725//------------------------------------------------------------------------------
1726
1727 inline const Matrix< real > &
1728 Calculator::Jm( const uint aIndex )
1729 {
1730 if ( aIndex != mJm->index() )
1731 {
1732 mJm->set_index( aIndex );
1733 mJm->matrix().matrix_data() =
1734 mMasterIntegration->dNdXi( aIndex ) * mXm ;
1735 }
1736
1737 return mJm->matrix() ;
1738 }
1739
1740//------------------------------------------------------------------------------
1741
1742 inline const Matrix< real > &
1743 Calculator::Js( const uint aIndex )
1744 {
1745 if ( aIndex != mJs->index() )
1746 {
1747 mJs->set_index( aIndex );
1748 mJs->matrix().matrix_data() =
1749 mSlaveIntegration->dNdXi( aIndex ) * mXs ;
1750 }
1751
1752 return mJs->matrix() ;
1753 }
1754
1755//------------------------------------------------------------------------------
1756
1757 inline const Matrix< real > &
1758 Calculator::invJ( const uint aIndex )
1759 {
1760 return ( this->*mFunInvJ )( aIndex );
1761 }
1762
1763//------------------------------------------------------------------------------
1764
1765 inline const Matrix< real > &
1766 Calculator::invJ2D3D( const uint aIndex )
1767 {
1768 if ( aIndex != mInvJ->index() )
1769 {
1770 // remember index
1771 mInvJ->set_index( aIndex );
1772
1773 // in 2D and 3D, we can directly use this value for dV
1774 mDetJIndex = aIndex ;
1775
1776 // compute inverse and remember determinant
1777 mDetJ = ( * mFunInvertJ )( this->J( aIndex ), mInvJ->matrix() );
1778 }
1779
1780 return mInvJ->matrix() ;
1781 }
1782
1783//------------------------------------------------------------------------------
1784
1785 inline const Matrix< real > &
1786 Calculator::invJaxsym( const uint aIndex )
1787 {
1788 if ( aIndex != mInvJ->index() )
1789 {
1790 // remember index
1791 mInvJ->set_index( aIndex );
1792
1793 // just invert the matrix but do not store the determinant
1794 ( * mFunInvertJ )( this->J( aIndex ), mInvJ->matrix() );
1795 }
1796
1797 return mInvJ->matrix() ;
1798 }
1799
1800//------------------------------------------------------------------------------
1801
1802 inline const Matrix< real > &
1803 Calculator::Nscalar( const uint aIndex )
1804 {
1805 return mDomainIntegration->N( aIndex );
1806 }
1807
1808//------------------------------------------------------------------------------
1809
1810 inline const Matrix< real > &
1811 Calculator::Nscalar_master( const uint aIndex )
1812 {
1813 return mMasterIntegration->N( aIndex );
1814 }
1815
1816//------------------------------------------------------------------------------
1817
1818 inline const Matrix< real > &
1819 Calculator::Nscalar_slave( const uint aIndex )
1820 {
1821 return mSlaveIntegration->N( aIndex );
1822 }
1823
1824//------------------------------------------------------------------------------
1825
1826 inline const Matrix< real > &
1827 Calculator::N2D( const uint aIndex )
1828 {
1829 if( mN->index() != aIndex )
1830 {
1831 // precomputed data
1832 const Vector< real > & tPhi = mDomainIntegration->phi( aIndex );
1833
1834 // remember the index
1835 mN->set_index( aIndex );
1836
1837 // get link to matrix
1838 Matrix< real > & tN = mN->matrix() ;
1839
1840 // initialize counter
1841 uint tCount = 0 ;
1842
1843 // populate matrix
1844 for( uint k=0; k<mNumberOfNodes; ++k )
1845 {
1846 tN( 0, tCount++ ) = tPhi( k );
1847 tN( 1, tCount++ ) = tPhi( k );
1848 }
1849 }
1850 return mN->matrix() ;
1851 }
1852
1853//------------------------------------------------------------------------------
1854
1855 inline const Matrix< real > &
1856 Calculator::N3D( const uint aIndex )
1857 {
1858 if( mN->index() != aIndex )
1859 {
1860 // precomputed data
1861 const Vector< real > & tPhi = mDomainIntegration->phi( aIndex );
1862
1863 // remember the index
1864 mN->set_index( aIndex );
1865
1866 // get link to matrix
1867 Matrix< real > & tN = mN->matrix() ;
1868
1869 // initialize counter
1870 uint tCount = 0 ;
1871
1872 // populate matrix
1873 for( uint k=0; k<mNumberOfNodes; ++k )
1874 {
1875 tN( 0, tCount++ ) = tPhi( k );
1876 tN( 1, tCount++ ) = tPhi( k );
1877 tN( 2, tCount++ ) = tPhi( k );
1878 }
1879 }
1880 return mN->matrix() ;
1881 }
1882
1883//------------------------------------------------------------------------------
1884
1885 inline const Matrix< real > &
1886 Calculator::Bscalar( const uint aIndex )
1887 {
1888 if ( aIndex != mB->index() )
1889 {
1890 // remember the index
1891 mB->set_index( aIndex );
1892
1893 mB->matrix() =
1894 this->invJ( aIndex ) * mDomainIntegration->dNdXi( aIndex );
1895 }
1896
1897 return mB->matrix() ;
1898 }
1899
1900
1901//------------------------------------------------------------------------------
1902
1903 inline const Matrix< real > &
1904 Calculator::Bscalar_master( const uint aIndex )
1905 {
1906
1907 if ( aIndex != mBm->index() )
1908 {
1909 // remember the index
1910 mBm->set_index( aIndex );
1911
1912 ( * mFunInvertJ )(this->Jm( aIndex ), mInvJm->matrix() ) ;
1913
1914 mBm->matrix() =mInvJm->matrix() * mMasterIntegration->dNdXi( aIndex );
1915 }
1916
1917 return mBm->matrix() ;
1918 }
1919
1920 inline const Matrix< real > &
1921 Calculator::Bscalar_slave( const uint aIndex )
1922 {
1923 if ( aIndex != mBs->index() )
1924 {
1925 // remember the index
1926 mBs->set_index( aIndex );
1927
1928 ( * mFunInvertJ )(this->Js( aIndex ), mInvJs->matrix() ) ;
1929
1930 mBs->matrix() = mInvJs->matrix() * mSlaveIntegration->dNdXi( aIndex );
1931 }
1932
1933 return mBs->matrix() ;
1934 }
1935
1936//------------------------------------------------------------------------------
1937
1938 inline const Matrix< real > &
1939 Calculator::Bplanestress( const uint aIndex )
1940 {
1941 if ( aIndex != mB->index() )
1942 {
1943 // remember the index
1944 mB->set_index( aIndex );
1945
1946 // compute derivative for scalar field
1947 Matrix< real > & tdN = mdN->matrix() ;
1948
1949 // compute derivatives
1950 tdN = this->invJ( aIndex ) * mDomainIntegration->dNdXi( aIndex );
1951
1952 // get matrix object
1953 Matrix< real > & tB = mB->matrix() ;
1954
1955 // initialize counter
1956 uint tCount = 0 ;
1957
1958 // populate data
1959 for( uint k=0; k<mNumberOfNodes; ++k )
1960 {
1961 tB( 0, tCount ) = tdN( 0, k );
1962 tB( 2, tCount++ ) = tdN( 1, k );
1963 tB( 1, tCount ) = tdN( 1, k );
1964 tB( 2, tCount++ ) = tdN( 0, k );
1965 }
1966 }
1967
1968 return mB->matrix() ;
1969 }
1970
1971//------------------------------------------------------------------------------
1972
1973 inline const Matrix< real > &
1974 Calculator::Bvoigt( const uint aIndex )
1975 {
1976 if ( aIndex != mB->index() )
1977 {
1978 // remember the index
1979 mB->set_index( aIndex );
1980
1981 // compute derivative for scalar field
1982 Matrix< real > & tdN = mdN->matrix() ;
1983
1984 // compute derivatives
1985 tdN = this->invJ( aIndex ) * mDomainIntegration->dNdXi( aIndex );
1986
1987 // get matrix object
1988 Matrix< real > & tB = mB->matrix() ;
1989
1990 // initialize counter
1991 uint tCount = 0 ;
1992
1993 // populate data
1994 for( uint k=0; k<mNumberOfNodes; ++k )
1995 {
1996 tB( 0, tCount ) = tdN( 0, k );
1997 tB( 4, tCount ) = tdN( 2, k );
1998 tB( 5, tCount ) = tdN( 1, k );
1999 ++tCount ;
2000
2001 tB( 1, tCount ) = tdN( 1, k );
2002 tB( 3, tCount ) = tdN( 2, k );
2003 tB( 5, tCount ) = tdN( 0, k );
2004 ++tCount ;
2005
2006 tB( 2, tCount ) = tdN( 2, k );
2007 tB( 3, tCount ) = tdN( 1, k );
2008 tB( 4, tCount ) = tdN( 0, k );
2009 ++tCount ;
2010 }
2011 }
2012
2013 return mB->matrix() ;
2014 }
2015
2016//------------------------------------------------------------------------------
2017
2018 inline const Matrix< real > &
2019 Calculator::N( const uint aIndex )
2020 {
2021 return ( this->*mFunN )( aIndex );
2022 }
2023
2024//------------------------------------------------------------------------------
2025
2026 inline const Matrix< real > &
2027 Calculator::Nm( const uint aIndex )
2028 {
2029 return ( this->*mFunNm )( aIndex );
2030 }
2031
2032//------------------------------------------------------------------------------
2033
2034 inline const Matrix< real > &
2035 Calculator::Ns( const uint aIndex )
2036 {
2037 return ( this->*mFunNs )( aIndex );
2038 }
2039
2040//------------------------------------------------------------------------------
2041
2042 inline const Matrix< real > &
2043 Calculator::Bm( const uint aIndex )
2044 {
2045 return ( this->*mFunBm )( aIndex );
2046 }
2047
2048//------------------------------------------------------------------------------
2049
2050 inline const Matrix< real > &
2051 Calculator::Bs( const uint aIndex )
2052 {
2053 return ( this->*mFunBs )( aIndex );
2054 }
2055
2056//------------------------------------------------------------------------------
2057
2058 inline real
2059 Calculator::node_interp( const uint aIndex, const Vector< real > & aNodeValues ) const
2060 {
2061 return dot( mDomainIntegration->phi( aIndex ), aNodeValues );
2062 }
2063
2064//------------------------------------------------------------------------------
2065 inline const Vector< real > &
2066 Calculator::Nvec( const uint aIndex ) const
2067 {
2068 return mDomainIntegration->phi( aIndex ) ;
2069 }
2070
2071//------------------------------------------------------------------------------
2072
2073 inline const Matrix< real > &
2074 Calculator::B( const uint aIndex )
2075 {
2076 return ( this->*mFunB )( aIndex );
2077 }
2078
2079//------------------------------------------------------------------------------
2080
2081 inline const Vector< real > &
2082 Calculator::normal( const uint aIndex )
2083 {
2084 return ( this->*mFunNormal) ( aIndex );
2085 }
2086
2087//------------------------------------------------------------------------------
2088
2089 inline real
2090 Calculator::dS( const uint aIndex )
2091 {
2092#if !defined( NDEBUG ) || defined( DEBUG )
2093 real adS = ( this->*mFundS )( aIndex );
2094 BELFEM_ASSERT( adS >= 0.0, "Negative Jacobian determinant" );
2095 return adS;
2096#else
2097 return ( this->*mFundS )( aIndex );
2098#endif
2099 }
2100
2101//------------------------------------------------------------------------------
2102
2103 inline real
2104 Calculator::dV( const uint aIndex )
2105 {
2106#if !defined( NDEBUG ) || defined( DEBUG )
2107 real adV = ( this->*mFundV)( aIndex );
2108 BELFEM_ERROR( adV >= 0.0, "Negative Jacobian determinant" );
2109 return adV;
2110#else
2111 return ( this->*mFundV)( aIndex );
2112#endif
2113 }
2114
2115//------------------------------------------------------------------------------
2116
2117 // tri and tet only!
2118 inline real
2119 Calculator::dV_tri6_tet10( const uint aIndex )
2120 {
2121 // we don't need to recompute this at every point for linear elements
2122 uint tIndex = mIsCurved ? aIndex : 0 ;
2123
2124 if( mDetJIndex != tIndex )
2125 {
2126 mDetJIndex = tIndex ;
2127 mDetJ = det( this->J( tIndex ) );
2128 }
2129 return mDetJ ;
2130 }
2131
2132 inline real
2133 Calculator::dV_ts( const uint aIndex )
2134 {
2135 // Formulation discriminator (not element-specific):
2136 // - Magnetic solve: an edge function is attached, and its
2137 // most recent C(aIndex) call has already cached the
2138 // Jacobian determinant at this integration point in
2139 // mDetJ. Reuse it.
2140 // - Thermal solve (or any Nedelec-free path): no edge
2141 // function; fall back to computing det(J) from the
2142 // scalar Lagrange shape on demand.
2143 // Used by PENTA6TS (where the edge function caches
2144 // thickness*surface as the element-level volume) and by
2145 // HEX8TS (where the edge function caches the per-IP
2146 // determinant of the 3D Jacobian); in both cases the
2147 // edge-function value is the correct dV weight.
2148 if (mEdgeFunction != nullptr)
2149 {
2150 return mEdgeFunction->det_J() ;
2151 }
2152 else
2153 {
2154 uint tIndex = mIsCurved ? aIndex : 0 ;
2155
2156 if( mDetJIndex != tIndex )
2157 {
2158 mDetJIndex = tIndex ;
2159 mDetJ = det( this->J( tIndex ) );
2160 }
2161 return mDetJ ;
2162 }
2163 }
2164
2165 inline real
2166 Calculator::dV_quad4ts( const uint aIndex )
2167 {
2168 // For any physically meaningful ( positive ) layer thickness,
2169 // QUAD4TS layer elements are wound CLOCKWISE by construction, so
2170 // their scalar Lagrange determinant is negative. A negative deck
2171 // thickness is currently accepted and inverts the stack, which
2172 // flips that sign; this function is correct either way, but the
2173 // reasoning below assumes the valid case. That winding is
2174 // forced, not accidental: the layer element's local edge node
2175 // order must follow the facet edge direction, because the edge
2176 // signs are read from the node order ( Element::compute_edge_
2177 // directions ) while the Nedelec tangent is read from the facet
2178 // ( EF_QUAD4TS::link ). Given the clockwise extrusion normal of
2179 // ThinShellFactory::process_nodes_line2, ( bottom0, bottom1,
2180 // top1, top0 ) is the only ordering that satisfies both, and it
2181 // is left-handed.
2182 //
2183 // The map is still a diffeomorphism, so the thermal assembly is
2184 // correct up to the volume weight: N is reference-space and
2185 // B = invJ * dNdXi returns true Cartesian gradients whichever
2186 // orientation the map has. Change of variables then asks for
2187 // | det J |, which is what this returns -- an exact weight, not
2188 // a clamp.
2189 //
2190 // The MAGNETIC solve reaches this function too, through its own
2191 // Calculator for the same block -- the dispatch is by element
2192 // type, so do NOT assume magnetic never enters here. It takes
2193 // dV_ts's edge-function branch, where det_J is
2194 // 0.25 * thickness * length: positive for any thickness that
2195 // reaches assembly, so the abs is the identity and magnetic
2196 // results are unchanged. A negative deck thickness would flip
2197 // that, but it cannot get this far -- Kernel::compute_element_
2198 // volumes forms facet_area * thickness per layer element and
2199 // hard-errors on a negative total before any assembly runs.
2200 //
2201 // Wraps dV_ts rather than duplicating its body ON PURPOSE. Both
2202 // functions share the mDetJ / mDetJIndex cache with invJ2D3D,
2203 // which the thermal kernels populate by calling B( k ) BEFORE
2204 // dV( k ). At k = 0 that leaves a valid cache entry, so a copy
2205 // of the body that only wrapped its own mDetJ = det( ... )
2206 // assignment would return the SIGNED cached value and never run
2207 // the abs. Taking it on the returned value is immune to that.
2208 return std::abs( this->dV_ts( aIndex ) );
2209 }
2210
2211//------------------------------------------------------------------------------
2212
2213 inline real
2214 Calculator::dV_hex( const uint aIndex )
2215 {
2216 if ( mEdgeFunction != nullptr )
2217 {
2218 mEdgeFunction->update_nabla( aIndex );
2219 return mEdgeFunction->det_J() ;
2220 }
2221 else
2222 {
2223
2224 if( mDetJIndex != aIndex )
2225 {
2226 mDetJIndex = aIndex ;
2227 mDetJ = det( this->J( aIndex ) );
2228 }
2229 return mDetJ ;
2230 }
2231 }
2232
2233//------------------------------------------------------------------------------
2234
2235 inline real
2236 Calculator::dV_axsymmx( const uint aIndex )
2237 {
2238 if( mDetJIndex != aIndex )
2239 {
2240 mDetJIndex = aIndex ;
2241
2242 // radius contribution
2243 mDetJ = mIsCurved ? dot(
2244 mDomainIntegration->phi( aIndex ).vector_data(), mX.col( 1 ) ) :
2245 dot( mLinearIntegration->phi( aIndex ).vector_data(), mXc.col( 1 ) );
2246
2247 mDetJ *= det( this->J( aIndex ) ) * 2.0 * constant::pi ;
2248 }
2249 return mDetJ ;
2250 }
2251
2252//------------------------------------------------------------------------------
2253
2254 inline real
2255 Calculator::dV_axsymmy( const uint aIndex )
2256 {
2257
2258 if( mDetJIndex != aIndex )
2259 {
2260 mDetJIndex = aIndex ;
2261
2262 // radius contribution
2263 mDetJ = mIsCurved ? dot(
2264 mDomainIntegration->phi( aIndex ).vector_data(), mX.col( 0 ) ) :
2265 dot( mLinearIntegration->phi( aIndex ).vector_data(), mXc.col( 0 ) );
2266
2267 mDetJ *= det( this->J( aIndex ) ) * 2.0 * constant::pi ;
2268 }
2269 return mDetJ ;
2270 }
2271
2272//------------------------------------------------------------------------------
2273
2274
2275 inline real
2276 Calculator::dS_cartesian ( const uint aIndex )
2277 {
2278 // compute the normal if it hasn't been computed so far
2279 ( this->*mFunNormal ) ( aIndex );
2280
2281 // return the surface increment
2282 // note that the surface increment is scaled reciprocally to the
2283 // integration weights.
2284 return mSurfaceIncrement ;
2285 }
2286
2287//------------------------------------------------------------------------------
2288
2289 inline real
2290 Calculator::dS_axsymmx( const uint aIndex )
2291 {
2292 // compute the normal if it hasn't been computed so far
2293 ( this->*mFunNormal ) ( aIndex );
2294
2295 // return the surface increment
2296 return dot( mMasterIntegration->phi( aIndex ).vector_data(),
2297 mXm.col( 1 ) )
2298 * mSurfaceIncrement * 2.0 * constant::pi ;
2299
2300 }
2301
2302//------------------------------------------------------------------------------
2303
2304 inline real
2305 Calculator::dS_axsymmy( const uint aIndex )
2306 {
2307 // compute the normal if it hasn't been computed so far
2308 ( this->*mFunNormal ) ( aIndex );
2309
2310 // return the surface increment
2311 return dot( mMasterIntegration->phi( aIndex ).vector_data(),
2312 mXm.col( 0 ) )
2313 * mSurfaceIncrement * 2.0 * constant::pi ;
2314
2315 }
2316
2317//------------------------------------------------------------------------------
2318
2319 inline const IntegrationData *
2320 Calculator::integration() const
2321 {
2322 return mDomainIntegration ;
2323 }
2324
2325//------------------------------------------------------------------------------
2326
2327 inline const IntegrationData *
2328 Calculator::master_integration() const
2329 {
2330 return mMasterIntegration ;
2331 }
2332
2333//------------------------------------------------------------------------------
2334
2335 inline const IntegrationData *
2336 Calculator::slave_integration() const
2337 {
2338 return mSlaveIntegration ;
2339 }
2340
2341//------------------------------------------------------------------------------
2342
2343 inline const IntegrationData *
2344 Calculator::volume_integration() const
2345 {
2346 return mMasterVolumeIntegration ;
2347 }
2348
2349//------------------------------------------------------------------------------
2350
2351 inline const IntegrationData *
2352 Calculator::volume_enrichment() const
2353 {
2354 return mVolumeEnrichment ;
2355 }
2356
2357//------------------------------------------------------------------------------
2358
2359 inline const IntegrationData *
2360 Calculator::sideset_enrichment() const
2361 {
2362 return mSideSetEnrichment ;
2363 }
2364
2365//------------------------------------------------------------------------------
2366
2367 inline Matrix< real > &
2368 Calculator::K()
2369 {
2370 return mK ;
2371 }
2372
2373//------------------------------------------------------------------------------
2374
2375 inline Matrix< real > &
2376 Calculator::M()
2377 {
2378 return mM ;
2379 }
2380
2381//------------------------------------------------------------------------------
2382
2383 inline Matrix< real > &
2384 Calculator::JN()
2385 {
2386 return mJN ;
2387 }
2388
2389//------------------------------------------------------------------------------
2390
2391 inline Vector< real > &
2392 Calculator::f()
2393 {
2394 return mf ;
2395 }
2396
2397//------------------------------------------------------------------------------
2398
2399 inline uint
2400 Calculator::num_intpoints() const
2401 {
2402 return mNumberOfIntegrationPoints ;
2403 }
2404
2405//------------------------------------------------------------------------------
2406
2407 inline real
2408 Calculator::timestep() const
2409 {
2410 return mTimestep ;
2411 }
2412
2413//------------------------------------------------------------------------------
2414 inline void
2415 Calculator::set_model_parameters( const Vector< real > & aParams )
2416 {
2417 mModelParameters = aParams ;
2418 }
2419
2420//------------------------------------------------------------------------------
2421
2422 inline const Vector< real > &
2423 Calculator::model_parameters() const
2424 {
2425 return mModelParameters ;
2426 }
2427
2428//------------------------------------------------------------------------------
2429
2430 inline Element *
2431 Calculator::element()
2432 {
2433 return mElement ;
2434 }
2435
2436//------------------------------------------------------------------------------
2437
2438 inline Group *
2439 Calculator::group()
2440 {
2441 return mGroup ;
2442 }
2443
2444//------------------------------------------------------------------------------
2445
2446 inline const Vector< real > &
2447 Calculator::nedelec_data_h()
2448 {
2449 return (this->*mFunNedelecDataH )();
2450 }
2451
2452//------------------------------------------------------------------------------
2453
2454 inline const Vector< real > &
2455 Calculator::nedelec_data_a()
2456 {
2457 return (this->*mFunNedelecDataA )();
2458 }
2459
2460//------------------------------------------------------------------------------
2461
2462 inline const Vector< real > &
2463 Calculator::nedelec_data_linear_h()
2464 {
2465 return this->nedelec_data_linear( "edge_h");
2466 }
2467
2468//------------------------------------------------------------------------------
2469
2470 inline const Vector< real > &
2471 Calculator::nedelec_data_linear_a()
2472 {
2473 return this->nedelec_data_linear( "edge_a");
2474 }
2475
2476//------------------------------------------------------------------------------
2477
2478 inline const Vector< real > &
2479 Calculator::nedelec_data_quadratic_h_2d()
2480 {
2481 return this->nedelec_data_quadratic_2d(
2482 "edge_h",
2483 "face_h",
2484 "nedelec_h");
2485 }
2486
2487//------------------------------------------------------------------------------
2488
2489 inline const Vector< real > &
2490 Calculator::nedelec_data_quadratic_h_3d()
2491 {
2492 return this->nedelec_data_quadratic_3d(
2493 "edge_h",
2494 "face_h",
2495 "nedelec_h");
2496 }
2497
2498//------------------------------------------------------------------------------
2499
2500 inline const Vector< real > &
2501 Calculator::nedelec_data_quadratic_a_3d()
2502 {
2503 return this->nedelec_data_quadratic_3d(
2504 "edge_a",
2505 "face_a",
2506 "nedelec_a");
2507 }
2508
2509//------------------------------------------------------------------------------
2510
2511 inline uint
2512 Calculator::integration_order() const
2513 {
2514 return mIntegrationOrder ;
2515 }
2516
2517//------------------------------------------------------------------------------
2518
2519 inline Mesh * Calculator::mesh()
2520 {
2521 return mMesh ;
2522 }
2523
2524//------------------------------------------------------------------------------
2525
2526 inline bool Calculator::element_is_linear() const
2527 {
2528 return mIsLinear ;
2529 }
2530
2531 inline real
2532 Calculator::bj_angle( const Vector< real > & b, const Vector< real > & j, real & norm_b, real & norm_j ) const
2533 {
2534 return ( this->*mFunBJAngle )( b, j, norm_b, norm_j );
2535 }
2536
2537 inline real
2538 Calculator::bn_angle( const Vector< real > & b, const Vector< real > & n, real & norm_b ) const
2539 {
2540 norm_b = norm( b ) ;
2541 BELFEM_ASSERT( std::abs( norm( n ) - 1.0 ) < 1e-6, "tape normal must be normalized" ) ;
2542
2543 // UNFOLDED angle between field and tape normal, [ 0, pi ]
2544 // ( 2026-08-16 ). theta < pi/2 iff b has a component
2545 // along +n; n is the outward normal of the mid-surface facet's
2546 // MASTER volume element ( master -> slave ), which is also the
2547 // layer-stack direction — a minus prefix on the thinshell
2548 // sidesets key flips both together. Measured jc(theta) tables
2549 // are asymmetric about pi/2 and consume this angle as is;
2550 // analytic laws that are even in theta ( ModifiedKim ) fold
2551 // internally by construction
2552 return norm_b < 1e-6 ? constant::pi*0.5 :
2553 std::acos( std::clamp( dot( n, b ) / norm_b, -1.0, 1.0 ) ) ;
2554 }
2555
2556 inline real
2557 Calculator::bj_angle_3d( const Vector< real > & b, const Vector< real > & j, real & norm_b, real & norm_j ) const
2558 {
2559 norm_b = norm( b ) ;
2560 norm_j = norm( j ) ;
2561
2562 // Angle between magnetic field and current density (beta)
2563 return norm_b < 1e-6 || norm_j < 1e-6 ?
2564 constant::pi*0.5 : std::acos( std::min(std::abs( dot( b, j ) / ( norm_b * norm_j ) ) , 1.0 )) ;
2565 }
2566
2567 inline real
2568 Calculator::bj_angle_2d( const Vector< real > & b, const Vector< real > & j, real & norm_b, real & norm_j ) const
2569 {
2570 norm_b = norm( b ) ;
2571 norm_j = norm( j ) ;
2572
2573 // Angle between magnetic field and current density (beta)
2574 // in 2D: current is always perpendicular to model plane
2575 return constant::pi*0.5 ;
2576 }
2577
2578 inline Vector< real > &
2579 calculator::MaxwellData::link_vector( const string & aLabel )
2580 {
2581 if ( ! mMaxwellCalculator->vector_exists( aLabel ) )
2582 {
2583 uint tNumDim = mMaxwellCalculator->mesh()->number_of_dimensions() ;
2584
2585 // for 2D calculations, j is always perpendicular to model plane
2586 // this is the only special case that we have to catch here
2587 if ( aLabel == "j" && tNumDim == 2 )
2588 {
2589 tNumDim = 1 ;
2590 }
2591
2592 VectorData * tData = mMaxwellCalculator->create_vector( aLabel, tNumDim );
2593 return tData->vector();
2594 }
2595 return mMaxwellCalculator->vector( aLabel );
2596 }
2597
2598 inline void
2599 calculator::MaxwellData::compute_x( const uint aIndex )
2600 {
2601 if ( ! this->is_current( MaxwellDataValue::x, aIndex ) )
2602 {
2603 ( this->*mFunX )( aIndex );
2604 this->set( MaxwellDataValue::x, aIndex );
2605 }
2606 }
2607
2608 inline void
2609 calculator::MaxwellData::compute_x_2d( const uint aIndex )
2610 {
2611 mCoords = mMaxwellCalculator->N( aIndex ) * mMaxwellCalculator->X();
2612 mX = mCoords( 0, 0 );
2613 mY = mCoords( 0, 1 );
2614 }
2615
2616 inline void
2617 calculator::MaxwellData::compute_x_3d( const uint aIndex )
2618 {
2619 mCoords = mMaxwellCalculator->N( aIndex ) * mMaxwellCalculator->X();
2620 mX = mCoords( 0, 0 );
2621 mY = mCoords( 0, 1 );
2622 mZ = mCoords( 0, 2 );
2623 }
2624
2625 inline const Vector< real > &
2626 calculator::MaxwellData::compute_j( const uint aIndex )
2627 {
2628 if ( aIndex == mLastIndex( static_cast< uint > ( MaxwellDataValue::j ) ) )
2629 {
2630 return mJ ;
2631 }
2632
2633 mLastIndex( static_cast< uint > ( MaxwellDataValue::j ) ) = aIndex ;
2634 mJ = mMaxwellCalculator->C( aIndex ) * mMaxwellCalculator->q();
2635
2636 return mJ ;
2637 }
2638
2639
2668 inline void
2669 compute_h_trace(
2670 Calculator * aNormalCalc,
2671 const bool aMaster,
2672 const bool aIsConductor,
2673 const Vector< real > & aPhi,
2674 Vector< real > & aScratch,
2675 Vector< real > & aH )
2676 {
2677 if ( aIsConductor )
2678 {
2679 const Vector< real > & q = aMaster ?
2680 aNormalCalc->nedelec_data_master_h() :
2681 aNormalCalc->nedelec_data_slave_h() ;
2682
2683 const Vector< real > & w = aNormalCalc->integration()->weights() ;
2684
2685 const uint tNumPoints = aNormalCalc->num_intpoints() ;
2686
2687 real tSum = 0.0 ;
2688 aH.fill( 0.0 );
2689
2690 for ( uint k = 0; k < tNumPoints; ++k )
2691 {
2692 aScratch = ( aMaster ? aNormalCalc->Em( k ) : aNormalCalc->Es( k ) ) * q ;
2693 aH += w( k ) * aScratch ;
2694 tSum += w( k );
2695 }
2696 aH /= tSum ;
2697 }
2698 else
2699 {
2700 aH = -1. * ( aMaster ? aNormalCalc->Bm( 0 ) : aNormalCalc->Bs( 0 ) ) * aPhi ;
2701 }
2702 }
2703
2712 inline const Vector< real > &
2713 compute_hn( Calculator * aCalc , const uint k )
2714 {
2715 // normal field
2716 Vector< real > & hn = aCalc->vector("hn");
2717
2718 // warning:
2719 // for linear elements, we only need to do this once
2720 BELFEM_ASSERT( aCalc->element_is_linear(), "compute_hn doesn't work for non-linear elements because it expects surface integration points but we integrate over the volume" );
2721
2722 if ( k == 0 )
2723 {
2724 // dofs from master element
2725 Vector< real > & phi_m = aCalc->vector("phi_m");
2726
2727 // dofs from slave element
2728 Vector< real > & phi_s = aCalc->vector("phi_s");
2729
2730 // get calculator for facet and the kind of each volume side
2731 bool tMasterIsConductor ;
2732 bool tSlaveIsConductor ;
2733
2734 Calculator * tCalc = aCalc->get_normal_calculator(
2735 phi_m, phi_s, tMasterIsConductor, tSlaveIsConductor );
2736
2737 // scratch for the per-point products
2738 Vector< real > & hk = aCalc->vector("hk");
2739
2740 // h-field on master side
2741 Vector< real > & hm = aCalc->vector("hm");
2742 compute_h_trace( tCalc, true, tMasterIsConductor, phi_m, hk, hm );
2743
2744 // h-field on slave side
2745 Vector< real > & hs = aCalc->vector("hs");
2746 compute_h_trace( tCalc, false, tSlaveIsConductor, phi_s, hk, hs );
2747
2748 // temporarily writing average into hn container
2749 hn = 0.5 * ( hm + hs ) ;
2750
2751 // normal vector
2752 Vector< real > & n = aCalc->vector("normal");
2753 n = tCalc->normal( k );
2754
2755 // normal component of h
2756 hn = dot( hn, n ) * n ;
2757 }
2758 return hn;
2759 }
2760
2761 inline void
2762 Calculator::link_element_default( Element * aElement )
2763 {
2764 mElement = aElement;
2765 }
2766
2767
2768
2769 inline void
2770 calculator::MaxwellData::reset()
2771 {
2772 mLastIndex.fill( BELFEM_UINT_MAX );
2773 mBeta = 0.5 * constant::pi ;
2774 mTClamped = false ;
2775 mRhoClamped = false ;
2776 mFrameCurrent = false ;
2777 }
2778
2779 inline bool
2780 calculator::MaxwellData::is_current( MaxwellDataValue aValue, const uint aIndex ) const
2781 {
2782 return mLastIndex( static_cast< uint >( aValue ) ) == aIndex ;
2783 }
2784
2785 inline void
2786 calculator::MaxwellData::set( MaxwellDataValue aValue, const uint aIndex )
2787 {
2788 mLastIndex( static_cast< uint >( aValue ) ) = aIndex ;
2789 }
2790
2791 inline const Vector< real > &
2792 calculator::MaxwellData::compute_h( const uint aIndex )
2793 {
2794 if ( ! this->is_current( MaxwellDataValue::H, aIndex ) )
2795 {
2796 this->set( MaxwellDataValue::H, aIndex );
2797 return ( this->*mFunH )( aIndex );
2798
2799 }
2800 return mH ;
2801 }
2802
2803 inline real
2804 calculator::MaxwellData::compute_mu( const uint aIndex )
2805 {
2806 if ( ! this->is_current( MaxwellDataValue::mu, aIndex ) )
2807 {
2808 this->set( MaxwellDataValue::mu, aIndex );
2809 mMu = ( this->*mFunMu )( aIndex );
2810 }
2811 return mMu ;
2812 }
2813
2814 inline real
2815 calculator::MaxwellData::compute_mu_0( const uint aIndex )
2816 {
2817 return constant::mu0 ;
2818 }
2819
2820 inline real
2821 calculator::MaxwellData::compute_mu_const( const uint aIndex )
2822 {
2823 return mMaterial->constant_property( MaterialProperty::mu );
2824 }
2825
2826 inline real
2827 calculator::MaxwellData::compute_mu_h( const uint aIndex )
2828 {
2829 if ( ! this->is_current( MaxwellDataValue::normH, aIndex ) )
2830 {
2831 mNormH = norm( this->compute_h( aIndex ));
2832 this->set( MaxwellDataValue::normH, aIndex );
2833 }
2834 return mMaterial->mu( mNormH );
2835 }
2836
2837 inline const Vector< real > &
2838 calculator::MaxwellData::compute_b( const uint aIndex )
2839 {
2840 if ( ! this->is_current( MaxwellDataValue::B, aIndex ) )
2841 {
2842 this->set( MaxwellDataValue::B, aIndex );
2843 return ( this->*mFunB )( aIndex );
2844 }
2845 return mB ;
2846 }
2847
2848 inline const Vector< real > &
2849 calculator::MaxwellData::compute_b_bulk( const uint aIndex )
2850 {
2851 mB = this->compute_mu( aIndex ) * this->compute_h( aIndex ) ;
2852 return mB ;
2853 }
2854
2855 inline const Vector< real > &
2856 calculator::MaxwellData::compute_b_ts( const uint aIndex )
2857 {
2858 this->compute_h( aIndex );
2859 real mu = this->compute_mu( aIndex );
2860 mBt = mu * mHt ;
2861 mBn = mu * mHn ;
2862 mB = mBt + mBn ;
2863 return mB ;
2864 }
2865
2866 inline real
2867 calculator::MaxwellData::compute_rho( const uint aIndex )
2868 {
2869 if ( ! this->is_current( MaxwellDataValue::rho, aIndex ) )
2870 {
2871 real tRho = ( this->*mFunRho )( aIndex );
2872
2873 // numerical guard on the power-law output: one clamped value
2874 // for all consumers ( K, Joule source, element mean ). While
2875 // clamped, the J-derivative is zero ( consistent tangent )
2876 mRhoClamped = ( tRho < gRhoMin ) || ( tRho > gRhoMax ) ;
2877 mRho = mRhoClamped ? std::clamp( tRho, gRhoMin, gRhoMax ) : tRho ;
2878
2879 this->set( MaxwellDataValue::rho, aIndex );
2880 }
2881 return mRho ;
2882 }
2883
2884
2885 inline real
2886 calculator::MaxwellData::compute_drhodj( const uint aIndex )
2887 {
2888 if ( ! this->is_current( MaxwellDataValue::drhodj, aIndex ) )
2889 {
2890 // make sure rho — and with it mRhoClamped — is current here
2891 this->compute_rho( aIndex );
2892
2893 // consistent tangent: a clamped rho has zero derivative
2894 mdRhodJ = mRhoClamped ? 0.0 : ( this->*mFundRhodJ )( aIndex );
2895 this->set( MaxwellDataValue::drhodj, aIndex );
2896 }
2897 return mdRhodJ ;
2898 }
2899
2900 inline real
2901 calculator::MaxwellData::compute_drhodT( const uint aIndex )
2902 {
2903 if ( ! this->is_current( MaxwellDataValue::drhodT, aIndex ) )
2904 {
2905 // make sure rho — and with it mRhoClamped — is current here
2906 this->compute_rho( aIndex );
2907
2908 // refresh mTClamped for THIS index ( cache-hit if the rho
2909 // law above already computed it ): a cached rho would
2910 // otherwise leave the flag at whatever index last touched
2911 // compute_T — same self-sufficiency contract as
2912 // compute_dcpdT ( Grok phase-3 hardening, 2026-08-13 )
2913 this->compute_T( aIndex );
2914
2915 // consistent tangent: a clamped rho has zero derivative, and
2916 // a clamped T iterate zeroes ALL dT derivatives — same
2917 // contract as compute_dcpdT / compute_dlambdadT ( Newton
2918 // must not push on a flat clamp; 2026-08-13 audit )
2919 mdRhodT = ( mTClamped || mRhoClamped ) ? 0.0 : ( this->*mFundRhodT )( aIndex );
2920 this->set( MaxwellDataValue::drhodT, aIndex );
2921 }
2922 return mdRhodT ;
2923 }
2924
2925
2926 inline real
2927 calculator::MaxwellData::compute_T( const uint aIndex )
2928 {
2929 return ( this->*mFunT )( aIndex );
2930 }
2931
2932 inline real
2933 calculator::MaxwellData::compute_cp( const uint aIndex )
2934 {
2935 if ( ! this->is_current( MaxwellDataValue::cp, aIndex ) )
2936 {
2937 mCp = mMaterial->cp( this->compute_T( aIndex ) );
2938 this->set( MaxwellDataValue::cp, aIndex );
2939 }
2940 return mCp ;
2941 }
2942
2943 inline real
2944 calculator::MaxwellData::compute_dcpdT( const uint aIndex )
2945 {
2946 if ( ! this->is_current( MaxwellDataValue::dcpdT, aIndex ) )
2947 {
2948 real T = this->compute_T( aIndex );
2949
2950 // consistent tangent: no T-sensitivity while T is clamped
2951 mdCpdT = mTClamped ? 0.0 : mMaterial->dcpdT( T );
2952 this->set( MaxwellDataValue::dcpdT, aIndex );
2953 }
2954 return mdCpdT ;
2955 }
2956
2957 inline real
2958 calculator::MaxwellData::compute_lambda( const uint aIndex )
2959 {
2960 if ( ! this->is_current( MaxwellDataValue::lambda, aIndex ) )
2961 {
2962 mLambda = ( this->*mFunLambda)( aIndex );
2963 this->set( MaxwellDataValue::lambda, aIndex );
2964 }
2965 return mLambda ;
2966 }
2967
2968 inline real
2969 calculator::MaxwellData::compute_dlambdadT( const uint aIndex )
2970 {
2971 if ( ! this->is_current( MaxwellDataValue::dlambdadT, aIndex ) )
2972 {
2973 // the variant computes T internally and updates mTClamped
2974 mdLambdadT = ( this->*mFundLambdadT)( aIndex );
2975
2976 // consistent tangent: no T-sensitivity while T is clamped
2977 if ( mTClamped )
2978 {
2979 mdLambdadT = 0.0 ;
2980 }
2981 this->set( MaxwellDataValue::dlambdadT, aIndex );
2982 }
2983 return mdLambdadT ;
2984 }
2985
2986 inline real
2987 calculator::MaxwellData::compute_T_fem( const uint aIndex )
2988 {
2989 if ( ! this->is_current( MaxwellDataValue::T, aIndex ) )
2990 {
2991 real T = dot( mThermalCalculator->Nvec( aIndex ) , mThermalCalculator->q() ) ;
2992
2993 // transient nonlinear iterates can swing outside the physical
2994 // table window; clamp on both edges and remember: while
2995 // clamped, all dT-derivatives are zero ( consistent tangent ).
2996 // A CONVERGED solution at a clamp is a modeling error — the
2997 // controller can query T_clamped() for a diagnostic.
2998 mTClamped = ( T < gTmin ) || ( T > mTmax ) ;
2999 mT = mTClamped ? std::clamp( T, gTmin, mTmax ) : T ;
3000
3001 this->set( MaxwellDataValue::T, aIndex );
3002 }
3003 return mT ;
3004 }
3005
3006 inline real
3007 calculator::MaxwellData::compute_T_const( const uint aIndex )
3008 {
3009 // same contract as compute_T_fem: a user-set bulk temperature
3010 // outside the material's table window is pinned, and the flag
3011 // zeroes the dT-derivatives ( consistent tangent )
3012 mTClamped = ( gTbulk < gTmin ) || ( gTbulk > mTmax );
3013 mT = mTClamped ? std::clamp( gTbulk, gTmin, mTmax ) : gTbulk ;
3014 return mT ;
3015 }
3016
3017 inline const Vector< real > &
3018 calculator::MaxwellData::compute_h_bulk_node( const uint aIndex )
3019 {
3020 mH = -1. * mMaxwellCalculator->B( aIndex ) * mMaxwellCalculator->q();
3021 return mH ;
3022 }
3023
3024 inline const Vector< real > &
3025 calculator::MaxwellData::compute_h_bulk_edge( const uint aIndex )
3026 {
3027 mH = mMaxwellCalculator->E( aIndex ) * mMaxwellCalculator->q();
3028 return mH ;
3029 }
3030
3031 inline const Vector< real > &
3032 calculator::MaxwellData::compute_h_ts_edge( const uint aIndex )
3033 {
3034 mHn = compute_hn( mMaxwellCalculator, aIndex );
3035 this->set( MaxwellDataValue::n, aIndex ); // compute_hn writes n
3036 mHt = mMaxwellCalculator->E( aIndex ) * mMaxwellCalculator->q();
3037
3038 mH = mHn+mHt ;
3039 return mH ;
3040 }
3041
3042 inline real
3043 calculator::MaxwellData::compute_lambda_bulk( const uint aIndex )
3044 {
3045 return mMaterial->lambda( this->compute_T( aIndex ) );
3046 }
3047
3048 inline real
3049 calculator::MaxwellData::compute_dlambdadT_bulk( const uint aIndex )
3050 {
3051 return mMaterial->dlambdadT( this->compute_T( aIndex ) );
3052 }
3053
3054 inline real
3055 calculator::MaxwellData::compute_lambda_metal( const uint aIndex )
3056 {
3057 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3058 {
3059 const Vector< real > & b = this->compute_b( aIndex );
3060 const Vector< real > & j = this->compute_j( aIndex );
3061
3062 // the bj_angle function updates NormB and normJ
3063 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3064
3065 this->set( MaxwellDataValue::normB, aIndex );
3066 this->set( MaxwellDataValue::normJ, aIndex );
3067 this->set( MaxwellDataValue::beta, aIndex );
3068 }
3069
3070 return mMaterial->lambda( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3071 }
3072
3073 inline real
3074 calculator::MaxwellData::compute_dlambdadT_metal( const uint aIndex )
3075 {
3076 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3077 {
3078 const Vector< real > & b = this->compute_b( aIndex );
3079 const Vector< real > & j = this->compute_j( aIndex );
3080
3081 // the bj_angle function updates NormB and normJ
3082 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3083
3084 this->set( MaxwellDataValue::normB, aIndex );
3085 this->set( MaxwellDataValue::normJ, aIndex );
3086 this->set( MaxwellDataValue::beta, aIndex );
3087 }
3088
3089 return mMaterial->dlambdadT( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3090 }
3091
3092
3093
3094 inline real
3095 calculator::MaxwellData::compute_rho_bulk( const uint aIndex )
3096 {
3097 return mMaterial->rho( this->compute_T( aIndex ) );
3098 }
3099
3100 inline real
3101 calculator::MaxwellData::compute_drhodT_bulk( const uint aIndex )
3102 {
3103 return mMaterial->drhodT( this->compute_T( aIndex ) );
3104 }
3105
3106 inline real
3107 calculator::MaxwellData::norm_b( const uint aIndex )
3108 {
3109 if ( ! this->is_current( MaxwellDataValue::normB, aIndex ) )
3110 {
3111 mNormB = norm( this->compute_b( aIndex ) );
3112 this->set( MaxwellDataValue::normB, aIndex );
3113 }
3114 return mNormB ;
3115 }
3116
3117 inline real
3118 calculator::MaxwellData::norm_j( const uint aIndex )
3119 {
3120 if ( ! this->is_current( MaxwellDataValue::normJ, aIndex ) )
3121 {
3122 mNormJ = norm( this->compute_j( aIndex ) );
3123 this->set( MaxwellDataValue::normJ, aIndex );
3124 }
3125 return mNormJ ;
3126 }
3127
3128 inline real
3129 calculator::MaxwellData::compute_rho_metal( const uint aIndex )
3130 {
3131 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3132 {
3133 const Vector< real > & b = this->compute_b( aIndex );
3134 const Vector< real > & j = this->compute_j( aIndex );
3135
3136 // the bj_angle function updates NormB and normJ
3137 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3138
3139 this->set( MaxwellDataValue::normB, aIndex );
3140 this->set( MaxwellDataValue::normJ, aIndex );
3141 this->set( MaxwellDataValue::beta, aIndex );
3142 }
3143
3144 return mMaterial->rho( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3145 }
3146
3147 // HTS drho/dT family ( T-leg ): preambles mirror the
3148 // corresponding compute_drhodb_* wrappers — bn_angle block for the
3149 // thin-shell variants, beta_dummy for bulk, compute_x for defect
3150
3151 inline real
3152 calculator::MaxwellData::compute_drhodT_powerlaw_ts( const uint aIndex )
3153 {
3154 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3155 {
3156 const Vector< real > & b = this->compute_b( aIndex );
3157
3158 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3159
3160 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3161
3162 this->set( MaxwellDataValue::beta, aIndex );
3163 this->set( MaxwellDataValue::normB, aIndex );
3164 }
3165
3166 return mMaterial->drho_powerlaw_dT(
3167 this->norm_j( aIndex ),
3168 this->compute_T( aIndex ),
3169 this->norm_b( aIndex ),
3170 mBeta );
3171 }
3172
3173 inline real
3174 calculator::MaxwellData::compute_drhodT_powerlaw_ts_defect( const uint aIndex )
3175 {
3176 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3177 {
3178 const Vector< real > & b = this->compute_b( aIndex );
3179
3180 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3181
3182 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3183
3184 this->set( MaxwellDataValue::beta, aIndex );
3185 this->set( MaxwellDataValue::normB, aIndex );
3186 }
3187
3188 this->compute_x( aIndex );
3189
3190 return mMaterial->drho_powerlaw_dT(
3191 this->norm_j( aIndex ),
3192 this->compute_T( aIndex ),
3193 this->norm_b( aIndex ),
3194 mBeta,
3195 mX, mY, mZ, mTime );
3196 }
3197
3198 inline real
3199 calculator::MaxwellData::compute_drhodT_piecewise_ts( const uint aIndex )
3200 {
3201 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3202 {
3203 const Vector< real > & b = this->compute_b( aIndex );
3204
3205 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3206
3207 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3208
3209 this->set( MaxwellDataValue::beta, aIndex );
3210 this->set( MaxwellDataValue::normB, aIndex );
3211 }
3212
3213 return mMaterial->drho_piecewise_dT(
3214 this->norm_j( aIndex ),
3215 this->compute_T( aIndex ),
3216 this->norm_b( aIndex ),
3217 mBeta );
3218 }
3219
3220 inline real
3221 calculator::MaxwellData::compute_drhodT_piecewise_ts_defect( const uint aIndex )
3222 {
3223 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3224 {
3225 const Vector< real > & b = this->compute_b( aIndex );
3226
3227 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3228
3229 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3230
3231 this->set( MaxwellDataValue::beta, aIndex );
3232 this->set( MaxwellDataValue::normB, aIndex );
3233 }
3234
3235 this->compute_x( aIndex );
3236
3237 return mMaterial->drho_piecewise_dT(
3238 this->norm_j( aIndex ),
3239 this->compute_T( aIndex ),
3240 this->norm_b( aIndex ),
3241 mBeta,
3242 mX, mY, mZ, mTime );
3243 }
3244
3245 inline real
3246 calculator::MaxwellData::compute_drhodT_powerlaw_bulk( const uint aIndex )
3247 {
3248 return mMaterial->drho_powerlaw_dT(
3249 this->norm_j( aIndex ),
3250 this->compute_T( aIndex ),
3251 this->norm_b( aIndex ),
3252 this->beta_dummy() );
3253 }
3254
3255 inline real
3256 calculator::MaxwellData::compute_drhodT_powerlaw_bulk_defect( const uint aIndex )
3257 {
3258 this->compute_x( aIndex );
3259
3260 return mMaterial->drho_powerlaw_dT(
3261 this->norm_j( aIndex ),
3262 this->compute_T( aIndex ),
3263 this->norm_b( aIndex ),
3264 this->beta_dummy(),
3265 mX, mY, mZ, mTime );
3266 }
3267
3268 inline real
3269 calculator::MaxwellData::compute_drhodT_piecewise_bulk( const uint aIndex )
3270 {
3271 return mMaterial->drho_piecewise_dT(
3272 this->norm_j( aIndex ),
3273 this->compute_T( aIndex ),
3274 this->norm_b( aIndex ),
3275 this->beta_dummy() );
3276 }
3277
3278 inline real
3279 calculator::MaxwellData::compute_drhodT_piecewise_bulk_defect( const uint aIndex )
3280 {
3281 this->compute_x( aIndex );
3282
3283 return mMaterial->drho_piecewise_dT(
3284 this->norm_j( aIndex ),
3285 this->compute_T( aIndex ),
3286 this->norm_b( aIndex ),
3287 this->beta_dummy(),
3288 mX, mY, mZ, mTime );
3289 }
3290
3291//------------------------------------------------------------------------------
3292// riva-law wrappers ( 2026-08-27 ), mirroring the piecewise family
3293//------------------------------------------------------------------------------
3294
3295 inline real
3296 calculator::MaxwellData::compute_rho_riva_ts( const uint aIndex )
3297 {
3298 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3299 {
3300 const Vector< real > & b = this->compute_b( aIndex );
3301
3302 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3303
3304 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3305
3306 this->set( MaxwellDataValue::beta, aIndex );
3307 this->set( MaxwellDataValue::normB, aIndex );
3308 }
3309
3310 return mMaterial->rho_riva(
3311 this->norm_j( aIndex ),
3312 this->compute_T( aIndex ),
3313 this->norm_b( aIndex ),
3314 mBeta );
3315 }
3316
3317 inline real
3318 calculator::MaxwellData::compute_rho_riva_ts_defect( const uint aIndex )
3319 {
3320 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3321 {
3322 const Vector< real > & b = this->compute_b( aIndex );
3323
3324 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3325
3326 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3327
3328 this->set( MaxwellDataValue::beta, aIndex );
3329 this->set( MaxwellDataValue::normB, aIndex );
3330 }
3331
3332 this->compute_x( aIndex );
3333
3334 return mMaterial->rho_riva(
3335 this->norm_j( aIndex ),
3336 this->compute_T( aIndex ),
3337 this->norm_b( aIndex ),
3338 mBeta,
3339 mX, mY, mZ, mTime );
3340 }
3341
3342 inline real
3343 calculator::MaxwellData::compute_rho_riva_bulk( const uint aIndex )
3344 {
3345 return mMaterial->rho_riva(
3346 this->norm_j( aIndex ),
3347 this->compute_T( aIndex ),
3348 this->norm_b( aIndex ),
3349 this->beta_dummy() );
3350 }
3351
3352 inline real
3353 calculator::MaxwellData::compute_rho_riva_bulk_defect( const uint aIndex )
3354 {
3355 this->compute_x( aIndex );
3356
3357 return mMaterial->rho_riva(
3358 this->norm_j( aIndex ),
3359 this->compute_T( aIndex ),
3360 this->norm_b( aIndex ),
3361 this->beta_dummy(),
3362 mX, mY, mZ, mTime );
3363 }
3364
3365 inline real
3366 calculator::MaxwellData::compute_drhodj_riva_ts( const uint aIndex )
3367 {
3368 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3369 {
3370 const Vector< real > & b = this->compute_b( aIndex );
3371
3372 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3373
3374 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3375
3376 this->set( MaxwellDataValue::beta, aIndex );
3377 this->set( MaxwellDataValue::normB, aIndex );
3378 }
3379
3380 return mMaterial->drho_riva_dJ(
3381 this->norm_j( aIndex ),
3382 this->compute_T( aIndex ),
3383 this->norm_b( aIndex ),
3384 mBeta );
3385 }
3386
3387 inline real
3388 calculator::MaxwellData::compute_drhodj_riva_ts_defect( const uint aIndex )
3389 {
3390 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3391 {
3392 const Vector< real > & b = this->compute_b( aIndex );
3393
3394 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3395
3396 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3397
3398 this->set( MaxwellDataValue::beta, aIndex );
3399 this->set( MaxwellDataValue::normB, aIndex );
3400 }
3401
3402 this->compute_x( aIndex );
3403
3404 return mMaterial->drho_riva_dJ(
3405 this->norm_j( aIndex ),
3406 this->compute_T( aIndex ),
3407 this->norm_b( aIndex ),
3408 mBeta,
3409 mX, mY, mZ, mTime );
3410 }
3411
3412 inline real
3413 calculator::MaxwellData::compute_drhodj_riva_bulk( const uint aIndex )
3414 {
3415 return mMaterial->drho_riva_dJ(
3416 this->norm_j( aIndex ),
3417 this->compute_T( aIndex ),
3418 this->norm_b( aIndex ),
3419 this->beta_dummy() );
3420 }
3421
3422 inline real
3423 calculator::MaxwellData::compute_drhodj_riva_bulk_defect( const uint aIndex )
3424 {
3425 this->compute_x( aIndex );
3426
3427 return mMaterial->drho_riva_dJ(
3428 this->norm_j( aIndex ),
3429 this->compute_T( aIndex ),
3430 this->norm_b( aIndex ),
3431 this->beta_dummy(),
3432 mX, mY, mZ, mTime );
3433 }
3434
3435 inline real
3436 calculator::MaxwellData::compute_drhodb_riva_ts( const uint aIndex )
3437 {
3438 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3439 {
3440 const Vector< real > & b = this->compute_b( aIndex );
3441
3442 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3443
3444 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3445
3446 this->set( MaxwellDataValue::beta, aIndex );
3447 this->set( MaxwellDataValue::normB, aIndex );
3448 }
3449
3450 return mMaterial->drho_riva_dB(
3451 this->norm_j( aIndex ),
3452 this->compute_T( aIndex ),
3453 this->norm_b( aIndex ),
3454 mBeta );
3455 }
3456
3457 inline real
3458 calculator::MaxwellData::compute_drhodb_riva_ts_defect( const uint aIndex )
3459 {
3460 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3461 {
3462 const Vector< real > & b = this->compute_b( aIndex );
3463
3464 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3465
3466 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3467
3468 this->set( MaxwellDataValue::beta, aIndex );
3469 this->set( MaxwellDataValue::normB, aIndex );
3470 }
3471
3472 this->compute_x( aIndex );
3473
3474 return mMaterial->drho_riva_dB(
3475 this->norm_j( aIndex ),
3476 this->compute_T( aIndex ),
3477 this->norm_b( aIndex ),
3478 mBeta,
3479 mX, mY, mZ, mTime );
3480 }
3481
3482 inline real
3483 calculator::MaxwellData::compute_drhodb_riva_bulk( const uint aIndex )
3484 {
3485 return mMaterial->drho_riva_dB(
3486 this->norm_j( aIndex ),
3487 this->compute_T( aIndex ),
3488 this->norm_b( aIndex ),
3489 this->beta_dummy() );
3490 }
3491
3492 inline real
3493 calculator::MaxwellData::compute_drhodb_riva_bulk_defect( const uint aIndex )
3494 {
3495 this->compute_x( aIndex );
3496
3497 return mMaterial->drho_riva_dB(
3498 this->norm_j( aIndex ),
3499 this->compute_T( aIndex ),
3500 this->norm_b( aIndex ),
3501 this->beta_dummy(),
3502 mX, mY, mZ, mTime );
3503 }
3504
3505 inline real
3506 calculator::MaxwellData::compute_drhodT_riva_ts( const uint aIndex )
3507 {
3508 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3509 {
3510 const Vector< real > & b = this->compute_b( aIndex );
3511
3512 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3513
3514 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3515
3516 this->set( MaxwellDataValue::beta, aIndex );
3517 this->set( MaxwellDataValue::normB, aIndex );
3518 }
3519
3520 return mMaterial->drho_riva_dT(
3521 this->norm_j( aIndex ),
3522 this->compute_T( aIndex ),
3523 this->norm_b( aIndex ),
3524 mBeta );
3525 }
3526
3527 inline real
3528 calculator::MaxwellData::compute_drhodT_riva_ts_defect( const uint aIndex )
3529 {
3530 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3531 {
3532 const Vector< real > & b = this->compute_b( aIndex );
3533
3534 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3535
3536 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3537
3538 this->set( MaxwellDataValue::beta, aIndex );
3539 this->set( MaxwellDataValue::normB, aIndex );
3540 }
3541
3542 this->compute_x( aIndex );
3543
3544 return mMaterial->drho_riva_dT(
3545 this->norm_j( aIndex ),
3546 this->compute_T( aIndex ),
3547 this->norm_b( aIndex ),
3548 mBeta,
3549 mX, mY, mZ, mTime );
3550 }
3551
3552 inline real
3553 calculator::MaxwellData::compute_drhodT_riva_bulk( const uint aIndex )
3554 {
3555 return mMaterial->drho_riva_dT(
3556 this->norm_j( aIndex ),
3557 this->compute_T( aIndex ),
3558 this->norm_b( aIndex ),
3559 this->beta_dummy() );
3560 }
3561
3562 inline real
3563 calculator::MaxwellData::compute_drhodT_riva_bulk_defect( const uint aIndex )
3564 {
3565 this->compute_x( aIndex );
3566
3567 return mMaterial->drho_riva_dT(
3568 this->norm_j( aIndex ),
3569 this->compute_T( aIndex ),
3570 this->norm_b( aIndex ),
3571 this->beta_dummy(),
3572 mX, mY, mZ, mTime );
3573 }
3574
3575 inline real
3576 calculator::MaxwellData::compute_drhodT_metal( const uint aIndex )
3577 {
3578 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3579 {
3580 const Vector< real > & b = this->compute_b( aIndex );
3581 const Vector< real > & j = this->compute_j( aIndex );
3582
3583 // the bj_angle function updates NormB and normJ
3584 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3585
3586 this->set( MaxwellDataValue::normB, aIndex );
3587 this->set( MaxwellDataValue::normJ, aIndex );
3588 this->set( MaxwellDataValue::beta, aIndex );
3589 }
3590
3591 return mMaterial->drhodT( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3592 }
3593
3594 inline real
3595 calculator::MaxwellData::compute_drhodb( const uint aIndex )
3596 {
3597 if ( ! this->is_current( MaxwellDataValue::drhodb, aIndex ) )
3598 {
3599 // make sure rho — and with it mRhoClamped — is current here
3600 this->compute_rho( aIndex );
3601
3602 // consistent tangent: a clamped rho has zero derivative
3603 mdRhodB = mRhoClamped ? 0.0 : ( this->*mFundRhodB )( aIndex );
3604 this->set( MaxwellDataValue::drhodb, aIndex );
3605 }
3606 return mdRhodB ;
3607 }
3608
3609 inline real
3610 calculator::MaxwellData::compute_drhodbeta( const uint aIndex )
3611 {
3612 if ( ! this->is_current( MaxwellDataValue::drhodbeta, aIndex ) )
3613 {
3614 // make sure rho — and with it mRhoClamped — is current here
3615 this->compute_rho( aIndex );
3616
3617 // consistent tangent: a clamped rho has zero derivative
3618 mdRhodBeta = mRhoClamped ? 0.0 : ( this->*mFundRhodBeta )( aIndex );
3619 this->set( MaxwellDataValue::drhodbeta, aIndex );
3620 }
3621 return mdRhodBeta ;
3622 }
3623
3624 inline real
3625 calculator::MaxwellData::compute_drhodb_metal( const uint aIndex )
3626 {
3627 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3628 {
3629 const Vector< real > & b = this->compute_b( aIndex );
3630 const Vector< real > & j = this->compute_j( aIndex );
3631
3632 // the bj_angle function updates NormB and normJ
3633 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3634
3635 this->set( MaxwellDataValue::normB, aIndex );
3636 this->set( MaxwellDataValue::normJ, aIndex );
3637 this->set( MaxwellDataValue::beta, aIndex );
3638 }
3639
3640 return mMaterial->drhodB( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3641 }
3642
3643 inline real
3644 calculator::MaxwellData::compute_drhodbeta_metal( const uint aIndex )
3645 {
3646 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3647 {
3648 const Vector< real > & b = this->compute_b( aIndex );
3649 const Vector< real > & j = this->compute_j( aIndex );
3650
3651 // the bj_angle function updates NormB and normJ
3652 mBeta = mMaxwellCalculator->bj_angle( b, j, mNormB, mNormJ );
3653
3654 this->set( MaxwellDataValue::normB, aIndex );
3655 this->set( MaxwellDataValue::normJ, aIndex );
3656 this->set( MaxwellDataValue::beta, aIndex );
3657 }
3658
3659 return mMaterial->drhodbeta( this->compute_T( aIndex ), this->norm_b( aIndex ), mBeta );
3660 }
3661
3662//------------------------------------------------------------------------------
3663// HTS drho/d|B| wrappers: each mirrors the preamble of its
3664// compute_rho_* counterpart exactly, so the cached beta / normB state is
3665// updated identically whichever of the pair runs first in an iterate.
3666//------------------------------------------------------------------------------
3667
3668 inline real
3669 calculator::MaxwellData::compute_drhodb_powerlaw_ts( const uint aIndex )
3670 {
3671 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3672 {
3673 const Vector< real > & b = this->compute_b( aIndex );
3674
3675 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3676
3677 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3678
3679 this->set( MaxwellDataValue::beta, aIndex );
3680 this->set( MaxwellDataValue::normB, aIndex );
3681 }
3682
3683 return mMaterial->drho_powerlaw_dB(
3684 this->norm_j( aIndex ),
3685 this->compute_T( aIndex ),
3686 this->norm_b( aIndex ),
3687 mBeta );
3688 }
3689
3690 inline real
3691 calculator::MaxwellData::compute_drhodb_powerlaw_ts_defect( const uint aIndex )
3692 {
3693 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3694 {
3695 const Vector< real > & b = this->compute_b( aIndex );
3696
3697 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3698
3699 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3700
3701 this->set( MaxwellDataValue::beta, aIndex );
3702 this->set( MaxwellDataValue::normB, aIndex );
3703 }
3704
3705 this->compute_x( aIndex );
3706
3707 return mMaterial->drho_powerlaw_dB(
3708 this->norm_j( aIndex ),
3709 this->compute_T( aIndex ),
3710 this->norm_b( aIndex ),
3711 mBeta,
3712 mX, mY, mZ, mTime );
3713 }
3714
3715 inline real
3716 calculator::MaxwellData::compute_drhodb_piecewise_ts( const uint aIndex )
3717 {
3718 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3719 {
3720 const Vector< real > & b = this->compute_b( aIndex );
3721
3722 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3723
3724 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3725
3726 this->set( MaxwellDataValue::beta, aIndex );
3727 this->set( MaxwellDataValue::normB, aIndex );
3728 }
3729
3730 return mMaterial->drho_piecewise_dB(
3731 this->norm_j( aIndex ),
3732 this->compute_T( aIndex ),
3733 this->norm_b( aIndex ),
3734 mBeta );
3735 }
3736
3737 inline real
3738 calculator::MaxwellData::compute_drhodb_piecewise_ts_defect( const uint aIndex )
3739 {
3740 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3741 {
3742 const Vector< real > & b = this->compute_b( aIndex );
3743
3744 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3745
3746 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3747
3748 this->set( MaxwellDataValue::beta, aIndex );
3749 this->set( MaxwellDataValue::normB, aIndex );
3750 }
3751
3752 this->compute_x( aIndex );
3753
3754 return mMaterial->drho_piecewise_dB(
3755 this->norm_j( aIndex ),
3756 this->compute_T( aIndex ),
3757 this->norm_b( aIndex ),
3758 mBeta,
3759 mX, mY, mZ, mTime );
3760 }
3761
3762 inline real
3763 calculator::MaxwellData::compute_drhodb_powerlaw_bulk( const uint aIndex )
3764 {
3765 return mMaterial->drho_powerlaw_dB(
3766 this->norm_j( aIndex ),
3767 this->compute_T( aIndex ),
3768 this->norm_b( aIndex ),
3769 this->beta_dummy() );
3770 }
3771
3772 inline real
3773 calculator::MaxwellData::compute_drhodb_powerlaw_bulk_defect( const uint aIndex )
3774 {
3775 this->compute_x( aIndex );
3776
3777 return mMaterial->drho_powerlaw_dB(
3778 this->norm_j( aIndex ),
3779 this->compute_T( aIndex ),
3780 this->norm_b( aIndex ),
3781 this->beta_dummy(),
3782 mX, mY, mZ, mTime );
3783 }
3784
3785 inline real
3786 calculator::MaxwellData::compute_drhodb_piecewise_bulk( const uint aIndex )
3787 {
3788 return mMaterial->drho_piecewise_dB(
3789 this->norm_j( aIndex ),
3790 this->compute_T( aIndex ),
3791 this->norm_b( aIndex ),
3792 this->beta_dummy() );
3793 }
3794
3795 inline real
3796 calculator::MaxwellData::compute_drhodb_piecewise_bulk_defect( const uint aIndex )
3797 {
3798 this->compute_x( aIndex );
3799
3800 return mMaterial->drho_piecewise_dB(
3801 this->norm_j( aIndex ),
3802 this->compute_T( aIndex ),
3803 this->norm_b( aIndex ),
3804 this->beta_dummy(),
3805 mX, mY, mZ, mTime );
3806 }
3807
3808 inline real
3809 calculator::MaxwellData::compute_rho_powerlaw_bulk( const uint aIndex )
3810 {
3811 return mMaterial->rho_powerlaw(
3812 this->norm_j( aIndex ),
3813 this->compute_T( aIndex ),
3814 this->norm_b( aIndex ),
3815 this->beta_dummy() );
3816 }
3817
3818 inline real
3819 calculator::MaxwellData::compute_rho_piecewise_bulk( const uint aIndex )
3820 {
3821 return mMaterial->rho_piecewise(
3822 this->norm_j( aIndex ),
3823 this->compute_T( aIndex ),
3824 this->norm_b( aIndex ),
3825 this->beta_dummy() );
3826 }
3827
3828 inline real
3829 calculator::MaxwellData::compute_drhodj_powerlaw_bulk( const uint aIndex )
3830 {
3831 return mMaterial->drho_powerlaw_dJ(
3832 this->norm_j( aIndex ),
3833 this->compute_T( aIndex ),
3834 this->norm_b( aIndex ),
3835 this->beta_dummy() );
3836 }
3837
3838 inline real
3839 calculator::MaxwellData::compute_drhodj_piecewise_bulk( const uint aIndex )
3840 {
3841 return mMaterial->drho_piecewise_dJ(
3842 this->norm_j( aIndex ),
3843 this->compute_T( aIndex ),
3844 this->norm_b( aIndex ),
3845 this->beta_dummy() );
3846 }
3847
3848 inline real
3849 calculator::MaxwellData::compute_rho_powerlaw_ts( const uint aIndex )
3850 {
3851
3852 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3853 {
3854 const Vector< real > & b = this->compute_b( aIndex );
3855
3856 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3857
3858 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3859
3860 this->set( MaxwellDataValue::beta, aIndex );
3861 this->set( MaxwellDataValue::normB, aIndex );
3862 }
3863
3864 return mMaterial->rho_powerlaw(
3865 this->norm_j( aIndex ),
3866 this->compute_T( aIndex ),
3867 this->norm_b( aIndex ),
3868 mBeta );
3869 }
3870
3871 inline real
3872 calculator::MaxwellData::compute_drhodj_powerlaw_ts( const uint aIndex )
3873 {
3874
3875 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3876 {
3877 const Vector< real > & b = this->compute_b( aIndex );
3878
3879 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3880
3881 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3882
3883 this->set( MaxwellDataValue::beta, aIndex );
3884 this->set( MaxwellDataValue::normB, aIndex );
3885 }
3886
3887 return mMaterial->drho_powerlaw_dJ(
3888 this->norm_j( aIndex ),
3889 this->compute_T( aIndex ),
3890 this->norm_b( aIndex ),
3891 mBeta );
3892 }
3893
3894 inline real
3895 calculator::MaxwellData::compute_rho_piecewise_ts( const uint aIndex )
3896 {
3897
3898 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3899 {
3900 const Vector< real > & b = this->compute_b( aIndex );
3901
3902 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3903
3904 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3905
3906 this->set( MaxwellDataValue::beta, aIndex );
3907 this->set( MaxwellDataValue::normB, aIndex );
3908 }
3909
3910 return mMaterial->rho_piecewise(
3911 this->norm_j( aIndex ),
3912 this->compute_T( aIndex ),
3913 this->norm_b( aIndex ),
3914 mBeta );
3915 }
3916
3917 inline real
3918 calculator::MaxwellData::compute_drhodj_piecewise_ts( const uint aIndex )
3919 {
3920
3921 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3922 {
3923 const Vector< real > & b = this->compute_b( aIndex );
3924
3925 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3926
3927 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3928
3929 this->set( MaxwellDataValue::beta, aIndex );
3930 this->set( MaxwellDataValue::normB, aIndex );
3931 }
3932
3933 return mMaterial->drho_piecewise_dJ(
3934 this->norm_j( aIndex ),
3935 this->compute_T( aIndex ),
3936 this->norm_b( aIndex ),
3937 mBeta );
3938 }
3939
3940 inline real
3941 calculator::MaxwellData::compute_rho_powerlaw_bulk_defect( const uint aIndex )
3942 {
3943 this->compute_x( aIndex );
3944
3945 return mMaterial->rho_powerlaw(
3946 this->norm_j( aIndex ),
3947 this->compute_T( aIndex ),
3948 this->norm_b( aIndex ),
3949 this->beta_dummy(),
3950 mX, mY, mZ, mTime );
3951 }
3952
3953 inline real
3954 calculator::MaxwellData::compute_drhodj_powerlaw_bulk_defect( const uint aIndex )
3955 {
3956 this->compute_x( aIndex );
3957
3958 return mMaterial->drho_powerlaw_dJ(
3959 this->norm_j( aIndex ),
3960 this->compute_T( aIndex ),
3961 this->norm_b( aIndex ),
3962 this->beta_dummy(),
3963 mX, mY, mZ, mTime );
3964 }
3965
3966 inline real
3967 calculator::MaxwellData::compute_rho_powerlaw_ts_defect( const uint aIndex )
3968 {
3969
3970 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3971 {
3972 const Vector< real > & b = this->compute_b( aIndex );
3973
3974 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
3975
3976 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
3977
3978 this->set( MaxwellDataValue::beta, aIndex );
3979 this->set( MaxwellDataValue::normB, aIndex );
3980 }
3981
3982 this->compute_x( aIndex );
3983
3984 return mMaterial->rho_powerlaw(
3985 this->norm_j( aIndex ),
3986 this->compute_T( aIndex ),
3987 this->norm_b( aIndex ),
3988 mBeta,
3989 mX, mY, mZ, mTime );
3990 }
3991
3992 inline real
3993 calculator::MaxwellData::compute_drhodj_powerlaw_ts_defect( const uint aIndex )
3994 {
3995
3996 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
3997 {
3998 const Vector< real > & b = this->compute_b( aIndex );
3999
4000 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
4001
4002 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
4003
4004 this->set( MaxwellDataValue::beta, aIndex );
4005 this->set( MaxwellDataValue::normB, aIndex );
4006 }
4007
4008 this->compute_x( aIndex );
4009
4010 return mMaterial->drho_powerlaw_dJ(
4011 this->norm_j( aIndex ),
4012 this->compute_T( aIndex ),
4013 this->norm_b( aIndex ),
4014 mBeta,
4015 mX, mY, mZ, mTime );
4016 }
4017
4018 inline real
4019 calculator::MaxwellData::compute_rho_piecewise_bulk_defect( const uint aIndex )
4020 {
4021 this->compute_x( aIndex );
4022 return mMaterial->rho_piecewise(
4023 this->norm_j( aIndex ),
4024 this->compute_T( aIndex ),
4025 this->norm_b( aIndex ),
4026 this->beta_dummy(),
4027 mX, mY, mZ, mTime );
4028 }
4029
4030 inline real
4031 calculator::MaxwellData::compute_drhodj_piecewise_bulk_defect( const uint aIndex )
4032 {
4033 this->compute_x( aIndex );
4034 return mMaterial->drho_piecewise_dJ(
4035 this->norm_j( aIndex ),
4036 this->compute_T( aIndex ),
4037 this->norm_b( aIndex ),
4038 this->beta_dummy(),
4039 mX, mY, mZ, mTime );
4040 }
4041
4042 inline real
4043 calculator::MaxwellData::compute_rho_piecewise_ts_defect( const uint aIndex )
4044 {
4045
4046 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
4047 {
4048 const Vector< real > & b = this->compute_b( aIndex );
4049
4050 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
4051
4052 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
4053
4054 this->set( MaxwellDataValue::beta, aIndex );
4055 this->set( MaxwellDataValue::normB, aIndex );
4056 }
4057
4058 this->compute_x( aIndex );
4059
4060 return mMaterial->rho_piecewise(
4061 this->norm_j( aIndex ),
4062 this->compute_T( aIndex ),
4063 this->norm_b( aIndex ),
4064 mBeta,
4065 mX, mY, mZ, mTime);
4066 }
4067
4068
4069 inline real
4070 calculator::MaxwellData::compute_drhodj_piecewise_ts_defect( const uint aIndex )
4071 {
4072
4073 if ( ! this->is_current( MaxwellDataValue::beta, aIndex ) )
4074 {
4075 const Vector< real > & b = this->compute_b( aIndex );
4076
4077 BELFEM_ASSERT( this->is_current( MaxwellDataValue::n, aIndex ), "normal vector has not been updated" );
4078
4079 mBeta = mMaxwellCalculator->bn_angle( b, mN, mNormB );
4080
4081 this->set( MaxwellDataValue::beta, aIndex );
4082 this->set( MaxwellDataValue::normB, aIndex );
4083 }
4084
4085 this->compute_x( aIndex );
4086
4087 return mMaterial->drho_piecewise_dJ(
4088 this->norm_j( aIndex ),
4089 this->compute_T( aIndex ),
4090 this->norm_b( aIndex ),
4091 mBeta,
4092 mX, mY, mZ, mTime );
4093 }
4094
4095 inline real
4096 calculator::MaxwellData::compute_dmudh( const uint aIndex )
4097 {
4098 if ( ! this->is_current( MaxwellDataValue::dmudh, aIndex ) )
4099 {
4100 this->set( MaxwellDataValue::dmudh, aIndex );
4101 return ( this->*mFundMudH ) ( aIndex );
4102 }
4103 return mdMudH ;
4104 }
4105
4106
4107 inline real
4108 calculator::MaxwellData::compute_dmu_zero( const uint aIndex )
4109 {
4110 mdMudH = 0.0 ;
4111 return mdMudH ;
4112 }
4113
4114 inline real
4115 calculator::MaxwellData::compute_dmu_material( const uint aIndex )
4116 {
4117 if ( ! this->is_current( MaxwellDataValue::normH, aIndex ) )
4118 {
4119 mNormH = norm( this->compute_h( aIndex ) );
4120 this->set( MaxwellDataValue::normH, aIndex );
4121 }
4122 mMaterial->dmudH( mNormH, mMu, mdMudH );
4123 this->set( MaxwellDataValue::mu, aIndex );
4124 return mdMudH ;
4125 }
4126
4127 inline real
4128 calculator::MaxwellData::compute_volumetric_heatload( const uint aIndex )
4129 {
4130 return ( this->*mFunHeat )( aIndex );
4131 }
4132
4133 inline real
4134 calculator::MaxwellData::compute_heatload_user( const uint aIndex )
4135 {
4136 this->compute_x( aIndex );
4137
4138 return mMaterial->volumetric_heatload( mX, mY, mZ, mTime );
4139 }
4140
4141 inline real
4142 calculator::MaxwellData::return_zero( const uint aIndex )
4143 {
4144 return 0.0 ;
4145 }
4146
4147 inline real
4148 calculator::MaxwellData::beta_dummy() const
4149 {
4150 // bulk HTS has no meaningful tape normal, so the material call
4151 // takes a dummy beta that must still be the reset() value.
4152 // beta conventions: ( b, n ) bn_angle for HTS thin shells,
4153 // ( b, j ) bj_angle for metals — never both in one material
4154 // ( enforced in the constructor )
4155 BELFEM_ASSERT( mBeta == 0.5 * constant::pi,
4156 "bulk-HTS dummy beta was polluted ( expected pi/2 )" );
4157 return mBeta ;
4158 }
4159
4160 inline real
4161 calculator::MaxwellData::density() const
4162 {
4163 BELFEM_ASSERT( ! std::isnan( mDensity ),
4164 "material %s provides neither density nor ref_density",
4165 mMaterial->label().c_str() );
4166 return mDensity ;
4167 }
4168
4169 inline bool
4170 calculator::MaxwellData::T_clamped() const
4171 {
4172 return mTClamped ;
4173 }
4174
4175 inline bool
4176 calculator::MaxwellData::rho_clamped() const
4177 {
4178 return mRhoClamped ;
4179 }
4180
4181 inline Calculator *
4182 calculator::MaxwellData::maxwell()
4183 {
4184 return mMaxwellCalculator ;
4185 }
4186
4187 inline Calculator *
4188 calculator::MaxwellData::thermal()
4189 {
4190 return mThermalCalculator ;
4191 }
4192
4193 inline calculator::MaxwellData *
4194 Calculator::maxwell()
4195 {
4196 BELFEM_ASSERT( mMaxwellData != nullptr, "MaxwellData has not been initialized" );
4197 return mMaxwellData ;
4198 }
4199
4200
4201 }
4202}
4203
4204#endif //BELFEM_CL_FEM_CALCULATOR_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_EF_EdgeFunction.hpp
cl_IF_IntegrationData.hpp
cl_Material.hpp
cl_Matrix.hpp
cl_Mesh.hpp
cl_Vector.hpp
belfem::Cell
Cell is a wrapper around the standard vector.
Definition cl_Cell.hpp:42
belfem::Map
Hash map (unordered key-value).
Definition cl_Map.hpp:75
belfem::Material
Base class for all materials in BELFEM.
Definition cl_Material.hpp:279
belfem::Mesh
Top-level container for all mesh entities.
Definition cl_Mesh.hpp:60
belfem::Vector::fill
void fill(const T &aValue)
write value into all entries of the vector
Definition cl_AR_Vector.hpp:226
belfem::fem::Calculator
Definition cl_FEM_Calculator.hpp:600
belfem::fem::Calculator::M
Matrix< real > & M()
return the mass matrix
Definition cl_FEM_Calculator.hpp:2376
belfem::fem::Calculator::master_integration
const IntegrationData * master_integration() const
Definition cl_FEM_Calculator.hpp:2328
belfem::fem::Calculator::nedelec_data_slave_h
const Vector< real > & nedelec_data_slave_h()
slave twin of nedelec_data_master_h(), into vector( "nedelec_h_s" )
Definition cl_FEM_Calculator.cpp:2201
belfem::fem::Calculator::slave_integration
const IntegrationData * slave_integration() const
Definition cl_FEM_Calculator.hpp:2336
belfem::fem::Calculator::Cs
const Matrix< real > & Cs(const uint aIndex)
curl interpolator on slave
Definition cl_FEM_Calculator.cpp:3096
belfem::fem::Calculator::timestep
real timestep() const
Definition cl_FEM_Calculator.hpp:2408
belfem::fem::Calculator::G
const Matrix< real > & G(const uint aIndex)
gradient interpolator for nedelec elements
Definition cl_FEM_Calculator.cpp:3032
belfem::fem::Calculator::link
void link(Group *aGroup)
Definition cl_FEM_Calculator.cpp:1766
belfem::fem::Calculator::Cm
const Matrix< real > & Cm(const uint aIndex)
curl interpolator on master
Definition cl_FEM_Calculator.cpp:3074
belfem::fem::Calculator::volume_is_conductor
bool volume_is_conductor(const mesh::Element *aVolume) const
true if the volume element sits on an h-conductor block ( DomainType::Conductor ),...
Definition cl_FEM_Calculator.cpp:2055
belfem::fem::Calculator::initialize_integration
void initialize_integration(const ElementType aElementType, const InterpolationType aInterpolationType)
Definition cl_FEM_Calculator.cpp:893
belfem::fem::Calculator::print_dofs
void print_dofs()
Definition cl_FEM_Calculator.cpp:2128
belfem::fem::Calculator::integration_order
uint integration_order() const
Definition cl_FEM_Calculator.hpp:2512
belfem::fem::Calculator::dS
real dS(const uint aIndex=0)
surface increment
Definition cl_FEM_Calculator.hpp:2090
belfem::fem::Calculator::Es
const Matrix< real > & Es(const uint aIndex)
edge interpolator on slave
Definition cl_FEM_Calculator.cpp:3085
belfem::fem::Calculator::node_interp
real node_interp(const uint aIndex, const Vector< real > &aNodeValues) const
node interpolation function for local values
Definition cl_FEM_Calculator.hpp:2059
belfem::fem::Calculator::allocate_memory
void allocate_memory()
Definition cl_FEM_Calculator.cpp:950
belfem::fem::Calculator::K
Matrix< real > & K()
return the stiffness matrix
Definition cl_FEM_Calculator.hpp:2368
belfem::fem::Calculator::matrix
Matrix< real > & matrix(const string &aLabel)
return a matrix object
Definition cl_FEM_Calculator.hpp:1680
belfem::fem::Calculator::Nvec
const Vector< real > & Nvec(const uint aIndex) const
node interpolator, but as vector
Definition cl_FEM_Calculator.hpp:2066
belfem::fem::Calculator::Xm
const Matrix< real > & Xm() const
node coordinates on master element
Definition cl_FEM_Calculator.hpp:1696
belfem::fem::Calculator::invJ
const Matrix< real > & invJ(const uint aIndex)
inverse of jacobian matrix
Definition cl_FEM_Calculator.hpp:1758
belfem::fem::Calculator::allocate
void allocate()
called by dof manager
Definition cl_FEM_Calculator.cpp:1204
belfem::fem::Calculator::link_maxwell
void link_maxwell(Kernel *aMaxwellKernel, Kernel *aThermalKernel=nullptr)
register the kernels for the maxwell data helper.
Definition cl_FEM_Calculator.cpp:1733
belfem::fem::Calculator::Nm
const Matrix< real > & Nm(const uint aIndex)
node interpolator master
Definition cl_FEM_Calculator.hpp:2027
belfem::fem::Calculator::Em
const Matrix< real > & Em(const uint aIndex)
edge interpolator on master
Definition cl_FEM_Calculator.cpp:3063
belfem::fem::Calculator::element_is_linear
bool element_is_linear() const
Definition cl_FEM_Calculator.hpp:2526
belfem::fem::Calculator::integration
const IntegrationData * integration() const
Definition cl_FEM_Calculator.hpp:2320
belfem::fem::Calculator::matrix_exists
bool matrix_exists(const string aLabel) const
Definition cl_FEM_Calculator.hpp:1672
belfem::fem::Calculator::material
const Material * material() const
Definition cl_FEM_Calculator.cpp:3107
belfem::fem::Calculator::set_integration_order
void set_integration_order(const uint aOrder)
Definition cl_FEM_Calculator.cpp:916
belfem::fem::Calculator::create_matrix
calculator::MatrixData * create_matrix(const string &aLabel, const uint aNumRows, const uint aNumCols)
Definition cl_FEM_Calculator.cpp:2094
belfem::fem::Calculator::Ns
const Matrix< real > & Ns(const uint aIndex)
node interpolator slave
Definition cl_FEM_Calculator.hpp:2035
belfem::fem::Calculator::qswap
Vector< real > & qswap()
return a swap vector for the dofs
Definition cl_FEM_Calculator.cpp:2987
belfem::fem::Calculator::init_qold_table
void init_qold_table()
Definition cl_FEM_Calculator.cpp:2922
belfem::fem::Calculator::normal
const Vector< real > & normal(const uint aIndex=0)
returns the normal of a surface
Definition cl_FEM_Calculator.hpp:2082
belfem::fem::Calculator::group
Group * group()
Definition cl_FEM_Calculator.hpp:2439
belfem::fem::Calculator::num_intpoints
uint num_intpoints() const
Definition cl_FEM_Calculator.hpp:2400
belfem::fem::Calculator::set_model_parameters
void set_model_parameters(const Vector< real > &aParams)
Definition cl_FEM_Calculator.hpp:2415
belfem::fem::Calculator::bn_angle
real bn_angle(const Vector< real > &b, const Vector< real > &n, real &norm_b) const
Definition cl_FEM_Calculator.hpp:2538
belfem::fem::Calculator::volume_integration
const IntegrationData * volume_integration() const
Definition cl_FEM_Calculator.hpp:2344
belfem::fem::Calculator::B
const Matrix< real > & B(const uint aIndex=0)
gradient operator
Definition cl_FEM_Calculator.hpp:2074
belfem::fem::Calculator::C
const Matrix< real > & C(const uint aIndex)
curl interpolator
Definition cl_FEM_Calculator.cpp:3020
belfem::fem::Calculator::E
const Matrix< real > & E(const uint aIndex)
edge interpolator
Definition cl_FEM_Calculator.cpp:3008
belfem::fem::Calculator::JN
Matrix< real > & JN()
return the Newton correction matrix
Definition cl_FEM_Calculator.hpp:2384
belfem::fem::Calculator::Xs
const Matrix< real > & Xs() const
node coordinates on slave element
Definition cl_FEM_Calculator.hpp:1704
belfem::fem::Calculator::Jm
const Matrix< real > & Jm(const uint aIndex)
jacobian matrix for master
Definition cl_FEM_Calculator.hpp:1728
belfem::fem::Calculator::element
Element * element()
Definition cl_FEM_Calculator.hpp:2431
belfem::fem::Calculator::print_local_dofs
void print_local_dofs()
Definition cl_FEM_Calculator.cpp:2136
belfem::fem::Calculator::node_data
const Vector< real > & node_data(const string &aNodeField)
return a vector object
Definition cl_FEM_Calculator.cpp:2884
belfem::fem::Calculator::volume_enrichment
const IntegrationData * volume_enrichment() const
Definition cl_FEM_Calculator.hpp:2352
belfem::fem::Calculator::nedelec_data_master_h
const Vector< real > & nedelec_data_master_h()
edge dofs of the master volume of the linked facet, read from the edge_h mesh field ( the live dof st...
Definition cl_FEM_Calculator.cpp:2171
belfem::fem::Calculator::bj_angle
real bj_angle(const Vector< real > &b, const Vector< real > &j, real &norm_b, real &norm_j) const
Definition cl_FEM_Calculator.hpp:2532
belfem::fem::Calculator::vector_exists
bool vector_exists(const string aLabel) const
Definition cl_FEM_Calculator.hpp:1657
belfem::fem::Calculator::num_nedelec_dofs
uint num_nedelec_dofs() const
Definition cl_FEM_Calculator.cpp:3044
belfem::fem::Calculator::vector
Vector< real > & vector(const string &aLabel)
return a vector object
Definition cl_FEM_Calculator.hpp:1665
belfem::fem::Calculator::J
const Matrix< real > & J(const uint aIndex)
jacobian matrix
Definition cl_FEM_Calculator.hpp:1712
belfem::fem::Calculator::qold
const Vector< real > & qold(const uint aStep=0)
return the dof vector at an old timestep
Definition cl_FEM_Calculator.cpp:2894
belfem::fem::Calculator::nedelec_data_h
const Vector< real > & nedelec_data_h()
Definition cl_FEM_Calculator.hpp:2447
belfem::fem::Calculator::X
const Matrix< real > & X() const
node coordinates on element
Definition cl_FEM_Calculator.hpp:1688
belfem::fem::Calculator::nedelec_data_a
const Vector< real > & nedelec_data_a()
Definition cl_FEM_Calculator.hpp:2455
belfem::fem::Calculator::edge_function
EdgeFunction * edge_function()
Definition cl_FEM_Calculator.hpp:1339
belfem::fem::Calculator::mesh
Mesh * mesh()
Definition cl_FEM_Calculator.hpp:2519
belfem::fem::Calculator::Bs
const Matrix< real > & Bs(const uint aIndex=0)
gradient operator slave
Definition cl_FEM_Calculator.hpp:2051
belfem::fem::Calculator::Calculator
Calculator(Group *aGroup, const ModelDimensionality aDimensionality)
default constructor
Definition cl_FEM_Calculator.cpp:825
belfem::fem::Calculator::Bm
const Matrix< real > & Bm(const uint aIndex=0)
gradient operator master
Definition cl_FEM_Calculator.hpp:2043
belfem::fem::Calculator::f
Vector< real > & f()
return the load vector
Definition cl_FEM_Calculator.hpp:2392
belfem::fem::Calculator::dV
real dV(const uint aIndex=0)
volume increment
Definition cl_FEM_Calculator.hpp:2104
belfem::fem::Calculator::sideset_enrichment
const IntegrationData * sideset_enrichment() const
Definition cl_FEM_Calculator.hpp:2360
belfem::fem::Calculator::Js
const Matrix< real > & Js(const uint aIndex)
jacobian matrix for slave
Definition cl_FEM_Calculator.hpp:1743
belfem::fem::Calculator::maxwell
calculator::MaxwellData * maxwell()
Definition cl_FEM_Calculator.hpp:4194
belfem::fem::Calculator::model_parameters
const Vector< real > & model_parameters() const
Definition cl_FEM_Calculator.hpp:2423
belfem::fem::Calculator::create_vector
calculator::VectorData * create_vector(const string &aLabel, const uint aSize, const EntityType aType=EntityType::UNDEFINED)
Definition cl_FEM_Calculator.cpp:2063
belfem::fem::Calculator::q
const Vector< real > & q()
return the dof vector at current timestep
Definition cl_FEM_Calculator.cpp:2995
belfem::fem::Calculator::N
const Matrix< real > & N(const uint aIndex)
node interpolator
Definition cl_FEM_Calculator.hpp:2019
belfem::fem::Calculator::get_normal_calculator
Calculator * get_normal_calculator(Vector< real > &aPhiM, Vector< real > &aPhiS, bool &aMasterIsConductor, bool &aSlaveIsConductor)
links the tape-sideset calculator to the facet of the linked layer element and returns it.
Definition cl_FEM_Calculator.cpp:1985
belfem::fem::EdgeFunction
the edge function base class
Definition cl_EF_EdgeFunction.hpp:32
belfem::fem::Element
Definition cl_FEM_Element.hpp:41
belfem::fem::Group
Definition cl_FEM_Group.hpp:45
belfem::fem::IntegrationData
Precomputed bundle of integration points, weights and shape-function values.
Definition cl_IF_IntegrationData.hpp:28
belfem::fem::IntegrationData::weights
const Vector< real > & weights() const
return the integration weights
Definition cl_IF_IntegrationData.hpp:266
belfem::fem::Kernel
Top-level orchestrator; owns the mesh, materials, boundary conditions and DOF managers.
Definition cl_FEM_Kernel.hpp:50
belfem::fem::calculator::MatrixData
Definition cl_FEM_Calculator.hpp:79
belfem::fem::calculator::MatrixData::set_index
void set_index(const uint aIndex)
Definition cl_FEM_Calculator.hpp:1636
belfem::fem::calculator::MatrixData::label
const string & label() const
Definition cl_FEM_Calculator.hpp:1630
belfem::fem::calculator::MatrixData::index
uint index() const
Definition cl_FEM_Calculator.hpp:1642
belfem::fem::calculator::MatrixData::matrix
Matrix< real > & matrix()
Definition cl_FEM_Calculator.hpp:1650
belfem::fem::calculator::MatrixData::MatrixData
MatrixData(const string &aLabel, const uint aNumRows, const uint aNumCols)
Definition cl_FEM_Calculator.cpp:56
belfem::fem::calculator::MaxwellData
Definition cl_FEM_Calculator.hpp:132
belfem::fem::calculator::MaxwellData::compute_lambda
real compute_lambda(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2958
belfem::fem::calculator::MaxwellData::compute_cp
real compute_cp(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2933
belfem::fem::calculator::MaxwellData::compute_dlambdadT
real compute_dlambdadT(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2969
belfem::fem::calculator::MaxwellData::compute_drhodbeta
real compute_drhodbeta(const uint aIndex)
Definition cl_FEM_Calculator.hpp:3610
belfem::fem::calculator::MaxwellData::compute_h
const Vector< real > & compute_h(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2792
belfem::fem::calculator::MaxwellData::compute_b
const Vector< real > & compute_b(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2838
belfem::fem::calculator::MaxwellData::T_clamped
bool T_clamped() const
Definition cl_FEM_Calculator.hpp:4170
belfem::fem::calculator::MaxwellData::maxwell
Calculator * maxwell()
Definition cl_FEM_Calculator.hpp:4182
belfem::fem::calculator::MaxwellData::density
real density() const
Definition cl_FEM_Calculator.hpp:4161
belfem::fem::calculator::MaxwellData::reset
void reset()
Definition cl_FEM_Calculator.hpp:2770
belfem::fem::calculator::MaxwellData::compute_drhodj
real compute_drhodj(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2886
belfem::fem::calculator::MaxwellData::compute_T
real compute_T(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2927
belfem::fem::calculator::MaxwellData::compute_drhodb
real compute_drhodb(const uint aIndex)
Definition cl_FEM_Calculator.hpp:3595
belfem::fem::calculator::MaxwellData::compute_dmudh
real compute_dmudh(const uint aIndex)
Definition cl_FEM_Calculator.hpp:4096
belfem::fem::calculator::MaxwellData::norm_b
real norm_b(const uint aIndex)
Definition cl_FEM_Calculator.hpp:3107
belfem::fem::calculator::MaxwellData::compute_j
const Vector< real > & compute_j(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2626
belfem::fem::calculator::MaxwellData::compute_x
void compute_x(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2599
belfem::fem::calculator::MaxwellData::compute_rho
real compute_rho(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2867
belfem::fem::calculator::MaxwellData::compute_dcpdT
real compute_dcpdT(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2944
belfem::fem::calculator::MaxwellData::rho_clamped
bool rho_clamped() const
Definition cl_FEM_Calculator.hpp:4176
belfem::fem::calculator::MaxwellData::compute_drhodT
real compute_drhodT(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2901
belfem::fem::calculator::MaxwellData::thermal
Calculator * thermal()
Definition cl_FEM_Calculator.hpp:4188
belfem::fem::calculator::MaxwellData::norm_j
real norm_j(const uint aIndex)
Definition cl_FEM_Calculator.hpp:3118
belfem::fem::calculator::MaxwellData::compute_mu
real compute_mu(const uint aIndex)
Definition cl_FEM_Calculator.hpp:2804
belfem::fem::calculator::MaxwellData::MaxwellData
MaxwellData(Calculator *aCalculator, Kernel *aMaxwellKernel, Kernel *aThermalKernel)
Definition cl_FEM_Calculator.cpp:66
belfem::fem::calculator::MaxwellData::compute_volumetric_heatload
real compute_volumetric_heatload(const uint aIndex)
artificial volumetric heat load [ W/m³ ] at integration point aIndex, from the material's heating plu...
Definition cl_FEM_Calculator.hpp:4128
belfem::fem::calculator::VectorData
Definition cl_FEM_Calculator.hpp:49
belfem::fem::calculator::VectorData::label
const string & label() const
Definition cl_FEM_Calculator.hpp:1600
belfem::fem::calculator::VectorData::vector
Vector< real > & vector()
Definition cl_FEM_Calculator.hpp:1618
belfem::fem::calculator::VectorData::entity_type
EntityType entity_type() const
Definition cl_FEM_Calculator.hpp:1624
belfem::fem::calculator::VectorData::VectorData
VectorData(const string &aLabel, const uint aSize, const EntityType aType)
Definition cl_FEM_Calculator.cpp:45
belfem::fem::calculator::VectorData::set_index
void set_index(const uint aIndex)
Definition cl_FEM_Calculator.hpp:1606
belfem::fem::calculator::VectorData::index
uint index() const
Definition cl_FEM_Calculator.hpp:1612
belfem::mesh::Facet
Definition cl_Facet.hpp:24
en_IWGs.hpp
fn_det.hpp
Determinant of a square matrix.
fn_dot.hpp
Inner products of vectors and matrices.
fn_inv2.hpp
Closed-form inverse of a 2x2 matrix.
fn_inv3.hpp
Closed-form inverse of a 3x3 matrix.
fn_norm.hpp
Euclidean length of a vector.
belfem::dot
auto dot(const Vector< T > &aA, const Vector< T > &aB) -> decltype(arma::dot(aA.vector_data(), aB.vector_data()))
Scalar product of two vectors.
Definition fn_AR_dot.hpp:24
belfem::norm
auto norm(const Vector< T > &aA) -> decltype(norm(aA.vector_data()))
Euclidean (L2) norm of a vector.
Definition fn_norm.hpp:56
belfem::constant::phi
const real phi
golden number
Definition constants.hpp:48
belfem::constant::mu0
const real mu0
magnetic constant in V*s/(A*m) ( recommended value ) http://physics.nist.gov/cgi-bin/cuu/Value?...
Definition constants.hpp:73
belfem::constant::pi
const real pi
circle number
Definition constants.hpp:41
belfem::fem::calculator
Definition cl_FEM_Calculator.cpp:42
belfem::fem::calculator::MaxwellDataValue
MaxwellDataValue
Definition cl_FEM_Calculator.hpp:106
belfem::fem::calculator::MaxwellDataValue::drhodbeta
@ drhodbeta
Definition cl_FEM_Calculator.hpp:127
belfem::fem::calculator::MaxwellDataValue::dlambdadT
@ dlambdadT
Definition cl_FEM_Calculator.hpp:122
belfem::fem::calculator::MaxwellDataValue::j
@ j
Definition cl_FEM_Calculator.hpp:110
belfem::fem::calculator::MaxwellDataValue::normJ
@ normJ
Definition cl_FEM_Calculator.hpp:117
belfem::fem::calculator::MaxwellDataValue::drhodT
@ drhodT
Definition cl_FEM_Calculator.hpp:120
belfem::fem::calculator::MaxwellDataValue::n
@ n
Definition cl_FEM_Calculator.hpp:118
belfem::fem::calculator::MaxwellDataValue::rho
@ rho
Definition cl_FEM_Calculator.hpp:111
belfem::fem::calculator::MaxwellDataValue::normH
@ normH
Definition cl_FEM_Calculator.hpp:115
belfem::fem::calculator::MaxwellDataValue::mu
@ mu
Definition cl_FEM_Calculator.hpp:124
belfem::fem::calculator::MaxwellDataValue::dmudh
@ dmudh
Definition cl_FEM_Calculator.hpp:125
belfem::fem::calculator::MaxwellDataValue::lambda
@ lambda
Definition cl_FEM_Calculator.hpp:113
belfem::fem::calculator::MaxwellDataValue::beta
@ beta
Definition cl_FEM_Calculator.hpp:114
belfem::fem::calculator::MaxwellDataValue::cp
@ cp
Definition cl_FEM_Calculator.hpp:112
belfem::fem::calculator::MaxwellDataValue::dcpdT
@ dcpdT
Definition cl_FEM_Calculator.hpp:121
belfem::fem::calculator::MaxwellDataValue::B
@ B
Definition cl_FEM_Calculator.hpp:109
belfem::fem::calculator::MaxwellDataValue::x
@ x
Definition cl_FEM_Calculator.hpp:123
belfem::fem::calculator::MaxwellDataValue::drhodj
@ drhodj
Definition cl_FEM_Calculator.hpp:119
belfem::fem::calculator::MaxwellDataValue::normB
@ normB
Definition cl_FEM_Calculator.hpp:116
belfem::fem::calculator::MaxwellDataValue::T
@ T
Definition cl_FEM_Calculator.hpp:107
belfem::fem::calculator::MaxwellDataValue::H
@ H
Definition cl_FEM_Calculator.hpp:108
belfem::fem::calculator::MaxwellDataValue::drhodb
@ drhodb
Definition cl_FEM_Calculator.hpp:126
belfem::fem
Definition cl_IFB_LINE3.hpp:21
belfem::fem::compute_hn
const Vector< real > & compute_hn(Calculator *aCalc, const uint k)
Computes the purely-normal magnetic field hn at the thin-shell facet: the average of the master and s...
Definition cl_FEM_Calculator.hpp:2713
belfem::fem::compute_h_trace
void compute_h_trace(Calculator *aNormalCalc, const bool aMaster, const bool aIsConductor, const Vector< real > &aPhi, Vector< real > &aScratch, Vector< real > &aH)
Full field trace of one volume side of a thin-shell facet, for linear elements ( the k = 0 contract o...
Definition cl_FEM_Calculator.hpp:2669
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::InterpolationType
InterpolationType
Definition Mesh_Enums.hpp:99
belfem::gTbulk
greal gTbulk
Definition globals.hpp:45
belfem::uint
unsigned int uint
Definition typedefs.hpp:30
belfem::ModelDimensionality
ModelDimensionality
Definition en_IWGs.hpp:20
belfem::gRhoMin
greal gRhoMin
minimim resistrivity in Ohm*m, default: 0
Definition globals.hpp:48
belfem::ElementType
ElementType
Element types.
Definition Mesh_Enums.hpp:27
belfem::gTmin
constexpr real gTmin
Definition typedefs.hpp:65
belfem::det
auto det(const T &A) -> decltype(arma::det(A))
Definition fn_AR_det.hpp:23
belfem::EntityType
EntityType
Definition Mesh_Enums.hpp:116
belfem::EntityType::UNDEFINED
@ UNDEFINED
Definition Mesh_Enums.hpp:124
belfem::MaterialProperty::q
@ q
Definition cl_Material.hpp:176
belfem::MaterialProperty::mu
@ mu
Definition cl_Material.hpp:159
belfem::gRhoMax
greal gRhoMax
maximim resistrivity in Ohm*m, default: 1e10
Definition globals.hpp:51
belfem::index_t
uint32_t index_t
Definition typedefs.hpp:52
belfem::real
double real
Definition typedefs.hpp:36
mesh.Element
Definition Element.py:1
mesh
Definition __init__.py:1
BELFEM_UINT_MAX
#define BELFEM_UINT_MAX
Definition typedefs.hpp:80
BELFEM_REAL_MAX
#define BELFEM_REAL_MAX
Definition typedefs.hpp:81
BELFEM_QUIET_NAN
#define BELFEM_QUIET_NAN
Definition typedefs.hpp:87