David Lambert/isoparametric shape functions

From J Wiki
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
)