AbstractConvergenceTester.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #ifndef ABSTRACTCONVERGENCETESTER_HPP_
00031 #define ABSTRACTCONVERGENCETESTER_HPP_
00032 
00033 #include "BidomainProblem.hpp"
00034 #include "MonodomainProblem.hpp"
00035 
00036 #include <petscvec.h>
00037 #include <vector>
00038 #include <iostream>
00039 #include <fstream>
00040 #include <cmath>
00041 
00042 #include "AbstractTetrahedralMesh.hpp"
00043 #include "TrianglesMeshReader.hpp"
00044 #include "AbstractCardiacCellFactory.hpp"
00045 #include "OutputFileHandler.hpp"
00046 #include "TrianglesMeshWriter.hpp"
00047 #include "PropagationPropertiesCalculator.hpp"
00048 #include "Hdf5DataReader.hpp"
00049 #include "GeneralPlaneStimulusCellFactory.hpp"
00050 #include "CuboidMeshConstructor.hpp"
00051 #include "ZeroStimulusCellFactory.hpp"
00052 #include "SimpleStimulus.hpp"
00053 #include "ConstBoundaryCondition.hpp"
00054 #include "StimulusBoundaryCondition.hpp"
00055 
00056 #include "Warnings.hpp"
00057 
00058 typedef enum StimulusType_
00059 {
00060     PLANE=0,
00061     QUARTER,
00062     NEUMANN
00063 } StimulusType;
00064 
00069 template <class CELL, unsigned DIM>
00070 class RampedQuarterStimulusCellFactory : public AbstractCardiacCellFactory<DIM>
00071 {
00072 private:
00074     std::vector< boost::shared_ptr<SimpleStimulus> > mpStimuli;
00076     double mMeshWidth;
00078     double mStepSize;
00082     unsigned mLevels;
00083 public:
00084 
00091     RampedQuarterStimulusCellFactory(double meshWidth, unsigned numElemAcross, double fullStim=-1e7)
00092         : AbstractCardiacCellFactory<DIM>(),
00093           mMeshWidth(meshWidth),
00094           mStepSize(meshWidth/numElemAcross),
00095           mLevels(numElemAcross/4)
00096     {
00097         assert(numElemAcross%4 == 0); //numElemAcross is supposed to be a multiple of 4
00098 
00099         for (unsigned level=0; level<mLevels; level++)
00100         {
00101             double this_stim = fullStim - (level*fullStim)/mLevels;
00102             //this_stim is full_stim at the zero level and would be zero at level=mLevels
00103             mpStimuli.push_back((boost::shared_ptr<SimpleStimulus>)new SimpleStimulus(this_stim, 0.5));
00104         }
00105     }
00106 
00107 
00113     AbstractCardiacCell* CreateCardiacCellForTissueNode(unsigned node)
00114     {
00115         double x = this->GetMesh()->GetNode(node)->GetPoint()[0];
00116         double d_level = x/mStepSize;
00117         unsigned level = (unsigned) d_level;
00118         assert(fabs(level-d_level) < DBL_MAX); //x ought to really be a multiple of the step size
00119 
00120         if (level < mLevels)
00121         {
00122             return new CELL(this->mpSolver, this->mpStimuli[level]);
00123         }
00124         else
00125         {
00126             return new CELL(this->mpSolver, this->mpZeroStimulus);
00127         }
00128     }
00129 };
00130 
00131 
00136 class AbstractUntemplatedConvergenceTester
00137 {
00138 protected:
00140     double mMeshWidth;
00141 public:
00143     double OdeTimeStep;
00145     double PdeTimeStep;
00150     unsigned MeshNum;
00151     double RelativeConvergenceCriterion; 
00152     double LastDifference; 
00153     double Apd90FirstQn; 
00154     double Apd90ThirdQn; 
00155     double ConductionVelocity;  
00159     bool PopulatedResult;
00163     bool FixedResult;
00167     bool UseAbsoluteStimulus;
00172     double AbsoluteStimulus;
00173     bool SimulateFullActionPotential; 
00174     bool Converged; 
00175     StimulusType Stimulus; 
00176     double NeumannStimulus; 
00178     AbstractUntemplatedConvergenceTester();
00179 
00185     virtual void Converge(std::string nameOfTest)=0;
00186 
00187     virtual ~AbstractUntemplatedConvergenceTester();
00188 };
00189 
00193 template<class CARDIAC_PROBLEM, unsigned DIM>
00194 void SetConductivities(CARDIAC_PROBLEM& rCardiacProblem);
00195 
00196 template<unsigned DIM>
00197 void SetConductivities(BidomainProblem<DIM>& rProblem)
00198 {
00199     c_vector<double, DIM> conductivities;
00200     for (unsigned i=0; i<DIM; i++)
00201     {
00202         conductivities[i] = 1.75;
00203     }
00204     HeartConfig::Instance()->SetIntracellularConductivities(conductivities);
00205 
00206     for (unsigned i=0; i<DIM; i++)
00207     {
00208         conductivities[i] = 7.0;
00209     }
00210     HeartConfig::Instance()->SetExtracellularConductivities(conductivities);
00211 }
00212 
00213 template<unsigned DIM>
00214 void SetConductivities(MonodomainProblem<DIM>& rProblem)
00215 {
00216     c_vector<double, DIM> conductivities;
00217     for (unsigned i=0; i<DIM; i++)
00218     {
00219         conductivities[i] = 1.75;
00220     }
00221     HeartConfig::Instance()->SetIntracellularConductivities(conductivities);
00222 }
00223 
00228 template<class CELL, class CARDIAC_PROBLEM, unsigned DIM, unsigned PROBLEM_DIM>
00229 class AbstractConvergenceTester : public AbstractUntemplatedConvergenceTester
00230 {
00231 public:
00236     void Converge(std::string nameOfTest)
00237     {
00238         std::cout << "=========================== Beginning Test...==================================\n";
00239         // Create the meshes on which the test will be based
00240         const std::string mesh_dir = "ConvergenceMesh";
00241         OutputFileHandler output_file_handler(mesh_dir);
00242         ReplicatableVector voltage_replicated;
00243 
00244         unsigned file_num=0;
00245 
00246         // Create a file for storing conduction velocity and AP data and write the header
00247         OutputFileHandler conv_info_handler("ConvergencePlots", false);
00248         out_stream p_conv_info_file;
00249         if (PetscTools::AmMaster())
00250         {
00251             p_conv_info_file = conv_info_handler.OpenOutputFile(nameOfTest+"_info.csv");
00252             (*p_conv_info_file) << "#Abcisa\t"
00253                                 << "l2-norm-full\t"
00254                                 << "l2-norm-onset\t"
00255                                 << "Max absolute err\t"
00256                                 << "APD90_1st_quad\t"
00257                                 << "APD90_3rd_quad\t"
00258                                 << "Conduction velocity (relative diffs)" << std::endl;
00259         }
00260         SetInitialConvergenceParameters();
00261 
00262         double prev_apd90_first_qn=0.0;
00263         double prev_apd90_third_qn=0.0;
00264         double prev_cond_velocity=0.0;
00265         std::vector<double> prev_voltage;
00266         std::vector<double> prev_times;
00267         PopulateStandardResult(prev_voltage, prev_times);
00268 
00269         do
00270         {
00271             CuboidMeshConstructor<DIM> constructor;
00272 
00273             //If the printing time step is too fine, then simulations become I/O bound without much improvement in accuracy
00274             double printing_step = this->PdeTimeStep;
00275 #define COVERAGE_IGNORE
00276             while (printing_step < 1.0e-4)
00277             {
00278                 printing_step *= 2.0;
00279                 std::cout<<"Warning: PrintingTimeStep increased\n";
00280             }
00281 #undef COVERAGE_IGNORE
00282             HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(this->OdeTimeStep, this->PdeTimeStep, printing_step);
00283 #define COVERAGE_IGNORE
00284             if (SimulateFullActionPotential)
00285             {
00286                 HeartConfig::Instance()->SetSimulationDuration(400.0);
00287             }
00288             else
00289             {
00290                 HeartConfig::Instance()->SetSimulationDuration(8.0);
00291             }
00292 #undef COVERAGE_IGNORE
00293             HeartConfig::Instance()->SetOutputDirectory("Convergence");
00294             HeartConfig::Instance()->SetOutputFilenamePrefix("Results");
00295 
00296             DistributedTetrahedralMesh<DIM, DIM> mesh;
00297             constructor.Construct(mesh, this->MeshNum, mMeshWidth);
00298 
00299             unsigned num_ele_across = (unsigned) SmallPow(2, this->MeshNum+2); // number of elements in each dimension
00300 
00301             AbstractCardiacCellFactory<DIM>* p_cell_factory=NULL;
00302 
00303             switch (this->Stimulus)
00304             {
00305                 case NEUMANN:
00306                 {
00307                     p_cell_factory = new ZeroStimulusCellFactory<CELL, DIM>();
00308                     break;
00309                 }
00310                 case PLANE:
00311                 {
00312                     if (this->UseAbsoluteStimulus)
00313                     {
00314                         #define COVERAGE_IGNORE
00315                         p_cell_factory = new GeneralPlaneStimulusCellFactory<CELL, DIM>(0, this->AbsoluteStimulus, true);
00316                         #undef COVERAGE_IGNORE
00317                     }
00318                     else
00319                     {
00320                         p_cell_factory = new GeneralPlaneStimulusCellFactory<CELL, DIM>(num_ele_across, constructor.GetWidth(), false, this->AbsoluteStimulus);
00321                     }
00322                     break;
00323                 }
00324                 case QUARTER:
00325                 {
00327                     p_cell_factory = new RampedQuarterStimulusCellFactory<CELL, DIM>(constructor.GetWidth(), num_ele_across, this->AbsoluteStimulus/10.0);
00328                     break;
00329                 }
00330             }
00331 
00332 
00333             CARDIAC_PROBLEM cardiac_problem(p_cell_factory);
00334             cardiac_problem.SetMesh(&mesh);
00335 
00336             // Calculate positions of nodes 1/4 and 3/4 through the mesh
00337             unsigned third_quadrant_node;
00338             unsigned first_quadrant_node;
00339             switch(DIM)
00340             {
00341                 case 1:
00342                 {
00343                     first_quadrant_node = (unsigned) (0.25*constructor.GetNumElements());
00344                     third_quadrant_node = (unsigned) (0.75*constructor.GetNumElements());
00345                     break;
00346                 }
00347                 case 2:
00348                 {
00349                     unsigned n= (unsigned) SmallPow (2, this->MeshNum+2);
00350                     first_quadrant_node =   (n+1)*(n/2)+  n/4 ;
00351                     third_quadrant_node =   (n+1)*(n/2)+3*n/4 ;
00352                     break;
00353                 }
00354                 case 3:
00355                 {
00356                     const unsigned first_quadrant_nodes_3d[5]={61, 362, 2452, 17960, 137296};
00357                     const unsigned third_quadrant_nodes_3d[5]={63, 366, 2460, 17976, 137328};
00358                     assert(this->PdeTimeStep<5);
00359                     first_quadrant_node = first_quadrant_nodes_3d[this->MeshNum];
00360                     third_quadrant_node = third_quadrant_nodes_3d[this->MeshNum];
00361                     break;
00362                 }
00363 
00364                 default:
00365                     NEVER_REACHED;
00366             }
00367 
00368             double mesh_width=constructor.GetWidth();
00369 
00370             // We only need the output of these two nodes
00371             std::vector<unsigned> nodes_to_be_output;
00372             nodes_to_be_output.push_back(first_quadrant_node);
00373             nodes_to_be_output.push_back(third_quadrant_node);
00374             cardiac_problem.SetOutputNodes(nodes_to_be_output);
00375 
00376             // The results of the tests were originally obtained with the following conductivity
00377             // values. After implementing fibre orientation the defaults changed. Here we set
00378             // the former ones to be used.
00379             SetConductivities(cardiac_problem);
00380 
00381             cardiac_problem.Initialise();
00382 
00383             boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,PROBLEM_DIM> > p_bcc(new BoundaryConditionsContainer<DIM,DIM,PROBLEM_DIM>);
00384             SimpleStimulus stim(NeumannStimulus, 0.5);
00385             if (Stimulus==NEUMANN)
00386             {
00387 
00388                 StimulusBoundaryCondition<DIM>* p_bc_stim = new StimulusBoundaryCondition<DIM>(&stim);
00389 
00390                 // get mesh
00391                 AbstractTetrahedralMesh<DIM, DIM> &r_mesh = cardiac_problem.rGetMesh();
00392                 // loop over boundary elements
00393                 typename AbstractTetrahedralMesh<DIM, DIM>::BoundaryElementIterator iter;
00394                 iter = r_mesh.GetBoundaryElementIteratorBegin();
00395                 while (iter != r_mesh.GetBoundaryElementIteratorEnd())
00396                 {
00397                     double x = ((*iter)->CalculateCentroid())[0];
00398                     if (x*x<=1e-10)
00399                     {
00400                         p_bcc->AddNeumannBoundaryCondition(*iter, p_bc_stim);
00401                     }
00402                     iter++;
00403                 }
00404                 // pass the bcc to the problem
00405                 cardiac_problem.SetBoundaryConditionsContainer(p_bcc);
00406             }
00407 
00408             DisplayRun();
00409             double time_before=MPI_Wtime();
00411             //cardiac_problem.SetWriteInfo();
00412 
00413             try
00414             {
00415                 cardiac_problem.Solve();
00416             }
00417             catch (Exception e)
00418             {
00419                 WARNING("This run threw an exception.  Check convergence results\n");
00420                 std::cout << e.GetMessage() << std::endl;
00421             }
00422 
00423             std::cout << "Time to solve = " << MPI_Wtime()-time_before << " seconds\n";
00424 
00425             OutputFileHandler results_handler("Convergence", false);
00426             Hdf5DataReader results_reader = cardiac_problem.GetDataReader();
00427 
00428             {
00429                 std::vector<double> transmembrane_potential = results_reader.GetVariableOverTime("V", third_quadrant_node);
00430                 std::vector<double> time_series = results_reader.GetUnlimitedDimensionValues();
00431 
00432                 OutputFileHandler plot_file_handler("ConvergencePlots", false);
00433                 if (PetscTools::AmMaster())
00434                 {
00435                     // Write out the time series for the node at third quadrant
00436                     {
00437                         std::stringstream plot_file_name_stream;
00438                         plot_file_name_stream<< nameOfTest << "_Third_quadrant_node_run_"<< file_num << ".csv";
00439                         out_stream p_plot_file = plot_file_handler.OpenOutputFile(plot_file_name_stream.str());
00440                         for (unsigned data_point = 0; data_point<time_series.size(); data_point++)
00441                         {
00442                             (*p_plot_file) << time_series[data_point] << "\t" << transmembrane_potential[data_point] << "\n";
00443                         }
00444                         p_plot_file->close();
00445                     }
00446                     // Write time series for first quadrant node
00447                     {
00448                         std::vector<double> transmembrane_potential_1qd=results_reader.GetVariableOverTime("V", first_quadrant_node);
00449                         std::vector<double> time_series_1qd = results_reader.GetUnlimitedDimensionValues();
00450                         std::stringstream plot_file_name_stream;
00451                         plot_file_name_stream<< nameOfTest << "_First_quadrant_node_run_"<< file_num << ".csv";
00452                         out_stream p_plot_file = plot_file_handler.OpenOutputFile(plot_file_name_stream.str());
00453                         for (unsigned data_point = 0; data_point<time_series.size(); data_point++)
00454                         {
00455                             (*p_plot_file) << time_series_1qd[data_point] << "\t" << transmembrane_potential_1qd[data_point] << "\n";
00456                         }
00457                         p_plot_file->close();
00458                     }
00459                 }
00460                 // calculate conduction velocity and APD90 error
00461                 PropagationPropertiesCalculator ppc(&results_reader);
00462 
00463                 try
00464                 {
00465                     #define COVERAGE_IGNORE
00466                     if (SimulateFullActionPotential)
00467                     {
00468                         Apd90FirstQn = ppc.CalculateActionPotentialDuration(90.0, first_quadrant_node);
00469                         Apd90ThirdQn = ppc.CalculateActionPotentialDuration(90.0, third_quadrant_node);
00470                     }
00471                     #undef COVERAGE_IGNORE
00472                     ConductionVelocity  = ppc.CalculateConductionVelocity(first_quadrant_node,third_quadrant_node,0.5*mesh_width);
00473                 }
00474                 catch (Exception e)
00475                 {
00476                     #define COVERAGE_IGNORE
00477                     std::cout << "Warning - this run threw an exception in calculating propagation.  Check convergence results\n";
00478                     std::cout << e.GetMessage() << std::endl;
00479                     #undef COVERAGE_IGNORE
00480                 }
00481                 double cond_velocity_error = 1e10;
00482                 double apd90_first_qn_error = 1e10;
00483                 double apd90_third_qn_error = 1e10;
00484 
00485                 if (this->PopulatedResult)
00486                 {
00487                     if (prev_cond_velocity != 0.0)
00488                     {
00489                         cond_velocity_error = fabs(ConductionVelocity - prev_cond_velocity) / prev_cond_velocity;
00490                     }
00491 #define COVERAGE_IGNORE
00492                     if (prev_apd90_first_qn != 0.0)
00493                     {
00494                         apd90_first_qn_error = fabs(Apd90FirstQn - prev_apd90_first_qn) / prev_apd90_first_qn;
00495                     }
00496                     if (prev_apd90_third_qn != 0.0)
00497                     {
00498                         apd90_third_qn_error = fabs(Apd90ThirdQn - prev_apd90_third_qn) / prev_apd90_third_qn;
00499                     }
00500                     if (apd90_first_qn_error == 0.0)
00501                     {
00502                         apd90_first_qn_error = DBL_EPSILON; //Avoid log zero on plot
00503                     }
00504                     if (apd90_third_qn_error == 0.0)
00505                     {
00506                         apd90_third_qn_error = DBL_EPSILON; //Avoid log zero on plot
00507                     }
00508 #undef COVERAGE_IGNORE
00509                     if (cond_velocity_error == 0.0)
00510                     {
00511                         cond_velocity_error = DBL_EPSILON; //Avoid log zero on plot
00512                     }
00513                 }
00514 
00515                 prev_cond_velocity = ConductionVelocity;
00516                 prev_apd90_first_qn = Apd90FirstQn;
00517                 prev_apd90_third_qn = Apd90ThirdQn;
00518 
00519                 // calculate l2norm
00520                 double max_abs_error = 0;
00521                 double sum_sq_abs_error =0;
00522                 double sum_sq_prev_voltage = 0;
00523                 double sum_sq_abs_error_full =0;
00524                 double sum_sq_prev_voltage_full = 0;
00525                 if (this->PopulatedResult)
00526                 {
00527                     //If the PDE step is varying then we'll have twice as much data now as we use to have
00528 
00529                     unsigned time_factor=(time_series.size()-1) / (prev_times.size()-1);
00530                     assert (time_factor == 1 || time_factor == 2 || time_factor == 8);
00531                     //Iterate over the shorter time series data
00532                     for (unsigned data_point = 0; data_point<prev_times.size(); data_point++)
00533                     {
00534                         unsigned this_data_point=time_factor*data_point;
00535 
00536                         assert(time_series[this_data_point] == prev_times[data_point]);
00537                         double abs_error = fabs(transmembrane_potential[this_data_point]-prev_voltage[data_point]);
00538                         max_abs_error = (abs_error > max_abs_error) ? abs_error : max_abs_error;
00539                         //Only do resolve the upstroke...
00540                         sum_sq_abs_error_full += abs_error*abs_error;
00541                         sum_sq_prev_voltage_full += prev_voltage[data_point] * prev_voltage[data_point];
00542                         if (time_series[this_data_point] <= 8.0)
00543                         {
00544                             //In most regular cases we always do this, since the simulation stops at ms
00545                             sum_sq_abs_error += abs_error*abs_error;
00546                             sum_sq_prev_voltage += prev_voltage[data_point] * prev_voltage[data_point];
00547                         }
00548                     }
00549 
00550                 }
00551                 if (!this->PopulatedResult || !FixedResult)
00552                 {
00553                     prev_voltage = transmembrane_potential;
00554                     prev_times = time_series;
00555                 }
00556 
00557                 if (this->PopulatedResult)
00558                 {
00559 
00560                     if (PetscTools::AmMaster())
00561                     {
00562                         (*p_conv_info_file) << std::setprecision(8)
00563                                             << Abscissa() << "\t"
00564                                             << sqrt(sum_sq_abs_error_full/sum_sq_prev_voltage_full) << "\t"
00565                                             << sqrt(sum_sq_abs_error/sum_sq_prev_voltage) << "\t"
00566                                             << max_abs_error << "\t"
00567                                             << Apd90FirstQn <<" "<< apd90_first_qn_error <<""<< "\t"
00568                                             << Apd90ThirdQn <<" "<< apd90_third_qn_error <<""<< "\t"
00569                                             << ConductionVelocity <<" "<< cond_velocity_error  <<""<< std::endl;
00570                     }
00571                     // convergence criterion
00572                     this->Converged = sum_sq_abs_error/sum_sq_prev_voltage<this->RelativeConvergenceCriterion;
00573                     this->LastDifference=sum_sq_abs_error/sum_sq_prev_voltage;
00574 #define COVERAGE_IGNORE
00575                     if (time_series.size() == 1u)
00576                     {
00577                         std::cout << "Failed after successful convergence - give up this convergence test\n";
00578                         break;
00579                     }
00580 #undef COVERAGE_IGNORE
00581 
00582                 }
00583 
00584                 if (!this->PopulatedResult)
00585                 {
00586                     this->PopulatedResult=true;
00587 
00588                 }
00589             }
00590 
00591             // Get ready for the next test by halving the time step
00592             if (!this->Converged)
00593             {
00594                 UpdateConvergenceParameters();
00595                 file_num++;
00596             }
00597             delete p_cell_factory;
00598         }
00599         while (!GiveUpConvergence() && !this->Converged);
00600 
00601 
00602         if (PetscTools::AmMaster())
00603         {
00604             p_conv_info_file->close();
00605 
00606             std::cout << "Results: " << std::endl;
00607             MPIABORTIFNON0(system, "cat " + conv_info_handler.GetOutputDirectoryFullPath() + nameOfTest + "_info.csv");
00608         }
00609 
00610     }
00611 
00612     void DisplayRun()
00613     {
00614         unsigned num_ele_across = (unsigned) SmallPow(2, this->MeshNum+2);// number of elements in each dimension
00615         double scaling = mMeshWidth/(double) num_ele_across;
00616 
00617         std::cout<<"================================================================================"<<std::endl;
00618         std::cout<<"Solving in "<<DIM<<"D\t";
00619         std::cout<<"Space step "<< scaling << " cm (mesh " << this->MeshNum << ")" << "\n";
00620         std::cout<<"PDE step "<<this->PdeTimeStep<<" ms"<<"\t";
00621         std::cout<<"ODE step "<<this->OdeTimeStep<<" ms"<<"\t";
00622         if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
00623         {
00624             std::cout<<"KSP absolute "<<HeartConfig::Instance()->GetAbsoluteTolerance()<<"\t";
00625         }
00626         else
00627         {
00628             std::cout<<"KSP relative "<<HeartConfig::Instance()->GetRelativeTolerance()<<"\t";
00629         }
00630         switch (this->Stimulus)
00631         {
00632             case PLANE:
00633             std::cout<<"Stimulus = Plane\n";
00634             break;
00635 
00636             case QUARTER:
00637             std::cout<<"Stimulus = Ramped first quarter\n";
00638             break;
00639 
00640             case NEUMANN:
00641             std::cout<<"Stimulus = Neumann\n";
00642             break;
00643 
00644         }
00645         EXPECT0(system, "date");//To keep track of what Nightly things are doing
00648         if (this->UseAbsoluteStimulus)
00649         {
00650             #define COVERAGE_IGNORE
00651             std::cout<<"Using absolute stimulus of "<<this->AbsoluteStimulus<<std::endl;
00652             #undef COVERAGE_IGNORE
00653         }
00654         std::cout << std::flush;
00655         //HeartEventHandler::Headings();
00656         //HeartEventHandler::Report();
00657 
00658     }
00659 
00660 public:
00661     virtual ~AbstractConvergenceTester() {}
00662 
00664     virtual void SetInitialConvergenceParameters()=0;
00666     virtual void UpdateConvergenceParameters()=0;
00668     virtual bool GiveUpConvergence()=0;
00670     virtual double Abscissa()=0;
00677     virtual void PopulateStandardResult(std::vector<double> &result, std::vector<double> &times)
00678     {
00679         assert(this->PopulatedResult==false);
00680         assert(result.size()==0);
00681         assert(times.size()==0);
00682     }
00683 
00687     bool IsConverged()
00688     {
00689         return Converged;
00690     }
00691 
00695     void SetMeshWidth(double meshWidth)
00696     {
00697         mMeshWidth=meshWidth;
00698     }
00699 };
00700 
00701 #endif /*ABSTRACTCONVERGENCETESTER_HPP_*/

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5