C++ Array1类代码示例

void SphSystemData2::computeMass() {
    Array1<Vector2D> points;
    TrianglePointGenerator pointsGenerator;
    BoundingBox2D sampleBound(
        Vector2D(-1.5*_kernelRadius, -1.5*_kernelRadius),
        Vector2D(1.5*_kernelRadius, 1.5*_kernelRadius));

    pointsGenerator.generate(sampleBound, _targetSpacing, &points);

    double maxNumberDensity = 0.0;
    SphStdKernel2 kernel(_kernelRadius);

    for (size_t i = 0; i < points.size(); ++i) {
        const Vector2D& point = points[i];
        double sum = 0.0;

        for (size_t j = 0; j < points.size(); ++j) {
            const Vector2D& neighborPoint = points[j];
            sum += kernel(neighborPoint.distanceTo(point));
        }

        maxNumberDensity = std::max(maxNumberDensity, sum);
    }

    JET_ASSERT(maxNumberDensity > 0);

    double newMass = _targetDensity / maxNumberDensity;

    ParticleSystemData2::setMass(newMass);
}

TEST(PointSimpleListSearcher3, ForEachNearbyPoint) {
    Array1<Vector3D> points = {
        Vector3D(0, 1, 3),
        Vector3D(2, 5, 4),
        Vector3D(-1, 3, 0)
    };

    PointSimpleListSearcher3 searcher;
    searcher.build(points.accessor());

    int cnt = 0;
    searcher.forEachNearbyPoint(
        Vector3D(0, 0, 0),
        std::sqrt(10.0),
        [&](size_t i, const Vector3D& pt) {
            EXPECT_TRUE(i == 0 || i == 2);

            if (i == 0) {
                EXPECT_EQ(points[0], pt);
            } else if (i == 2) {
                EXPECT_EQ(points[2], pt);
            }

            ++cnt;
        });

    EXPECT_EQ(2, cnt);
}

TEST(PointParallelHashGridSearcher3, CopyConstructor) {
    Array1<Vector3D> points = {
        Vector3D(0, 1, 3),
        Vector3D(2, 5, 4),
        Vector3D(-1, 3, 0)
    };

    PointParallelHashGridSearcher3 searcher(4, 4, 4, std::sqrt(10));
    searcher.build(points.accessor());

    PointParallelHashGridSearcher3 searcher2(searcher);
    int cnt = 0;
    searcher2.forEachNearbyPoint(
        Vector3D(0, 0, 0),
        std::sqrt(10.0),
        [&](size_t i, const Vector3D& pt) {
            EXPECT_TRUE(i == 0 || i == 2);

            if (i == 0) {
                EXPECT_EQ(points[0], pt);
            } else if (i == 2) {
                EXPECT_EQ(points[2], pt);
            }

            ++cnt;
        });
    EXPECT_EQ(2, cnt);
}

void Force_fitter::fit_force(Array1 <doublevar> & r, 
                             Array1 <doublevar> & bare_force,
                             Array1 <doublevar> & fit_force) {
  int ndim=bare_force.GetDim(0);

  fit_force=bare_force; return;

  fit_force.Resize(ndim);
  fit_force=0;

  if(r(0) < cut) {
    for(int d=0; d< ndim; d++) {
      doublevar basis=0;
      for(int i=0; i< nexp; i++) {
        basis+=coeff(i)*pow(r(0),i+m);
      }
      fit_force(d)=bare_force(d)*basis;
    }
  }
  else {
    fit_force=bare_force;
  }


}

JET_END_TEST_F

JET_BEGIN_TEST_F(FlipSolver2, Rotation) {
    // Build solver
    auto solver = FlipSolver2::builder()
        .withResolution({10, 10})
        .withDomainSizeX(1.0)
        .makeShared();

    solver->setGravity({0, 0});
    solver->setPressureSolver(nullptr);

    // Build emitter
    auto box = Sphere2::builder()
        .withCenter({0.5, 0.5})
        .withRadius(0.4)
        .makeShared();

    auto emitter = VolumeParticleEmitter2::builder()
        .withSurface(box)
        .withSpacing(1.0 / 20.0)
        .withIsOneShot(true)
        .makeShared();

    solver->setParticleEmitter(emitter);

    Array1<double> r;

    for (Frame frame; frame.index < 360; ++frame) {
        auto x = solver->particleSystemData()->positions();
        auto v = solver->particleSystemData()->velocities();
        r.resize(x.size());
        for (size_t i = 0; i < x.size(); ++i) {
            r[i] = (x[i] - Vector2D(0.5, 0.5)).length();
        }

        solver->update(frame);

        if (frame.index == 0) {
            x = solver->particleSystemData()->positions();
            v = solver->particleSystemData()->velocities();
            for (size_t i = 0; i < x.size(); ++i) {
                Vector2D rp = x[i] - Vector2D(0.5, 0.5);
                v[i].x = rp.y;
                v[i].y = -rp.x;
            }
        } else {
            for (size_t i = 0; i < x.size(); ++i) {
                Vector2D rp = x[i] - Vector2D(0.5, 0.5);
                if (rp.lengthSquared() > 0.0) {
                    double scale = r[i] / rp.length();
                    x[i] = scale * rp + Vector2D(0.5, 0.5);
                }
            }
        }

        saveParticleDataXy(solver->particleSystemData(), frame.index);
    }
}

//----------------------------------------------------------------------
doublevar dot(const Array1 <doublevar> & a,  const Array1 <doublevar> & b){
  int dim=a.GetSize();
  assert(a.GetSize()==b.GetSize());
  doublevar sum=0;
  for(int i=0;i<dim;i++)
    sum+=a(i)*b(i);
  return sum;
}

void copy_array (Array1& source, Array2& dest) {
  assert(std::equal(source.shape(),source.shape()+source.num_dimensions(),
                    dest.shape()));
  // Dispatch to the proper function
  typedef typename Array1::element element_type;
  copy_dispatch<element_type>::
    copy_array(source.begin(),source.end(),dest.begin());
}

TEST(PointKdTreeSearcher3, ForEachNearbyPointEmpty) {
    Array1<Vector3D> points;

    PointKdTreeSearcher3 searcher;
    searcher.build(points.accessor());

    searcher.forEachNearbyPoint(Vector3D(0, 0, 0), std::sqrt(10.0),
                                [](size_t, const Vector3D&) {});
}

log_value<dcomplex>
TransposeInverseMatrix(const Array2 < complex <doublevar> > & a, 
                       Array2 < complex <doublevar> > & a1, 
                       const int n)
{
  Array2 <complex <doublevar> >  temp(n,n);
  Array1 <int> indx(n);
  doublevar d;

  // a(i,j) first index i is row index (convention)
  // elements of column vectors are stored contiguous in memory in C style arrays
  // a(i) refers to a column vector

  // calculate the inverse of the transposed matrix because this
  // allows to pass a column vector to lubksb() instead of a row

  // put the transposed matrix in temp
  //cout << "temp " << endl;
  for(int i=0;i<n;++i)
  {
    for(int j=0;j<n;++j)
    {
      temp(i,j)=a(i,j);
      a1(i,j)=complex <doublevar> (0.0,0.0);
    }
    a1(i,i)=complex <doublevar> (1.0,0.0);
  }

  //cout << "ludcmp" << endl;
  //if the matrix is singular, the determinant is zero.
  if(ludcmp(temp,n,indx,d)==0) { 
#ifdef SUPERDEBUG
    cout << "log_value<dcomplex>TransposeInverseMatrix:zero determinant " << endl;
#endif
    return dcomplex(0.0,0.0);
  }

  //cout << "lubksb" << endl;

  for(int j=0;j<n;++j)
  {
    // get column vector
    Array1 <complex <doublevar> > yy;//(a1(j));
    yy.refer(a1(j));
    lubksb(temp,n,indx,yy);
  }


  //complex <doublevar> det(d,0);
  log_value<dcomplex> det=dcomplex(d,0);
  for(int j=0;j<n;++j) {
    det *= temp(j,j);
  }
  return det;
}

    void exportStates(Array1<double>& x, Array1<double>& y) const
    {
        x.resize(positions.size());
        y.resize(positions.size());

        for (size_t i = 0; i < positions.size(); ++i)
        {
            x[i] = positions[i].x;
            y[i] = positions[i].y;
        }
    }

void Properties_final_average::showSummary(ostream & os, Array1 <Average_generator*> avg_gen) { 
  int nwf=avg.GetDim(1);
  if(avg_gen.GetDim(0)%nwf != 0 ) error("Something wrong in the showSummary translation");
  int navg=avg_gen.GetDim(0)/nwf;
  Array2 <Average_generator*> avg2(nwf,navg);
  int count=0;
  for(int w=0; w< nwf; w++) {
    for(int i=0; i< navg; i++) avg2(w,i)=avg_gen(count++);
  }
  showSummary(os, avg2);
}

int HEG_system::enforcePbc(Array1 <doublevar> & pos, Array1 <int> & nshifted) {
  assert(pos.GetDim(0) >=3);
  int shifted=0;
  nshifted.Resize(3);
  nshifted=0;

  int nshift=0;
  for(int i=0; i< 3; i++) {
    int shouldcontinue=1;
    while(shouldcontinue) {

      //Whether we're past the origin side
      doublevar tooshort=0;
      for(int j=0; j<3;j++) tooshort+=normVec(i,j)*(pos(j)-origin(j));

      //Whether we're past the lattice vector
      doublevar toofar=0;
      for(int j=0; j< 3; j++) toofar+=normVec(i,j)*(pos(j)-corners(i,j));

      //the 1e-12 seems to help avoid numerical problems, esp 
      //when integrating over a grid(which tends to hit the edges)
      if(tooshort < -1e-12) {
	//cout <<"tooshort " << tooshort << endl;
        for(int j=0; j< 3; j++) pos(j)+=latVec(i,j);
        shifted=1;
        nshifted(i)+=1;
        nshift++;
      }
      else if(toofar > 1e-12) {
	//cout << "toofar " << toofar << endl;
        for(int j=0; j< 3; j++) pos(j)-=latVec(i,j);
        shifted=1;
	// JK: here was +1, which works fine for real k-points that
	// correspond to standing waves. For general (complex) k-points the
	// wavefunction is "directional", however.
        nshifted(i)-=1;

        nshift++;
      }
      else {
        shouldcontinue=0;
      }
      if(nshift > 1000)
	
        error("Did over 1000 shifts and we're still out of the simulation cell."
	            "  There's probably something wrong.  Position : ", pos(i));
    }
  }

  return shifted;

}

void Molecular_system::setIonPos(int ion, Array1 <doublevar> & r)
{
  assert(r.GetDim(0) >=3);
  Array1 <doublevar> temp(r.GetDim(0));
  temp=r;
  if(use_bounding_box) {
    bounding_box.enforcePbc(temp);
  }


  for(int i=0; i< 3; i++) {

    ions.r(i,ion)=temp(i);
  }
}

void Center_set::assignBasis(Array1 <CBasis_function *> basisfunc)
{
  int totbasis=basisfunc.GetDim(0);
  basis.Resize(ncenters, totbasis);

  for(int i=0; i< ncenters; i++)
  {
    int found=0;
    for(int j=0; j < totbasis; j++)
    {
      if(basisfunc(j)->label() == labels[i])
      {
        appendbasis(j,i);
        found=1;
        //break;
      }

    }
    if(!found)
    {
      cout << "\n****WARNING*****\n"
      << "Couldn't find basis for atom " << labels[i]
      << endl;
    }
  }
}

//----------------------------------------------------------------------
//Keeping this separate from calcLoc for efficiency reasons..
//It's most likely faster to add to a double than fill matrices..
//We also don't need to calculate the ion-ion interaction for this
void Molecular_system::separatedLocal(Sample_point * sample,Array1 <doublevar> & onebody,
    Array2 <doublevar> & twobody)
{
  int nions=sample->ionSize();
  int nelectrons=sample->electronSize();
  onebody.Resize(nelectrons);
  twobody.Resize(nelectrons,nelectrons);
  onebody=0.0; twobody=0.0;
  Array1 <doublevar> R(5);
  sample->updateEIDist();
  sample->updateEEDist();

  for(int e=0; e< nelectrons; e++) {
    for(int i=0; i < nions; i++) {
      sample->getEIDist(e,i, R);
      onebody(e)-=sample->getIonCharge(i)/R(0);
    }
  }
  Array1 <doublevar> R2(5);
  for(int i=0; i< nelectrons; i++) {
    for(int j=i+1; j<nelectrons; j++) {
      sample->getEEDist(i,j,R2);
      twobody(i,j)+= 1/R2(0);
    }
  }

  Array1 <doublevar> pos(3);
  for(int e=0; e< nelectrons; e++) {
    sample->getElectronPos(e,pos);
    for(int d=0; d< 3; d++) 
      onebody(e)-=electric_field(d)*pos(d);
  }
}