updatemf.cpp

OFELI is an object oriented library of C++ classes for development of finite element codes. Its main

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

                                    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 "UpdateMF.h"

void UpdateMF(Mesh &ms0, Mesh &ms1, double omega, std::complex<double> &volt,
              Vect<std::complex<double> > &b0, Vect<std::complex<double> > &b1)
//------------------------------------------------------------------------------
//      Update magnetic field by taking into account boundary conditions
//------------------------------------------------------------------------------
{
   std::complex<double> a(0,omega);
   Element *el;
   Side    *sd;
   int i;

   std::complex<double> yy;
   for (ms0.TopElement(); (el=ms0.GetElement());) {
      EC2D1T3 eq(el,b0);
      yy += eq.IntegND();
   }
   std::complex<double> xx;
   for (ms1.TopElement(); (el=ms1.GetElement());) {
      EC2D1T3 eq(el,b1);
      xx += eq.IntegMF();
   }

// Calculate area of vacuum
   double aa=0.;
   for (ms0.TopSide(); (sd=ms0.GetSide());) {
      EC2D1T3 eq(sd);
      aa += eq.VacuumArea();
   }
   for (ms1.TopSide(); (sd=ms1.GetSide());) {
      EC2D1T3 eq(sd);
      aa -= eq.VacuumArea();
   }

   for (i=0; i<b0.Size(); i++)
      b0[i] *= volt/(a*(xx+MU0*aa)+yy);
   for (i=0; i<b1.Size(); i++)
      b1[i] *= volt/(a*(xx+MU0*aa)+yy);
}