mulmats.c

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

	    for(j = nnsize - 1; j >= 0; j--) {
	      mat[i][offset + j] = nmat[i][j];
	    }
	  }
	  else {
#if DPCOMP == ON
	    fprintf(stdout, "Source mat, nnbr to nc\n");
	    dumpMat(nmat, nsize, nnsize);
#endif
	    find_flux_density_row(mat, nmat, i, nnsize, nsize, 0, offset,
				  nc_pc, nnbr_pc, nc_dummy, nnbr_dummy);
	  }
	}
	/* Get the row of the big matrix associated with this nnbr. */
	for(noffset = 0, l = -1; l < nc->numnbrs; l++) { /* lp on nc's nbrs */
	  if(l < 0) nnnbr = nc;
	  else nnnbr = nc->nbrs[l];
	  if(NEAR(nnnbr, nj, nk, nl)) {  /* Note, near to nc!! */
	    if(nnbr == nnnbr) m = -1;
	    else { /* Find this nnnbr's position in nnbr's list */
	      for(m=0; m < nnbr->numnbrs; m++) {
		if(nnbr->nbrs[m] == nnnbr) break;
	      }
	      ASSERT(m < nnbr->numnbrs);
	    }
	    nnnsize = nnbr->directnumeles[m+1];
	    nmat = nnbr->directmats[m+1];
	    ASSERT(nnbr->directnumeles[m+1] == nnnbr->directnumeles[0]);
	    nnnbr_pc = nnnbr->chgs; /* panels in nnnbr */
	    nnnbr_dummy = nnnbr->nbr_is_dummy[0];
#if CHKDUM == ON
	    chkDummyList(nnnbr_pc, nnnbr_dummy, nnnsize);
#endif
	    for(i = nnsize - 1; i >= 0; i--) { /* loop on panels in nnbr */
	      if(nnbr_dummy[i]) continue;
	      if(nnbr_pc[i]->surf->type != DIELEC) {
		for(j = nnnsize - 1; j >= 0; j--) {
		  mat[offset + i][noffset+j] = nmat[i][j];
		}
	      }
	      else {
#if DPCOMP == ON
		fprintf(stdout, "Source mat, nnnbr to nnbr\n");
		dumpMat(nmat, nnsize, nnnsize);
#endif
		find_flux_density_row(mat, nmat, i, nnnsize, nnsize, offset,
				      noffset, nnbr_pc, nnnbr_pc, nnbr_dummy,
				      nnnbr_dummy);
	      }
	    }
	    noffset += nnnsize;
	  }
	}
	offset += nnsize;
      }
    }

    /* set up the local is_dummy vector for the rows/cols of mat */
    /* THIS COULD BE AVOIDED BY USING CUBE is_dummy's INSIDE invert() */
    if(big_mat_size < offset) {	/* allocate only if larger array needed */
      CALLOC(is_dummy, offset, int, ON, AMSC);
    }
    /* dump sections of the dummy vector in order cubes appear in nbr lst */
    /* (use fragment of Jacob's loop above) */
    nnnsize = noffset = nc->directnumeles[0];
    nc_dummy = nc->nbr_is_dummy[0];
    for(i = nnnsize - 1; i >= 0; i--) {
      is_dummy[i] = nc_dummy[i];
    }
    for(l = 0; l < nc->numnbrs; l++) {
      nnnbr = nc->nbrs[l];
      if(NEAR(nnnbr, nj, nk, nl)) {
	nnnsize = nnnbr->directnumeles[0];
	nc_dummy = nnnbr->nbr_is_dummy[0];
	for(i = nnnsize - 1; i >= 0; i--) {
	  is_dummy[i + noffset] = nc_dummy[i];
	}
	noffset += nnnsize;
      }
    }

    /* The big Matrix is filled in, invert it and get the preconditioner. */
#if DPCOMP == ON
    fprintf(stdout, "Before compression\n");
    dumpMat(mat, offset, offset);
#endif
    nnnsize = compressMat(mat, offset, is_dummy, BOTH);
#if DPCOMP == ON
    fprintf(stdout, "After compression\n");
    dumpMat(mat, nnnsize, nnnsize);
#endif
    invert(mat, nnnsize, NULL);
    expandMat(mat, offset, nnnsize, is_dummy, BOTH);
#if DPCOMP == ON
    fprintf(stdout, "After expansion\n");
    dumpMat(mat, offset, offset);
#endif

    /* Copy out the preconditioner to the saved matrices. */
    for(i = nsize - 1; i >= 0; i--) {
      for(j = nsize - 1; j >= 0; j--) {
	nc->precondmats[0][i][j] = mat[i][j];
      }
    }
    offset = nsize;
    for(k=0; k < nc->numnbrs; k++) {
      nnbr = nc->nbrs[k];
      if(NEAR(nnbr, nj, nk, nl)) {
	nnsize = nc->directnumeles[k+1];
	nmat = nc->precondmats[k+1];
	for(i = nsize - 1; i >= 0; i--) {
	  for(j = nnsize - 1; j >= 0; j--) {
	    nmat[i][j] = mat[i][offset + j];
	  }
	}
	offset += nnsize;
      }
      else nc->precondmats[k+1] = NULL;
    }
  }
}

/*
  finds a row of flux density coeffs from three potential coeff rows
  - to_mat[eval_row][] is the destination row; from_mat[eval_row][]
    initially contains the potential coefficients for evals at the 
    center of eval_panels[eval_row] (unless NUMDPT == 2, is garbage then)
  - the eval panels are scaned until eval_panels[eval_row]'s
    dummies are found and the corresponding two rows are identified
  - the divided differences built with entries in the same columns in
    these three rows replace the to_mat[eval_row][] entries:
    to_mat[eval_row][j] = a1*from_mat[eval_row][j]
              + a2*from_mat[pos_dum_row][j] + a3*from_mat[neg_dum_row][j]
  - if a dummy panel is not found in the panel list, its row is generated
    using explicit calcp() calls (shouldn't happen much)
  - global flags used here
    NUMDPT = number of divided diff points, 2 or 3
    SKIPDQ = ON=>don't do cancellation-prone add-subtract of identical
      influence of DIELEC/BOTH panels' charges on dummy panel pot. evals
*/
find_flux_density_row(to_mat, from_mat, eval_row, n_chg, n_eval, row_offset,
		      col_offset, eval_panels, chg_panels, eval_is_dummy, 
		      chg_is_dummy)
double **to_mat, **from_mat;
int eval_row, n_chg, n_eval, row_offset, col_offset;
charge **eval_panels, **chg_panels;
int *eval_is_dummy, *chg_is_dummy;
{
  int dindex, j;
  double factor, calcp();
  charge *dp;
  surface *surf = eval_panels[eval_row]->surf;

  /* do divided difference w/ three rows to get dielectric row */
#if NUMDPT == 3
  /* - dielectric panel row first */
  factor = -(surf->outer_perm + surf->inner_perm)/
      (eval_panels[eval_row]->pos_dummy->area);
#if DPDDIF == ON
  fprintf(stdout, "Center row, factor = %g\n", factor);
#endif
  for(j = n_chg - 1; j >= 0; j--) { /* loop on columns */
    if(!chg_is_dummy[j]) 
	to_mat[row_offset + eval_row][col_offset + j] 
	    = from_mat[eval_row][j]*factor;
#if DPDDIF == ON
    fprintf(stdout, " %.16e", from_mat[eval_row][j]);
#endif				/* #if DPDDIF == ON */
  }
#endif				/* #if NUMDPT == 3 */
  /* - do positive dummy row */
  /*   first find the dummy row */
  dindex = -1;
  dp = eval_panels[eval_row]->pos_dummy; /* get dummy panel from eval panel */
  for(j = n_eval - 1; j >= 0; j--) {
    if(!eval_is_dummy[j]) continue;
    if(dp == eval_panels[j]) {
      dindex = j;
      break;
    }
  }
  if(dindex != -1) { /* dummy row found */
#if NUMDPT == 3
    factor = surf->outer_perm/eval_panels[dindex]->area;
#else
    /* this is the only factor required for two dummy rows in two point case */
    factor = (surf->inner_perm - surf->outer_perm)
	/(eval_panels[eval_row]->neg_dummy->area 
	  + eval_panels[eval_row]->pos_dummy->area);
#endif
#if DPDDIF == ON
    fprintf(stdout, "\nPos dummy row, factor = %g\n", factor);
#endif
    for(j = n_chg - 1; j >= 0; j--) {
#if SKIPDQ == ON
      if(chg_panels[j]->index == eval_panels[eval_row]->index) {
	to_mat[row_offset + eval_row][col_offset + j] = 0.0;
	continue;
      }
#endif
      if(!chg_is_dummy[j])
#if NUMDPT == 3
	  to_mat[row_offset + eval_row][col_offset + j] 
	      += from_mat[dindex][j]*factor;
#else				/* make sure to overwrite possible garbage */
	  to_mat[row_offset + eval_row][col_offset + j] 
	      = -from_mat[dindex][j]*factor;
#endif
#if DPDDIF == ON
      fprintf(stdout, " %.16e (%d)", from_mat[dindex][j],chg_panels[j]->index);
#endif
    }
  }
  else {		/* dummy row out of cube => build it w/calcp */
#if NUMDPT == 3
    factor = surf->outer_perm/dp->area;
#else
    /* this is the only factor required for two dummy rows in two point case */
    factor = (surf->inner_perm - surf->outer_perm)
	/(eval_panels[eval_row]->neg_dummy->area 
	  + eval_panels[eval_row]->pos_dummy->area);
#endif
#if DPDDIF == ON
    fprintf(stdout, "\nPos dummy calcp row, factor = %g\n", factor);
#else
    fprintf(stderr, "\nolmulMatPrecond: building pos. dummy row\n");
#endif
    for(j = n_chg - 1; j >= 0; j--) {
#if SKIPQD == ON
      if(chg_panels[j]->index == eval_panels[eval_row]->index) {
	to_mat[row_offset + eval_row][col_offset + j] = 0.0;
	continue;
      }
#endif
      if(!chg_is_dummy[j]) {
#if NUMDPT == 3
	to_mat[row_offset + eval_row][col_offset + j] 
	    += calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL)*factor;
#else
	to_mat[row_offset + eval_row][col_offset + j] 
	    = -calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL)*factor;
#endif
#if DPDDIF == ON
	fprintf(stdout, " %.16e (%d)",
		calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL),
		chg_panels[j]->index);
      }
      else {
	fprintf(stdout, " dummy");
#endif
      }
    }
  }
  /* - do negative dummy row */
  /*   first find the dummy row */
  dindex = -1;
  dp = eval_panels[eval_row]->neg_dummy; /* get dummy panel from eval panel */
  for(j = n_eval - 1; j >= 0; j--) {
    if(!eval_is_dummy[j]) continue;
    if(dp == eval_panels[j]) {
      dindex = j;
      break;
    }
  }
  if(dindex != -1) { /* dummy row found */
#if NUMDPT == 3
    factor = surf->inner_perm/eval_panels[dindex]->area;
#endif
#if DPDDIF == ON
    fprintf(stdout, "\nNeg dummy row, factor = %g\n", factor);
#endif
    for(j = n_chg - 1; j >= 0; j--) {
#if SKIPQD == ON
      if(chg_panels[j]->index == eval_panels[eval_row]->index) continue;
#endif
      if(!chg_is_dummy[j])
	  to_mat[row_offset + eval_row][col_offset + j] 
	      += from_mat[dindex][j]*factor;
#if DPDDIF == ON
      fprintf(stdout, " %.16e (%d)", from_mat[dindex][j],chg_panels[j]->index);
#endif
    }
  }
  else {		/* dummy row out of cube => build it w/calcp */
    factor = surf->inner_perm/dp->area;
#if DPDDIF == ON
    fprintf(stdout, "\nNeg dummy calcp row, factor = %g\n", factor);
#else
    fprintf(stderr, "olmulMatPrecond: building neg. dummy row\n");
#endif
    for(j = n_chg - 1; j >= 0; j--) {
#if SKIPQD == ON
      if(chg_panels[j]->index == eval_panels[eval_row]->index) continue;
#endif
      if(!chg_is_dummy[j]) {
	to_mat[row_offset + eval_row][col_offset + j] 
	    += calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL)*factor;
#if DPDDIF == ON
	fprintf(stdout, " %.16e (%d)",
		calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL),
		chg_panels[j]->index);
      }
      else {
	fprintf(stdout, " dummy");
#endif
      }
    }
  }
#if NUMDPT == 2
  /* - do row entry due to panel contribution 
     - entry only necessary if eval panel is in chg panel list */
  /*   search for the eval panel in the charge panel list */
  dp = NULL;
  for(j = n_chg - 1; j >= 0; j--) {
    if(!chg_is_dummy[j]) {
      if(eval_panels[eval_row] == chg_panels[j]) {
	dp = eval_panels[eval_row];
	break;
      }
    }
  }
  /*   set entry if eval panel found in chg panel list 
       - this is an overwrite; contributions of other rows should cancel */
  if(dp != NULL) {
    to_mat[row_offset + eval_row][col_offset + j]
	= -(2*M_PI*(surf->inner_perm + surf->outer_perm)
	    /eval_panels[eval_row]->area);
  }
#endif

#if DPDDIF == ON
  fprintf(stdout, "\nDivided difference row (%d)\n", 
	  eval_panels[eval_row]->index);
  for(j = n_chg - 1; j >= 0; j--) {
    fprintf(stdout, " %.16e (%d)", 
	    to_mat[row_offset + eval_row][col_offset + j], 
	    chg_panels[j]->index);
  }
  fprintf(stdout, "\n\n");
#endif
}


/* 
MulMatUp computes the multipole to multipole or charge to
multipole matrices that map to a parent's multipole coeffs from its
children's multipoles or charges. Note that only one set of
multipole to multipole matrices is computed per level by exploiting the
uniform break-up of three-space (ie many shifts have similar geometries).  
*/
mulMatUp(sys) 
ssystem *sys; 
{
cube *nextc, *kid;
int i, j, numterms, depth, order = sys->order;
double **multimats[8];