mulsetup.c

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

    /* check if current charge is same as those already in the cube `nextc' */
    for(compflag = FALSE, i = (nextc->upnumeles[0] - 1); i >= 0; i--) {
      compq = nextc->chgs[i];
      if((compq->x == nextq->x) &&
         (compq->y == nextq->y) && (compq->z == nextq->z)) {
	fprintf(stderr, "placeq: Warning, removing identical");
	if(compq->shape == 3) fprintf(stderr, " triangular");
	else if(compq->shape == 4) fprintf(stderr, " quadrilateral");
	else fprintf(stderr, " illegal-shape");
	fprintf(stderr, " panel\n  rmved ctr = (%g %g %g) surf = `%s'", 
		compq->x, compq->y, compq->z, hack_path(compq->surf->name));
	fprintf(stderr, " trans = (%g %g %g)\n", compq->surf->trans[0],
		compq->surf->trans[1], compq->surf->trans[2]);
	fprintf(stderr, "  saved ctr = (%g %g %g) surf = `%s'", 
		nextq->x, nextq->y, nextq->z, hack_path(nextq->surf->name));
	fprintf(stderr, " trans = (%g %g %g)\n", nextq->surf->trans[0],
		nextq->surf->trans[1], nextq->surf->trans[2]);
        /* Remove charge from linked list. */
        for(compq = charges; compq->next != nextq; compq = compq->next) {};
        compq->next = nextq->next;
        nextq = compq;
        compflag = TRUE;
      }
    }
    if(compflag == FALSE) nextc->chgs[nextc->upnumeles[0]++] = nextq;

  }
  return(depth);
}
      
/* 
Place the charges in the cube structure. 
*/
static placeqold(sys, charges)
ssystem *sys;
charge *charges;
{
double length;
double minx, maxx, miny, maxy, minz, maxz, tilelength();
int depth=sys->depth;
int j, k, l, xindex, yindex, zindex, side = sys->side;
int totalq;
charge *nextq;
cube *nextc, *****cubes = sys->cubes;

  nextq = charges;
  minx = maxx = nextq->x;
  miny = maxy = nextq->y;
  minz = maxz = nextq->z;

/* Figure out the length of domain and total number of charges. */
  for(totalq=1, nextq = nextq->next; nextq != NULL;
      totalq++, nextq = nextq->next) {
    maxx = MAX(nextq->x, maxx);
    minx = MIN(nextq->x, minx);
    maxy = MAX(nextq->y, maxy);
    miny = MIN(nextq->y, miny);
    maxz = MAX(nextq->z, maxz);
    minz = MIN(nextq->z, minz);
  }

/* Create the vectors for storing the charges and coefficients. */
  CALLOC(sys->q, totalq+1, double, ON, AMSC);
  CALLOC(sys->p, totalq+1, double, ON, AMSC);

/* Determine side lengths. Make sure domain is eps bigger than problem. */
  length = MAX((maxx - minx), (maxy - miny));
  length = MAX((maxz - minz), length);
  length = (1.01 * length) / side;
  if(length == 0) {
    fprintf(stderr,"placeq: All the lengths in the problem are zero\n");
    exit(1);
  }
  sys->length = length;
  sys->minx = minx;
  sys->miny = miny;
  sys->minz = minz;

/* Count the number of charges per cube. */
  for(nextq = charges; nextq != NULL; nextq = nextq->next) {
    xindex = (nextq->x - minx) / length;
    yindex = (nextq->y - miny) / length;
    zindex = (nextq->z - minz) / length;
    nextc = cubes[depth][xindex][yindex][zindex];
    if(nextc == NULL) {
      CALLOC(nextc, 1, cube, ON, AMSC);
      cubes[depth][xindex][yindex][zindex] = nextc;
      nextc->upnumvects = 1;
      CALLOC(nextc->upnumeles, 1, int, ON, AMSC);
      nextc->upnumeles[0] = 1;
    }
    else {
      nextc->upnumeles[0]++;
    }
  }

/* Allocate space for the charges. */
  for(j=0; j < side; j++) {
    for(k=0; k < side; k++) {
      for(l=0; l < side; l++) {
        nextc = sys->cubes[depth][j][k][l];
        if(nextc != NULL) {  /* Only fill out nonempty cubes. */
	/* Allocate for the charge ptrs, and get q vector pointer. */
	  CALLOC(nextc->chgs, nextc->upnumeles[0], charge*, ON, AMSC);
	  CALLOC(nextc->upnumeles, 1, int, ON, AMSC);
	/* Zero the numchgs to use as index. */
	  nextc->upnumeles[0] = 0;
	}
      }
    }
  }

/* Put the charges in the cube and check to make sure they are not too big. */
  for(totalq = 0, nextq = charges; nextq != NULL; nextq = nextq->next) {
    if(tilelength(nextq) > length) {
/*
      printf("WARNING: A Tile is Larger than the cube containing it.\n");
      printf("Cube length = %g tile length = %g\n", length, tilelength(nextq));
      disfchg(nextq);
*/
    }
    xindex = (nextq->x - minx) / length;
    yindex = (nextq->y - miny) / length;
    zindex = (nextq->z - minz) / length;
    nextc = cubes[depth][xindex][yindex][zindex];
    nextc->chgs[nextc->upnumeles[0]++] = nextq;
  }
}

/*
GetRelations allocates parents links the children. 
*/
getrelations(sys)
ssystem *sys;
{
cube *nextc, *parent, *****cubes = sys->cubes;
int i, j, k, l, side;
  for(i = sys->depth, side = sys->side; i >= 0; i--, side /= 2) {
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
          nextc = cubes[i][j][k][l];
	  if(nextc != NULL) {
	/* Get the parents and children pointers of nonempty cubes. */
	    if(i < sys->depth) {
	      nextc->numkids = 8; /* all cubes, even empties, are counted */
	      CALLOC(nextc->kids, nextc->numkids, cube*, ON, AMSC);
	      nextc->kids[0] = cubes[i+1][2*j][2*k][2*l]; /* empties get */
	      nextc->kids[1] = cubes[i+1][2*j][2*k][2*l+1]; /* null pointers */
	      nextc->kids[2] = cubes[i+1][2*j][2*k+1][2*l];
	      nextc->kids[3] = cubes[i+1][2*j][2*k+1][2*l+1];
	      nextc->kids[4] = cubes[i+1][2*j+1][2*k][2*l];
	      nextc->kids[5] = cubes[i+1][2*j+1][2*k][2*l+1];
	      nextc->kids[6] = cubes[i+1][2*j+1][2*k+1][2*l];
	      nextc->kids[7] = cubes[i+1][2*j+1][2*k+1][2*l+1];
	    }
	    if(i > 0) {
	      parent = cubes[i-1][j/2][k/2][l/2];
	      if(parent == NULL) {
		CALLOC(parent, 1, cube, ON, AMSC);
		cubes[i-1][j/2][k/2][l/2] = parent;
	      }
	      nextc->parent = parent;
	    }
	  }
	}
      }
    }
  }
}

/*
Set the position coordinates of the cubes.
*/
setPosition(sys)
ssystem *sys;
{
int i, j, k, l;
int side = sys->side;
double length = sys->length;
cube *nextc;

/* Mark the position of the lowest level cubes. */
  for(i=sys->depth; i >= 0; i--, side /= 2, length *= 2.0) {
    for(j=0; j < side; j++) {
      for(k=0; k < side; k++) {
	for(l=0; l < side; l++) {
	  nextc = sys->cubes[i][j][k][l];
	  if(nextc != NULL) {
	    nextc->x = length * ((double) j + 0.5) + sys->minx;
	    nextc->y = length * ((double) k + 0.5) + sys->miny;
	    nextc->z = length * ((double) l + 0.5) + sys->minz;
	    nextc->level = i;
	    nextc->j = j;
	    nextc->k = k;
	    nextc->l = l;
	  }
	}
      }
    }
  }
}

/*
Recursive routine to give indexes to the charges so that those in each 
cube are contiguous. In addition, insure that the charges in each parent 
at each level are numbered contiguously.  This is used to support a 
psuedo-adaptive scheme.  Also get the pointer to the appropriate section
of the charge and potential vector.  Uses the eval vector for the potential
coeffs at the lowest level.  Also index the lowest level cubes.
*/
static indexkid(sys, dad, pqindex, pcindex)
ssystem *sys;
cube *dad;
int *pqindex, *pcindex;
{
  int i;
  
  if(dad != NULL) {
    if((dad->numkids == 0) && (dad->upnumvects > 0)) {
      CALLOC(dad->upvects, 1, double*, ON, AMSC);
      CALLOC(dad->nbr_is_dummy, 1, int*, ON, AMSC);
      dad->upvects[0] = &(sys->q[*pqindex]);
      dad->eval = &(sys->p[*pqindex]); /* changed from local to eval 17Feb90 */
      dad->nbr_is_dummy[0] = &(sys->is_dummy[*pqindex]);
      dad->is_dielec = &(sys->is_dielec[*pqindex]);
      dad->index = (*pcindex)++;
      for(i=0; i < dad->upnumeles[0]; i++) {
	(dad->chgs[i])->index = (*pqindex)++;
      }
    }
    else {
      for(i=0; i < dad->numkids; i++) {
	indexkid(sys, dad->kids[i], pqindex, pcindex);
      }
    }
  }
}


/* 
SetExact marks as exact those cubes containing fewer than numterms
number of charges.  If the number of charges in the kids is less than
numterms, the box is marked as exact and the charges are copied up.
In addition, the vector of local expansion coeffs is set to the
potential vector.  Otherwise, the number of nonzero kids is counted
and put in upnumvects as usual.  
*/
/* added 30Mar91: provisions for loc_exact and mul_exact */
setExact(sys, numterms)
ssystem *sys;
int numterms;
{
int i, j, k, l, m, n;
int side = sys->side;
int depth = sys->depth;
int numchgs, num_eval_pnts, first, multerms();
cube *nc, *nkid, *****cubes = sys->cubes;
int all_mul_exact, all_loc_exact, p, num_real_panels;

  for(i=depth; i > 0; 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(i == depth) {
	      ASSERT(nc->upnumvects != 0);
	      /* count the number of true panels in this cube */
	      num_real_panels = 0;
	      for(p = 0; p < nc->upnumeles[0]; p++) {
		if(!nc->chgs[p]->dummy) num_real_panels++;
	      }
	      if(num_real_panels <= numterms) {
		nc->mul_exact = TRUE;
		nc->multisize = nc->upnumeles[0];
	      }
	      else {
		nc->mul_exact = FALSE;
		nc->multisize = multerms(sys->order);
	      }
	      if(nc->upnumeles[0] <= numterms) {
		nc->loc_exact = TRUE;
		nc->localsize = nc->upnumeles[0];
	      }
	      else {
		nc->loc_exact = FALSE;
		nc->localsize = multerms(sys->order); 
	      }
	    }
	    else {  
	      /* Count the number of charges and nonempty kids. */
	      all_loc_exact = all_mul_exact = TRUE;
	      num_eval_pnts = numchgs = nc->upnumvects = 0;
	      for(m = 0; m < nc->numkids; m++) {
		nkid = nc->kids[m];
		if(nkid != NULL) {
		  nc->upnumvects += 1;
		  if(nkid->mul_exact == FALSE) all_mul_exact = FALSE;
		  else {
		    num_eval_pnts += nkid->upnumeles[0];
		    for(p = 0; p < nkid->upnumeles[0]; p++) {
		      if(!nkid->chgs[p]->dummy) numchgs++;
		    }
		  }
		  if(nkid->loc_exact == FALSE) all_loc_exact = FALSE;
		}
	      }
	      /* If all nonempty kids exact, # chgs <= # terms, mark exact, 
		 copy chgs, and promote pointers to charge and potential.  
		 Note EXPLOITS special ordering of the pot and charge vectors.
		 */
	      if(!all_mul_exact || (numchgs > numterms)) { /* multi req'd */
		nc->mul_exact = FALSE;
		nc->multisize = multerms(sys->order);
	      }
	      else if(all_mul_exact && (numchgs <= numterms)) { 
		nc->mul_exact = TRUE;
		nc->upnumvects = 1;
		CALLOC(nc->upvects, 1, double*, ON, AMSC);
		CALLOC(nc->upnumeles, 1, int, ON, AMSC);
		nc->upnumeles[0] = num_eval_pnts; /* was numchgs 30Mar91 */
		nc->multisize = num_eval_pnts; /* was numchgs */
		CALLOC(nc->chgs, num_eval_pnts, charge*, ON, AMSC);
		num_eval_pnts = 0;
		for(m=0, first=TRUE; m < nc->numkids; m++) {
		  nkid = nc->kids[m]; 
		  if(nkid != NULL) {
		    if(first == TRUE) {
		      nc->upvects[0] = nkid->upvects[0];
		      if(nc->nbr_is_dummy == NULL)
			  CALLOC(nc->nbr_is_dummy, 1, int*, ON, AMSC);
		      nc->nbr_is_dummy[0] = nkid->nbr_is_dummy[0];
		      first = FALSE;
		    }
		    for(n=0; n < nkid->upnumeles[0]; n++) {
		      nc->chgs[num_eval_pnts++] = nkid->chgs[n];
		    }
		  }
		}
	      }

	      /* do the same for local expansion */
	      /* if local exact, must be multi exact => no promotion reqd */
	      if(!all_loc_exact || (num_eval_pnts > numterms)) { /* le req'd */
		nc->loc_exact = FALSE;
		nc->localsize = multerms(sys->order);
	      }
	      else if(all_loc_exact && (num_eval_pnts <= numterms)) { 
		nc->loc_exact = TRUE;
		nc->localsize = num_eval_pnts;
	      }
	    }
	  }
	}
      }
    }
  }
}