mulmats_old.c

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

#if DPPCMP == ON
	fprintf(stdout, "Self mat, nc = (%d %d %d), nsize = %d\n",
		nc->j, nc->k, nc->l, nsize);
	dumpMat(nc->directmats[0], nsize, nsize);
#endif
      }
      else {
#if DPCOMP == ON
	fprintf(stdout, "Source mat, nc to nc\n");
	dumpMat(nc->directmats[0], nsize, nsize);
#endif
	find_flux_density_row(mat, nc->directmats[0], i, nsize, nsize, 0, 0,
			      nc_pc, nc_pc, nc_dummy, nc_dummy);
      }	
    }
    offset = nsize;
    for(k=0; k < nc->numnbrs; k++) { /* loop on neighbors of nc */
      nnbr = nc->nbrs[k];
      if(NEAR(nnbr, nj, nk, nl)) {
	nnsize = nc->directnumeles[k+1];
	nmat = nc->directmats[k+1];
	ASSERT(nc->directnumeles[k+1] == nnbr->directnumeles[0]);
	nnbr_dummy = nnbr->nbr_is_dummy[0];
	nnbr_pc = nnbr->chgs;
#if CHKDUM == ON
	chkDummyList(nnbr_pc, nnbr_dummy, nnsize);
#endif
	for(i = nsize - 1; i >= 0; i--) {
	  if(nc_dummy[i]) continue;
	  if(nc_pc[i]->surf->type != DIELEC) {
	    for(j = nnsize - 1; j >= 0; j--) {
	      mat[i][offset + j] = nmat[i][j];
	    }
#if DPPCMP == ON
	    fprintf(stdout, 
      "Nbr mat, nc = (%d %d %d), nnbr = (%d %d %d), nsize = %d, nnsize = %d\n",
	        nc->j, nc->k, nc->l, nnbr->j, nnbr->k, nnbr->l, nsize, nnsize);
	    dumpMat(nmat, nsize, nnsize);
#endif
	  }
	  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];
		}
#if DPPCMP == ON
		fprintf(stdout, 
  "Nbr mat, nnbr = (%d %d %d), nnnbr = (%d %d %d), nnsize = %d, nnnsize = %d\n",
			nnbr->j, nnbr->k, nnbr->l, 
			nnnbr->j, nnnbr->k, nnnbr->l, nnsize, nnnsize);
		dumpMat(nmat, nnsize, nnnsize);
#endif
	      }
	      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 contians the potential coefficients for evals at the 
    center of eval_panels[eval_row]
  - 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)
*/
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 */
  /* - 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, " %.3g", from_mat[eval_row][j]);
#endif
  }
  /* - 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 */
    factor = surf->outer_perm/eval_panels[dindex]->area;
#if DPDDIF == ON
    fprintf(stdout, "\nPos dummy row, factor = %g\n", factor);
#endif
    for(j = n_chg - 1; j >= 0; j--) {
      if(!chg_is_dummy[j])
	  to_mat[row_offset + eval_row][col_offset + j] 
	      += from_mat[dindex][j]*factor;
#if DPDDIF == ON
      fprintf(stdout, " %.3g", from_mat[dindex][j]);
#endif
    }
  }
  else {		/* dummy row out of cube => build it w/calcp */
    factor = surf->outer_perm/dp->area;
#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(!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, " %.3g",
		calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL));
      }
      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 */
    factor = surf->inner_perm/eval_panels[dindex]->area;
#if DPDDIF == ON
    fprintf(stdout, "\nNeg dummy row, factor = %g\n", factor);
#endif
    for(j = n_chg - 1; j >= 0; j--) {
      if(!chg_is_dummy[j])
	  to_mat[row_offset + eval_row][col_offset + j] 
	      += from_mat[dindex][j]*factor;
#if DPDDIF == ON
      fprintf(stdout, " %.3g", from_mat[dindex][j]);
#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(!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, " %.3g",
		calcp(chg_panels[j], dp->x, dp->y, dp->z, NULL));
      }
      else {
	fprintf(stdout, " dummy");
#endif
      }
    }
  }
#if DPDDIF == ON
  fprintf(stdout, "\nDivided difference row\n");
  for(j = n_chg - 1; j >= 0; j--) {
    fprintf(stdout, " %.3g", to_mat[row_offset + eval_row][col_offset + j]);
  }
  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];

#if OFF == ON			/* OFF == OFF produces the M2M error!! */
    for(i=0; i < 8; i++) multimats[i] = NULL;
#endif

  numterms = multerms(order);