muldisplay.c

很经典的电磁计算（电容计算）软件 MIT90年代开发

/*
Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
All rights reserved.

This Agreement gives you, the LICENSEE, certain rights and obligations.
By using the software, you indicate that you have read, understood, and
will comply with the terms.

Permission to use, copy and modify for internal, noncommercial purposes
is hereby granted.  Any distribution of this program or any part thereof
is strictly prohibited without prior written consent of M.I.T.

Title to copyright to this software and to any associated documentation
shall at all times remain with M.I.T. and LICENSEE agrees to preserve
same.  LICENSEE agrees not to make any copies except for LICENSEE'S
internal noncommercial use, or to use separately any portion of this
software without prior written consent of M.I.T.  LICENSEE agrees to
place the appropriate copyright notice on any such copies.

Nothing in this Agreement shall be construed as conferring rights to use
in advertising, publicity or otherwise any trademark or the name of
"Massachusetts Institute of Technology" or "M.I.T."

M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.  By
way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
M.I.T. shall not be held liable for any liability nor for any direct,
indirect or consequential damages with respect to any claim by LICENSEE
or any third party on account of or arising from this Agreement or use
of this software.
*/

#include "mulGlobal.h"
#include "zbufGlobal.h"

disExtrasimpcube(pc)
cube *pc;
{
  printf("cubes[%d][%d][%d][%d]\n", pc->level, pc->j, pc->k, pc->l);
}

disExParsimpcube(pc)
cube *pc;
{
  cube *pa = pc->parent;
  printf("cubes[%d][%d][%d][%d], ", pc->level, pc->j, pc->k, pc->l);
  printf("parent = cubes[%d][%d][%d][%d]\n", pa->level, pa->j, pa->k, pa->l);
}

dissimpcube(pc)
cube *pc;
{
int i;
  printf("cube center: x=%g y=%g z=%g\n", pc->x, pc->y, pc->z);
  printf("index=%d dindex=%d level=%d loc_exact=%d mul_exact=%d numkids=%d\n",
	 pc->index, pc->dindex, pc->level,
	 pc->loc_exact, pc->mul_exact, pc->numkids);
  printf("numnbrs=%d upnumvects=%d directnumvects=%d downnumvects=%d\n",
	 pc->numnbrs, pc->upnumvects, pc->directnumvects, pc->downnumvects);
}

discube(pc)
cube *pc;
{
int i;
  printf("cube center: x=%g y=%g z=%g\n", pc->x, pc->y, pc->z);
  printf("index=%d dindex=%d level=%d loc_exact=%d mul_exact=%d numkids=%d\n",
	 pc->index, pc->dindex, pc->level,
	 pc->loc_exact, pc->mul_exact, pc->numkids);
  printf("numnbrs=%d upnumvects=%d directnumvects=%d downnumvects=%d\n",
	 pc->numnbrs, pc->upnumvects, pc->directnumvects, pc->downnumvects);
  if(pc->directnumvects > 0) {
    printf("num of elements in ");
    for(i=0; i < pc->directnumvects; i++) {
      printf("v%d = %d ", i, pc->directnumeles[i]);
    }
    printf("\nchgs\n");
    for(i=0; i < pc->directnumeles[0]; i++) {
      dischg(pc->chgs[i]);
    }
  }
  if(pc->downnumvects > 0) {
    printf("num of down elements in ");
    for(i=0; i < pc->downnumvects; i++) {
      printf("v%d = %d ", i, pc->downnumeles[i]);
    }
  }
}

disupcube(pc)
cube *pc;
{


}

disdirectcube(pc)
cube *pc;
{
int i;
  for(i=0; i < pc->directnumvects; i++) {
    printf("matrix %d\n", i);
    dismat(pc->directmats[i], pc->directnumeles[0], pc->directnumeles[i]);
    if(i==0) {
      printf("lu factored matrix\n");
      dismat(pc->directlu, pc->directnumeles[0], pc->directnumeles[i]);
    }
  }
}


dissys(sys)
ssystem *sys;
{
int i, j, k, l, side;
  printf("side=%d depth=%d order=%d\n",
	 sys->side, sys->depth, sys->order);
  printf("Cube corner is x=%g y=%g z=%g\n", sys->minx, sys->miny, sys->minz);
  printf("Cube side length= %g\n", sys->length);
  printf("Printing all the cubes\n");
  for(i = 0, side = 1; i <= sys->depth; i++, side *= 2) {
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  fprintf(stdout, "\ncubes[%d][%d][%d][%d]\n", i, j, k, l);
	  dissimpcube(&(sys->cubes[i][j][k][l]));
/*	  disdirectcube(&(sys->cubes[i][j][k][l])); */
	}
      }
    }
  }
}



dismat(mat, rows, cols)
double **mat;
int rows, cols;
{
int i,j;
  if(cols != 0) {
    for(i=0; i < rows; i++) {
      printf("\n i=%d\n", i);
      for(j=0; j < cols; j++) {
        printf("%d %g  ", j, mat[i][j]);
        if(((j+1) % 5) == 0) printf("\n");
      }
    }
    printf("\n");
  }
}


disvect(v, size)
double *v;
int size;
{
int i;
  for(i=0; i < size; i++) {
    printf("i=%d %g ", i, v[i]);
    if(((i+1) % 5) == 0) printf("\n");
  }
  printf("\n");
}

dischg(pq)
charge *pq;
{
  printf("cond=%d index=%d\n", pq->cond, pq->index);
}

disallchg(pq) 
charge *pq;
{
charge *nq;
  for(nq = pq; nq != NULL; nq = nq->next) disfchg(pq);
}

disfchg(pq) 
charge *pq;
{
/*
  printf("Cond=%d Corners\n", pq->cond);
  printf("x0=%g y0=%g z0=%g\n", pq->x0, pq->y0, pq->z0);
  printf("x1=%g y1=%g z1=%g\n", pq->x1, pq->y1, pq->z1);
  printf("x2=%g y2=%g z2=%g\n", pq->x2, pq->y2, pq->z2);
  printf("x3=%g y3=%g z3=%g\n", pq->x3, pq->y3, pq->z3);
  printf("Center\n");
  printf("x=%g y=%g z=%g\n", pq->x, pq->y, pq->z);
*/
}

/*
  dumps a rows x cols matrix of doubles; assumes indices from zero 
*/
void dumpMat(mat, rows, cols)
int rows, cols;
double **mat;
{
  int i, j;
  for(i = 0; i < rows; i++) {
    fprintf(stdout, "    row%d ", i);
    for(j = 0; j < cols; j++) {
      if(mat[i][j] < 0.0) fprintf(stdout, "%.5e ", mat[i][j]);
      else fprintf(stdout, " %.5e ", mat[i][j]);
    }
    fprintf(stdout, "\n");
  }
}

/*
  dumps a rows x cols matrix of doubles; assumes indices from zero 
*/
void dumpCorners(fp, mat, rows, cols)
int rows, cols;
double **mat;
FILE *fp;
{
  int i, j;
  for(i = 0; i < rows; i++) {
    fprintf(fp, "  corner%d ", i);
    for(j = 0; j < cols; j++) {
      if(mat[i][j] < 0.0) fprintf(fp, "%.5e ", mat[i][j]);
      else fprintf(fp, " %.5e ", mat[i][j]);
    }
    fprintf(fp, "\n");
  }
}

/*
  dumps a vector of itegers along side a vector of doubles, index from zero
*/
void dumpVecs(dblvec, intvec, size)
double *dblvec;
int *intvec, size;
{
  int i;

  for(i = 0; i < size; i++) {
    fprintf(stdout, "%d %d %g\n", i, intvec[i], dblvec[i]);
  }
}

/*
  dumps the relative coordinates of an array of charges or evaluation pnts
*/
void dumpChgs(chgs, numchgs, x, y, z)
int numchgs;
double x, y, z;
charge **chgs;
{
  int i;
  double rho, cosA, beta;
  for(i = 0; i < numchgs; i++) {
    xyz2sphere(chgs[i]->x, chgs[i]->y, chgs[i]->z,
	       x, y, z, &rho, &cosA, &beta);
    fprintf(stdout, "    %d %d ", chgs[i]->index, chgs[i]->cond);
    if(rho < 0) fprintf(stdout, "(%.5e ", rho);
    else fprintf(stdout, "( %.5e ", rho);
    if(cosA < 0) fprintf(stdout, "%.5e ", cosA);
    else fprintf(stdout, " %.5e ", cosA);
    if(beta < 0) fprintf(stdout, "%.5e) ", beta);
    else fprintf(stdout, " %.5e) ", beta);
    if(x < 0) fprintf(stdout, "(%.5e ", chgs[i]->x);
    else fprintf(stdout, "( %.5e ", chgs[i]->x);
    if(y < 0) fprintf(stdout, "%.5e ", chgs[i]->y);
    else fprintf(stdout, " %.5e ", chgs[i]->y);
    if(z < 0) fprintf(stdout, "%.5e)\n", chgs[i]->z);
    else fprintf(stdout, " %.5e)\n", chgs[i]->z);
  }
}

/*
  dumps the relative coordinates of an array of charges or evaluation pnts
  - also dumps dummy bit
*/
void dumpChgsWDummy(chgs, numchgs, is_dummy, x, y, z)
int numchgs, *is_dummy;
double x, y, z;
charge **chgs;
{
  int i;
  double rho, cosA, beta;
  for(i = 0; i < numchgs; i++) {
    xyz2sphere(chgs[i]->x, chgs[i]->y, chgs[i]->z,
	       x, y, z, &rho, &cosA, &beta);
    fprintf(stdout, "    %d %d(%d) %d ", chgs[i]->index, is_dummy[i], 
	    chgs[i]->dummy, chgs[i]->cond);
    if(rho < 0) fprintf(stdout, "(%.5e ", rho);
    else fprintf(stdout, "( %.5e ", rho);
    if(cosA < 0) fprintf(stdout, "%.5e ", cosA);
    else fprintf(stdout, " %.5e ", cosA);
    if(beta < 0) fprintf(stdout, "%.5e) ", beta);
    else fprintf(stdout, " %.5e) ", beta);
    if(x < 0) fprintf(stdout, "(%.5e ", chgs[i]->x);
    else fprintf(stdout, "( %.5e ", chgs[i]->x);
    if(y < 0) fprintf(stdout, "%.5e ", chgs[i]->y);
    else fprintf(stdout, " %.5e ", chgs[i]->y);
    if(z < 0) fprintf(stdout, "%.5e)\n", chgs[i]->z);
    else fprintf(stdout, " %.5e)\n", chgs[i]->z);
  }
}

/*
  display the matrix built for a given charge to multipole transformation
*/
void dispQ2M(mat, chgs, numchgs, x, y, z, order)
int numchgs, order;
double **mat, x, y, z;
charge **chgs;
{
  fprintf(stdout, "\nQ2M MATRIX: cube at (%.5e %.5e %.5e)\n", x, y, z);
  dumpMat(mat, multerms(order), numchgs);
  fprintf(stdout, 
	  "    CHARGES IN CUBE # cond (rho_i cos(alpha_i) beta_i) (x y z):\n");
  dumpChgs(chgs, numchgs, x, y, z);
}

/*
  display the matrix built for a given multipole to local transformation
*/
void dispM2L(mat, x, y, z, xp, yp, zp, order)
int order;
double **mat, x, y, z, xp, yp, zp;
{
  fprintf(stdout, 
   "\nM2L MATRIX: multi at (%.5e %.5e %.5e) -> local at (%.5e %.5e %.5e)\n",
	  x, y, z, xp, yp, zp);
  dumpMat(mat, multerms(order), multerms(order));
}

/*
  display the matrix built for a given charge to local transformation
*/
void dispQ2L(mat, chgs, numchgs, x, y, z, order)
int numchgs, order;
double **mat, x, y, z;
charge **chgs;
{
  fprintf(stdout, "\nQ2L MATRIX: cube at (%.5e %.5e %.5e)\n", x, y, z);
  dumpMat(mat, multerms(order), numchgs);
  fprintf(stdout, 
	  "    CHARGES IN CUBE # cond (rho_i cos(alpha_i) beta_i) (x y z):\n");
  dumpChgs(chgs, numchgs, x, y, z);
}

/*
  display the matrix built for a given charge to potential transformation
*/
void dispQ2P(mat, chgs, numchgs, is_dummy, pchgs, numpchgs)
int numchgs, numpchgs, *is_dummy;
double **mat;
charge **chgs, **pchgs;
{
  fprintf(stdout, "\nQ2P MATRIX:\n");
  dumpMat(mat, numpchgs, numchgs);
  fprintf(stdout, 
	  "    PANELS IN CUBE # dummy(real) cond (rho_i cos(alpha_i) beta_i) (x y z):\n");
  dumpChgsWDummy(chgs, numchgs, is_dummy, 0.0, 0.0, 0.0);
  fprintf(stdout, 
	  "    EVALS IN CUBE # cond (rho_i cos(alpha_i) beta_i) (x y z):\n");
  dumpChgs(pchgs, numpchgs, 0.0, 0.0, 0.0);
}

/*
  display the matrix built for a given charge to potential transformation
*/
void dispQ2PDiag(mat, chgs, numchgs, is_dummy)
int numchgs, *is_dummy;
double **mat;
charge **chgs;
{
  fprintf(stdout, "\nQ2PDiag MATRIX:\n");
  dumpMat(mat, numchgs, numchgs);
  fprintf(stdout, 
	  "    PANELS IN CUBE # dummy(real) cond (rho_i cos(alpha_i) beta_i) (x y z):\n");
  dumpChgsWDummy(chgs, numchgs, is_dummy, 0.0, 0.0, 0.0);
}

/*
  display the matrix built for a given multipole to multipole transformation
*/
void dispM2M(mat, x, y, z, xp, yp, zp, order)
int order;
double **mat, x, y, z, xp, yp, zp;
{
  fprintf(stdout, 
      "\nM2M MATRIX: cube at (%.5e %.5e %.5e) shifted to (%.5e %.5e %.5e)\n", 
	  x, y, z, xp, yp, zp);
  dumpMat(mat, multerms(order), multerms(order));
}

/*
  display the matrix built for a given local to local transformation
*/
void dispL2L(mat, x, y, z, xp, yp, zp, order)
int order;
double **mat, x, y, z, xp, yp, zp;
{
  fprintf(stdout, 
      "\nL2L MATRIX: cube at (%.5e %.5e %.5e) shifted to (%.5e %.5e %.5e)\n", 
	  x, y, z, xp, yp, zp);
  dumpMat(mat, multerms(order), multerms(order));
}

/*
  display the matrix built for a given multipole to potential transformation
*/
void dispM2P(mat, x, y, z, chgs, numchgs, order)
int numchgs, order;
double **mat, x, y, z;
charge **chgs;
{
  fprintf(stdout, "\nM2P MATRIX: cube at (%.5e %.5e %.5e)\n", x, y, z);
  dumpMat(mat, numchgs, multerms(order));
  fprintf(stdout, 
	  "    EVAL PNTS IN CUBE # cond (rho_i, cos(alpha_i), beta_i):\n");
  dumpChgs(chgs, numchgs, x, y, z);
}

/*
  display the matrix built for a given local to potential transformation
*/
void dispL2P(mat, x, y, z, chgs, numchgs, order)
int numchgs, order;
double **mat, x, y, z;
charge **chgs;
{
  fprintf(stdout, "\nL2P MATRIX: cube at (%.5e %.5e %.5e)\n", x, y, z);
  dumpMat(mat, numchgs, multerms(order));
  fprintf(stdout, 
	  "    EVAL PNTS IN CUBE # cond (rho_i, cos(alpha_i), beta_i):\n");
  dumpChgs(chgs, numchgs, x, y, z);
}

/*
  displays upward pass and moment vectors associated with a cube - debug only
*/
void dumpUpVecs(pc)
cube *pc;
{
  int i, j;
  fprintf(stdout, 
    "\nUPWARD PASS/MOMENT VECTORS, LEVEL %d CUBE AT (%.5e %.5e %.5e):\n",
	  pc->level, pc->x, pc->y, pc->z);
  for(i = 0; i < pc->upnumvects; i++) {
    fprintf(stdout, "%d", i);
    for(j = 0; j < pc->upnumeles[i]; j++) {
      if(pc->upvects[i][j] < 0.0) 
	  fprintf(stdout, " %.5e", pc->upvects[i][j]);
      else fprintf(stdout, "  %.5e", pc->upvects[i][j]);
    }
    fprintf(stdout, "\n");
  }
  fprintf(stdout, "M");
  for(j = 0; j < pc->multisize; j++) {
    if(pc->multi[j] < 0.0) fprintf(stdout, " %.5e", pc->multi[j]);
    else fprintf(stdout, "  %.5e", pc->multi[j]);
  }
  fprintf(stdout, "\n");
}

/*
  displays the upward pass vectors for the eight level 1 cubes - debug only
*/
void dumpLevOneUpVecs(sys)
ssystem *sys;
{
  int i, j, k;
  cube *****cubes = sys->cubes;
  for(i = 0; i < 2; i++) {
    for(j = 0; j < 2; j++) {
      for(k = 0; k < 2; k++) {
	if(cubes[1][i][j][k] != NULL) dumpUpVecs(cubes[1][i][j][k]);
      }
    }
  }
}

/*
  checks a cube (direct, local or eval) list for bad cube structs - debug only
  -- doesn't quite do this - always uses direct list for one thing
*/
void chkList(sys, listtype)
ssystem *sys;
int listtype;			/* DIRECT, LOCAL or EVAL */
{
  int cnt[BUFSIZ];		/* # of cubes processed by level */
  int depth = sys->depth;
  int lev, nn;
  int i, j, k;
  cube *nc;
  for(i = 0; i <= depth; i++) cnt[i] = 0;
  nc = sys->directlist;
  while(nc != NULL) {
    /* check number and level of neighbors */
    lev = nc->level;
    nn = nc->numnbrs;
    for(i = 0; i < nn; i++) {
      if(lev != ((nc->nbrs)[i])->level) {
	fprintf(stderr, "chkList: level %d cube has a level %d nbr\n", lev,
		((nc->nbrs)[i])->level);
	fprintf(stderr, " ok cubes ");
	for(j = 0; j <= depth; j++) fprintf(stderr, "lev%d: %d ", j, cnt[j]);
	fprintf(stderr, "\n");
	exit(1);
      }
    }
    /* check number of kids */
    if(lev == depth && nc->numkids != 0) {
      fprintf(stderr, "chkList: level %d cube has children\n", lev);
      fprintf(stderr, " ok cubes ");
      for(j = 0; j <= depth; j++) fprintf(stderr, "lev%d: %d ", j, cnt[j]);
      fprintf(stderr, "\n");
      exit(1);
    }
    /* if lowest level, check status of eval and direct vects */
    if(lev == depth) {
      if(nc->dindex == 0 || nc->directnumeles == NULL) {
	fprintf(stderr, "chkList: level %d cube has bad direct info\n", lev);
	fprintf(stderr, " ok cubes ");
	for(j = 0; j <= depth; j++) fprintf(stderr, "lev%d: %d ", j, cnt[j]);