Different result of the same problem using different meshes

General Discussion
1

meshEDDY.zip (3.5 MB)
eddyGMSH.zip (2.4 MB)

Hi everybody!
I have developed in the past times a code that solves eddy current problems with freefem++ (MUMPS solver). Solving for meshes directly developed in FF++ I can get results that can be considered correct (I have benchmarked them against other results), by the way passing to meshes built in gmsh it becomes impossible to recover correct results (I mean on the very same problem and geometry, changing only the mesher!!!)

Is somebody aware of which could be the reason? I am using Edge03d elements to solve the problem, are there any compatibility issues between edge elements and meshes generated in gmsh?

For sake of being clearer I’m attaching the code I’m using with the adopted meshes.

I want already to thank a lot all the people will help me in trying to solve this issue!

yours,

Marco

P.S. if you comment line 6 and de-comment line 7 you get a working code.

    load "msh3"
load "tetgen"

include "MeshSurface.idp"

mesh3 Th("eddyGMSH.mesh");
//mesh3 Th("meshEDDY.mesh");

plot(Th, cmm = "Volume mesh from gmsh", wait = 1);

cout<<"number of cells = "<<Th.nt<<endl;
cout<<"Mesh regions = "<<regions(Th)<<endl;
cout<<"mesh labels = "<<labels(Th)<<endl;

int cable = Th(0.0,0.0,1.0).region;
cout<<"Cable identifier = "<<cable<<endl;
int jacket = Th(0.0,0.0,0.89).region;
cout<<"Jacket identifier = "<<jacket<<endl;
int out = Th(0.0,0.0,2.5).region;
cout<<"Air identifier = "<<out<<endl;
real vc = int3d(Th,cable)(1.);
cout<<"cable volume = "<<vc<<" [m^3]"<<endl;
real vj = int3d(Th,jacket)(1.);
cout<<"jacket volume = "<<vj<<" [m^3]"<<endl;

///////////////////////////////////////////////
//////////// DEFINITION OF FESPACES ///////////
///////////////////////////////////////////////

func Pk = Edge03d;
func Ps = P0;
func Pj = [P1dc,P1dc,P1dc];

fespace Vh(Th,Pk);
fespace Vhp(Th,Ps);
fespace Vhj(Th,Pj);

///////////////////////////////////////////////
/// DEFINITION OF SOME RELEVANT PROPERTIES ///
//////////////////////////////////////////////

real tempo;
real J0 = 1.0e4;
func Jmod = J0*tempo*(tempo<=1.0)+J0*(tempo>1.0 && tempo<=1.5)+J0*exp(-(tempo-1.5))*(tempo>1.5);
real mu = 4*pi*1.0e-7;
real nu = 1.0/mu;
Vhp sigma = 1.0*(region == cable)+1.12e7*(region == jacket)+1.0*(region == out);

/*
savemesh(Th,"permatlab.mesh");
ffSaveVh(Th,Vhp,"Vhpermatlab.txt");
ffSaveData(sigma,"sigmapermatlab.txt");
*/


//////////////////////////////////////////////
//////////// TRANSIENT INFORMATION ///////////
/////////////////////////////////////////////

Vh [A0x,A0y,A0z] = [0.,0.,0.];
real dt; //defined... the value will be discussed later on
real tend = 2.5;
Vhj [Jx,Jy,Jz];

/////////////////////////////////////////////
/////////////////// MACROS //////////////////
////////////////////////////////////////////

macro Curl(Ax,Ay,Az) [dy(Az)-dz(Ay), dz(Ax)-dx(Az), dx(Ay)-dy(Ax)]//EOM

/////////////////////////////////////////////
///////////// PROBLEM DEFINITION ////////////
////////////////////////////////////////////

varf vAmp([Ax,Ay,Az],[Wx,Wy,Wz]) = int3d(Th)(nu*Curl(Ax,Ay,Az)'*Curl(Wx,Wy,Wz))+int3d(Th)(sigma*[Ax,Ay,Az]'*[Wx,Wy,Wz]/dt);
varf vAmprhs([Ax,Ay,Az],[Wx,Wy,Wz]) = int3d(Th)(sigma*[A0x,A0y,A0z]'*[Wx,Wy,Wz]/dt)+int3d(Th)([Jx,Jy,Jz]'*[Wx,Wy,Wz]);

matrix AA;
real[int] rhs(Vh.ndof);

tempo = 0.0;

///////////////////////////////////////////
///////////////// TRANSIENT ///////////////
//////////////////////////////////////////
ofstream fa("AtimeGMSHandTETGEN.dat");

while (tempo<=tend)
{
  // now need to impose the value of dt:
  if (tempo<1.5)
  {
    dt = 0.1;
  }
  else
  {
    dt = 0.01;
  }

  tempo = tempo + dt;
  [Jx,Jy,Jz] = [Jmod*(z/sqrt(x^2+z^2))*(region==cable),0.,Jmod*(-x/sqrt(x^2+z^2))*(region==cable)];

  AA = vAmp(Vh,Vh);
  rhs = vAmprhs(0,Vh);
  set(AA,solver = sparsesolver);
  Vh [Ax,Ay,Az];
  Ax[] = AA^-1*rhs;

  {
    ofstream fa("AtimeGMSHandTETGEN.dat",append);
    fa<<tempo<<" "<<Ax(0.0,0.0,0.88)<<" "<<Ay(0.0,0.0,0.88)<<" "<<Az(0.0,0.0,0.88)<<" "<<Ax(0.0,0.0,0.92)<<" "<<Ay(0.0,0.0,0.92)<<" "<<Az(0.0,0.0,0.92)<<endl;
  }

  cout<<"Done at time = "<<tempo<<" [s]"<<endl;

  [A0x,A0y,A0z] = [Ax,Ay,Az];

}