-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathStokes_2.edp~
49 lines (35 loc) · 921 Bytes
/
Stokes_2.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
38
39
40
41
42
43
44
45
46
47
48
// border b1(t=0,1){x=t;y=0;label=1;};
// border b2(t=0,1){x=1-t;y=t;label=2;};
// border b3(t=0,1){x=0 ;y=1-t;label=3;};
// mesh Th = buildmesh(b1(1) + b2(1) + b3(1));
// savemesh(Th, "trian_ref.msh");
// plot(Th);
border b1(t=0,1){x=3-t;y=1-t;label=1;};
border b2(t=0,1){x=2+t;y=0;label=2;};
border b3(t=0,1){x=3 ;y=t;label=3;};
mesh Th = buildmesh(b1(1) + b2(1) + b3(1));
savemesh(Th, "trian_ref.msh");
plot(Th);
//mesh Th = readmesh("carre.msh");
real nu=0.1;
real dt=0.1;
real alpha=1./dt;
real eps = 10e-6;
varf NSA(u,v) = int2d(Th) ( (dx(u)*dx(v) + dy(u)*dy(v)) );
// varf NSA(u,v) = int2d(Th) ( (u*v) );
fespace Uhvec(Th,P2);
matrix A=NSA(Uhvec,Uhvec);
real nbr = 0;
for(int i = 0; i < A.n; ++i){
for(int j = 0; j < A.m; ++j){
if( (abs(A(i,j)) < 10e6) ){
//nbr += A(i,j);
cout << A(i,j) << "\t";
}
}
cout << endl;
//cout << nbr << endl;
}
//{ofstream f("mat.dat");
//f << A;
//}