Trying to reconstruct the mesh results in an apparent infinite loop and eventual
use of all memory ( I had to wait minutes for system to come back to write this ).
Before it died, I caught a few stacks in the debugger and one suggested something
about bad faces. The code is below. Any thoughts? Thanks.
#0 0x00007ffff38ec4c7 in tetgenmesh::insertpoint(double*, tetgenmesh::triface*, tetgenmesh::face*, tetgenmesh::face*, tetgenmesh::insertvertexflags*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#1 0x00007ffff39069fb in tetgenmesh::split_tetrahedron(tetgenmesh::triface*, double*, int, int, tetgenmesh::insertvertexflags&) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#2 0x00007ffff39076b1 in tetgenmesh::repairbadtets(double, int) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#3 0x00007ffff3908931 in tetgenmesh::delaunayrefinement() ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#4 0x00007ffff3921718 in tetrahedralize(tetgenbehavior*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#5 0x00007ffff39236a5 in tetrahedralize(char*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#6 0x00007ffff388acba in ReconstructionRefine_tetgen (
switch_tetgen=switch_tetgen@entry=0x1541aa0 “raAQ”, Th3=…, nbhole=@0x7fffffffa9a0: 0,
tabhole=, nbregion=@0x7fffffffa9a4: 1, tabregion=,
nbfacecl=@0x7fffffffa9a8: 0, tabfacecl=0x15419d0, tsizevol=0x199bc10) at tetgen.cpp:1446
#7 0x00007ffff389ea53 in ReconstructionRefine_Op::operator() (this=0x1555ae0, stack=0x1538280)
at tetgen.cpp:2159
#0 0x00007ffff3963c5f in insphere(double*, double*, double*, double*, double*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#1 0x00007ffff38c3044 in tetgenmesh::insphere_s(double*, double*, double*, double*, double*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#2 0x00007ffff38fe8e4 in tetgenmesh::lawsonflip3d(tetgenmesh::flipconstraints*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#3 0x00007ffff3906b68 in tetgenmesh::split_tetrahedron(tetgenmesh::triface*, double*, int, int, tetgenmesh::insertvertexflags&) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#4 0x00007ffff39076b1 in tetgenmesh::repairbadtets(double, int) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#5 0x00007ffff3908931 in tetgenmesh::delaunayrefinement() ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#6 0x00007ffff3921718 in tetrahedralize(tetgenbehavior*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#7 0x00007ffff39236a5 in tetrahedralize(char*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#8 0x00007ffff388acba in ReconstructionRefine_tetgen (
switch_tetgen=switch_tetgen@entry=0x1541aa0 “raAQ”, Th3=…, nbhole=@0x7fffffffa9a0: 0,
tabhole=, nbregion=@0x7fffffffa9a4: 1, tabregion=,
nbfacecl=@0x7fffffffa9a8: 0, tabfacecl=0x15419d0, tsizevol=0x199bc10) at tetgen.cpp:1446
#9 0x00007ffff389ea53 in ReconstructionRefine_Op::operator() (this=0x1555ae0, stack=0x1538280)
at tetgen.cpp:2159
#0 0x00007ffff38c1c57 in tetgenmesh::memorypool::alloc() ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#1 0x00007ffff38c8170 in tetgenmesh::flippush(tetgenmesh::badface*&, tetgenmesh::triface*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#2 0x00007ffff38f08dc in tetgenmesh::insertpoint(double*, tetgenmesh::triface*, tetgenmesh::face*, tetgenmesh::face*, tetgenmesh::insertvertexflags*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#3 0x00007ffff39069fb in tetgenmesh::split_tetrahedron(tetgenmesh::triface*, double*, int, int, tetgenmesh::insertvertexflags&) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#4 0x00007ffff39076b1 in tetgenmesh::repairbadtets(double, int) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#5 0x00007ffff3908931 in tetgenmesh::delaunayrefinement() ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#6 0x00007ffff3921718 in tetrahedralize(tetgenbehavior*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) () from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#7 0x00007ffff39236a5 in tetrahedralize(char*, tetgenio*, tetgenio*, tetgenio*, tetgenio*) ()
from /home/ubuntu/dev/freefem/install/lib/ff++/4.12/lib/tetgen.so
#8 0x00007ffff388acba in ReconstructionRefine_tetgen (
///home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$ vi Poisson-cube-ballon.edp
//marchywka@happy:/home/ubuntu/dev/freefem/FreeFem-sources-master/examples/3d$
load “msh3”
load “tetgen”
load “medit”
load “mshmet”
//
real volumetet; // use in tetg.
meshS Thx0,Thx1, Thy0, Thy1, Thz0, Thz1;
int nx=10,ny=10,nz=10;
func meshS mycube ( real x0, real x1, real y0,real y1, real z0, real z1, int orient, int lb)
{
meshS ThHex;
// first build the 6 faces of the hex.
//real x0=-1,x1=1;
//real y0=-1.1,y1=1.1;
//real z0=-1.2,z1=1.2;
// a volume of on tet.
volumetet= (x1-x0)(y1-y0)(z1-z0)/ (nxnynz) /6.;
mesh Thx = square(ny,nz,[y0+(y1-y0)*x,z0+(z1-z0)*y]);
mesh Thy = square(nx,nz,[x0+(x1-x0)*x,z0+(z1-z0)*y]);
mesh Thz = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
int[int] refz=[0,5]; // bas
int[int] refZ=[0,6]; // haut
int[int] refy=[0,3]; // devant
int[int] refY=[0,4]; // derriere
int[int] refx=[0,1]; // gauche
int[int] refX=[0,2]; // droite
int f1=-1;
int f2=1;
if (orient==1) { f1=1; f2=-1; }
int[int] l0=[0,lb+0];
//meshS
Thx0 = movemesh23(Thx,transfo=[x0,x,y],label=l0,orientation=f1,region=refx,removeduplicate=false);
++l0[1];
//meshS
Thx1 = movemesh23(Thx,transfo=[x1,x,y],label=l0,orientation=f2,region=refX,removeduplicate=false);
++l0[1];
//meshS
Thy0 = movemesh23(Thy,transfo=[x,y0,y],label=l0,orientation=f2,region=refy,removeduplicate=false);
++l0[1];
//meshS
Thy1 = movemesh23(Thy,transfo=[x,y1,y],label=l0,orientation=f1,region=refY,removeduplicate=false);
++l0[1];
//meshS
Thz0 = movemesh23(Thz,transfo=[x,y,z0],label=l0,orientation=f1,region=refz,removeduplicate=false);
++l0[1];
//meshS
Thz1 = movemesh23(Thz,transfo=[x,y,z1],label=l0,orientation=f2,region=refZ,removeduplicate=false);
cout<<“l0=”<<l0<<endl;
//medit(" — ", Thx0,Thx1,Thy0,Thy1,Thz0,Thz1);
ThHex = Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;
return ThHex;
}
real x0=-1,x1=1; real yy0=-1.1,yy1=.1; real z0=-1.2,z1=.2;
real xa0=-.5; real ya0=-.8; real za0=-1;
real xb0=-.1; real yb0=-.8; real zb0=-1;
real sz=.25;
meshS ThHex = mycube(x0,x1,yy0,yy1,z0,z1,1,0);
real vol=volumetet;
meshS Thsrc=Thx0;
nx=10; ny=10; nz=10;
meshS Tha = mycube(xa0,xa0+sz,ya0,ya0+sz,za0,za0+sz,-1,10);
meshS Thb = mycube(xb0,xb0+sz,yb0,yb0+sz,zb0,zb0+sz,-1,20);
meshS Thc= Tha+Thb+ThHex;
//medit(" asdf ",Thc);
real[int] domaine = [-.9,-1,-1.1,.5,vol];
mesh3 Th = tetg(Thc,switch=“pqaAAYYQ”,nbofregions=1,regionlist=domaine);
// Tetrahelize the interior of the cube with tetgen
//medit(“tetg”,Th,wait=1);
for( int i=0; i<1; ++i)
{
cout<<" ++++++++++++++++++++++++++++++++++++ "<<i<<endl;
fespace Xh(Th,P23d);
Xh psi, wx,dxpsi,dypsi;
macro grad(u) [dx(u),dy(u),dz(u)]// def of grad operator
// Solve
solve potential(psi, wx)
= int3d(Th)(
grad(psi)'grad(wx)) // scalar product
//+int2d(Thsrc) (2wx)
+int2d(Th,1) (.200wx)
+int2d(Th,2) (-.200wx)
//+on(Thl,psi=1)
//+on(102,psi=1)
+on(1,psi=0)
//+on(1,psi=0)
//+on(1,psi=1)
//+on(3,dx(psi)=0)
//+on(4,psi=1)
//+on(5,psi=1)
//+on(6,psi=1)
;
dxpsi=dx(psi);
dypsi=dy(psi);
{
ofstream ff(“psi.txt”);
for (int ix = 0; ix < Th.nv; ix++)
{
// this does produce dups, maybe confusing everything…
for (int j = 0; j < 3; j++) // foo[ix]+=Thx[ix][j]*Thx[ix][j];
//ff << Th[ix][j].x << " “<< Th[ix][j].y << " " << psi[Xh(ix,j)] << endl;
ff << Th[ix][j].x << " “<< Th[ix][j].y << " "
<<Th[ix][j].z<<” "
<< dxpsi[Xh(ix,j)]<<” “<<dypsi[Xh(ix,j)]<<” "<<psi[Xh(ix,j)] << endl;
//ff << Th[ix][0].x << " " << Th[ix][0].y << " " << psi[Xh(ix,0)] << “\n\n\n”;
}
} // scoping
//real[int] viso=[1.2,1.15,1.1,1.02,1.01];
//real[int] viso=[2,1.1,1.05,1.01,1,.999,.99,.9];
//real[int] viso=[1.7,1.1,1,.999,.99,.9];
real[int] viso=[0,.01,.05,.1,.2,.3];
plot(psi, viso=viso, value=true, wait=1);
//plot(psi, value=true, wait=1);
//real[int] xxx=[dx(psi),dy(psi)];
//plot(xxx, value=true, wait=1);
Xh hh; real lerr=.001;
//hh=mshmet(Th,(dx(psi)*dx(psi)+dy(psi)*dy(psi)+dz(psi)*dz(psi)),normalization=1,aniso=0,nbregul=1,hmin=Th.measure/100,hmax=Th.measure/10,err=lerr);
hh=mshmet(Th,psi,normalization=1,aniso=0,nbregul=1,hmin=Th.measure/100,hmax=Th.measure/10,err=lerr);
cout<<" +++++++++ RECO ++++++++++++++++++++ "<<i<<endl;
// zero regions?
Th=tetgreconstruction(Th,switch=“raAQ”,nbofregions=1,regionlist=domaine,sizeofvolume=hhhhhh/6.);
medit(“tetg”,Th,wait=1);
} //i