oned_init.c

#include "array.h"
#include "poisson.h"

void oned_init(dmatrix *a, dmatrix *b, dmatrix *f, Decomposition oned)
{
  /* initialize all local matrices */

   int     p1_global_x, p1_global_y;
   int     p4_global_x, p4_global_y;
   double  offset;
   double rho;                             /* charge density */

   /* allocate arrays with ghostpoints  */
   new_dmatrix(a, oned.lx+2, oned.ly+2);
   new_dmatrix(b, oned.lx+2, oned.ly+2);
   new_dmatrix(f, oned.lx+2, oned.ly+2);

   /* 
    * set quadrupol sources
    */

   offset = (1 - QUAD_DIST)/2;
   rho    = (double) oned.gx * oned.gy;   /* total charge is constant */

   /* get global values */
   p1_global_x = (int) (offset * oned.gx) + 1;
   p1_global_y = (int) (offset * oned.gy) + 1;
   p4_global_x = (int) ((offset + QUAD_DIST) * oned.gx) + 1;
   p4_global_y = (int) ((offset + QUAD_DIST) * oned.gy) + 1;

   /* find local processor for P1 and P2 and set values */
   if ( (oned.llc_x <= p1_global_x) && (p1_global_x < oned.llc_x + oned.lx) )
   {
      INDEX(*f, p1_global_x - oned.llc_x + 1, p1_global_y) =  rho;   /* P1 */
      INDEX(*f, p1_global_x - oned.llc_x + 1, p4_global_y) = -rho;   /* P2 */
   }

   /* and now for P3 and P4 */
   if ( (oned.llc_x <= p4_global_x) && (p4_global_x < oned.llc_x + oned.lx) )
   {
      INDEX(*f, p4_global_x - oned.llc_x + 1, p1_global_y) = -rho;   /* P3 */
      INDEX(*f, p4_global_x - oned.llc_x + 1, p4_global_y) = +rho;   /* P4 */
   }
}