Generalized eigenvalue problem in a mesh

General Discussion
1

Dear FreeFem community,

I hope you are all doing well.

I am trying to solve the following generalized eigenvalue problem on a simplicial mesh \mathcal{T}_h: for each K \in \mathcal{T}_h, find \lambda \in \mathbb{R} and u_h \in \mathbb{P}_k(K)/\mathbb{R} such that

(\nabla \cdot (\mathcal{A} \nabla u_h), \nabla \cdot (\mathcal{A} \nabla v_h))_{0,K} = \lambda\, (\mathcal{A} \nabla u_h, \mathcal{A} \nabla v_h)_{0,K} \quad \forall\, v_h \in \mathbb{P}_k(K)/\mathbb{R}.

Here, \mathcal{A} is a differentiable, elliptic tensor. For each K \in \mathcal{T}_h, I want to compute the local eigenvalues, and extract the maximum one called \lambda_K, which I plan to use in a stabilized method. It is very important for me to associate the maximum eigenvalue to its corresponding element K, since the stabilization may differ from one triangle to another.

Below is my current code for the case \mathcal{A} := I_d:

//----Loop in each triangle
Rh fespace(calP, P0);
Rh lambdaK;
/*
mesh[int] Th(NbTriangles);
for(int K = 0; K < NbTriangles; K++) {

    real coordX0, coordX1, coordX2;
    real coordY0, coordY1, coordY2;

    for(int i = 0; i< NbVertices; i++){
        if(elemNodesGlobal(K,0) == i){
            coordX0 = calP(i).x;
            coordY0 = calP(i).y;
        }
        if(elemNodesGlobal(K,1) == i){
            coordX1 = calP(i).x;
            coordY1 = calP(i).y;
        }
        if(elemNodesGlobal(K,2) == i){
            coordX2 = calP(i).x;
            coordY2 = calP(i).y;
        }
    }

    real[int, int] nodes = [[coordX0, coordY0], [coordX1,coordY1], [coordX2, coordY2]];

    border f0(t=0,1){P.x=nodes(2,0)*t + nodes(1,0)*(1-t);P.y=nodes(2,1)*t + nodes(1,1)*(1-t); label=11;};
    border f1(t=0,1){P.x=nodes(0,0)*t + nodes(2,0)*(1-t);P.y=nodes(0,1)*t + nodes(2,1)*(1-t); label=22;};
    border f2(t=0,1){P.x=nodes(1,0)*t + nodes(0,0)*(1-t);P.y=nodes(1,1)*t + nodes(0,1)*(1-t); label=33;};

    Th[K] = buildmesh(f0(1) + f1(1) + f2(1));

    int NbTrianglesLocal = Th[K].nt;
    fespace VhK(Th[K], Pk);
    int ndof = VhK.ndof;

    VhK wh, v;

    varf a(wh,v)= int2d(Th[K])( (dxx(wh)+dyy(wh)) * (dxx(v)+dyy(v)) );
    varf b(wh,v)= int2d(Th[K])( dx(wh)*dx(v) + dy(wh)*dy(v) );

    matrix A = a(VhK, VhK);
    matrix B = b(VhK, VhK);
   
    int fixedDOF = 0; // We have to fix one DOF to avoid singularity of B

    A(fixedDOF, fixedDOF) = 1e30; 
    B(fixedDOF, fixedDOF) = 1e-30; 

    for(int j = 0; j < ndof; j++) {
        if(j != fixedDOF) {
            A(fixedDOF, j) = 0;  
            A(j, fixedDOF) = 0;  
            B(fixedDOF, j) = 0;
            B(j, fixedDOF) = 0;
        }
    }


    // Using LAPACK to compute the eigenvalues and eigenvectors

    // Creating full matrices for LAPACK
    real[int,int] fullA(ndof, ndof);
    real[int,int] fullB(ndof, ndof);
    
    for(int i = 0; i < ndof; i++) {
        for(int j = 0; j < ndof; j++) {
            fullA(i,j) = A(i,j);
            fullB(i,j) = B(i,j);
        }
    }

    real[int] nev(ndof );
    complex[int] cnev(ndof );
    complex[int,int] cneV(ndof, ndof );
    dggev(fullA, fullB, cnev, nev, cneV);
   
    // Convert complex eigenvalues to real and sort them
    real[int] realEigenvalues(ndof);
    for(int i = 0; i < ndof; i++) {
        realEigenvalues[i] = real(cnev[i]);
    }

    real[int] sortedEigenvalues;
    int[int] indices = sortEigenvaluesDescending(realEigenvalues, sortedEigenvalues);

    lambdaK[][K] = sortedEigenvalues[1]; 

}
*/
//cout << "lambdaK: " << lambdaK[] << endl;

I find this implementation quite complicated, and I would like to take advantage of some of FreeFem’s built-in features to simplify and clean up the code. Is there a better or more elegant way to do this?

2

Here are a couple of things that could be improved:

  1. no need to use buildmesh, use trunc instead
  2. don’t impose penalisation by hand but use setBC instead
  3. don’t define the fespace and the varf in the for loop
  4. don’t use LAPACK directly but switch to SLEPc -eps_type lapack so you don’t have to deal with the conversion to a dense matrix yourself
3

Dear PRJ,
I’ll give it a try and get back with the changes. Thank you for your help.

4

Dear PRJ,

I’m not finding any information about setBC in the documentation.

My other problem is with the use of trunc. Do you mean that I should use it inside the for loop for each K \in \mathcal{T}_h? So I assume you’re saying I have to use trunc(Th, 1, split=1), but how can I know that this corresponds to the mesh element K \in \mathcal{T}_h?

5

You can look for examples using setBC, it’s pretty straightforward: it sets the boundary conditions.
Use trunc(Th, nuTriangle == k).

6

I’m trying some modifications that you suggested:

load "PETSc"

int coarseMesh = 3;
mesh calP = square(2^coarseMesh , 2^coarseMesh); // Structured mesh
fespace Vh(calP, P2); // FEM space for test functions and solution
fespace Rh(calP, P0); // Space for constant by triangle

//(Bi)linear forms
Vh wh, v; 
varf a(wh,v)= int2d(calP)( (dxx(wh)+dyy(wh)) * (dxx(v)+dyy(v)) );
varf b(wh,v)= int2d(calP)( dx(wh)*dx(v) + dy(wh)*dy(v) );


Rh lambdaK;

for(int K = 0; K < calP.nt; K++) {
    mesh ThK = trunc(calP , nuTriangle == K); //Extracting the K-th triangle
    fespace VhK(ThK, P2);
    int ndofK = VhK.ndof;

    matrix<real> A = a(VhK, VhK); 
    matrix<real> B = b(VhK, VhK);
    //We fix the first DoF
    VhK fixdof = 0; // zero function
    fixdof[][0] = 1; // fixing the first DoF to 1

    setBC(A, fixdof[], -2); //Applying the condition
    setBC(B, fixdof[], -2);

    Mat AA(A);
    Mat BB(B); 

    //Solve generalized eigenvalue problem
    real[int] nv;
    real[int, int] nV(1, 1);
    int nev = EPSSolve(AA, BB, sparams = "-eps_type lapack -eps_view_values -eps_view_vectors", values = nv, array = nV);

    //cout << "lambda_K:" << nv.max << endl;
    lambdaK[][K] = nv.max; // Store the maximum eigenvalue for triangle K
}

cout << "lambda_K vector" << lambdaK[] << endl;

I’m running some tests right now. The code seems fine, but I’m encountering a few issues. I’m not sure how to define the fespace VhK outside the loop. Can you please help me find the problems?

There is also and odd message when I run with coarseMesh=0

:

lambda_K vector2
         48      48
times: compile 0.108887s, execution 0.006459s,  mpirank:0
 CodeAlloc : nb ptr  4539,  size :603696 mpirank: 0
 WARNING NUMBER bad SearchMethod cas in 2d: 14 int 3d 0(essai d2: 0 3d: 0 )
Ok: Normal End

Thanks in advance, and I apologize if my questions are very basic, I’m new to FreeFEM++.

7

can you please help me find the problems?

What are the problems?

8

You mentioned that the fespace can be declared outside the for loop to improve the implementation. How can I do that?

Also, regarding the message:

lambda_K vector2
         48      48
times: compile 0.108887s, execution 0.006459s,  mpirank:0
 CodeAlloc : nb ptr  4539,  size :603696 mpirank: 0
 WARNING NUMBER bad SearchMethod cas in 2d: 14 int 3d 0(essai d2: 0 3d: 0 )
Ok: Normal End

Should I be worried about that?

9

How can I do that?

By moving the declaration outside of the loop?

Should I be worried about that?

No.

10

Ok, thanks for your help.