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
|
import graph3;
import palette;
import animate;
settings.tex="pdflatex";
settings.render=0;
settings.prc=false;
unitsize(1cm);
animation a;
currentprojection=perspective(-20,-18,18);
currentlight=light(1,1,10);
int n=26;
real L=2.5;
real dx=2*L/n;
real CFL=0.125;
real dt=CFL*dx^2;
real[][] Z=new real[n][n];
real[][] W=new real[n][n];
guide initcond1=shift((-1,-1))*scale(0.5)*unitcircle;
guide initcond2=shift((0.5,0))*yscale(1.2)*unitcircle;
real f(pair p) {return (inside(initcond1,p)||inside(initcond2,p)) ? 2 : 0;}
//Initialize
for(int i=0; i < n; ++i)
for (int j=0; j < n; ++j)
Z[i][j]=f((-L,-L)+(2*L/n)*(i,j));
real f(pair t) {
int i=round((n/2)*(t.x/L+1));
int j=round((n/2)*(t.y/L+1));
if(i > n-1) i=n-1;
if(j > n-1) j=n-1;
return Z[i][j];
}
surface sf;
void advanceZ(int iter=20) {
for(int k=0; k < iter; ++k) {
for(int i=0; i < n; ++i)
for(int j=0; j < n; ++j) W[i][j]=0;
for(int i=1; i < n-1; ++i)
for(int j=1; j< n-1; ++j)
W[i][j]=Z[i][j]+(dt/dx^2)*(Z[i+1][j]+Z[i-1][j]+Z[i][j-1]+Z[i][j+1]
-4*Z[i][j]);
for(int i=0; i < n; ++i)
for(int j=0; j < n; ++j)
Z[i][j]=W[i][j];
};
}
pen[] Pal=Rainbow(96);
int endframe=40;
for(int fr=0; fr < endframe; ++fr) {
if(fr == 0) {// smoothing of initial state; no Spline, but full grid
advanceZ(3);
sf=surface(f,(-L,-L),(L,L),nx=n);
} else // use Spline and fewer grid points to save memory
sf=surface(f,(-L,-L),(L,L),nx=round(n/2),Spline);
sf.colors(palette(sf.map(zpart),Pal[0:round(48*max(sf).z)]));
draw(sf);
a.add();
erase();
advanceZ(30);
};
label(a.pdf(delay=400,"controls,loop"));
shipout(bbox(3mm,darkblue+3bp+miterjoin,FillDraw(fillpen=paleblue)),"pdf");
|