# Eliminate degrees of freedom using periodic fespace

**URL:** <https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238>\
**Category:** General Discussion\
**Created:** [March 9, 2026, 11:45am UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238 "2026-03-09T11:45:18Z")\
**Posts on this page:** 9\
**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 9, 2026, 11:45am UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/1 "2026-03-09T11:45:18Z")

</div>

`current line = 163`  
`Assertion fail : (kkk++ < 10)`  
`line :1141, in file lgfem.cpp`  
`Segmentation fault (core dumped)`

If we use the parameterised space (`\theta, \phi)` to solve PDEs on a sphere, then the boundaries at `\phi=0` and `\phi=\pi` should really be collapsed to single points – we may try various tricks and boundary conditions to, but none of them seems entirely satisfactory.

Anyway, I am trying to implement this using periodic boundary conditions. However, it seems that the maximum number of periodic pairs is **10**. Is this really the case in FreeFEM?

```auto
int m = 8;

// --------------------------------------------------
// geometry arrays
// --------------------------------------------------

real[int] xxb(m+1), yyb(m+1);
real[int] xxt(m+1), yyt(m+1);

int[int] nbot(m), ntop(m);
real[int] lbot(m), ltop(m);

macro dist(ax,ay,bx,by) sqrt(square(ax-bx)+square(ay-by)) //EOM

// bottom boundary
for(int i=0;i<=m;i++){
xxb[i] = 2pii/m;
yyb[i] = phi1;
}

for(int i=0;i<m;i++){
nbot[i] = 1;
lbot[i] = dist(xxb[i],yyb[i],xxb[i+1],yyb[i+1]);
}

border Bot(t=0,1;i){
x = xxb[i]*(1-t) + xxb[i+1]*t;
y = yyb[i];
label = i+1;
}

// top boundary
for(int i=0;i<=m;i++){
xxt[i] = 2pi(m-i)/m;
yyt[i] = phi2;
}

for(int i=0;i<m;i++){
ntop[i] = 1;
ltop[i] = dist(xxt[i],yyt[i],xxt[i+1],yyt[i+1]);
}

border Top(t=0,1;i){
x = xxt[i]*(1-t) + xxt[i+1]*t;
y = yyt[i];
label = 101 + i;
}

// side boundaries
border cr(t=phi1,phi2){
x = 2*pi;
y = t;
label = 1001;
};

border cl(t=phi2,phi1){
x = 0;
y = t;
label = 1002;
};

// build mesh
mesh Th = buildmesh(Bot(nbot) + cr(m/2) + Top(ntop) + cl(m/2));

plot(Th,wait=1);

// --------------------------------------------------
// periodic parametrisation
// --------------------------------------------------

macro PERIOBOT(k)
[k,abs(x-xxb[k-1])/lbot[k-1]] //EOM

macro PERIOTOP(k)
[100+k,abs(xxt[k-1]-x)/ltop[k-1]] //EOM

macro BOTPAIR(k) PERIOBOT(k-1),PERIOBOT(k) //EOM
macro TOPPAIR(k) PERIOTOP(k-1),PERIOTOP(k) //EOM

// --------------------------------------------------
// periodic list
// --------------------------------------------------
func perio = [
[1001,y],
[1002,y],
BOTPAIR(2), BOTPAIR(3), BOTPAIR(4), BOTPAIR(5), BOTPAIR(6),
BOTPAIR(7), BOTPAIR(8),
TOPPAIR(2), TOPPAIR(3), TOPPAIR(4), TOPPAIR(5), TOPPAIR(6),
TOPPAIR(7), TOPPAIR(8)];

//////////////////////////////////////////////////////////////////////
// Finite element spaces with periodicity
//////////////////////////////////////////////////////////////////////

func Pk=[P1,P1];
fespace Rh(Th, Pk, periodic=perio);

```

**The code looks ugly, but it works up to m=8 at least.** The following error appears when set m=9:

current line = 130  
Assertion fail : (kkk++ \< 10)  
line :1141, in file lgfem.cpp  
Segmentation fault (core dumped)

**Doe someone have a better idea to collapse the degrees of freedom, or general ideas to deal with singular boundaries ? Thank you very much in advance!**

By the way, years ago when I worked on a software package called FEPG (Finite Element Programmer Generator), the periodic conditions were just handled simply as constraints. If all constraints are eliminated during the assembly of the global matrix, then periodic conditions can be handled quite conveniently at that stage. However, I don’t think this is the way FreeFEM or most other FEM software packages deal with constraints?

---

<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 9, 2026, 6:18pm UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/2 "2026-03-09T18:18:12Z")

</div>

Seems this is a good idea:

first, pairing the two line segments with the corresponding nodes.

then shift one element and pairing then again as shown below:

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

This should collapse all the nodes to one single node. However, it seem FreeFEM needs to pair the edge elements as well… anyone know how to implement this idea?

---

<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 10, 2026, 10:12am UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/3 "2026-03-10T10:12:10Z")

</div>

The difficulty with spherical coordinates is the singularities at the poles, that lead to the problem of merging the degrees of freedom.  
I think it is interesting to use other types of meshes without singularity.  
There is in particular the mesh built from an icosahedron that can be obtained as

```auto
include "MeshSurface.idp"

real radius=1.;
int nbsplit=5;
int orient=1;
meshS Th=Sphere20(radius,nbsplit,orient);

plot(Th);

```

---

<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 10, 2026, 10:29am UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/4 "2026-03-10T10:29:52Z")

</div>

Thank you a lot. The mesh looks really nice, and the surfaceFEM works well on this mesh.

However, I am still looking for a formulation on the parameter space…

It seems that the singularity is not a problem, because the Gaussian quadrature points never touch the pole.

I am surprised that collapsing the nodes can’t be implemented in FreeFEM.

---

<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 10, 2026, 11:18am UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/5 "2026-03-10T11:18:22Z")

</div>

About constraints, there are mainly two ways to deal with them:  
– Use a formulation with Lagrange multipliers. It adds unknowns to the system but it is generally handable.  
or  
– Modify by hand the matrix and right-hand side of the system. For this there is the tool setBC() that is useful.

About merging nodes and eliminating the degrees of freedom, I think that a way to do it could be

1. Start from a meshS Th0 transported from the plane coordinate system. Then the poles appear several times with different node numbers (if necessary one can use a transport map that does not exactly close, so that the image is approximately the sphere minus a meridian).  
From this write data (by I/O commands) to a file .mesh to the define a new meshS where the redundant nodes are deleted, and as a consequence the node numbering is changed. During this operation we have to keep track of the o2n mapping (old node number)–\>(new node number).
2. Load the file .mesh to a new meshS Th. Then if you have a matrix and rhs corresponding to a varational formulation on Th0, you can build the new matrix and new rhs corresponding to Th by using the o2n mapping, thus removing the redundant degrees of freedom.

A maybe more direct method is as follows by Lagrange multiplier. Consider that your have the system AX=b, where X is the vector of degrees of freedom. Assume that there is a “pole” corresponding to n degrees of freedom. The list of indices for these degrees of freedom is I(0),I(1),\ldots I(n-1). We want to add the n-1 constraints X\_{I(k)}=X\_{I(n-1)} for k=0,\ldots,n-2.  
We add n-1 new unknowns (the multipliers), leading to the matrix  
M=\pmatrix{A & B^t\\ B & 0}  
for solving the system M(X,Y)=(b,0).  
The matrix B is for 0\leq i\leq n-2 and j in the range of indices of X  
B\_{ij}=\delta\_{I(i),j}-\delta\_{I(n-1),j}  
with \delta the Kronecker symbol.  
Then solving the new system means to write the variational formulation with the constraints on the unknown, for only the test functions that satisfy the constraints.

---

<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 10, 2026, 1:47pm UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/6 "2026-03-10T13:47:13Z")

</div>

Yes, I thought about the Lagrange multiplier method, and enforcing a “zero gradient along the boundary“ in the weak form may lead to constant values on the boundary.

I quickly tried the penalty approach, but haven’t tried the Lagrange multiplier yet…

Ideally, I want to remove the degrees of freedom. One reason is that I am also looking at an eigenvalue problem corresponding to the Hodge Laplacian on the surface. The eigenvalue 0 corresponds to the three Killling vector fields, see animation here: [SurfaceFluids](https://yongxingwang.github.io/surfacefluids/)

> **[SurfaceFluids](https://yongxingwang.github.io/surfacefluids/)**
>
> Yongxing’s personal page.

I am not sure whether the Lagrange multiplier would change the eigenvalue problem, but I will give it a try. I really appreciate your help.

---

<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 10, 2026, 2:32pm UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/7 "2026-03-10T14:32:23Z")

</div>

If you want to compute the eigenvalues of AX=\lambda X with constraints, you have to write M(X,Y)=\lambda J(X,Y), with  
J=\pmatrix{Id & 0 \\ 0 & 0}.

---

<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 10, 2026, 4:19pm UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/8 "2026-03-10T16:19:14Z")

</div>

@fb77, I am now trying an eigenvalue problem using the Lagrange multiplier method – the new composite fespace makes the assembly very convenient, **however…**

```auto
mesh Th = square(m,m) ;

int[int] labs = [1,3];
meshL ThL = extract(Th, label=labs);

func Pk=[P1,P1];
fespace Rh(Th, Pk, periodic=perio);
Rh [u1, u2];

fespace Lh(ThL, Pk);
Lh [a1,a2];

fespace RLh = Rh * Lh;
varf Hodge(<[u1,u2],[a1, a2]>,<[uh1,uh2], [ah1, ah2]>) = …;

varf bf(<[u1,u2],[a1, a2]>,<[uh1,uh2], [ah1, ah2]>) = … ;

matrix A = Hodge(RLh,RLh); 
matrix B = bf(RLh,RLh); 

int k=EigenValue(A,B,sym=true,sigma=sigma,value=ev,vector=eu1,tol=1e-16,maxit=0,ncv=30);

```

However, do you know how to define **vector** in the above EigenValue ( )?

Usually we do

```auto
real sigma = 1.e-12;
int nev=10; // number of computed eigen valeu close to sigma

real[int] ev(nev); // to store nev eigein value

Rh[int] [eu1,eu2] (nev); // to store nev eigen vector

```

now we have do something like

`RLh[int] [eu1, eu2, ea1, ea2] (nev)`

but this expression is illegal in FreeFEM, because `RLh=Rh * Lh` , any idea how to resolve this issue?

---

<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 10, 2026, 5:32pm UTC](https://community.freefem.org/t/eliminate-degrees-of-freedom-using-periodic-fespace/4238/9 "2026-03-10T17:32:55Z")

</div>

You should use the array of dofs instead of the finite element functions. I think it is not possible to do it with `EigenValue()`. However it is possible if you use `EPSSolve()` with PETSc. Then instead of `vector=` you put `array=`.  
There is an example in

> [@Generalized eigenvalue problem solving with Lagrange multipliers with EPSSOLVE](https://community.freefem.org/t/generalized-eigenvalue-problem-solving-with-lagrange-multipliers-with-epssolve/3650):
>
> Hello, I am facing a problem with a generalized eigenvalues problem using EPSSolve in FreeFEM. I want to solve a variational formulation with constraints, which I handle using Lagrange multipliers. This leads to a matrix problem of the form AX = \lambda BX. (B has some lines of zeros but it shouldn’t be a problem to solve this matricial equation) I can solve this problem in MATLAB, but I cannot get it to work with the EPSSOLVE function in FreeFEM. If I specify the parameters for EPSSOLVE as s…

For the case of composite spaces, you can do as  
`real[int,int] Vectab(ndof,nev);`  
then you put `array=Vectab` in `EPSSolve()`.  
Finally you define

```auto
Rh [uu1,uu2];
Lh [aa1,aa2];
for (int i=0;i<nev;i++){
 [uu1[],aa1[]]=Vectab(:,i);
 //here you can use the eigenvector [uu1,uu2] as you like
 //the Lagrange multiplier [aa1,aa2] is useless
}

```

There is also this example

> [@Fascinating theoretical challenge: Simulating EM eigenmodes in a tetrahedron](https://community.freefem.org/t/fascinating-theoretical-challenge-simulating-em-eigenmodes-in-a-tetrahedron/4049/10):
>
> Dear Frodo, I have written a version with Laplacian and boundary condition E\times n =0 [EMtet.edp](https://community.freefem.org/uploads/short-url/AobXDpCZsR0nsGFhn5FxXqcv0sg.edp) (7.8 KB) I get for the 10 first eigenvalues 78.88776433 79.02651853 79.02651853 153.9973688 163.2695866 163.7718604 163.7718604 174.020952 175.1561748 197.8465549 It seems that there is a first eigenvalue with triple multiplicity, then an eigenvalue 154 (your fundamental mode ?), then again an eigenvalue with triple multiplicity. You have to run the code EMtet.edp to plot (buil…