libfreefem-dev 3.5.8-5.1+b2 / usr / include / fem / femSolver.hpp

// Emacs will be in -*- Mode: c++ -*-
//
// ********** DO NOT REMOVE THIS BANNER **********
//
// SUMMARY: Language for a Finite Element Method
// RELEASE: 
// AUTHORS:  C. Prud'homme
// ORG    :          
// E-MAIL :  prudhomm@users.sourceforge.net
//
// ORIG-DATE:     June-94
// LAST-MOD: 12-Jul-01 at 10:01:47 by 
//
// DESCRIPTION:
//   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 2 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, write to the Free Software
//   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
// DESCRIP-END.
//
#ifndef __FEM_H
#ifdef __GNUG__
#pragma interface
#endif
#define __FEM_H 1

#define nhowmax 20

#include <femCommon.hpp>

#include <femMisc.hpp>

namespace fem
{
DECLARE_CLASS( FEM );
DECLARE_CLASS( femMesh );

/*!
  \class FEM
  \brief this class drives the resolution of the pde using the Finite Element Method.
  
  @author Christophe Prud'homme <Christophe.Prudhomme@ann.jussieu.fr>
  @see femMesh
  @version #$Id: femSolver.hpp,v 1.2 2001/07/12 14:11:57 prudhomm Exp $#
*/
class FEM
{
public:

  /** Typedefs
   */
  //@{
  DECLARE_TYPE( femMesh::femPoint, femPoint );
  DECLARE_TYPE( femMesh::femTriangle, femTriangle );
  
  DECLARE_TYPE( creal*, cmatptr );
  DECLARE_TYPE( float*, matptr );
  //@}
  
  /** Constructors, destructor and methods
   */ 
  //@{

  //! default constructor
  FEM( femMeshPtr = 0, int quadra = 0 );

  //! destructor
  ~FEM();

  /**
     \brief solve the PDE
     \param param contain all the possible data for computation
     \param how defines if the P1 quadrature
   */
  float solvePDE( fcts *param, int how);

  /**
   */
  creal deriv (int m, creal * f, int ksolv, int i);

  /**
   */
  creal convect(creal* f, creal* u1, creal* u2, float dt, int i);

  /**
   */
  creal rhsConvect(creal* f, creal* u1, creal* u2, float dt, int i);

  /**
   */
  creal fctval(creal* f,float x,float y);

  /**
     @return the ngt of a femTriangle to which belongs vertex k
   */
  int getregion(int k);
  
  creal gfemuser( creal what, creal *f, int i);
  creal P1ctoP1(creal* f, int i);
  creal prodscalar(creal* f, creal* g);
  
  creal ginteg(int, int, int, creal*, creal*, int);
  creal binteg(int, int, int, creal*, creal*, int);
  
  void initvarmat ( int how, int flagcomplexe,int N, fcts* param);
  void assemble ( int how, int flagcomplexe,int N, int k, creal* a, creal* b, fcts* param);
  void solvevarpde(int N,  fcts* param, int how);

  //@}

private:

  float norm( float, float ) const;
  
  int doedge (void);
  int Tconvect (int k, double u1k, double u2k, double x, double y, double *mu, double *nu);
  int xtoX (creal * u1, creal * u2, float *dt, float *x, float *y, int *k);
  int barycoor (float x, float y, int k, float *xl0, float *xl1, float *xl2);
  void rhsPDE(Acvect & fw, Acvect & f, Acvect & g);
  int buildarea();
  void connectiv (void);
  
  /** Acmat routines
   */
  cmat id(cvect& x);
  
  void pdemat (Acmat & a, Acmat & alpha,
               Acmat & rho11, Acmat & rho12, Acmat & rho21, Acmat & rho22,
               Acmat & u1, Acmat & u2, Acmat & beta);
  float gaussband (Acmat & a, Acvect & x, long n, long bdth, int first, float eps);
  float pdeian (Acmat & a, Acvect & u, Acvect & f, Acvect & g, Acvect & u0,
                Acmat & alpha, Acmat & rho11, Acmat & rho12, Acmat & rho21, Acmat & rho22,
                Acmat & u1, Acmat & u2, Acmat & beta, int factorize);


  /** Float routines
   */
  float id(float x);

  /**
     Computes the right hand side of the linear system of the Laplace equation
     $fw(j)=\int_\Omega f w^j + \int_{ng<>0} g w^j + penal|_{u0!=0} u0^j$ 
     if quadra = i then f is P^i, i=0 or 1
     OTHER INPUT next,area,q,me,ng,nt,penal
  */
  void rhsPDE(float*  fw, float*  f, float*  g);

  /**
     Factorise (first!=0) and/or solve    Ay = x  with result in x 
     LU is stored in A ; returns the value of the smallest pivot  
     all pivots less than eps are put to eps 
     a[i][j] is stored in a[i-j+bdth][j]=a[n*(i-j+bdth)+j) 
     where -bdwth <= i-j <= bdth
  */
  float gaussband (float*  a, float*  x, long n, long bdthl, int first, float eps);

  void pdemat (float*  a, float*  alpha,
               float*  rho11, float*  rho12, float*  rho21, float*  rho22,
               float*  u1, float*  u2, float*  beta);
  /**
    Solves  alpha u + (u1,u2)grad u- div(mat[rho] grad u) = f in ã, 
    u|_{ng<0} = u0, 
    beta u + mat[rho]¦u/¦n|_{ng<>0} = g
    where f is P^quadra (quadra=0 or 1) 
  */
  float pdeian(float*  a, float*  u, float*  f, float*  g, float*  u0, 
               float*  alpha, float*  rho11, float*  rho12, float*  rho21, float*  rho22,
               float*  u1, float*  u2, float*  beta, int factorize);

  /** creal(complex) routines
   */
  void rhsPDE(creal*  fw, creal*  f, creal*  g);
  
  float gaussband (creal*  a, creal*  x, long n, long bdthl, int first, float eps);
  
  float pdeian(creal*  a, creal*  u, creal*  f, creal*  g, creal*  u0, 
               creal*  alpha, creal*  rho11, creal*  rho12, creal*  rho21, creal*  rho22,
               creal*  u1, creal*  u2, creal*  beta, int factorize);

  void pdemat (creal*  a, creal*  alpha, 
               creal*  rho11, creal*  rho12, creal*  rho21, creal*  rho22,
               creal*  u1, creal*  u2, creal*  beta);
  
  /**
     integration routines
   */
  creal  binteg_t (int k, int ref1, int ref2, int ref3, creal* f, creal* g);
  creal  ginteg_t (int k, creal* f, creal* g);

public:

  float* normlx;
  float* normly;

private:
  
  femMeshPtr __mesh;
  int __quadra;
  
  int              ns, nt;
  femPoint     *q;
  femTriangle  *me;
  int            *ng;
  int            *ngt;

  int             bug;
  int             nhow;
  int             nhow1;
  int             nhow2;

  struct
  {
    int             N;
    int             nhow;
  }whatmat[nhowmax];
  
  float          *area;
  AAcmat          a2;
  creal          *a1c[nhowmax];
  creal          *gconvect;
  int             rhsQuadra;
  float          *a1[nhowmax];
  int            *triauptr;
  int            *triaunb;
  long            bdth;

  int            *Tleft;
  int            *Tright;
  int            *lowV;
  int            *highV;
  int            *listHead;
  femTriangle       *edgeT;
  int             ne;

};
}
#endif /* __FEM_H */