tfem_p1.mac - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
 (HTM) git clone git://src.adamsgaard.dk/pism
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tfem_p1.mac (1307B)
       ---
            1 /* -*- mode: maxima -*- */
            2 
            3 /* Derivation of the Jacobian of the map from the reference P1 element
            4 to physical elements embedded in a Q1 element aligned with the x and y
            5 axes. */
            6 
            7 /* Coordinates of the nodes of a Q1 rectangular element aligned with the coordinate system. */
            8 pts : [
            9   [-dx/2, -dy/2],
           10   [ dx/2, -dy/2],
           11   [ dx/2,  dy/2],
           12   [-dx/2,  dy/2]]$
           13 
           14 /* Short-cuts for nodes */
           15 A : pts[1]$
           16 B : pts[2]$
           17 C : pts[3]$
           18 D : pts[4]$
           19 
           20 /* Four triangles embedded in a rectangular element */
           21 tris : [
           22   [A, B, D],
           23   [B, C, A],
           24   [C, D, B],
           25   [D, A, C]
           26 ]$
           27 
           28 /* P1 shape functions, using the reference element defined above */
           29 chi[1](xi,eta) := 1 - xi - eta$
           30 chi[2](xi,eta) := xi$
           31 chi[3](xi,eta) := eta$
           32 
           33 /* Map from the reference element. */
           34 x(xi,eta) := 'sum(x[j]*chi[j](xi,eta), j, 1, 3)$
           35 y(xi,eta) := 'sum(y[j]*chi[j](xi,eta), j, 1, 3)$
           36 
           37 J : matrix(
           38   ['diff(x(xi,eta), xi), 'diff(y(xi,eta), xi)],
           39   ['diff(x(xi,eta), eta),'diff(y(xi,eta), eta)])$
           40 
           41 /* Compute the Jacobian of the mapping from the reference element */
           42 J_s : ratsimp(ev(J, nouns))$
           43 
           44 eqs[k] := block([T],
           45   T : tris[k],
           46   [
           47   x[1] = T[1][1],
           48   y[1] = T[1][2],
           49   x[2] = T[2][1],
           50   y[2] = T[2][2],
           51   x[3] = T[3][1],
           52   y[3] = T[3][2]
           53   ])$
           54 
           55 J[k] := ev(J_s, eqs[k]);
           56 
           57 J_det[k] := determinant(J[k]);
           58 
           59 J_inv[k] := invert(J[k]);
           60 
           61 J[1];
           62 J[2];
           63 J[3];
           64