boundaries.cpp

来自mit的fdtd开放性源代码meep

}

void fields_chunk::zero_metal(field_type ft) {
  for (int i=0;i<num_zeroes[ft];i++) *(zeroes[ft][i]) = 0.0;
}

void fields::find_metals() {
  for (int i=0;i<num_chunks;i++)
    if (chunks[i]->is_mine()) {
      const volume vi = chunks[i]->v;
      FOR_FIELD_TYPES(ft) {
        delete[] chunks[i]->zeroes[ft];
        // First electric components...
        chunks[i]->num_zeroes[ft] = 0;
        DOCMP FOR_COMPONENTS(c)
          if (type(c) == ft && chunks[i]->f[c][cmp])
	    LOOP_OVER_VOL_OWNED(vi, c, n) 
	      if (IVEC_LOOP_AT_BOUNDARY) { // todo: just loop over boundaries
		IVEC_LOOP_ILOC(vi, here);
		if (on_metal_boundary(here))
		  chunks[i]->num_zeroes[ft]++;
	      }
        typedef double *double_ptr;
        chunks[i]->zeroes[ft] = new double_ptr[chunks[i]->num_zeroes[ft]];
        int num = 0;
        DOCMP FOR_COMPONENTS(c)
          if (type(c) == ft && chunks[i]->f[c][cmp])
	    LOOP_OVER_VOL_OWNED(vi, c, n)
	      if (IVEC_LOOP_AT_BOUNDARY) { // todo: just loop over boundaries
		IVEC_LOOP_ILOC(vi, here);
		if (on_metal_boundary(here))
		  chunks[i]->zeroes[ft][num++] = &(chunks[i]->f[c][cmp][n]);
	      }
      }
    }
}

void fields::connect_the_chunks() {
  int *nc[NUM_FIELD_TYPES][3][2];
  FOR_FIELD_TYPES(f)
    for (int ip=0;ip<3;ip++)
      for (int io=0;io<2;io++) {
	nc[f][ip][io] = new int[num_chunks];
	for (int i=0;i<num_chunks;i++) nc[f][ip][io][i] = 0;
      }

  /* For some of the chunks, H==B, and we definitely don't need to
     send B between two such chunks.   We'll still send B when
     the recipient has H != B, since the recipient needs to get B
     from somewhere (although it could get it locally, in principle).
     When the sender has H != B, we'll skip sending B (we'll only send H)
     since we need to get the correct curl H in the E update.  This is
     a bit subtle since the non-owned B may be different from H even
     on an H==B chunk (true?), but since we don't use the non-owned B
     for anything(?) it shouldn't matter. */
  int *B_redundant = new int[num_chunks*2 * 5];
  for (int i = 0; i < num_chunks; ++i)
    FOR_H_AND_B(hc,bc)
      B_redundant[5*(num_chunks+i) + bc-Bx] 
      = chunks[i]->f[hc][0] == chunks[i]->f[bc][0];
  and_to_all(B_redundant + 5*num_chunks, B_redundant, 5*num_chunks);

  for (int i=0;i<num_chunks;i++) {
    // First count the border elements...
    const volume vi = chunks[i]->v;
    FOR_FIELD_TYPES(ft)
      for (int ip=0;ip<3;ip++)
	for (int j=0;j<num_chunks;j++)
	  comm_sizes[ft][ip][j+i*num_chunks] = 0;
    FOR_COMPONENTS(corig) {
      if (have_component(corig))
	LOOP_OVER_VOL_NOTOWNED(vi, corig, n) {
	  IVEC_LOOP_ILOC(vi, here);
	  component c = corig;
	  // We're looking at a border element...
	  complex<double> thephase;
	  if (locate_component_point(&c,&here,&thephase)
	      && !on_metal_boundary(here))
	    for (int j=0;j<num_chunks;j++) {
	      if ((chunks[i]->is_mine() || chunks[j]->is_mine())
		  && chunks[j]->v.owns(here)
		  && !(is_B(corig) && is_B(c) &&
		       B_redundant[5*i+corig-Bx] && B_redundant[5*j+c-Bx])) {
		const int pair = j+i*num_chunks;
		const connect_phase ip = thephase == 1.0 ? CONNECT_COPY 
		  : (thephase == -1.0 ? CONNECT_NEGATE : CONNECT_PHASE);
		{
		  const int nn = is_real?1:2;
		  nc[type(corig)][ip][Incoming][i] += nn;
		  nc[type(c)][ip][Outgoing][j] += nn;
		  comm_sizes[type(c)][ip][pair] += nn;
		}
		if (is_electric(corig)) {
		  int common_pols = 0;
		  for (polarization *pi = chunks[i]->pol; pi; pi = pi->next)
		    for (polarization *pj = chunks[j]->pol; pj; pj = pj->next)
		      if (pi->pb->get_identifier() == pj->pb->get_identifier())
			common_pols += 1;
		  const int nn = (is_real?1:2) * common_pols * 2;
		  nc[P_stuff][ip][Incoming][i] += nn;
		  nc[P_stuff][ip][Outgoing][j] += nn;
		  comm_sizes[P_stuff][ip][pair] += nn;
		  // Note above that the factor of two in 2*nn comes from
		  // the fact that we have two polarization arrays, pol and
		  // olpol.
		}
	      } // if is_mine and owns...
	    } // loop over j chunks
        } // LOOP_OVER_VOL_NOTOWNED
    } // FOR_COMPONENTS

    // Allocating comm blocks as we go...
    FOR_FIELD_TYPES(ft)
      for (int j=0;j<num_chunks;j++) {
        delete[] comm_blocks[ft][j+i*num_chunks];
        comm_blocks[ft][j+i*num_chunks] =
          new double[comm_size_tot(ft, j+i*num_chunks)];
      }
  } // loop over i chunks

  /* Note that the ordering of the connections arrays must be kept
     consistent with the fields::step_boundaries.  In particular, we
     must set up the connections array so that all of the connections
     for process i come before all of the connections for process i'
     for i < i' */

  // wh stores the current indices in the connections array(s)
  int *wh[NUM_FIELD_TYPES][3][2];

  /* Now allocate the connection arrays... this is still slightly
     wasteful (by a factor of 2) because we allocate for chunks we
     don't own if we have a connection with them. Removing this waste
     would mean a bunch more is_mine() checks below. */
  FOR_FIELD_TYPES(f)
    for (int ip=0;ip<3;ip++) {
      for (int io=0;io<2;io++) {
	for (int i=0;i<num_chunks;i++)
	  chunks[i]->alloc_extra_connections(field_type(f),
					     connect_phase(ip),
					     in_or_out(io),
					     nc[f][ip][io][i]);
	delete[] nc[f][ip][io];
	wh[f][ip][io] = new int[num_chunks];
      }
      for (int i=0;i<num_chunks;i++) wh[f][ip][Outgoing][i] = 0;
    }

  // Next start setting up the connections...
  
  for (int i=0;i<num_chunks;i++) {
    const volume vi = chunks[i]->v;
    
    // initialize wh[f][ip][Incoming][j] to sum of comm_sizes for jj < j
    FOR_FIELD_TYPES(f)
      for (int ip=0;ip<3;ip++) {
	wh[f][ip][Incoming][0] = 0;
	for (int j = 1; j < num_chunks; ++j)
	  wh[f][ip][Incoming][j] = wh[f][ip][Incoming][j-1]
	    + comm_sizes[f][ip][(j-1)+i*num_chunks];
      }

    FOR_COMPONENTS(corig)
      if (have_component(corig))
	LOOP_OVER_VOL_NOTOWNED(vi, corig, n) {
  	  IVEC_LOOP_ILOC(vi, here);
	  component c = corig;
	  // We're looking at a border element...
	  complex<double> thephase;
	  if (locate_component_point(&c,&here,&thephase)
	      && !on_metal_boundary(here))
	    for (int j=0;j<num_chunks;j++) {
	      if ((chunks[i]->is_mine() || chunks[j]->is_mine())
		  && chunks[j]->v.owns(here)
		  && !(is_B(corig) && is_B(c) &&
		       B_redundant[5*i+corig-Bx] && B_redundant[5*j+c-Bx])) {
		const connect_phase ip = thephase == 1.0 ? CONNECT_COPY 
		  : (thephase == -1.0 ? CONNECT_NEGATE : CONNECT_PHASE);
		const int m = chunks[j]->v.index(c, here);
		const int f = type(c);

		if (ip == CONNECT_PHASE)
		  chunks[i]->connection_phases[f][wh[f][ip][Incoming][j]/2] =
		    thephase;
		DOCMP {
		  chunks[i]->connections[f][ip][Incoming]
		    [wh[f][ip][Incoming][j]++] = chunks[i]->f[corig][cmp] + n;
		  chunks[j]->connections[f][ip][Outgoing]
		    [wh[f][ip][Outgoing][j]++] = chunks[j]->f[c][cmp] + m;
		}
		
		if (is_electric(corig)) {
		  const int f = P_stuff;
		  for (int ipol = 0; ipol < 2; ++ipol) // pol then olpol
		    for (polarization *pi = 
			   ipol ? chunks[i]->olpol : chunks[i]->pol;
			 pi; pi = pi->next)
		      for (polarization *pj = 
			     ipol ? chunks[j]->olpol : chunks[j]->pol;
			   pj; pj = pj->next)
			if (pi->pb->get_identifier()
			    == pj->pb->get_identifier()) {
			  if (ip == CONNECT_PHASE)
			    chunks[i]->connection_phases[f]
			      [wh[f][ip][Incoming][j]/2] = thephase;
			  DOCMP {
			    chunks[i]->connections[f][ip][Incoming]
			      [wh[f][ip][Incoming][j]++] = pi->P[corig][cmp]+n;
			    chunks[j]->connections[f][ip][Outgoing]
			      [wh[f][ip][Outgoing][j]++] = pj->P[c][cmp]+m;
			  }
			}
		} // is_electric(corig)
	      } // if is_mine and owns...
	    } // loop over j chunks
      } // LOOP_OVER_VOL_NOTOWNED
  } // loop over i chunks
  FOR_FIELD_TYPES(f)
    for (int ip=0;ip<3;ip++)
      for (int io=0;io<2;io++)
	delete[] wh[f][ip][io];
  delete[] B_redundant;
}

void fields_chunk::alloc_extra_connections(field_type f, connect_phase ip,
					   in_or_out io, int num) {
  if (num == 0) return; // No need to go to any bother...
  const int tot = num_connections[f][ip][io] + num;
  if (io == Incoming && ip == CONNECT_PHASE) {
    delete[] connection_phases[f];
    connection_phases[f] = new complex<double>[tot];
  }
  typedef double *double_ptr;
  double **conn = new double_ptr[tot];
  if (!conn) abort("Out of memory!\n");
  delete[] connections[f][ip][io];
  connections[f][ip][io] = conn;
  num_connections[f][ip][io] = tot;
}

} // namespace meep