Plot segvaults after several mesh adaptations

General Discussion
1

I can’t upload the example since new users can’t upload apparently.
so source included at bottom.

Attached file appears to die in plot after several adaptations
with the following stack trace,
FreeFem++ FreeFem++ - version 4.12 (Thu Jan 5 20:06:10 EST 2023 - git no git) 64bits

Warning: too few vertices to split all internal edges
We lost 6 edges. Sorry!
– adaptmesh statistics: Nb triangles 17246 , h min 0.0496506 , h max 3.2041
area = 398 , M area = 28179.4 , M area/( |Khat| nt) 3.77348
infiny-regularity: min 0.251023 max 24.7795
anisomax 13.2528, beta max = 1.71978 min 0.061045
– mesh: Nb of Triangles = 17246, Nb of Vertices 9000
– Solve :
min 4.13289e-32 max 1
number of required edges : 0
Warning: too few vertices to split all internal edges
We lost 8 edges. Sorry!

Program received signal SIGSEGV, Segmentation fault.
bamg::Vertex::Smoothing (this=0x1c3b580, Th=…, BTh=…,
tstart=@0x1a4d2e0: 0x19182f8, omega=omega@entry=1)
at …/bamglib/Mesh2.cpp:1417
1417 tria->SetUnMarkUnSwap(1);
(gdb) bt
#0 bamg::Vertex::Smoothing (this=0x1c3b580, Th=…, BTh=…,
tstart=@0x1a4d2e0: 0x19182f8, omega=omega@entry=1)
at …/bamglib/Mesh2.cpp:1417
#1 0x0000000000bad8f9 in bamg::Triangles::SmoothingVertex (
this=this@entry=0x1579db0, nbiter=nbiter@entry=3, omega=omega@entry=1)
at …/bamglib/Mesh2.cpp:4460
#2 0x0000000000aba595 in Adaptation::operator() (this=,
stack=0x1539290) at lgmesh.cpp:974
#3 0x0000000000d3861e in E_F_F0F0<Fem2D::Mesh const**, Fem2D::Mesh const**, Fem2D::Mesh const*>::operator() (this=0x1551860, s=0x1539290)
at AFunction.hpp:1021
#4 0x000000000095b284 in ListOfInst::operator() (this=0x1548ca0, s=0x1539290)
at AFunction2.cpp:794
#5 0x00000000008e2ea6 in FFor (s=0x1539290, i0=, i1=0x1548a60,
i2=0x1548aa0, ins=0x1548ca0) at AFunction.cpp:284
#6 0x00000000008d8bc2 in E_F0_CFunc4::operator() (this=,
s=) at ./…/fflib/AFunction.hpp:908
#7 0x000000000095b284 in ListOfInst::operator() (this=0x15494f0, s=0x1539290)
at AFunction2.cpp:794
#8 0x000000000095b284 in ListOfInst::operator() (this=0x153a7d0, s=0x1539290)
at AFunction2.cpp:794
#9 0x00000000008d2f81 in CListOfInst::eval (this=0x7fffffffbeb8, s=0x1539290)
at ./…/fflib/AFunction.hpp:1486
#10 lgparse () at lg.ypp:360
#11 0x00000000008d7bda in Compile () at lg.ypp:853
#12 0x00000000008d8604 in mainff (argc=2, argv=0x7fffffffd388) at lg.ypp:1027
#13 0x00007ffff5568840 in __libc_start_main (
main=0x8c1480 <main(int, char**)>, argc=2, argv=0x7fffffffd388,
init=, fini=, rtld_fini=,
stack_end=0x7fffffffd378) at …/csu/libc-start.c:291
#14 0x00000000008c5209 in _start ()
(gdb)

load “medit”
load “msh3”
// mesh Th=Th2; // myqs(-1,0,1,1,-2,-1,5,4,1,0,1,1,0);
real x0=-1; real y0=0; real w=1; real h=1;
real ax0=-10; real ay0=-12; real aw=20; real ah=20 ;
real bx0=1; real by0=0; real bw=1; real bh=1;

int m=10;
int p=10;
int n=30;
// outside
border F01(t=0, 1) {x=ax0+(awt); y=ay0; }
border F11(t=0, 1) {x=ax0; y=ay0+ah-aht; }
border F10(t=0, 1) {x=ax0+aw; y=ay0+aht; }
border F00(t=0, 1) {x=ax0+aw-awt; y=ay0+ah; }
// first
border D01(t=0, 1) {x=x0+(wt); y=y0; }
border D11(t=0, 1) {x=x0; y=y0+h-ht; }
border D10(t=0, 1) {x=x0+w; y=y0+ht; }
border D00(t=0, 1) {x=x0+w-wt; y=y0+h; }

border E01(t=0, 1) {x=bx0+(bwt); y=by0; }
border E11(t=0, 1) {x=bx0; y=by0+bh-bht; }
border E10(t=0, 1) {x=bx0+bw; y=by0+bht; }
border E00(t=0, 1) {x=bx0+bw-bwt; y=by0+bh; }

mesh Th2=buildmesh(
F01(n)+F10(n)+F00(n)+F11(n)+
E01(-p)+E10(-p)+E00(-p)+E11(-p)+
D01(-m)+D10(-m)+D00(-m)+D11(-m));

mesh Th=Th2; // myqs(-1,0,1,1,-2,-1,5,4,1,0,1,1,0);
real tol=.5;
real cut=.01;

for(int i=0; i<10; ++i)
{
//fespace Xh(Th,RT03d);
//fespace Xh(Th,RT0);
fespace Xh(Th,P1);
Xh psi, wx;
macro grad(u) [dx(u),dy(u)]// def of grad operator
// Solve
solve potential(psi, wx)
= int2d(Th)(
grad(psi)'*grad(wx)) // scalar product
//+ on(E00,psi=0) + on(E10,psi=0) + on(E01,psi=0) + on(E11,psi=0)

  • on(D00,psi=0) + on(D10,psi=0) + on(D01,psi=0) + on(D11,psi=1)
  • on(F00,psi=1) + on(F10,psi=1) + on(F01,psi=1) + on(F11,psi=1)
    //)
    //+ on(C, psi = [uinfty1,uinfty2]'*[y,-x])
    //+ on(S, psi=-lift) // to get a correct value
    ;
    //plot(psi, wait=1);
    plot(Th, wait=0);

// this segfaults after a few iterations.
// changing hmin to 1/100 instead of 1/1000 seems to work better
Th = adaptmesh(Th, [dx(psi), dy(psi)], splitpbedge=1, abserror=0, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./10);
tol=tol/2;
cut=cut/2;
} // i

cout<<" done "<<endl;

//plot(Th2,wait=1);

2

Try changing the nbvx parameter. For example:

Th = adaptmesh(Th,u1,u2,...,nbvx = 1e5);
3

Thanks. I can play with it but generally a segfault is a bug .
If it is just a parameter value, it looks like the mesh is getting too fine
for the “smoother” whatever that does. For now I just disabled plotting
and dumped to a txt file and display in R.

4

The suggested change appears to work but the mesh is pretty dense
everywhere not just around the edges. Before it was getting coarse
in the space between squares and border. Since I’m running this
at home I’d like as few extra elements as possible.

5

Have you tried setting abserror=1?

6

Thanks, that seems to work with 1e4 and the grid looks more organizaed at 1e4 or 1e5.
Is the segfault just a pathological case or is it worth looking at? I was not going to
look at the code unless there is some point.

Th = adaptmesh(Th, [dx(psi), dy(psi)], splitpbedge=1, nbvx=1e5, abserror=1, cutoff=cut, err=tol, inquire=0, ratio=1.5, hmin=1./10);

7

Glad that helped. The segfault is probably just a pathological case.