libtrilinos-intrepid-dev 12.12.1-5 / usr / include / trilinos / Intrepid_AdaptiveSparseGridDef.hpp

// @HEADER
// ************************************************************************
//
//                           Intrepid Package
//                 Copyright (2007) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
//                    Denis Ridzal  (dridzal@sandia.gov), or
//                    Kara Peterson (kjpeter@sandia.gov)
//
// ************************************************************************
// @HEADER

/** \file   Intrepid_AdaptiveSparseGridDef.hpp
    \brief  Definition file for the Intrepid::AdaptiveSparseGrid class.
    \author Created by D. Kouri and D. Ridzal.
*/

namespace Intrepid {

template<class Scalar, class UserVector>
bool AdaptiveSparseGrid<Scalar,UserVector>::isAdmissible(
               std::vector<int> index, 
	       int direction, 
	       std::set<std::vector<int> > inOldIndex,
	       AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
  /*
    Check if inOldIndex remains admissible if index is added, i.e.
       index-ek in inOldIndex for all k=1,...,dim.
  */
  int dimension = problem_data.getDimension();
  for (int i=0; i<dimension; i++) {
    if (index[i]>1 && i!=direction) {
      index[i]--;
      if (!inOldIndex.count(index)) {
	return false;
      }
      index[i]++;
    }
  }
  return true;
}

template<class Scalar, class UserVector>
void AdaptiveSparseGrid<Scalar,UserVector>::build_diffRule(
	       CubatureTensorSorted<Scalar> & outRule, 
	       std::vector<int> index,
	       AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
	       
  int numPoints     = 0;
  int dimension     = problem_data.getDimension();
  bool isNormalized = problem_data.isNormalized();
  std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
  std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);

  for (int i=0; i<dimension; i++) {
    // Compute 1D rules
    numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
    CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
 
    if (numPoints!=1) { // Compute differential rule
      numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
      CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
      diffRule1.update(-1.0,rule1,1.0);
    }
    // Build Tensor Rule
    outRule = kron_prod<Scalar>(outRule,diffRule1); 
  }
}

template<class Scalar, class UserVector>
void AdaptiveSparseGrid<Scalar,UserVector>::build_diffRule(
	       CubatureTensorSorted<Scalar> & outRule, 
	       std::vector<int> index, int dimension,
	       std::vector<EIntrepidBurkardt> rule1D, 
	       std::vector<EIntrepidGrowth> growth1D,
	       bool isNormalized) {
	       
  int numPoints = 0;
  for (int i=0; i<dimension; i++) {
    // Compute 1D rules
    numPoints = growthRule1D(index[i],growth1D[i],rule1D[i]);
    CubatureLineSorted<Scalar> diffRule1(rule1D[i],numPoints,isNormalized);
 
    if (numPoints!=1) { // Differential Rule
      numPoints = growthRule1D(index[i]-1,growth1D[i],rule1D[i]);
      CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
      diffRule1.update(-1.0,rule1,1.0);
    }
    // Build Tensor Rule
    outRule = kron_prod<Scalar>(outRule,diffRule1); 
  }
}
 
// Update Index Set - no knowledge of active or old indices
template<class Scalar, class UserVector>
Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
   	         typename std::multimap<Scalar,std::vector<int> >  & indexSet, 
		 UserVector & integralValue,
		 AdaptiveSparseGridInterface<Scalar,UserVector> & problem_data) {
  
  int dimension = problem_data.getDimension();
  std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
  std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
  
  // Copy Multimap into a Set for ease of use
  typename std::multimap<Scalar,std::vector<int> >::iterator it;
  std::set<std::vector<int> > oldSet;
  std::set<std::vector<int> >::iterator it1(oldSet.begin());
  for (it=indexSet.begin(); it!=indexSet.end(); it++) {
    oldSet.insert(it1,it->second);
    it1++;
  }
  indexSet.clear();
  
  // Find Possible Active Points
  int flag = 1;
  std::vector<int> index(dimension,0);
  typename std::multimap<Scalar,std::vector<int> > activeSet;
  for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
    index = *it1;
    for (int i=0; i<dimension; i++) {
      index[i]++;
      flag = (int)(!oldSet.count(index));
      index[i]--;
      if (flag) {
	activeSet.insert(std::pair<Scalar,std::vector<int> >(1.0,index));
	oldSet.erase(it1);
	break;
      }
    }
  }

  // Compute local and global error indicators for active set
  typename std::multimap<Scalar,std::vector<int> >::iterator it2;
  Scalar eta = 0.0;
  Scalar G   = 0.0;
  Teuchos::RCP<UserVector> s = integralValue.Create();
  for (it2=activeSet.begin(); it2!=activeSet.end(); it2++) {
    // Build Differential Quarature Rule
    index = it2->second;
    CubatureTensorSorted<Scalar> diffRule(0,dimension);
    build_diffRule(diffRule,index,problem_data);
    
    // Apply Rule to function
    problem_data.eval_cubature(*s,diffRule);
    
    // Update local error indicator and index set
    G  = problem_data.error_indicator(*s);
    activeSet.erase(it2);
    activeSet.insert(it2,std::pair<Scalar,std::vector<int> >(G,index));
    eta += G;
  }

  // Refine Sparse Grid
  eta = refine_grid(activeSet,oldSet,integralValue,eta,
		    dimension,rule1D,growth1D);

  // Insert New Active and Old Index sets into indexSet
  indexSet.insert(activeSet.begin(),activeSet.end());
  for (it1=oldSet.begin(); it1!=oldSet.end(); it1++) {
    index = *it1;
    indexSet.insert(std::pair<Scalar,std::vector<int> >(-1.0,index));
  }
    
  return eta;
}

// Update index set and output integral
template<class Scalar, class UserVector> 
Scalar AdaptiveSparseGrid<Scalar,UserVector>::refine_grid(
          typename std::multimap<Scalar,std::vector<int> > & activeIndex, 
      	  std::set<std::vector<int> > & oldIndex, 
	  UserVector & integralValue,
	  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;
      }
    }
    else { // Max Level Exceeded 
      //std::cout << "Maximum Level Exceeded" << std::endl;
    }
    index[k]--;
  }
  return eta;
}

// Update index set and output integral/sparse grid
template<class Scalar, class UserVector> 
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;
}

template<class Scalar, class UserVector>
void AdaptiveSparseGrid<Scalar,UserVector>::buildSparseGrid(
	      CubatureTensorSorted<Scalar> & output,
	      int dimension, int maxlevel,
	      std::vector<EIntrepidBurkardt> rule1D, 
	      std::vector<EIntrepidGrowth> growth1D,
	      bool isNormalized) {

  if (dimension == 2) {
    std::vector<int> index(dimension,0);
    for (int i=0; i<maxlevel; i++) {
      for (int j=0; j<maxlevel; j++) {
	if(i+j+dimension <= maxlevel+dimension-1) {
	  index[0] = i+1; index[1] = j+1; 
	  CubatureTensorSorted<Scalar> diffRule(0,dimension);
	  build_diffRule(diffRule,index,dimension,rule1D,growth1D,isNormalized);
	  output.update(1.0,diffRule,1.0);
	}
      }
    } 
  }
  else if (dimension == 3) {    
    std::vector<int> index(dimension,0);
    for (int i=0; i<maxlevel; i++) {
      for (int j=0; j<maxlevel; j++) {
	for (int k=0; k<maxlevel; k++) {
	  if(i+j+k+dimension <= maxlevel+dimension-1) {
	    index[0] = i+1; index[1] = j+1; index[2] = k+1;
	    CubatureTensorSorted<Scalar> diffRule(0,dimension);
	    build_diffRule(diffRule,index,dimension,rule1D,
			  growth1D,isNormalized);
	    output.update(1.0,diffRule,1.0);
	  }
	}
      }
    } 
  }
  else if (dimension == 4) {
    std::vector<int> index(dimension,0);
    for (int i=0; i<maxlevel; i++) {
      for (int j=0; j<maxlevel; j++) {
	for (int k=0; k<maxlevel; k++) {
	  for (int l=0; l<maxlevel; l++) {
	    if(i+j+k+l+dimension <= maxlevel+dimension-1) {
	      index[0] = i+1; index[1] = j+1; index[2] = k+1; index[3] = l+1;
	      CubatureTensorSorted<Scalar> diffRule(0,dimension);
	      build_diffRule(diffRule,index,dimension,rule1D,
			    growth1D,isNormalized);
	      output.update(1.0,diffRule,1.0);
	    }
	  }
	}
      } 
    }
  }
  else if (dimension == 5) {
    std::vector<int> index(dimension,0);
    for (int i=0; i<maxlevel; i++) {
      for (int j=0; j<maxlevel; j++) {
	for (int k=0; k<maxlevel; k++) {
	  for (int l=0; l<maxlevel; l++) {
	    for (int m=0; m<maxlevel; m++) {
	      if(i+j+k+l+m+dimension <= maxlevel+dimension-1) {
		index[0] = i+1; index[1] = j+1; index[2] = k+1; 
		index[3] = l+1; index[4] = m+1;
		CubatureTensorSorted<Scalar> diffRule(0,dimension);
		build_diffRule(diffRule,index,dimension,rule1D,
			      growth1D,isNormalized);
		output.update(1.0,diffRule,1.0);
	      }
	    }
	  }
	} 
      }
    }
  }
  else 
    std::cout << "Dimension Must Be Less Than 5\n";
}

} // end Intrepid namespace