/
tutomesh1d.edp
37 lines (33 loc) · 1.06 KB
/
tutomesh1d.edp
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
load "msh3" // bizarre
int nn= 4;
int[int] ll=[10,11];
meshL Th= segment(nn,[x*nn,0,0],label=ll,region=1); //
// print information of Th
// bug label corrige ??
cout << " Element : "<<endl;
for(int k=0; k< Th.nt; ++k)
{
cout << Th[k][0] << " "<< Th[k][1] << ": "
<< Th[k][0].x << " " << Th[k][0].label << " , "
<< Th[k][1].x << " " << Th[k][1].label << " " << Th[k].region << endl;
}
for(int i=0; i< Th.nv; ++i)
{
cout << Th(i).x << " " << Th(i).y << " " << Th(i).y << " " << Th(i).label << endl;
}
for(int i=0; i< Th.nbe; ++i)
{
cout << " be " << i << " " << Th.be(i) << " " << Th.be(i).label <<endl;
}
// verif integration ..
func real track() { cout << " Track " << x << " N " << N << " Tl= " << Tl << endl; return 1.;}
fespace Vh(Th,P1);
varf va(u,v) = int1d(Th)(track()*dx(u)*dx(v)) + int0d(Th,10)(track()*u*v) - int0d(Th,10)(track()*v);
matrix A = va(Vh,Vh);
verbosity= 10;
real[int] b= va(0,Vh);
cout << " A = " << A << endl;
cout << " b = " << b << endl;
assert(abs(b.sum+1)< 1e-7);
assert(abs(A(0,0)-2)<1e-7);
assert(abs(A(nn,nn)-1)<1e-7);