# Bug with Petsc in building a local matrix with edge elements 

**URL:** <https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706>\
**Category:** General Discussion\
**Created:** [January 23, 2025, 10:06am UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706 "2025-01-23T10:06:51Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ctain](https://avatars.discourse-cdn.com/v4/letter/c/dbc845/32.png) [@ctain](https://community.freefem.org/u/ctain)\
**Post date:** [January 23, 2025, 10:06am UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706/1 "2025-01-23T10:06:51Z")

</div>

Dear all,

I encounter a bug when assembling a matrix with Nedelec elements. The problem does not show up with 4 procs but with 8 and more.

Thank you for your help

//=================================================  
load “msh3”  
load “medit”  
load “PETSc-complex”  
//load “PETSc”  
load “mmg”  
load “tetgen”  
include “MeshSurface.idp”  
macro dimension()3// EOM // 2D or 3D  
include “macro\_ddm.idp” // additional DDM functions

//==================================================  
// parameters (physical, numerics)  
//==================================================  
int R = 5;  
real kappa = 4.;  
real nodePerXi = 1;  
int nbsegul = kappa\*nodePerXi;  
real hsize= 1./nbsegul;

// applied field  
real Bax = 0.;  
real Bay = 0.;  
real Baz = 0.9;

real omega = 1.;  
tgv = 1e30;  
real dt = 0.1;  
real idt=1./dt;

//==================================================  
// files for savings  
//==================================================

ofstream fout(“./a\_echo”);  
fout.precision(12);

//==================================================  
// macros  
//==================================================  
macro div(A)(dx(A) + dy(A#y) + dz(A#z))//EOM  
macro Grad(u)([dx(u), dy(u), dz(u)])//

//==================================================  
//finite element space  
//==================================================

func PkMagEdge = Edge03d;  
func PkMagRT = RT03d;

//===============================================================  
// global mesh  
//===============================================================

//---------- sphere de référence  
meshS Thsph0 = Sphere(R,hsize,1,1);  
real[int] domain = [0.,0.,0.,1,hsize^3/6];  
mesh3 Th=tetg(Thsph0,switch=“paAAQYY”,nbofregions=1,regionlist=domain);  
if(mpirank==0) medit(“sphere”, Th);

// if(mpirank==0)  
// {  
// cout \<\< " ========== sphere.hmin = " \<\< Thsph.hmin \<\< endl;  
// cout \<\< " ========== sphere.hmax = " \<\< Thsph.hmax \<\< endl;  
// cout \<\< " ========== sphere.nt = " \<\< Thsph.nt \<\< endl;  
// fout \<\< " ========== sphere.hmin = " \<\< Thsph.hmin \<\< endl;  
// fout \<\< " ========== sphere.hmax = " \<\< Thsph.hmax \<\< endl;  
// fout \<\< " ========== sphere.nt = " \<\< Thsph.nt \<\< endl;  
// }

//===============================================================  
//local FE spaces  
//===============================================================

buildDmesh(Th)

//local FE spaces  
fespace VhMagEdge(Th, PkMagEdge);  
fespace VhMagRT(Th, PkMagRT);

if(mpirank==0)  
{  
cout \<\< " \*\*\*\*\*\*\*\*\*\* local mesh on proc = " \<\< mpirank \<\< endl;  
cout \<\< " ========== ntetra = " \<\< Th.nt \<\< endl;  
cout \<\< " ========== dof magField = " \<\< VhMagEdge.ndof \<\< endl;  
cout \<\< " ========== dof vectPot = " \<\< VhMagRT.ndof \<\< endl;

```
fout << " ********** local mesh on proc = " << mpirank << endl;
fout << " ========== ntetra = " << Th.nt << endl;
fout << " ========== dof magField = " << VhMagEdge.ndof << endl;
fout << " ========== dof vectPot = " << VhMagRT.ndof << endl;

```

}

//===============================================================  
//local matrices  
//local weak forms  
//===============================================================

//local matrices  
Mat Asys11, Asys22;

{  
macro def(i)[i, i#y, i#z]//EOM // vector field definition  
macro init(i)[i, i, i]//EOM // vector field initialization  
createMat(Th, Asys11, PkMagEdge);

// macro necessaire pour la construction de Asys22 (bug PETSc)  
macro ThPostProcessD(D) {  
VhMagRT def(u), def(v);  
varf onG(def(u), def(v)) = on(-111111, u = 10 + x - z, uy = 100 + y - x, uz = 1000 + z - y);  
v = D;  
u = onG(0, VhMagRT);  
for [j, dj : u] dj = abs(dj) \> 1e-2 ? 0.0 : 1.0;  
def(u) = [u, uy, uz];  
D = u;  
} //EOM  
createMat(Th, Asys22, PkMagRT);  
}

// local weak forms

varf vSys11([sigma, sigmay, sigmaz],[ki, kiy, kiz]) =  
int3d(Th, qfV = qfV2)([sigma, sigmay, sigmaz]'\*[ki, kiy, kiz])

- on(1, sigma = Bax, sigmay = Bay, sigmaz = Baz);

varf vSys22([A, Ay, Az],[v, vy, vz]) =  
int3d(Th, qfV = qfV2)(idt\*[A, Ay, Az]'\*[v, vy, vz])

- int3d(Th, qfV = qfV2)(omega\*div(A)\*div(v))
- on(1, A = 0., Ay = 0., Az = 0.); // [A, Ay, Az] = 0 ==\> A.n = 0 here

// local matrices  
Asys11 = vSys11(VhMagEdge, VhMagEdge); // bug ici avec 16 proc mais 4 proc ok  
//Asys22 = vSys22(VhMagRT, VhMagRT);

---

<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:** [January 23, 2025, 10:07am UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706/2 "2025-01-23T10:07:21Z")

</div>

Can’t copy/paste your code =\> can’t help you.

---

<div class="post-metadata">

**Author:** ![ctain](https://avatars.discourse-cdn.com/v4/letter/c/dbc845/32.png) [@ctain](https://community.freefem.org/u/ctain)\
**Post date:** [January 23, 2025, 2:20pm UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706/3 "2025-01-23T14:20:18Z")

</div>

Here is the code:

> **(attachments)**
>
> [Zcode.edp](https://community.freefem.org/uploads/short-url/wK6N63GFYNZpmUNqyCuUCFmcgA1.edp) (4.04 KB)

---

<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:** [January 24, 2025, 11:39am UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706/4 "2025-01-24T11:39:05Z")

</div>

I see two wrong things:

1. you should broadcast your mesh from processor 0 after the mesh adaptation step because there are no guarantees that the mesh will be the same among all processes
2. for edge elements, you should not define `init` and `part` but instead `PkPart`, see [FreeFem-sources/examples/hpddm/maxwell-3d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/maxwell-3d-PETSc.edp#L29)

With both fixes, the code runs to completion:

```auto
$ ff-mpirun -n 16 Zcode.edp -v 0 -version
[...]
********** local mesh on proc = 0
 ========== ntetra = 29545
 ========== dof magField = 37680
 ========== dof vectPot = 61270
WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 2 unused database options. They are:
Option left: name:-nw value: Zcode.edp source: command line
Option left: name:-v value: 0 source: command line

```

---

<div class="post-metadata">

**Author:** ![ctain](https://avatars.discourse-cdn.com/v4/letter/c/dbc845/32.png) [@ctain](https://community.freefem.org/u/ctain)\
**Post date:** [January 24, 2025, 3:21pm UTC](https://community.freefem.org/t/bug-with-petsc-in-building-a-local-matrix-with-edge-elements/3706/5 "2025-01-24T15:21:43Z")

</div>

Thank you very much for your help