Source: k2dgrflow/galarsolver2d.h

/***************************************************************************
                          galarsolver2d.h  -  description
                             -------------------
    begin                : Tue May 29 2001
    copyright            : (C) 2001 by benny
    email                : bm@cage.rug.ac.be
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef GALARSOLVER2D_H
#define GALARSOLVER2D_H


/**class with the general outline of a problem
Every problem wich solves a FEM problem should be
derived from this abstract class. With this we solve generally the problem :

    Problem definition :
   div [-Adif(x,y)*grad u(x,y) - Acon(x,y)*u(x,y)] + Asou(x,y)*u(x,y) = Fsou(x,y) + div Fdiv(x,y)

                                                                                          u(x,y) = Gdir(x,y)  on Dirichlet boundary

      -Adif(x,y)*grad u(x,y) - Acon(x,y)*u(x,y) - Grob(x,y)*u(x,y) = Gneu(x,y)  on Neumann/Robin boundary

  *@author benny
  */

//matrix package  - a C package for LU decomposition
/* add constant to make sure the complex function clog does not conflict with iostream clog */
#define __BENNY_CLOG
extern "C" {
   #include "meschach/matrix.h"
   #include "meschach/matrix2.h"
   #include "meschach/sparse.h"
   #include "meschach/sparse2.h"
   #include "meschach/iter.h"
}

class DomainNet2d;
class TopoElement2d;
class TopoNode2d;
class TopoNodeDir2d;
class TopoNodeNeu2d;
class TopoNodeRobin2d;

class GalarSolver2d {
public:
     /** construct the solution class by giving at a minimum the domain on which to solve the equation
       */
	GalarSolver2d(DomainNet2d* );
	virtual ~GalarSolver2d();
    /** solve the problem and return a pointer to the solution.
         beware : when the problem is deallocated, also the memory where the solution resides is freed */
    virtual VEC* solution()=0;
    /** construct the global stiffness matrix, by adding the correct local stiffness matrices
         If done before, init global stiffness first ! */
    SPMAT* assembleStiffness();
    /** construct the global force vector, by adding the correct local force vectorterms per node
         If done before, init global force first !*/
    VEC* assembleForce();
    /** return the present value of the sparse matrix for node I (row) and J (column)
      */
    double stiffnessMatrixIJ(int globalI, int globalJ) {return sp_get_val(stiffness , globalI, globalJ) ;} ;
    /** return the present value of the force vector
      */
    double forceVectorJ(int globalJ ) {return force->ve[globalJ] ;} ;
    /** init the Stiffness quickly, necessary before you do a second assemble
      */
    void initStiffness();
    /** init the Force quickly, necessary before you do a second assemble
      */
    void initForce();
    /** return a pointer to the net on wich we solve a problem
      */
    DomainNet2d* net2d() {return net ;} ;
    /** print solution (unknowns) to a certain file, given by its name, eg /home/me/k2dgrflow/data.u
      */
    void printSolFile(const char *filename);
    /** print stiffnessmatrix A from A x = f to a certain file, given by its name, eg /home/me/k2dgrflow/data.a
      */
    void printStiffness(const char *filename);
    /** print forcevecor f from A x = f to a certain file, given by its name, eg /home/me/k2dgrflow/data.f
      */
    void printForce(const char *filename);

    /**a static method to easely read in vectors from files
      */
    static VEC *readVEC(FILE *fp, VEC *vec);

protected :
     /** pointer to the net on the domain
       */
     DomainNet2d *net;
     /** pointer to the sparse global stiffness matrix
       */
     SPMAT *stiffness;
     /** pointer to the global force vector
       */
     VEC *force;
     /** pointer to the vector with the unknowns. contains the solution when the problem is solved
       */
     VEC *unknowns;

      /** protected virtual methods that have to be implemented to solve a specific Galarking problem    */
      /** calculate the localstiffness matrix with a specific element
        */
      virtual void localStiffness(TopoElement2d* el) =0 ;
      /** calculate the local force vector with a specific element
        */
      virtual void localForce(TopoElement2d* el)=0 ;
      /** calculate the diffusion term for local node I,J in an element
        */
      virtual double diffusion(TopoElement2d* el, unsigned int localI, unsigned int localJ)=0 ;
      /** calculate the convection term for local node I,J in an element
        */
      virtual double convection(TopoElement2d* el, unsigned int localI, unsigned int localJ)=0 ;
      /** calculate the robin term for local node I,J in an element
        */
      virtual double robin(TopoElement2d* el, unsigned int localI, unsigned int localJ)=0 ;
      /** calculate the head term for local node I,J in an element
        */
      virtual double headTerm(TopoElement2d* el, unsigned int localI, unsigned int localJ)=0 ;
      /** calculate the source term for local node J in an element, part of the force
        */
      virtual double sourceTerm(TopoElement2d* el, unsigned int localJ)=0 ;
      /** calculate the neumann term for local node J in an element, part of the force
        */
      virtual double neumannTerm(TopoElement2d* el, unsigned int localJ)=0 ;
      /** return the gdir function in a node
        */
      virtual double gDir(TopoNodeDir2d* nod)=0 ;
      /** return the gNeu function in a node
        */
      virtual double gNeu(TopoNodeNeu2d* nod)=0 ;
      /** return the gRob function in a node
        */
      virtual double gRob(TopoNodeRobin2d* nod)=0 ;

private :

      /** some internal functions to speed things up
        */
      double sp_get_val_dir(SPMAT *A, int  i, int j) ;
      /** some internal functions to speed things up
        */
      double sp_set_val_dir(SPMAT *A, int  i, int j, double val) ;
      /** some internal functions to speed things up
        */
      double sp_set_val_dir_expand(SPMAT *A, int  i, int j, double val, SPROW * r, int idx) ;
};

inline double GalarSolver2d::sp_get_val_dir(SPMAT	*A, int i, int j)
{ //return value of sparse matrix without checking anything ! quicker
  SPROW	*r = A->row+i;
  int	idx( sprow_idx(r,j) );
  if ( idx < 0 )
     return 0.0;
   /* else */
   return r->elt[idx].val;
}
/* sp_set_val _dir-- sets the (i,j) entry of the sparse matrix A */
inline double GalarSolver2d::sp_set_val_dir(SPMAT *A, int i,int j, double val)
{
   SPROW	*r;

   r = A->row+i;
   int idx( sprow_idx(r,j));
   /* printf("sp_set_val: idx = %d\n",idx); */
   if ( idx >= 0 )
   {	r->elt[idx].val = val;	return val;	}
   /* else , not inlined*/
   return sp_set_val_dir_expand(A,i,j,val,r,idx);
}

#endif