An Object Finite Element Library

SteklovPoincare2DBE.h

Go to the documentation of this file.

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

                                    O  F  E  L  I

                           Object  Finite  Element  Library

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

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

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

                       Definition of class SteklovPoincare2DBE
              for boundary elements for the 2-D Steklov Poincare solver
                          using P0 collocation elements

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


#ifndef __STEKLOV_POINCARE_2DBE_H
#define __STEKLOV_POINCARE_2DBE_H

#include "Mesh.h"
#include "DMatrix.h"
#include "Precond.h"
#include "GMRes.h"
#include "Vect.h"
#include "Triang3.h"

namespace OFELI {

#ifndef DOXYGEN_SHOULD_SKIP_THIS
struct ErrorInSteklovPoincare2DBE {
   void Message(const char *file, size_t line, int code);
};
#endif /* DOXYGEN_SHOULD_SKIP_THIS */


class SteklovPoincare2DBE {

public :

    SteklovPoincare2DBE(bool ext=false)
   { 
      _ext = -1;
      if (ext)
         _ext = 1; 
   }

    SteklovPoincare2DBE(const Mesh &mesh, bool ext=false);

    SteklovPoincare2DBE(const Mesh &mesh, const SideVect<double> &g, Vect<double> &b, bool ext=false);

    SteklovPoincare2DBE(const Mesh &mesh, const NodeVect<double> &g, Vect<double> &b, bool ext=false);

    SteklovPoincare2DBE(const Mesh &mesh, const SideVect<double> &g, SideVect<double> &b, bool ext=false);

    SteklovPoincare2DBE(const Mesh &mesh, const SideVect<double> &g, NodeVect<double> &b, bool ext=false);

    ~SteklovPoincare2DBE() { }

    void setMesh(const Mesh &mesh, bool ext=false);

    void Solve();

    int Solve(Vect<double> &b, const SideVect<double> &g);

    int Solve(Vect<double> &b, const NodeVect<double> &g);

private :
   int                        _ext;
   const Mesh                 *_theMesh;
   size_t                     _nb_eq;
   Vect<Point<double> >       _nn, _ttg, _center;
   Vect<double>               _length;
   double                     _h;
   DMatrix<double>            _A;
   ErrorInSteklovPoincare2DBE _e;

// Return integral of the Green function over the side sd at point x 
   inline double _single_layer(size_t j, const Point<double> &z) const
   {
      double t = z.NNorm();
      if (t < OFELI_EPSMCH)
         return -0.25/OFELI_PI*_h*_I1(0.5*_h);
      else
         return -0.125/OFELI_PI*_h*_I2(0.25*_h*_h,z*_ttg(j),t);
   }

// Return integral of the normal derivative of the Green function over the side sd at point x
   inline double _double_layer(size_t j, const Point<double> &z) const
   {
      double t = _nn(j)*z;
      if (fabs(t) < OFELI_EPSMCH)
         return 0;
      else
         return -0.25/OFELI_PI*t*_I3(0.25*_h*_h,z*_ttg(j),z.NNorm()); 
   }

// Return integral of log|a*t| over (-1,1)
   inline double _I1(double a) const { return 2*log(a)-2; }

// Return integral of log(a*t^2+b*t+c) over (-1,1) when 4*a*c-b^2>0.
   inline double _I2(double a, double b, double c) const
   {
      double d = sqrt(4*a*c-b*b);
      return 0.5*(log(a-b+c)*(2*a-b)+log(a+b+c)*(2*a+b)+2*(atan((2*a+b)/d)+atan((2*a-b)/d))*sqrt(4*a*c-b*b)-8*a)/a;
   }

// Return integral of 1/(a*t^2+b*t+c) over (-1,1) when 4*a*c-b^2>0.
   inline double _I3(double a, double b, double c) const
   {
      double d = sqrt(4*a*c-b*b);
      return 2*(atan((2*a-b)/d)+atan((2*a+b)/d))/d;
   }

// Calculate the normal to sides and their lengths
   void _util();
};

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani