//  Discussion:
//
//    This example solves equations related to a model of a microwave cooker.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/freefem++/microwave.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 February 2015
//
//  Author:
//
//    An example from the FreeFem++ manual.
//    Modifications and comments by John Burkardt.
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//
cout << "\n";
cout << "microwave\n";
cout << "  FreeFem++ version\n";
cout << "  Simulation of a microwave food processor.\n";
//
//  These parameters are used to define the boundaries.
//
real a = 20.0;
real b = 20.0;
real c = 15.0;
real d = 8.0;
real e = 2.0;
real l = 12.0;
real f = 2.0;
real g = 2.0;
//
//  Define the line segments A0, A1, A2, A3, A4, A5 of the external boundary.
//  Segment A4 (label = 2 ) is where the nonzero wave will be emitted.
//
//  Y
//  |
// 20  +--------A2------------+
//  |  |                      |
//  | A3                      |
//  |  |                      |
//  | A4                     A1
//  |  |                      |
//  | A5                      |
//  |  |                      |
//  0  +--------A0------------+
//  |
//  +--0---------------------20----X
//
border a0 (t=0,1) { x=a*t;       y=0.0;             label=1; }
border a1 (t=1,2) { x=a;         y=b*(t-1.0);       label=1; }
border a2 (t=2,3) { x=a*(3.0-t); y=b;               label=1; }
border a3 (t=3,4) { x=0.0;       y=b-(b-c)*(t-3.0); label=1; }
border a4 (t=4,5) { x=0.0;       y=c-(c-d)*(t-4.0); label=2; }
border a5 (t=5,6) { x=0.0;       y=d*(6.0-t);       label=1; }
//
//  Define the line segments B0, B1, B2, B3, of the internal interface.
//
//  Y
//  |
// 20  +----------------------+
//  |  |                      |
//  |  |                      |
// 14  |                +B2+  |
//  |  |               B3  |  |
//  |  |                | B1  |
//  2  |                +B0+  |
//  |  |                      |
//  0  +----------------------+
//  |
//  +--0---------------16-18-20----X
//
border b0 (t=0,1) { x=a-f+e*(t-1.0); y=g;                  label=3; }
border b1 (t=1,4) { x=a-f;           y=g+l*(t-1.0)/3.0;    label=3; }
border b2 (t=4,5) { x=a-f-e*(t-4.0); y=l+g;                label=3; }
border b3 (t=5,8) { x=a-e-f;         y= l+g-l*(t-5.0)/3.0; label=3; }
//
//  Compute a mesh by filling in the line segments with a given number of
//  points.  The value of N controls how fine the lines are.
//
int n = 2;
mesh Th = buildmesh ( 
    a0 ( 10 * n ) + a1 ( 10 * n) + a2 ( 10 * n ) + a3 ( 10 * n ) 
  + a4 ( 10 * n ) + a5 ( 10 * n )
  + b0 ( 5 * n ) + b1 ( 10 * n ) + b2 ( 5 * n ) + b3 ( 10 * n ) );
//
//  Plot the mesh.
//
plot ( Th, ps = "microwave_mesh.ps", wait = true );
//
//  Use a finite element space of piecewise linear functions on mesh Th.
//
fespace Vh ( Th, P1 );
//
//  The mesh is naturally divided into two regions.
//  We can request a "label" for a region by specifying a point inside that region.
//  The command 
//    "meat = Th ( xm, ym ).region" 
//  returns, as the value "meat", the label that FreeFem++ uses to identify
//  the portion of the mesh that contains (xm,ym), which happens to be the
//  region we think of as the meat.  Similarly, we can choose a point (xa,ya)
//  that lies in the other region, to get the label for the "air" portion
//  of our problem.  
//
real xm = 17.0;
real ym =  8.0;
real meat = Th ( xm, ym ).region;

real xa = 0.01;
real ya = 0.01;
real air = Th ( xa, ya ).region;
//
//  Now that we have the numeric labels "meat" and "air", we can
//  define R, a finite element interpolant function, which will be
//  1 in the meat region and 0 in the air region:
//
Vh R = ( region - air ) / ( meat - air );
//
//  Our unknown electromagnetic function is complex valued.
//  U1 will be our solution function, and V1 our test function.
//
Vh<complex> u1;
Vh<complex> v1;
//
//  Solve for the electromagnetic function U1.
//  The boundary condition is mostly zero, but over the boundary segment
//  with label 2, there is a real-valued sinusoidal boundary condition.
//
solve muwave ( u1, v1 ) = 
    int2d(Th) ( u1 * v1 * ( 1.0 + R ) 
                - ( dx(u1)*dx(v1)+dy(u1)*dy(v1)) * ( 1.0 - 0.5i) )
  + on ( 1, u1 = 0.0 ) 
  + on ( 2, u1 = sin(pi*(y-c)/(c-d)) );
//
//  Define the real and imaginary parts of U1, so we can plot them,
//  and so we can write a formula for the magnitude of the signal.
//
Vh u1r = real ( u1 );
Vh u1i = imag ( u1 );
//
//  Plot the real and imaginary portions of the solution.
//
plot ( u1r, wait = true, ps = "microwave_r.ps", fill = true );
plot ( u1i, wait = true, ps = "microwave_i.ps", fill = true );
//
//  Now we use the microwave solution to solve the second problem,
//  for the temperature distribution over the region.
//
fespace Uh ( Th, P1 ); 
//
//  U2 is our solution, and V2 is our test function.
//
Uh u2;
Uh v2;
//
//  The right hand side of the equations is FF, which is proportional to
//  the square of the signal intensity.  However, FF is only nonzero within
//  the meat region, so we multiply by R, the function that is 0 outside 
//  the meat.
//
Uh ff = 100000.0 * ( u1r^2 + u1i^2 ) * R;
//
//  Now we can solve the usual heat equation over the whole microwave
//  region, using a right hand side that is squelched outside the meat region.
//
solve temperature ( u2, v2 ) = 
    int2d(Th) ( dx(u2)* dx(v2) + dy(u2)* dy(v2) )
  - int2d(Th) ( ff * v2 ) 
  + on ( 1, u2 = 0.0 )
  + on ( 2, u2 = 0.0 );
//
//  Plot the temperature.
//
plot ( u2, wait = true, ps = "microwave_t.ps", fill = true );
//
//  Terminate.
//
cout << "\n";
cout << "microwave\n";
cout << "  Normal end of execution.\n";

exit ( 0 );