# Generalized eigenvalue problem in a mesh

**URL:** <https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022>\
**Category:** General Discussion\
**Created:** [July 27, 2025, 7:02pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022 "2025-07-27T19:02:18Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 27, 2025, 7:02pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/1 "2025-07-27T19:02:18Z")

</div>

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:

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 28, 2025, 6:35am UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/2 "2025-07-28T06:35:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 28, 2025, 12:12pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/3 "2025-07-28T12:12:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 28, 2025, 2:01pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/4 "2025-07-28T14:01:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 28, 2025, 2:22pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/5 "2025-07-28T14:22:18Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 28, 2025, 3:34pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/6 "2025-07-28T15:34:57Z")

</div>

I’m trying some modifications that you suggested:

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

:

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 29, 2025, 5:06am UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/7 "2025-07-29T05:06:10Z")

</div>

> can you please help me find the problems?

What are the problems?

---

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 29, 2025, 12:15pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/8 "2025-07-29T12:15:43Z")

</div>

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:

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [July 29, 2025, 12:23pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/9 "2025-07-29T12:23:10Z")

</div>

> How can I do that?

By moving the declaration outside of the loop?

> Should I be worried about that?

No.

---

<div class="post-metadata">

**Author:** ![juanpaca](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/juanpaca/32/2807_2.png) [@juanpaca](https://community.freefem.org/u/juanpaca)\
**Post date:** [July 29, 2025, 12:28pm UTC](https://community.freefem.org/t/generalized-eigenvalue-problem-in-a-mesh/4022/10 "2025-07-29T12:28:58Z")

</div>

Ok, thanks for your help.