C++ Basis::getBasis方法代码示例

// Matrix and Residual Fills
bool
Pitchfork_FiniteElementProblem::evaluate(FillType f,
                     const Epetra_Vector* soln,
                     Epetra_Vector* tmp_rhs,
                     Epetra_RowMatrix* tmp_matrix,
                     double jac_coeff,
                     double mass_coeff)
{
  flag = f;

  // Set the incoming linear objects
  if (flag == F_ONLY) {
    rhs = tmp_rhs;
  } else if (flag == MATRIX_ONLY) {
    A = dynamic_cast<Epetra_CrsMatrix*> (tmp_matrix);
    assert(A != NULL);
  } else if (flag == ALL) {
    rhs = tmp_rhs;
    A = dynamic_cast<Epetra_CrsMatrix*> (tmp_matrix);
    assert(A != NULL);
  } else {
    std::cout << "ERROR: Pitchfork_FiniteElementProblem::fillMatrix() - FillType flag is broken" << std::endl;
    throw;
  }

  // Create the overlapped solution and position vectors
  Epetra_Vector u(*OverlapMap);
  Epetra_Vector x(*OverlapMap);

  // Export Solution to Overlap vector
  u.Import(*soln, *Importer, Insert);

  // Declare required variables
  int i,j,ierr;
  int OverlapNumMyElements = OverlapMap->NumMyElements();

  int OverlapMinMyGID;
  if (MyPID==0) OverlapMinMyGID = StandardMap->MinMyGID();
  else OverlapMinMyGID = StandardMap->MinMyGID()-1;

  int row, column;
  double jac;
  double xx[2];
  double uu[2];
  Basis basis;

  // Create the nodal coordinates
  double Length=2.0;
  double dx=Length/((double) NumGlobalElements-1);
  for (i=0; i < OverlapNumMyElements; i++) {
    x[i]=-1.0 + dx*((double) OverlapMinMyGID+i);
  }

  // Zero out the objects that will be filled
  if ((flag == MATRIX_ONLY) || (flag == ALL)) {
    i = A->PutScalar(0.0);
    assert(i == 0);
  }
  if ((flag == F_ONLY)    || (flag == ALL)) {
    i = rhs->PutScalar(0.0);
    assert(i == 0);
  }

  // Loop Over # of Finite Elements on Processor
  for (int ne=0; ne < OverlapNumMyElements-1; ne++) {

    // Loop Over Gauss Points
    for(int gp=0; gp < 2; gp++) {
      // Get the solution and coordinates at the nodes
      xx[0]=x[ne];
      xx[1]=x[ne+1];
      uu[0]=u[ne];
      uu[1]=u[ne+1];
      // Calculate the basis function at the gauss point
      basis.getBasis(gp, xx, uu);

      // Loop over Nodes in Element
      for (i=0; i< 2; i++) {
    row=OverlapMap->GID(ne+i);
    //printf("Proc=%d GlobalRow=%d LocalRow=%d Owned=%d\n",
    //     MyPID, row, ne+i,StandardMap.MyGID(row));
    if (StandardMap->MyGID(row)) {
      if ((flag == F_ONLY)    || (flag == ALL)) {
        (*rhs)[StandardMap->LID(OverlapMap->GID(ne+i))]+=
          +basis.wt*basis.dx
          *((-1.0/(basis.dx*basis.dx))*basis.duu*
        basis.dphide[i]-source_term(basis.uu)*basis.phi[i]);
      }
    }
    // Loop over Trial Functions
    if ((flag == MATRIX_ONLY) || (flag == ALL)) {
      for(j=0;j < 2; j++) {
        if (StandardMap->MyGID(row)) {
          column=OverlapMap->GID(ne+j);
          jac=jac_coeff*basis.wt*basis.dx*
        ((-1.0/(basis.dx*basis.dx))*basis.dphide[j]*basis.dphide[i]
         -source_deriv(basis.uu)*basis.phi[j]*basis.phi[i]) +
        mass_coeff*basis.wt*basis.dx*basis.phi[j]*basis.phi[i];
          ierr=A->SumIntoGlobalValues(row, 1, &jac, &column);
//.........这里部分代码省略.........

// A fill specialized to the single node at the coupling interface
void 
ConvDiff_PDE::computeHeatFlux( const Epetra_Vector * soln )
{
  
  int numDep = depProblems.size();

  // Create the overlapped solution and position vectors
  Epetra_Vector u(*OverlapMap);
  Epetra_Vector uold(*OverlapMap);
  std::vector<Epetra_Vector*> dep(numDep);
  for( int i = 0; i < numDep; ++i)
    dep[i] = new Epetra_Vector(*OverlapMap);

  Epetra_Vector xvec(*OverlapMap);

  // Export Solution to Overlap vector
  // If the vector to be used in the fill is already in the Overlap form,
  // we simply need to map on-processor from column-space indices to
  // OverlapMap indices. Note that the old solution is simply fixed data that
  // needs to be sent to an OverlapMap (ghosted) vector.  The conditional
  // treatment for the current soution vector arises from use of
  // FD coloring in parallel.

  uold.Import(*oldSolution, *Importer, Insert);

  for( int i = 0; i < numDep; ++i )
    (*dep[i]).Import(*( (*(depSolutions.find(depProblems[i]))).second ), *Importer, Insert);

  xvec.Import(*xptr, *Importer, Insert);

  if( NULL == soln )
    u.Import(*initialSolution, *Importer, Insert);
  else
    u.Import(*soln, *Importer, Insert);

  // Declare required variables
  int row;
  double * xx = new double[2];
  double * uu = new double[2]; 
  double * uuold = new double[2];
  std::vector<double*> ddep(numDep);
  for( int i = 0; i < numDep; ++i)
    ddep[i] = new double[2];

  Basis basis;

  // Bundle up the dependent variables in the way needed for computing
  // the source terms of each reaction
  map<string, double*> depVars;
  depVars.insert( pair< std::string, double*>(getName(), uu) );
  for( int i = 0; i < numDep; ++i )
    depVars.insert( pair<string, double*>(myManager->getProblemName(depProblems[i]), ddep[i]) );

  myFlux = 0.0;

  // Loop Over Gauss Points
  for( int gp = 0; gp < 2; ++gp ) 
  {
    // Get the solution and coordinates at the nodes 
    xx[0]=xvec[interface_elem];
    xx[1]=xvec[interface_elem+1];
    uu[0] = u[interface_elem];
    uu[1] = u[interface_elem+1];
    uuold[0] = uold[interface_elem];
    uuold[1] = uold[interface_elem+1];
    for( int i = 0; i < numDep; ++i ) 
    {
      ddep[i][0] = (*dep[i])[interface_elem];
      ddep[i][1] = (*dep[i])[interface_elem+1];
    }

    // Calculate the basis function and variables at the gauss points
    basis.getBasis(gp, xx, uu, uuold, ddep);

    row = OverlapMap->GID( interface_elem + local_node );

    if( StandardMap->MyGID(row) ) 
    {
      myFlux += 
        + basis.wt * basis.dx
        * ( peclet * (basis.duu / basis.dx) * basis.phi[local_node] 
        +   kappa * (1.0/(basis.dx*basis.dx)) * basis.duu * basis.dphide[local_node] );
    }
  }

  // Sync up processors to be safe
  Comm->Barrier();
 
  // Cleanup
  for( int i = 0; i < numDep; ++i)
  {
    delete [] ddep[i];
    delete     dep[i];
  }

  delete [] xx    ;
  delete [] uu    ;
  delete [] uuold ;

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

//.........这里部分代码省略.........
                    (myManager->getProblemName(depProblems[i]), ddep[i]) );
  // Do a check on this fill
//  map<string, double*>::iterator iter;
//  for( iter = depVars.begin(); iter != depVars.end(); iter++)
//    std::cout << "Inserted ... " << iter->first << "\t" << iter->second << std::endl;
//  std::cout << "--------------------------------------------------" << std::endl;
//  for( iter = depVars.begin(); iter != depVars.end(); iter++)
//	  std::cout << iter->first << "\t" << (iter->second)[0] << ", " 
//               << (iter->second)[1] << std::endl;
//  std::cout << "--------------------------------------------------" << std::endl;

  // Zero out the objects that will be filled
  if ( fillMatrix ) A->PutScalar(0.0);
  if ( fillF ) rhs->PutScalar(0.0);

  // Loop Over # of Finite Elements on Processor
  for (int ne=0; ne < OverlapNumMyNodes-1; ne++) 
  {
    // Loop Over Gauss Points
    for(int gp=0; gp < 2; gp++) 
    {
      // Get the solution and coordinates at the nodes 
      xx[0]=xvec[ne];
      xx[1]=xvec[ne+1];
      uu[0] = u[ne];
      uu[1] = u[ne+1];
      uuold[0] = uold[ne];
      uuold[1] = uold[ne+1];
      for( int i = 0; i<numDep; i++ ) {
        ddep[i][0] = (*dep[i])[ne];
        ddep[i][1] = (*dep[i])[ne+1];
      }
      // Calculate the basis function and variables at the gauss points
      basis.getBasis(gp, xx, uu, uuold, ddep);

      // Loop over Nodes in Element
      for (int i=0; i< 2; i++) {
	row=OverlapMap->GID(ne+i);
	if (StandardMap->MyGID(row)) {
	  if ( fillF ) {

            // First do time derivative and diffusion operator
	    (*rhs)[StandardMap->LID(OverlapMap->GID(ne+i))]+=
	      +basis.wt*basis.dx
	      *((basis.uu - basis.uuold)/dt * basis.phi[i] 
              +(1.0/(basis.dx*basis.dx))*diffCoef*basis.duu*basis.dphide[i]);

            // Then do source term contributions
	    //
            for( srcTermIter = SrcTermWeight.begin(); 
                          srcTermIter != srcTermEnd; srcTermIter++) {
              HMX_PDE &srcTermProb = 
                       dynamic_cast<HMX_PDE&>(
			       myManager->getProblem((*srcTermIter).first) );
              srcTermProb.computeSourceTerm(2, depVars, srcTerm);
              (*rhs)[StandardMap->LID(OverlapMap->GID(ne+i))]+=
                +basis.wt*basis.dx
                *( basis.phi[i] * ( - (*srcTermIter).second * srcTerm[i] ));
            }
	    //
	  }
	}
	// Loop over Trial Functions
	if ( fillMatrix ) {
		/*
	  for(j=0;j < 2; j++) {

// Matrix and Residual Fills
bool FiniteElementProblem::
evaluate(FillType f, 
	 const Epetra_Vector* soln, 
	 Epetra_Vector* tmp_rhs, 
	 Epetra_RowMatrix* tmp_matrix,
	 NOX::Epetra::Interface::Required::FillType fillFlag)
{
  flag = f;

  // Set the incoming linear objects
  if (flag == F_ONLY) {
    rhs = tmp_rhs;
  } else if (flag == MATRIX_ONLY) {
    //A = dynamic_cast<Epetra_CrsMatrix*> (tmp_matrix);
  } else if (flag == ALL) { 
    rhs = tmp_rhs;
    //A = dynamic_cast<Epetra_CrsMatrix*> (tmp_matrix);
  } else {
    cout << "ERROR: FiniteElementProblem::fillMatrix() - FillType flag is broken" << endl;
    throw;
  }

  // Create the overlapped solution and position vectors
  Epetra_Vector u(*OverlapMap);
  Epetra_Vector x(*OverlapMap);

  // Export Solution to Overlap vector
  // If the vector to be used in the fill is already in the Overlap form,
  // we simply need to map on-processor from column-space indices to
  // OverlapMap indices.
  if( fillFlag == NOX::Epetra::Interface::Required::FD_Res)
  {
    u.Export(*soln, *ColumnToOverlapImporter, Insert);
  }
  // If the incoming vector used in the fill has not been scattered to 
  // OverlapMap (ghosted) entities, do it now.  This potentially involves
  // interprocessor communication unlike the previous condition.
  else
    u.Import(*soln, *Importer, Insert);

  // Declare required variables
  int i,j,ierr;
  int OverlapNumMyElements = OverlapMap->NumMyElements();

  int OverlapMinMyGID;
  if (MyPID==0) OverlapMinMyGID = StandardMap->MinMyGID();
  else OverlapMinMyGID = StandardMap->MinMyGID()-1;

  int row, column;
  double factor=1000.0;
  double jac;
  double xx[2];
  double uu[2];
  Basis basis;

  // Create the nodal coordinates
  double Length=1.0;
  double dx=Length/((double) NumGlobalElements-1);
  for (i=0; i < OverlapNumMyElements; i++) {
    x[i]=dx*((double) OverlapMinMyGID+i);
  }
  
  // Zero out the objects that will be filled
  if ((flag == MATRIX_ONLY) || (flag == ALL)) i=A->PutScalar(0.0);
  if ((flag == F_ONLY)    || (flag == ALL)) i=rhs->PutScalar(0.0);

  // Loop Over # of Finite Elements on Processor
  for (int ne=0; ne < OverlapNumMyElements-1; ne++) {
    
    // Loop Over Gauss Points
    for(int gp=0; gp < 2; gp++) {
      // Get the solution and coordinates at the nodes 
      xx[0]=x[ne];
      xx[1]=x[ne+1];
      uu[0]=u[ne];
      uu[1]=u[ne+1];
      // Calculate the basis function at the gauss point
      basis.getBasis(gp, xx, uu);
	            
      // Loop over Nodes in Element
      for (i=0; i< 2; i++) {
	row=OverlapMap->GID(ne+i);
	if (StandardMap->MyGID(row)) {
	  if ((flag == F_ONLY)    || (flag == ALL)) {
	    (*rhs)[StandardMap->LID(OverlapMap->GID(ne+i))]+=
	      +basis.wt*basis.dx
	      *((1.0/(basis.dx*basis.dx))*basis.duu*
		basis.dphide[i]+factor*basis.uu*basis.uu*basis.phi[i]);
	  }
	}
	// Loop over Trial Functions
	if ((flag == MATRIX_ONLY) || (flag == ALL)) {
	  for(j=0;j < 2; j++) {
	    if (StandardMap->MyGID(row)) {
	      column=OverlapMap->GID(ne+j);
	      jac=basis.wt*basis.dx*((1.0/(basis.dx*basis.dx))*
				     basis.dphide[j]*basis.dphide[i]
				     +2.0*factor*basis.uu*basis.phi[j]*
				     basis.phi[i]);  
//.........这里部分代码省略.........