Error with restrict() and loaded mesh

General Discussion
1

Hello,

I am having problems using restrict() with a mesh loaded with loadDmesh(). Here is a MWE that returns an error when either of the last 2 lines is not commented.
Could you please tell me if I am doing something wrong, or if there is a problem with my installation?

load “PETSc”
int[int] n2oSaved;
int[int] n2oLoaded;
{
macro dimension()2//
include “macro_ddm.idp”
macro def(i)[i, i#B]//
macro init(i)[i, i]//
int[int] n2oTh;
macro ThN2O()n2oTh//
mesh Th = square(40, 40);
mesh ThGlobal = Th;
Mat A;
createMat(Th, A, [P1b,P1b])
saveDmesh(Th, “dump-2d”)
savemesh(ThGlobal, “dump-2d-G.mesh”);
fespace XhC(Th, P1b);
fespace VhC(Th, [P1b,P1b]);
fespace XhCGlob(ThGlobal, P1b);
fespace VhCGlob(ThGlobal, [P1b,P1b]);
int[int] restX = restrict(XhC, XhCGlob, n2oTh); // Ok
int[int] restV = restrict(VhC, VhCGlob, n2oTh); // Ok
}
{
macro dimension()2//
include “macro_ddm.idp”
macro def(i)[i, i#B]//
macro init(i)[i, i]//
int[int] n2oTh;
macro ThN2O()n2oTh//
mesh Th;
fespace XhC(Th, P1b);
fespace VhC(Th, [P1b,P1b]);
loadDmesh(Th, “dump-2d”)
mesh ThGlobal;
ThGlobal = Th;
fespace XhCGlob(ThGlobal, P1b);
fespace VhCGlob(ThGlobal, [P1b,P1b]);
Mat A;
createMat(Th, A, [P1b,P1b])
//int[int] restX = restrict(XhC, XhCGlob, n2oTh); // Error here
//int[int] restV = restrict(VhC, VhCGlob, n2oTh); // Error here
}

The error message looks like this:

Assertion fail : (kf >= 0 && kf <> neF)
line :415, in file lgmat.cpp

Thank you!

2

Please make your code runnable when pasting a MWE. Your snippet needs to be edited, it cannot be run directly (issue with the double quotes). Anyway, the fix is:

--- original.edp	2024-07-06 07:03:39
+++ fixed.edp	2024-07-06 07:02:39
@@ -13,7 +13,7 @@ savemesh(ThGlobal, "dump-2d-G.mesh");
 Mat A;
 createMat(Th, A, [P1b,P1b])
 saveDmesh(Th, "dump-2d")
-savemesh(ThGlobal, "dump-2d-G.mesh");
+savemesh(ThGlobal, "dump-2d-G.msh");
 fespace XhC(Th, P1b);
 fespace VhC(Th, [P1b,P1b]);
 fespace XhCGlob(ThGlobal, P1b);
@@ -32,12 +32,11 @@ mesh ThGlobal;
 fespace XhC(Th, P1b);
 fespace VhC(Th, [P1b,P1b]);
 loadDmesh(Th, "dump-2d")
-mesh ThGlobal;
-ThGlobal = Th;
+mesh ThGlobal = readmesh("dump-2d-G.msh");
 fespace XhCGlob(ThGlobal, P1b);
 fespace VhCGlob(ThGlobal, [P1b,P1b]);
 Mat A;
 createMat(Th, A, [P1b,P1b])
-//int[int] restX = restrict(XhC, XhCGlob, n2oTh); // Error here
-//int[int] restV = restrict(VhC, VhCGlob, n2oTh); // Error here
+int[int] restX = restrict(XhC, XhCGlob, n2oTh); // Error here
+int[int] restV = restrict(VhC, VhCGlob, n2oTh); // Error here
 }
3

Thank you very much.