David Lambert/isoparametric shape functions lab
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