#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 */ } }