/usr/share/doc/freefem++/examples/examples++-tutorial/schwarz-gc.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 | // Schwarz without overlapping (Shur complenement Neumann -> Dirichet)
verbosity=2;
real cpu=clock();
int inside = 2;
int outside = 1;
border Gamma1(t=1,2){x=t;y=0;label=outside;};
border Gamma2(t=0,1){x=2;y=t;label=outside;};
border Gamma3(t=2,0){x=t ;y=1;label=outside;};
border GammaInside(t=1,0){x = 1-t; y = t;label=inside;};
border GammaArc(t=pi/2, 2*pi){ x= cos(t); y = sin(t);label=outside;};
int n=4;
mesh Th1 = buildmesh( Gamma1(5*n) + Gamma2(5*n) + GammaInside(5*n) + Gamma3(5*n));
mesh Th2 = buildmesh ( GammaInside(-5*n) + GammaArc(25*n) );
plot(Th1,Th2);
fespace Vh1(Th1,P1), Vh2(Th2,P1);
Vh1 u1=0,v1; Vh2 u2,v2;
Vh1 lambda=0; // take $\lambda \in V_{h1}$
/*
Remark, to day is not possible to
defined a function just on a border, so the $\ lambda $
function is defined on the all domain $\Omega_1$
by:
*/
int i=0; // for factorization optimization
problem Pb2(u2,v2,init=i,solver=Cholesky) =
int2d(Th2)( dx(u2)*dx(v2)+dy(u2)*dy(v2) )
+ int2d(Th2)( -v2)
+ int1d(Th2,inside)(-lambda*v2) + on(outside,u2= 0 ) ;
problem Pb1(u1,v1,init=i,solver=Cholesky) =
int2d(Th1)( dx(u1)*dx(v1)+dy(u1)*dy(v1) )
+ int2d(Th1)( -v1)
+ int1d(Th1,inside)(+lambda*v1) + on(outside,u1 = 0 ) ;
varf b(u2,v2,solver=CG) =int1d(Th1,inside)(u2*v2);
matrix B= b(Vh1,Vh1,solver=CG);
// $\lambda \longrightarrow \int_{\Gamma_i }(u_1-u_2) v_{1}$
func real[int] BoundaryProblem(real[int] &l)
{
lambda[]=l;
Pb1;
Pb2;
i++;
v1=-(u1-u2);
lambda[]=B*v1[];
return lambda[] ;
};
/*
Remark, the difference between the tow notation \ttCC{v1} and \ttCC{v1[]} is:
\ttCC{v1} is the finite element function and \ttCC{v1[]}
is the vector in the canonical basis of the finite element function \ttCC{v1} .
*/
Vh1 p=0,q=0;
// solve the problem with Conjugue Gradient
LinearCG(BoundaryProblem,p[],q[],eps=1.e-6,nbiter=100);
// compute the final solution, because CG works with increment
BoundaryProblem(p[]); // solve again to have right u1,u2
cout << " -- CPU time schwarz-gc:" << clock()-cpu << endl;
plot(u1,u2);
|