Error encountered when constructing the stiffness matrix with PETSc

Why do I encounter an error at A = laplace(W1, W1); when executing the following PETSc parallel code, but the problem disappears when I remove rhoe? Additionally, I experience the same issue when I set rhoe to 1. I’m not sure why this is happening.

include “set_material.edp”

        varf laplace(phi,phiv, eps=1.0e-8)=
            int3d(Th)( Grad(phi)'*Grad(phiv)/rhoe )
            +on(labHead,phi=voltage/2)
            +on(labTail,phi=-voltage/2)
            +on(labBx0,phi=0)
            +on(labBx1,phi=0)
            +on(labBz0,phi=0)
            +on(labBz1,phi=0)
            +on(labBy0,phi=0)
            +on(labBy1,phi=0);

        real time01 = clock();
        Mat A ;
        // MatCreate(Th,A,P13d);
        buildMat(Th, getARGV("-split", 1), A, P13d, mpiCommWorld)  
        real time02 = clock();
        A= laplace(W1, W1); 
        real[int] b = laplace(0, W1); 

        real time03 = clock();
        set(A, sparams = "-pc_type gamg -ksp_type gmres");
        real time04 = clock();
        phi[] = A^-1 * b;

Isn’t it simply that rhoe=0? Try rhoe=1; befor you construc the matrix.

I have tried setting rhoe = 1, but I still get the same error.

What is the error? Please share an actual usable code.

test16.edp (2.3 KB)
When I use a single core, it runs normally, but when I use multiple cores, it errors out.

I have already figured out where this issue is coming from, thank you for your answer. I apologize, but I have another question: how can I create a mesh that cuts a smaller cylinder out of a larger cylinder, like in this picture? I am looking forward to your reply.

Use Gmsh, there are many ressources available online.

Thank you very much for your answers, and I will also take a close look at the Gmsh software. In the question two above, when I changed the code from mesh3 Th4 = cube(N/2, N/3, N, [x, 2+y, z], label = LL4, region=4); to mesh3 Th4 = cube(N, N/3, N, [x, 2+y, z], label = LL4, region=4);, it became possible to run normally. Could you explain why this is the case?

If I have a slightly more complex mesh, like this, is there any convenient way to solve this problem? Looking forward to your reply.
test16.edp (2.5 KB)

What is the problem?

When I run this code, the following error occurs:
Exec error : Div by 0
– number :1
Exec error : Div by 0
– number :1
err code 8 , mpirank 0

This is only part of the error message. The line number should be print just before, so that you can look back at your file and fix the division by 0.