Gmsh - errror element type

ftrillaudp · October 26, 2022, 12:08pm

Dear all,

I am starting using FreeFem++ benchmarking some code. I have created a 3D ring reusing the magnetostatic example provided in FreeFem++ website. This example in 3D does not work producing the mesh with the lastest version of Gmsh 4.11. The software complains about the element type not available:

“Element of type 12 is not considered in Freefem++”

I am joining the gmsh *.geo and the *edp.

Best,

Frederic

GEO file:

DefineConstant[
coilInnerRadius = {0.104, Name “Input/0Geometry/0Coil inner radius (m)”},
coilWidth = {0.016, Name “Input/0Geometry/1Coil width (m)”},
coilThickness = {0.046, Name “Input/0Geometry/2Coil thickness (m)”}
];
airRadius = 2*(coilInnerRadius+coilThickness);

airMeshSize = 4;
coilMeshSize = 4;
lc_1 = coilWidth / coilMeshSize;
lc_2 = airRadius / airMeshSize;

// — Constant definition for regions —
coilID = 1;
airID = 2;
airBoundaryID = 3;

/// — Geometry and Meshing ----
SetFactory(“OpenCASCADE”);
General.ExpertMode = 1;

/// Geometry options:
Geometry.Color.Lines = Black;
Geometry.LineWidth = 1;
Geometry.PointSize = 1;

/// Optimize mesh:
Mesh.Optimize = 1;

// Save mesh as binary file (to save storage): cannot be used with ElmerGrid, choose 0
Mesh.Binary = 0;

/// — Build geometry —
Cylinder(1) = {0., 0., (-1)0.5coilThickness, 0., 0., coilThickness, coilInnerRadius};
Cylinder(2) = {0., 0., (-1)0.5coilThickness, 0., 0., coilThickness, coilInnerRadius+coilWidth};
Sphere(3) = {0., 0., 0., airRadius};

BooleanDifference (4) = { Volume{2}; Delete ;}{ Volume{1}; Delete ;};
BooleanDifference (5) = { Volume{3}; Delete ;}{ Volume{4}; };

surfaces_1() = Boundary{Volume{4};};
MeshSize{ PointsOf{ Volume{4}; } } = lc_1;
surfaces_2() = Boundary{Volume{5};};
MeshSize{ PointsOf{ Surface{surfaces_2(0)}; } } = lc_2;

/// Remove duplicate geometrical elementary entities:
Coherence;

/// Remove duplicate nodes of the mesh:
Coherence Mesh;

// Definition of volumes and boundaries //
Physical Volume(“coil”, coilID) = {4};
Color Red {Volume {4};}
Physical Volume(“air”, airID) = {5};
Color Blue {Volume {5};}
// For ElmerGrid, we cannot have the same submesh in different Physical Surface.
Physical Surface(“airBoundary”, airBoundaryID) = {surfaces_2(0)};

/// — Miscellaneous information —
cpuTime = Cpu;
memoryUsage = Memory;
Printf(“CPU time”, cpuTime);
Printf(“Memory usage”, memoryUsage);

EDP file:

load “gmsh”

//Parameters
real Mu0 = 4.pi1.e-7; //Vacuum magnetic permeability
real MuC = 1.25e-6; //Copper magnetic permeability
real J0 = 500.e4; //Current density

//Mesh
real R = 2.;
real H = 0.046;
real Rint = 0.104;
real Rext = 0.104+0.016;

int Coil = 10;
int BoxWall = 1;

mesh3 Th = gmshload3(“benchmark1.msh”);

//FESpaces
func Pk = P2;
fespace Ah(Th, [Pk, Pk, Pk]);
Ah [Ax, Ay, Az] = [0, 0, 0];

//Macro
macro Curl(Ax, Ay, Az) [dy(Az)-dz(Ay), dz(Ax)-dx(Az), dx(Ay)-dy(Ax)] //
macro Divergence(Ax, Ay, Az) (dx(Ax) + dy(Ay) + dz(Az)) //

//Current
func r = sqrt(x^2+y^2);
Ah [Jx, Jy, Jz] = J0*[
cos(atan2(y, x)+pi/2.) * (r >= Rint)(r <= Rext)(z >= -H/2.)(z <= H/2.),
0,
sin(atan2(y, x)+pi/2.) * (r >= Rint)(r <= Rext)(z >= -H/2.)(z <= H/2.)
];

//Problem
func Mu = Mu0 + (MuC-Mu0)*(region==Coil);
func Nu = 1. / Mu;
varf vLaplacian ([Ax, Ay, Az], [AAx, AAy, AAz])
= int3d(Th)(
Nu * Curl(Ax, Ay, Az)’ * Curl(AAx, AAy, AAz)
+ (1./Mu) * Divergence(Ax, Ay, Az) * Divergence(AAx, AAy, AAz)
)
- int3d(Th)(
[Jx, Jy, Jz]’ * [AAx, AAy, AAz]
)
+ on(BoxWall, Ax=0, Ay=0, Az=0)
;

matrix Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver);
real[int] LaplacianBoundary = vLaplacian(0, Ah);
Ax = Laplacian^-1 * LaplacianBoundary;

//Magnetic induction
Ah [Bx, By, Bz];
[Bx, By, Bz] = Curl(Ax, Ay, Az);

//Magnetic field
Ah [Hx, Hy, Hz];
[Hx, Hy, Hz] = Nu * [Bx, By, Bz];

prj · October 26, 2022, 12:49pm

What is the format of benchmark1.msh (ASCII or binary, Gmsh format version)?

ftrillaudp · October 26, 2022, 2:00pm

It is ASCII. I thought that it was necessary. But you are right the error disappear using binary format for the mesh.

I tried with the basic magnetostatic example but it does not seem to run in 3D. IN 2D it runs fine.

I am running gmsh 4.11 and FreeFem++ 4.11 on ubuntu 20.04LTS.

I have attached the code downloaded from the magnetostatic example for a try. The geo file is renamed benchmark1.geo to create the mesh in Gmsh.

Best,

Frederic

load "gmsh"

//Parameters
real Mu0 = 4.*pi*1.e-7;	//Vacuum magnetic permeability
real MuC = 1.25e-6;		//Copper magnetic permeability
real J0 = 500.e4;		//Current density

//Mesh
real R = 2.;
real H = 1.;
real Rint = 0.1;
real Rext = 0.2;

int Coil = 10;
int BoxWall = 1;

mesh3 Th = gmshload3("benchmark1.msh");

//FESpaces
func Pk = P2;
fespace Ah(Th, [Pk, Pk, Pk]);
Ah [Ax, Ay, Az] = [0, 0, 0];

//Macro
macro Curl(Ax, Ay, Az) [dy(Az)-dz(Ay), dz(Ax)-dx(Az), dx(Ay)-dy(Ax)] //
macro Divergence(Ax, Ay, Az) (dx(Ax) + dy(Ay) + dz(Az)) //

//Current
func r = sqrt(x^2+y^2);
Ah [Jx, Jy, Jz] = J0*[
	cos(atan2(y, x)+pi/2.) * (r >= Rint)*(r <= Rext)*(z >= -H/2.)*(z <= H/2.),
	0,
	sin(atan2(y, x)+pi/2.) * (r >= Rint)*(r <= Rext)*(z >= -H/2.)*(z <= H/2.)
	];

//Problem
func Mu = Mu0 + (MuC-Mu0)*(region==Coil);
func Nu = 1. / Mu;
varf vLaplacian ([Ax, Ay, Az], [AAx, AAy, AAz])
	= int3d(Th)(
		  Nu * Curl(Ax, Ay, Az)' * Curl(AAx, AAy, AAz)
		+ (1./Mu) * Divergence(Ax, Ay, Az) * Divergence(AAx, AAy, AAz)
	)
	- int3d(Th)(
		  [Jx, Jy, Jz]' * [AAx, AAy, AAz]
	)
	+ on(BoxWall, Ax=0, Ay=0, Az=0)
	;

matrix<real> Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver);
real[int] LaplacianBoundary = vLaplacian(0, Ah);
Ax[] = Laplacian^-1 * LaplacianBoundary;

//Magnetic induction
Ah [Bx, By, Bz];
[Bx, By, Bz] = Curl(Ax, Ay, Az);

//Magnetic field
Ah [Hx, Hy, Hz];
[Hx, Hy, Hz] = Nu * [Bx, By, Bz];

plot(Bz);

gmsh:

//////////////////
///Optimization///
//////////////////
Mesh.Optimize = 1;
Mesh.Binary = 1;

////////////////
///Parameters///	include in the Gmsh GUI
////////////////
DefineConstant
[
	R = {2., Name "Geometry/01 «Infinite» boundary radius"},
	LM = {1., Name "Geometry/02 Magnet length"},
	RMext = {0.2, Name "Geometry/03 Magnet external radius"},
	RMint = {0.1, Name "Geometry/04 Magnet internal radius"},
	h = {1., Name "Mesh/01 Characteristic length for boundary"},
	hM = {0.05, Name "Mesh/02 Characteristic length for magnet"}
];

//////////////////////
///Global variables///
//////////////////////
Wall[] = {};
MagnetUp[] = {};
MagnetDown[] = {};
MagnetMid[] = {};

///////////////////
///Global sphere///
///////////////////
p = newp;
Point(p+0) = {0, 0, 0, h};
Point(p+1) = {R, 0, 0, h};
Point(p+2) = {0, R, 0, h};
Point(p+3) = {-R, 0, 0, h};
Point(p+4) = {0, -R, 0, h};

l = newl;
Circle(l+0) = {p+1, p+0, p+2};
Circle(l+1) = {p+2, p+0, p+3};
Circle(l+2) = {p+3, p+0, p+4};
Circle(l+3) = {p+4, p+0, p+1};

D1[] = Duplicata {
	Line{l+0, l+1, l+2, l+3};
};
//Printf("D1 :", D1[0], D1[1], D1[2], D1[3]);

Rotate { {1, 0, 0}, {0, 0, 0}, Pi/2.} {
	Line{D1[0], D1[1], D1[2], D1[3]};
}

D2[] = Duplicata {
	Line{l+0, l+1, l+2, l+3};
};
//Printf("D2 :", D2[0], D2[1], D2[2], D2[3]);

Rotate { {0, 1, 0}, {0, 0, 0}, Pi/2.} {
	Line{D2[0], D2[1], D2[2], D2[3]};
}

ll =newll;
Line Loop(ll+0) = {l+0, D2[1], -D1[0]};
Line Loop(ll+1) = {l+0, D1[3], -D2[0]};
Line Loop(ll+2) = {l+1, D1[2], D2[0]};
Line Loop(ll+3) = {l+1, -D1[1], -D2[1]};
Line Loop(ll+4) = {l+2, D2[3], -D1[2]};
Line Loop(ll+5) = {l+2, -D2[2], D1[1]};
Line Loop(ll+6) = {l+3, -D1[3], -D2[3]};
Line Loop(ll+7) = {l+3, D1[0], D2[2]};


s = news;
Ruled Surface(s+0) = {ll+0};	Wall[0] = s+0;
Ruled Surface(s+1) = {ll+1};	Wall[1] = s+1;
Ruled Surface(s+2) = {ll+2};	Wall[2] = s+2;
Ruled Surface(s+3) = {ll+3};	Wall[3] = s+3;
Ruled Surface(s+4) = {ll+4};	Wall[4] = s+4;
Ruled Surface(s+5) = {ll+5};	Wall[5] = s+5;
Ruled Surface(s+6) = {ll+6};	Wall[6] = s+6;
Ruled Surface(s+7) = {ll+7};	Wall[7] = s+7;

////////////
///Magnet///
////////////
pext = newp;
Point(pext+0) = {0, 0, -LM/2., hM};
Point(pext+1) = {RMext, 0, -LM/2., hM};
Point(pext+2) = {0, RMext, -LM/2., hM};
Point(pext+3) = {-RMext, 0, -LM/2., hM};
Point(pext+4) = {0, -RMext, -LM/2., hM};

Point(pext+5) = {0, 0, LM/2., hM};
Point(pext+6) = {RMext, 0, LM/2., hM};
Point(pext+7) = {0, RMext, LM/2., hM};
Point(pext+8) = {-RMext, 0, LM/2., hM};
Point(pext+9) = {0, -RMext, LM/2., hM};

lext = newl;
Circle(lext+0) = {pext+1, pext+0, pext+2};
Circle(lext+1) = {pext+2, pext+0, pext+3};
Circle(lext+2) = {pext+3, pext+0, pext+4};
Circle(lext+3) = {pext+4, pext+0, pext+1};

Circle(lext+4) = {pext+6, pext+5, pext+7};
Circle(lext+5) = {pext+7, pext+5, pext+8};
Circle(lext+6) = {pext+8, pext+5, pext+9};
Circle(lext+7) = {pext+9, pext+5, pext+6};

Line(lext+8) = {pext+1, pext+6};
Line(lext+9) = {pext+2, pext+7};
Line(lext+10) = {pext+3, pext+8};
Line(lext+11) = {pext+4, pext+9};

llext = newll;
Line Loop(llext+0) = {lext+0, lext+9, -(lext+4), -(lext+8)};
Line Loop(llext+1) = {lext+1, lext+10 ,-(lext+5), -(lext+9)};
Line Loop(llext+2) = {lext+2, lext+11, -(lext+6), -(lext+10)};
Line Loop(llext+3) = {lext+3, lext+8, -(lext+7), -(lext+11)};

sext = news;
Ruled Surface(sext+0) = {llext+0};	MagnetMid[0] = sext+0;
Ruled Surface(sext+1) = {llext+1};	MagnetMid[1] = sext+1;
Ruled Surface(sext+2) = {llext+2};	MagnetMid[2] = sext+2;
Ruled Surface(sext+3) = {llext+3};	MagnetMid[3] = sext+3;

Rotate { {0, 1, 0}, {0, 0, 0}, Pi/2. } {
	Surface{sext+0, sext+1, sext+2, sext+3};
}
Rotate { {0, 0, 1}, {0, 0, 0}, Pi/2. } {
	Surface{sext+0, sext+1, sext+2, sext+3};
}

pint = newp;
Point(pint+0) = {0, 0, -LM/2., hM};
Point(pint+1) = {RMint, 0, -LM/2., hM};
Point(pint+2) = {0, RMint, -LM/2., hM};
Point(pint+3) = {-RMint, 0, -LM/2., hM};
Point(pint+4) = {0, -RMint, -LM/2., hM};

Point(pint+5) = {0, 0, LM/2., hM};
Point(pint+6) = {RMint, 0, LM/2., hM};
Point(pint+7) = {0, RMint, LM/2., hM};
Point(pint+8) = {-RMint, 0, LM/2., hM};
Point(pint+9) = {0, -RMint, LM/2., hM};

lint = newl;
Circle(lint+0) = {pint+1, pint+0, pint+2};
Circle(lint+1) = {pint+2, pint+0, pint+3};
Circle(lint+2) = {pint+3, pint+0, pint+4};
Circle(lint+3) = {pint+4, pint+0, pint+1};

Circle(lint+4) = {pint+6, pint+5, pint+7};
Circle(lint+5) = {pint+7, pint+5, pint+8};
Circle(lint+6) = {pint+8, pint+5, pint+9};
Circle(lint+7) = {pint+9, pint+5, pint+6};

Line(lint+8) = {pint+1, pint+6};
Line(lint+9) = {pint+2, pint+7};
Line(lint+10) = {pint+3, pint+8};
Line(lint+11) = {pint+4, pint+9};

llint = newll;
Line Loop(llint+0) = {lint+0, lint+9, -(lint+4), -(lint+8)};
Line Loop(llint+1) = {lint+1, lint+10 ,-(lint+5), -(lint+9)};
Line Loop(llint+2) = {lint+2, lint+11, -(lint+6), -(lint+10)};
Line Loop(llint+3) = {lint+3, lint+8, -(lint+7), -(lint+11)};

sint = news;
Ruled Surface(sint+0) = {llint+0};	MagnetMid[4] = sint+0;
Ruled Surface(sint+1) = {llint+1};	MagnetMid[5] = sint+1;
Ruled Surface(sint+2) = {llint+2};	MagnetMid[6] = sint+2;
Ruled Surface(sint+3) = {llint+3};	MagnetMid[7] = sint+3;

Rotate { {0, 1, 0}, {0, 0, 0}, Pi/2. } {
	Surface{sint+0, sint+1, sint+2, sint+3};
}
Rotate { {0, 0, 1}, {0, 0, 0}, Pi/2. } {
	Surface{sint+0, sint+1, sint+2, sint+3};
}

l = newl;
Line(l+0) = {pext+1, pint+1};
Line(l+1) = {pext+2, pint+2};
Line(l+2) = {pext+3, pint+3};
Line(l+3) = {pext+4, pint+4};
Line(l+4) = {pext+6, pint+6};
Line(l+5) = {pext+7, pint+7};
Line(l+6) = {pext+8, pint+8};
Line(l+7) = {pext+9, pint+9};

ll = newll;
Line Lopp(ll+0) = {lext+0, l+1, -(lint+0), -(l+0)};
Line Lopp(ll+1) = {lext+1, l+2, -(lint+1), -(l+1)};
Line Lopp(ll+2) = {lext+2, l+3, -(lint+2), -(l+2)};
Line Lopp(ll+3) = {lext+3, l+0, -(lint+3), -(l+3)};

Line Lopp(ll+4) = {lext+4, l+5, -(lint+4), -(l+4)};
Line Lopp(ll+5) = {lext+5, l+6, -(lint+5), -(l+5)};
Line Lopp(ll+6) = {lext+6, l+7, -(lint+6), -(l+6)};
Line Lopp(ll+7) = {lext+7, l+4, -(lint+7), -(l+7)};

s = news;
Ruled Surface(s+0) = {ll+0};	MagnetDown[0] = s+0;
Ruled Surface(s+1) = {ll+1};	MagnetDown[1] = s+1;
Ruled Surface(s+2) = {ll+2};	MagnetDown[2] = s+2;
Ruled Surface(s+3) = {ll+3};	MagnetDown[3] = s+3;

Ruled Surface(s+4) = {ll+4};	MagnetUp[0] = s+4;
Ruled Surface(s+5) = {ll+5};	MagnetUp[1] = s+5;
Ruled Surface(s+6) = {ll+6};	MagnetUp[2] = s+6;
Ruled Surface(s+7) = {ll+7};	MagnetUp[3] = s+7;

///////////////////////
///Physical surfaces///
///////////////////////
WALL = 1;
MAGNETUP = 2;
MAGNETDOWN = 3;
MAGNETMID = 4;

Physical Surface(WALL) = {Wall[]};
Physical Surface(MAGNETUP) = {MagnetUp[]};
Physical Surface(MAGNETDOWN) = {MagnetDown[]};
Physical Surface(MAGNETMID) = {MagnetMid[]};

////////////
///Volume///
////////////
sl = newsl;
Surface Loop(sl+0) = {Wall[]};
Surface Loop(sl+1) = {MagnetDown[], MagnetUp[], MagnetMid[]};

v = newv;
Volume(v+0) = {sl+1};
Volume(v+1) = {sl+0, (sl+1)};

//////////////////////
///Physical volumes///
//////////////////////
MAGNETVOLUME = 10;
VACUUMVOLUME = 20;

Physical Volume(MAGNETVOLUME) = {v+0};
Physical Volume(VACUUMVOLUME) = {v+1};

frederichecht · October 6, 2023, 2:37pm

just save the file as mesh file with:

/Applications/Gmsh.app/Contents/MacOS/gmsh benchmark1.msh -format mesh -save .msh -format mesh -save
and read the benchmark1.mesh
with
mesh3 Th(“benchmark1.mesh”);

Topic		Replies	Views
Gmsh: Error reading the mesh General Discussion	1	1028	April 3, 2020
Error reading Gmsh mesh in FreeFem++ General Discussion	3	314	January 11, 2022
Problem import 3d mesh from gmsh General Discussion	2	2334	February 19, 2020
Error reading 3D gmsh mesh General Discussion	0	644	May 30, 2020
Problem reading gmsh mesh in FreeFem General Discussion	0	475	May 18, 2021

Gmsh - errror element type

Related Topics