Gradient for surface FEM

General Discussion
1

Hi,

There are examples of surface FEM for a scalar function, such as LapEigenBeltrami.edp, however I haven’t found an example of “expressing the gradient of a vector on a surface mesh ThS”.

For a scalar u, the surface gradient is computed as follows:

meshS ThS=square3(nx,ny,[torex,torey,torez],removeduplicate=true) ;

fespace Vh(ThS,P1);

Vh u.

macro Grad3(u) [dx(u),dy(u),dz(u)]  // EOM

real sigma = -1;
varf aS(u,v) = int2d(Th)(Grad3(u)'Grad3(v)- sigma uv);
varf mS(u,v) = int2d(Th)( uv);

But how to express the gradient of a vector (u1, u2) (intrinsic form: two components in the local coordinates)?

fespace VhS(Th, [P1, P1] );
Vh [u1, u2] .

My week form in the local coordinate evolves the computation of

image

where

image

with g_ij being the metric tensor, g^ij being its inverse, and " ; “ denotes the covariant derivative. Alternatively, the non-symmetric form may also be used if it’s easier:image

How to express these expressions on the surface mesh?

For the following code,

fespace VhS(Th, [P1, P1] );
Vh [u1, u2] .

If we express u1 and u2 using the shape function \phi, the derivative ultimately acts on the shape function \phi (x, y, z). The relevant question is: how does FreeFEM compute dx(\phi)?

In some FEM software, it’s only meaningful to define the derivatives with respect to the local coordinates. Can I say

macro Grad3(u) [dx(u),dy(u),dz(u)]  // EOM

has not been projected to the surface yet and the user has to do the projection?

If my assumption is true, then it seems the weak form in exmaple LapEigenBeltrami.edp has a problem: (Grad3(u)'Grad3(v) includes the normal component, which is wrong?

2

Hi,
For a scalar u defined on a meshS surface, the gradient \nabla u defined as [dx(u),dy(u),dz(u)] is a 3d vector which is tangent to the meshS. In means that interpreted as a element of the tangent space, it is the gradient in the Riemanian sense. In particular, if \bf x is on the surface and \bf y also an the surface but approaching \bf x, one has u(\bf y)\simeq u(\bf x)+\nabla u\cdot (\bf y-\bf x).
This is true indendently of any choice of local coordinates.

In the case of a vector field, I am not sure on how to proceed.

A first approach is to define a 3d vector with the space [P1,P1,P1]. I think that you would like to describe a tangent vector field, thus when defining this 3d vector field you have to built it so that it is tangent to the surface.
Then you can use the operators dx(), dy(), dz() on each of the 3 components of the vector field to build up the gradient. You have to keep in mind that only the 3d vector [dx(),dy(),dz()] is geometrically meaningful (indeed dx() gives the x component of the gradient, dy() the y component, dz() the z component).

A second approach is to consider 2 scalar coordinates by decomposing your vector field in a local basis of the tangent space. But then you have to know what is the gradient of this local basis. It could be computed by the first approach.

3

Thank you a lot @fb77 for your comments!
If the vector is defined by an extrinsic form u = (ux, uy, uz) in the Euclidean space, maybe the projection operator P = I - n x n can be used to build the surface gradient ,

func Pk=[P1,P1,P1];
fespace Vh(ThS, Pk);

Vh [u1, u2, u3],[uh1, uh2, uh3];

macro Nor [N.x, N.y, N.z] //
macro vect(u) [u#1, u#2, u#3]//
macro P(u) (vect(u) - (vect(u)'*Nor)*Nor ) //

Because u has to be projected to the tangent space first, then take the gradient, and then project to the tangent space again, i.e.: something like
P (Grad(Pu))

The macro is a bit messy…. and I haven’t got it right yet.

Moreover, I am looking for an intrinsic form with u = u^i \partial_i being expressed locally. It seems that the integration in the parameter space together with metric tensor is more convenient than the surface FEM.

4

The normal to the surface is [Ns.x,Ns.y,Ns.z].
If your metric is attached to the scalar product inherited from R^3, there is hope that your operator D(u) can be expressed in terms of the freefem gradient applied to the 3 components [u1,u2,u3].
Then if [u1,u2,u3] is your unknown and needs to satisfy u\cdot N_s=0, you could use a Lagrange multiplier formulation.
It leads to look for u1,u2,u3 scalars defined on ThS, and also p scalar defined on ThS such that (for example), denoting v1,v2,v3 and q the test functions
\int_S \nabla[u1,u2,u3]:\nabla[v1,v2,v3]+\int_S p [v1,v2,v3]\cdot N_s+\int_S q[u1,u2,u3]\cdot N_s-\int_S f\cdot [v1,v2,v3]=0

Concerning the expression of your “intrinsic” operators, consider \xi_1,\xi_2 coordinates for your surface and \frac{\partial x}{\partial\xi_i} the associate two vectors (i=1,2), that span the tangent space (x stands for the position in R^3).

As a start consider a scalar P. Then \frac{\partial P}{\partial\xi_i}=\nabla P\cdot\frac{\partial x}{\partial\xi_i} for i=1,2. Denote by G the 3\times 2 matrix defined by G_{ij}=\frac{\partial x_i}{\partial\xi_j}.
Then [\frac{\partial P}{\partial {\xi_1}},\frac{\partial P}{\partial {\xi_2}}]=G^t\nabla P
Next for a tangent vector u=[u1,u2,u3] you can apply the previous formula for each component u1,u2,u3 and you get an expression of the partial derivatives \frac{\partial u_i}{\partial\xi_j} in terms of the (matrix) gradient of u=[u1,u2,u3].
Then you need to express your two local coordinates of u in terms of the vector [u1,u2,u3], hopefully obtained via the scalar product of [u1,u2,u3] with the vectors \frac{\partial x}{\partial\xi_i} for i=1,2, and differentiale the result with respect to \xi_i …

5

@fb77 @frederichecht

Dear François and Frederic,

Thank you very much for your previous help. Can I ask you another question about implementation of the surface FEM.

It seems we cannot avoid computing the gradient of the normal n in order to compute the surface gradient of a vector?

image

where u = (ux, uy, uz) is a 3D vector in the ambient Euclidean space, PPreformatted text is the projection operator.

However, since 1. the normal n is piece-wise constant (if I am rignt), its gradient is zero.
In addition, 2. the surface mesh is also not isoparametric (flat triangles),

Can I say it’s really hard to implement the surface FEM (for vector PDEs) in FreeFEM, or you may have a good example?

Although, It’s OK for special cases, such scalar PDEs, or a sphere where P\nabla n P = P. Do we use Lagrange multiplier for the case of general curves instead of computing the gradient of the normal?

6

Do you mean you are interested in the case when there is no analytic description of the surface you discretize ? Hence you have no formula for the gradient of the normal?

7

Yes, for a general surface which evolves.

8

I understand now. What I said above was with the idea of having an analytic description of the surface, for which you can write down a local basis of the tangent plane and differentiate it analytically.
If you have just the vertices coordinates this is more difficult. Ideally one should use finite difference type formulas, and a good theoretical framework to do that is discrete geometry. See for example the book (in French)
https://perso.math.u-pem.fr/romon.pascal/gdd/index.html
But maybe you can do something less “intrinsic” but more simple. For example, you could define a P1 approximation of the normal Ns by interpolation, then differentiate this P1 normal with the dx,dy,dz operators. With this strategy you could maybe have some way to approximate the differential operators you need.

1 Like