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, arXiv
3. M. Grbac, I. Ivec, M. Vrdoljak: Classical Optimal Designs for Stationary Diffusion with Multiple Phases, arXiv
4. A. Alphonse, P. Kunštek, M. Vrdoljak: Optimal design with uncertainties: a risk-averse approach, arXiv

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);	
}

The following code implements a homogenization procedure for a stochastic diffusion process. The stochastic forcing is represented through a truncated spectral expansion, and the associated deterministic auxiliary problems are solved on a fixed two-dimensional domain. The effective (homogenized) material coefficients are computed by solving a sequence of elliptic boundary value problems and by assembling the corresponding second-order moment tensors. An iterative optimization loop is then employed to determine the optimal homogenized conductivity under a prescribed volume constraint.

load  "iovtk"
load  "medit"
load  "UMFPACK64"

include  "updateCond2D.idp"


//---------------------------------------------------
// Output - name
//---------------------------------------------------

string name="primjerStochastic";
ofstream data("podaci/"+name+".txt");
ofstream geo("podaci/"+name+"geo.txt");
ofstream check("podaci/"+name+"check.txt");
int npdata = data.precision(12);
int npcheck = data.precision(12);


//---------------------------------------------------
// Parameters
//---------------------------------------------------

real alpha=1;  				// physical parameters
real beta=2;

int n=12;
int nn = n*n+1;
int D = nn+(nn-1)*nn/2;


//---------------------------------------------------
// Mesh - domain DD
//--------------------------------------------------- 
border boundary1(t=0,1) { x=t; y=0; label=0;}
border boundary2(t=0,1) { x=1; y=t; label=0;}
border boundary3(t=0,1) { x=1-t; y=1; label=0;}
border boundary4(t=0,1) { x=0; y=1-t; label=0;}

mesh DD = buildmesh(boundary1(10)+boundary2(10)+boundary3(10)+boundary4(10)); 
real Nmesh = 0.01;
DD = adaptmesh (DD, Nmesh, IsMetric=1); 
DD = trunc(DD,1,split=2);

int NbTrianglesDD = DD.nt;
real measureDD = DD.measure;
int NbVerticesDD = DD.nv;


//---------------------------------------------------
// Level set function - first initialization
//---------------------------------------------------
fespace VHF(DD,P2);
fespace VH(DD,P1);
fespace VH0(DD,P0);
VH0 h = hTriangle;
real hmax = h[].max;


//---------------------------------------------------
// Functions right-hand sides
//---------------------------------------------------
func truncate = (abs(x-0.5)<0.499 && abs(y-0.5)<0.499) ? 1 : 0;


//---------------------------------------------------
// Main state equation
//---------------------------------------------------
varf state(u,v) =
int2d(DD)(
  a11*dx(u)*dx(v)
+ a12*(dx(u)*dy(v)+dy(u)*dx(v))
+ a22*dy(u)*dy(v)
)
+ on(0,u=0.0);

varf RHS(u,v) =
int2d(DD)(f*v)
+ on(0,u=0.0);

// create RHS
real[int][int] bU(nn); 


//---------------------------------------------------
// Homogenization part 
// .................................
// Initial value for homogenization 
//---------------------------------------------------
// set all to alpha*I
theta=1; a11=alpha; a12=0; a22=alpha; 
thetaOld=theta;

//---------------------------------------------------
// Create u_i for homogenization
//---------------------------------------------------
matrix A = state(VH, VH, init=0);
for (int i=0; i < nn; ++i) {
    f = ff[i];
    real[int] b = RHS(0, VH);
    uu[i][] = A^-1 * b;
    // u = uu[i];
}

//---------------------------------------------------
// Create M for homogenization
//---------------------------------------------------
cout << "\n\nCreate M for homogenization\n\n" << endl;
cout << "\n\nSending to python\n\n" << endl;

real[int] dd(D);
{
    ofstream file("f.txt");
    file.precision(16);
    file << nn-1 << endl << D << endl;

    real ispis = -int2d(DD)(ff[0] * uu[0]);
    file << ispis << endl;

    for (int iter = 1; iter < nn; ++iter) {
        ispis = -sqrt(lambda[iter])
                * (int2d(DD)(ff[0] * uu[iter] + ff[iter] * uu[0]));
        file << ispis << endl;
    }

    for (int iter1 = 1; iter1 < nn; ++iter1)
        for (int iter2 = iter1; iter2 < nn; ++iter2) {
            ispis = -sqrt(lambda[iter1]) * sqrt(lambda[iter2])
                    * (int2d(DD)(ff[iter1] * uu[iter2]
                              + ff[iter2] * uu[iter1]));
            file << ispis << endl;
        }
}

exec("python3 skripta02.py");

{
    count = 0;
    ifstream file("Python_Out.txt");

    for (int iter = 1; iter < nn; ++iter) {
        file >> dd[count];
        count++;
    }

    for (int iter1 = 1; iter1 < nn; ++iter1)
        for (int iter2 = iter1; iter2 < nn; ++iter2) {
            file >> dd[count];
            count++;
        }
}


exec("rm f.txt");
exec("rm Python_Out.txt");

M11 = -(dx(uu[0]) * dx(uu[0]));
M12 = -((dx(uu[0]) * dy(uu[0]) + dy(uu[0]) * dx(uu[0])) / 2);
M22 = -(dy(uu[0]) * dy(uu[0]));

count = 0;
for (int iter = 1; iter < nn; ++iter) {
    M11 = M11 - 2 * sqrt(lambda[iter]) * dd[count]
          * (dx(uu[0]) * dx(uu[iter]));
    M12 = M12 - 2 * sqrt(lambda[iter]) * dd[count]
          * ((dx(uu[0]) * dy(uu[iter])
            + dy(uu[0]) * dx(uu[iter])) / 2);
    M22 = M22 - 2 * sqrt(lambda[iter]) * dd[count]
          * (dy(uu[0]) * dy(uu[iter]));
    count++;
}

for (int iter1 = 1; iter1 < nn; ++iter1)
    for (int iter2 = iter1; iter2 < nn; ++iter2) {
        M11 = M11 - 2 * sqrt(lambda[iter1]) * sqrt(lambda[iter2]) * dd[count]
              * (dx(uu[iter1]) * dx(uu[iter2]));
        M12 = M12 - 2 * sqrt(lambda[iter1]) * sqrt(lambda[iter2]) * dd[count]
              * ((dx(uu[iter1]) * dy(uu[iter2])
                + dy(uu[iter1]) * dx(uu[iter2])) / 2);
        M22 = M22 - 2 * sqrt(lambda[iter1]) * sqrt(lambda[iter2]) * dd[count]
              * (dy(uu[iter1]) * dy(uu[iter2]));
        count++;
    }