psn_common.c

Finite Volume Poisson PDE Solver

/*psn_lib/psn_common.c*/
/***********************************************************************
    Finite Volume Poisson PDE Solver: C-Library & Matlab Toolbox
    Implements numerical solution of Poisson PDE
    in 2D  Cartesian and Cylindrical coordinates

    Copyright (C) 2004 Igor Kaufman
    Copyright (C) 2004 Lancaster University, UK

    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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


    Author's email: [e-mail]

    APPLICATION    : COMMON ROUTINES
    VERSION        : 1.0
************************************************************************/
/*17.01.04*/
#include "psn.h"
#ifdef __MEX__
  #include "mex.h"
  #define psn_alloc(n,sz)  mxCalloc((n),(sz))
  #define psn_free(a)  mxFree(a)
#else
  #define psn_alloc(n,sz)  malloc((n)*(sz))
  #define psn_free(a)  free(a)
#endif

PSN_VECTOR psn_vector_create(size_t n)
{
   return (PSN_VECTOR)psn_alloc(n,sizeof(double));
}

void psn_vector_free(PSN_VECTOR a)
{
  if (a) psn_free(a);
}

PSN_MATRIX psn_matrix_create(size_t row_num, size_t col_num)
{
  size_t i, mn;
  double *a;
  PSN_MATRIX tmp =(PSN_MATRIX)psn_alloc(col_num,sizeof(double*));
#ifndef __MEX__
  for (i=0;i<col_num;i++) {
    tmp[i]=psn_vector_create(row_num);
  }
#else
  mn=col_num*row_num;
  a=(double*)psn_alloc(mn,sizeof(double));
  for (i=0;i<col_num;i++) {
    tmp[i]=&a[i*row_num];
  }
#endif
  return tmp;
}

#ifdef __MEX__
PSN_MATRIX psn_matrix_create_from_mxArray(const mxArray *ar)
{
  size_t i;
  size_t row_num=mxGetM(ar);
  size_t col_num=mxGetN(ar);
  double *data=mxGetPr(ar);
  PSN_MATRIX tmp=(PSN_MATRIX)psn_alloc(col_num,sizeof(double*));
  for (i=0;i<col_num;i++) {
    tmp[i]=&data[i*row_num];
  }
  return tmp;
}

int psn_matrix_free_from_mxArray(PSN_MATRIX a)
{
  if (a) psn_free(a);
  return 0;
}

#endif

void psn_matrix_free(PSN_MATRIX a, size_t col_num)
{
#ifndef __MEX__
  size_t i;
  for (i=0;i<col_num;i++) {
    psn_vector_free(a[i]);
  }
#endif
  if (a) psn_free(a);
}

void psn_matrix_get_rowcol(
           const size_t N,
           size_t num,
           int mode,   //1-row, 0-col
           PSN_MATRIX mat,
           PSN_VECTOR vec)
{
   size_t i;
   switch (mode) {
      case smCol:
        for(i=0;i<=N+1;i++)
         psn_vector_setval(vec,i,psn_matrix_getval(mat,i,num));
        break;
      case smRow:
        for(i=0;i<=N+1;i++)
          psn_vector_setval(vec,i,psn_matrix_getval(mat,num,i));
        break;
   }

}

void psn_matrix_set_rowcol(const size_t N,
           size_t num,
           int mode,
           PSN_MATRIX mat,
           PSN_VECTOR vec)
{
   size_t i;
   switch (mode)   {
     case smCol:
       for(i=0;i<=N+1;i++)
          psn_matrix_setval(mat,i,num,psn_vector_getval(vec,i));
       break;
     case smRow:
       for(i=0;i<=N+1;i++)
          psn_matrix_setval(mat,num,i,psn_vector_getval(vec,i));
       break;
   }
}

double norm_vector(size_t N, PSN_VECTOR oldvec, PSN_VECTOR newvec) {
  size_t i;
  double tmp=0, delta;
  for (i=1;i<=N;i++) {
     delta=fabs(psn_vector_getval(oldvec,i)-psn_vector_getval(newvec,i));
     tmp=max(tmp,delta);
  }
  return tmp;
}