C++ RCP::ExtractView方法代码示例

void
PeridigmNS::InterfaceData::WriteExodusOutput(int timeStep, const float & timeValue, Teuchos::RCP<Epetra_Vector> x, Teuchos::RCP<Epetra_Vector> y){

  int error_int = 0;

  int CPU_word_size = 0;
  int IO_word_size = 0;
  float version = 0;
  std::string outputFileNameStr = filename.str();
  std::vector<char> writable(outputFileNameStr.size() + 1);
  std::copy(outputFileNameStr.begin(), outputFileNameStr.end(), writable.begin());

  exoid = ex_open(&writable[0], EX_WRITE, &CPU_word_size, &IO_word_size, &version);

  error_int = ex_put_time(exoid, timeStep, &timeValue);
  TEUCHOS_TEST_FOR_EXCEPTION(error_int,std::logic_error, "ex_put_time(): Failure");

  float * quadValues = new float[numQuads];
  float * triValues = new float[numTris];

  // populate the quad values
  int quadIndex = 0;
  int triIndex = 0;
  for(int i=0;i<numOwnedPoints;++i){
    if(interfaceNodesMap->ElementSize(i)==4){
      quadValues[quadIndex] = (*interfaceAperture)[i];
      quadIndex++;
    }
    else if(interfaceNodesMap->ElementSize(i)==3){
      triValues[triIndex] = (*interfaceAperture)[i];
      triIndex++;
    }
    else{
      TEUCHOS_TEST_FOR_EXCEPTION(true,std::invalid_argument,"size of this element is not recognized: " << interfaceNodesMap->ElementSize(i));
    }
  }

  int blockIndex = 0;
  const int varIndex = 1;
  blockIndex++;
  if(numQuads > 0){
    error_int = ex_put_elem_var(exoid, timeStep, varIndex, blockIndex, numQuads, &quadValues[0]);
    TEUCHOS_TEST_FOR_EXCEPTION(error_int,std::logic_error,"Failure ex_put_elem_var(): ");
  }
  blockIndex++;
  if(numTris > 0){
    error_int = ex_put_elem_var(exoid, timeStep, varIndex, blockIndex, numTris, &triValues[0]);
    TEUCHOS_TEST_FOR_EXCEPTION(error_int,std::logic_error,"Failure ex_put_elem_var(): ");
  }

  delete [] quadValues;
  delete [] triValues;

  // update the apertures...
  // import the mothership vectors x and y to the overlap epetra vectors
  Teuchos::RCP<const Epetra_Import> importer = Teuchos::rcp(new Epetra_Import(*elemOverlapMap, x->Map()));

  Teuchos::RCP<Epetra_Vector> xOverlap = Teuchos::rcp(new Epetra_Vector(*elemOverlapMap,true));
  xOverlap->Import(*x,*importer,Insert);
  Teuchos::RCP<Epetra_Vector> yOverlap = Teuchos::rcp(new Epetra_Vector(*elemOverlapMap,true));
  yOverlap->Import(*y,*importer,Insert);

  double *xValues;
  xOverlap->ExtractView( &xValues );
  double *yValues;
  yOverlap->ExtractView( &yValues );

  double xLeft=0,yLeft=0,zLeft=0,xRight=0,yRight=0,zRight=0;
  double XLeft=0,YLeft=0,ZLeft=0,XRight=0,YRight=0,ZRight=0;
  double X=0,Y=0;
  double dx=0,dy=0,dz=0,dX=0,dY=0,dZ=0;
  int elemIndexLeft=-1,elemIndexRight=-1,GIDLeft=-1,GIDRight=-1;

  for(int i=0;i<numOwnedPoints;++i){
    GIDLeft = elementLeft[i];
    GIDRight = elementRight[i];

    elemIndexLeft = xOverlap->Map().FirstPointInElement(elemOverlapMap->LID(GIDLeft));
    elemIndexRight = xOverlap->Map().FirstPointInElement(elemOverlapMap->LID(GIDRight));

    xLeft = xValues[elemIndexLeft+0];
    yLeft = xValues[elemIndexLeft+1];
    zLeft = xValues[elemIndexLeft+2];
    xRight = xValues[elemIndexRight+0];
    yRight = xValues[elemIndexRight+1];
    zRight = xValues[elemIndexRight+2];

    XLeft = yValues[elemIndexLeft+0];
    YLeft = yValues[elemIndexLeft+1];
    ZLeft = yValues[elemIndexLeft+2];
    XRight = yValues[elemIndexRight+0];
    YRight = yValues[elemIndexRight+1];
    ZRight = yValues[elemIndexRight+2];

    dx = xRight - xLeft;
    dy = yRight - yLeft;
    dz = zRight - zLeft;

    dX = XRight - XLeft;
    dY = YRight - YLeft;
//.........这里部分代码省略.........

int Ifpack_SORa::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
  if(!IsComputed_) return -1;
  Time_.ResetStartTime();
  bool initial_guess_is_zero=false;
  const int lclNumRows = W_->NumMyRows();
  const int NumVectors = X.NumVectors();
  Epetra_MultiVector Temp(A_->RowMatrixRowMap(),NumVectors);

  double omega=GetOmega();

  // need to create an auxiliary vector, Xcopy
  Teuchos::RCP<const Epetra_MultiVector> Xcopy;
  if (X.Pointers()[0] == Y.Pointers()[0]){
    Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
    // Since the user didn't give us anything better, our initial guess is zero.
    Y.Scale(0.0);
    initial_guess_is_zero=true;
  }
  else
    Xcopy = Teuchos::rcp( &X, false );

  Teuchos::RCP< Epetra_MultiVector > T2;
  // Note: Assuming that the matrix has an importer.  Ifpack_PointRelaxation doesn't do this, but given that
  // I have a CrsMatrix, I'm probably OK.
  // Note: This is the lazy man's version sacrificing a few extra flops for avoiding if statements to determine
  // if things are on or off processor.
  // Note: T2 must be zero'd out
  if (IsParallel_ && W_->Importer())  T2 = Teuchos::rcp( new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors,true));
  else T2 = Teuchos::rcp( new Epetra_MultiVector(A_->RowMatrixRowMap(),NumVectors,true));

  // Pointer grabs
  int* rowptr,*colind;
  double *values;
  double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
  T2->ExtractView(&t2_ptr);
  Y.ExtractView(&y_ptr);
  Temp.ExtractView(&t_ptr);
  Xcopy->ExtractView(&x_ptr);
  Wdiag_->ExtractView(&d_ptr);
  IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));


  for(int i=0; i<NumSweeps_; i++){
    // Calculate b-Ax
    if(!initial_guess_is_zero  || i > 0) {
      A_->Apply(Y,Temp);
      Temp.Update(1.0,*Xcopy,-1.0);
    }
    else
      Temp.Update(1.0,*Xcopy,0.0);

    // Note: The off-processor entries of T2 never get touched (they're always zero) and the other entries are updated
    // in this sweep before they are used, so we don't need to reset T2 to zero here.

    // Do backsolve & update
    // x = x  + W^{-1} (b - A x)
    for(int j=0; j<lclNumRows; j++){
      double diag=d_ptr[j];
      for (int m=0 ; m<NumVectors; m++) {
        double dtmp=0.0;
        // Note: Since the diagonal is in the matrix, we need to zero that entry of T2 here to make sure it doesn't contribute.
        t2_ptr[m][j]=0.0;
        for(int k=rowptr[j];k<rowptr[j+1];k++){
          dtmp+= values[k]*t2_ptr[m][colind[k]];
        }
        // Yes, we need to update both of these.
        t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;
        y_ptr[m][j] += omega*t2_ptr[m][j];
      }
    }
  }

  // Counter update
  NumApplyInverse_++;
  ApplyInverseTime_ += Time_.ElapsedTime();
  return 0;
}

void PeridigmNS::Block::createMapsFromGlobalMaps(Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarBondMap,
                                                 Teuchos::RCP<const Epetra_Vector>   globalBlockIds,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalNeighborhoodData,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalContactNeighborhoodData)
{
  double* globalBlockIdsPtr;
  globalBlockIds->ExtractView(&globalBlockIdsPtr);

  // Create a list of all the on-processor elements that are part of this block

  vector<int> IDs;
  IDs.reserve(globalOverlapScalarPointMap->NumMyElements()); // upper bound
  vector<int> bondIDs;
  bondIDs.reserve(globalOverlapScalarPointMap->NumMyElements());
  vector<int> bondElementSize;
  bondElementSize.reserve(globalOwnedScalarPointMap->NumMyElements());

  for(int iLID=0 ; iLID<globalOwnedScalarPointMap->NumMyElements() ; ++iLID){
    if(globalBlockIdsPtr[iLID] == blockID) {
      int globalID = globalOwnedScalarPointMap->GID(iLID);
      IDs.push_back(globalID);
    }
  }

  // Record the size of these elements in the bond map
  // Note that if an element has no bonds, it has no entry in the bondMap
  // So, the bond map and the scalar map can have a different number of entries (different local IDs)

  for(int iLID=0 ; iLID<globalOwnedScalarBondMap->NumMyElements() ; ++iLID){
    int globalID = globalOwnedScalarBondMap->GID(iLID);
    int localID = globalOwnedScalarPointMap->LID(globalID);
    if(globalBlockIdsPtr[localID] == blockID){
      bondIDs.push_back(globalID);
      bondElementSize.push_back(globalOwnedScalarBondMap->ElementSize(iLID));
    }
  }

  // Create the owned scalar point map, the owned vector point map, and the owned scalar bond map

  int numGlobalElements = -1;
  int numMyElements = IDs.size();
  int* myGlobalElements = 0;
  if(numMyElements > 0)
    myGlobalElements = &IDs.at(0);
  int elementSize = 1;
  int indexBase = 0;
  ownedScalarPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  elementSize = 3;
  ownedVectorPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  numMyElements = bondElementSize.size();
  myGlobalElements = 0;
  int* elementSizeList = 0;
  if(numMyElements > 0){
    myGlobalElements = &bondIDs.at(0);
    elementSizeList = &bondElementSize.at(0);
  }
  ownedScalarBondMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSizeList, indexBase, globalOwnedScalarPointMap->Comm()));

  // Create a list of nodes that need to be ghosted (both across material boundaries and across processor boundaries)
  set<int> ghosts;

  // Check the neighborhood list for things that need to be ghosted
  int* const globalNeighborhoodList = globalNeighborhoodData->NeighborhoodList();
  int globalNeighborhoodListIndex = 0;
  for(int iLID=0 ; iLID<globalNeighborhoodData->NumOwnedPoints() ; ++iLID){
    int numNeighbors = globalNeighborhoodList[globalNeighborhoodListIndex++];
    if(globalBlockIdsPtr[iLID] == blockID) {
      for(int i=0 ; i<numNeighbors ; ++i){
        int neighborGlobalID = globalOverlapScalarPointMap->GID( globalNeighborhoodList[globalNeighborhoodListIndex + i] );
        ghosts.insert(neighborGlobalID);
      }
    }
    globalNeighborhoodListIndex += numNeighbors;
  }

  // Check the contact neighborhood list for things that need to be ghosted
  if(!globalContactNeighborhoodData.is_null()){
    int* const globalContactNeighborhoodList = globalContactNeighborhoodData->NeighborhoodList();
    int globalContactNeighborhoodListIndex = 0;
    for(int iLID=0 ; iLID<globalContactNeighborhoodData->NumOwnedPoints() ; ++iLID){
      int numNeighbors = globalContactNeighborhoodList[globalContactNeighborhoodListIndex++];
      if(globalBlockIdsPtr[iLID] == blockID) {
        for(int i=0 ; i<numNeighbors ; ++i){
          int neighborGlobalID = globalOverlapScalarPointMap->GID( globalContactNeighborhoodList[globalContactNeighborhoodListIndex + i] );
          ghosts.insert(neighborGlobalID);
        }
      }
      globalContactNeighborhoodListIndex += numNeighbors;
    }
  }

//.........这里部分代码省略.........