Preprints:

1. P. Kunštek, M. Vrdoljak: Stability of critical shapes in two-phase conductivity problems
2. M. Erceg, P. Kunštek, M. Vrdoljak: Hybrid Optimization Techniques for Multi-State Optimal Design Problems

Codes:

The following code was used to test numerical stability using the shape derivative gradient method. It is straightforward to verify that if the analysis in [1] demonstrates that the solutions are local extrema, the same result is valid for the previous method. For further details, please refer to [1].

// ball 2d 
load  "iovtk" 
ofstream file("AB.txt");
int fileNP = file.precision(7);
//---------------------------------------------------
// parameters
//---------------------------------------------------
real alpha = 1.0;
real beta = 2.0;


// max number of iterations
int kmax=400; 
// coefficinet for movemesh
real coef=4; // determining size of "jump"

real lambda = -0.62;

// radii 
real rI=0.76;
real rB=1.0;

//right-hand side 
func ff=6-6*sqrt(x^2+y^2);
func dxff=-6*x/sqrt(x^2+y^2);
func dyff=-6*y/sqrt(x^2+y^2);
func truncate=(((x^2+y^2)<0.9999*rB^2)? 1:0);
func hole=(((x^2+y^2)<0.1)? (0.0):(min((x^2+y^2)-0.1,1.0) ));

//---------------------------------------------------
// Meshes
//---------------------------------------------------
border OuterBoundary (t=0,2*pi) { x=rB*cos(t); y=rB*sin(t); label=2; }
border Interface (t=0,2*pi) { x=(rI+0.1*cos(12*t))*cos(t); y=(rI+0.1*cos(12*t))*sin(t); label=1;}

mesh Omega = buildmesh (OuterBoundary(60)+Interface(300));

Omega = adaptmesh (Omega, 1./25., IsMetric=1);
//---------------------------------------------------
// main loop
//---------------------------------------------------
for ( int k=0 ; k<kmax ; ++k)
{
	fespace Ph(Omega,P0);
	Ph reg=region;
	cout<<reg(0,0)<<" "<<reg(.99,0)<<endl;
	Ph aa = beta*(reg)+alpha*(1-reg);
	
	//------------------------------------------------
	fespace Vh(Omega,P2);
	Vh u,phi,theta1,theta2,psi1,psi2;
	problem PDEeq (u,phi)=
	int2d(Omega)(aa*(dx(u)*dx(phi)+dy(u)*dy(phi)))
	-int2d(Omega) (ff*phi)
	+ on (2, u=0);
	
	problem Direction (theta1,theta2,psi1,psi2) = 
	int2d(Omega) (aa*(dx(theta1)*dx(psi1)+dx(theta2)*dx(psi2)+dy(theta1)*dy(psi1)+dy(theta2)*dy(psi2)))
	-int2d(Omega) ((aa*(dx(u)^2*dx(psi1) - dx(psi1)*dy(u)^2 - dx(u)^2*dy(psi2) + dy(u)^2*dy(psi2) + 2*dx(u)*dx(psi2)*dy(u) + 2*dx(u)*dy(u)*dy(psi1))+2*dx(psi1)*ff*u + 2*dy(psi2)*ff*u + 2*dxff*psi1*u + 2*dyff*psi2*u))
	-int2d(Omega) (lambda*(1-reg)*(dx(psi1)+dy(psi2)))
	+ on (2, theta1=0, theta2=0 );
	PDEeq; // obtain u
	Direction; // obtain direction theta
	
	//------------------------------------------------
	// determining step
	//------------------------------------------------
	real minT00 =	checkmovemesh (Omega,[x,y])/5.0;
	while(1) 
	{
		real minTT = checkmovemesh (Omega,[x+coef*truncate*theta1,y+coef*truncate*theta2]);	
		if ( minT00 < minTT ) break ; // now movemesh can work
		coef/=10;
	}
	cout<<"determing step size ended\n stepsize:"<<coef<<endl;
	real volA=int2d(Omega)(1-reg);
	real volB=int2d(Omega)(reg);
	real L=int2d(Omega)(ff*u-lambda*reg);
	real er=sqrt(int1d(Omega,1)(theta1^2+theta2^2));
	if (er<1e-4) 
	coef=min(1.0,500*er);
	file<<k<<" "<<L<<" "<<er<<endl;
	savevtk("numerika/AB"+k+".vtk",Omega,aa,[theta1,theta2],dataname="kaj theta");
	
	// moving phaze alpha and beta
	Omega = movemesh (Omega,[x+coef*truncate*theta1,y+coef*truncate*theta2]);
	Omega = adaptmesh (Omega,1./25.,IsMetric=1);	
}