#ifndef MESH_ #define MESH_ #include #include #include #include #include #include "sutil.h" #include "zeroin.h" /////////////////////// Simplified verison of mesh1D is defined in this header //////////////////////// //*****************************************// // Classes implemented in this header file // //*****************************************// class mesh; // Basic class class mesh1D;// 1D grid class MeshJoin; // Small class for root-finding routine typedef int tint; // For position in mesh we use tint clss, which is at the moment just int. // We use class name for flexibility in implementation. //*********************************************************// // Additional nonstandard classes used in this header file // //*********************************************************// class intpar;// Simple class for exchanging data between functions and meshes when doing linear interpolation //************************************************// // Declaration of basic class for holding grids // //************************************************// class mesh{ protected: int N, N0; // size, size of the allocated memory (might be larger than size) double *om; // grid points double *delta; // contains weights for interpolation double *dh; // contains weights for integration static const int dN = 10; // when searching ordered table, searching is done first between a0 and a0+dN... protected: mesh(): N(0),N0(0),om(NULL),delta(NULL),dh(NULL){}; // constructor is made protected such that mesh can not be instantiated ~mesh(){}; // This class is used only as base class public: // OPERATORS double& operator[](int i) {Assert(i 1) { k=(khi+klo) >> 1; if (om[k] > x) khi=k; else klo=k; } return klo; } inline int mesh::find(double x) const { // If nothing is known about the point x if (x<=om[0]) return 0; if (x>=om[N-2]) return N-2; int ai0=0, ai1=N-1; return bisection(x,ai0,ai1,N); } inline int mesh::find(double x, int& dummy) const { return find(x);} inline int mesh::find_(double x, int& ai0) const { // This is most often called searching routine // It is used for searching table in increasing order if (x=N-2) return ai0; // Needs to check for the end of the table if (x=N-2) return ai1; // Again checks for the end of the table ai0 = ai1+1; // makes another step if (ai1+dNom[bi]) return bi; // Checks weather x is still in [bi:bi+1] if (bi<=0) return bi; // Checks start of the table if (x>om[bi-1]) return bi-1; // Checks weather x is in [bi-1:bi] if (bi-2<=ai0) return ai0; // If [bi-2:bi-1] is first interval (equal to [0:1]) we are done if (x>om[bi-2]) return bi-2; // Checks interbal [bi-2:bi-1] if (bi-dN>ai0){ // Bisection only between [bi-dN:ai0] if (x>om[bi-dN]){ int ai1, ai00; ai00 = bi-dN; return bisection(x,ai00,ai1,bi-1); } bi=bi-dN+1; } else bi-=1; { // If everything else failed, search everywhere between [ai0:bi] int ai1; return bisection (x,ai0,ai1,bi); } } inline int mesh::findBoth(const double x, int& ix) const { // This routine is uded when point is usually close to the previous point but not // on left or right side of it if (x==om[ix]) return ix; if (x>om[ix]) return find_(x,ix); else return _find(x,ix); } inline intpar mesh::Interp(const double x) const { int i = find(x); #ifdef _DEBUG double p = (x-om[i])*delta[i]; if ((p>1.0 || p<0) && i!=0 && i!=(N-2)) {std::cerr<<"Variable p="<