mulsetup.c

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

/*
Find all the nearest neighbors.
At the bottom level, get neighbors due to a parents being exact.
*/
static getnbrs(sys)
ssystem *sys;
{
cube *nc, *np, *****cubes = sys->cubes;
int depth = sys->depth;
int i, j, k, l, m, n, p, side, es;
int numnbrs;

/* Return if depth = 0, no neighbors. */
  if(depth == 0) return;

/*
At the every level, get the nearest nbrs combined with nbrs due to parents
being exact.
*/
  /* exactness for local expansion is checked - nbrs used only in dwnwd pass */
  for(i = 1, side = 2; i <= depth; i++, side *= 2) {
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nc = cubes[i][j][k][l];
	  if(nc != NULL) {
	    /* Find sidelength of exact cube. */
	    for(es=1, np=nc->parent; np->loc_exact==TRUE; 
		np = np->parent, es *= 2); /* exact -> loc_exact 1Apr91 */
	  
	    /* Stack up the nearest nbrs plus nbrs in exact cube. */
	    numnbrs = 0;
	    for(m = MIN((j-NNBRS), es * (j/es));
		m < MAX((j+NNBRS+1), es * (1 + (j / es))); m++) {
	      for(n = MIN((k-NNBRS), es * (k/es));
		  n < MAX((k+NNBRS+1), es * (1 + (k/es))); n++) {
		for(p = MIN((l-NNBRS), es * (l/es));
		    p < MAX((l+NNBRS+1), es * (1+(l/es))); p++) {
		  if( (m >= 0) && (n >= 0) && (p >= 0)
		     && (m < side) && (n < side) && (p < side)
		     && ((m != j) || (n != k) || (p != l))
		     && (cubes[i][m][n][p] != NULL)) {
		    cstack[numnbrs++] = cubes[i][m][n][p];
		  }
		}
	      }
	    }
	    nc->numnbrs = numnbrs;
	    if(nc->numnbrs > 0) CALLOC(nc->nbrs, numnbrs, cube*, ON, AMSC);
	    for(m=numnbrs-1; m >= 0; m--) nc->nbrs[m] = cstack[m];
	  }
	}
      }
    }
  }
}

/*
  returns the number of charges in the lowest level cubes contained in "cp"
*/
int cntDwnwdChg(cp, depth)
int depth;			/* number of lowest level */
cube *cp;
{
  int i;
  int cnt;
  cube *kidc;

  if(cp->level == depth) return(cp->upnumeles[0]);
  else for(i = 0; i < cp->numkids; i++) 
      cnt += cntDwnwdChg(cp->kids[i], depth);
  return(cnt);
}

/* 
Set up the links between cubes requiring multi work on each level, one
for the cubes requiring local expansion work, one for the cubes requiring
direct methods and one for cubes with potential evaluation points. 
Note, upnumvects and exact must be set!!!
*/
static linkcubes(sys)
ssystem *sys;
{
  cube *nc, **plnc, **pdnc, **pmnc, *****cubes = sys->cubes;
  int i, j, k, l, cnt = 0;
  int dindex, side, depth=sys->depth, numterms=multerms(sys->order);

  /* Allocate the vector of heads of cubelists. */
  CALLOC(sys->multilist, sys->depth+1, cube*, ON, AMSC);
  CALLOC(sys->locallist, sys->depth+1, cube*, ON, AMSC);

  pdnc = &(sys->directlist);
  for(dindex = 1, i=0, side = 1; i <= sys->depth; i++, side *= 2) {
    pmnc = &(sys->multilist[i]);
    plnc = &(sys->locallist[i]);
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nc = cubes[i][j][k][l];
	  if(nc != NULL) {
	    /* Do the multi expansion if the cube is not treated exactly. */
	    if(i > 1) {		/* no multis over the root cube and its kids */
	      if(nc->mul_exact == FALSE) { /* exact -> mul_exact 1Apr91 */
		CALLOC(nc->multi, numterms, double, ON, AMSC);
		*pmnc = nc;
		pmnc = &(nc->mnext);
	      }
	    }
	    
	    /* Do the local expansion on a cube if it has chgs inside, and it's
	     not exact and not the root (lev 0) nor one of its kids (lev 1). */
	    if(i > 1) {		/* no locals with level 0 or 1 */
	      if(nc->loc_exact == FALSE) { /* exact -> loc_exact 1Apr91 */
		*plnc = nc;
		plnc = &(nc->lnext);
		CALLOC(nc->local, numterms, double, ON, AMSC);
	      }
	    }

	    /* Add to direct list if at bot level and not empty. */
	    if(i == depth) { 
	      *pdnc = nc;  /* For the direct piece, note an index. */
	      pdnc = &(nc->dnext);
	      nc->dindex = dindex++;
	    }
	  }
	}
      }
    }
  }
}

/*
Determine maximum number of chgs contained in a single cube.
*/
static setMaxq(sys)
ssystem *sys;
{
  int i, j, k, l, side, p, kids_are_exact, all_null, depth = sys->depth;
  int mul_maxq, mul_maxlq, loc_maxq, loc_maxlq, num_chgs, real_panel_cnt;
  cube *nc, *****cubes = sys->cubes;

  mul_maxq = mul_maxlq = loc_maxq = loc_maxlq = 0;
  for(i = 1, side = 2; i <= depth; i++, side *= 2) {
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nc = cubes[i][j][k][l];
	  if(nc != NULL) {
	    if(nc->mul_exact == TRUE) {
	      num_chgs = 0;
	      for(p = 0; p < nc->upnumeles[0]; p++) {
		if(!nc->nbr_is_dummy[0][p]) num_chgs++;
	      }
	      mul_maxq = MAX(mul_maxq, num_chgs);
	      if(i == depth) mul_maxlq = MAX(mul_maxlq, num_chgs); 
	    }
	    if(nc->loc_exact == TRUE) {
	      loc_maxq = MAX(loc_maxq, nc->upnumeles[0]);
	      if(i == depth) loc_maxlq = MAX(loc_maxlq,nc->upnumeles[0]); 
	    }
	  }
	}
      }
    }
  }
  sys->loc_maxq = loc_maxq;	/* max evaluation points, over all loc_exact */
  sys->loc_maxlq = loc_maxlq;	/* max evaluation pnts, over lowest level */
  sys->mul_maxq = mul_maxq;	/* max panels, over all mul_exact cubes */
  sys->mul_maxlq = mul_maxlq;	/* max panels, over lowest level cubes */

  /* find the maximum #panels in all non-exact cubes w/exact (or no) kids */
  sys->max_panel = 0;
  for(j = 2; j <= depth; j++) {
    for(nc = sys->multilist[j]; nc != NULL; nc = nc->mnext) {
      if(nc->level == depth) {
	real_panel_cnt = 0;
	for(i = 0; i < nc->upnumeles[0]; i++) {
	  if(!nc->nbr_is_dummy[0][i]) real_panel_cnt++;
	}
	sys->max_panel = MAX(sys->max_panel, real_panel_cnt);
      }
      else {
	kids_are_exact = all_null = TRUE;
	real_panel_cnt = 0;
	for(i = 0; i < nc->numkids && kids_are_exact; i++) {
	  if(nc->kids[i] == NULL) continue;
	  all_null = FALSE;
	  if(!nc->kids[i]->mul_exact) kids_are_exact = FALSE;
	  else {		/* count real panels */
	    for(l = 0; l < nc->kids[i]->upnumeles[0]; l++) {
	      if(!((nc->kids[i]->nbr_is_dummy[0])[l])) real_panel_cnt++;
	    }
	  }
	}
	if(kids_are_exact && !all_null) {
	  sys->max_panel = MAX(sys->max_panel, real_panel_cnt);
	}
      }
    }
  }

  /* find the maximum #eval points in all non-exact cubes w/exact children */
  sys->max_eval_pnt = 0;
  for(j = 2; j <= depth; j++) {
    for(nc = sys->locallist[j]; nc != NULL; nc = nc->lnext) {
      if(nc->level == depth) {
	sys->max_eval_pnt = MAX(sys->max_eval_pnt, nc->upnumeles[0]);
      }
      else {
	kids_are_exact = all_null = TRUE;
	real_panel_cnt = 0;
	for(i = 0; i < nc->numkids && kids_are_exact; i++) {
	  if(nc->kids[i] == NULL) continue;
	  all_null = FALSE;
	  if(!nc->kids[i]->loc_exact) kids_are_exact = FALSE;
	  else real_panel_cnt += nc->kids[i]->upnumeles[0];
	}
      }
      if(kids_are_exact && !all_null)
	  sys->max_eval_pnt = MAX(sys->max_eval_pnt, real_panel_cnt);
    }
  }
}

/* 
  markup sets the flag to "flag" in the child and its nearest nbrs
*/
static markUp(child, flag)
cube *child;
int flag;
{
  int i,j;
  cube *nc, *np;

  child->flag = flag;
  for(i = 0; i < child->numnbrs; i++) {
    child->nbrs[i]->flag = flag;
  }
}

/* 
  forms the true interaction list (see also comment at mulMatEval())
   for cube "child", excluding only empty cubes
  -interaction list pointer is saved in the interList cube struct field
*/
static getInter(child)
cube *child;
{
  int i, j, vects, usekids, lc, jc, kc, ln, jn, kn;
  int numnbr = (child->parent)->numnbrs; /* number of neighbors */
  cube **nbrc = (child->parent)->nbrs; /* list of neighbor pointers */
  cube *sib;			/* pointer to sibling (same level as child) */
  cube **pstack = &(cstack[0]); /* temporary storage pointer */

  /* mark the child cube and all its neighbors */
  markUp(child, TRUE);

  /* unmarked children of child's parent's neighbors become the ilist */
  for(i = 0; i < numnbr; i++) { /* loop on neighbors */
    /* Check nbr's kids for a marked kid. */
    for(usekids = FALSE, j = 0; j < nbrc[i]->numkids; j++) { 
      sib = (nbrc[i]->kids)[j];
      if((sib != NULL) && (sib->flag == TRUE)) { usekids = TRUE; break; };
    }
    /* Use nbr if no kids marked. */
    /* ...and it's really not a 1st nrst nbr of the parent 
       - this stops parent-sized cubes from getting into the ilist
         when they have empty child-sized cubes that are 2nd or 1st
	 nrst nbrs of the child cube 
       - should work with NNBRS = 1 (never allows parent-sized in list)
         and NNBRS > 2 (but cannot allow greater than parent-sized)
       (29May90) */
#if ON == ON
    lc = (child->parent)->l;
    jc = (child->parent)->j;
    kc = (child->parent)->k;
    ln = nbrc[i]->l;
    jn = nbrc[i]->j;
    kn = nbrc[i]->k;
    if((RADINTER == ON) && (usekids == FALSE) &&
       ((lc-1 != ln && lc+1 != ln && lc != ln)
       || (jc-1 != jn && jc+1 != jn && jc != jn)
       || (kc-1 != kn && kc+1 != kn && kc != kn))) {  
      *pstack = nbrc[i];
      pstack++;
    }
#else				/* USE THIS PART FOR TESTING ONLY */
    if(RADINTER && (usekids == FALSE)) { /* PRODUCES INCORRECT ILISTS!!! */
      *pstack = nbrc[i];
      pstack++;
    }
#endif
    else for(j = 0; j < nbrc[i]->numkids; j++) { /* use nbr's kids. */
      sib = (nbrc[i]->kids)[j];	/* get sib of child cube of interest */
      if((sib != NULL) && (sib->flag == FALSE)) { 
	*pstack = sib;
	pstack++;
      }
    }
  }

  /* clear all the flags */
  markUp(child, FALSE);

  /* allocate and save the interaction list */
  child->interSize = vects = pstack - &(cstack[0]);
  if(vects > 0) CALLOC(child->interList, vects, cube*, ON, AMSC);
  for(j = 0; j < vects; j++) child->interList[j] = cstack[j];

  return(vects);		/* return number of interaction elements */
}

/*
  generates explicit, true interaction lists for all non-empty cubes w/lev > 1
*/
static getAllInter(sys)
ssystem *sys;
{
  int i, j, k, l, side, depth = sys->depth;
  cube *nc, *****cubes = sys->cubes;
  for(i = 2, side = 4; i <= depth; i++, side *= 2) {
    for(j=0; j < side; j++) {	/* loop through all cubes at levels > 1 */
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nc = cubes[i][j][k][l];
	  if(nc != NULL) getInter(nc);
	}
      }
    }
  }
}

/* 
  inits the dummy and dielec mask vectors; used to tell which entries to skip
  - mask vectors are redundant (could read flags in charge struct) 
  - done for speed in potential eval loop
*/
static set_vector_masks(sys)
ssystem *sys;
{
  int i;
  cube *cp;

  for(cp = sys->directlist; cp != NULL; cp = cp->dnext) {
    for(i = 0; i < cp->upnumeles[0]; i++) {
      if(!(cp->nbr_is_dummy[0][i] = cp->chgs[i]->dummy))
	  cp->is_dielec[i] = 
	      (cp->chgs[i]->surf->type == DIELEC 
	       || cp->chgs[i]->surf->type == BOTH);
    }
  }

}