1load "Curvature" 2load "isoline" 3cout << " tersca " << Tresca(1,2,3,0,0,0) << endl;; 4real R= 1; 5real meshsize= 0.1; 6border DC(t=pi/2,-pi/2) { x=R*cos(t); y=R*sin(t); label=1;} 7border axe(t=R,-R) {x=0; y=t; label=2;} 8 9mesh Th=buildmesh(DC(-R*pi/meshsize)+axe(2*R/meshsize)); 10 11fespace Vh(Th,P1); 12Vh ca,c,cc,ccc;; 13c[]=curvature(Th,1); 14real cmean= int1d(Th,1,qforder=1)(c)/int1d(Th,1,qforder=1)(1.); 15cout<< " cmean = " << cmean << " == " << 1/R << endl; 16plot(c, wait=1,cmm=" curvature" ); 17verbosity=1000; 18ca[]=raxicurvature(Th,1); 19verbosity=1; 20plot(ca, wait=1,cmm="axi curvature" ); 21real s=int1d(Th,1,qforder=1)(x*2*pi); 22cout << " s = " << s << " " << 4*pi*R*R << endl; 23real cmeana= int1d(Th,1,qforder=1)(2*pi*ca)/s; 24solve AAA(cc,ccc)= int2d(Th)(cc*ccc*1e-10)+ int1d(Th,1)(x*cc*ccc)- int1d(Th,1)(ca*ccc); 25//cc = ca/max(x,0.01); 26plot(cc,wait=1); 27cout<< " cmeana = " << cmeana << " == " << 2/R << endl; 28assert(abs(cmean-1/R)< 0.05/R); 29 30int[int] ll=[1,2]; 31real[int,int] b12(1,3); 32real l12=extractborder(Th,ll,b12); 33cout << " size b12 = " << b12.n << " x " << b12.m << endl; 34border BB(t=0,1){ P=Curve(b12,t);label=3;} 35plot(BB(100),wait=1);