mulsetup.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"

cube *cstack[1024];		/* Stack used in several routines. */

static int placeq();
static getnbrs();
static indexkid();
static linkcubes();
static setMaxq();
static getAllInter();
static set_vector_masks();

/*
  sets up the partitioning of space and room for charges and expansions
*/
ssystem *mulInit(autom, depth, order, charges)
int autom;			/* ON => choose depth automatically */
int depth, order;
charge *charges;
{
  ssystem *sys;
  int qindex=1, cindex=1;

  CALLOC(sys, 1, ssystem, ON, AMSC);
  sys->depth = depth;		/* overwritten below if autom = ON */
  sys->order = order;
  
  sys->depth = placeq(autom, sys, charges); /* create cubes, put in charges */

  getrelations(sys);		/* Get all the prnts and kids for each cube. */

  setPosition(sys);		/* Figures out position of cube center. */
  indexkid(sys, sys->cubes[0][0][0][0], &qindex, &cindex); 
				/* Index chgs and cubes. */

#if ADAPT == ON
  setExact(sys, multerms(sys->order)); /* Note cubes to be done exactly and
					   determine the number of nonempty
					   kids for each cube. */
#else
  setExact(sys, 0);		/* Note cubes to be done exactly and
					   determine the number of nonempty
					   kids for each cube. */
#endif

  getnbrs(sys);			/* Get all the nearest neighbors. At bot level
				   add as nearest nbrs cubes in exact block. */
  linkcubes(sys);		/* Make linked-lists of direct, multis, and
				   locals to do at each level. */
  set_vector_masks(sys);	/* set up sys->is_dummy and sys->is_dielec */
  setMaxq(sys);                 /* Calculates the max # chgs in cubes treated
				   exactly, and over lowest level cubes. */
  getAllInter(sys);		/* Get the interaction lists at all levels. */

  return(sys);
}

/*
  places the charges using best number of levels for adaptive algorithm
  - returns number of levels used
  - can be called with flag == OFF to set depth = sys->depth
  - uses entire panel list, dieletric and dummy panels are included
     (excluding BOTH or DIELELC type panels if NUMDPT == 2, ie 
     two-point flux-density-differences are being use ),
     so the constraining type of exactness is usually
     local expansion exactness---using this method the switch to multipole
     exactness is done automatically if there are no dielectric panels
  - when the method without dielectric panels is set up, these loops will
     need to ignore all dielectric panels
  - this routine is still called to set automatic levels if ADAPT is OFF,
     ie even when the calculation is not adaptive, so results can be compared
*/
static int placeq(flag, sys, charges)
int flag;			/* ON => set depth automatically */
ssystem *sys;
charge *charges;
{
  int i, j, k, l, side, totalq, isexact, multerms(), depth;
  int xindex, yindex, zindex, limit = multerms(sys->order), compflag;
  int exact_cubes_this_level, cubes_this_level;
  double length0, length, exact_ratio;
  double minx, maxx, miny, maxy, minz, maxz, tilelength(), maxTileLength;
  charge *nextq, *compq;
  cube *****cubes, *nextc;
  char *hack_path();

  /* Figure out the length of lev 0 cube and total number of charges. */
  nextq = charges;
  minx = maxx = nextq->x;
  miny = maxy = nextq->y;
  minz = maxz = nextq->z;

  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);
  }

  sys->minx = minx;
  sys->miny = miny;
  sys->minz = minz;

  /* Make sure cube isn't smaller than a tile. */
  for(maxTileLength = 0.0, nextq = charges;
      nextq != NULL; nextq = nextq->next) {
    length = tilelength(nextq);
    maxTileLength = MAX(maxTileLength, length);
  }
  if((maxx - minx) < maxTileLength) {
    maxx += 0.5 * maxTileLength;
    minx -= 0.5 * maxTileLength;
  }
  if((maxy - miny) < maxTileLength) {
    maxy += 0.5 * maxTileLength;
    miny -= 0.5 * maxTileLength;
  }
  if((maxz - minz) < maxTileLength) {
    maxz += 0.5 * maxTileLength;
    minz -= 0.5 * maxTileLength;
  }

  /* (see below for test for panel size vs cube size) */

  length0 = MAX((maxx - minx), (maxy - miny));
  length0 = MAX((maxz - minz), length0);

  /* 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);

  /* set up mask vector: is_dummy[i] = TRUE => panel i is a dummy */
  CALLOC(sys->is_dummy, totalq+1, int, ON, AMSC);

  /* set up mask vector: is_dielec[i] = TRUE => panel i is on DIELEC or BOTH */
  CALLOC(sys->is_dielec, totalq+1, int, ON, AMSC);

  if(flag == ON) {		/* set depth of partitions automatically */
    /* alloc spine of cube pntr array - leave enough room for depth = MAXDEP */
    CALLOC(cubes, MAXDEP+1, cube****, ON, AMSC); 

    /* allocate for levels 0, 1, and 2 (always used) */
    for(side = 1, i=0; i <= 2; side *= 2, i++) {
      CALLOC(cubes[i], side, cube***, ON, AMSC);
      for(j=0; j < side; j++) {
	CALLOC(cubes[i][j], side, cube**, ON, AMSC);
	for(k=0; k < side; k++) {
	  CALLOC(cubes[i][j][k], side, cube*, ON, AMSC);
	}
      }
    }
    /* side /= 2; */

    /* for each level > 2: allocate for full cubes, count panels in each 
       - quit loop if all lowest level cubes are exact */
    for(isexact = FALSE; isexact == FALSE; side *= 2, i++) {

      if(i > MAXDEP) {
	fprintf(stderr, 
		"placeq: out of cube pntr space - increase MAXDEP == %d\n", 
		MAXDEP);
	exit(1);
      }

      length = (1.01 * length0)/side;

      CALLOC(cubes[i], side, cube***, OFF, AMSC);
      if(cubes[i] == NULL) {
	fprintf(stderr, "placeq: %d levels set up\n", i-1);
	exit(1);
      }
      for(j=0; j < side; j++) {
	CALLOC(cubes[i][j], side, cube**, OFF, AMSC);
	if(cubes[i][j] == NULL) {
	  fprintf(stderr, "placeq: %d levels set up\n", i-1);
	  exit(1);
	}
	for(k=0; k < side; k++) {
	  CALLOC(cubes[i][j][k], side, cube*, OFF, AMSC);
	  if(cubes[i][j][k] == NULL) {
	    fprintf(stderr, "placeq: %d levels set up\n", i-1);
	    exit(1);
	  }
	}
      }

      /* Count the number of charges per cube and allocate if needed */
      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[i][xindex][yindex][zindex];
	if(nextc == NULL) {
	  CALLOC(nextc, 1, cube, OFF, AMSC);
	  if(nextc == NULL) {
	    fprintf(stderr, "placeq: %d levels set up\n", i-1);
	    exit(1);
	  }
	  cubes[i][xindex][yindex][zindex] = nextc;
	  nextc->upnumvects = 1;
	  CALLOC(nextc->upnumeles, 1, int, OFF, AMSC);
	  if(nextc->upnumeles == NULL) {
	    fprintf(stderr, "placeq: %d levels set up\n", i-1);
	    exit(1);
	  }
	  nextc->upnumeles[0] = 1;
	}
	else {
	  nextc->upnumeles[0]++;
	}
      }

      /* if the current lowest level is not exact, loop back until it is */
      /*    check for exactness of this level, get cube statistics */
      isexact = TRUE;
      cubes_this_level = 0;
      exact_cubes_this_level = 0;
      for(j = 0; j < side; j++) {
	for(k = 0; k < side; k++) {
	  for(l = 0; l < side; l++) {
	    if(cubes[i][j][k][l] != NULL) {
	      if(cubes[i][j][k][l]->upnumeles[0] > limit) {
		isexact = FALSE;
	      }
	      else exact_cubes_this_level++;
	      cubes_this_level++;
	    }
	  }
	}
      }

      /*    decide whether to go down another level by checking exact/ttl */
      exact_ratio = (double)exact_cubes_this_level/(double)cubes_this_level;
      if(exact_ratio > EXRTSH) 
	  isexact = TRUE;      	/* set up to terminate level build loop */
      /* fprintf(stdout, "Level %d, %g%% exact\n", i, exact_ratio*100.0); */

      /* clean up cube structs if need to go down another level */
      if(isexact == FALSE) {
	for(j = 0; j < side; j++) {
	  for(k = 0; k < side; k++) {
	    for(l = 0; l < side; l++) {
	      if(cubes[i][j][k][l] != NULL) {
		cubes[i][j][k][l]->upnumeles[0] = 0;
		cubes[i][j][k][l]->upnumvects = 0;
	      }
	    }
	  }
	}
      }
    }
    depth = i - 1;		/* the automatically set depth */
    side /= 2;
  }
  else {			/* old code - uses sys->depth for depth */
    /* Allocate the cubes, note calloc used because zeros everything. */
    depth = sys->depth;
    CALLOC(cubes, sys->depth+1, cube****, ON, AMSC);
    for(side = 1, i=0; i <= depth; side *= 2, i++) {
      CALLOC(cubes[i], side, cube***, ON, AMSC);
      for(j=0; j < side; j++) {
	CALLOC(cubes[i][j], side, cube**, ON, AMSC);
	for(k=0; k < side; k++) {
	  CALLOC(cubes[i][j][k], side, cube*, ON, AMSC);
	}
      }
    }
    side /= 2;
    length = (1.01 * length0)/side;

    /* 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]++;
      }
    }
  }
  sys->length = length;
  sys->side = side;
  sys->cubes = cubes;

  /* 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 cubes; check to make sure they are not too big. */
  for(nextq = charges; nextq != NULL; nextq = nextq->next) {
#if 1 == 0
    if(tilelength(nextq) > length) {
      fprintf(stderr,
	      "\nplaceq: Warning, a panel is larger than the cube supposedly containing it\n");
      fprintf(stderr,"  cube length = %g panel length = %g\n", 
	      length, tilelength(nextq));
      /* disfchg(nextq); */
    }
#endif
    xindex = (nextq->x - minx) / length;
    yindex = (nextq->y - miny) / length;
    zindex = (nextq->z - minz) / length;
    nextc = cubes[depth][xindex][yindex][zindex];