# Gradient for surface FEM

**URL:** https://community.freefem.org/t/gradient-for-surface-fem/4233
**Category:** General Discussion
**Created:** [March 2, 2026, 10:57am UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233 "2026-03-02T10:57:05Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [March 2, 2026, 10:57am UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/1 "2026-03-02T10:57:05Z")

</div>

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:

```auto
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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/7/7c87c0f434e10c6ccbfd08c404fae4034f87470d.png)

where

![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/1b8a329d5434f307c82b41c2cbee96dffac57eee.png)

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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/cd23e4fb636cce2705358d40d9818581108b6a9e.png)

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

```auto
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?

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 2, 2026, 5:57pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/2 "2026-03-02T17:57:11Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [March 2, 2026, 8:11pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/3 "2026-03-02T20:11:22Z")

</div>

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 ,

```auto
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.

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 3, 2026, 12:56pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/4 "2026-03-03T12:56:32Z")

</div>

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 …

---

<div class="post-metadata">

### Author: ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [March 14, 2026, 11:38am UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/5 "2026-03-14T11:38:18Z")

</div>

@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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/19c7ba8a0d41a10ae81d9e87128cd8004e9f4193.png)

where `u = (ux, uy, uz)` is a 3D vector in the ambient Euclidean space, P`Preformatted 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?**

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 14, 2026, 2:49pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/6 "2026-03-14T14:49:21Z")

</div>

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?

---

<div class="post-metadata">

### Author: ![yongxing](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/yongxing/32/3632_2.png) [@yongxing](https://community.freefem.org/u/yongxing)
#### Post date: [March 14, 2026, 4:00pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/7 "2026-03-14T16:00:39Z")

</div>

Yes, for a general surface which evolves.

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [March 14, 2026, 9:29pm UTC](https://community.freefem.org/t/gradient-for-surface-fem/4233/8 "2026-03-14T21:29:44Z")

</div>

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](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.
