David Lambert/isoparametric shape functions
Jump to navigation
Jump to search
Note 'how to use the interpolants iXXX' Rows of y are the coordinates (one point per line) in the element system (xi, eta, nu). These are the locations to interpolate the field. Rank of y is 2. The number of columns of y is the dimensionality and should match the first digit of interplant name. Anticipated that each row of y is a point on [_1, 1]. Rows of x are the known field values at the nodes, in, of course, the correct order. Beware, this node numbering scheme is unusual. The tally of columns of x is the tally of nodes in the fundamental element and must match the third digit of the interpolant name. These are the real world field values. On interpolating coordinates, note that x, y, and z are separate fields. Example: Interpolate two fields in one dimension. They are linear and the field values at the _1 end of the element are _4 and 2, and at the 1 end of the element the values are _8 and 3. The element is i122---1 dimensional, 2 "corners", and 2 nodes. Now we want to interpolate the field at five uniformly spaced positions. [ LINEAR_ELEMENT_COORDINATES =: ,.(%~i:)2x _1 _1r2 0 1r2 1 NB. interpolate 2 fields at these 5 points (_4 _8,:2 3) i122 LINEAR_ELEMENT_COORDINATES _4 2 _5 2.25 _6 2.5 _7 2.75 _8 3 To create a graded mesh use a double interpolation. Start with a uniform grid y argument on [_1, 1]. This is easy to compute. Use the x argument to transform onto the interval [-1, 1] with grading of a higher order interpolant. Verb prepare38z is especially nice for this. Use a second interpolation to map the coordinate block to the real world position. A "corners only" transfom might be useful. (X,Y,:Z) i388 0.2 Prepare38z node_grid (NX, NY, NZ) ) Note 'FEA' Here we develop a general method to generate isoparametric interpolants. Example for 4 node quadrilateral element. The interpolant is the dot product of the four shape function values evaluated at some position within the element against the known nodal values. The sum of four linear functions of two variables (xi, eta) is 1 at each of four nodes. Let the base element have nodal coordinates (xi, eta) of +/-1. f2 f3 (1,1) +---------------+ | f(x,y) | | (x,y) | | * (0,0) | | * | | | | | |(_1,_1) | +---------------+ f0 f1 determine n0(xi,eta), ..., n3(xi,eta). n0(-1,-1) = 1, n0(all other corners) is 0. n1( 1,-1) = 1, n1(all other corners) is 0. ... Choose a shape function: Use shape functions C0 + C1*xi + C2*eta + C3*xi*eta . Given (xi,eta) as the vector y form a vector of the coefficients of the constants (1, xi, eta, and their product) This method extends to higher order interpolants having more nodes or to other dimensions. ) example=: 3 :0 require'viewmat' GRID =: |.,~"0/~(%~i:)100 SADDLE =: 1.0 2.0 2.2 0.7 i244 GRID viewmat SADDLE (<./ , >./) , SADDLE ) Note 'Some elemental information' Node order 1D: 0 2 1 0 1 2 3 cubic 2D: 2 7 3 5 8 6 Node 8 at origin, Node 3 at (1,1) 0 4 1 c d e f bicubic (f) 8 9 a b 4 5 6 7 0 1 2 3 Names for interpolants, shape functions and constants: n means shape function n1 all linear combinations n2 all quadratic combinations n2Face n2 without center n3 all cubic combinations n38k 3D, corners + mid-edge 'k' = Num_j_,26}.Alpha_j_ i249: 2 dimensions, 4 corners (quadrilateral), 9 nodes, or C244: C constants for 2 dimensions, 4 corners (quadrilateral), 4 nodes 3D At z = _1 z = 0 z = 1 2 b 3 e o f 6 j 7 9 k a m q n h p i 0 8 1 c l d 4 g 5 terms included in shape funcions n38k includes corners and mid-edge nodes |1|x|y|z|xx|xy|xz|yy|yz|zz|xxy|xyy|xyz|xxz|xzz|yyz|yzz|xxyz|xyyz|xyzz| 1 1 x ] y z xy */"1@:((2 combinations 3)&{) xz yz xx *: yy zz xyz (1&, * */) xxyz xyyz xyzz xxy 0 0 1 *:@:(0&{) * 1 2&{ xxz 0 0 2 yyx 1 1 0 *:@:(1&{) * 0 2&{ yyz 1 1 2 zzx 2 2 0 *:@:(2&{) * 0 1&{ zzy 2 2 1 The shape functions for 27 node brick are given by (;:'1 z zz'),&.>/(;:'1 y yy'),&.>/(;:'1 x xx') ;L:1{(;:'1 z zz');(;:'1 y yy');<(;:'1 x xx') which leads directly to the product tables. ) combinations=: 4 : 0 if. x e. 0 1 do. z=.<((x!y),x)$ i.y else. t=. |.(<.@-:)^:(i.<. 2^.x)x z=.({.t) ([:(,.&.><@;\.)/ >:@-~[\i.@]) ({.t)+y-x for_j. 2[\t do. z=.([ ;@:(<"1@[ (,"1 ({.j)+])&.> ])&.> <@;\.({&.><)~ (1+({.j)-~{:"1)&.>) z if. 2|{:j do. z=.(i.1+y-x)(,.>:)&.> <@;\.z end. end. end. ;z ) mp=: +/ .*~~ NB. ($: |:) : (+/ .*) NB. A Atranspose : matrix product A B identity=: =@:i. NB. generate identity matrix NB. master nodes N1=: 4 A.-.3x#.inv i.3^1 N2=: 336330 A.-.3x#.inv i.3^2 NB. 336330 (-: A.) 8 2 6 0 5 7 1 3 4 N3=: 267337661061030402017459663x A.<:3x#.inv i.3^3 NB. 267337661061030402017459663x (-: A.) 0 18 6 24 2 20 8 26 9 3 21 15 1 19 7 25 11 5 23 17 12 10 4 22 16 14 13 NB. form of shape functions Shape=: 1 : '[: , [: *// ^/&(i.m)' n1 =: 2 Shape NB. all linear combinations n2 =: 3 Shape NB. all quadratic combinations n3 =: 4 Shape NB. all cubic combinations n2Face=: }:@:n2 NB. exclude center n38k =: 1 , ] , */"1@:((2 combinations 3)&{) , *: , (1&, * */) , ,@:(*:@:|. (*"0 1) (2 combinations 3)&{) NB. corners and mid-edge nodes NB. constants Const=: ("1)([:`(x:inv)`(identity@:#)`%.`)(`:6) C122=: n1 Const 2{.N1 C123=: n2 Const 3{.N1 C124=: n3 Const ,.3x%~i:3j3 NB. nodes simply in a row. C244=: n1 Const 4{.N2 NB. serendipity C248=: n2Face Const 8{.N2 NB. serendipity C249=: n2 Const 9{.N2 NB. non-serendipity C24g=: n3 Const <:2r3*(|."1)4#.inv i.4^2 NB. bicubic, nodes simply layed out i._4 4 C388=: n1 Const 8{.N3 C38k=: n38k Const 36bk{.N3 C38q=: n2Face Const 36bq{.N3 NB. 'q' = (8 corners + 12 edges + 6 faces) { Num_j_ , 26 }. Alpha_j_ C38r=: n2 Const 36br{.N3 NB. 'r' = (8 corners + 12 edges + 6 faces + 1 center) { Num_j_ , 26 }. Alpha_j_ C38z=: n3 Const <:2r3*(|."1)4#.inv i.4^3 NB. tricubic, nodes simply layed out i._4 4 NB. z stands for 64 Note'interpolants' rows of y are the coordinates in the element system (xi, eta, nu) where to interpolate the field. y''s rank is 2. The number of columns is the dimensionality. (first digit of interplant name) rows of y become the transformed output. rows of x are the known field values at the nodes. The tally of columns of x is the tally of nodes. (third digit of the interpolant name) Real world coordinates. Suppose the element has 2 nodes linear. We need to evaluate at 5 sites. xi, the sites at which to evaluate the field, shall be ,.(%~i:)2 and x is _4 at one end, _8 at the other end. while y runs from 2 to 3. (_4 _8,:2 3)i122 ,.(%~i:)2x _4 2 _5 2.25 _6 2.5 _7 2.75 _8 3 [N=:(prepare38z 1r5)i38z grid 18 2 2 NB. generate nodes in elemental system ) Weights=: 2 :'(m mp~ v)"_ 1' Interp=: 2 :'(mp (m mp~ v))"_ 1 f.' i122=: C122 Interp n1 i123=: C123 Interp n2 i124=: C124 Interp n3 i244=: C244 Interp n1 NB. corners only 4 nodes i248=: C248 Interp n2Face NB. corners and edges 8 nodes i249=: C249 Interp n2 NB. corners, edges, and face 9 nodes i24g=: C24g Interp n3 NB. bicubic 16 nodes i388=: C388 Interp n1 NB. verified corners 8 nodes i38k=: C38k Interp n38k NB. verified corners, edges 20 nodes i38q=: C38q Interp n2Face NB. verified corners, edges, faces 26 nodes i38r=: C38r Interp n2 NB. tested corners, edges, faces, center 27 nodes i38z=: C38z Interp n3 NB. verified tricubic 64 nodes NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. NB. mesh generation Monad=: : [: Dyad =: [: : odometer=: #: i.@:(*/) NB. node_grid generates a block of nodes spanning the element coordinate system [_1, 1] node_grid=: ([: <: odometer@:, (*"1) 2&%@:<:)&.:(|."1)@:(2&>.) Monad NB. node_grid Nx [Ny [Nz]] element_grid=: (([: - 2 ^ #) ]\ [: , (1 ,: # $ 2:) ];._3 i.@:|.)@:(2&>.) Monad NB. mesh NX NY ... creates a mesh entity ELEMENTS ; NODES on [_1, 1]... mesh=: (element_grid@:] ; node_grid) Monad Note 'prepare verbs, and adverbs' The prepare verbs generate x arguments for the interpolants that map onto _1 1 . As input are the locations of the midedge nodes in the elemental system on (_1, 1) maintaining rectilinearity, were that input as expected. ) NODES38R=: ".;._2 [ noun define _1 _1 _1 NB. 0 corner 1 _1 _1 NB. 1 corner _1 1 _1 NB. 2 corner 1 1 _1 NB. 3 corner _1 _1 1 NB. 4 corner 1 _1 1 NB. 5 corner _1 1 1 NB. 6 corner 1 1 1 NB. 7 corner 0 _1 _1 NB. 8 bottom mid-edge _1 0 _1 NB. 9 bottom mid-edge 1 0 _1 NB. 10 bottom mid-edge 0 1 _1 NB. 11 bottom mid-edge _1 _1 0 NB. 12 middle mid-edge 1 _1 0 NB. 13 middle mid-edge _1 1 0 NB. 14 middle mid-edge 1 1 0 NB. 15 middle mid-edge 0 _1 1 NB. 16 top mid-edge _1 0 1 NB. 17 top mid-edge 1 0 1 NB. 18 top mid-edge 0 1 1 NB. 19 top mid-edge 0 0 _1 NB. 20 face center 0 _1 0 NB. 21 face center _1 0 0 NB. 22 face center 1 0 0 NB. 23 face center 0 1 0 NB. 24 face center 0 0 1 NB. 25 face center 0 0 0 NB. 26 body center ) prepare38r=: 3 :0 'ux uy uz'=. 3 {. 0 ,~ y AN=. Num_j_,26}.Alpha_j_ ix=. AN i.'8kblqogpj' iy=. AN i.'9kamqnhpi' iz=. AN i.'eofmqncld' result=. NODES38R result=. ( ux) [`([: < ix ; 0:)`]} result result=. ( uy) [`([: < iy ; 1:)`]} result result=. ( uz) [`([: < iz ; 2:)`]} result |: result ) prepare38z=: 3 :0 NB. u[xyz] are the amount to stretch the interior nodes. NB. greater than one third opens the mesh interior, NB. less than one third squeezes the interior. NB. values not far from one third will cause overlapping mesh. NB. prepare38z generates the "x" argument that redistributes nodes NB. within the bounding element box for i38z . 'ux uy uz'=. 3 {. ,&(3 # 1r3) y nx=.1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 ny=.4 5 6 7 20 21 22 23 36 37 38 39 52 53 54 55 nz=. 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 px=. nx + 4 ^ 0 py=. ny + 4 ^ 1 pz=. nz + 4 ^ 2 result=. node_grid 4 4 4 result=. (-ux) [`([: < nx ; 0:)`]} result result=. ( ux) [`([: < px ; 0:)`]} result result=. (-uy) [`([: < ny ; 1:)`]} result result=. ( uy) [`([: < py ; 1:)`]} result result=. (-uz) [`([: < nz ; 2:)`]} result result=. ( uz) [`([: < pz ; 2:)`]} result |: result ) prepare123=: 0 1 4 {"1 [: }: prepare248 prepare124=: 1 4 {. prepare24g prepare248=: 0 1 2 3 8 9 10 11 {"1 [: }: prepare38r prepare249=: 0 1 2 3 8 9 10 11 20 {"1 [: }: prepare38r prepare24g=: (2 16 {. prepare38z) Monad prepare38k=: _ 20 {. prepare38r prepare38q=: _ 26 {. prepare38r Prepare123=: adverb def 'prepare123 m' Prepare124=: adverb def 'prepare124 m' Prepare248=: adverb def 'prepare248 m' Prepare249=: adverb def 'prepare249 m' Prepare24g=: adverb def 'prepare24g m' Prepare38k=: adverb def 'prepare38k m' Prepare38q=: adverb def 'prepare38q m' Prepare38r=: adverb def 'prepare38r m' Prepare38z=: adverb def 'prepare38z m' NB. mm=: 1 :'m % 1000' NB. um=: mm mm NB. ((0 6 0 6 mm,:0 0 250 250 um) ; 0.3) i244`i248`prepare248 Mesh 61 6 NB. x m Mesh y Mesh=: adverb define : prepare_grade=. (2{m)`:6 apply_grade=. (1{m)`:6 position=. (0{m)`:6 NODE_TALLIES=. y GRADATION=. 1 {:: x FINAL_POSITION=. 0 {:: x 'E N'=. mesh NODE_TALLIES N=. (prepare_grade GRADATION) apply_grade N E ; FINAL_POSITION position N ) spherical_octant=: 3 :0 NB. y controls node tally and should have length 0, 1, 3 or 4: NB. 0{y longitude 45 through 90 degrees NB. 1{y longitude 0 through 45 degrees NB. 2{y determines spacing latitude 0 to 45 NB. 3{y determines the outer shell through-thickness node tally (the outer shell delta radius) NB. the spherical surface has a dimple, certainly due to the triquadratic approximation. NB. after merging the mesh, move the nodes with distance from origin > 1 or a bit less, NB. and move them onto the surface. NB. roughly, nodes=.. ((* %@:mp)nodes) indexes} nodes if. 3 = # y do. y=. 2 1 1 # y elseif. (# y) -.@:e. 1 4 do. y=. 8 end. 'n0 n1 n2 n3'=. y NB. to do: are the elements right handed? 0 1 3 2 4 5 7 6 NB. merge the mesh NB. check with ANSYS NB. return the mesh instead of filing it. NB. mesh a spherical octant, radius 1. NB. warning! overwrites file /tmp/a and more round=. 1&$: :([ * [: <. 1r2 + *inv) NB. round y to nearest x create=. (1!:2 boxopen)~ append=. (1!:3 boxopen)~ d3edges=. [: (,:"2) 1&({::) {~ [: ~. [: /:~"1 [: ,/ (12 2$0 1 0 2 1 3 2 3 0 4 1 5 2 6 3 7 4 5 4 6 5 7 6 7) ({"_ 1) 0&({::) literate=. LF ,~ (([ }. (] ([: ,/ ,.&LF)"2)^:(-@:[))~ (0 <. -.@:(#@:$)))@:": NB. as a dyad x is x of ": algebraic_negative=. (=&'_')`(,:&'-')} comma=. #!.','~ (1 j. 0 ,~ 2 </\' '&=) ansys_elements=. [: literate [: ('E,' , comma)"1 [: ": 0 1 3 2 4 5 7 6 {"1 >: ansys_nodes=. [: algebraic_negative@:literate [: ('N,,' , comma)"1 ": export_ansys=. ansys_nodes@:(1&{::) , LF , ansys_elements@:(0&{::) NB. spherical to Cartesian. NB. The spherical system is as the earth system rather than that of a physics text. NB. That is, latitude 0 is the equitorial great circle. TAU=. 2p1 NB. tauday.com deg=. *&(TAU % 360) Deg=. 1 :'deg m' NB. this pretty adverb is not so useful here sin=. 1&o. cos=. 2&o. arcsin=. _1&o. hav=. [: *: [: sin -: NB. haversine radius=: 1: lat=. {. lon=. {: phis=: ,&lat NB. latitude extracted from x and y lambdas=: ,&lon NB. longitudes ... extracted from x and y NB. (LAT1, LON1) arclen (LAT2, LON2) in radians arclen=. (2 * radius * [: arcsin [: %: hav@:([: -/ phis) + (*&:cos/@:phis) * hav@:([: -/ lambdas))~ NB. optimizing f=:([: mp 2 -/\ (0 0,0 90,:90 0)Deg&(arclen"1))@:deg NB. found 35.2644 45 define angle of radius through center of octant 0 0,0 90,:90 0 eks=. * (cos@:lat * cos@:lon) why=. * (cos@:lat * sin@:lon) zee=. * sin@:lat cart=. 1&$: :(eks, why, zee)"1 0 1 NB. Cartesian triple results from RADIUS cart LATITUDE , LONGITUDE avg=. # *inv +/ NB. four/dollar and five/percent keys do not work NB. PING is center of the spherical triangle and came from optimization of NB.([: mp 2 -/\ (3 2$0 0 0 1.57079632679489656 1.57079632679489656 0)&(arclen"1))@:deg PING=. 35.2644 45 NB. hexahedral corners NB. Cartesian coordinates Spherical coordinates C=. (_8 {. 7 # 1r2) (cart"0 1 deg) _2 [\ 0 0 0 0 0 90 0 45 90 0 45 0 45 90 , PING D=. (DR=. 8 $ 1r2 1) (cart"0 1 deg) DS=. _2 [\ 0 0 0 0 0 45 0 45 45 0 45 0 , PING , PING E=. (ER=. 8 $ 1r2 1) (cart"0 1 deg) ES=. _2 [\ 0 45 0 45 0 90 0 90 , PING , PING , 45 90 45 90 F=. (FR=. 4 # 1r2 1) (cart"0 1 deg) FS=. _2 [\ 90 0 45 0 45 90 , PING , 90 0 45 0 45 90 , PING NB. append mid-edges as the mean of their segment endpoints enabling i38k for the final placement...? NB. Average the spherical coordinates for edges on the surface. edges=. 12 2$0 1 0 2 1 3 2 3 0 4 1 5 2 6 3 7 4 5 4 6 5 7 6 7 NB. C doesn't change D=. D , > (1 1 -:"1 edges { DR)} (avg"2 edges { D) ,:&:{ 1 (cart"1 deg) avg"2 edges { DS NB. D , avg"2 edges { D NB. E=. E , > (1 1 -:"1 edges { ER)} (avg"2 edges { E) ,:&:{ 1 (cart"1 deg) avg"2 edges { ES F=. F , avg"2 edges { F NB. perhaps some of these curves should be great circles. D=. (m2=.(* %@:mp)avg(_2{edges){D) _2}D E=. (m1=.(* %@:mp)avg(_2{edges){E) _2}E F=. (m2,:m1)_2 _1}F i=. _4 F=. (cart@:deg 90*3r4 0) i}F i=. _3 F=. (cart@:deg 90*3r4 1) i}F f=. algebraic_negative@:literate d=. f@:d3edges m012=. ((|: C) ; 0) i388`i38k`prepare38k Mesh n0, n1, n2 m312=. ((|: D) ; 0) i38k`i38k`prepare38k Mesh n3, n1, n2 m302=. ((|: E) ; 0) i38k`i38k`prepare38k Mesh n3, n0, n2 m013=. ((|: F) ; 0) i38k`i38k`prepare38k Mesh n0, n1, n3 NB. gnuplot> splot'a'w points pointtype 5 NB. '/tmp/a' create f 1 {:: m012 NB. '/tmp/a' append f 1 {:: m312 NB. '/tmp/a' append f 1 {:: m302 NB. '/tmp/a' append f 1 {:: m013 NB. gnuplot> splot'a'u 1:2:3 w l '/tmp/a'create d m012 '/tmp/a'append LF,LF '/tmp/a'append d m312 '/tmp/a'append LF,LF '/tmp/a'append d m302 '/tmp/a'append LF,LF '/tmp/a'append d m013 NB. edit a lot then read into prep7 NB.'/tmp/b' create export_ansys m012 NB.'/tmp/b' append export_ansys m312 NB.'/tmp/b' append export_ansys m302 NB.'/tmp/b' append export_ansys m013 )