// thermal_design.edp
//
//  Discussion:
//
//    A thermal design problem.  The wall is a composite material
//    made of one cheap and one expensive material.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/thermal_design/thermal_design.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    17 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 << "thermal_design\n";
cout << "  FreeFem++ version\n";

real Lx = 4;
real Ly = 5;
real lr = 1;
real lx = 0.3;
real Tr = 20;
real Text = -5.0;
real hy = 0.1;
real hx = 0.1;
real mu = 0.2;
real kappaa = 10.0;
real kappac = 1.0;
real kappae = 0.008;

real V = lx * Ly;
real Ve = mu * V;
real xem = Lx + hx;
real deltaxe = Ve / ( Ly - 2.0 * hy );
//
//  Define boundaries.
//
border c1 ( t = 0, Lx )               { x = t; y = 0.0; }
border c1w ( t = Lx, Lx + lx )        { x = t; y = 0.0; }
border cext ( t = 0, Ly )             { x = Lx + lx; y = t; }
border c3w ( t = Lx+lx, Lx )          { x = t; y = Ly; }
border c3 ( t = Lx, 0.0 )             { x = t; y = Ly; }
border c4 ( t = Ly, Ly/2+lr/2 )       { x = 0; y = t; }
border cr ( t = Ly/2+lr/2,Ly/2-lr/2 ) { x = 0; y = t; }
border c5 ( t = Ly/2-lr/2,0.0 )       { x = 0; y = t; }
border c0 ( t = 0, Ly )               { x = Lx; y = t; }

border ce1 ( t = xem, xem + deltaxe ) { x = t; y = hy; }
border ce2 ( t = hy, Ly - hy )        { x = xem+deltaxe; y = t; }
border ce3 ( t = xem+deltaxe, xem )   { x = t; y = Ly-hy; }
border ce4 ( t = Ly-hy, hy )          { x = xem; y = t; }
//
//  Generate the mesh.
//
mesh Th = buildmesh (
  c1(30) + c1w(10) + cext(100) 
+ c3w(10) + c3(30) + c4(15) + cr(15) + c5(15) + c0(100) 
+ ce1(15) + ce2(500) + ce3(15) + ce4(500) );

plot ( Th, wait = true, ps = "thermal_design_mesh.ps" );
//
//  Generate the second mesh.
//
mesh Th2 = buildmesh (
  c1(30) + c0(100) + c3(30) + c4(15) + cr(15) + c5(15) );

plot ( Th, wait = true, ps = "thermal_design_mesh2.ps" );
//
//  Pick points that can identify the two regions.
//
real region1 = Th ( Lx/2, Ly/2).region;
real region2 = Th ( xem+deltaxe/2,Ly/2).region;
//
//  Define the finite element space.
//
fespace Vh ( Th, P1 );
Vh kappa;
Vh th;
Vh vh;
//
//  Define a piecewise linear function KAPPA which has the values
//
//  KAPPAA in region1,
//  KAPPAE in region2,
//  KAPPAC elsewhere.
//
kappa = kappac + ( kappaa - kappac ) * ( region == region1 )
               + ( kappae - kappac ) * ( region == region2 );

plot ( kappa, nbiso = 60, fill = true, value = true, 
  ps = "thermal_design_kappa.ps" );
//
//  Define the weak form of the thermal design problem.
//
problem thermaldesign ( th, vh ) =
  int2d ( Th ) ( kappa * dx ( th ) * dx ( vh )
               + kappa * dy ( th ) * dy ( vh ) )
+ on ( cext, th = Text )
+ on ( cr, th = Tr );
//
//  Solve the problem.
//
thermaldesign;
//
//  Define colors to be associated with contour plot levels.
//
//  Since the name "colorhsv" occurs nowhere else, is this a
//  recognized keyword with a default value?
//
real[int] colorhsv = [
  4.0 / 6.0, 1.0, 0.5,
  4.0 / 6.0, 1.0, 1.0,
  5.0 / 6.0, 1.0, 1.0,
  1.0,       1.0, 1.0,
  1.0,       0.5, 1.0 ];
//
//  Define the levels of interest for a contour plot.
//
real[int] viso(26);
for ( int i = 0; i < viso.n; i++ )
{
  viso[i] = i - 5;
}
//
//  Plot T.
//
plot ( th, viso = viso(0:viso.n-1), value = true, fill = true, 
  ps = "thermal_design_tfield.ps" );
//
//  Compute performance indicators defined on the second mesh:
//  J1 = mean value of T;
//  J2 = minimum value of T;
//  J3 = standard deviation of T.
//
fespace Vh2 ( Th2, P1 );
Vh2 th2 = th;
real J1;
real J2;
real J3;

J1 = int2d ( Th2 ) ( th2 ) / Th2.area;
J2 = th2[].min;
J3 = sqrt ( int2d ( Th2 ) ( ( th2 - J1 )^2 ) / Th2.area );

cout << "  Mean temperature = " << J1 << "\n";
cout << "  Minimum temperature = " << J2 << "\n";
cout << "  Standard deviation << " << J3 << "\n";
//
//  Terminate.
//
cout << "\n";
cout << "thermal_design:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );