// heat.edp
//
//  Discussion:
//
//    Time dependent heat equation with inhomogeneous Dirichlet
//    and Neumann flux boundary conditions.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/heat/heat.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    14 June 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
//
//  Parameters:
//
//    real KAPPA, the thermal conductivity.
//
//    real THETA, the initial and boundary temperature.
//
real theta = 20.0;
real kappa = 0.1;
real dt = 0.5;
//
//  Boundary.
//
border Gamma1 ( t = 0.0, 2.0 * pi ) { x = cos ( t ); y = sin ( t ); }
border Gamma2 ( t = 0.0, 1.0      ) { x = -0.5 + t; y = - 0.3; }
//
//  Define the mesh.
//
mesh Th = buildmesh ( Gamma1 ( 100 ) + Gamma2 ( 60 ) );
plot ( Th, wait = true, ps = "heat_mesh.ps" );
//
//  Define the finite element space.
//
fespace Vh ( Th, P1 );
Vh uh;
Vh uold;
Vh vh;
//
//  The solution is taken to be constant at the initial time,
//  and at the "previous" timestep.
//
uh = theta;
uold = theta;
//
//  Define the weak form of the heat equation.
//
problem heatstep ( uh, vh, solver = LU ) = 
  int2d ( Th ) ( uh * vh / dt )
- int2d ( Th ) ( uold * vh / dt )
+ int2d ( Th ) ( kappa * dx ( uh ) * dx ( vh )
               + kappa * dy ( uh ) * dy ( vh ) )
+ int1d ( Th, Gamma1 ) ( kappa * vh )
+ on ( Gamma2, uh = theta );
//
//  Perform 6 timesteps.
//
for ( int it = 1; it <= 6; it++ )
{
  heatstep;
  plot ( uh, nbiso = 50, fill = true, value = true, wait = true, 
    ps = "heat_uh_" + it + ".ps" );
  uold = uh;
}
//
//  Terminate.
//
cout << "\n";
cout << "heat:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );