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

                                    O  F  E  L  I

                            Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2002 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"

#ifndef LARGE
#define LARGE 1.e20
#endif

using namespace OFELI;

void ScaledMF(Mesh &ms, SkMatrix<std::complex<double> > &a, Vect<std::complex<double> > &b, 
              std::complex<double> &current, double omega, int flag)
//------------------------------------------------------------------------------
//              Calculate scaled magnetic field
//------------------------------------------------------------------------------
{
   Element *el;
   Node *nd;
   for (ms.TopElement(); (el=ms.GetElement());) {
      EC2D1T3 equa(el);
      equa.Magnetic(omega,1.);
      equa.Electric();
      a.Assembly(el,equa.A());
   }

   for (ms.TopNode(); (nd=ms.GetNode());) {
      int m = nd->DOF(1);
      if (nd->Code(1)==1) {
        a.Set(m,m,a(m,m)*LARGE);
        b[m-1] = a(m,m)*current;
      }
      if (nd->Code(1)==2) {
        a.Set(m,m,a(m,m)*LARGE);
        b[m-1] = 0;
      }
   }

   a.Factor();
   a.Solve(b);
}