amroc/clawpack/ClpIntegratorBase.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2002 Ralf Deiterding
00004 // Brandenburgische Universitaet Cottbus
00005 //
00006 // Copyright (C) 2003-2007 California Institute of Technology
00007 // Ralf Deiterding, ralf@cacr.caltech.edu
00008 
00009 #ifndef AMROC_CLP_INTEGRATORBASE_H
00010 #define AMROC_CLP_INTEGRATORBASE_H
00011 
00019 #include "Integrator.h"
00020 #include "Timing.h"
00021 #include "FixupBase.h"
00022 
00023 #include <iostream>
00024 #include <string>
00025 #include <cfloat>
00026 
00027 #define FACTOR 1.0e5
00028 
00029 #define f_init_common FORTRAN_NAME(combl, COMBL)
00030 extern "C" {
00031    void f_init_common();
00032 }
00033 
00042 template <class VectorType, class AuxVectorType, int dim>
00043 class ClpIntegratorBase : public Integrator<VectorType,dim> {
00044   typedef typename VectorType::InternalDataType DataType;
00045   typedef Integrator<VectorType,dim> base;
00046   typedef ClpIntegratorBase<VectorType,AuxVectorType,dim> self;
00047 
00048 public:
00049   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00050   typedef typename base::vec_grid_data_type vec_grid_data_type;
00051   typedef GridData<VectorType,minus_1<dim>::dim>    ld_vec_grid_data_type;
00052   typedef GridData<AuxVectorType,minus_1<dim>::dim> ld_aux_grid_data_type;
00053   typedef generic_fortran_func generic_func_type;
00054 
00055   typedef void (*check_1_func_type) ( FI(1,VectorType), BI, const INTEGER& meqn, 
00056                                       const INTEGER& mout, INTEGER& result );
00057   typedef void (*check_2_func_type) ( FI(2,VectorType), BI, const INTEGER& meqn, 
00058                                       const INTEGER& mout, INTEGER& result );
00059   typedef void (*check_3_func_type) ( FI(3,VectorType), BI, const INTEGER& meqn, 
00060                                       const INTEGER& mout, INTEGER& result );
00061 
00062   typedef void (*aux_1_func_type) ( const INTEGER& maxmx,
00063                                     const INTEGER& meqn, const INTEGER& mbc,
00064                                     const INTEGER& ibx,
00065                                     const INTEGER& mx, VectorType q[],
00066                                     DOUBLE aux[], const INTEGER& maux,
00067                                     const DOUBLE& cornx,
00068                                     const DOUBLE& dx,
00069                                     const DOUBLE& t, const DOUBLE& dt );
00070   typedef void (*aux_2_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,
00071                                     const INTEGER& meqn, const INTEGER& mbc,
00072                                     const INTEGER& ibx, const INTEGER& iby,
00073                                     const INTEGER& mx, const INTEGER& my, VectorType q[],
00074                                     DOUBLE aux[], const INTEGER& maux,
00075                                     const DOUBLE& cornx, const DOUBLE& corny,
00076                                     const DOUBLE& dx, const DOUBLE& dy,
00077                                     const DOUBLE& t, const DOUBLE& dt );
00078   typedef void (*aux_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy, const INTEGER& maxmz,
00079                                     const INTEGER& meqn, const INTEGER& mbc,
00080                                     const INTEGER& ibx, const INTEGER& iby, const INTEGER& ibz,
00081                                     const INTEGER& mx, const INTEGER& my, const INTEGER& mz, VectorType q[],
00082                                     DOUBLE aux[], const INTEGER& maux,
00083                                     const DOUBLE& cornx, const DOUBLE& corny, const DOUBLE& cornz,
00084                                     const DOUBLE& dx, const DOUBLE& dy, const DOUBLE& dz,
00085                                     const DOUBLE& t, const DOUBLE& dt );
00086 
00087   typedef void (*source_1_func_type) ( const INTEGER& maxmx,
00088                                        const INTEGER& meqn, const INTEGER& mbc,
00089                                        const INTEGER& ibx,
00090                                        const INTEGER& mx, 
00091                                        VectorType q[],
00092                                        const DOUBLE aux[], const INTEGER& maux,
00093                                        const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );
00094   typedef void (*source_2_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,
00095                                        const INTEGER& meqn, const INTEGER& mbc,
00096                                        const INTEGER& ibx, const INTEGER& iby,
00097                                        const INTEGER& mx, const INTEGER& my, 
00098                                        VectorType q[],
00099                                        const DOUBLE aux[], const INTEGER& maux,
00100                                        const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );
00101   typedef void (*source_3_func_type) ( const INTEGER& maxmx, const INTEGER& maxmy,  const INTEGER& maxmz,
00102                                        const INTEGER& meqn, const INTEGER& mbc,
00103                                        const INTEGER& ibx, const INTEGER& iby, const INTEGER& ibz,
00104                                        const INTEGER& mx, const INTEGER& my, const INTEGER& mz, 
00105                                        VectorType q[],
00106                                        const DOUBLE aux[], const INTEGER& maux,
00107                                        const DOUBLE& t, const DOUBLE& dt, const INTEGER& ibnd );
00108 
00109   ClpIntegratorBase(const int equ, const int wv) : 
00110     base(), _EqUsed(equ), _Waves(wv), f_chk(0), f_aux(0), f_src(0)
00111   { ConsInit(); }
00112   ClpIntegratorBase(const int equ, const int wv, generic_func_type check, 
00113                     generic_func_type aux, generic_func_type source) : 
00114     base(), _EqUsed(equ), _Waves(wv), f_chk(check), f_aux(aux), f_src(source)
00115   { ConsInit(); }
00116 
00117   void ConsInit() {
00118     limiter = new int[NWaves()];
00119     mthlim = new int[NWaves()];
00120     for (int i=0; i<NWaves(); i++) {
00121       limiter[i] = -1; mthlim[i] = -1;
00122     }
00123     _name = "";
00124     for (int d=0; d<7; d++) method[d] = 0;
00125   }    
00126 
00127   ~ClpIntegratorBase() {
00128     if (limiter) delete [] limiter;
00129     if (mthlim) delete [] mthlim;
00130   }
00131 
00132   //******************************************************************************
00133   // Abstract class interface
00134   //******************************************************************************
00135   virtual void ComputeRP(ld_vec_grid_data_type& LeftState, ld_aux_grid_data_type& auxl, 
00136                          ld_vec_grid_data_type& RightState, ld_aux_grid_data_type& auxr, 
00137                          ld_vec_grid_data_type& NewState, const int& direction) = 0;
00138   virtual double ComputeGrid(vec_grid_data_type& StateVec, 
00139                              vec_grid_data_type* Flux[], DataType* aux,
00140                              const double& t, const double& dt, const int& mpass) = 0;
00141   //******************************************************************************
00142 
00143   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00144     ControlDevice LocCtrl = Ctrl.getSubDevice(prefix+"ClawpackIntegrator");
00145     _name = prefix+"ClawpackIntegrator";
00146     char MethodName[16];
00147     register int d;
00148     for (d=0; d<7; d++) {
00149       sprintf(MethodName,"Method(%d)",d+1);
00150       RegisterAt(LocCtrl,MethodName,method[d]);
00151     }
00152     char LimiterName[16];
00153     for (d=0; d<NWaves(); d++) {
00154       sprintf(LimiterName,"Limiter(%d)",d+1);
00155       RegisterAt(LocCtrl,LimiterName,limiter[d]);
00156     }
00157   }
00158   virtual void register_at(ControlDevice& Ctrl) {
00159     register_at(Ctrl, "");
00160   }
00161    
00162   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00163     for (int i=0; i<NWaves(); i++) 
00164       if (limiter[i] > 0) mthlim[i] = limiter[i];
00165       else mthlim[i] = limiter[0];
00166     if (NAux()!=0 && NAux()!=AuxVectorType::Length())
00167       method[6] = 0;
00168     base::SetMaxIntegratorPasses(NMaxPass());
00169     base::_abort = (NCheck()>=10 ? DAGHTrue : DAGHFalse);
00170     f_init_common();
00171     base::SetupData(gh,ghosts);
00172   }
00173 
00174   virtual void finish() {
00175     if (limiter) {
00176       delete [] limiter;
00177       limiter = (int *) 0;
00178     }
00179     if (mthlim) {
00180       delete [] mthlim;
00181       mthlim = (int *) 0;
00182     }
00183     base::finish();
00184   } 
00185 
00186   void CalculateRP(ld_vec_grid_data_type &LeftState, BBox& bbl, double tl, double dtl, 
00187                    ld_vec_grid_data_type &RightState, BBox& bbr, double tr, double dtr, 
00188                    ld_vec_grid_data_type &NewState, const int& direction) {
00189     
00190     //------------------------------------------------------------------                      
00191     // Calculate AuxArray and source term on left and right state.
00192     //------------------------------------------------------------------                      
00193     vec_grid_data_type left(bbl, LeftState.data());
00194     vec_grid_data_type right(bbr, RightState.data());
00195     AllocError::SetTexts("ClpIntegrator","calculate_Riemann_Problem(): aux arrays");
00196     DataType* auxldata = new DataType[LeftState.bbox().size()*NAux()];
00197     DataType* auxrdata = new DataType[RightState.bbox().size()*NAux()];
00198 
00199     // Set aux arrays.
00200     if (NAux()>0) {
00201       SetAux(auxldata,left,tl,dtl);
00202       SetAux(auxrdata,right,tr,dtr);
00203     }      
00204 
00205     // Apply the source term, if necessary.
00206     START_WATCH;
00207     if (NSrc() > 0) {
00208       // Integrate source term on left and right data to the same time level.
00209       int srcbnd = 1;
00210       if (tl < tr) {
00211         SetSrc(auxldata,left,tl,tr-tl,srcbnd);
00212         tl += tr-tl;
00213       }
00214       if (tl > tr) {
00215         SetSrc(auxrdata,right,tr,tl-tr,srcbnd);
00216         tr += tl-tr;
00217       }
00218       
00219       if (NSrc() > 1) {
00220         double dt = Min(dtl,dtr);
00221         if (NSrc()==2) {
00222           SetSrc(auxldata,left,tl,dt/2.0,srcbnd);
00223           SetSrc(auxrdata,right,tr,dt/2.0,srcbnd);
00224         }
00225         if (NSrc()==3) {
00226           SetSrc(auxldata,left,tl,dt,srcbnd);
00227           SetSrc(auxrdata,right,tr,dt,srcbnd);
00228         }
00229         if (NAux()>0) {
00230           SetAux(auxldata,left,tl,dtl);
00231           SetAux(auxrdata,right,tr,dtr);
00232         }      
00233       }   
00234     }
00235     END_WATCH_SOURCE_INTEGRATION;
00236     VectorType* data_dummy;
00237     left.deallocate(data_dummy);
00238     right.deallocate(data_dummy);
00239     ld_aux_grid_data_type auxl(LeftState.bbox(), (AuxVectorType*) auxldata);
00240     ld_aux_grid_data_type auxr(RightState.bbox(), (AuxVectorType*) auxrdata);
00241 
00242     // Solve the Riemann problems. 
00243     ComputeRP(LeftState, auxl, RightState, auxr, NewState, direction);
00244 
00245     AuxVectorType* aux_data_dummy;
00246     auxl.deallocate(aux_data_dummy);
00247     auxr.deallocate(aux_data_dummy);
00248     delete [] auxldata; delete [] auxrdata;
00249   }
00250 
00251   virtual double CalculateGrid(vec_grid_data_type& NewStateVec, 
00252                                vec_grid_data_type& OldStateVec,
00253                                vec_grid_data_type* Flux[],
00254                                const int& level, const double& t, 
00255                                const double& dt, const int& mpass) {
00256 
00257     double cfl = 0.0;
00258 #ifdef DEBUG_PRINT_INTEGRATOR
00259     ( comm_service::log() << "on: " << OldStateVec.bbox() ).flush();
00260 #endif
00261 
00262     if (NMaxPass()==0 || mpass==1) 
00263       NewStateVec.copy(OldStateVec);
00264 
00265     if (dt/t >= DBL_EPSILON*FACTOR) {
00266       AllocError::SetTexts("ClpIntegrator","calculate_time_step(): aux");
00267       DataType* aux = new DataType[OldStateVec.bbox().size()*NAux()];     
00268       if (NAux() > 0) SetAux(aux,NewStateVec,t,dt);
00269 
00270 #ifdef DEBUG_AMR
00271       if (NCheck()%10>1 && (NSrc()==2 || NSrc()==3))
00272         base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before SrcInt");
00273 #endif
00274 
00275       // The source term has to be applied also on ghostcells, otherwise boundary conditions
00276       // would not be consistent. See CLAWPACK-Manual Note #5 "Source Terms".
00277       if (NSrc()>1 && (NMaxPass()==0 || mpass==1)) {
00278         START_WATCH;
00279         int srcbnd = 1;
00280         if (NSrc()==2) SetSrc(aux,NewStateVec,t,dt/2.0,srcbnd);
00281         if (NSrc()==3) SetSrc(aux,NewStateVec,t,dt,srcbnd);
00282         if (NAux()>0) SetAux(aux,NewStateVec,t,dt);
00283         if (NCheck()>=0) 
00284           if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0) {
00285             ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
00286             delete [] aux;
00287             if (cfl>0) cfl = std::max(cfl, 2.0);
00288             return cfl;
00289           }
00290         END_WATCH_SOURCE_INTEGRATION;
00291       }
00292 
00293       if (NCheck()>0 && (NMaxPass()==0 || mpass==1))
00294         base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::before f_step");
00295 
00296       // Update NewStateVec. 
00297       cfl = ComputeGrid(NewStateVec, Flux, aux, t, dt, mpass);
00298 
00299       if (NCheck()>=0) 
00300         // The negation captures cfl=NaN!
00301         if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0 || !(cfl>0)) {
00302           ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
00303           delete [] aux;
00304           if (cfl>0) cfl = std::max(cfl, 2.0);
00305           return cfl;
00306         }
00307 
00308 #ifdef DEBUG_AMR
00309       if (NCheck()%10>1 && (NMaxPass()==0 || mpass==NMaxPass()))
00310         base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after f_step");
00311 #endif
00312 
00313       // The source term is not applied to ghostcells, because they are overwritten anyway by
00314       // the application of boundary conditions in the next step.
00315       if (NSrc()>0 &&  NSrc()<3 && (NMaxPass()==0 || mpass==NMaxPass())) {
00316         START_WATCH;
00317         int srcbnd = 0;
00318         if (NAux()>0) SetAux(aux,NewStateVec,t,dt);
00319         if (NSrc()==2) SetSrc(aux, NewStateVec, t+dt/2.0, dt/2.0, srcbnd);
00320         if (NSrc()==1) SetSrc(aux, NewStateVec, t, dt, srcbnd);
00321         if (NCheck()>=0) 
00322           if (ControlGrid(NewStateVec,level,NewStateVec.bbox(),t,0)==0) {
00323             ResetGrid(NewStateVec,OldStateVec,Flux,cfl,level,t,mpass);
00324             delete [] aux;
00325             if (cfl>0) cfl = std::max(cfl, 2.0);
00326             return cfl;
00327           }
00328         END_WATCH_SOURCE_INTEGRATION
00329       }
00330 
00331 #ifdef DEBUG_AMR
00332       if (NCheck()%10>1 && (NSrc()==2 || NSrc()==1) && (mpass==0 || mpass==NMaxPass()))
00333         base::CheckGrid(NewStateVec, level, NewStateVec.bbox(), t, "CalculateGrid::after SrcInt");
00334 #endif
00335 
00336       delete [] aux;
00337     }
00338     else {
00339       VectorType zero(0.0);
00340       for (int d=0; d<dim; d++)
00341         Flux[d]->equals(zero);
00342 #ifdef DEBUG_PRINT_INTEGRATOR
00343       ( comm_service::log() <<  "dt < eps*1.0e5  " ).flush();
00344 #endif
00345     }
00346 
00347 #ifdef DEBUG_PRINT_INTEGRATOR
00348       ( comm_service::log() << "  CFL = " << cfl <<
00349         "   t = " << t << "   dt = " << dt << "\n" ).flush();
00350 #endif
00351     
00352     return cfl;
00353   }
00354 
00355   virtual void AllocGridFluxes(const BBox &bb, vec_grid_data_type**& Flux) {
00356     Flux = new vec_grid_data_type* [2*base::Dim()]; 
00357     for (register int d=0; d<2*base::Dim(); d++) 
00358       Flux[d] = new vec_grid_data_type(bb); 
00359   }
00360 
00361   virtual void DeAllocGridFluxes(vec_grid_data_type**& Flux) {
00362     if (Flux) {
00363       for (register int d=0; d<2*base::Dim(); d++) 
00364         if (Flux[d]) delete Flux[d];
00365       delete [] Flux;
00366     }
00367   }
00368     
00369   virtual void ResetGridFluxes(vec_grid_data_type**& Flux) {
00370     if (Flux) {
00371       VectorType zero(0.0);
00372       for (register int d=0; d<2*base::Dim(); d++) 
00373         if (Flux[d]) 
00374           Flux[d]->equals(zero);
00375     }
00376   }
00377 
00378   virtual int ControlGrid(vec_grid_data_type& StateVec, const int& level, 
00379                           const BBox& where, const double& time, const int verbose) {
00380     int result = 1;
00381     if (!f_chk) return result;
00382     if (dim == 1) 
00383       ((check_1_func_type) f_chk)(FA(1,StateVec),BOUNDING_BOX(where),base::NEquations(),
00384                                   verbose,result);
00385     else if (dim == 2) 
00386       ((check_2_func_type) f_chk)(FA(2,StateVec),BOUNDING_BOX(where),base::NEquations(),
00387                                   verbose,result);
00388     else if (dim == 3) 
00389       ((check_3_func_type) f_chk)(FA(3,StateVec),BOUNDING_BOX(where),base::NEquations(),
00390                                   verbose,result);
00391     return result; 
00392   }  
00393 
00394   void ResetGrid(vec_grid_data_type& StateVec, vec_grid_data_type& RecoverStateVec,
00395                  vec_grid_data_type* Flux[], double& cfl, const int& level, 
00396                  const double& t, const int& mpass) {
00397     ResetGridFluxes(Flux);
00398     StateVec.copy(RecoverStateVec);
00399     if (cfl>0) cfl = std::max(cfl, 2.0);
00400     //#ifdef DEBUG_PRINT_INTEGRATOR
00401 #ifdef DEBUG_PRINT
00402     ( comm_service::log() << "Data in " << StateVec.bbox() 
00403       << " recovered (Level " << level << "  Time = " << t << "  mpass = " 
00404       << mpass << "). Using CFL=" << cfl << "." << std::endl ).flush();
00405 #endif
00406   }
00407 
00408   void SetAux(DataType* aux, vec_grid_data_type& StateVec, const double& t, const double& dt) {
00409     if (NAux() <=0 || !f_aux) return;    
00410     DCoords dx = base::GH().worldStep(StateVec.stepsize());
00411     DCoords lbcorner = base::GH().worldCoords(StateVec.lower(),StateVec.stepsize());
00412     int mx[3], ibx[3];
00413     DimExtBBox(StateVec.bbox(), mx, ibx);
00414     if (dim == 1)
00415       ((aux_1_func_type) f_aux)(AA(1,mx),base::NEquations(),base::NGhosts(),AA(1,ibx),AA(1,mx),
00416                                 StateVec.data(),aux,NAux(),AA(1,lbcorner()),AA(1,dx()),t,dt);   
00417     else if (dim == 2)
00418       ((aux_2_func_type) f_aux)(AA(2,mx),base::NEquations(),base::NGhosts(),AA(2,ibx),AA(2,mx),
00419                                 StateVec.data(),aux,NAux(),AA(2,lbcorner()),AA(2,dx()),t,dt); 
00420     else if (dim == 3)
00421       ((aux_3_func_type) f_aux)(AA(3,mx),base::NEquations(),base::NGhosts(),AA(3,ibx),AA(3,mx),
00422                                 StateVec.data(),aux,NAux(),AA(3,lbcorner()),AA(3,dx()),t,dt); 
00423   }
00424 
00425   virtual void SetSrc(DataType* aux, vec_grid_data_type& StateVec, const double t, 
00426                        const double dt, const int& bnd) {
00427     if (NSrc() <=0 || !f_src) return;    
00428     DCoords dx = base::GH().worldStep(StateVec.stepsize());
00429     DCoords lbcorner = base::GH().worldCoords(StateVec.lower(),StateVec.stepsize());
00430     int mx[3], ibx[3];
00431     DimExtBBox(StateVec.bbox(), mx, ibx);
00432     if (dim == 1)
00433       ((source_1_func_type) f_src)(AA(1,mx),base::NEquations(),base::NGhosts(),AA(1,ibx),AA(1,mx),
00434                                    StateVec.data(),aux,NAux(),t,dt,bnd); 
00435     else if (dim == 2)
00436       ((source_2_func_type) f_src)(AA(2,mx),base::NEquations(),base::NGhosts(),AA(2,ibx),AA(2,mx),
00437                                    StateVec.data(),aux,NAux(),t,dt,bnd); 
00438     else if (dim == 3)
00439       ((source_3_func_type) f_src)(AA(3,mx),base::NEquations(),base::NGhosts(),AA(3,ibx),AA(3,mx),
00440                                    StateVec.data(),aux,NAux(),t,dt,bnd); 
00441   }
00442 
00443   void DimExtBBox(const BBox& bb, int mx[], int ibx[]) {
00444     for (int d=0; d<dim; d++) 
00445       if (bb.extents(d)-2*base::NGhosts()<1) 
00446         { mx[d] = bb.extents(d); ibx[d] = 0; }
00447       else
00448         { mx[d] = bb.extents(d)-2*base::NGhosts(); ibx[d] = 1; }
00449   }
00450 
00451   virtual int NMethodOrder() const { return (std::max(method[1],2)); }
00452   inline const int& NCheck() const { return method[3]; }
00453   inline int NSrc() const { return (method[4]%10); }
00454   inline bool DimSplit() const { return (method[2]<0); }
00455   inline int NMaxPass() const { return (method[2]<-1 && dim>1 ? dim : 0); }
00456   inline const int& NAux() const { return method[6]; }
00457   inline void SetEqUsed(const int& eq) { _EqUsed = eq; }
00458   inline const int& NEqUsed() const { return _EqUsed; }
00459   inline void SetWaves(const int& wv) { _Waves = wv; }
00460   inline const int& NWaves() const { return _Waves; }  
00461 
00462   inline void SetCheckFunc(generic_func_type check) { f_chk = check; }
00463   generic_func_type GetCheckFunc() const { return f_chk; }
00464   inline void SetAuxFunc(generic_func_type aux) { f_aux = aux; }
00465   generic_func_type GetAuxFunc() const { return f_aux; }
00466   inline void SetSrcFunc(generic_func_type src) { f_src = src; }
00467   generic_func_type GetSrcFunc() const { return f_src; }
00468 
00469   std::ostream& debug(std::ostream& out, const self& P) {
00470     out << "ClpIntegrator registered at " << P._name << "\n";
00471     out << "   Methods:";
00472     for (int i=0; i<7; i++)
00473       out << " " << P.method[i];
00474     out << "   Limiter:";
00475     for (int j=0; j<P.NWaves(); j++)
00476       out << " " << P.limiter[j];
00477     out << "\n";
00478     return out;
00479   }
00480 
00481 protected:
00482   int _EqUsed, _Waves;
00483   int method[7];
00484   int *limiter, *mthlim;
00485   std::string _name;
00486   generic_func_type f_chk, f_aux, f_src;
00487   int order;
00488 };
00489 
00490 
00491 template <class VectorType, class AuxVectorType, int dim> class ClpIntegrator;
00492 
00493 #endif