Trunc assertion fail ie==nbe

General Discussion
1

Hi guys,

once more I have a question regarding meshing in FreeFEM.

My goal is to build a hourglass-esque mesh with equidistant layering in z (this is later important to impose the BCs).
Layer-by-layer I want to add the layers consecutively to the already existing. My approach was to First build the full geometry and then trunc it for the calculations (to reduce comp. effort and keep the geometry).
The hourglass geoemtry is also build from single BUILDLAYERS (with a certain height) that originiate from a buildmesh which is defined by a BORDER with a certain radius (depending on the z level).

The TRUNC yields the problem :confused:

Here are the important code snippets:
real RR =5;
real HH =12.5;
real zmax=20;
real SLheight=0.4; //HH/SLheight = INT! UND zmax/SLheight=INT!
int nn= zmax/SLheight;
real rrnn;

for (int ii=0; ii<zmax/SLheight;ii=ii+1){ //zmax : max height, SLheight --> layer height
	real rrnn1=RR-ii*SLheight*RR/HH; //RR max radius of hourglass, HH hourglass height
	real rrnn2=ii*SLheight*RR/HH-(zmax-HH)*RR/HH;
	cout << rrnn1 << " "<< rrnn2 << endl;
	rrnn=max(rrnn1,rrnn2); //dont go to 0 radius, instead overlap in hourglass
	cout << rrnn << endl;
	border a(t=0,2*pi){x=(rrnn)*cos(t); y=(rrnn)*sin(t); label=-(ii+1);}
	nn=floor(2*pi*rrnn/SLheight); //keep elements roughly same size in xy &z direction
	cout << " n border " << nn << endl;
	mesh aa=buildmesh(a(nn));
	int[int] ll= labels(aa);		
	cout << "labels aa:" << ll << endl;	
	int[int] lup=[0,ii+1];
	int[int] ldown=[0,200+ii+1];
	mesh3 temp=buildlayers(aa,1,zbound=[ii*SLheight,(ii+1)*SLheight], labelup=lup,labeldown=ldown);
	int[int] kk= labels(temp);		
	cout << "lables temp: " <<kk << endl;	
	ACone=ACone+temp; //assemble cone step by step --> until hourglass is finished
}

nn=zmax/SLheight;	
ACone=Th;
for (int ii=0;ii<nn;ii=ii+1){
	cout << ii << endl;
	curL=ii+1;
	verbosity=99;
	Th=trunc(ACone,z<=((ii+1)*SLheight+0.01),flabel=67,removeduplicate=1,cleanmesh=1);
	verbosity=0;
	real Asurf=int2d(Th,curL) (1);
	cout << ii << " A [mm^2] " << Asurf << endl;
	plot(Th);
}

I hope someone of you knows what my problem is!

Best regards,
Jonas

2

Your mesh is nonconforming at the interface between two layers. Impossible to do in FreeFEM.

3

If you want to day this

the simple way is to build a piece of the hourglass using tetgen or
build a 2d of a vertical section

like this:

// example to build a mesh a cone
load “msh3”
load “medit”
// cone using buildlayers with a triangle
real RR=1,HH=1;
border Taxe(t=0,HH){x=t;y=0;label=0;};
border Hypo(t=1,0){x=HHt;y=RRt;label=1;};
border Vert(t=0,RR){x=HH;y=t;label=2;};

int nn=10;
real h= 1./nn;
mesh Th2= trunc(square(nn,nn,flags=2),x+y<1); //buildmesh( Taxe(HHnn)+ Hypo(sqrt(HHHH+RR*RR)nn) + Vert(RRnn) ) ;
plot(Th2,wait=1);
//medit(“circle”,ThT);

int MaxLayersT=(int(2piRR/h)/4)4;
func zminT = 0;
func zmaxT = 2pi;
func fx= ycos(z);// / max( abs(cos(z) ), abs(sin(z)));
func fy= ysin(z);// / max( abs(cos(z) ), abs(sin(z)));
func fz= x;
int[int] r1T=[0,0], r2T=[0,0,2,2];
int[int] r4T=[0,2];
mesh3 Th3T=buildlayers(Th2,coef= max(.01,y/max(x,0.4) ), MaxLayersT,zbound=[zminT,zmaxT],transfo=[fx,fy,fz],facemerge=1, region=r1T, labelmid=r2T);
medit(“cone”,Th3T,wait=1);
// FFCS: testing 3d plots
plot(Th3T,cmm=“cone”);

4

Thank you both for your feedback.

After reading “prj”'s answer in the morning I decided to “buildlayers” a cylinder and “trunc” the mesh in a manner that the hourglass geometry is obtaind. (by defining a function that yields 1 inside hourglass and 0 outside).

Now my code works! Thank you very much!

Best regards,

Jonas