quad.h

Go to the documentation of this file.

/* =====================================================================================
 
                      The UQ Toolkit (UQTk) version 3.1.5
                          Copyright (2024) NTESS
                        https://www.sandia.gov/UQToolkit/
                        https://github.com/sandialabs/UQTk
 
     Copyright 2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
     Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government
     retains certain rights in this software.
 
     This file is part of The UQ Toolkit (UQTk)
 
     UQTk is open source software: you can redistribute it and/or modify
     it under the terms of BSD 3-Clause License
 
     UQTk is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     BSD 3 Clause License for more details.
 
     You should have received a copy of the BSD 3 Clause License
     along with UQTk. If not, see https://choosealicense.com/licenses/bsd-3-clause/.
 
     Questions? Contact the UQTk Developers at https://github.com/sandialabs/UQTk/discussions
     Sandia National Laboratories, Livermore, CA, USA
===================================================================================== */
 
#ifndef QUAD_H_SEEN
#define QUAD_H_SEEN
 
#include "Array1D.h"
#include "Array2D.h"
 
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <sstream>
using namespace std; // needed for python string conversion
 
#define QD_MAX 20000000
 
class Quad {
public:
  Quad(char *grid_type, char *fs_type,int ndim,int param,double alpha=0.0, double betta=1.0);
 
    Quad(Array1D<string>& grid_types, char *fs_type, Array1D<int>& param,Array1D<double>& alphas, Array1D<double>& bettas);
 
  Quad() {};
  ~Quad() {};
 
  void init();
 
  void SetAlpha(double alpha){this->alpha_=alpha; }
  void SetBeta(double betta){this->beta_=betta; }
 
  void SetDomain(Array1D<double>& aa, Array1D<double>& bb);
  void SetDomain(Array1D<double>& aa);
  void GetDomain(Array1D<double>& aa, Array1D<double>& bb) const {aa=aa_; bb=bb_;}
  void GetDomain(Array1D<double>& aa) const {aa=aa_;}
 
  void SetRule(Array2D<double>& q, Array1D<double>& w);
  void SetRule(Array2D<double>& q, Array1D<double>& w, Array2D<int>& ind){this->SetRule(q,w);}
  //void SetRule(Array2D<double>& q, Array1D<double>& w, Array2D<int>& ind, int param);
 
  void SetRule();
 
  void GetRule(Array2D<double>& q, Array1D<double>& w);
  void GetRule(Array2D<double>& q, Array1D<double>& w, Array2D<int>& ind){this->GetRule(q,w);}
 
  void SetQdpts(Array2D<double>& q){  rule_.qdpts=q; return;}
  void SetWghts(Array1D<double>& w){  rule_.wghts=w; return;}
  // void SetIndices(Array2D<int>& ind){  rule_.indices=ind; return;}
 
  void GetQdpts(Array2D<double>& q){ q=rule_.qdpts;  return;}
  void GetWghts(Array1D<double>& w){  w=rule_.wghts; return;}
  //void GetIndices(Array2D<int>& ind){  ind=rule_.indices; return;}
 
  void SetLevel(int param) {nlevel_=param; return;}
 
  void nextLevel();
 
  int GetNQ() {return rule_.qdpts.XSize(); }
 
  void SetVerbosity(int verbosity) { quadverbose_ = verbosity; }
 
 
 
 private:
 
  Quad(const Quad &) {};
 
  Array1D<double> aa_;
  Array1D<double> bb_;
 
  int quadverbose_;
 
  double alpha_;
  double beta_;
  Array1D<double> alphas_;
  Array1D<double> betas_;
 
  typedef struct
  {
    Array2D<double> qdpts;
    Array1D<double> wghts;
  } QuadRule;
 
  QuadRule rule_;
 
  int ndim_;
 
  int nlevel_;
 
  int maxlevel_;
  Array1D<int> param_;
 
  Array2D<int> npts_all, npts_1_all;
  Array2D<QuadRule> qr_all, qr_1_all;
 
  void MultiplyTwoRules(QuadRule *rule1,QuadRule *rule2,QuadRule *rule_prod);
  void MultiplyManyRules(int nrules, QuadRule *rules, QuadRule *rule_prod);
  void MultiplyManyRules_(int nrules, QuadRule *rules, QuadRule *rule_prod);
  void AddTwoRules(QuadRule *rule1,QuadRule *rule2,QuadRule *rule_sum);
  void SubtractTwoRules(QuadRule *rule1,QuadRule *rule2,QuadRule *rule_sum);
 
  void create1DRule(string gridtype,Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
 
  void create1DRule_CC(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_LU(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_HG(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr);
  void create1DRule_NC(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_NCO(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_CCO(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_JB(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_LG(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr);
  void create1DRule_SW(Array1D<double>& qdpts,Array1D<double>& wghts,int ngr);
  void create1DRule_pdf(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
  void create1DRule_GP3(Array1D<double>& qdpts,Array1D<double>& wghts, int ngr, double a, double b);
 
  void getMultiIndexLevel(Array2D<int>& multiIndexLevel, int level, int ndim);
 
  int growth_rule_;
  Array1D<int> growth_rules_;
 
  string grid_type_;
  Array1D<string> grid_types_;
 
  string fs_type_;
 
  void compressRule(QuadRule *rule);
 
 
 
};
#endif /* QUAD_H_SEEN */