femsolver.hpp

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

// 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 :  [e-mail]
//
// 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 <[e-mail]>
  @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