/usr/share/doc/freefem++/examples/examples++-load/cmaes-VarIneq.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 | load "ff-cmaes"
int NN = 10;
mesh Th = square(NN,NN);
func f1=1.;
func f2=-1.;
func g1=0.;
func g2=0.1;
int iter=0;
int nadapt=1;
real starttol=1e-10,bctol=6.e-12;
fespace Vh(Th,P1);
Vh ou1,ou2;
real pena=100;
for(int al=0;al<nadapt;++al)
{
varf BVF(v,w) = int2d(Th)(0.5*dx(v)*dx(w) + 0.5*dy(v)*dy(w));
varf LVF1(v,w) = int2d(Th)(f1*w);
varf LVF2(v,w) = int2d(Th)(f2*w);
matrix A = BVF(Vh,Vh);
real[int] b1 = LVF1(0,Vh) , b2 = LVF2(0,Vh);
varf Vbord(v,w) = on(1,2,3,4,v=1);
//real[int] bord = Vbord(0,Vh);
//real[int] in = bord ? 0 : 1;
Vh In,Bord;
Bord[] = Vbord(0,Vh,tgv=1);
In[] = Bord[] ? 0:1;
Vh gh1=Bord*g1,gh2=Bord*g2;
//Function which create a vector of the search space type from
//two finite element functions
func int FEFToSSP(real[int] &fef1,real[int] &fef2,real[int] &ssp)
{
int kX=0;
for(int i=0;i<Vh.ndof;++i)
{
if(In[][i])
{
ssp[kX] = fef1[i];
ssp[kX+In[].sum] = fef2[i];
++kX;
}
}
return 1;
}
//Function spliting a vector from the search space and fills
//two finite element functions with it
func int SSPToFEF(real[int] &fef1,real[int] &fef2,real[int] &ssp)
{
int kX=0;
for(int i=0;i<Vh.ndof;++i)
{
if(In[][i])
{
fef1[i] = ssp[kX];
fef2[i] = ssp[kX+In[].sum];
++kX;
}
else
{
fef1[i] = gh1[][i];
fef2[i] = gh2[][i];
}
}
return 1;
}
func real IneqC(real[int] &X)
{
real[int] constraints(In[].sum);
for(int i=0;i<In[].sum;++i)
{
constraints[i] = X[i] - X[i+In[].sum];
constraints[i] = constraints[i] <= 0 ? 0. : constraints[i];
}
return constraints.l2;
}
func real J(real[int] &X)
{
Vh u1,u2;
SSPToFEF(u1[],u2[],X);
iter++;
real[int] Au1 = A*u1[], Au2 = A*u2[];
Au1 -= b1;
Au2 -= b2;
real val = u1[]'*Au1 + u2[]'*Au2;
val += pena * IneqC(X);
if(iter%200==199)
plot(u1,u2,nbiso=30,fill=1,dim=3,cmm="adapt level " + al + " - iteration " + iter + " - J = " + val,value=1);
return val ;
}
real[int] start(2*In[].sum);
if(al==0)
{
start(0:In[].sum-1) = 0.;
start(In[].sum:2*In[].sum-1) = 0.1;
}
else FEFToSSP(ou1[],ou2[],start);
real mini = cmaes(J,start,stopMaxFunEval=1000*(al+1),stopTolX=1.e-1/(10*(al+1)),initialStdDev=(0.025/(pow(100.,al))));
Vh best1,best2;
SSPToFEF(best1[],best2[],start);
Th = adaptmesh(Th,best1,best2);
ou1 = best1;
ou2 = best2;
}
|