Question: Calculus of variations for a deformable surface. how to solve efficiently a Differential geometry problem.

Hello everyone.

I am not really sure how to solve the following problem:

Given a parametriaztion of a surface as an interpolation in terms of given nodes as:
X(ξ,η)= Σ Νi(ξ,η) Xi, it is possible to calculate tangent and normal vectors as:
α1(ξ,η)= Σ Νi,ξ(ξ,η) Xi     α2(ξ,η)= Σ Νi,η(ξ,η) Xi  and a3(ξ,η)= cross(a1,a2) and their derivatives
α1,ξ α1,η α2,ξ α2,η etc.
For a variation of Xi i.e.  a variation of X can be calculated δX(ξ,η)= Σ Νi(ξ,η) δXi in terms of 
intermediate quantities that are as simple as  Σ Νi(ξ,η) for the position vector but get quite complicate
for the normal vector (δα3 that depends on a cross product) of given vectors  and for higher order derivatives,
or for the second variation. These variations can be also expressed with a directional derivative in the direction 
X,Y, or Z of a node i with initial position Xi.

To calculate those more complicate derivations I try to express them in terms of known less complicate ones.
Symbolizing r and s the first and second variation, and α,β firs and second derivative withh respect to ξ οr h or ξ,η.

Regarding the variation of the normal vector "a3_t" as known (i.e. a3_t,r)
and the variation of its norm "a3_n" (i.e. a3_n,r) as known too we can get the variation of
the of unit normal using the following commands:

a3 := a3_t(a, s, r)*(1/a3_n(a, s, r)): (vector-scalar operation) (1)
diff(a3, r):

the second variation
diff(a3, r):
diff(diff(a3, r), s) (2)

and even variations of its derivatives 
diff(diff(diff(a3, a), r), s): (3) 

That yield the results of the posted picture.

The problem is that if we want to calculate in the same way the variation of the derivative 
of expressions that involve the vector-vector operations i.e  variation of diff(a3_t, a)
where a3_t = CrossProduct(a1(a,r,s),a2(a,r,s)) it is not possible to do that as in (1)
we cannot substitute CrossProduct or dot product operation  with simple multiplication.

If we use the same strategy in VectorCalculus package we can get some results:
using the following commands: 

a1 := PositionVector([a1x(a, s, r), a1y(a, s, r), a1z(a, s, r)], cartesian[x, y, z])
a2 := PositionVector([a2x(a, s, r), a2y(a, s, r), a2z(a, s, r)], cartesian[x, y, z])
a3T := CrossProduct(a1, a2)
diff(a3T, r)

but then (as it is shown in the upload picture 2) the calculations are performed component-wise
are very lengthy so it is impossible to get second derivative or second derivation.

So I would like to ask is it possible to declare a1 as a vector but not give its components a1x, a1y etc explicitely but instead
declare it as a vector valued function of (a, r, s). so that it yields simpler expressions 
as "a3_t,r= a1,r(x)a2 + a1(x)a2,r".

I would also like to know if those calculations can be made simpler with and another package and if so is there an appropriate example
or a book that explains how to do that.

I would really appreciate  any idea.
Because this is a difficult problem Maple can make it very easy.

Thank you in advance.

upload1.mw