% 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)).