David Lambert/isoparametric shape functions lab

From J Wiki
Jump to navigation Jump to search

Change the file loaded by load'/home/lambertdw/src/j/isoparametric_shape_functions.ijs' to your file created from the isoparametric_shape_functions page.

LABTITLE=: 'Isoparametric Shape Functions'

NB. =========================================================
Lab Section Overview

 Here we develop a general method to generate isoparametric interpolants.  These can interpolate data on non-uniform grids.  The interpolants will be of linear through cubic order in each dimension.  The dimensions range from 1 to 3; the interpolants are first through third order.

 Isoparametric representation means to use the same nodes for both geometry and displacements.  This property is essentially non-essential here because I plan to use these functions to create finite element meshes.  Meshing requires a match between faces containing the same nodes in the interpolated values.

 Merging meshes requires that co-located nodes are within a small floating point tolerance of each other.  Use the serendipity elements for mesh frames.  Serendipity elements have no internal nodes.  The cubic nodes and complete-term interpolants here are not serendipity elements.

)
PREPARE

load'/home/lambertdw/src/j/isoparametric_shape_functions.ijs'

PREPARE

NB. =========================================================
Lab Section Data: interpolating a field

 We wish to esimate the temperature at Whiteman Air Force Base.  Uninspired, but this works.  We know the location of Whiteman AFB.  We also know the location and temperature at 4 surrounding cities which form a convex quadrilateral containing Whiteman AFB.  The locations are approximate and in graph paper units.  The master element is non-dimensional.  Springfield was choses as the origin to simplify reading the graph paper.

 Location         X      Y      T
  KC               _13    25     33
  Springfield        0     0     36
  JC                12    18     36
  Columbia          10    23     34
  AFB               _3    20    ????
)

DATA=: noun define  NB.3 fields aligned by column
   -13    25      33
     0     0      36
    12    18      36
    10    23      34
    -3    20
)

NB. The interpolant requires rows of same field
FIELDS=: |: 0&".;._2 DATA

(,.'XYT');(0 _1&}. ; _ _1&{.)FIELDS

NB. =========================================================
Lab Section Interpolation verb: interpolating a field

 Flat Earthers can love this 2D problem.  The known points form a quadrilateral containing the AFB.  Thus it is an interpolation (opposed to extrapolation) problem.  Extrapolated results should always be cross-checked.

 The interpolant (i) is 2d (2), has 4 corners (4), and 4 nodes (4).  i244 is appropriate.

 The order of the data is important, and this module uses an unusual node order.  This affects the order in which to specify the data.  The quadrilateral has nodes

    2     3


    0     1

 The cities where all fields are known can be in the order (looking at a map)
 KC, Springfield, Columbia, JC.

)

NB. =========================================================
Lab Section Solution: interpolating a field

 The Interpolate2d conjunction manages the data and computes the interpolated temperature at Whiteman AFB.

)

NB. Coordinates of Known Temperatures
[COKT=: ((<0 ;0 1 3 2),(<1;0 1 3 2)){FIELDS

NB. coordinates of desired temperature
[CODT=: , 2 _1{.FIELDS

NB. Known Temperatures
[KT=: _1 4 {. FIELDS

NB. Estimate the temperature at Whiteman AFB, the solution!
KT COKT Interpolate2d i244 CODT
NB. (actual temperature was 37 degrees F.)

NB. =========================================================
Lab Section Master coordinates: interpolating a field

 The master coordinates can also be found.  The element contains the location if the master coordinates are all <: 1 in absolute value

)

NB. Parenthesize the verb for clarity
NB. as a dyad yields the interpolated value.
KT (COKT Interpolate2d i244) CODT

NB. as a monad yield the master coordinates.
(COKT Interpolate2d i244) CODT

NB. =========================================================
Lab Section 4 node square master element

 Example for 4 node quadrilateral element.  A colorful picture appears at the end of this lab.

 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 within the element.  Name these shape functions n_0 through n_3.  The n_i satisfy the Kronecker property that it is 1 at node i and 0 at the other nodes.

 Let the base element have nodal coordinates (xi, eta) of +/-1.

           f2              f3 (1,1)
           N2──────────────N3
           │   f(x,y)      │
           │   (x,y)       │
           │  *     (0,0)  │
           │       +       │
           │               │
           │               │
           │(_1,_1)        │
           N0──────────────N1
           f0              f1

)

NB. =========================================================
Lab Section Choose shape functions

 The quadrilateral is planar, so there are two variables.
 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.
 n2(-1, 1) = 1, n1(all other corners) is 0.
 n3( 1, 1) = 1, n1(all other corners) is 0.

 Choose terms for the shape functions.  Each function is 1 at a node and 0 at 3 other nodes.  There are 4 equations to solve and we need 4 terms.  The product of two linear functions is the usual choice.  (a + b xi) ( c + d eta) == 1
   
 ni has the form C0 + C1*xi + C2*eta + C3*xi*eta .
 C_i are determined because given that the function values f_i at the nodes are known.

   This method extends to higher order interpolants having more nodes or to other dimensions.

)

NB. =========================================================
Lab Section Node order

 This module uses non-standard node ordering for the elements.  Reordering is necessary to communicate with most other external programs.  For the n-dimensional linear elements use only the corner nodes.  Use additional nodes, mid-edge or mid face, for various quadratic interpolants.  The cubics use a different order.

   Node order
   1D:

   0   2   1


   2D:

   2   7   3

   5   8   6   Node 8 at origin, Node 3 at (1,1)

   0   4   1

   3D
   At z = _1           z = 1            z = 0
   2   b   3           6   j   7        e   o   f

   9   k   a           h   p   i        m   q   n

   0   8   1           4   g   5        c   l   d
)

NB. =========================================================
Lab Section 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 nodes
   i249:   2 dimensions, 4 corners (quadrilateral), 9 nodes, or
   C244:   C constants for   2 dimensions, 4 corners (quadrilateral), 4 nodes

)

NB. =========================================================
Lab Section Confusing table

 This table shows reasonable choices of terms and how to access them.  I imagine it was created with editing applied to {4#<'1xyz' used to generate symbolic multi-dimensional polynomial terms.  Shown are j accessors.

   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
)

NB. Terms of the 27 node brick, i38r
<"2 ;L:1{(;:'1 z zz');(;:'1 y yy');<(;:'1 x xx')

36br   NB. "r" identifier has 27 nodes.
       NB. (Incidental fun with base 0, also true that 27 = 0br)

NB. =========================================================
Lab Section Node order of cubic interpolants

 The 4 node line has nodes numbered 0 1 2 3.
 Nodes 0 and 3 are the ends.

 The 16 node face has nodes numbered i. _4 4
 0 3 12 15 are the corners.

 The 64 node brick in nodal order i. _4 _4 4
 0 3 12 15 48 51 60 63 are the corners.

)

i. _4 _4 4

NB. =========================================================
Lab Section Meshing tips.

 It happens to be quite convenient for construction and debugging that if N are the 27 known values to be interpolated of an "r" element, then 8&{.N are the corners for an "8" element, 20&{.N are the corners and mid-side nodes for a "k" element, and 26&{.N are in correct order to include the faces of a "q" element as well.

 Some other crazy scheme for the cubics would be difficult to remember, and would not help at all with the exploitable convenience for the lower order elements.

)

NB. =========================================================
Lab Section Interpolants

 Rows of y are the coordinates in the element system (xi, eta, nu) where to interpolate the field.  The number of columns matches dimensionality of the element, 1, 2, or 3.

 Rows of x are the known values at the nodes to interpolate.  Hence the number of columns equals the number of nodes in the element, and of course the data order must match the node order in the element.

 Both x and y should have rank 2.  There can be many points, and many fields can be interpolated at once.

 Extrapolation---specifying the master coordinates outside of the interval from _1 to 1 is a bad idea, except for the 2 node line element i122 which is fully linear.

 In following example i122 evaluates 8 fields at xi = 0, 1, 2, 3, and 4.  The known values of the first field are 0 1, of the second field are 2 3, etcetera.

)
(i.8 2) i122 i.5 1

NB. =========================================================
Lab Section Interpolation output explanation.

 The first row are the eight field values interpolated at xi=0.  xi=0 is the element midpoint, and the result is indeed the average of the field values.  The second output row are the values at xi=1, which is the location of the second node of the master element.

 Read by column see a single field value interpolated at each node.  The output shape is ((rows of y) , (rows of x))

 For meshing, y can be an easily produced uniform grid on _1 to 1 (see the odometer essay), and the known field values x can be real world coordinates X,Y,Z onto which the grid should be mapped.  To map to a curved space use an element having mid-edge nodes.  For a graded mesh I start with a uniform grid, which I map onto _1 1 having midside nodes repositioned along the line between the corners.  Transposing the result becomes a non-uniform grid on _1 to 1 which I then map onto world coordinates.

)

NB. =========================================================
Lab Section Method Summary, neglecting critical details.

 Choose some nodes.  Choose some shape functions having terms consisting of varioius combinations of the local coordinates.  Solve the linear systems of equations that each shape function is 1 at "its" node and 0 at the others.

)

NB. =========================================================
Lab Section Interpolant verbs

 Let F be the number of fields and N be the number of points.
 The shapes input and -> output are shown

 Line segments
  (F,2) i122 (N,1) -> (N,F)   Ends
  (F,3) i123 (N,1) -> (N,F)   Ends and midpoint
  (F,4) i124 (N,1) -> (N,F)   End point point End
 Quadrilaterals
  (F, 4) i244 (N,2) -> (N,F)  Corners
  (F, 8) i248 (N,2) -> (N,F)  Corners and midedge
  (F, 9) i249 (N,2) -> (N,F)  Corners and midedge and face center
  (F,16) i24g (N,2) -> (N,F)  y=_1 by 4, y=_1r3 by 4, y=1r3 by 4, y=1
 Hexahedrons
  (F, 8) i388 (N,3) -> (N,F)  Corners
  (F,20) i38k (N,3) -> (N,F)  Corners and midedge
  (F,26) i38q (N,3) -> (N,F)  Corners and midedge and face centers
  (F,27) i38r (N,3) -> (N,F)  Corners, midedges, face centers, body center
  (F,64) i38z (N,3) -> (N,F)  i._4 _4 4.

 The verb prepare_38z constructs a 3 field x argument for i38z to squeeze a mesh toward the center of the cube (positive arguments less than 1r3) or squeeze toward the edges.

)

require'viewmat'
GRID =: |.,~"0/~(%~i:)100
SADDLE =: 1.0 2.0 2.2 0.7 i244 GRID
viewmat SADDLE
(<./ , >./) , SADDLE