boundaries.cpp

麻省理工的计算光子晶体的程序

    for (int sn=0;sn<S.multiplicity();sn++) {
      const ivec here=S.transform(*there,sn);
      if (v.owns(here)) {
        *there = here;
        *phase *= S.phase_shift(*c,sn);
        *c = direction_component(*c,
                                 S.transform(component_direction(*c),sn).d);
        return true;
      }
    }
  return false;
}

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 (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)) {
		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 code in boundaries.cpp.  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)) {
		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];
}

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