-
Notifications
You must be signed in to change notification settings - Fork 197
/
Copy pathexample20.edp
35 lines (26 loc) · 1.07 KB
/
example20.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
// example20.edp
// original file: membrane.edp
real theta = 4.*pi/3.;
real a = 2.,b = 1.; // semimajor and semiminor axes
func z=x; // elevation of Gamma1 boundary
border Gamma1(t=0, theta) { x = a * cos(t); y = b*sin(t); }
border Gamma2(t=theta, 2*pi) { x = a * cos(t); y = b*sin(t); }
mesh Th = buildmesh( Gamma1(100) + Gamma2(50) ); // construct mesh
fespace Vh(Th,P2); // P2 conforming triangular FEM
Vh phi,w, f=1; // phi is shape function, w is test function, f is load
// problem definition
solve Laplace(phi, w) = int2d(Th)( dx(phi)*dx(w) + dy(phi)*dy(w) )
- int2d(Th)( f*w ) + on(Gamma1, phi=z);
plot(phi, wait=true, fill=true, ps="solution20.eps"); //Plot solution
plot(Th, wait=true, ps="mesh20.eps"); //Plot mesh
// to build a gnuplot data file
{ ofstream ff("graph20.txt");
for (int i=0;i<Th.nt;i++){
for (int j=0; j <3; j++){
ff<<Th[i][j].x << " "<< Th[i][j].y<< " "<<phi[][Vh(i,j)]<<endl;
}
ff<<Th[i][0].x << " "<< Th[i][0].y<< " "<<phi[][Vh(i,0)]<<endl
<<endl<<endl;
}
}
// savemesh(Th,"Th.msh");