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