C++ Array::push_back方法代码示例

//Neumann BCs
void
Albany::PNPProblem::constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
{

   // Note: we only enter this function if sidesets are defined in the mesh file
   // i.e. meshSpecs.ssNames.size() > 0
   
    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;

   // Check to make sure that Neumann BCs are given in the input file

   if(!nbcUtils.haveBCSpecified(this->params)) {
      return;
   }

   // Construct BC evaluators for all side sets and names
   // Note that the string index sets up the equation offset, so ordering is important
   //
   // Currently we aren't exactly doing this right.  I think to do this
   // correctly we need different neumann evaluators for each DOF (velocity,
   // pressure, temperature, flux) since velocity is a vector and the 
   // others are scalars.  The dof_names stuff is only used
   // for robin conditions, so at this point, as long as we don't enable
   // robin conditions, this should work.
   
   std::vector<std::string> nbcNames;
   Teuchos::RCP< Teuchos::Array<std::string> > dof_names =
     Teuchos::rcp(new Teuchos::Array<std::string>);
   // TODO: arbitraty numSpecies
   Teuchos::Array<Teuchos::Array<int> > offsets;
   int idx = 0;
     nbcNames.push_back("C1");
     offsets.push_back(Teuchos::Array<int>(1,idx++));
     if (numSpecies>=2) {
       nbcNames.push_back("C2");
       offsets.push_back(Teuchos::Array<int>(1,idx++));
     }
     if (numSpecies==3) {
       nbcNames.push_back("C3");
       offsets.push_back(Teuchos::Array<int>(1,idx++));
     }
     nbcNames.push_back("Phi");
     offsets.push_back(Teuchos::Array<int>(1,idx++));
     dof_names->push_back("Concentration");
     dof_names->push_back("Potential");

   // Construct BC evaluators for all possible names of conditions
   // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both
   std::vector<std::string> condNames; //dudx, dudy, dudz, dudn, basal 


   nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, nbcNames,
                                           Teuchos::arcp(dof_names),
                                           true, 0, condNames, offsets, dl,
                                           this->params, this->paramLib);
}

void
Tpetra::Utils::generateMatrix(const Teuchos::RCP<Teuchos::ParameterList> &plist,
                              const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
                              const Teuchos::RCP<Node> &node,
                              Teuchos::RCP< Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &A)
{
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    TEUCHOS_TEST_FOR_EXCEPTION( plist == Teuchos::null, std::runtime_error,
                                "Tpetra::Utils::generateMatrix(): ParameterList is null.");
    TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::isParameterType<std::string>(*plist,"mat_type") == false, std::runtime_error,
                                "Tpetra::Utils::generateMatrix(): ParameterList did not contain string parameter ""mat_type"".");
    std::string mat_type = plist->get<std::string>("mat_type");
    if (mat_type == "Lap3D") {
        // 3D Laplacian, grid is a cube with dimension gridSize x gridSize x gridSize
        const GlobalOrdinal gridSize = as<GlobalOrdinal>(plist->get<int>("gridSize",100));
        const GlobalOrdinal gS2 = gridSize*gridSize;
        const GlobalOrdinal numRows = gS2*gridSize;
        Teuchos::RCP<Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap;
        rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(as<global_size_t>(numRows),as<GlobalOrdinal>(0),comm,GloballyDistributed,node));
        A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,7,Tpetra::StaticProfile));
        // fill matrix, one row at a time
        Teuchos::Array<GlobalOrdinal> neighbors;
        Teuchos::Array<Scalar> values(7, -ST::one());
        values[0] = (Scalar)6;
        for (GlobalOrdinal r = rowMap->getMinGlobalIndex(); r <= rowMap->getMaxGlobalIndex(); ++r) {
            neighbors.clear();
            neighbors.push_back(r); // add diagonal
            GlobalOrdinal ixy, iz, ix, iy;  // (x,y,z) coords and index in xy plane
            ixy = r%gS2;
            iz = (r - ixy)/gS2;
            ix = ixy%gridSize;
            iy = (ixy - ix)/gridSize;
            //
            if ( ix != 0 )          neighbors.push_back( r-1 );
            if ( ix != gridSize-1 ) neighbors.push_back( r+1 );
            if ( iy != 0 )          neighbors.push_back( r-gridSize );
            if ( iy != gridSize-1 ) neighbors.push_back( r+gridSize );
            if ( iz != 0 )          neighbors.push_back( r-gS2 );
            if ( iz != gridSize-1 ) neighbors.push_back( r+gS2 );
            A->insertGlobalValues( r, neighbors(), values(0,neighbors.size()) );
        }
        A->fillComplete();
    }
    else {
        TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
                                    "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
    }
}

Teuchos::Array<int> getMyBlockLIDs(
    const Teuchos::ParameterList &blockingParams,
    const Albany::AbstractDiscretization &disc)
{
  Teuchos::Array<int> result;

  const std::string nodeSetLabel = blockingParams.get<std::string>("Node Set");
  const int dofRank = blockingParams.get<int>("Dof");

  const Albany::NodeSetList &nodeSets = disc.getNodeSets();
  const Albany::NodeSetList::const_iterator it = nodeSets.find(nodeSetLabel);
  TEUCHOS_TEST_FOR_EXCEPT(it == nodeSets.end());
  {
    typedef Albany::NodeSetList::mapped_type NodeSetEntryList;
    const NodeSetEntryList &nodeEntries = it->second;

    for (NodeSetEntryList::const_iterator jt = nodeEntries.begin(); jt != nodeEntries.end(); ++jt) {
      typedef NodeSetEntryList::value_type NodeEntryList;
      const NodeEntryList &entries = *jt;
      result.push_back(entries[dofRank]);
    }
  }

  return result;
}

Teuchos::Array<bool> createResponseTable(
    int count,
    const std::string selectionType,
    int index,
    const Teuchos::ArrayView<const int> &list)
{
  Teuchos::Array<bool> result;

  if (count > 0) {
    if (selectionType == "All") {
      result.resize(count, true);
    } else if (selectionType == "Last") {
      result = createResponseTableFromIndex(count - 1, count);
    } else if (selectionType == "AllButLast") {
      result.reserve(count);
      result.resize(count - 1, true);
      result.push_back(false);
    } else if (selectionType == "Index") {
      result = createResponseTableFromIndex(index, count);
    } else if (selectionType == "List") {
      result.resize(count, false);
      for (Teuchos::ArrayView<const int>::const_iterator it = list.begin(), it_end = list.end(); it != it_end; ++it) {
        result.at(*it) = true;
      }
    } else {
      TEUCHOS_TEST_FOR_EXCEPT(false);
    }
  }

  return result;
}

Teuchos::Array<GO>
Albany::NodeGIDsSolutionCullingStrategy::
selectedGIDsT(Teuchos::RCP<const Tpetra_Map> sourceMapT) const
{
  Teuchos::Array<GO> result;
  {
    Teuchos::Array<GO> mySelectedGIDs;

    // Subract 1 to convert exodus GIDs to our GIDs
    for (int i=0; i<nodeGIDs_.size(); i++)
      if (sourceMapT->isNodeGlobalElement(nodeGIDs_[i] -1) ) mySelectedGIDs.push_back(nodeGIDs_[i] - 1);

    Teuchos::RCP<const Teuchos::Comm<int> >commT = sourceMapT->getComm(); 
    {
      GO selectedGIDCount;
      {
        GO mySelectedGIDCount = mySelectedGIDs.size();
        Teuchos::reduceAll<LO, GO>(*commT, Teuchos::REDUCE_SUM, 1, &mySelectedGIDCount, &selectedGIDCount); 
      }
      result.resize(selectedGIDCount);
    }

    const int ierr = Tpetra::GatherAllV(
        commT,
        mySelectedGIDs.getRawPtr(), mySelectedGIDs.size(),
        result.getRawPtr(), result.size());
    TEUCHOS_ASSERT(ierr == 0);
  }

  std::sort(result.begin(), result.end());

  return result;
}

  Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
  createMeshMap (const LO& blockSize, const Tpetra::Map<LO,GO,Node>& pointMap)
  {
    typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
    typedef Tpetra::Map<LO,GO,Node> map_type;

    //calculate mesh GIDs
    Teuchos::ArrayView<const GO> pointGids = pointMap.getNodeElementList();
    Teuchos::Array<GO> meshGids;
    GO indexBase = pointMap.getIndexBase();

    // Use hash table to track whether we've encountered this GID previously.  This will happen
    // when striding through the point DOFs in a block.  It should not happen otherwise.
    // I don't use sort/make unique because I don't want to change the ordering.
    meshGids.reserve(pointGids.size());
    Tpetra::Details::HashTable<GO,int> hashTable(pointGids.size());
    for (int i=0; i<pointGids.size(); ++i) {
      GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
      if (hashTable.get(meshGid) == -1) {
        hashTable.add(meshGid,1);   //(key,value)
        meshGids.push_back(meshGid);
      }
    }

    Teuchos::RCP<const map_type> meshMap = Teuchos::rcp( new map_type(TOT::invalid(), meshGids(), 0, pointMap.getComm()) );
    return meshMap;

  }

// ============================================================================
Teuchos::RCP<Recti::Domain::Abstract>
Recti::Domain::Factory::
buildPolygon() const
{
    const Teuchos::ParameterList & pointsList =
            pList_.sublist("Vertices");
    
    Teuchos::Array<Point> vertices;
    
    Teuchos::ParameterList::ConstIterator k;
    for ( k=pointsList.begin(); k!=pointsList.end(); ++k )
    {
        Teuchos::ParameterEntry e = pointsList.entry(k);
        TEUCHOS_ASSERT( e.isList() );
        
        std::string label = pointsList.name(k);
        const Teuchos::ParameterList & coordinates =
                pointsList.sublist(label);
        
        Point X = Teuchos::tuple( coordinates.get<double>("x"),
                                  coordinates.get<double>("y"),
                                  0.0
                                );
        vertices.push_back( X );
    }
    
    return Teuchos::rcp( new Polygon( vertices ) );
}

Teuchos::Array<std::string>
SolverFactory<Scalar, MV, OP>::supportedSolverNames () const
{
  typedef std::vector<std::string>::const_iterator iter_type;
  Teuchos::Array<std::string> names;

  {
    std::vector<std::string> aliases = Details::solverNameAliases ();
    for (iter_type iter = aliases.begin (); iter != aliases.end (); ++iter) {
      names.push_back (*iter);
    }
  }
  {
    std::vector<std::string> canonicalNames = Details::canonicalSolverNames ();
    for (iter_type iter = canonicalNames.begin (); iter != canonicalNames.end (); ++iter) {
      names.push_back (*iter);
    }
  }
  return names;
}

void EricksonManufacturedSolution::getConstraints(FieldContainer<double> &physicalPoints, 
						  FieldContainer<double> &unitNormals,
						  vector<map<int,FieldContainer<double > > > &constraintCoeffs,
						  vector<FieldContainer<double > > &constraintValues){
    
  int numCells = physicalPoints.dimension(0);
  int numPoints = physicalPoints.dimension(1);
  int spaceDim = physicalPoints.dimension(2);       
  map<int,FieldContainer<double> > outflowConstraint;
  FieldContainer<double> uCoeffs(numCells,numPoints);
  FieldContainer<double> beta_sigmaCoeffs(numCells,numPoints);
  FieldContainer<double> outflowValues(numCells,numPoints);
  double tol = 1e-12;
  // default to no constraints, apply on outflow only
  uCoeffs.initialize(0.0);
  beta_sigmaCoeffs.initialize(0.0);
  outflowValues.initialize(0.0);

  Teuchos::Array<int> pointDimensions;
  pointDimensions.push_back(2);    
  for (int cellIndex=0;cellIndex<numCells;cellIndex++){
    for (int pointIndex=0;pointIndex<numPoints;pointIndex++){
      FieldContainer<double> physicalPoint(pointDimensions,
                                           &physicalPoints(cellIndex,pointIndex,0));
      FieldContainer<double> unitNormal(pointDimensions,
                                        &unitNormals(cellIndex,pointIndex,0));

      double x = physicalPoints(cellIndex,pointIndex,0);
      double y = physicalPoints(cellIndex,pointIndex,1);

      double beta_n = _beta_x*unitNormals(cellIndex,pointIndex,0)+_beta_y*unitNormals(cellIndex,pointIndex,1);
      
      if ( abs(x-1.0) < tol ) { // if on outflow boundary
	TEUCHOS_TEST_FOR_EXCEPTION(beta_n < 0,std::invalid_argument,"Inflow condition on boundary");
	
	// this combo isolates sigma_n
	//uCoeffs(cellIndex,pointIndex) = 1.0;
	uCoeffs(cellIndex,pointIndex) = beta_n;
	beta_sigmaCoeffs(cellIndex,pointIndex) = -1.0;	    
	double beta_n_u_minus_sigma_n = solutionValue(ConfusionBilinearForm::BETA_N_U_MINUS_SIGMA_HAT, physicalPoint, unitNormal);
	double u_hat = solutionValue(ConfusionBilinearForm::U_HAT, physicalPoint, unitNormal);
	outflowValues(cellIndex,pointIndex) = beta_n*u_hat - beta_n_u_minus_sigma_n; // sigma_n
      }	
    }
  }
  outflowConstraint[ConfusionBilinearForm::U_HAT] = uCoeffs;
  outflowConstraint[ConfusionBilinearForm::BETA_N_U_MINUS_SIGMA_HAT] = beta_sigmaCoeffs;	        
  if (!_useWallBC){
    constraintCoeffs.push_back(outflowConstraint); // only one constraint on outflow
    constraintValues.push_back(outflowValues); // only one constraint on outflow
  }
}

  /*! \brief Get the parts belonging to this rank
   *  \param numParts on return, set to the number of parts assigned to rank.
   *  \param parts on return, pointer (view) to the parts assigned to rank
   */
  void getMyPartsView(part_t &numParts, part_t *&parts)
  {
    bool useAlg = true;

    // Check first whether this algorithm answers getMyPartsView.
    if (mapping_algorithm != Teuchos::null) {
      try {
        mapping_algorithm->getMyPartsView(numParts, parts);
      }
      catch (NotImplemented &e) {
        // It is OK if the algorithm did not implement this method;
        // we'll get the information from the solution below.
        useAlg = false;
      }
      Z2_FORWARD_EXCEPTIONS;
    }

    if (!useAlg) {  

      // Algorithm did not implement this method.

      // Did the algorithm register a result with the solution?
      if ((partsForRank==Teuchos::null) && (rankForPart==Teuchos::null)) {
        numParts = 0;
        parts = NULL;
        throw std::runtime_error("No mapping solution available.");
      }
  
      if (partsForRank == Teuchos::null) {
        // Need to create the array since we haven't created it before.
        Teuchos::Array<part_t> tmp;

        part_t cnt = 0;
        for (typename rankmap_t::iterator it = rankForPart->begin();
             it != rankForPart->end(); it++) {
          if (it->second == myRank) {
            tmp.push_back(it->first); 
            cnt++;
          }
        }
        if (cnt)
          partsForRank = arcp(&tmp[0], 0, cnt, true);
      }

      numParts = partsForRank.size();
      if (numParts)
        parts = partsForRank.getRawPtr();
      else 
        parts = NULL;
    }
  }

void EricksonManufacturedSolution::imposeBC(int varID, FieldContainer<double> &physicalPoints,
					    FieldContainer<double> &unitNormals,
					    FieldContainer<double> &dirichletValues,
					    FieldContainer<bool> &imposeHere) {
  int numCells = physicalPoints.dimension(0);
  int numPoints = physicalPoints.dimension(1);
  int spaceDim = physicalPoints.dimension(2);
  
  TEUCHOS_TEST_FOR_EXCEPTION( ( spaceDim != 2  ),
		      std::invalid_argument,
		      "ConfusionBC expects spaceDim==2.");  
  
  TEUCHOS_TEST_FOR_EXCEPTION( ( dirichletValues.dimension(0) != numCells ) 
		      || ( dirichletValues.dimension(1) != numPoints ) 
		      || ( dirichletValues.rank() != 2  ),
		      std::invalid_argument,
		      "dirichletValues dimensions should be (numCells,numPoints).");
  TEUCHOS_TEST_FOR_EXCEPTION( ( imposeHere.dimension(0) != numCells ) 
		      || ( imposeHere.dimension(1) != numPoints ) 
		      || ( imposeHere.rank() != 2  ),
		      std::invalid_argument,
		      "imposeHere dimensions should be (numCells,numPoints).");
  Teuchos::Array<int> pointDimensions;
  pointDimensions.push_back(2);
  for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
      FieldContainer<double> physicalPoint(pointDimensions,
                                           &physicalPoints(cellIndex,ptIndex,0));
      FieldContainer<double> unitNormal(pointDimensions,
                                        &unitNormals(cellIndex,ptIndex,0));
      double beta_n = unitNormals(cellIndex,ptIndex,0)*_beta_x + unitNormals(cellIndex,ptIndex,1)*_beta_y;
      double x = physicalPoint(cellIndex,ptIndex,0);
      double y = physicalPoint(cellIndex,ptIndex,1);
      imposeHere(cellIndex,ptIndex) = false;
      //      if (varID==ConfusionBilinearForm::BETA_N_U_MINUS_SIGMA_HAT) {
      //	dirichletValues(cellIndex,ptIndex) = solutionValue(varID, physicalPoint, unitNormal);
      if (varID==ConfusionBilinearForm::U_HAT) {
	dirichletValues(cellIndex,ptIndex) = solutionValue(varID, physicalPoint);
	if ( abs(x-1.0) > 1e-12) { // if not the outflow (pts on boundary already)
	  imposeHere(cellIndex,ptIndex) = true;
	}

	// wall boundary 
	if (abs(x-1.0)<1e-12 && _useWallBC){
	  dirichletValues(cellIndex,ptIndex) = solutionValue(varID, physicalPoint);
	  imposeHere(cellIndex,ptIndex) = true;
	}
      } 
    }
  }
}

void epetraFromThyra(
  const Teuchos::RCP<const Epetra_Comm> &comm,
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(epetraVectorFromThyra(comm, *it));
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(epetraMultiVectorFromThyra(comm, *jt));
    }
    sensitivities.push_back(sens);
  }
}

Teuchos::Array< int >
splitStringOfIntsWithCommas(std::string data)
{
  Teuchos::Array< int > result;
  size_t current = 0;
  while (current < data.size())
  {
    size_t next = data.find(",", current);
    if (next == std::string::npos) next = data.size();
    result.push_back(atoi(data.substr(current, next-current).c_str()));
    current = next + 1;
  }
  return result;
}

void tpetraFromThyra(
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Tpetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(Teuchos::nonnull(*it) ? ConverterT::getConstTpetraVector(*it) : Teuchos::null);
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
  ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(Teuchos::nonnull(*jt) ? ConverterT::getConstTpetraMultiVector(*jt) : Teuchos::null);
    }
    sensitivities.push_back(sens);
  }
}

void ExactSolution::solutionValues(FieldContainer<double> &values, int trialID,
                                   FieldContainer<double> &physicalPoints) {
  int numCells = physicalPoints.dimension(0);
  int numPoints = physicalPoints.dimension(1);
  int spaceDim = physicalPoints.dimension(2);
  
  Teuchos::Array<int> pointDimensions;
  pointDimensions.push_back(spaceDim);
  
//  cout << "ExactSolution: physicalPoints:\n" << physicalPoints;
  
  for (int cellIndex=0; cellIndex<numCells; cellIndex++) {
    for (int ptIndex=0; ptIndex<numPoints; ptIndex++) {
      FieldContainer<double> point(pointDimensions,&physicalPoints(cellIndex,ptIndex,0));
      double value = solutionValue(trialID, point);
      values(cellIndex,ptIndex) = value;
    }
  }
}