muldisplay.c

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

  char unit[BUFSIZ], name[BUFSIZ], *padName(), *spaces(), cond_name[BUFSIZ];
  char *getConductorName();
  extern NAME *start_name;	/* NAME structs giving conductor names */
  Name *cname;
  extern ITER *kill_num_list, *kinp_num_list;
  extern double iter_tol;

  first_offd = TRUE;
  minoffd = capmat[1][1];	/* this entry is always present */
				/* - in the 1 cond case, assign is still ok */

  /* set up symetrized matrix storage */
  CALLOC(sym_mat, numconds+1, double*, ON, AMSC);
  for(i=1; i <= numconds; i++)  {
    CALLOC(sym_mat[i], numconds+1, double, ON, AMSC);
  }

  /* get the smallest and largest (absolute value) symmetrized elements */
  /* check for non-M-matrix symmetrized capacitance matrix */
  for(i = 1; i <= numconds; i++) {

    /* skip conductors removed from input */
    if(want_this_iter(kinp_num_list, i)) continue;

    i_killed = want_this_iter(kill_num_list, i);

    if(capmat[i][i] <= 0.0 && !i_killed) {
      fprintf(stderr, "\nmksCapDump: Warning - capacitance matrix has non-positive diagonal\n  row %d\n", i+1);
    }
    maxdiag = MAX(maxdiag, fabs(capmat[i][i]));
    rowttl = 0.0;
    for(j = 1; j <= numconds; j++) {

      /* skip conductors removed from input */
      if(want_this_iter(kinp_num_list, j)) continue;

      if(j == i) {
	sym_mat[i][i] = capmat[i][i];
	continue;
      }

      /* if this column was not calculated and neither was the column
         with the same number as the current row, then symetrized mat has
	 no entry at [i][j], [j][i] */
      j_killed = want_this_iter(kill_num_list, j);
      if(i_killed && j_killed) continue;

      /* if this column was calculated but column with the same number
         as the current row wasnt, then symmetrized mat has unaveraged entry 
	 at [i][j], [j][i] */
      else if(i_killed && !j_killed) mat_entry = capmat[i][j];

      /* if this column was not calculated but column with the same number
         as the current row was, then symmetrized mat has unaveraged entry 
	 at [i][j], [j][i] */
      else if(!i_killed && j_killed) mat_entry = capmat[j][i];

      /* if this column was calculated and column with the same number
         as the current row was also, then symmetrized mat has averaged entry 
	 at [i][j], [j][i] */
      else mat_entry = (capmat[i][j] + capmat[j][i])/2.0;

      rowttl += mat_entry;
      if(mat_entry >= 0.0) {
	fprintf(stderr, "\nmksCapDump: Warning - capacitance matrix has non-negative off-diagonals\n  row %d col %d\n", i, j);
      }
      if(fabs(mat_entry) != 0.0) {
	if(first_offd) {
	  minoffd = fabs(mat_entry);
	  first_offd = FALSE;
	}
	else minoffd = MIN(minoffd, fabs(mat_entry));
      }

      sym_mat[i][j] = mat_entry;
    }
    if(rowttl + capmat[i][i] <= 0.0 && !i_killed) {
      fprintf(stderr, "\nmksCapDump: Warning - capacitance matrix is not strictly diagonally dominant\n  due to row %d\n", i);
    }
  }

  /* figure the units to use for the matrix entries 
     - set up so smallest element is between 0.1 and 10 */
  if(minoffd*FPIEPS*relperm*scale > 10.0) toobig = TRUE;
  else toobig = FALSE;
  if(minoffd*FPIEPS*relperm*scale < 0.1) toosmall = TRUE;
  else toosmall = FALSE;
  while(toobig == TRUE || toosmall == TRUE) {
    if(toobig == TRUE) {
      scale *= 1e-3;
      if(minoffd*FPIEPS*relperm*scale <= 10.0) break;
    }
    if(toosmall == TRUE) {
      scale *= 1e+3;
      if(minoffd*FPIEPS*relperm*scale >= 0.1) break;
    }
  }

  /* get the appropriate unit string */
  if(scale == 1e-18) strcpy(unit, "exa");
  else if(scale == 1e-15) strcpy(unit, "peta");
  else if(scale == 1e-12) strcpy(unit, "tera");
  else if(scale == 1e-9) strcpy(unit, "giga");
  else if(scale == 1e-6) strcpy(unit, "mega");
  else if(scale == 1e-3) strcpy(unit, "kilo");
  else if(scale == 1.0) strcpy(unit, "");
  else if(scale == 1e+3) strcpy(unit, "milli");
  else if(scale == 1e+6) strcpy(unit, "micro");
  else if(scale == 1e+9) strcpy(unit, "nano");
  else if(scale == 1e+12) strcpy(unit, "pico");
  else if(scale == 1e+15) strcpy(unit, "femto");
  else if(scale == 1e+18) strcpy(unit, "atto");
  else sprintf(unit, "every unit is %g ", 1/scale);

  /* get the length of the longest name */
  maxlen = 0;
  for(i = 1; i <= numconds; i++) {
    maxlen = MAX(strlen(getConductorName(i, name_list)), maxlen);
  }

  /* print the matrix */
  sigfig = 2+log10(1.0/iter_tol);	/* get no. significant figs to prnt */
  colwidth = sigfig+6;		/* field width for cap mat columns */
  if(ITRDAT == OFF) fprintf(stdout, "\n");
  if(kill_num_list != NULL) 
      fprintf(stdout, "\nPARTIAL CAPACITANCE MATRIX, %sfarads\n", unit);
  else fprintf(stdout, "\nCAPACITANCE MATRIX, %sfarads\n", unit);
  if(numconds < 10) fprintf(stdout, "%s", spaces(unit, maxlen+2));
  else if(numconds < 100) fprintf(stdout, "%s", spaces(unit, maxlen+3));
  else fprintf(stdout, "%s", spaces(unit, maxlen+4));
  for(j = 1; j <= numconds; j++) { /* column headings */
    if(want_this_iter(kinp_num_list, j)) continue;
    sprintf(name, "%d ", j);
    sprintf(unit, "%%%ds", colwidth+1);
    fprintf(stdout, unit, name);
  }
  fprintf(stdout, "\n");
  for(i = 1; i <= numconds; i++) { /* rows */

    /* skip conductors removed from input */
    if(want_this_iter(kinp_num_list, i)) continue;

    sprintf(unit, "%d", i);

    strcpy(cond_name, getConductorName(i, name_list));

    if(numconds < 10)
	fprintf(stdout, "%s %1s", padName(name, cond_name, maxlen), unit);
    else if(numconds < 100)
	fprintf(stdout, "%s %2s", padName(name, cond_name, maxlen), unit);
    else
	fprintf(stdout, "%s %3s", padName(name, cond_name, maxlen), unit);

    for(j = 1; j <= numconds; j++) {

      /* skip conductors removed from input */
      if(want_this_iter(kinp_num_list, j)) continue;

      if(want_this_iter(kill_num_list, i) 
	 && want_this_iter(kill_num_list, j)) {
	/* print a blank if capacitance was not calculated */
	fprintf(stdout, "%s", spaces(unit, colwidth+1));
      }
      else {
	sprintf(unit, " %%%d.%dg", colwidth, sigfig);
	fprintf(stdout, unit, scale*FPIEPS*relperm*sym_mat[j][i]);
      }
    }
    fprintf(stdout, "\n");
  }
}

/*
  dumps brief information about multipole set up
*/
void dumpMulSet(sy, numLev, order)
ssystem *sy;
int numLev, order;
{
  int numcubes, numsides, i, multerms();

  for(numcubes = 1, i = 0; i < numLev; numcubes *= 8, i++);
  for(numsides = 1, i = 0; i < numLev; numsides *= 2, i++);

  fprintf(stdout, "\nMULTIPOLE SETUP SUMMARY\n");
  fprintf(stdout, "Level 0 cube extremal coordinates\n");
  fprintf(stdout, "  x: %g to %g\n", 
	  sy->minx, sy->minx + numsides * (sy->length));
  fprintf(stdout, "  y: %g to %g\n", 
	  sy->miny, sy->miny + numsides * (sy->length));
  fprintf(stdout, "  z: %g to %g\n", 
	  sy->minz, sy->minz + numsides * (sy->length));
  fprintf(stdout, "Level %d (lowest level) cubes\n  %d total\n", 
	  numLev, numcubes);
  fprintf(stdout, 
	  "  side length = %g\n  maximum number of panels in each = %d\n",
	  sy->length, sy->mul_maxlq);
  fprintf(stdout, "  maximum number of evaluation points in each = %d\n",
	  sy->loc_maxlq);
  fprintf(stdout, 
	  "Maximum number of panels treated with a multipole expansion = %d\n",
	  sy->max_panel);
  fprintf(stdout, 
  "Maximum number of evaluation points treated with a local expansion = %d\n",
	  sy->max_eval_pnt);
  fprintf(stdout, 
	  "Maximum number of panels treated exactly = %d (limit = %d)\n",
	  sy->mul_maxq, multerms(order));
  fprintf(stdout, 
   "Maximum number of evaluation points treated exactly = %d (limit = %d)\n",
	  sy->loc_maxq, multerms(order));
}

/*
  dump all the preconditioner matrices in the direct list cubes as one
  big matrix for matlab read in "Ctil"
  - figures preconditioner by multiplying it by columns of the inverse
    and dumping results one column at a time
  - figures the unpreconditioned matrix (using calcp) and dumps it to "P"
  - type determines which of the matrices to dump
*/
void dump_preconditioner(sys, chglist, type)
charge *chglist;
ssystem *sys;
int type;			/* 1=>dump P and C; 2=>P only; 3=>C only */
{
  int num_panels, i, j;
  charge *pp, *pi;
  cube *cp;
  double calcp();
  FILE *fp;

  /* find the number of panels */
  for(num_panels = 0, pp = chglist; pp != NULL; pp = pp->next, num_panels++);

  /* open the output file */
  if((fp = fopen("prec.mat","w")) == NULL) {
    fprintf(stderr, "dump_preconditioner: can't open `prec.mat'\n");
    exit(1);
  }
  
  if(type == 1 || type == 3) {
    fprintf(stdout, "\nDumping preconditioner to `prec.mat' as `Ctil'\n");
    /* dump the preconditioner one column at a time */
    /* - savemat arg "type" can be used to make rowwise dumps
         type = x0xx  =>  columnwise dumps 
         type = x1xx  =>  rowwise dumps (see matlab manual) */
    for(j = 1; j < num_panels+1; j++) {
      /* load I col into q and zero p */
      for(i = 0; i < num_panels+1; i++) {
	if(i == j) sys->q[i] = 1.0;
	else sys->q[i] = 0.0;
      }
      /* figure the column of C in p (xfered to q after calculation) */
      mulPrecond(sys, PRECOND);
      /* dump the preconditioner column */
      if(j == 1) savemat_mod(fp, 1000, "Ctil", num_panels, num_panels, 0, 
			     &(sys->q[1]), (double *)NULL, 0, num_panels);
      else savemat_mod(fp, 1000, "Ctil", num_panels, num_panels, 0, 
		       &(sys->q[1]), (double *)NULL, 1, num_panels);

    }
  }

  if(type == 1 || type == 2) {
    fprintf(stdout, "\nDumping pot. coeff. mat. to `prec.mat' as `P'\n");
    /* dump the P matrix - DANGER: does n^2 calcp() calls,
       but storage is only n */
    /* are q indices from 1?? */
    for(j = 1; j < num_panels+1; j++) {
      /* find the panel with the current index (select column) */
      for(pp = chglist; pp != NULL; pp = pp->next) {
	if(pp->index == j) break;
      }
      if(pp == NULL) {
	fprintf(stderr, "dump_preconditioner: no charge with index %d\n", j);
	exit(1);
      }
      for(i = 0; i < num_panels+1; i++) sys->p[i] = 0.0;
      /* find the column---influence of q_j on potentials at each q_i */
      for(i = 1, pi = chglist; pi != NULL; i++, pi = pi->next) {
	sys->p[pi->index] = calcp(pp, pi->x, pi->y, pi->z, NULL);
      }
      /* dump the column */
      if(j == 1) savemat_mod(fp, 1000, "P", num_panels, num_panels, 0, 
			     &(sys->p[1]), (double *)NULL, 0, num_panels);
      else savemat_mod(fp, 1000, "P", num_panels, num_panels, 0, 
		       &(sys->p[1]), (double *)NULL, 1, num_panels);

    }
  }

  fclose(fp);
}

/*
  do an n^2/2 check to see if any panels have the same center points
  (debug only)
*/
int has_duplicate_panels(fp, chglst)
charge *chglst;
FILE *fp;
{
  int no_duplicates;
  charge *cp, *cpinner;

  no_duplicates = TRUE;
  for(cp = chglst; cp != NULL; cp = cp->next) {
    for(cpinner = cp->next; cpinner != NULL; cpinner = cpinner->next) {
      if(cp->x == cpinner->x && cp->y == cpinner->y && cp->z == cpinner->z) {
	no_duplicates = FALSE;

	if(cp->surf->type == CONDTR) fprintf(fp, "Panels %d(CONDTR)",
						cp->index);
	if(cp->surf->type == DIELEC) fprintf(fp, "Panels %d(DIELEC)",
						cp->index);
	if(cp->surf->type == BOTH) fprintf(fp, "Panels %d(BOTH)",
					      cp->index);

	if(cpinner->surf->type == CONDTR) fprintf(fp, " and %d(CONDTR)",
						cpinner->index);
	if(cpinner->surf->type == DIELEC) fprintf(fp, " and %d(DIELEC)",
						cpinner->index);
	if(cpinner->surf->type == BOTH) fprintf(fp, " and %d(BOTH)",
					      cpinner->index);

	fprintf(fp, " both have center (%.3g %.3g %.3g)\n",
		cp->x, cp->y, cp->z);
      }
    }
  }

  if(no_duplicates) return(FALSE);
  else return(TRUE);
}

#if DSQ2PD == ON
/*
  dump the condensed matrix for matlab use
*/
void dumpQ2PDiag(nextc)
cube *nextc;
{
  int i, j, ind, pos_d, neg_d;
  double temp[BUFSIZ], temp_mat[100][100], **rmat;
  double pos_fact, neg_fact, panel_fact;
  FILE *fp, *fopen();

  /* dump the diag matrix to a matlab file along with its dummy vector */
  /*   must complie on sobolev with /usr/local/matlab/loadsave/savemat.o */
  if((fp = fopen("Q2PDiag.mat", "w")) == NULL) {
    fprintf(stderr, "dumpQ2PDiag: can't open `Q2PDiag.mat' to write\n");
    exit(1);
  }
  if(sizeof(temp) < nextc->upnumeles[0]*nextc->upnumeles[0]) {
    fprintf(stderr, "dumpQ2PDiag: temporary arrays not big enough\n");
    exit(1);
  }

  /* incorporate the electric field evaluation into the matrix */
  /* if a row corresponds to a flux density difference, alter it */
  rmat = nextc->directmats[0];
  for(i = 0; i < nextc->upnumeles[0]; i++) {
    if(nextc->chgs[i]->dummy) {
      for(j = 0; j < nextc->upnumeles[0]; j++) temp_mat[i][j] = 0.0;
      continue;
    }

    if(nextc->chgs[i]->surf->type == CONDTR) {
      for(j = 0; j < nextc->upnumeles[0]; j++) {
	temp_mat[i][j] = rmat[i][j];
      }
    }
    else {

      pos_fact 
	  = nextc->chgs[i]->surf->outer_perm/nextc->chgs[i]->pos_dummy->area;
      pos_d = nextc->chgs[i]->pos_dummy->index - 1;
      neg_fact 
	  = nextc->chgs[i]->surf->inner_perm/nextc->chgs[i]->neg_dummy->area;
      neg_d = nextc->chgs[i]->neg_dummy->index - 1;
      panel_fact = pos_fact + neg_fact;

      for(j = 0; j < nextc->upnumeles[0]; j++) {
	temp_mat[i][j] = pos_fact*rmat[pos_d][j] - panel_fact*rmat[i][j]
	    + neg_fact*rmat[neg_d][j];
      }
    }
  }
  
  /* flatten the Q2PDiag matrix */
  for(j = ind = 0; j < nextc->upnumeles[0]; j++) {
    for(i = 0; i < nextc->upnumeles[0]; i++) {
      temp[ind++] = temp_mat[i][j];
    }
  }
  savemat(fp, 1000, "A", nextc->upnumeles[0], nextc->upnumeles[0],
	  0, temp, (double *)NULL);

  /* make the is_dummy vector a vector of doubles */
  for(i = 0; i < nextc->upnumeles[0]; i++) 
      temp[i] = (double)(nextc->nbr_is_dummy[0][i]);
  savemat(fp, 1000, "is_dummy", nextc->upnumeles[0], 1,
	  0, temp, (double *)NULL);

  /* make a vector with 0 => CONDTR 1 => DIELEC 2 => BOTH -1 => dummy */
  for(i = 0; i < nextc->upnumeles[0]; i++) {
    if(nextc->chgs[i]->dummy) temp[i] = -1.0;
    else temp[i] = (double)(nextc->chgs[i]->surf->type);
  }
  savemat(fp, 1000, "surf_type", nextc->upnumeles[0], 1,
	  0, temp, (double *)NULL);

  fclose(fp);
  fprintf(stdout, "Dumped Q2PDiag matrix to `Q2PDiag.mat'\n");
}

#endif

#if 1 == 0
/*
  for debug only on SPARC2 -- NaN trap error handler (see man matherr)
*/
int matherr(exc)
struct exception *exc;
{
  fprintf(stderr, "matherr: ");

  if(exc->type == DOMAIN) fprintf(stderr, "DOMAIN error in");
  else if(exc->type == SING) fprintf(stderr, "SING error in");
  else if(exc->type == OVERFLOW) fprintf(stderr, "OVERFLOW error in");
  else if(exc->type == UNDERFLOW) fprintf(stderr, "UNDERFLOW error in");

  fprintf(stderr, " %s\n", exc->name);

  fprintf(stderr, " args: %g %g returning: %g\n", exc->arg1, exc->arg2,
	  exc->retval);
  exit(1);
}
#endif

/*
  debug only - check a vector to make sure it has zeros in dummy entries
*/
int chkDummy(vector, is_dummy, size)
double *vector;
int *is_dummy;
int size;
{
  int i, first = TRUE;

  for(i = 0; i < size; i++) {
    if(is_dummy[i] && vector[i] != 0.0) {
      if(first) {
	first = FALSE;
	fprintf(stderr, "\nchkDummy: entries should be 0.0: %d", i);
      }
      else fprintf(stderr, " %d", i);
    }
  }
  if(!first) fprintf(stderr, "\n");
}

/*
  debug only - print message if dummy list doesn't match panel list
*/
void chkDummyList(panels, is_dummy, n_chgs)
charge **panels;
int *is_dummy;
int n_chgs;
{
  int i;
  int first = TRUE;
  
  for(i = 0; i < n_chgs; i++) {
    if(is_dummy[i] && !panels[i]->dummy || !is_dummy[i] && panels[i]->dummy) {
      if(first) {
	first = FALSE;
	fprintf(stderr, "chkDummyList: inconsistent dummy list entries:\n");
      }
      fprintf(stderr, " %d is_dummy = %d, panel->dummy = %d\n", i,
	      is_dummy[i], panels[i]->dummy);
    }
  }

}

/*
  print the conductor names to a file
*/
void dumpCondNames(fp, name_list)
FILE *fp;
Name *name_list;
{ 
  int i;
  char *last_alias();
  Name *cur_name;

  fprintf(fp, "CONDUCTOR NAMES\n");
  for(cur_name = name_list, i = 0; cur_name != NULL; 
      cur_name = cur_name->next, i++) {
    fprintf(fp, "  %d `%s'\n", i+1, last_alias(cur_name));
  }
}

/*
  debug only: dump state of internal conductor name list
*/
int dumpNameList(name_list)
Name *name_list;
{
  Name *cur_name, *cur_alias;

  for(cur_name = name_list; cur_name != NULL; cur_name = cur_name->next) {
    fprintf(stdout, "`%s'\n", cur_name->name);
    for(cur_alias = cur_name->alias_list; cur_alias != NULL; 
	cur_alias = cur_alias->next) {
      fprintf(stdout, "  `%s'\n", cur_alias->name);
    }
  }
  return(TRUE);
}