simple_dif.edp

//Définition du domaine:

border b(t=0,1){x=t;y=0;label=1;};
border d(t=0,1){x=1;y=t;label=1;};
border h(t=0,1){x=1-t;y=1-3*t/4;label=1;};
border g(t=0,1){x=0;y=(1-t)/4;label=1;};
//border cirh(t=0,2*pi){x=0.4+0.1*cos(t);y=0.4+0.1*sin(t);label=1;};

mesh Th = buildmesh(d(50)+h(50)+g(50/2)+b(50));

/*
border aa(t=1,0)  { x=0; y=t ;label=1;};        //  left beam
border b(t=0,2) { x=t; y=0 ;label=1;};        //  bottom of beam
border c(t=0,1)  { x=2; y=t ;label=1;};       //  rigth beam
border d(t=0,2) { x=2-t; y=1; label=1;};//  top beam 
border cirh(t=0,2*pi){x=0.8+0.1*cos(t);y=0.8+0.1*sin(t);label=1;};
mesh Th = buildmesh(aa(50)+b(50)+c(50)+d(50)+cirh(-50));
*/
//mesh Th =square(20,20);
//Th= adaptmesh(Th,1./30.,IsMetric=1,nbvx=10000);
plot(Th,wait=1);
fespace Vh(Th,P1);

//Variables :
Vh M,M0,v;
real mu,k,dt;
dt=0.1;
k=0.0003;
mu=0.500;
M0=0.1;

//Le problème de l'activateur:
problem Activator(M,v)=
int2d(Th)(M*v/dt)-
int2d(Th)(M0*v/dt)+//+?
int2d(Th)(mu*(dx(M)*dx(v)+dy(M)*dy(v)))+
int2d(Th)(k*M0*v)+on(1,M=0);

real nStep = 200;

for(int i=1;i<nStep-2;i++){

Activator;

M0=M;

plot(M,fill=1,value=1,wait=0);	
}