C++ UserVector::Update方法代码示例

void AdaptiveSparseGridInterface<Scalar,UserVector>::eval_cubature(
	       UserVector & output, 
	       CubatureTensorSorted<Scalar> & cubRule) {

  //int dimf      = 0;                      // Dimension of the integrand
  Scalar weight = 0.0;
  std::vector<Scalar> point(dimension_,(Scalar)0.0);
  //std::vector<Scalar> f(1,0.0);
  Teuchos::RCP<UserVector> f = output.Create(); output.Update(-1.0,output);

  typename std::map<std::vector<Scalar>,int>::iterator it;
  for (it=cubRule.begin(); it!=cubRule.end(); it++) {
    // Evaluate Function
    point.assign((it->first).begin(),(it->first).end()); // Extract point
    f->Update(-1.0,*f);
    eval_integrand(*f,point);     // Evaluate Integrand at point
 
    // Update integral
    weight = cubRule.getWeight(it->second);
    output.Update(weight,*f);
  }
}

void AdaptiveSparseGridInterface<Scalar,UserVector>::init(UserVector & output) {
  std::vector<int> index(dimension_,1);
  CubatureTensorSorted<Scalar> cubRule(
            dimension_,index,rule1D_,growth1D_,isNormalized_);
  
  // Evaluate the initial contribution to the integral
  initialDiff_ = 1.0;
  output.Update(-1.0,output);
  eval_cubature(output,cubRule);

  // Compute the initial error indicator
  initialDiff_ = error_indicator(output);
  if (fabs(initialDiff_)<INTREPID_TOL) 
    initialDiff_ = 1.0;
}

Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
	typename std::multimap<Scalar,std::vector<int> > & activeIndex, 
	std::set<std::vector<int> > & oldIndex, 
	UserVector & integralValue,
	CubatureTensorSorted<Scalar> & cubRule,
	Scalar globalErrorIndicator,
	AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {

  TEUCHOS_TEST_FOR_EXCEPTION((activeIndex.empty()),std::out_of_range,
              ">>> ERROR (AdaptiveSparseGrid): Active Index set is empty.");  

  int dimension = problem_data.getDimension();
  std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
  std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);

  // Initialize Flags
  bool maxLevelFlag     = true;
  bool isAdmissibleFlag = true;

  // Initialize Cubature Rule
  Teuchos::RCP<UserVector> s = integralValue.Create();

  // Initialize iterator at end of inOldIndex
  std::set<std::vector<int> >::iterator it1(oldIndex.end());  

  // Initialize iterator at end of inActiveIndex
  typename std::multimap<Scalar,std::vector<int> >::iterator it;

  // Obtain Global Error Indicator as sum of key values of inActiveIndex
  Scalar eta = globalErrorIndicator;

  // Select Index to refine
  it = activeIndex.end(); it--;        // Decrement to position of final value 
  Scalar G               = it->first;  // Largest Error Indicator is at End
  eta                   -= G;          // Update global error indicator
  std::vector<int> index = it->second; // Get Corresponding index
  activeIndex.erase(it);               // Erase Index from active index set
  // Insert Index into old index set
  oldIndex.insert(it1,index); it1 = oldIndex.end(); 
  
  // Refinement process
  for (int k=0; k<dimension; k++) {
    index[k]++; // index + ek
    // Check Max Level
    maxLevelFlag = problem_data.max_level(index);
    if (maxLevelFlag) {
      // Check Admissibility
      isAdmissibleFlag = isAdmissible(index,k,oldIndex,problem_data);
      if (isAdmissibleFlag) { // If admissible
	// Build Differential Quarature Rule
	CubatureTensorSorted<Scalar> diffRule(0,dimension);
	build_diffRule(diffRule,index,problem_data);
	
	// Apply Rule to function
	problem_data.eval_cubature(*s,diffRule);
	
	// Update integral value
	integralValue.Update(*s);
	
	// Update local error indicator and index set
	G  = problem_data.error_indicator(*s); 	
	if (activeIndex.end()!=activeIndex.begin()) 
	  activeIndex.insert(activeIndex.end()--,
			   std::pair<Scalar,std::vector<int> >(G,index));
	else
	  activeIndex.insert(std::pair<Scalar,std::vector<int> >(G,index));
		
	// Update global error indicators
	eta += G;

	// Update adapted quadrature rule nodes and weights
	cubRule.update(1.0,diffRule,1.0);
      }
    }
    else { // Max Level Exceeded 
      //std::cout << "Maximum Level Exceeded" << std::endl;
    }
    index[k]--;
  }
  return eta;
}