/*==============================================================================

               A Finite Element Code to Solve Helmholtz Equation
                      in a Bounded Domain using OFELI

  ------------------------------------------------------------------------------

   Copyright (C) 1998 - 2004 Rachid Touzani

   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; Version 2 of the License.

   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, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================*/

#include "OFELI.h"
#include "Electromagnetics.h"
#include "User.h"
using namespace OFELI;

int main(int argc, char *argv[])
{
   Mesh     ms;
   Element  *el;
   Side     *sd;
   ifstream mf, bcf;

   if (argc < 2) {
     cout << "\nUsage:  helmholtz <parameter_file>\n";
     return 0;
   }

   IPF data("Helmholtz - 1.0",argv[1]);
   int output_flag = data.Output();
   int bc_flag = data.BC();
   double wave_nb = data.DoublePar(1);

   if (output_flag) {
     cout << endl << endl;
     cout << "    *******************************************************\n";
     cout << "    *                      H e l m h o l t z              *\n";
     cout << "    *            Helmholtz equation in a bounded domain   *\n";
     cout << "    *******************************************************\n\n\n";
     cout << "=====================================================================\n\n";
     cout << "               A Finite Element Code for Helmholtz equation\n";
     cout << "                        in a 2-D bounded domain\n\n";
     cout << "            Helmholtz uses OFELI Library of Finite Element Classes\n\n";
     cout << "                           V E R S I O N   1.0\n\n";
     cout << "                     Copyright R. Touzani, 1999\n\n";
     cout << "=====================================================================\n\n";
   }

//----------
// Read data
//----------

// Read Mesh data
   if (output_flag > 1)
     cout << "Reading mesh data ...\n";
   ms.Get(data.MeshFile());
   int nb_dof = 1;
   wave_nb = data.DoublePar(1);

   if (output_flag > 1)
     cout << ms;
   User ud(ms);

// Declare problem data (matrix, rhs, boundary conditions, body forces)
   if (output_flag > 1)
     cout << "Allocating memory for matrix and R.H.S. ...\n";
   SkMatrix<std::complex<double> > a(ms);
   Vect<std::complex<double> > b(ms.NbDOF());

// Read boundary conditions, body and boundary forces
   if (output_flag > 1)
     cout << "Reading boundary conditions ...\n";
   Vect<std::complex<double> > bc(ms.NbDOF());
   if (!bc_flag)
     ud.SetDBC(bc);

   NodeVect<std::complex<double> > uf(ms,nb_dof);

// Read in boundary conditions and body forces
   ud.SetDBC(bc);

// Loop over elements
// ------------------

   for (ms.TopElement(); (el=ms.GetElement());) {
      HelmholtzBT3 eq(el);
      eq.LHS(wave_nb);
      a.Assembly(el,eq.A());
   }

// Loop over sides
// ---------------

   for (ms.TopSide(); (sd=ms.GetSide());) {
      HelmholtzBT3 eq(sd);
      eq.BoundaryRHS(ud);
      b.Assembly(sd,eq.b());
   }

   a.Prescribe(ms,b,bc);
   a.Factor();
   a.Solve(b);

   uf.FromVect(b,1,"Magnetic Field",0);
   if (output_flag > 0)
     cout << uf;

   void error(const Mesh &ms, User &ud, const Vect<complex<double> > &u);
   error(ms,ud,b);

   return 0;
}