/usr/share/doc/freefem++/examples/examples++-load/convect_dervieux.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | bool fast=true;
load "mat_dervieux"; // external module in C++ must be loaded
border a(t=0, 2*pi){ x = cos(t); y = sin(t); }
mesh th = buildmesh(a(100));
fespace Vh(th,P1);
Vh vh,vold,u1 = y, u2 = -x;
Vh v = exp(-10*((x-0.3)^2 +(y-0.3)^2)), vWall=0, rhs =0;
real dt = 0.025;
// qf1pTlump means mass lumping is used
problem FVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
- int2d(th,qft=qf1pTlump)(vold*vh/dt)
+ int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh)
+ rhs[] ;
matrix A;
MatUpWind1(A,th,vold,[u1,u2]);
if(fast)
{
varf vFVM(v,vh) = int2d(th,qft=qf1pTlump)(v*vh/dt)
- int1d(th,a)(((u1*N.x+u2*N.y)<0)*(u1*N.x+u2*N.y)*vWall*vh) ;
real[int] rhs0=vFVM(0,Vh);
matrix M=vFVM(Vh,Vh,solver=CG);
A=-A+M;
for ( real t=0; t< pi ; t+=dt)
{
vold=v;
rhs[]=rhs0;
rhs[] += A * vold[] ;
v[]= M^-1*rhs[];
plot(v,wait=0);
}
}
else
for ( real t=0; t< pi ; t+=dt){
vold=v;
rhs[] = A * vold[] ;
FVM;
plot(v,wait=0);
};
|