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

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