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
)