testaverage.cpp

强大的并行科学计算软件包

// -*- C++ -*-
// ACL:license
// ----------------------------------------------------------------------
// This software and ancillary information (herein called "SOFTWARE")
// called POOMA (Parallel Object-Oriented Methods and Applications) is
// made available under the terms described here.  The SOFTWARE has been
// approved for release with associated LA-CC Number LA-CC-98-65.
// 
// Unless otherwise indicated, this SOFTWARE has been authored by an
// employee or employees of the University of California, operator of the
// Los Alamos National Laboratory under Contract No.  W-7405-ENG-36 with
// the U.S. Department of Energy.  The U.S. Government has rights to use,
// reproduce, and distribute this SOFTWARE, and to allow others to do so.
// The public may copy and use this SOFTWARE, FOR NONCOMMERCIAL USE ONLY,
// without charge, provided that this Notice and any statement of
// authorship are reproduced on all copies.  Neither the Government nor
// the University makes any warranty, express or implied, or assumes any
// liability or responsibility for the use of this SOFTWARE.
// 
// If SOFTWARE is modified to produce derivative works, such modified
// SOFTWARE should be clearly marked, so as not to confuse it with the
// version available from LANL.
// 
// For more information about POOMA, send e-mail to [e-mail],
// or visit the POOMA web page at http://www.acl.lanl.gov/pooma/.
// ----------------------------------------------------------------------
// ACL:license

#include "Pooma/Fields.h"

#include <iostream>
using namespace std;

int main(int argc, char *argv[])
{
  Pooma::initialize(argc,argv);

  // Create the physical domains (1D, 2D, and 3D):
  const int nVerts[3] = {4, 4, 4};
  Interval<1> dom1(nVerts[0]);
  Interval<2> dom2(dom1, nVerts[1]);
  Interval<3> dom3(dom2, nVerts[2]);

  int i, j, k, d;
  double epsilon = 1.0e-6;
  double summer = 0.0;
  int nCellsOut, nVertsOut;

  // Flag global passing/failure of all tests:
  bool passed = true;

  // ///////////////////////////////////////////////////////////////////////////
  // 1D

  // Create the mesh:
  typedef RectilinearMesh<1> Mesh1_t;
  Mesh1_t::PointType_t origin1(0.0);
  Vector< 1, Array<1, Mesh1_t::AxisType_t> > spacings1;
  for (d = 0; d < 1; d++) {
    spacings1(d).initialize(dom1[d].size()-1);
    for (i = 0; i < dom1[d].size()-1; i++) {
      spacings1(d)(i) = origin1(d) + (d + 2)*(i + 1);
    }
  }
  Mesh1_t mesh1(dom1, origin1, spacings1);
  
  // Create the geometries:
  DiscreteGeometry<Vert, Mesh1_t> geom1v(mesh1, GuardLayers<1>(1));
  DiscreteGeometry<Cell, Mesh1_t> geom1c(mesh1, GuardLayers<1>(1));
  
  // Make the fields:
  Field<DiscreteGeometry<Vert, Mesh1_t>, Vector<1> > f1v(geom1v);
  Field<DiscreteGeometry<Cell, Mesh1_t>, Vector<1> > f1c(geom1c);
  Field<DiscreteGeometry<Vert, Mesh1_t>, double> w1v(geom1v);
  Field<DiscreteGeometry<Cell, Mesh1_t>, double> w1c(geom1c);
    
  // Set the (constant) Field and weight-Field values:
  f1v(f1v.totalDomain()) = Vector<1>(1.0);
  w1v = 2.0;
  f1c(f1c.totalDomain()) = Vector<1>(1.0);
  w1c = 2.0;
      
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Unweighted Average
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  Interval<1> outCells1 = average<Cell>(f1v).physicalDomain();
  nCellsOut = outCells1.size();
  Vector<1> sum1(0.0);
  for (i = outCells1[0].first(); i <= outCells1[0].last(); i++) {
    sum1 += average<Cell>(f1v)(i);
  }
  sum1 /= nCellsOut;
  summer = dot((sum1 - Vector<1>(1.0)), (sum1 - Vector<1>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Cell>(1D Vert position field) != const=1.0) ; " 
	      << "sum1 = " << sum1 << std::endl;
    passed = false;
  }

  // Cell-->Vert
  Interval<1> outVerts1 = average<Vert>(f1c).physicalDomain();
  nVertsOut = outVerts1.size();
  sum1 = Vector<1>(0.0);
  for (i = outVerts1[0].first(); i <= outVerts1[0].last(); i++) {
    sum1 += average<Vert>(f1c)(i);
  }
  sum1 /= nVertsOut;
  summer = dot((sum1 - Vector<1>(1.0)), (sum1 - Vector<1>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Vert>(1D Cell position field) != const=1.0) ; " 
	      << "sum1 = " << sum1 << std::endl;
    passed = false;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Weighted Average:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  outCells1 = average<Cell>(f1v, w1v).physicalDomain();
  nCellsOut = outCells1.size();
  sum1 = Vector<1>(0.0);
  for (i = outCells1[0].first(); i <= outCells1[0].last(); i++) {
    sum1 += average<Cell>(f1v, w1v)(i);
  }
  sum1 /= nCellsOut;
  summer = dot((sum1 - Vector<1>(1.0)), (sum1 - Vector<1>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Cell>(1D constant field) != const="
	      << 1.0 << ") ; sum1 = " << sum1 << std::endl;
    passed = false;
  }
  
  // Cell-->Vert
  outVerts1 = average<Vert>(f1c, w1c).physicalDomain();
  nVertsOut = outVerts1.size();
  sum1 = Vector<1>(0.0);
  for (i = outVerts1[0].first(); i <= outVerts1[0].last(); i++) {
    sum1 += average<Vert>(f1c, w1c)(i);
  }
  sum1 /= nVertsOut;
  summer = dot((sum1 - Vector<1>(1.0)), (sum1 - Vector<1>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Vert>(1D Cell constant field) != const="
	      << 1.0 << ") ; sum1 = " << sum1 << std::endl;
    passed = false;
  }

  // ///////////////////////////////////////////////////////////////////////////
  // 2D

  // Create the mesh:
  typedef RectilinearMesh<2> Mesh2_t;
  Mesh2_t::PointType_t origin2(0.0);
  Vector< 2, Array<1, Mesh2_t::AxisType_t> > spacings2;
  for (d = 0; d < 2; d++) {
    spacings2(d).initialize(dom2[d].size()-1);
    for (i = 0; i < dom2[d].size()-1; i++) {
      spacings2(d)(i) = origin2(d) + (d + 2)*(i + 1);
    }
  }
  Mesh2_t mesh2(dom2, origin2, spacings2);

  // Create the geometries:
  DiscreteGeometry<Vert, Mesh2_t> geom2v(mesh2, GuardLayers<2>(1));
  DiscreteGeometry<Cell, Mesh2_t> geom2c(mesh2, GuardLayers<2>(1));
  
  // Make the fields:
  Field<DiscreteGeometry<Vert, Mesh2_t>, Vector<2> > f2v(geom2v);
  Field<DiscreteGeometry<Cell, Mesh2_t>, Vector<2> > f2c(geom2c);
  Field<DiscreteGeometry<Vert, Mesh2_t>, double> w2v(geom2v);
  Field<DiscreteGeometry<Cell, Mesh2_t>, double> w2c(geom2c);
    
  // Set the (constant) Field and weight-Field values:
  f2v(f2v.totalDomain()) = Vector<2>(1.0);
  w2v = 2.0;
  f2c(f2c.totalDomain()) = Vector<2>(1.0);
  w2c = 2.0;

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Unweighted Average
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  Interval<2> outCells2 = average<Cell>(f2v).physicalDomain();
  nCellsOut = outCells2.size();
  summer = 0.0;
  Vector<2> sum2(0.0);
  for (i = outCells2[0].first(); i <= outCells2[0].last(); i++) {
    for (j = outCells2[1].first(); j <= outCells2[1].last(); j++) {
      sum2 += average<Cell>(f2v)(i,j);
    }
  }
  sum2 /= nCellsOut;
  summer = dot((sum2 - Vector<2>(1.0)), (sum2 - Vector<2>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Cell>(2D Vert position field) != const=1.0) ; " 
	      << "sum2 = " << sum2 << std::endl;
    passed = false;
  }

  // Cell-->Vert
  Interval<2> outVerts2 = average<Vert>(f2c).physicalDomain();
  nVertsOut = outVerts2.size();
  sum2 = Vector<2>(0.0);
  for (i = outVerts2[0].first(); i <= outVerts2[0].last(); i++) {
    for (j = outVerts2[1].first(); j <= outVerts2[1].last(); j++) {
      sum2 += average<Vert>(f2c)(i,j);
    }
  }
  sum2 /= nVertsOut;
  summer = dot((sum2 - Vector<2>(1.0)), (sum2 - Vector<2>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Vert>(2D Cert position field) != const=1.0) ; " 
	      << "sum2 = " << sum2 << std::endl;
    passed = false;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Weighted Average:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  outCells2 = average<Cell>(f2v, w2v).physicalDomain();
  nCellsOut = outCells2.size();
  sum2 = Vector<2>(0.0);
  for (i = outCells2[0].first(); i <= outCells2[0].last(); i++) {
    for (j = outCells2[1].first(); j <= outCells2[1].last(); j++) {
      sum2 += average<Cell>(f2v, w2v)(i,j);
    }
  }
  sum2 /= nCellsOut;
  summer = dot((sum2 - Vector<2>(1.0)), (sum2 - Vector<2>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Cell>(2D constant field) != const="
	      << 1.0 << ") ; sum2 = " << sum2 << std::endl;
    passed = false;
  }

  // Cell-->Vert
  outVerts2 = average<Vert>(f2c, w2c).physicalDomain();
  nVertsOut = outVerts2.size();
  sum2 = Vector<2>(0.0);
  for (i = outVerts2[0].first(); i <= outVerts2[0].last(); i++) {
    for (j = outVerts2[1].first(); j <= outVerts2[1].last(); j++) {
      sum2 += average<Vert>(f2c, w2c)(i,j);
    }
  }
  sum2 /= nVertsOut;
  summer = dot((sum2 - Vector<2>(1.0)), (sum2 - Vector<2>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Vert>(2D Cell constant field) != const="
	      << 1.0 << ") ; sum2 = " << sum2 << std::endl;
    passed = false;
  }


  // ///////////////////////////////////////////////////////////////////////////
  // 3D

  // Create the mesh
  typedef RectilinearMesh<3> Mesh3_t;
  Mesh3_t::PointType_t origin3(0.0);
  Vector< 3, Array<1, Mesh3_t::AxisType_t> > spacings3;
  for (d = 0; d < 3; d++) {
    spacings3(d).initialize(dom3[d].size()-1);
    for (i = 0; i < dom3[d].size()-1; i++) {
      spacings3(d)(i) = origin3(d) + (d + 2)*(i + 1);
    }
  }
  Mesh3_t mesh3(dom3, origin3, spacings3);
  
  // Create the geometries:
  DiscreteGeometry<Vert, Mesh3_t> geom3v(mesh3, GuardLayers<3>(1));
  DiscreteGeometry<Cell, Mesh3_t> geom3c(mesh3, GuardLayers<3>(1));
  
  // Make the fields:
  Field<DiscreteGeometry<Vert, Mesh3_t>, Vector<3> > f3v(geom3v);
  Field<DiscreteGeometry<Cell, Mesh3_t>, Vector<3> > f3c(geom3c);
  Field<DiscreteGeometry<Vert, Mesh3_t>, double> w3v(geom3v);
  Field<DiscreteGeometry<Cell, Mesh3_t>, double> w3c(geom3c);

  // Set the (constant) Field and weight-Field values:
  f3v(f3v.totalDomain()) = Vector<3>(1.0);
  w3v = 2.0;
  f3c(f3c.totalDomain()) = Vector<3>(1.0);
  w3c = 2.0;
      
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Unweighted Average
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  Interval<3> outCells3 = average<Cell>(f3v).physicalDomain();
  nCellsOut = outCells3.size();
  Vector<3> sum3(0.0);
  for (i = outCells3[0].first(); i <= outCells3[0].last(); i++) {
    for (j = outCells3[1].first(); j <= outCells3[1].last(); j++) {
      for (k = outCells3[2].first(); k <= outCells3[2].last(); k++) {
	sum3 += average<Cell>(f3v)(i,j,k);
      }
    }
  }
  sum3 /= nCellsOut;
  summer = dot((sum3 - Vector<3>(1.0)), (sum3 - Vector<3>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Cell>(3D Vert position field) != const=1.0) ; " 
	      << "sum3 = " << sum3 << std::endl;
    passed = false;
  }

  // Cell-->Vert
  Interval<3> outVerts3 = average<Vert>(f3c).physicalDomain();
  nVertsOut = outVerts3.size();
  sum3 = Vector<3>(0.0);
  for (i = outVerts3[0].first(); i <= outVerts3[0].last(); i++) {
    for (j = outVerts3[1].first(); j <= outVerts3[1].last(); j++) {
      for (k = outVerts3[2].first(); k <= outVerts3[2].last(); k++) {
	sum3 += average<Vert>(f3c)(i,j,k);
      }
    }
  }
  sum3 /= nVertsOut;
  summer = dot((sum3 - Vector<3>(1.0)), (sum3 - Vector<3>(1.0)));
  if (summer > epsilon) {
    std::cout << "average<Vert>(3D Cell position field) != const=1.0) ; " 
	      << "sum3 = " << sum3 << std::endl;
    passed = false;
  }

  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  // Weighted Average:
  // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  // Vert-->Cell
  outCells3 = average<Cell>(f3v, w3v).physicalDomain();
  nCellsOut = outCells3.size();
  sum3 = Vector<3>(0.0);
  for (i = outCells3[0].first(); i <= outCells3[0].last(); i++) {
    for (j = outCells3[1].first(); j <= outCells3[1].last(); j++) {
      for (k = outCells3[2].first(); k <= outCells3[2].last(); k++) {
	sum3 += average<Cell>(f3v, w3v)(i,j,k);
      }
    }
  }
  sum3 /= nCellsOut;
  summer = dot((sum3 - Vector<3>(1.0)), (sum3 - Vector<3>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Cell>(3D constant field) != const="
	      << 1.0 << ") ; sum3 = " << sum3 << std::endl;
    passed = false;
  }

  // Cell-->Vert
  outVerts3 = average<Vert>(f3c, w3c).physicalDomain();
  nVertsOut = outVerts3.size();
  sum3 = Vector<3>(0.0);
  for (i = outVerts3[0].first(); i <= outVerts3[0].last(); i++) {
    for (j = outVerts3[1].first(); j <= outVerts3[1].last(); j++) {
      for (k = outVerts3[2].first(); k <= outVerts3[2].last(); k++) {
	sum3 += average<Vert>(f3c, w3c)(i,j,k);
      }
    }
  }
  sum3 /= nVertsOut;
  summer = dot((sum3 - Vector<3>(1.0)), (sum3 - Vector<3>(1.0)));
  if (summer > epsilon) {
    std::cout << "weighted average<Vert>(3D Cell constant field) != const="
	      << 1.0 << ") ; sum3 = " << sum3 << std::endl;
    passed = false;
  }

  std::cout << ( (passed) ? "PASSED" : "FAILED" ) 
    << " TestAverage" << std::endl;;

  Pooma::finalize();
  return 0;
}

// ACL:rcsinfo
// ----------------------------------------------------------------------
// $RCSfile: TestAverage.cpp,v $   $Author: swhaney $
// $Revision: 1.6 $   $Date: 1999/08/20 22:11:05 $
// ----------------------------------------------------------------------
// ACL:rcsinfo