Electrostatic stripline analysis segfaults during adaptation

General Discussion
1

Sorry for the large blob of text, but as a new user I cannot upload attachments.
The code which has been included calculates the impedance of stripline electrostatically. It is a conversion of this , VERY old code, which I found on the signal integrity mailing list.

It was relatively straightforward to bring it up to modern freefem. As you can see I have implemented adaptation based on the Poisson example.

There is one problem, for NAdapt = 12 , Freefem segfaults.

I suspect that the ever increasing mesh density is simply causing an issue of some sort, perhaps a too small triangle, which in turn causes the segfault.
However I was hoping someone could make sure that I’m not making some subtle error in the code.
If the code is correct, is there some way to gauge if the adaptation process should stop before it segfaults ? :grinning:

Thank you !

load "iovtk"

// permittivity constant of free space
real eps0 = 8.854e-12;

// permeability constant of free space
real mu = 4*pi*1E-7;

// total effective width of the stripline
real width=300;

// metal thickness
real thick=2.5;

// width of the center conductor
real traceWidth=70;

// total height of the stripline
real H=90;

// height of the strip above the bottom ground
// if this is not = H/2 then the stripline is asymmetric
real hd=45;

// relative dielectric constant of the substrate
real eps = 2.33;

// speed of light in the material
real vc = 1.0 / sqrt(eps0 * eps * mu);

// the direction of the edges matters or the interior will be meshed

/*
  Triangulation keywords assume that the domain is defined as being on
  the left (resp right) of its oriented parameterized boundary
*/

// meshing constant
// this is the number of mesh points along the edge

// on the side of the metal trace this number must be very small or
// the edge of the trace will mesh too much !


int m1 = 5;
int m2 = 40;

// outer boundary
// bottom
border B01(t=0, width) { x = t; y = 0.0; }
// right
border B02(t=0, H) {x = width; y = t; }
// top
border B03(t=0, width) { x=width-t; y=H; }
// left
border B04(t=0, H) {x=0.0; y=H-t;}

// bottom edge
border C01(t=0, traceWidth) { x=t + 0.5*width - 0.5*traceWidth; y = hd; }
// right edge
border C02(t=0, thick) { x=0.5*width + 0.5*traceWidth; y=t + hd; }
// top edge
border C03(t=0, traceWidth) { x = -t + traceWidth + 0.5*width - 0.5*traceWidth; y = hd+thick; }
// left  edge
border C04(t=0, thick) { x = 0.5*width - 0.5*traceWidth; y = -t + hd + thick; }

mesh Th = buildmesh(C01(-m2) + C02(-m1) + C03(-m2) + C04(-m1) + B01(m2) + B02(m2) + B03(m2) + B04(m2));

// plot(Th);

fespace Vh(Th, P1);

Vh u, v;

macro grad(u) [ dx(u), dy(u) ] //

func f = 0.0;
func g = 1.0;

// program segfaults for NAdapt = 12
// strangely the Poisson example does NOT.

int NAdapt = 10;

/*

Remember that in Poisson's equation , the equation solves for the
Voltage in the solution space, NOT the electric field.

del^2 V = rho/epsilon

where

E = -del V
*/

problem Poisson(u, v, solver=CG, eps=-1.e-6)
  = int2d(Th)(grad(u)' * grad(v))
  + on(C01, C02, C03, C04, u=g)
  + on(B01, B02, B03, B04, u=f)
  ;

real error = 0.1;
real coef = 0.1^(1./5.);

for (int i = 0; i < NAdapt; i++){
	// Solve
	Poisson;
	
	// Plot
	plot(Th, u);
        // plot(u);
	
	// Adaptmesh
	Th = adaptmesh(Th, u, inquire=1, err=error);
	error = error * coef;
	
	// TO BUILD IMG ONLY
	string DataName = "Poisson";
	int[int] Order = [1];
	// savevtk("Result_"+i+".vtu", Th, u, dataname=DataName, order=Order);
	// END
 }


real energy = eps0*eps* int2d(Th) ((dx(u))*(dx(u)) + (dy(u))*(dy(u)));
cout << "energy = " << energy << endl;

real impenerg = 1/(vc * energy);
cout << "Impedance (in ohms) calculated from energy:" << impenerg << endl;