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

 void Jacobian_2(Teuchos::RCP<Intrepid::FieldContainer<Real> > & jac,
                 const Teuchos::RCP<const Intrepid::FieldContainer<Real> > & u_coeff,
                 const Teuchos::RCP<const Intrepid::FieldContainer<Real> > & z_coeff = Teuchos::null,
                 const Teuchos::RCP<const std::vector<Real> > & z_param = Teuchos::null) {
   if ( z_coeff != Teuchos::null ) {
     // GET DIMENSIONS
     int c = u_coeff->dimension(0);
     int f = basisPtr_->getCardinality();
     // INITIALIZE JACOBIAN
     jac = Teuchos::rcp(new Intrepid::FieldContainer<Real>(c, f, f));
     // ADD CONTROL TERM
     Intrepid::FunctionSpaceTools::integrate<Real>(*jac, *(fe_vol_->N()), *(fe_vol_->NdetJ()), Intrepid::COMP_CPP, false);
     // APPLY DIRICHLET CONDITIONS
     int numSideSets = bdryCellLocIds_.size();
     if (numSideSets > 0) {
       for (int i = 0; i < numSideSets; ++i) {
         int numLocalSideIds = bdryCellLocIds_[i].size();
         for (int j = 0; j < numLocalSideIds; ++j) {
           int numCellsSide = bdryCellLocIds_[i][j].size();
           int numBdryDofs = fidx_[j].size();
           for (int k = 0; k < numCellsSide; ++k) {
             int cidx = bdryCellLocIds_[i][j][k];
             for (int l = 0; l < numBdryDofs; ++l) {
               //std::cout << "\n   j=" << j << "  l=" << l << "  " << fidx[j][l];
               for (int m = 0; m < f; ++m) {
                 (*jac)(cidx,fidx_[j][l],m) = static_cast<Real>(0);
               }
             }
           }
         }
       }
     }
   }
   else {
     throw Exception::Zero(">>> (PDE_Poisson::Jacobian_2): Jacobian is zero.");
   }
 }

/*----------------------------------------------------------------------*
 | project nodes master to slave along slave cont. normal field         |
 *----------------------------------------------------------------------*/
bool MOERTEL::Interface::ProjectNodes_MastertoSlave_NormalField()
{ 
  if (!IsComplete())
  {
	std::stringstream oss;
    oss << "***ERR*** MOERTEL::Interface::ProjectNodes_MastertoSlave_NormalField:\n"
         << "***ERR*** Complete() not called on interface " << Id() << "\n"
         << "***ERR*** file/line: " << __FILE__ << "/" << __LINE__ << "\n";
	throw ReportError(oss);
  }
  if (!lComm()) return true;
  
  int mside = MortarSide();
  int sside = OtherSide(mside);

  // iterate over all nodes of the master side and project those belonging to me
  std::map<int,Teuchos::RCP<MOERTEL::Node> >::iterator mcurr;
  for (mcurr=rnode_[mside].begin(); mcurr!=rnode_[mside].end(); ++mcurr)
  {
	Teuchos::RCP<MOERTEL::Node> mnode = mcurr->second;
    if (NodePID(mnode->Id()) != lComm()->MyPID())
      continue;
    
    const double* mx = mnode->X();
    double mindist = 1.0e+20;
	Teuchos::RCP<MOERTEL::Node> closenode = Teuchos::null;
    
    // find a node on the slave side that is closest to me
	std::map<int,Teuchos::RCP<MOERTEL::Node> >::iterator scurr;
    for (scurr=rnode_[sside].begin(); scurr!=rnode_[sside].end(); ++scurr)
    {
	  Teuchos::RCP<MOERTEL::Node> snode = scurr->second;
      const double* sx = snode->X();
      
      // build distance | snode->X() - mnode->X() |
      double dist = 0.0;
      for (int i=0; i<3; ++i) dist += (mx[i]-sx[i])*(mx[i]-sx[i]);
      dist = sqrt(dist);
      if (dist < mindist)
      {
        mindist = dist;
        closenode = snode;
      }
    } 
    if (closenode == Teuchos::null)
    {
	  std::stringstream oss;
      oss << "***ERR*** MOERTEL::Interface::ProjectNodes_MastertoSlave_NormalField:\n"
           << "***ERR*** Weired: for master node " << mnode->Id() << " no closest master node found\n"
           << "***ERR*** file/line: " << __FILE__ << "/" << __LINE__ << "\n";
		throw ReportError(oss);
    }

#if 0
    cout << "snode     " << *mnode;
    cout << "closenode " << *closenode;
#endif    
    
    // get segments attached to closest node closenode
    int  nseg = closenode->Nseg();
    MOERTEL::Segment** segs = closenode->Segments(); 
    
    // create a projection operator
    MOERTEL::Projector projector(IsOneDimensional(),OutLevel());
    
    // loop these segments and find best projection
    double bestdist[2];
    const double tol = 0.2;
	double bestgap = 0.0, gap;
    bestdist[0] = bestdist[1] = 1.0e+20;
    MOERTEL::Segment* bestseg = NULL;
    for (int i=0; i<nseg; ++i)
    {
      // project the master node on the slave segment along the segments interpolated normal field
      double xi[2]; xi[0] = xi[1] = 0.0;
      projector.ProjectNodetoSegment_SegmentNormal(*mnode,*(segs[i]),xi,gap);
      
      // check whether xi is better then previous projections
      if (IsOneDimensional())
      {
        if (abs(xi[0]) < abs(bestdist[0])) 
	{ 
	  bestdist[0] = xi[0];
	  bestdist[1] = xi[1];
	  bestseg = segs[i];
	  bestgap = gap;
	}
      }
      else
      {
        double third = 1./3.;
        // it's 'inside' with some tolerance
        if ( xi[0]<=1.+tol && xi[1]<=abs(1.-xi[0])+tol && xi[0]>=0.-tol && xi[1]>=0.-tol )
        {
          // it's better in both directions
          if ( sqrt((xi[0]-third)*(xi[0]-third)) < sqrt((bestdist[0]-third)*(bestdist[0]-third)) &&
               sqrt((xi[1]-third)*(xi[1]-third)) < sqrt((bestdist[1]-third)*(bestdist[1]-third)) )
//.........这里部分代码省略.........