//  Discussion:
//
//    Set up a region with four internal boxes where the solution jumps.
//
//    Export the data, using ffmatlib, for graphic processing by MATLAB or Octave.
//
//  Licensing:
//
//    Copyright (C) 2018 Chloros2 <chloros2@gmx.de>
//
//    This program is free software: you can redistribute it and/or modify it
//    under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hopeC that it will be useful, but
//    WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this program.  If not, see
//    <https://www.gnu.org/licenses/>.
//
//  Modified:
//
//    2018-05-19
//
// Author: 
//
//   Chloros2 <chloros2@gmx.de>
//

include "ffmatlib.idp"

cout << endl;
cout << "region:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Set up a mesh on an interesting region, with four internal boxes" << endl;
cout << "  where the solution jumps." << endl;
cout << "  Using the ffmatlib interface, export the data for graphic" << endl;
cout << "  processing by MATLAB/Octave." << endl;

real[int] XX = [-1,1,-1,1];
real[int] YY = [-1,-1,1,1];

real[int] PHI = [ 0.0, 0.5, 0.6, 0.7 ];

int[int] CC = [ 1, 2, 3, 4 ];

real r0 = 0.2;
real ra = 3.0;
border outer ( t = 0, 2 * pi ) 
  { x = ra * cos ( t ); y = ra * sin ( t ); label = 123; }

border cycloides(t=0,2*pi;i)
  { x = r0 * ( 3.0 * cos ( t + PHI[i] ) + 0.3 * cos ( 3.0 * t ) ) + XX[i];
    y = r0 * ( 3.0 * sin ( t + PHI[i] ) - 0.3 * sin ( 3.0 * t ) ) + YY[i];
    label = CC[i]; }

int[int] NN = 40 * [ -1, 1, 1, 1 ];

mesh Th = buildmesh ( cycloides(NN) + outer(100) );

plot ( Th, wait = false,
  ps = "region_mesh.ps" );
//
//  Use piecewise constant space.
//
fespace Vh ( Th, P0 );

int reg1 = Th(XX[2],YY[2]).region;
int reg2 = Th(XX[1],YY[1]).region;

Vh u = 2.0 * ( region == reg1 ) - 1.0 * ( region == reg2 );

plot ( u, fill = true, wait = false,
  ps = "region_values.ps" );
//
//  Write out data for later processing by MATLAB/Octave.
//
savemesh ( Th, "region.msh" );
ffSaveVh ( Th, Vh, "region_vh.txt" );
ffSaveData ( u, "region_data.txt" );
//
//  Terminate.
//
cout << "\n";
cout << "region:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );