//  poisson_ell.edp
//
//  Discussion:
//
//    Solve the Poisson equation over the L-shaped region.
//
//    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 hope 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 << "poisson_ell():" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Solve the Poisson equation over the L-shaped region." << endl;
cout << "  Using the ffmatlib interface, export the data for graphic" << endl;
cout << "  processing by MATLAB/Octave." << endl;

verbosity = 0;
border a ( t = 0.0, 1.0 ) { x=t;y=0;label=1; };
border b ( t = 0.0, 0.5 ) { x=1;y=t;label=1; };
border c ( t = 0.0, 0.5 ) { x=1-t;y=0.5;label=1; };
border d ( t = 0.5, 1.0 ) { x=0.5;y=t;label=1; };
border e ( t = 0.5, 1.0 ) { x=1-t;y=1;label=1; };
border f ( t = 0.0, 1.0 ) { x=0;y=1-t;label=1; };

int nn = 40;
mesh Th = buildmesh ( a(nn) + b(.5*nn) + c(.5*nn) +d(.5*nn) + e(.5*nn) + f(nn));

fespace Vh ( Th, P2 );

Vh u;
Vh v;

func rho = 1.0;

solve Poisson(u,v,solver=LU) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
                           - int2d(Th)(rho*v)+on(a,b,c,d,e,f,u=0);

plot ( Th, u, 
  ps = "poisson_ell.ps" );

cout << endl;
cout << "NbBoundaryElements:\t" <<  Th.nbe << endl;
cout << "NbTriangles:\t\t" <<  Th.nt << endl;
cout << "NbVertices:\t\t" <<  Th.nv << endl;
cout << "nDoF:\t\t\t" << Vh.ndof << endl;
cout << endl;

real lambda = 1.0;
Vh qx, qy;
qx = -lambda*dx(u);
qy = -lambda*dy(u);
//
//  Write out data for later processing by MATLAB/Octave.
//
savemesh ( Th, "poisson_ell.msh" );
ffSaveVh ( Th, Vh, "poisson_ell_vh.txt" );
ffSaveData3 ( u, qx, qy, "poisson_ell_data.txt" );
//
//  Terminate.
//
cout << "\n";
cout << "poisson_ell():\n";
cout << "  Normal end of execution.\n";

exit ( 0 );