Microsoft PowerPoint - FEM3D-C.ppt [互換モード]

Transcription

1 3D-FEM n C: Stead State Heat Conducton Kengo aajma Informaton echnolog Center Programmng for Parallel Computng ( ) Semnar on Advanced Computng (66-009)

2 FEM3D 2 3D Stead-State Heat Conducton Z =0@Z= ma Z Q 0 Heat Generaton Unform thermal conductvt HEX meshes cubes X Y Z cubes n each drecton Boundar Condtons X Y X Y =0@Z= ma Heat Gen. Rate s a functon of locaton (cell center: c c ) Q QVOL C C

3 FEM3D 3 3D Stead-State Heat Conducton Q 0 Hgher temperature at nodes far from the orgn. Heat Gen. Rate s a functon of locaton (cell center: c c ) Q C C move

4 FEM3D Fnte-Element Procedures Governng Equatons Galern Method: Wea Form Element-b-Element Integraton Element Matr Global Matr Boundar Condtons Lnear Solver

5 FEM3D 5 FEM Procedures: Program Intalaton Control Data ode Connectvt of Elements (: ode# E: Elem#) Intalaton of Arras (Global/Element Matrces) Element-Global Matr Mappng (Inde Item) Generaton of Matr Element-b-Element Operatons (do cel= E) Element matrces Accumulaton to global matr Boundar Condtons Lnear Solver Conjugate Gradent Method

6 FEM3D 6 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program

7 FEM3D 7 Etenson to 2D Prob.: rangles 三角形要素 rangles can handle arbtrarl shaped object Lnear trangular elements provde low accurac therefore the are not used n practcal applcatons.

8 FEM3D Etenson to 2D Prob.: Quadrlaterals 四角形要素 Formulaton of quad. elements s possble f same shape functons n D elements are appled along X- and Y- as. More accurate than trangles Each edge must be parallel wth X- and Y- as. Smlar to FDM 3 hs tpe of elements cannot be consdered. 2

9 FEM3D 9 Isoparametrc Element (/3) Each element s mapped to square element [±±] on natural/local coordnate () Components of global coordnate sstem of each node () for certan nds of elements are defned b shape functons [] on natural/local coordnate sstem where shape functons [] are also used for nterpolaton of dependent varables.

10 0 Isoparametrc Element (2/3) 2 3 Coordnate of each node: ( ) ( 2 2 ) ( 3 3 ) ( ) emperature at each node: FEM3D ) ( ) ( ) (

11 FEM3D Isoparametrc Element (3/3) Hgher-order shape functon can handle curved lnes/surfaces. atural coordnate sstem Sub-Parametrc Super-Parametrc

12 2 Shape Fn s on 2D atural Coord. (/3) Polnomal shape functons on squares of natural coordnate: Coeffcents are calculated as follows: FEM3D

13 Shape Fn s on 2D atural Coord. (2/3) FEM3D

14 Shape Fn s on 2D atural Coord. (2/3) FEM3D

15 5 Shape Fn s on 2D atural Coord. (3/3) s defned as follows accordng to : ) ( ) ( ) ( ) ( 3 2 Shape functons : Also nown as b-lnear nterpolaton Calculate at each node FEM3D

16 FEM3D 6 Etenson to 3D Problems etrahedron/etrahedra( 四面体 ): rangles n 2D can handle arbtrar shape objects Lnear elements are generall less accurate not practcal Hgher-order tetrahedral elements are wdel used. In ths class tr-lnear heahedral elements (soparametrc) are used( 六面体要素 )

17 7 Shape Fn s: 3D atural/local Coord. ) ( ) ( ) ( ) ( 3 2 ) ( ) ( ) ( ) ( FEM3D ) ( ) ( ) ( ) (

18 FEM3D Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program

19 Galern Method (/3) Governng Equaton for 3D Stead State Heat Conducton Problems (unform ): { Followng ntegral equaton s obtaned at each element b Galern method where [] s are also weghtng functons: dv Q V Q FEM3D 9 Dstrbuton of temperature n each element (matr form) : emperature at each node

20 Galern Method (2/3) Green s heorem (3D) Appl ths to the st 3-parts of the equaton wth 2 nd - order dff. (surface ntegraton terms are gnored): Consder the followng terms: { { { { V V S dv B A B A B A ds n B A dv B B B A dv dv V V FEM3D 20

21 FEM3D 2 Galern Method (3/3) Fnall followng equaton s obtaned b consderng heat generaton term : Q dv Q dv 0 V hs s called wea form( 弱形式 ). Orgnal PDE conssts of terms wth 2 nd -order dff. but ths wea form onl ncludes st -order dff b Green s theorem. Requrements for shape functons are weaer n wea form. Lnear functons can descrbe effects of 2 nd -order dfferentaton. Same as D problem V

22 Wea Form wth B.C.: on each elem. dv dv dv V V V e ) ( ) ( ) ( ) ( e e e f dv Q f V e ) ( FEM3D 22

23 FEM3D 23 Element Matr: j j j

24 FEM3D 2 Element Matr: j j j dv ( e) V V dv V dv j V j j j dv

25 FEM3D 25 et Stage: Integraton

26 FEM3D 26 Methods for umercal Integraton rapeodal Rule Smpson s Rule Gaussan Quadrature (or Gauss-Legendre) accurate Values of functons at fnte numbers of sample ponts are utled: X X 2 f m d w f

27 FEM3D 27 Gaussan Quadrature n D more accurate than Smpson s rule Smpson s f Gauss ( ) sn( ) X h S X X 2 X 3 0 X 2 X 3 2 X h 3 2 X X 3 X f X f X f X S h 2 2 / f d W f f h d

28 FEM3D 2 Gaussan Quadrature ガウスの積分公式 On normaled natural (or local) coordnate sstem [-+]( 自然座標系 局所座標系 ) Can appromate up to (2m-)-th order of functons b m quadrature ponts (m=2 s enough for quadratc shape functons). f m m m m d w f w w w / w.00 5/ 9 =- =0 = = =+

29 FEM3D 29 Gaussan Quadrature can be easl etended to 2D & 3D I m n W W j f ( j ) j f ( ) dd mn: number of quadrature ponts n -drecton ( j ): Coordnates of Quad s W W j : Weghtng Factor

30 30 hs confguraton s wdel used. In 2D problem ntegraton s done usng values of f at quad. ponts. Gaussan Quadrature ガウスの積分公式

31 3 hs confguraton s wdel used. In 2D problem ntegraton s done usng values of f at quad. ponts. Gaussan Quadrature ガウスの積分公式 I f ( ) dd m n j W ( ).0.0 f ( ).0.0 f ( ).0.0 f ( ).0.0 f ( ) W j f j

32 et Stage: Integraton FEM3D 32 3D atural/local Coordnate (): Gaussan Quadrature j W W W ) ( j j j M j L f W W W d d d f I ) ( ) ( LM: number of quadrature ponts n -drecton : Weghtng Factor : Coordnates of Quad s

33 Partal Dff. on atural Coord. (/) Accordng to formulae: ) ( ) ( ) ( can be easl derved accordng to defntons. are requred for computatons. FEM3D 33

34 Partal Dff. on atural Coord. (2/) In matr form: J J : Jacob matr Jacoban FEM3D 3

35 Partal Dff. on atural Coord. (3/) Components of Jacoban: J J J J J J J J J FEM3D 35

36 Partal Dff. on atural Coord. (/) Partal dfferentaton on global coordnate sstem s ntroduced as follows (wth nverse of Jacoban matr (3 3)) J FEM3D 36

37 Integraton on Element V j j j V j j j j dv dv j j FEM3D 37

38 Integraton on atural Coord. d d d J ddd dv j j j j j j V j j j det j j FEM3D 3

39 Gaussan Quadrature j j M j L f W W W d d d f I ) ( ) ( d d d J j j j det FEM3D 39 j j

40 FEM3D 0 Remanng Procedures Element matrces have been formed. Accumulaton to Global Matr Implementaton of Boundar Condtons Solvng Lnear Equatons Detals of mplementaton wll be dscussed n classes later than net wee through eplanaton of programs

41 Accumulaton: Local -> Global Matrces () () 3 () 2 () () () 3 () 2 () () () 3 () 2 () () 3 () 33 () 32 () 3 () 2 () 23 () 22 () 2 () () 3 () 2 () () () () { ]{ [ f f f f f (2) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) 3 (2) 33 (2) 32 (2) 3 (2) 2 (2) 23 (2) 22 (2) 2 (2) (2) 3 (2) 2 (2) (2) (2) (2) { ]{ [ f f f f f { ]{ [ B B B B B B D X X X X D X X X X D X X X X X X D X X X X D X X X X D F K FEM3D

42 FEM3D 2 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Eercse Runnng the Code Data Structure Overvew of the Program

43 FEM3D 3 Eercse Develop a program and calculate area of the followng quadrlateral usng Gaussan Quadrature. V 3 2 :(.0.0) 2:(.0 2.0) 3:( ) :(2.0.0) I V dv det J dd

44 FEM3D ps (/2) Calculate Jacoban Appl Gaussan Quadrature (n=2) I m W W j f ( j ) f ( ) dd mplct REAL* (A-HO-Z) real* W(2) real* POI(2) W()=.0d0 W(2)=.0d0 POI()= d0 POI(2)= d0 n j SUM= 0.d0 do jp= 2 do p= 2 FC = F(POI(p)POI(jp)) SUM= SUM + W (p)*w (jp)*fc enddo enddo

45 FEM3D 5 ps (2/2) J J det ) ( ) ( ) ( ) ( 3 2

46 FEM3D 6 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program

47 FEM3D 7 3D Stead-State Heat Conducton Z =0@Z= ma Z Q 0 Heat Generaton Unform thermal conductvt HEX meshes cubes X Y Z cubes n each drecton Boundar Condtons X Y X Y =0@Z= ma Heat Gen. Rate s a functon of locaton (cell center: c c ) Q QVOL C C

48 FEM3D Cop/Installaton 3D-FEM Code >$ cd <$E-OP> >$ cp /home03/sengon/documents/class_eps/f/fem3d.tar. >$ cp /home03/sengon/documents/class_eps/c/fem3d.tar. >$ tar vf fem3d.tar >$ cd fem3d >$ ls run src Install >$ cd <$E-OP>/fem3d/src >$ mae >$ ls../run/sol sol Install of Mesh Generator >$ cd <$E-OP>/fem3d/run >$ g95 O3 mgcube.f o mgcube

49 FEM3D 9 Operatons Startng from Grd Generaton to Computaton Flenames are fed mgcube mesh generator cube.0 mesh fle sol FEM Solver IPU.DA Control Data test.np for vsualaton

50 FEM3D 50 Mesh Generaton Z >$ cd <$fem>/fem3d/run >$./mgcube =0@Z= ma X Y Z umber of Elem s n each drecton eample Z >$ ls cube.0 confrmaton cube.0 X Y X Y

51 FEM3D 5 Control Fle: IPU.DA IPU.DA cube.0 fname 2000 IER.0.0 COD QVOL.0e-0 RESID fname: IER: COD: QVOL: RESID: Q ame of Mesh Fle Ma. Iteratons for CG hermal Conductvt Heat Generaton Rate Crtera for Convergence of CG QVOL C C Q 0

52 FEM3D 52 >$ cd <$E-OP>/fem3d/run >$./sol Runnng >$ ls test.np Confrmaton test.np

53 FEM3D 53 ParaVew Openng fles Dsplang fgures Savng mage fles

54 FEM3D 5 UCD Format (/3) Unstructured Cell Data 要素の種類点線三角形四角形四面体角錐三角柱六面体二次要素線 2 三角形 2 四角形 2 四面体 2 角錐 2 三角柱 2 六面体 2 キーワード pt lne tr quad tet pr prsm he lne2 tr2 quad2 tet2 pr2 prsm2 he2

55 FEM3D 55 UCD Format (2/3) Orgnall for AVS mcroavs Etenson of the UCD fle s np here are two tpes of formats. Onl old tpe can be read b ParaVew.

56 FEM3D 56 UCD Format (3/3): Old Format ( 全節点数 ) ( 全要素数 ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( モデルのデータ数 ) ( 節点番号 ) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 要素データ成分 のラベル )( 単位 ) ( 要素データ成分 2 のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( 節点のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2 のラベル )( 単位 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ ) ( 節点データ 2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)

57 FEM3D 57 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program

58 FEM3D 5 odes Overvew of Mesh Fle: cube.0 numberng starts from ode # (How man nodes?) ode ID Coordnates Elements Element # Element pe Element ID Materal ID Connectvt ode Groups Group # ode # n each group Group ame odes n each group

59 FEM3D 59 Mesh Generator: mgcube.f (/5) mplct REAL* (A-HO-Z) real(nd=) dmenson(::) allocatable :: X Y real(nd=) dmenson(::) allocatable :: X0 Y0 character(len=0) :: GRIDFILE HHH nteger dmenson(::) allocatable :: IW!C!C !C II.!C !C=== wrte (**) 'XYZ' read (**) XYZ XP= X + YP= Y + ZP= Z + DX=.d0 IODO= XP*YP*ZP ICELO= X *Y *Z IBODO= XP*YP allocate (IW(IODO)) IW= 0 ode number n X-drecton ode number n Y-drecton ode number n Z-drecton otal number of odes otal number of Elements umber of nodes n XY-plane

60 FEM3D 60!C=== Mesh Generator: mgcube.f (2/5) cou= 0 b = do = ZP do j= YP = cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 2 do = ZP j= do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 = do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = = ZP do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo odes on X=Xmn surface are stored n IW(b) where b= YP*ZP X Z =0@Z= ma Y X Z Y

61 FEM3D 6!C=== Mesh Generator: mgcube.f (2/5) cou= 0 b = do = ZP do j= YP = cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 2 do = ZP j= do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 = do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = = ZP do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo odes on X=Xmn surface are stored n IW(b) where b= YP*ZP X Z =0@Z= ma Y X Z Y

62 FEM3D 62 Mesh Generator: mgcube.f (2/5)!C=== cou= 0 b = do = ZP do j= YP = cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 2 do = ZP j= do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 = do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = = ZP do j= YP do = XP cou= cou + = (-)*IBODO + (j-)*xp + IW(coub)= enddo enddo X Z =0@Z= ma Y X odes on Z=Zmn surface are stored n IW(b3) where b= XP*YP odes on Z=Zma surface are stored n IW(b) where b= XP*YP Z Y

63 FEM3D 63 Mesh Generator: mgcube.f (3/5)!C!C !C GeoFEM data!c !C=== = 0 wrte (**) 'GeoFEM grdfle name?' GRIDFILE= 'cube.0' open (2 fle= GRIDFILE status='unnown'form='formatted') wrte(2'(00)') IODO cou= 0 do = ZP do j= YP do = XP XX= dfloat(-)*dx YY= dfloat(j-)*dx ZZ= dfloat(-)*dx cou= cou + wrte (2'(03(pe6.6))') cou XX YY ZZ enddo enddo enddo wrte(2'(0)') ICELO IELMYPL= 36 wrte(2'(00)') (IELMYPL =ICELO) ode # ode ID Coordnates Element # Element pe (not n use)

64 FEM3D 6 Eample of cube.0 (X=Y=Z=) ode Element-pe 25 =5*5* E E E E E E E E E E E E E E E E E E E E E E E E E E E+00 ( 途中省略 ) E E E E E E E E E E E E E E E+00 6 =** move

65 FEM3D 65 Mesh Generator: mgcube.f (/5) cou= 0 mat= do = Z do j= Y do = X cou= cou + n = (-)*IBODO + (j-)*xp + n2 = n + n3 = n2 + XP n = n3 - n5 = n + IBODO n6 = n2 + IBODO n7 = n3 + IBODO n = n + IBODO wrte (2'(00)') cou mat n n2 n3 n & & n5 n6 n7 n enddo enddo enddo mat: Materal ID (=)

66 FEM3D 66 Eample of cube.0 (X=Y=Z=) Element Connectvt (...)

67 FEM3D 67 X=Y=Z= XP=YP=ZP=5 ICELO= 6 IODO= 25 IBODO= 25 = = = = =

68 FEM3D 6 Mesh Generator: mgcube.f (5/5) IGO= IB= YP*ZP IB2= XP*ZP + IB IB3= XP*YP + IB2 IB= XP*YP + IB3 wrte (2'(00)') IGO wrte (2'(00)') IB IB2 IB3 IB HHH= 'Xmn' wrte (2'(a0)') HHH wrte (2'(00)') (IW() =YP*ZP) HHH= 'Ymn' wrte (2'(a0)') HHH wrte (2'(00)') (IW(2) =XP*ZP) HHH= 'Zmn' wrte (2'(a0)') HHH wrte (2'(00)') (IW(3) =XP*YP) HHH= 'Zma' wrte (2'(a0)') HHH wrte (2'(00)') (IW() =XP*YP) IGO Group # (XmnYmnZmnZma) IB Accumulated # ( 以下略 )!C=== stop end deallocate (IW) close (2)

69 FEM3D 69 Eample of cube.0 (X=Y=Z=) Group Info. Xmn Ymn Zmn Zma (not n use after ths pont)

70 FEM3D 70 Mesh Generaton Bg echncal & Research Issue Complcated Geometr Large Scale Parallelaton s dffcult Commercal Mesh Generator FEMAP Interface to CAD Data Format move

71 FEM3D 7 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program

72 FEM3D 72 FEM Procedures: Program Intalaton Control Data ode Connectvt of Elements (: ode# E: Elem#) Intalaton of Arras (Global/Element Matrces) Element-Global Matr Mappng (Inde Item) Generaton of Matr Element-b-Element Operatons (do cel= E) Element matrces Accumulaton to global matr Boundar Condtons Lnear Solver Conjugate Gradent Method

73 FEM3D 73 test man nput_cntl nput of control data nput_grd nput of mesh nfo fnd_node searchng nodes mat_con0 connectvt of matr msor sortng mat_con connectvt of matr mat_ass_man coeffcent matr jacob Jacoban mat_ass_bc boundar condtons Structure of heat3d solve control of lnear solver output_ucd vsualaton cg CG solver

74 FEM3D 7 Man Part /** program heat3d **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPU_CL(); etern vod IPU_GRID(); etern vod MA_CO0(); etern vod MA_CO(); etern vod MA_ASS_MAI(); etern vod MA_ASS_BC(); etern vod SOLVE(); etern vod OUPU_UCD(); nt man() { IPU_CL(); IPU_GRID(); MA_CO0(); MA_CO(); MA_ASS_MAI(); MA_ASS_BC() ; SOLVE(); OUPU_UCD() ;

75 FEM3D 75 Global Varables: pfem_utl.h (/3) ame pe Se I/O Defnton fname C [0] I ame of mesh fle P I I # ode ICELO I I # Element ODGRPtot I I # ode Group XYZ R [][3] I ode Coordnates ICELOD I [ICELO][] I Element Connectvt ODGRP_IDEX I [ODGRPtot+] I # ode n each ode Group ODGRP_IEM I [ODGRP_IDEX[ ODGRPO+]] I ode ID n each ode Group ODGRP_AME C0 [ODGRP_IDEX[ ODGRPO+]] I ame of odegroup LU I O # on-zero Off-Dagonals at each node PLU I O # on-zero Off-Dagonals D R [] O Dagonal Bloc of Global Matr BX R [] O RHS Unnown Vector

76 FEM3D 76 Global Varables: pfem_utl.h (2/3) ame pe Se I/O Defnton AMA R [PLU] O on-zero Off-Dagonal Components of Global Matr ndelu I [+] O # on-zero Off-Dagonal Components temlu I [PLU] O ILU I [] O IALU I [][LU] O Column ID of on-zero Off-Dagonal Components umber of on-zero Off-Dagonal Components at Each ode Column ID of on-zero Off-Dagonal Components at Each ode IWKX I [][2] O Wor Arras IER IERactual I I umber of CG Iteratons (MAX Actual) RESID R I Convergence Crtera (fed as.e-) pfemiarra I [00] O Integer Parameter Arra pfemrarra R [00] O Real Parameter Arra

77 FEM3D 77 Global Varables: pfem_utl.h (3/3) ame pe Se I/O Defnton Oth R I = 0.25 PQ PE P R [2][2][] O ~ at each Gaussan Quad. Pont POS WEI R [2] O COL COL2 I [00] O Wor arras for sortng Coordnates Weghtng Factor at each Gaussan Quad. Pont SHAPE R [2][2][2][] O (=~) at each Gaussan Quad Pont PX PY PZ R [2][2][2][] O ~ at each Gaussan Quad. Pont DEJ R [2][2][2] O Determnant of Jacoban Matr at each Gaussan Quad. Pont COD QVOL R I hermal Conductvt Heat Generaton Rate Q QVOL C C Q 0

78 FEM3D 7 IPU_CL: Control Data /** ** IPU_CL **/ #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" /** **/ vod IPU_CL() { FILE *fp; f( (fp=fopen("ipu.da""r")) == ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); fscanf(fp"%s"fname); fscanf(fp"%d"&ier); fscanf(fp "%lf %lf" &COD &QVOL); fscanf(fp "%lf" &RESID); fclose(fp); pfemrarra[0]= RESID; pfemiarra[0]= IER; IPU.DA cube.0 fname 2000 IER.0.0 COD QVOL.0e-0 RESID

79 FEM3D 79 IPU_GRID (/2) #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" #nclude "allocate.h" vod IPU_GRID() { FILE *fp; nt jnncelse; nt YPEIMA; /** **/ f( (fp=fopen(fname"r")) == ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); ODE fscanf(fp"%d"&); P=; XYZ=(KREAL**)allocate_matr(seof(KREAL)3); for(=0;<;++){ for(j=0;j<3;j++){ XYZ[][j]=0.0; for(=0;<;++){ fscanf(fp"%d %lf %lf %lf"&&xyz[][0]&xyz[][]&xyz[][2]);

80 FEM3D 0 allocate deallocate #nclude <stdo.h> #nclude <stdlb.h> vod* allocate_vector(nt sent m) { vod *a; f ( ( a=(vod * )malloc( m * se ) ) == ULL ) { fprntf(stdout"error:memor does not enough! n vector n"); et(); return a; vod deallocate_vector(vod *a) { free( a ); vod** allocate_matr(nt sent mnt n) { vod **aa; nt ; f ( ( aa=(vod ** )malloc( m * seof(vod*) ) ) == ULL ) { fprntf(stdout"error:memor does not enough! aa n matr n"); et(); f ( ( aa[0]=(vod * )malloc( m * n * se ) ) == ULL ) { fprntf(stdout"error:memor does not enough! n matr n"); et(); for(=;<m;++) aa[]=(char*)aa[-]+se*n; return aa; vod deallocate_matr(vod **aa) { free( aa ); Same nterface wth FORRA

81 FEM3D IPU_GRID (/2) /** **/ ELEME fscanf(fp"%d"&icelo); ICELOD=(KI**)allocate_matr(seof(KI)ICELO); for(=0;<icelo;++) fscanf(fp"%d"&ype); ICELOD[][j]: ode ID startng from. Element ID starts from 0. /** **/ for(cel=0;cel<icelo;cel++){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&ima &ICELOD[cel][0]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] &ICELOD[cel][]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ODGRP_IDEX=(KI* )allocate_vector(seof(ki)odgrptot+); ODGRP_AME =(CHAR0*)allocate_vector(seof(CHAR0)ODGRPtot); for(=0;<odgrptot+;++) ODGRP_IDEX[]=0; for(=0;<odgrptot;++) fscanf(fp"%d"&odgrp_idex[+]); nn=odgrp_idex[odgrptot]; ODGRP_IEM=(KI*)allocate_vector(seof(KI)nn); for(=0;<odgrptot;++){ S= ODGRP_IDEX[]; E= ODGRP_IDEX[+]; fscanf(fp"%s"odgrp_ame[].name); nn= E - S; f( nn!= 0 ){ for(=s;<e;++) fscanf(fp"%d"&odgrp_iem[]); ode Group: ode ID s start from fclose(fp);

82 FEM3D 2 test man nput_cntl nput of control data nput_grd nput of mesh nfo fnd_node searchng nodes mat_con0 connectvt of matr msor sortng mat_con connectvt of matr mat_ass_man coeffcent matr jacob Jacoban mat_ass_bc boundar condtons Structure of heat3d solve control of lnear solver output_ucd vsualaton cg CG solver

83 FEM3D 3 Global Varables: pfem_utl.h (/3) ame pe Se I/O Defnton fname C [0] I ame of mesh fle P I I # ode ICELO I I # Element ODGRPtot I I # ode Group XYZ R [][3] I ode Coordnates ICELOD I [ICELO][] I Element Connectvt ODGRP_IDEX I [ODGRPtot+] I # ode n each ode Group ODGRP_IEM I [ODGRP_IDEX[ ODGRPO+]] I ode ID n each ode Group ODGRP_AME C0 [ODGRP_IDEX[ ODGRPO+]] I ame of odegroup LU I O # on-zero Off-Dagonals at each node PLU I O # on-zero Off-Dagonals D R [] O Dagonal Bloc of Global Matr BX R [] O RHS Unnown Vector

84 FEM3D Global Varables: pfem_utl.h (2/3) ame pe Se I/O Defnton AMA R [PLU] O on-zero Off-Dagonal Components of Global Matr ndelu I [+] O # on-zero Off-Dagonal Components temlu I [PLU] O ILU I [] O IALU I [][LU] O Column ID of on-zero Off-Dagonal Components umber of on-zero Off-Dagonal Components at Each ode Column ID of on-zero Off-Dagonal Components at Each ode IWKX I [][2] O Wor Arras IER IERactual I I umber of CG Iteratons (MAX Actual) RESID R I Convergence Crtera (fed as.e-) pfemiarra I [00] O Integer Parameter Arra pfemrarra R [00] O Real Parameter Arra

85 FEM3D 5 Global Varables: pfem_utl.h (3/3) ame pe Se I/O Defnton Oth R I = 0.25 PQ PE P R [2][2][] O ~ at each Gaussan Quad. Pont POS WEI R [2] O COL COL2 I [00] O Wor arras for sortng Coordnates Weghtng Factor at each Gaussan Quad. Pont SHAPE R [2][2][2][] O (=~) at each Gaussan Quad Pont PX PY PZ R [2][2][2][] O ~ at each Gaussan Quad. Pont DEJ R [2][2][2] O Determnant of Jacoban Matr at each Gaussan Quad. Pont COD QVOL R I hermal Conductvt Heat Generaton Rate Q QVOL C C Q 0

86 FEM3D 6 owards Matr Assemblng In D t was eas to obtan nformaton related to nde and tem. 2 non-ero off-dagonals for each node ID of non-ero off-dagonal : + - where s node ID In 3D stuaton s more complcated: umber of non-ero off-dagonal components s between 7 and 26 for the current target problem More complcated for real problems. Generall there are no nformaton related to number of non-ero off-dagonal components beforehand. move

87 FEM3D 7 owards Matr Assemblng In D t was eas to obtan nformaton related to nde and tem. 2 non-ero off-dagonals for each node ID of non-ero off-dagonal : + - where s node ID In 3D stuaton s more complcated: umber of non-ero off-dagonal components s between 7 and 26 for the current target problem More complcated for real problems. Generall there are no nformaton related to number of non-ero off-dagonal components beforehand. Count number of non-ero off-dagonals usng arras: ILU[] IALU[][LU]

88 FEM3D Man Part /** program heat3d **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPU_CL(); etern vod IPU_GRID(); etern vod MA_CO0(); etern vod MA_CO(); etern vod MA_ASS_MAI(); etern vod MA_ASS_BC(); etern vod SOLVE(); etern vod OUPU_UCD(); nt man() { IPU_CL(); IPU_GRID(); MA_CO0(); MA_CO(); MA_ASS_MAI(); MA_ASS_BC() ; SOLVE(); OUPU_UCD() ; MA_CO0: generates IU IALU MA_CO: generates nde tem ode ID startng from

89 FEM3D 9 MA_CO0: Overvew do cel= ICELO generate ILU IALU accordng to nodes of he. elements (FID_ODE) enddo

90 FEM3D 90 Generatng Connectvt of Matr MA_CO0 (/) /** ** MA_CO0 **/ #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h etern FILE *fp_log; /*** eternal functons ***/ etern vod msor(nt* nt* nt); /*** statc functuons ***/ statc vod FID_S_ODE (ntnt); vod MA_CO0() { nt jceln; nt nn2n3nn5n6n7n; nt ; LU= 26; ILU=(KI* )allocate_vector(seof(ki)); IALU=(KI**)allocate_matr(seof(KI)LU); for(=0;<;++) ILU[]=0; for(=0;<;++) for(j=0;j<lu;j++) IALU[][j]=0; LU: umber of mamum number of connected nodes to each node (number of upper/lower non-ero off-dagonal nodes) In the current problem geometr s rather smple. herefore we can specf LU n ths wa. If t s not clear -> r more fleble mplementaton

91 FEM3D 9 Generatng Connectvt of Matr MA_CO0 (/) /** ** MA_CO0 **/ #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h etern FILE *fp_log; /*** eternal functons ***/ etern vod msor(nt* nt* nt); /*** statc functuons ***/ statc vod FID_S_ODE (ntnt); vod MA_CO0() { nt jceln; nt nn2n3nn5n6n7n; nt ; Arra Se Descrpton ILU IALU [] [][LU] umber of connected nodes to each node (lower/upper) Correspondng connected node ID (column ID) LU= 26; ILU=(KI* )allocate_vector(seof(ki)); IALU=(KI**)allocate_matr(seof(KI)LU); for(=0;<;++) ILU[]=0; for(=0;<;++) for(j=0;j<lu;j++) IALU[][j]=0;

92 FEM3D 92 Generatng Connectvt of Matr MA_CO0 (2/): Startng from for( cel=0;cel< ICELO;cel++){ n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; FID_S_ODE (nn2); FID_S_ODE (nn3); FID_S_ODE (nn); FID_S_ODE (nn5); FID_S_ODE (nn6); FID_S_ODE (nn7); FID_S_ODE (nn); FID_S_ODE (n2n); FID_S_ODE (n2n3); FID_S_ODE (n2n); FID_S_ODE (n2n5); FID_S_ODE (n2n6); FID_S_ODE (n2n7); FID_S_ODE (n2n); FID_S_ODE (n3n); FID_S_ODE (n3n2); FID_S_ODE (n3n); FID_S_ODE (n3n5); FID_S_ODE (n3n6); FID_S_ODE (n3n7); FID_S_ODE (n3n);

93 FEM3D 93 FID_S_ODE: Search Connectvt ILUIALU: Automatc Search /*** *** FID_S_ODE **/ statc vod FID_S_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; Arra Se Descrpton ILU IALU [] [][LU] umber of connected nodes to each node (lower/upper) Correspondng connected node ID (column ID)

94 FEM3D 9 FID_S_ODE: Search Connectvt ILUIALU: Automatc Search /*** *** FID_S_ODE **/ statc vod FID_S_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; If the target node s alread ncluded n IALU proceed to net par of nodes return;

95 FEM3D 95 FID_S_ODE: Search Connectvt ILUIALU: Automatc Search /*** *** FID_S_ODE **/ statc vod FID_S_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; If the target node s O ncluded n IALU store the node n IALU and add to ILU

96 FEM3D 96 Generatng Connectvt of Matr MA_CO0 (3/) FID_S_ODE (nn); FID_S_ODE (nn2); FID_S_ODE (nn3); FID_S_ODE (nn5); FID_S_ODE (nn6); FID_S_ODE (nn7); FID_S_ODE (nn); FID_S_ODE (n5n); FID_S_ODE (n5n2); FID_S_ODE (n5n3); FID_S_ODE (n5n); FID_S_ODE (n5n6); FID_S_ODE (n5n7); FID_S_ODE (n5n); FID_S_ODE (n6n); FID_S_ODE (n6n2); FID_S_ODE (n6n3); FID_S_ODE (n6n); FID_S_ODE (n6n5); FID_S_ODE (n6n7); FID_S_ODE (n6n); FID_S_ODE (n7n); FID_S_ODE (n7n2); FID_S_ODE (n7n3); FID_S_ODE (n7n); FID_S_ODE (n7n5); FID_S_ODE (n7n6); FID_S_ODE (n7n);

97 FEM3D 97 Generatng Connectvt of Matr MA_CO0 (/) FID_S_ODE (nn); FID_S_ODE (nn2); FID_S_ODE (nn3); FID_S_ODE (nn); FID_S_ODE (nn5); FID_S_ODE (nn6); FID_S_ODE (nn7); for(n=0;n<;n++){ =ILU[n]; for (=0;<;++){ COL[]=IALU[n][]; Sort IALU[][] n ascendng order b bubble sortng for less than 00 components. msor(colcol2); for(=;>0;--){ IALU[n][-]= COL[COL2[-]-];

98 FEM3D 9 MA_CO: CRS format #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MA_CO() { nt ; ndelu=(ki*)allocate_vector(seof(ki)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(ki*)allocate_vector(seof(ki)plu); for(=0;<;++){ for(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); C nde[ ] nde[0] 0 FORRA nde(0) 0 0 ILU[ ] nde( ) ILU( )

99 FEM3D 99 MA_CO: CRS format #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MA_CO() { nt ; ndelu=(ki*)allocate_vector(seof(ki)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(ki*)allocate_vector(seof(ki)plu); for(=0;<;++){ for(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); PLU=ndeLU[] Se of arra: temlu otal number of non-ero offdagonal blocs MA_CO

100 FEM3D 00 MA_CO: CRS format #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MA_CO() { nt ; ndelu=(ki*)allocate_vector(seof(ki)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(ki*)allocate_vector(seof(ki)plu); for(=0;<;++){ for(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); temlu store node ID startng from 0 MA_CO

101 FEM3D 0 MA_CO: CRS format #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MA_CO() { nt ; ndelu=(ki*)allocate_vector(seof(ki)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(ki*)allocate_vector(seof(ki)plu); for(=0;<;++){ for(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); ot requred an more MA_CO

102 FEM3D 02 Man Part /** program heat3d **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPU_CL(); etern vod IPU_GRID(); etern vod MA_CO0(); etern vod MA_CO(); etern vod MA_ASS_MAI(); etern vod MA_ASS_BC(); etern vod SOLVE(); etern vod OUPU_UCD(); nt man() { IPU_CL(); IPU_GRID(); MA_CO0(); MA_CO(); MA_ASS_MAI(); MA_ASS_BC() ; SOLVE(); OUPU_UCD() ;

103 FEM3D 03 MA_ASS_MAI: Overvew do pn= 2 Gaussan Quad. ponts n -drecton do jpn= 2 Gaussan Quad. ponts n -drecton do pn= 2 Gaussan Quad. Ponte n -drecton Defne Shape Functon at Gaussan Quad. Ponts (-ponts) Its dervatve on natural/local coordnate s also defned. enddo enddo enddo do cel= ICELO Loop for Element Jacoban and dervatve on global coordnate of shape functons at Gaussan Quad. Ponts are defned accordng to coordnates of nodes.(jacobi) do e= Local ode ID do je= Local ode ID Global ode ID: p jp Address of A pjp n temlu : j e do pn= 2 Gaussan Quad. ponts n -drecton do jpn= 2 Gaussan Quad. ponts n -drecton do pn= 2 Gaussan Quad. ponts n -drecton ntegraton on each element coeffcents of element matrces accumulaton to global matr enddo enddo enddo enddo enddo enddo e

104 FEM3D 0 #nclude <stdo.h> #nclude <math.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; etern vod JACOBI(); vod MA_ASS_MAI() { nt ; nt pjpp; nt pnjpnpn; nt cel; nt eje; nt SE; nt nn2n3nn5n6n7n; double SH; double QPQMEPEMPM; double XX2X3XX5X6X7X; double YY2Y3YY5Y6Y7Y; double ZZ2Z3ZZ5Z6Z7Z; double PXPYPZPXjPYjPZj; double COD0 QV0 QVC COEFj; double coef; KI nodlocal[]; MA_ASS_MAI (/6) AMA=(KREAL*) allocate_vector(seof(kreal)plu); B =(KREAL*) allocate_vector(seof(kreal) ); D =(KREAL*) allocate_vector(seof(kreal)); X =(KREAL*) allocate_vector(seof(kreal)); on-zero Off-Dagonal components (coef. matr) RHS vector Dagonal components (coef. matr) Unnowns for(=0;<plu;++) AMA[]=0.0; for(=0;< ;++) B[]=0.0; for(=0;< ;++) D[]=0.0; for(=0;< ;++) X[]=0.0; WEI[0]= e0; WEI[]= e0; POS[0]= e0; POS[]= e0;

105 FEM3D 05 MA_ASS_MAI (/6) #nclude <stdo.h> #nclude <math.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; etern vod JACOBI(); vod MA_ASS_MAI() { nt ; nt pjpp; nt pnjpnpn; nt cel; nt eje; nt SE; nt nn2n3nn5n6n7n; double SH; double QPQMEPEMPM; double XX2X3XX5X6X7X; double YY2Y3YY5Y6Y7Y; double ZZ2Z3ZZ5Z6Z7Z; double PXPYPZPXjPYjPZj; double COD0 QV0 QVC COEFj; double coef; KI nodlocal[]; AMA=(KREAL*) allocate_vector(seof(kreal)plu); B =(KREAL*) allocate_vector(seof(kreal) ); D =(KREAL*) allocate_vector(seof(kreal)); X =(KREAL*) allocate_vector(seof(kreal)); for(=0;<plu;++) AMA[]=0.0; for(=0;< ;++) B[]=0.0; for(=0;< ;++) D[]=0.0; for(=0;< ;++) X[]=0.0; WEI[0]= e0; WEI[]= e0; POS[0]= e0; POS[]= e0; POS: WEI: Quad. Pont Weghtng Factor

106 FEM3D 06 MA_ASS_MAI (2/6) /*** II. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b EA P - st-order dervatve of shape functon b ZE ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; P=.e0 + POS[p]; M=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * M; SHAPE[p][jp][p][]= Oth * QP * EM * M; SHAPE[p][jp][p][2]= Oth * QP * EP * M; SHAPE[p][jp][p][3]= Oth * QM * EP * M; SHAPE[p][jp][p][]= Oth * QM * EM * P; SHAPE[p][jp][p][5]= Oth * QP * EM * P; SHAPE[p][jp][p][6]= Oth * QP * EP * P; SHAPE[p][jp][p][7]= Oth * QM * EP * P;

107 FEM3D 07 MA_ASS_MAI (2/6) /*** II. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b EA P - st-order dervatve of shape functon b ZE ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; P=.e0 + POS[p]; M=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * M; SHAPE[p][jp][p][]= Oth * QP * EM * M; SHAPE[p][jp][p][2]= Oth * QP * EP * M; SHAPE[p][jp][p][3]= Oth * QM * EP * M; SHAPE[p][jp][p][]= Oth * QM * EM * P; SHAPE[p][jp][p][5]= Oth * QP * EM * P; SHAPE[p][jp][p][6]= Oth * QP * EP * P; SHAPE[p][jp][p][7]= Oth * QM * EP * P; QP EP P QM j EM j j M

108 FEM3D 0 MA_ASS_MAI (2/6) /*** II. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b EA P - st-order dervatve of shape functon b ZE ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; P=.e0 + POS[p]; M=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * M; SHAPE[p][jp][p][]= Oth * QP * EM * M; SHAPE[p][jp][p][2]= Oth * QP * EP * M; SHAPE[p][jp][p][3]= Oth * QM * EP * M; SHAPE[p][jp][p][]= Oth * QM * EM * P; SHAPE[p][jp][p][5]= Oth * QP * EM * P; SHAPE[p][jp][p][6]= Oth * QP * EP * P; SHAPE[p][jp][p][7]= Oth * QM * EP * P;

109 FEM3D 09 MA_ASS_MAI (2/6) /*** II. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b EA P - st-order dervatve of shape functon b ZE ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; P=.e0 + POS[p]; M=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * M; SHAPE[p][jp][p][]= Oth * QP * EM * M; SHAPE[p][jp][p][2]= Oth * QP * EP * M; SHAPE[p][jp][p][3]= Oth * QM * EP * M; SHAPE[p][jp][p][]= Oth * QM * EM * P; SHAPE[p][jp][p][5]= Oth * QP * EM * P; SHAPE[p][jp][p][6]= Oth * QP * EP * P; SHAPE[p][jp][p][7]= Oth * QM * EP * P; ) ( ) ( ) ( ) ( 3 2 ) ( ) ( ) ( ) ( 7 6 5

110 FEM3D 0 MA_ASS_MAI (3/6) PQ[jp][p][0]= - Oth * EM * M; PQ[jp][p][]= + Oth * EM * M; PQ[jp][p][2]= + Oth * EP * M; PQ[jp][p][3]= - Oth * EP * M; PQ[jp][p][]= - Oth * EM * P; PQ[jp][p][5]= + Oth * EM * P; PQ[jp][p][6]= + Oth * EP * P; PQ[jp][p][7]= - Oth * EP * P; PE[p][p][0]= - Oth * QM * M; PE[p][p][]= - Oth * QP * M; PE[p][p][2]= + Oth * QP * M; PE[p][p][3]= + Oth * QM * M; PE[p][p][]= - Oth * QM * P; PE[p][p][5]= - Oth * QP * P; PE[p][p][6]= + Oth * QP * P; PE[p][p][7]= + Oth * QM * P; P[p][jp][0]= - Oth * QM * EM; P[p][jp][]= - Oth * QP * EM; P[p][jp][2]= - Oth * QP * EP; P[p][jp][3]= - Oth * QM * EP; P[p][jp][]= + Oth * QM * EM; P[p][jp][5]= + Oth * QP * EM; P[p][jp][6]= + Oth * QP * EP; P[p][jp][7]= + Oth * QM * EP; for( cel=0;cel< ICELO;cel++){ COD0= COD; n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; PQ( j ) l ( j ) PE( ) l ( j ) P ( j) l ( j ) ( j ) 2 ( j ) 3 ( j ) 3 ( j ) ( j ) j j j Frst Order Dervatve of Shape Functons at j

111 FEM3D MA_ASS_MAI (3/6) PQ[jp][p][0]= - Oth * EM * M; PQ[jp][p][]= + Oth * EM * M; PQ[jp][p][2]= + Oth * EP * M; PQ[jp][p][3]= - Oth * EP * M; PQ[jp][p][]= - Oth * EM * P; PQ[jp][p][5]= + Oth * EM * P; PQ[jp][p][6]= + Oth * EP * P; PQ[jp][p][7]= - Oth * EP * P; PE[p][p][0]= - Oth * QM * M; PE[p][p][]= - Oth * QP * M; PE[p][p][2]= + Oth * QP * M; PE[p][p][3]= + Oth * QM * M; PE[p][p][]= - Oth * QM * P; PE[p][p][5]= - Oth * QP * P; PE[p][p][6]= + Oth * QP * P; PE[p][p][7]= + Oth * QM * P; P[p][jp][0]= - Oth * QM * EM; P[p][jp][]= - Oth * QP * EM; P[p][jp][2]= - Oth * QP * EP; P[p][jp][3]= - Oth * QM * EP; P[p][jp][]= + Oth * QM * EM; P[p][jp][5]= + Oth * QP * EM; P[p][jp][6]= + Oth * QP * EP; P[p][jp][7]= + Oth * QM * EP; for( cel=0;cel< ICELO;cel++){ COD0= COD; n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7];

112 FEM3D 2 MA_ASS_MAI (/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; ode ID (Global) QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; JACOBI(DEJ PQ PE P PX PY PZ X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z);

113 FEM3D 3 MA_ASS_MAI (/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; X-Coordnates of nodes Y-Coordnates of nodes QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; Z-Coordnates of nodes Coordnates: ode ID - JACOBI(DEJ PQ PE P PX PY PZ X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z);

114 FEM3D MA_ASS_MAI (/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; X-Coordnates of nodes Y-Coordnates of nodes Coordnates: ode ID QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; JACOBI(DEJ PQ PE P PX PY PZ X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z); Q QVOL C C Q 0 Heat Gen. Rate s a functon of locaton (cell center: c c )

115 FEM3D 5 MA_ASS_MAI (/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; Coordnates: ode ID QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; QVC C C JACOBI(DEJ PQ PE P PX PY PZ X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z); Q QVOL C C Q 0

116 FEM3D 6 MA_ASS_MAI (/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; JACOBI(DEJ PQ PE P PX PY PZ X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z);

117 FEM3D 7 JACOBI (/) #nclude <stdo.h> #nclude <math.h> #nclude "precson.h" #nclude "allocate.h" /*** *** JACOBI ***/ vod JACOBI( KREAL DEJ[2][2][2] KREAL PQ[2][2][]KREAL PE[2][2][]KREAL P[2][2][] KREAL PX[2][2][2][]KREAL PY[2][2][2][]KREAL PZ[2][2][2][] KREAL XKREAL X2KREAL X3KREAL XKREAL X5KREAL X6KREAL X7KREAL X KREAL YKREAL Y2KREAL Y3KREAL YKREAL Y5KREAL Y6KREAL Y7KREAL Y KREAL ZKREAL Z2KREAL Z3KREAL ZKREAL Z5KREAL Z6KREAL Z7KREAL Z) { /** calculates JACOBIA & IVERSE JACOBIA d/d d/d & d/d **/ nt pjpp; double dxdqdydqdzdqdxdedydedzdedxddyddzd; double coef; double aa2a3a2a22a23a3a32a33; for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ PX[p][jp][p][0]=0.0; PX[p][jp][p][]=0.0; PX[p][jp][p][2]=0.0; PX[p][jp][p][3]=0.0; PX[p][jp][p][]=0.0; PX[p][jp][p][5]=0.0; PX[p][jp][p][6]=0.0; PX[p][jp][p][7]=0.0; l l l Input l ~ Output l l l det J Values at each Gaussan Quad. Ponts: [p][jp][p] l l l

118 FEM3D Partal Dff. on atural Coord. (/) Accordng to formulae: ) ( ) ( ) ( can be easl derved accordng to defntons. are requred for computatons.

119 FEM3D 9 Partal Dff. on atural Coord. (2/) In matr form: J J J J J J J J J J J J : Jacob matr Jacoban

120 FEM3D 20 Partal Dff. on atural Coord. (3/) Components of Jacoban: J J J J J J J J J

121 FEM3D 2 JACOBI (2/) PY[p][jp][p][0]=0.0; PY[p][jp][p][]=0.0; PY[p][jp][p][2]=0.0; PY[p][jp][p][3]=0.0; PY[p][jp][p][]=0.0; PY[p][jp][p][5]=0.0; PY[p][jp][p][6]=0.0; PY[p][jp][p][7]=0.0; PZ[p][jp][p][0]=0.0; PZ[p][jp][p][]=0.0; PZ[p][jp][p][2]=0.0; PZ[p][jp][p][3]=0.0; PZ[p][jp][p][]=0.0; PZ[p][jp][p][5]=0.0; PZ[p][jp][p][6]=0.0; PZ[p][jp][p][7]=0.0; J J J J 2 3 J J J J J J /** **/ DEERMIA of the JACOBIA dxdq = PQ[jp][p][0]*X + PQ[jp][p][]*X2 + PQ[jp][p][2]*X3 + PQ[jp][p][3]*X + PQ[jp][p][]*X5 + PQ[jp][p][5]*X6 + PQ[jp][p][6]*X7 + PQ[jp][p][7]*X; dydq = PQ[jp][p][0]*Y + PQ[jp][p][]*Y2 + PQ[jp][p][2]*Y3 + PQ[jp][p][3]*Y + PQ[jp][p][]*Y5 + PQ[jp][p][5]*Y6 + PQ[jp][p][6]*Y7 + PQ[jp][p][7]*Y; dzdq = PQ[jp][p][0]*Z + PQ[jp][p][]*Z2 + PQ[jp][p][2]*Z3 + PQ[jp][p][3]*Z + PQ[jp][p][]*Z5 + PQ[jp][p][5]*Z6 + PQ[jp][p][6]*Z7 + PQ[jp][p][7]*Z; dxde = PE[p][p][0]*X + PE[p][p][]*X2 + PE[p][p][2]*X3 + PE[p][p][3]*X + PE[p][p][]*X5 + PE[p][p][5]*X6 + PE[p][p][6]*X7 + PE[p][p][7]*X; dxdq J dydq J dzdq J 2 3

122 FEM3D 22 /** IVERSE JACOBIA **/ JACOBI (3/) dyde = PE[p][p][0]*Y + PE[p][p][]*Y2 + PE[p][p][2]*Y3 + PE[p][p][3]*Y + PE[p][p][]*Y5 + PE[p][p][5]*Y6 + PE[p][p][6]*Y7 + PE[p][p][7]*Y; dzde = PE[p][p][0]*Z + PE[p][p][]*Z2 + PE[p][p][2]*Z3 + PE[p][p][3]*Z + PE[p][p][]*Z5 + PE[p][p][5]*Z6 + PE[p][p][6]*Z7 + PE[p][p][7]*Z; dxd = P[p][jp][0]*X + P[p][jp][]*X2 + P[p][jp][2]*X3 + P[p][jp][3]*X + P[p][jp][]*X5 + P[p][jp][5]*X6 + P[p][jp][6]*X7 + P[p][jp][7]*X; dyd = P[p][jp][0]*Y + P[p][jp][]*Y2 + P[p][jp][2]*Y3 + P[p][jp][3]*Y + P[p][jp][]*Y5 + P[p][jp][5]*Y6 + P[p][jp][6]*Y7 + P[p][jp][7]*Y; dzd = P[p][jp][0]*Z + P[p][jp][]*Z2 + P[p][jp][2]*Z3 + P[p][jp][3]*Z + P[p][jp][]*Z5 + P[p][jp][5]*Z6 + P[p][jp][6]*Z7 + P[p][jp][7]*Z; DEJ[p][jp][p]= dxdq*(dyde*dzd-dzde*dyd) + dydq*(dzde*dxd-dxde*dzd) + dzdq*(dxde*dyd-dyde*dxd); coef=.0 / DEJ[p][jp][p]; a= coef * ( dyde*dzd - dzde*dyd ); a2= coef * ( dzdq*dyd - dydq*dzd ); a3= coef * ( dydq*dzde - dzdq*dyde ); J J J J 2 3 J J J J J J a2= coef * ( dzde*dxd - dxde*dzd ); a22= coef * ( dxdq*dzd - dzdq*dxd ); a23= coef * ( dzdq*dxde - dxdq*dzde ); a3= coef * ( dxde*dyd - dyde*dxd ); a32= coef * ( dydq*dxd - dxdq*dyd ); a33= coef * ( dxdq*dyde - dydq*dxde ); DEJ[p][jp][p]=fabs(DEJ[p][jp][p]);

123 FEM3D 23 Partal Dff. on atural Coord. (/) Partal dfferentaton on global coordnate sstem s ntroduced as follows (wth nverse of Jacoban matr (3 3)) J

124 FEM3D 2 /** IVERSE JACOBIA JACOBI (3/) dyde = PE[p][p][0]*Y + PE[p][p][]*Y2 + PE[p][p][2]*Y3 + PE[p][p][3]*Y + PE[p][p][]*Y5 + PE[p][p][5]*Y6 + PE[p][p][6]*Y7 + PE[p][p][7]*Y; dzde = PE[p][p][0]*Z + PE[p][p][]*Z2 + PE[p][p][2]*Z3 + PE[p][p][3]*Z + PE[p][p][]*Z5 + PE[p][p][5]*Z6 + PE[p][p][6]*Z7 + PE[p][p][7]*Z; dxd = P[p][jp][0]*X + P[p][jp][]*X2 + P[p][jp][2]*X3 + P[p][jp][3]*X + P[p][jp][]*X5 + P[p][jp][5]*X6 + P[p][jp][6]*X7 + P[p][jp][7]*X; dyd = P[p][jp][0]*Y + P[p][jp][]*Y2 + P[p][jp][2]*Y3 + P[p][jp][3]*Y + P[p][jp][]*Y5 + P[p][jp][5]*Y6 + P[p][jp][6]*Y7 + P[p][jp][7]*Y; dzd = P[p][jp][0]*Z + P[p][jp][]*Z2 + P[p][jp][2]*Z3 + P[p][jp][3]*Z + P[p][jp][]*Z5 + P[p][jp][5]*Z6 + P[p][jp][6]*Z7 + P[p][jp][7]*Z; DEJ[p][jp][p]= dxdq*(dyde*dzd-dzde*dyd) + dydq*(dzde*dxd-dxde*dzd) + dzdq*(dxde*dyd-dyde*dxd); **/ coef=.0 / DEJ[p][jp][p]; a= coef * ( dyde*dzd - dzde*dyd ); a2= coef * ( dzdq*dyd - dydq*dzd ); a3= coef * ( dydq*dzde - dzdq*dyde ); a2= coef * ( dzde*dxd - dxde*dzd ); a22= coef * ( dxdq*dzd - dzdq*dxd ); a23= coef * ( dzdq*dxde - dxdq*dzde ); a3= coef * ( dxde*dyd - dyde*dxd ); a32= coef * ( dydq*dxd - dxdq*dyd ); a33= coef * ( dxdq*dyde - dydq*dxde ); J a a a 2 3 a a a a a a DEJ[p][jp][p]=fabs(DEJ[p][jp][p]);

125 FEM3D 25 JACOBI (/) /** set the d/dx d/dy & d/dz components **/ PX[p][jp][p][0]=a*PQ[jp][p][0]+a2*PE[p][p][0]+a3*P[p][jp][0]; PX[p][jp][p][]=a*PQ[jp][p][]+a2*PE[p][p][]+a3*P[p][jp][]; PX[p][jp][p][2]=a*PQ[jp][p][2]+a2*PE[p][p][2]+a3*P[p][jp][2]; PX[p][jp][p][3]=a*PQ[jp][p][3]+a2*PE[p][p][3]+a3*P[p][jp][3]; PX[p][jp][p][]=a*PQ[jp][p][]+a2*PE[p][p][]+a3*P[p][jp][]; PX[p][jp][p][5]=a*PQ[jp][p][5]+a2*PE[p][p][5]+a3*P[p][jp][5]; PX[p][jp][p][6]=a*PQ[jp][p][6]+a2*PE[p][p][6]+a3*P[p][jp][6]; PX[p][jp][p][7]=a*PQ[jp][p][7]+a2*PE[p][p][7]+a3*P[p][jp][7]; PY[p][jp][p][0]=a2*PQ[jp][p][0]+a22*PE[p][p][0]+a23*P[p][jp][0]; PY[p][jp][p][]=a2*PQ[jp][p][]+a22*PE[p][p][]+a23*P[p][jp][]; PY[p][jp][p][2]=a2*PQ[jp][p][2]+a22*PE[p][p][2]+a23*P[p][jp][2]; PY[p][jp][p][3]=a2*PQ[jp][p][3]+a22*PE[p][p][3]+a23*P[p][jp][3]; PY[p][jp][p][]=a2*PQ[jp][p][]+a22*PE[p][p][]+a23*P[p][jp][]; PY[p][jp][p][5]=a2*PQ[jp][p][5]+a22*PE[p][p][5]+a23*P[p][jp][5]; PY[p][jp][p][6]=a2*PQ[jp][p][6]+a22*PE[p][p][6]+a23*P[p][jp][6]; PY[p][jp][p][7]=a2*PQ[jp][p][7]+a22*PE[p][p][7]+a23*P[p][jp][7]; PZ[p][jp][p][0]=a3*PQ[jp][p][0]+a32*PE[p][p][0]+a33*P[p][jp][0]; PZ[p][jp][p][]=a3*PQ[jp][p][]+a32*PE[p][p][]+a33*P[p][jp][]; PZ[p][jp][p][2]=a3*PQ[jp][p][2]+a32*PE[p][p][2]+a33*P[p][jp][2]; PZ[p][jp][p][3]=a3*PQ[jp][p][3]+a32*PE[p][p][3]+a33*P[p][jp][3]; PZ[p][jp][p][]=a3*PQ[jp][p][]+a32*PE[p][p][]+a33*P[p][jp][]; PZ[p][jp][p][5]=a3*PQ[jp][p][5]+a32*PE[p][p][5]+a33*P[p][jp][5]; PZ[p][jp][p][6]=a3*PQ[jp][p][6]+a32*PE[p][p][6]+a33*P[p][jp][6]; PZ[p][jp][p][7]=a3*PQ[jp][p][7]+a32*PE[p][p][7]+a33*P[p][jp][7]; a a a a a a a a a

126 FEM3D 26 MA_ASS_MAI (5/6) /** COSRUC the GLOBAL MARIX **/ for(e=0;e<;e++){ p=nodlocal[e]; for(je=0;je<;je++){ jp=nodlocal[je]; =-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( =S;<E;++){ f( temlu[] == jp- ){ =; brea; on-zero Off-Dagonal Bloc n Global Mat A p jp : address n temlu p= nodlocal[e] jp= nodlocal[je] ode ID (pjp) startng from 2

127 FEM3D 27 Element Matr: j j j

128 FEM3D 2 MA_ASS_MAI (5/6) /** COSRUC the GLOBAL MARIX **/ for(e=0;e<;e++){ p=nodlocal[e]; for(je=0;je<;je++){ jp=nodlocal[je]; =-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( =S;<E;++){ f( temlu[] == jp- ){ =; brea; Element Matr( e ~j e ): Local ID Global Matr ( p ~j p ): Global ID :address n temlu startng from 0 : startng from 0 pjp: startng from

129 FEM3D 29 MA_ASS_MAI (6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(dej[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMA[]+= COEFj; j j j det J ddd

130 FEM3D 30 MA_ASS_MAI (6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(dej[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMA[]+= COEFj; d d d J j j j det j j M j L f W W W d d d f I ) ( ) (