%	vecalg.s
%
%	Purpose:	Vector algebra package and PDE vector notation
%	Copyright:	R.A. van Engelen, Leiden University, 1997

include linalg.
	% import dot- and cross-product operators

banner("Vector Algebra Package").

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	Define grad, div, curl, and lapl
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prefix grad has latex_name("\nabla").
prefix div has latex_name("\nabla\cdot").
prefix curl has latex_name("\nabla\times").
prefix lapl has latex_name("\nabla^2").

% Note: Vectors `coordinates' and `coefficients' define a coordinate system.
%       Scalar operator `diff(X, Y)' represents differentiation of X wrt. Y.
%       In def. of `grad X' vectors `coordinates' and `coefficients' are mapped.

grad X := diff(X, coordinates) / coefficients.

div X := diff(X, coordinates : 1) / coefficients : 1.
(div [X1, X2] :=
( ( diff(coefficients : 2 * X1, coordinates : 1)
+ diff(coefficients : 1 * X2, coordinates : 2)
)
/ (coefficients : 1 * coefficients : 2)
))
.
(div [X1, X2, X3] :=
( ( diff(( coefficients : 2
* coefficients : 3
* X1
), coordinates : 1)
+ diff(( coefficients : 3
* coefficients : 1
* X2
), coordinates : 2)
+ diff(( coefficients : 1
* coefficients : 2
* X3
), coordinates : 3)
)
/ (coefficients : 1 * coefficients : 2 * coefficients : 3)
))
.
curl _G2509 := error("Illegal use of curl"). % Scalar argument? (curl [X1, X2, X3] :=
[ ( ( diff(coefficients : 3 * X3, coordinates : 2)
- diff(coefficients : 2 * X2, coordinates : 3)
)
/ (coefficients : 2 * coefficients : 3)
)
, ( ( diff(coefficients : 1 * X1, coordinates : 3)
- diff(coefficients : 3 * X3, coordinates : 1)
)
/ (coefficients : 3 * coefficients : 1)
)
, ( ( diff(coefficients : 2 * X2, coordinates : 1)
- diff(coefficients : 1 * X1, coordinates : 2)
)
/ (coefficients : 1 * coefficients : 2)
) ])
.
lapl X := div grad X. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Vector notation % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Note: `nabla' represents the gradient symbol. It is used to formulate a PDE % in full vector notation with dyadic operators `+', `*', and `^'. nabla is inert. nabla_square is inert. nabla_plus(_G2509) is inert. nabla_times(_G2509) is inert. nabla_times_plus(_G2509, _G2510) is inert. nabla_square_plus(_G2509) is inert. nabla_square_times(_G2509) is inert. nabla_square_times_plus(_G2509, _G2510) is inert. nabla_dot(_G2509) is inert. nabla_cross(_G2509) is inert. nabla_times_dot(_G2509) is inert. nabla_times_cross(_G2509) is inert. % Extend definition of addition nabla + [X | Xs] := nabla_plus([X | Xs]). nabla_square + X := nabla_square_plus(X). nabla_plus([Y | Ys]) + [X | Xs] := nabla_plus([X | Xs] + [Y | Ys]). nabla_plus([Y | Ys]) + X := nabla_plus(X + [Y | Ys]). nabla_plus(Y) + [X | Xs] := nabla_plus([X | Xs] + Y). nabla_plus(Y) + X := nabla_plus(X + Y). nabla_square_plus(Y) + X := nabla_square_plus(X + Y). nabla_times(Y) + [X | Xs] := nabla_times_plus(Y, [X | Xs]). nabla_times(Y) + X := nabla_times_plus(Y, X). nabla_square_times(Y) + X := nabla_square_times_plus(Y, X). (nabla_times_plus(Y, [X | Xs]) + [Z | Zs] :=
nabla_times_plus(Y, [X | Xs] + [Z | Zs]))
.
nabla_times_plus(Y, [X | Xs]) + Z := nabla_times_plus(Y, [X | Xs] + Z). nabla_times_plus(Y, X) + [Z | Zs] := nabla_times_plus(Y, X + [Z | Zs]). nabla_times_plus(Y, X) + Z := nabla_times_plus(Y, X + Z). nabla_square_times_plus(Y, X) + Z := nabla_square_times_plus(Y, X + Z). [X | Xs] + nabla := nabla_plus([X | Xs]). X + nabla_square := nabla_square_plus(X). [X | Xs] + nabla_plus([Y | Ys]) := nabla_plus([X | Xs] + [Y | Ys]). X + nabla_plus([Y | Ys]) := nabla_plus(X + [Y | Ys]). [X | Xs] + nabla_plus(Y) := nabla_plus([X | Xs] + Y). X + nabla_plus(Y) := nabla_plus(X + Y). X + nabla_square_plus(Y) := nabla_square_plus(X + Y). [X | Xs] + nabla_times(Y) := nabla_times_plus(Y, [X | Xs]). X + nabla_times(Y) := nabla_times_plus(Y, X). X + nabla_square_times(Y) := nabla_square_times_plus(Y, X). ([Z | Zs] + nabla_times_plus(Y, [X | Xs]) :=
nabla_times_plus(Y, [X | Xs] + [Z | Zs]))
.
Z + nabla_times_plus(Y, [X | Xs]) := nabla_times_plus(Y, [X | Xs] + Z). [Z | Zs] + nabla_times_plus(Y, X) := nabla_times_plus(Y, X + [Z | Zs]). Z + nabla_times_plus(Y, X) := nabla_times_plus(Y, X + Z). Z + nabla_square_times_plus(Y, X) := nabla_times_plus(Y, X + Z). % Extend definition of product nabla * X := grad X. nabla_square * X := lapl X. nabla_plus([X | Xs]) * [Y | Ys] := [X | Xs] * [Y | Ys] + grad [Y | Ys]. nabla_plus([X | Xs]) * Y := [X | Xs] * Y + grad Y. nabla_plus(X) * Y := X * Y + grad Y. nabla_times(X) * Y := X * grad Y. nabla_times_plus(Y, [X | Xs]) * Z := [X | Xs] * Z + Y * grad Z. nabla_times_plus(Y, X) * Z := X * Z + Y * grad Z. nabla_square_plus(X) * Y := X * Y + lapl Y. nabla_square_times(X) * Y := X * lapl Y. nabla_square_times_plus(Y, X) * Z := X * Z + Y * lapl Z. nabla_dot([X | Xs]) * Y := [X | Xs] .* grad Y. nabla_cross([X | Xs]) * Y := [X | Xs] #* grad Y. nabla_times_dot(X) * Y := div (X * grad Y). nabla_times_cross(X) * Y := div (X * grad Y). X * nabla := nabla_times(X). X * nabla_square := nabla_square_times(X). X * nabla_times(Y) := nabla_times(X * Y). X * nabla_square_times(Y) := nabla_square_times(X * Y). % Extend definition of power nabla ^ 2 := nabla_square. % Extend definition of dot product nabla .* nabla := nabla_square. nabla .* nabla_times(X) := nabla_times_dot(X). nabla .* [X | Xs] := div [X | Xs]. nabla_plus([X | Xs]) .* [Y | Ys] := [X | Xs] .* [Y | Ys] + div [Y | Ys]. nabla_plus(X) .* Y := X * Y + div Y. nabla_times(X) .* [Y | Ys] := X * div [Y | Ys]. nabla_times(X) .* Y := X * div Y. (nabla_times_plus(Y, [X | Xs]) .* [Z | Zs] :=
[X | Xs] .* [Z | Zs] + Y * div [Z | Zs])
.
nabla_times_plus(Y, X) .* Z := X * Z + Y * div Z. [X | Xs] .* nabla := nabla_dot([X | Xs]). % Extend definition of cross product nabla #* nabla_times(X) := nabla_times_cross(X). nabla #* [X | Xs] := curl [X | Xs]. (nabla_plus([X | Xs]) #* [Y | Ys] :=
[X | Xs] #* [Y | Ys] + curl [Y | Ys])
.
nabla_times(X) #* [Y | Ys] := X * curl [Y | Ys]. [X | Xs] #* nabla := nabla_cross([X | Xs]). % Library of commonly used curvilinear grids (coordinate_system("1D Cartesian") :=
proc((coordinates := [x]; coefficients := [1]; return nil)))
.
(coordinate_system("2D Cartesian") :=
proc((coordinates := [x, y]; coefficients := [1, 1]; return nil)))
.
(coordinate_system("3D Cartesian") :=
proc(( coordinates := [x, y, z]
; coefficients := [1, 1, 1]
; return nil
)))
.
(coordinate_system("Polar") :=
proc((coordinates := [r, theta]; coefficients := [1, r]; return nil)))
.
coordinate_system(X) := error("Unknown coordinate system " // sprint(X)).