muldisplay.c

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

	fprintf(stderr, "\n");
	exit(1);
      }
      if(nc->evalnumvects == 0 && listtype == EVAL) {
	fprintf(stderr, "chkList: level %d cube has no eval info\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);
      }
    }
    cnt[lev]++;
    if(listtype == DIRECT) nc = nc->dnext;
    else if(listtype == LOCAL) nc = nc->lnext;
    else if(listtype == EVAL) nc = nc->enext;
    else {
      fprintf(stderr, "chkList: bad flag\n");
      exit(1);
    }
  }
  if(listtype == DIRECT) fprintf(stdout, "\nDirect ");
  else if(listtype == LOCAL) fprintf(stdout, "\nLocal ");
  else fprintf(stdout, "\nEval ");
  fprintf(stdout, "list ok: ");
  for(j = 0; j <= depth; j++) fprintf(stdout, "lev%d: %d ", j, cnt[j]);
  fprintf(stdout, "\n\n");
}

/*
  chks a cube for bad cube struct (direct, local or eval) entries - debug only
*/
void chkCube(sys, nc, listtype)
ssystem *sys;
cube *nc;
int listtype;			/* DIRECT, LOCAL or EVAL */
{
  int depth = sys->depth;
  int lev, nn;
  int i, j, k;
  if(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(stdout, "chkCube: level %d cube has a level %d nbr\n", lev,
		((nc->nbrs)[i])->level);
/*	exit(1);*/
      }
    }
    /* check number of kids */
    if(lev == depth && nc->numkids != 0) {
      fprintf(stdout, "chkCube: level %d cube has children\n", lev);
/*      exit(1);*/
    }
    /* if lowest level, check status of eval and direct vects */
    if(lev == depth) {
      if(nc->dindex == 0) {
	fprintf(stdout, "chkCube: level %d cube has zero direct index\n", lev);
/*	exit(1);*/
      }
      if(nc->directnumeles == NULL) {
	fprintf(stdout, 
		"chkCube: level %d cube has NULL directnumeles\n", lev);
/*	exit(1);*/
      }
      if(nc->evalnumvects == 0 && listtype == EVAL) {
	fprintf(stdout, "chkCube: level %d cube has no eval info\n", lev);
/*	exit(1);*/
      }
      if(nc->eval == NULL && listtype == EVAL) {
	fprintf(stdout, "chkCube: level %d cube has no eval pntr\n", lev);
      }
    }
  }
}

/*
  checks the lowest level cubes for trouble using chkCube - debug only
*/
void chkLowLev(sys, listtype)
ssystem *sys;
int listtype;			/* DIRECT, LOCAL or EVAL */
{
  int i, j, k, l, side, depth = sys->depth, cnt = 0;
  cube *nc, *****cubes = sys->cubes;
  for(i = 1, side = 1; i <= depth; i++, side *= 2);
  for(j=0; j < side; j++) {	/* loop through all cubes at level depth */
    for(k=0; k < side; k++) {
      for(l=0; l < side; l++) {
	nc = cubes[depth][j][k][l];
	if(nc != NULL) {
	  chkCube(sys, nc, listtype);
	  cnt++;
	}
      }
    }
  }
  fprintf(stdout,"Total lowest level (level %d) cubes checked = %d\n", 
	  depth, cnt);
}

/*
  dump the contents of a face struct
*/
void dump_face(fp, fac)
face *fac;
FILE *fp;
{
  int i, j;
  face **behind = fac->behind;

  fprintf(fp, "Face %d, %d sides, depth %d, mark %d, greylev %d\n", 
	  fac->index, fac->numsides, fac->depth, fac->mark, fac->greylev);
  fprintf(fp, "  plane: n = (%g %g %g) rhs = %g\n",
	  fac->normal[0], fac->normal[1], fac->normal[2], fac->rhs);
  fprintf(fp, "  behind [depth(index)]:");
  for(i = 0; i < fac->numbehind; i++) {
    fprintf(fp, " %d(%d)", behind[i]->depth, behind[i]->index);
    if(i % 10 == 0 && i != 0) fprintf(fp, "\n");
  }
  i--;
  if(!(i % 10 && i != 0)) fprintf(fp, "\n");
  fprintf(fp, " Corners\n");
  dumpCorners(fp, fac->c, fac->numsides, 3);
}  

/*
  core display routine used below
*/
void dumpSynCore1(str, depth, fcnt, excnt, emcnt, tcnt)
int depth, *fcnt, *excnt, *emcnt, *tcnt;
char *str;
{
  int i;
  fprintf(stdout, "%-13s", str);
  for(i = 0; i <= depth; i++) {
    sprintf(str, "%d/%d/%d/%d ", fcnt[i], excnt[i], emcnt[i], tcnt[i]);
    if(i < 2) fprintf(stdout, "%8s", str);
    else if(i == 2) fprintf(stdout, "%12s", str);
    else if(i == 3) fprintf(stdout, "%16s", str);
    else if(i == 4) fprintf(stdout, "%20s", str);
    else if(i == 5) fprintf(stdout, "%24s", str);
    else fprintf(stdout, "%28s", str);
  }
  fprintf(stdout, "\n");
}
/*
  core display rtn used below
*/
dumpSynCore2(str, depth, cnt)
int depth, *cnt;
char *str;
{
  int i;

  fprintf(stdout, "%-13s", str);
  for(i = 0; i <= depth; i++) {
    sprintf(str, "%d    ", cnt[i]);
    if(i < 2) fprintf(stdout, "%8s", str);
    else if(i == 2) fprintf(stdout, "%12s", str);
    else if(i == 3) fprintf(stdout, "%16s", str);
    else if(i == 4) fprintf(stdout, "%20s", str);
    else if(i == 5) fprintf(stdout, "%24s", str);
    else fprintf(stdout, "%28s", str);
  }
  fprintf(stdout, "\n");
}

/*
  displays number of exact, full, empty and total cubes per level in
  all cubes, and eval, direct, multi and local lists
*/
void dumpSynop(sys)
ssystem *sys;
{
  int i, j, k, l, side, depth = sys->depth, lev;
  int excnt[BUFSIZ], fcnt[BUFSIZ], emcnt[BUFSIZ], tcnt[BUFSIZ];
  extern int *multicnt, *localcnt;
  char str[BUFSIZ];
  cube *****cubes = sys->cubes, *nc;

  for(i = 0; i <= depth; i++) excnt[i] = fcnt[i] = emcnt[i] = tcnt[i] = 0;

  fprintf(stdout, 
	  "\nCUBE AND EXPANSION SYNOPSIS (full/mul_exact/empty/ttl):\n");
  fprintf(stdout, "             ");
  for(i = 0; i <= depth; i++) {
    sprintf(str, "level%d ", i);
    if(i < 2) fprintf(stdout, "%8s", str);
    else if(i == 2) fprintf(stdout, "%12s", str);
    else if(i == 3) fprintf(stdout, "%16s", str);
    else if(i == 4) fprintf(stdout, "%20s", str);
    else if(i == 5) fprintf(stdout, "%24s", str);
    else fprintf(stdout, "%28s", str);
  }
  fprintf(stdout, "\n");
  /* dump cube usage by level */
  for(i = 0, side = 1; i <= depth; i++, side *= 2) {
    for(j=0; j < side; j++) {	/* loop through all cubes at levels >= 0 */
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nc = cubes[i][j][k][l];
	  tcnt[i]++;
	  if(nc != NULL) {
	    lev = nc->level;
	    fcnt[i]++;
	    if(nc->mul_exact == TRUE) excnt[i]++;
	  }
	  else emcnt[i]++;
	}
      }
    }
  }
  sprintf(str, "All cubes");
  dumpSynCore1(str, depth, fcnt, excnt, emcnt, tcnt);
  
  for(i = 0; i <= depth; i++) excnt[i] = fcnt[i] = emcnt[i] = tcnt[i] = 0;
  /* dump cube direct list by level */
  for(nc = sys->directlist; nc != NULL; nc = nc->dnext) {
    lev = nc->level;
    tcnt[lev]++;
    if(nc->upnumvects > 0) fcnt[lev]++;
    else emcnt[lev]++;
    if(nc->mul_exact == TRUE) excnt[lev]++;
  }
  sprintf(str, "Direct list");
  dumpSynCore1(str, depth, fcnt, excnt, emcnt, tcnt);

  for(i = 0; i <= depth; i++) excnt[i] = fcnt[i] = emcnt[i] = tcnt[i] = 0;
  /* dump cube local list by level */
  for(i = 0; i <= depth; i++) {
    for(nc = sys->locallist[i]; nc != NULL; nc = nc->lnext) {
      lev = nc->level;
      tcnt[lev]++;
      if(nc->upnumvects > 0) fcnt[lev]++;
      else emcnt[lev]++;
      if(nc->mul_exact == TRUE) excnt[lev]++;
    }
  }
  sprintf(str, "Local list");
  dumpSynCore1(str, depth, fcnt, excnt, emcnt, tcnt);
    
  for(i = 0; i <= depth; i++) excnt[i] = fcnt[i] = emcnt[i] = tcnt[i] = 0;
  /* dump cube multipole list by level */
  for(i = 0; i <= depth; i++) {
    for(nc = sys->multilist[i]; nc != NULL; nc = nc->mnext) {
      lev = nc->level;
      tcnt[lev]++;
      if(nc->upnumvects > 0) fcnt[lev]++;
      else emcnt[lev]++;
      if(nc->mul_exact == TRUE) excnt[lev]++;
    }
  }
  sprintf(str, "Multi list");
  dumpSynCore1(str, depth, fcnt, excnt, emcnt, tcnt);


  sprintf(str, "Multis built");
  dumpSynCore2(str, depth, multicnt);

  sprintf(str, "Locals built");
  dumpSynCore2(str, depth, localcnt);

}

/*
  dumps the Gaussian unit (statcoulombs/meter^2) charge densities on panels
*/
void dumpChgDen(fp, q, chglist)
double *q;
charge *chglist;
FILE *fp;
{
  charge *panel;

  for(panel = chglist; panel != NULL; panel = panel->next) {
    if(panel->dummy) continue;
    fprintf(fp, "%d\tq/A = %.5e %g", panel->index,
	    q[panel->index]/panel->area, panel->area);
    if(panel->surf->type == CONDTR) fprintf(fp, " CONDTR");
    if(panel->surf->type == DIELEC) fprintf(fp, " DIELEC");
    if(panel->surf->type == BOTH) fprintf(fp, " BOTH");
    fprintf(fp, " (%.3g %.3g %.3g)", panel->x, panel->y, panel->z);
    fprintf(fp, " cond = %d\n", panel->cond);
  }
  fflush(fp);
}

/*
  like dumpMat but different formating and row labels (for dumpMatBldCnts)
*/
void dumpMatCnts(mat, depth, type)
int **mat, depth;
char *type;
{
  int i, j;
  char str[BUFSIZ];

  fprintf(stdout,
	  "\n%s MATRIX BUILD TOTALS (row = from cube, col = to cube):\n", 
	  type);

  for(i = 0; i <= depth; i++) {
    sprintf(str, " to %d ", i);
    if(i == 0) fprintf(stdout, "%13s", str);
    else if(i < 10) fprintf(stdout, "%6s", str);
    else fprintf(stdout, "%7s", str);
  }
  fprintf(stdout, "\n");

  for(i = 0; i <= depth; i++) {
    sprintf(str, "from %d ", i);
    fprintf(stdout, "%-7s", str); /* print row label */
    for(j = 0; j <= depth; j++) {
      sprintf(str, "%d ", mat[i][j]);
      if(j < 10) fprintf(stdout, "%6s", str);
      else fprintf(stdout, "%7s", str);
    }
    fprintf(stdout, "\n");
  }

}

/*
  display matrix build count totals
*/
void dumpMatBldCnts(sys)
ssystem *sys;
{
  int i;
  char type[BUFSIZ];
  extern int **Q2Mcnt, **Q2Lcnt, **Q2Pcnt, **L2Lcnt;
  extern int **M2Mcnt, **M2Lcnt, **M2Pcnt, **L2Pcnt, **Q2PDcnt;

  sprintf(type, "Q2M");
  dumpMatCnts(Q2Mcnt, sys->depth, type);

  sprintf(type, "Q2L");
  dumpMatCnts(Q2Lcnt, sys->depth, type);

  sprintf(type, "Q2P");
  dumpMatCnts(Q2Pcnt, sys->depth, type);

  sprintf(type, "M2M");
  dumpMatCnts(M2Mcnt, sys->depth, type);

  sprintf(type, "M2L");
  dumpMatCnts(M2Lcnt, sys->depth, type);

  sprintf(type, "M2P");
  dumpMatCnts(M2Pcnt, sys->depth, type);

  sprintf(type, "L2L");
  dumpMatCnts(L2Lcnt, sys->depth, type);

  sprintf(type, "L2P");
  dumpMatCnts(L2Pcnt, sys->depth, type);

  sprintf(type, "Q2PDiag");
  dumpMatCnts(Q2PDcnt, sys->depth, type);

}

/* 
  dumps state of important compile flags
*/
void dumpConfig(fp, name)
char *name;
FILE *fp;
{
  int size = -1;		/* for '#define MAXITER size' case */

  fprintf(fp, "\n%s CONFIGURATION FLAGS:\n", name);

  fprintf(fp, " DISCRETIZATION CONFIGURATION\n");

  fprintf(fp, "   WRMETH");
  if(WRMETH == COLLOC)
      fprintf(fp, " == COLLOC (point collocation)\n");
  else if(WRMETH == SUBDOM)
      fprintf(fp, " == SUBDOM (not implemented - do collocation)\n");
  else if(WRMETH == GALKIN)
      fprintf(fp, " == GALKIN (not implemented - do collocation)\n");
  fprintf(fp, "   ELTYPE");
  if(ELTYPE == CONST)
      fprintf(fp, " == CONST (constant panel densities)\n");
  else if(ELTYPE == AFFINE)
      fprintf(fp, " == AFFINE (not implemented - use constant)\n");
  else if(ELTYPE == QUADRA)
      fprintf(fp, " == QUADRA (not implemented - use constant)\n");

  fprintf(fp, " MULTIPOLE CONFIGURATION\n");

  fprintf(fp, "   DNTYPE");
  if(DNTYPE == NOLOCL) 
      fprintf(fp, " == NOLOCL (no locals in dwnwd pass)\n");
  else if(DNTYPE == NOSHFT) 
      fprintf(fp, " == NOSHFT (no local2local shift dwnwd pass)\n");
  else if(DNTYPE == GRENGD) 
      fprintf(fp, " == GRENGD (full Greengard dwnwd pass)\n");
  fprintf(fp, "   MULTI");
  if(MULTI == ON) fprintf(fp, " == ON (include multipole part of P*q)\n");
  else fprintf(fp, " == OFF (don't use multipole part of P*q)\n");
  fprintf(fp, "   RADINTER");
  if(RADINTER == ON) 
      fprintf(fp," == ON (allow parent level interaction list entries)\n");
  else 
   fprintf(fp," == OFF (use only cube level interaction list entries)\n");
  fprintf(fp, "   NNBRS == %d (max distance to a nrst neighbor)\n", NNBRS);
  fprintf(fp, "   ADAPT");
  if(ADAPT == ON) 
      fprintf(fp, " == ON (adaptive - no expansions in exact cubes)\n");
  else fprintf(fp, " == OFF (not adaptive - expansions in all cubes)\n");
  fprintf(fp, "   OPCNT");
  if(OPCNT == ON) 
      fprintf(fp, " == ON (count P*q ops - exit after mat build)\n");
  else fprintf(fp, " == OFF (no P*q op count - iterate to convergence)\n");

  fprintf(fp, "   MAXDEP");
  fprintf(fp, 
	  " == %d (assume no more than %d partitioning levels are needed)\n",
	  MAXDEP, MAXDEP);

  fprintf(fp, "   NUMDPT");
  fprintf(fp, 
	  " == %d (do %d potential evaluations for each dielectric panel)\n",
	  NUMDPT, NUMDPT);

  fprintf(fp, " LINEAR SYSTEM SOLUTION CONFIGURATION\n");

  fprintf(fp, "   ITRTYP");
  if(ITRTYP == GCR)
      fprintf(fp, " == GCR (generalized conjugate residuals)\n");
  else if(ITRTYP == GMRES)
      fprintf(fp, " == GMRES (generalized minimum residuals)\n");
  else fprintf(fp, " == %d (not implemented - use GCR)\n", ITRTYP);

  fprintf(fp, "   PRECOND");
  if(PRECOND == BD) {
    fprintf(fp, 
	    " == BD (use block diagonal preconditioner)\n");
  }
  else if(PRECOND == OL) {
    fprintf(fp, 
	    " == OL (use overlap preconditioner)\n");
  }
  else if(PRECOND == NONE) {
    fprintf(fp, 
	    " == NONE (no preconditioner)\n");
  }
  else fprintf(fp, " == %d (not implemented - use BD, OL or NONE)\n", PRECOND);

  fprintf(fp, "   DIRSOL");
  if(DIRSOL == ON) 
      fprintf(fp, " == ON (do the whole calculation directly)\n");
  else fprintf(fp, " == OFF (do the calculation iteratively)\n");

  fprintf(fp, "   EXPGCR");
  if(EXPGCR == ON) 
      fprintf(fp, " == ON (do all P*q's explicitly w/full matrix)\n");
  else fprintf(fp, " == OFF (do P*q's with multipole)\n");

  fprintf(fp, "   MAXITER");
  if(MAXITER < 0) {
    fprintf(fp, " == size (for n panel system, do at most n iterations)\n");
  }
  else fprintf(fp, " == %d (stop after %d iterations if not converged)\n", 
	  MAXITER, MAXITER);

  fprintf(fp, "   EXRTSH");
  fprintf(fp, 
	  " == %g (exact/ttl cubes > %g on lowest level => stop refinement)\n",
	  EXRTSH, EXRTSH);
}


/*
  pads a string on the right up to a given length, truncates if too long
*/
char *padName(tostr, frstr, len)
char *tostr, *frstr;
int len;
{
  int i;

  for(i = 0; frstr[i] != '\0'; i++) tostr[i] = frstr[i];
  if(i > len) tostr[len] = '\0';		/* truncate */
  else {			/* pad */
    for(; i < len; i++) tostr[i] = ' ';
    tostr[len] = '\0';
  }
  return(tostr);
}

/*
  returns a string of spaces (doesn't stdio have this somewhere?)
*/
char *spaces(str, num)
char *str;
int num;
{
  int i;

  for(i = 0; i < num; i++) str[i] = ' ';
  str[num] = '\0';
  return(str);
}
    
/*
  prints the capacitance matrix with symmetrized (averaged) off-diagonals
  - mks units (Farads) are used
  - some attempt to scale (eg pF, nF, uF etc) is made
  - also checks for M-matrix sign pattern, diag dominance
*/
void mksCapDump(capmat, numconds, relperm, name_list)
double **capmat, relperm;
int numconds;
Name **name_list;
{
  int i, j, toobig, toosmall, maxlen, sigfig, colwidth, i_killed, j_killed;
  int first_offd;
  double maxdiag = 0.0, minoffd, rowttl, rowdiag, scale = 1.0, **sym_mat;
  double mat_entry;