Level-set with parmmg3d

General Discussion
1

Dear FreeFem users,

I am trying to use parmmg3d from FreeFEM in level-set mode, but I am not sure whether this functionality is actually supported by the current FreeFEM interface to ParMmg.

From reading the parmmg.cpp plugin source in FreeFEM, it seems that the parameters

  • iso
  • ls

do exist and are forwarded to the ParMmg API. However, in practice, the call does not seem to modify the mesh. So, is ParMmg level-set mode really supported through parmmg3d in FreeFEM?

In a minimal test, if I run ParMmg externally from the shell with something like

parmmg_O3 -in mesh.mesh -sol phi.sol -ls 0 -out out.mesh

the level-set mode seems to be available.

But when I call it from FreeFEM with parmmg3d(..., metric=..., ls=0), the mesh does not change, and I also get the message

## Error: MMG5_scale_scalarMetric: at least 1 wrong metric.

Here is a minimal example:

load "mshmet"
load "PETSc"
load "parmmg"
load "medit"
macro dimension()3// EOM
include "macro_ddm.idp"

func real Heav(real x, real beta, real eta)
{
  return (tanh(beta*eta) + tanh(beta*(x-eta))) /
         (tanh(beta*eta) + tanh(beta*(1.0-eta)));
}

int[int] fforder = [1];
real err = 0.05;
int n = 15;

mesh3 Th3=cube(n,n,n);
mesh3 ThBackup = Th3;

Mat A;
int[int] n2o;
NewMacro Th3N2O() n2o EndMacro
DmeshCreate(Th3);

MatCreate(Th3, A, P1);
fespace Vh(Th3, P1);
Vh u;
plotDmesh(Th3, cmm = "Th3");

Vh uh = sqrt((x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2) - 0.3;
// Smooth
uh[] += 0.5;
uh = Heav(uh, 40, 0.5) - 0.5;

fespace VhBackup(ThBackup, P1);
VhBackup h, hReduced;
uh[] .*= A.D;
int[int] rest = restrict(Vh, VhBackup, n2o);
for[i, v : rest] h[][v] = uh[][i];
mpiAllReduce(h[], hReduced[], mpiCommWorld, mpiSUM); // centralize distributed solution

// real[int] met = mshmet(ThBackup, hReduced, aniso = 0, hmin = 1.0e-3, hmax = 2.0e-1, err = err);
ThBackup = parmmg3d(ThBackup, metric = uh[], hgrad = 2.3, verbose=3, ls=0);
Th3 = ThBackup;
Mat Adapt;
MatCreate(Th3, Adapt, P1);
A = Adapt;
cout << "mpirank # " << mpirank << " : Th3.nv " << Th3.nv << ", Th3.nt " << Th3.nt << endl;
plotDmesh(Th3, cmm = "Th3Adapt");

Thank you