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

typename boost::disable_if< boost::is_same<T, Albany::AbstractSTKFieldContainer::ScalarFieldType>, void >::type
Albany::GenericSTKFieldContainer<Interleaved>::fillVectorHelper(Epetra_Vector& soln,
    T* solution_field,
    const Teuchos::RCP<Epetra_Map>& node_map,
    const stk::mesh::Bucket& bucket, int offset) {

  // Fill the result vector
  // Create a multidimensional array view of the
  // solution field data for this bucket of nodes.
  // The array is two dimensional ( Cartesian X NumberNodes )
  // and indexed by ( 0..2 , 0..NumberNodes-1 )

  BucketArray<T> solution_array(*solution_field, bucket);

  const int num_vec_components = solution_array.dimension(0);
  const int num_nodes_in_bucket = solution_array.dimension(1);

  stk::mesh::BulkData& mesh = solution_field->get_mesh();

  for(std::size_t i = 0; i < num_nodes_in_bucket; i++)  {

    const GO node_gid = mesh.identifier(bucket[i]) - 1;
#ifdef ALBANY_64BIT_INT
    int node_lid = node_map->LID(static_cast<long long int>(node_gid));
#else
    //      const unsigned node_gid = bucket[i].identifier();
    int node_lid = node_map->LID(node_gid);
#endif

    for(std::size_t j = 0; j < num_vec_components; j++)

      soln[getDOF(node_lid, offset + j)] = solution_array(j, i);

  }
}

void
DirichletField<PHAL::AlbanyTraits::Residual, Traits>::
evaluateFields(typename Traits::EvalData dirichletWorkset) {

#if defined(ALBANY_EPETRA)
  const Albany::NodalDOFManager& fieldDofManager = dirichletWorkset.disc->getDOFManager(this->field_name);
  Teuchos::RCP<const Epetra_Map> fiedNodeMap = dirichletWorkset.disc->getNodeMap(this->field_name);
  const std::vector<GO>& nsNodesGIDs = dirichletWorkset.disc->getNodeSetGIDs().find(this->nodeSetID)->second;
#endif
  Teuchos::RCP<Tpetra_Vector> pvecT = dirichletWorkset.distParamLib->get(this->field_name)->vector();
  Teuchos::ArrayRCP<const ST> pT = pvecT->get1dView();
  Teuchos::RCP<Tpetra_Vector> fT = dirichletWorkset.fT;
  Teuchos::RCP<const Tpetra_Vector> xT = dirichletWorkset.xT;
  Teuchos::ArrayRCP<const ST> xT_constView = xT->get1dView();
  Teuchos::ArrayRCP<ST> fT_nonconstView = fT->get1dViewNonConst();

  const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
  for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
      int lunk = nsNodes[inode][this->offset];
#if defined(ALBANY_EPETRA)
      GO node_gid = nsNodesGIDs[inode];
      int lfield = fieldDofManager.getLocalDOF(fiedNodeMap->LID(node_gid),this->offset);
      fT_nonconstView[lunk] = xT_constView[lunk] - pT[lfield];
#else
      fT_nonconstView[lunk] = xT_constView[lunk] - pT[lunk];
#endif
  }
}

Teuchos::RCP<PeridigmNS::NeighborhoodData> PeridigmNS::Block::createNeighborhoodDataFromGlobalNeighborhoodData(Teuchos::RCP<const Epetra_BlockMap> globalOverlapScalarPointMap,
                                                                                                               Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalNeighborhoodData)
{
  int numOwnedPoints = ownedScalarPointMap->NumMyElements();
  int* ownedPointGlobalIDs = ownedScalarPointMap->MyGlobalElements();

  vector<int> ownedIDs(numOwnedPoints);
  vector<int> neighborhoodList;
  vector<int> neighborhoodPtr(numOwnedPoints);

  int* const globalNeighborhoodList = globalNeighborhoodData->NeighborhoodList();
  int* const globalNeighborhoodPtr = globalNeighborhoodData->NeighborhoodPtr();

  // Create the neighborhoodList and neighborhoodPtr for this block.
  // All the IDs in the neighborhoodList and neighborhoodPtr are local IDs into 
  // the block-specific overlap map.

  for(int i=0 ; i<numOwnedPoints ; ++i){
    neighborhoodPtr[i] = (int)(neighborhoodList.size());
    int globalID = ownedPointGlobalIDs[i];
    ownedIDs[i] = overlapScalarPointMap->LID(globalID);
    int globalNeighborhoodListIndex = globalNeighborhoodPtr[globalOverlapScalarPointMap->LID(globalID)];
    int numNeighbors = globalNeighborhoodList[globalNeighborhoodListIndex++];
    neighborhoodList.push_back(numNeighbors);
    for(int j=0 ; j<numNeighbors ; ++j){
      int globalNeighborID = globalOverlapScalarPointMap->GID(globalNeighborhoodList[globalNeighborhoodListIndex++]);
      neighborhoodList.push_back( overlapScalarPointMap->LID(globalNeighborID) );
    }
  }

  // create the NeighborhoodData for this block

  Teuchos::RCP<PeridigmNS::NeighborhoodData> blockNeighborhoodData = Teuchos::rcp(new PeridigmNS::NeighborhoodData);
  blockNeighborhoodData->SetNumOwned(ownedIDs.size());
  if(ownedIDs.size() > 0){
    memcpy(blockNeighborhoodData->OwnedIDs(), 
           &ownedIDs.at(0),
           ownedIDs.size()*sizeof(int));
  }
  if(neighborhoodPtr.size() > 0){
    memcpy(blockNeighborhoodData->NeighborhoodPtr(), 
           &neighborhoodPtr.at(0),
           neighborhoodPtr.size()*sizeof(int));
  }
  blockNeighborhoodData->SetNeighborhoodListSize(neighborhoodList.size());
  if(neighborhoodList.size() > 0){
    memcpy(blockNeighborhoodData->NeighborhoodList(),
           &neighborhoodList.at(0),
           neighborhoodList.size()*sizeof(int));
  }

  return blockNeighborhoodData;
}

// The solve is done in the felix_driver_run function, and the solution is passed back to Glimmer-CISM 
// IK, 12/3/13: time_inc_yr and cur_time_yr are not used here... 
void felix_driver_run(FelixToGlimmer * ftg_ptr, double& cur_time_yr, double time_inc_yr)
{
    //IK, 12/9/13: how come FancyOStream prints an all processors??    
    Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());

    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) {
      std::cout << "In felix_driver_run, cur_time, time_inc = " << cur_time_yr 
                << "   " << time_inc_yr << std::endl;
    }
    
    // ---------------------------------------------
    // get u and v velocity solution from Glimmer-CISM 
    // IK, 11/26/13: need to concatenate these into a single solve for initial condition for Albany/FELIX solve 
    // IK, 3/14/14: moved this step to felix_driver_run from felix_driver init, since we still want to grab and u and v velocities for CISM if the mesh hasn't changed, 
    // in which case only felix_driver_run will be called, not felix_driver_init.   
    // ---------------------------------------------
    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) 
      std::cout << "In felix_driver_run: grabbing pointers to u and v velocities in CISM..." << std::endl; 
    uVel_ptr = ftg_ptr ->getDoubleVar("uvel", "velocity"); 
    vVel_ptr = ftg_ptr ->getDoubleVar("vvel", "velocity"); 

    // ---------------------------------------------
    // Set restart solution to the one passed from CISM
    // IK, 3/14/14: moved this from felix_driver_init to felix_driver_run.  
    // ---------------------------------------------
    
    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) 
      std::cout << "In felix_driver_run: setting initial condition from CISM..." << std::endl;
    //Check what kind of ordering you have in the solution & create solutionField object.
    interleavedOrdering = meshStruct->getInterleavedOrdering();
    Albany::AbstractSTKFieldContainer::VectorFieldType* solutionField;
    if(interleavedOrdering)
      solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<true> >(meshStruct->getFieldContainer())->getSolutionField();
    else
      solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<false> >(meshStruct->getFieldContainer())->getSolutionField();

     //Create vector used to renumber nodes on each processor from the Albany convention (horizontal levels first) to the CISM convention (vertical layers first)
     nNodes2D = (global_ewn + 1)*(global_nsn+1); //number global nodes in the domain in 2D 
     nNodesProc2D = (nsn-2*nhalo+1)*(ewn-2*nhalo+1); //number of nodes on each processor in 2D  
     cismToAlbanyNodeNumberMap.resize(upn*nNodesProc2D);
     for (int j=0; j<nsn-2*nhalo+1;j++) { 
       for (int i=0; i<ewn-2*nhalo+1; i++) {
         for (int k=0; k<upn; k++) { 
           int index = k+upn*i + j*(ewn-2*nhalo+1)*upn; 
           cismToAlbanyNodeNumberMap[index] = k*nNodes2D + global_node_id_owned_map_Ptr[i+j*(ewn-2*nhalo+1)]; 
           //if (mpiComm->MyPID() == 0) 
           //  std::cout << "index: " << index << ", cismToAlbanyNodeNumberMap: " << cismToAlbanyNodeNumberMap[index] << std::endl; 
          }
        }
      }

     //The way it worked out, uVel_ptr and vVel_ptr have more nodes than the nodes in the mesh passed to Albany/CISM for the solve.  In particular, 
     //there is 1 row of halo elements in uVel_ptr and vVel_ptr.  To account for this, we copy uVel_ptr and vVel_ptr into std::vectors, which do not have the halo elements. 
     std::vector<double> uvel_vec(upn*nNodesProc2D); 
     std::vector<double> vvel_vec(upn*nNodesProc2D); 
     int counter1 = 0; 
     int counter2 = 0; 
     int local_nodeID; 
     for (int j=0; j<nsn-1; j++) {
       for (int i=0; i<ewn-1; i++) { 
         for (int k=0; k<upn; k++) {
           if (j >= nhalo-1 & j < nsn-nhalo) {
             if (i >= nhalo-1 & i < ewn-nhalo) {
#ifdef CISM_USE_EPETRA 
               local_nodeID = node_map->LID(cismToAlbanyNodeNumberMap[counter1]); 
#else
               local_nodeID = node_map->getLocalElement(cismToAlbanyNodeNumberMap[counter1]);
#endif
               uvel_vec[counter1] = uVel_ptr[counter2]; 
               vvel_vec[counter1] = vVel_ptr[counter2]; 
               counter1++;
            }
            }
            counter2++; 
         }
        }
     }
     //Loop over all the elements to find which nodes are active.  For the active nodes, copy uvel and vvel from CISM into Albany solution array to 
     //use as initial condition.
     //NOTE: there is some inefficiency here by looping over all the elements.  TO DO? pass only active nodes from Albany-CISM to improve this? 
     double velScale = seconds_per_year*vel_scaling_param;  
     for (int i=0; i<nElementsActive; i++) {
       for (int j=0; j<8; j++) {
        int node_GID =  global_element_conn_active_Ptr[i + nElementsActive*j]; //node_GID is 1-based
#ifdef CISM_USE_EPETRA      
        int node_LID =  node_map->LID(node_GID); //node_LID is 0-based
#else
        int node_LID =  node_map->getLocalElement(node_GID); //node_LID is 0-based
#endif
        stk::mesh::Entity node = meshStruct->bulkData->get_entity(stk::topology::NODE_RANK, node_GID);
        double* sol = stk::mesh::field_data(*solutionField, node);
        //IK, 3/18/14: added division by velScale to convert uvel and vvel from dimensionless to having units of m/year (the Albany units)  
        sol[0] = uvel_vec[node_LID]/velScale;
        sol[1] = vvel_vec[node_LID]/velScale;
      }
    }
    // ---------------------------------------------------------------------------------------------------
    // Solve 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    int il = (ordering == 0) * (j / (lElemColumnShift / 3))
        + (ordering == 1) * (j % (elemLayerShift / 3));
    int gId = il * elemColumnShift + elemLayerShift * indexToTriangleID[ib];
    int lId = il * lElemColumnShift + elemLayerShift * ib;
    for (int iTetra = 0; iTetra < 3; iTetra++) {
      stk::mesh::Entity elem = meshStruct->bulkData->get_entity(stk::topology::ELEMENT_RANK, ++gId);
      double* temperature = stk::mesh::field_data(*temperature_field, elem);
      temperature[0] = temperatureOnTetra[lId++];
    }
  }

  meshStruct->setHasRestartSolution(!first_time_step);

  if (!first_time_step) {
    meshStruct->setRestartDataTime(
        paramList->sublist("Problem").get("Homotopy Restart Step", 1.));
    double homotopy =
        paramList->sublist("Problem").sublist("FELIX Viscosity").get(
            "Glen's Law Homotopy Parameter", 1.0);
    if (meshStruct->restartDataTime() == homotopy) {
      paramList->sublist("Problem").set("Solution Method", "Steady");
      paramList->sublist("Piro").set("Solver Type", "NOX");
    }
  }



  if(!keptMesh) {
    albanyApp->createDiscretization();
    albanyApp->finalSetUp(paramList); //, albanyApp->getDiscretization()->getSolutionFieldT());
  }
  else
    albanyApp->getDiscretization()->updateMesh();

#ifdef MPAS_USE_EPETRA
  solver = slvrfctry->createThyraSolverAndGetAlbanyApp(albanyApp, mpiComm, mpiComm, Teuchos::null, false);
#else
   solver = slvrfctry->createAndGetAlbanyAppT(albanyApp, mpiCommT, mpiCommT, Teuchos::null, false);
#endif

  Teuchos::ParameterList solveParams;
  solveParams.set("Compute Sensitivities", false);

  Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > thyraResponses;
  Teuchos::Array<
      Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > thyraSensitivities;
  Piro::PerformSolveBase(*solver, solveParams, thyraResponses,
      thyraSensitivities);
#ifdef MPAS_USE_EPETRA
  const Epetra_Map& overlapMap(
      *albanyApp->getDiscretization()->getOverlapMap());
  Epetra_Import import(overlapMap, *albanyApp->getDiscretization()->getMap());
  Epetra_Vector solution(overlapMap);
  solution.Import(*albanyApp->getDiscretization()->getSolutionField(), import,
      Insert);
#else
  Teuchos::RCP<const Tpetra_Map> overlapMap = albanyApp->getDiscretization()->getOverlapMapT();
  Teuchos::RCP<Tpetra_Import> import = Teuchos::rcp(new Tpetra_Import(albanyApp->getDiscretization()->getMapT(), overlapMap));
  Teuchos::RCP<Tpetra_Vector> solution = Teuchos::rcp(new Tpetra_Vector(overlapMap));
  solution->doImport(*albanyApp->getDiscretization()->getSolutionFieldT(), *import, Tpetra::INSERT);
  Teuchos::ArrayRCP<const ST> solution_constView = solution->get1dView();
#endif


  for (UInt j = 0; j < numVertices3D; ++j) {
    int ib = (ordering == 0) * (j % lVertexColumnShift)
        + (ordering == 1) * (j / vertexLayerShift);
    int il = (ordering == 0) * (j / lVertexColumnShift)
        + (ordering == 1) * (j % vertexLayerShift);
    int gId = il * vertexColumnShift + vertexLayerShift * indexToVertexID[ib];

    int lId0, lId1;

#ifdef MPAS_USE_EPETRA
    if (interleavedOrdering) {
      lId0 = overlapMap.LID(2 * gId);
      lId1 = lId0 + 1;
    } else {
      lId0 = overlapMap.LID(gId);
      lId1 = lId0 + numVertices3D;
    }
    velocityOnVertices[j] = solution[lId0];
    velocityOnVertices[j + numVertices3D] = solution[lId1];
#else
    if (interleavedOrdering) {
      lId0 = overlapMap->getLocalElement(2 * gId);
      lId1 = lId0 + 1;
    } else {
      lId0 = overlapMap->getLocalElement(gId);
      lId1 = lId0 + numVertices3D;
    }
    velocityOnVertices[j] = solution_constView[lId0];
    velocityOnVertices[j + numVertices3D] = solution_constView[lId1];
#endif
  }

  keptMesh = true;

  //UInt componentGlobalLength = (nLayers+1)*nGlobalVertices; //mesh3DPtr->numGlobalVertices();
}