//  molding.edp
//
//  Discussion:
//
//    Simulation of injection molding.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/molding/molding.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 May 2015
//
//  Author:
//
//    Florian De Vuyst
//
//  Reference:
//
//    Florian De Vuyst,
//    Numerical modeling of transport problems using freefem++ software -
//    with examples in biology, CFD, traffic flow and energy transfer,
//    HAL id: cel-00842234
//    https://cel.archives-ouvertes.fr/cel-00842234
//
cout << "\n";
cout << "molding:\n";
cout << "  FreeFem++ version\n";
cout << "  Fill a hollow mold with a liquid.\n";

real[int] PA(2); PA = [0, 0];
real[int] PB(2); PB = [9, 0];
real[int] PC(2); PC = [9, 9];
real[int] PD(2); PD = [8, 9];
real[int] PE(2); PE = [0, 9];
real [int ] PF(2); PF = [0, 6];
real [int ] PG(2); PG = [0,4];
real [int ] PH(2); PH = [1, 1];
real [int ] PK(2); PK = [1,8];
real [int ] PL(2); PL = [8,8];
real [int ] PM(2); PM = [8,5];
real [int ] PO(2); PO = [5,5];
real [int ] PP(2); PP = [5,4];
real [int ] PQ(2); PQ = [8,4];
real [int ] PR(2); PR = [8,1];

border c1(t=0,1) {x =(1-t)*PA[0]+t*PB[0]; y =(1-t)*PA[1]+t*PB[1];}
border c2(t=0,1) {x =(1-t)*PB[0]+t*PC[0]; y =(1-t)*PB[1]+t*PC[1];}
border c3(t=0,1) {x =(1-t)*PC[0]+t*PD[0]; y =(1-t)*PC[1]+t*PD[1];}
border c4(t=0,1) {x =(1-t)*PD[0]+t*PE[0]; y =(1-t)*PD[1]+t*PE[1];}
border c5(t=0,1) {x =(1-t)*PE[0]+t*PF[0]; y =(1-t)*PE[1]+t*PF[1];}
border c6(t=0,1) {x =(1-t)*PF[0]+t*PG[0]; y =(1-t)*PF[1]+t*PG[1];}
border c7(t=0,1) {x =(1-t)*PG[0]+t*PA[0]; y =(1-t)*PG[1]+t*PA[1];}
//
border c8(t=0,1){x=(1-t)*PH[0]+t*PK[0]; y=(1-t)*PH[1]+t*PK[1];}
border c9(t=0,1){x=(1-t)*PK[0]+t*PL[0]; y=(1-t)*PK[1]+t*PL[1];}
border c10(t=0,1){x=(1-t)*PL[0]+t*PM[0]; y=(1-t)*PL[1]+t*PM[1];}
border c11(t=0,1){x=(1-t)*PM[0]+t*PO[0]; y=(1-t)*PM[1]+t*PO[1];}
border c12(t=0,1){x=(1-t)*PO[0]+t*PP[0]; y=(1-t)*PO[1]+t*PP[1];}
border c13(t=0,1){x=(1-t)*PP[0]+t*PQ[0]; y=(1-t)*PP[1]+t*PQ[1];}
border c14(t=0,1){x=(1-t)*PQ[0]+t*PR[0]; y=(1-t)*PQ[1]+t*PR[1];}
border c15(t=0,1){x=(1-t)*PR[0]+t*PH[0]; y=(1-t)*PR[1]+t*PH[1];}
//
//  Define the mesh, Th.
//
int m = 2;

mesh Th = buildmesh ( c1(80*m)+c2(50*m)+c3(10*m)+c4(40*m)
  + c5(30*m)+c6(10*m)+c7(30*m)
  + c8(80*m)+c9(40*m)+c10(25*m)+c11(20*m)
  + c12(5*m)+c13(20*m)+c14(25*m)+c15(30*m) );
//
//  Plot the mesh.
//
plot ( Th, ps = "molding_mesh.ps" );
//
real gravity = 9.81;
real rhom = 1.0;
real rhop = 20.0;
real num = 0.1;
real nup = 10;
real dt = 1.0;
//
//  Define the finite element spaces.
//
fespace Uh(Th, P1b);
fespace Vh(Th, P1);
fespace Wh(Th, P1b);
Uh u1, u2, u1old, u2old, u1h, u2h;
Vh p, ph;
Wh u1view, u2view;
Wh rho, rhoold, rho2, mu, muold;
//
// Initialization
//
rho = rhom;
mu = rhom * num;
rho = rho + (rhop-rhom)*((x<0.2)*(y>PG[1])*(y<PF[1]));
mu = mu + (rhop*nup-rhom*num)*((x<0.2)*(y>PG[1])*(y<PF[1]));
u1 = 0.0;
u1old = u1;
u2 = 0.0;
u2old = u2;
//
plot ( rho, nbiso = 60, fill = true );
//
//  Define the weak form of the Navier-Stokes equations.
//
problem NavierStokes ( [u1, u2, p], [u1h, u2h, ph] ) =
    int2d (Th) (rho*u1*u1h/dt)
  - int2d (Th) (rhoold*convect([u1old,u2old], -dt, u1old)*u1h/dt)
  + int2d (Th) (mu*dx(u1)*dx(u1h)+mu*dy(u1)*dy(u1h))
  + int2d (Th) (dx(p)*u1h)
  + int2d (Th) (rho*u2*u2h/dt)
  - int2d (Th) (rhoold*convect([u1old,u2old], -dt, u2old)*u2h/dt)
  + int2d (Th) (mu*dx(u2)*dx(u2h)+mu*dy(u2)*dy(u2h))
  + int2d (Th) (dy(p)*u2h)
  + int2d (Th) (rho*gravity*u2h)
  + int2d (Th) ((dx(u1)+dy(u2))*ph+0.000001*p*ph)
  + on (c6, u1=0.3, u2=0)
  + on (c1,c2,c4,c5,c7,c8,c9,c10,c11,c12,c13,c14,c15, u1=0, u2=0);
//
real[int] viso(1);
viso = [5, 10, 15];

for ( int it = 0; it < 40; it++ ) 
{
  for ( int sit = 0; sit < 5; sit++ )
  {
//
//  Solve the Navier-Stokes equations.
//
    NavierStokes;
//
//  Convect the density and dynamic viscosity.
//
    u1old = u1;
    u2old = u2;
    rhoold = rho;
    muold = mu;
    rho = convect ( [ u1old, u2old ], -dt, rho );
    mu = convect ( [ u1old, u2old ], -dt, mu );
  }
// Th = adaptmesh(Th, rho, u1, u2); rho = rho; mu = mu; u1=u1; u2=u2;
  rho2 = rhop - rho;

  plot ( rho, fill = true, nbiso = 3, viso = viso, value = true, 
    ps = "molding_rho_it"+it+".ps");

// plot ( rho2, nbiso = 60, fill = true, grey = true, bw = tru,
//   ps = "molding_rho_it"+it+".ps" );
}
//
//  Terminate.
//
cout << "\n";
cout << "molding:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );