/usr/share/doc/freefem++/examples/examples++-tutorial/AdaptResidualErrorIndicator.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 | // macro the get the current mesh size
// parameter
// in: Th the mesh
// Vh P1 fespace on Th
// out :
// h: the Vh finite element finite set to the current mesh size
macro MeshSizecomputation(Th,Vh,h)
{ /* Th mesh
Vh P1 finite element space
h the P1 mesh size value */
real[int] count(Th.nv);
/* mesh size (lenEdge = $\int_e 1$) */
varf vmeshsizen(u,v)=intalledges(Th,qfnbpE=1)(v);
/* number of edge / par vertex */
varf vedgecount(u,v)=intalledges(Th,qfnbpE=1)(v/lenEdge);
/*
computation of the mesh size
----------------------------- */
count=vedgecount(0,Vh);
h[]=0.;
h[]=vmeshsizen(0,Vh);
cout << " count min = "<< count.min << " " << count.max << endl;
h[]=h[]./count;
cout << " -- bound meshsize = " <<h[].min << " " << h[].max << endl;
} // end of macro MeshSizecomputation
// macro to remesh according the de residual indicator
// in:
// Th the mesh
// Ph P0 fespace on Th
// Vh P1 fespace on Th
// vindicator the varf of to evaluate the indicator to ^2
// coef on etameam ..
// ------
macro ReMeshIndicator(Th,Ph,Vh,vindicator,coef)
{
Vh h=0;
/*evalutate the mesh size */
MeshSizecomputation(Th,Vh,h);
Ph etak;
etak[]=vindicator(0,Ph);
cout << " global Eta : " << sqrt(etak[].sum) << " ......... " << Th.nv<< endl;
etak[]=sqrt(etak[]);
plot(etak,ps="arei-etak.eps",fill=1,value=1);
real etastar= coef*(etak[].sum/etak[].n);
cout << " etastar = " << etastar << " sum=" << etak[].sum << " " << endl;
/* here etaK is discontinous
we use the P1 L2 projection with mass lumping . */
Vh fn,sigma;
varf veta(unused,v)=int2d(Th)(etak*v);
varf vun(unused,v)=int2d(Th)(1*v);
fn[] = veta(0,Vh);
sigma[]= vun(0,Vh);
fn[]= fn[]./ sigma[];
fn = max(min(fn/etastar,3.),0.3333) ;
/* new mesh size */
h = h / fn ;
/* plot(h,wait=1); */
Th=adaptmesh(Th,IsMetric=1,h,splitpbedge=1,nbvx=10000);
}
// the classical L space problem.
// mesh definition
border ba(t=0,1.0){x=t; y=0; label=1;};
border bb(t=0,0.5){x=1; y=t; label=2;};
border bc(t=0,0.5){x=1-t; y=0.5;label=3;};
border bd(t=0.5,1){x=0.5; y=t; label=4;};
border be(t=0.5,1){x=1-t; y=1; label=5;};
border bf(t=0.0,1){x=0; y=1-t;label=6;};
mesh Th = buildmesh (ba(6) + bb(4) + bc(4) +bd(4) + be(4) + bf(6));
// FE space definition ---
fespace Vh(Th,P1); // for the mesh size
fespace Ph(Th,P0); // for the indicator
real hinit=0.2; //
Vh h=hinit; // the FE fonction for the mesh size
// to build a mesh with a given mesh size : meshsize
Th=adaptmesh(Th,h,IsMetric=1,splitpbedge=1,nbvx=10000);
plot(Th,wait=1,ps="RRI-Th-init.eps");
Vh u,v;
func f=(x-y);
problem Poisson(u,v) =
int2d(Th,qforder=5)( u*v*1.0e-10+ dx(u)*dx(v) + dy(u)*dy(v))
- int2d(Th,qforder=5)( f*v);
varf indicator2(unused,chiK) =
intalledges(Th)(chiK*lenEdge*square(jump(N.x*dx(u)+N.y*dy(u))))
+int2d(Th)(chiK*square(hTriangle*(f+dxx(u)+dyy(u))) );
for (int i=0;i< 10;i++)
{
u=u;
Poisson;
plot(Th,u,wait=1);
real cc=0.7;
if(i>5) cc=1;
if(i<9)
ReMeshIndicator(Th,Ph,Vh,indicator2,cc);
plot(u,Th,wait=1,ps="arei-Thu.eps",value=1);
}
|