3D-FEM n C: Stead State Heat Conducton Kengo aajma Informaton echnolog Center Programmng for Parallel Computng (66-2057) Semnar on Advanced Computng (66-009)
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
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
FEM3D Fnte-Element Procedures Governng Equatons Galern Method: Wea Form Element-b-Element Integraton Element Matr Global Matr Boundar Condtons Lnear Solver
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
FEM3D 6 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program
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.
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
FEM3D 9 Isoparametrc Element (/3) Each element s mapped to square element [±±] on natural/local coordnate () 3 2 + 3 - + - 2 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.
0 Isoparametrc Element (2/3) 2 3 Coordnate of each node: ( ) ( 2 2 ) ( 3 3 ) ( ) emperature at each node: 2 3 3 - + + - 2 FEM3D ) ( ) ( ) (
FEM3D Isoparametrc Element (3/3) + - + - Hgher-order shape functon can handle curved lnes/surfaces. atural coordnate sstem Sub-Parametrc Super-Parametrc
2 Shape Fn s on 2D atural Coord. (/3) Polnomal shape functons on squares of natural coordnate: Coeffcents are calculated as follows: - + + - 3 2 FEM3D 3 2 3 2 3 3 2 2 3 2 3 2
Shape Fn s on 2D atural Coord. (2/3) - + + - 3 2 FEM3D 3 3 2 3 2 3 2 3 2 3 2 3 2 3 2
Shape Fn s on 2D atural Coord. (2/3) - + + - 3 2 FEM3D 3 2 3 2 3 2 3 2 3 2 3 2 3 2 2 3
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 - + + - 3 2 FEM3D 3 3 2 2
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( 六面体要素 )
7 Shape Fn s: 3D atural/local Coord. ) ( ) ( ) ( ) ( 3 2 ) ( ) ( ) ( ) ( 7 6 5 FEM3D 2 3 5 6 7 ) ( ) ( ) ( ) (
FEM3D Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program
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: 0 2 2 2 2 2 2 dv Q V 0 2 2 2 2 2 2 Q FEM3D 9 Dstrbuton of temperature n each element (matr form) : emperature at each node
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 2 2 2 2 2 2 dv dv V V FEM3D 20
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
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
FEM3D 23 Element Matr: j j j 7 5 6 3 2
FEM3D 2 Element Matr: j j j 7 5 6 3 2 dv ( e) V V dv V dv j V j j j dv
FEM3D 25 et Stage: Integraton
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
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. 0023 2 3-0 + S h 2 2 / 2 2 0 f d W f 0.5773502692 0. 997 f h d
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 2 3 0.00 0.00 w 2.00 0.577350 w w / 9 0.77597 w.00 5/ 9 =- =0 =-0.577350 =+0.577350 =+
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 hs confguraton s wdel used. In 2D problem ntegraton s done usng values of f at quad. ponts. Gaussan Quadrature ガウスの積分公式
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.57735 0.57735).0.0 f ( 0.57735 0.57735).0.0 f ( 0.57735 0.57735).0.0 f ( 0.57735 0.57735) W j f j
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
Partal Dff. on atural Coord. (/) Accordng to formulae: ) ( ) ( ) ( can be easl derved accordng to defntons. are requred for computatons. FEM3D 33
Partal Dff. on atural Coord. (2/) In matr form: J J : Jacob matr Jacoban FEM3D 3
Partal Dff. on atural Coord. (3/) Components of Jacoban: J J J J J J J J J 33 32 3 23 22 2 3 2 FEM3D 35
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
Integraton on Element V j j j V j j j j dv dv 2 3 5 6 7 j j FEM3D 37
Integraton on atural Coord. d d d J ddd dv j j j j j j V j j j det 2 3 5 6 7 j j FEM3D 3
Gaussan Quadrature 2 3 5 6 7 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
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
Accumulaton: Local -> Global Matrces 2 5 6 3 2 2 3 2 2 3 () () 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 6 5 3 2 6 5 3 2 6 5 3 2 { ]{ [ 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
FEM3D 2 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Eercse Runnng the Code Data Structure Overvew of the Program
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:(3.0 5.0) :(2.0.0) I V dv det J dd
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()= -0.5773502692d0 POI(2)= +0.5773502692d0 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
FEM3D 5 ps (2/2) J J det ) ( ) ( ) ( ) ( 3 2
FEM3D 6 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program
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
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
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
FEM3D 50 Mesh Generaton Z >$ cd <$fem>/fem3d/run >$./mgcube =0@Z= ma X Y Z umber of Elem s n each drecton 202020 eample Z >$ ls cube.0 confrmaton cube.0 X Y X Y
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
FEM3D 52 >$ cd <$E-OP>/fem3d/run >$./sol Runnng >$ ls test.np Confrmaton test.np
FEM3D 53 ParaVew http://www.paravew.org/ Openng fles Dsplang fgures Savng mage fles http://nl.cc.u-too.ac.jp/class/howtouseparavewe.pdf http://nl.cc.u-too.ac.jp/class/howtouseparavewj.pdf
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
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.
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)
FEM3D 57 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program
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
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
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
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
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
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)
FEM3D 6 Eample of cube.0 (X=Y=Z=) ode Element-pe 25 =5*5*5 0.000000E+00 0.000000E+00 0.000000E+00 2.000000E+00 0.000000E+00 0.000000E+00 3 2.000000E+00 0.000000E+00 0.000000E+00 3.000000E+00 0.000000E+00 0.000000E+00 5.000000E+00 0.000000E+00 0.000000E+00 6 0.000000E+00.000000E+00 0.000000E+00 7.000000E+00.000000E+00 0.000000E+00 2.000000E+00.000000E+00 0.000000E+00 9 3.000000E+00.000000E+00 0.000000E+00 ( 途中省略 ) 2 0.000000E+00.000000E+00.000000E+00 22.000000E+00.000000E+00.000000E+00 23 2.000000E+00.000000E+00.000000E+00 2 3.000000E+00.000000E+00.000000E+00 25.000000E+00.000000E+00.000000E+00 6 =** 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 move
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 (=) 7 5 6 3 2
FEM3D 66 Eample of cube.0 (X=Y=Z=) Element Connectvt 2 7 6 26 27 32 3 2 2 3 7 27 2 33 32 3 3 9 2 29 3 33 5 0 9 29 30 35 3 5 6 7 2 3 32 37 36 6 7 3 2 32 33 3 37 7 9 3 33 3 39 3 9 0 5 3 35 0 39 9 2 7 6 36 37 2 0 2 3 7 37 3 3 2 3 9 3 39 3 2 5 20 9 39 0 5 3 6 7 22 2 2 7 6 (...) 2 62 63 6 67 7 93 92 3 63 6 69 6 9 9 93 6 65 70 69 9 90 95 9 5 66 67 72 7 9 92 97 96 6 67 6 73 72 92 93 9 97 7 6 69 7 73 93 9 99 9 69 70 75 7 9 95 00 99 9 76 77 2 0 02 07 06 50 77 7 3 2 02 03 0 07 5 7 79 3 03 0 09 0 52 79 0 5 0 05 0 09 53 2 7 6 06 07 2 5 2 3 7 07 0 3 2 55 3 9 0 09 3 56 5 90 9 09 0 5 57 6 7 92 9 2 7 6 5 7 93 92 2 3 7 59 9 9 93 3 9 60 9 90 95 9 5 20 9 6 9 92 97 96 6 7 22 2 62 92 93 9 97 7 23 22 63 93 9 99 9 9 2 23 6 9 95 00 99 9 20 25 2
FEM3D 67 X=Y=Z= XP=YP=ZP=5 ICELO= 6 IODO= 25 IBODO= 25 = 2 9 2 5 6 25 20 7 9 6 0 =2 9 37 36 2 2 3 7 9 20 33 3 35 36 32 50 5 32 33 3 3 35 =3 7 69 62 6 37 75 70 57 5 59 56 60 2 3 5 26 27 2 29 30 5 52 53 5 55 =5 2 9 6 25 20 99 9 6 00 95 = 53 2 07 0 09 06 0 9 50 5 52 0 02 03 0 05 6 53 7 2 3 5 9 50 5 52 76 77 7 79 0
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)
FEM3D 69 Eample of cube.0 (X=Y=Z=) Group Info. Xmn Ymn Zmn Zma 25 50 75 00 6 6 2 26 3 36 6 5 56 6 66 7 76 6 9 96 0 06 6 2 2 3 5 26 27 2 29 30 5 52 53 5 55 76 77 7 79 0 0 02 03 0 05 2 3 5 6 7 9 0 2 3 5 6 7 9 20 2 22 23 2 25 0 02 03 0 05 06 07 0 09 0 2 3 5 6 7 9 20 2 22 23 2 25 (not n use after ths pont)
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
FEM3D 7 Formulaton of 3D Element 3D Heat Equatons Galern Method Element Matrces Runnng the Code Data Structure Overvew of the Program
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
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
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() ;
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
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
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
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
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]);
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
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);
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
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
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
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
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
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]
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() { 3 5 6 7 9 9 0 2 5 6 5 6 7 2 3 2 3 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
FEM3D 9 MA_CO0: Overvew do cel= ICELO generate ILU IALU accordng to nodes of he. elements (FID_ODE) enddo 7 5 6 3 3 5 6 7 9 9 0 2 2 5 6 5 6 7 2 3 2 3
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
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;
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); 7 5 6 3 2 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);
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)
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; 3 5 6 7 9 9 0 2 5 6 5 6 7 2 3 2 3
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. 3 5 6 7 9 9 0 2 5 6 5 6 7 2 3 2 3
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); 7 5 6 3 2 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);
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[-]-];
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( )
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
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
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
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() ;
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
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]=.0000000000e0; WEI[]=.0000000000e0; POS[0]= -0.5773502692e0; POS[]= 0.5773502692e0;
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]=.0000000000e0; WEI[]=.0000000000e0; POS[0]= -0.5773502692e0; POS[]= 0.5773502692e0; POS: WEI: Quad. Pont Weghtng Factor
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;
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
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; 7 5 6 2 3
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
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
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]; 7 5 6 3 2
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) 7 5 6 2 3 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);
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 7 5 6 2 3 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);
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 2 5 6 Coordnates: ode ID - 7 3 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 )
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-][]; 2 5 6 Coordnates: ode ID - 7 3 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
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);
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
FEM3D Partal Dff. on atural Coord. (/) Accordng to formulae: ) ( ) ( ) ( can be easl derved accordng to defntons. are requred for computatons.
FEM3D 9 Partal Dff. on atural Coord. (2/) In matr form: J 33 32 3 23 22 2 3 2 J J J J J J J J J J J : Jacob matr Jacoban
FEM3D 20 Partal Dff. on atural Coord. (3/) Components of Jacoban: J J J J J J J J J 33 32 3 23 22 2 3 2
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 2 22 32 J J J 3 23 33 /** **/ 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
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 2 22 32 J J J 3 23 33 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]);
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
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 2 22 32 a a a 3 23 33 DEJ[p][jp][p]=fabs(DEJ[p][jp][p]);
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 33 32 3 23 22 2 3 2
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] 7 5 6 3 ode ID (pjp) startng from 2
FEM3D 27 Element Matr: j j j 7 5 6 3 2
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 7 5 6 3 2
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
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 ) ( ) (