Eliminate degrees of freedom using periodic fespace

General Discussion
1

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?

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?

2

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
image1308×872 145 KB

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?

3

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

include "MeshSurface.idp"

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

plot(Th);
1 Like
4

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.

5

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.

6

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

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.

7

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

1 Like
8

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

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

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?

9

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

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

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

1 Like