AbstractCardiacProblem.cpp

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 #include "AbstractCardiacProblem.hpp"
00030 
00031 #include "GenericMeshReader.hpp"
00032 #include "Exception.hpp"
00033 #include "HeartConfig.hpp"
00034 #include "HeartEventHandler.hpp"
00035 #include "TimeStepper.hpp"
00036 #include "PetscTools.hpp"
00037 #include "DistributedVector.hpp"
00038 #include "ProgressReporter.hpp"
00039 #include "LinearSystem.hpp"
00040 #include "PostProcessingWriter.hpp"
00041 #include "Hdf5ToMeshalyzerConverter.hpp"
00042 #include "Hdf5ToCmguiConverter.hpp"
00043 #include "Hdf5ToVtkConverter.hpp"
00044 
00045 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00046 AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AbstractCardiacProblem(
00047             AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory)
00048     : mMeshFilename(""),               // i.e. undefined
00049       mUseMatrixBasedRhsAssembly(true),
00050       mAllocatedMemoryForMesh(false),
00051       mWriteInfo(false),
00052       mPrintOutput(true),
00053       mpCardiacTissue(NULL),
00054       mpSolver(NULL),
00055       mpCellFactory(pCellFactory),
00056       mpMesh(NULL),
00057       mSolution(NULL),
00058       mCurrentTime(0.0),
00059       mpTimeAdaptivityController(NULL),
00060       mpWriter(NULL)
00061 {
00062     assert(mNodesToOutput.empty());
00063     if (!mpCellFactory)
00064     {
00065         EXCEPTION("AbstractCardiacProblem: Please supply a cell factory pointer to your cardiac problem constructor.");
00066     }
00067     HeartEventHandler::BeginEvent(HeartEventHandler::EVERYTHING);
00068 }
00069 
00070 
00071 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00072 AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AbstractCardiacProblem()
00073     // it doesn't really matter what we initialise these to, as they'll be overwritten by
00074     // the serialization methods
00075     : mMeshFilename(""),
00076       mUseMatrixBasedRhsAssembly(true),
00077       mAllocatedMemoryForMesh(false), // Handled by AbstractCardiacTissue
00078       mWriteInfo(false),
00079       mPrintOutput(true),
00080       mVoltageColumnId(UINT_MAX),
00081       mTimeColumnId(UINT_MAX),
00082       mNodeColumnId(UINT_MAX),
00083       mpCardiacTissue(NULL),
00084       mpSolver(NULL),
00085       mpCellFactory(NULL),
00086       mpMesh(NULL),
00087       mSolution(NULL),
00088       mCurrentTime(0.0),
00089       mpTimeAdaptivityController(NULL),
00090       mpWriter(NULL)
00091 {
00092 }
00093 
00094 
00095 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00096 AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::~AbstractCardiacProblem()
00097 {
00098     delete mpCardiacTissue;
00099     if (mSolution)
00100     {
00101         VecDestroy(mSolution);
00102     }
00103 
00104     if (mAllocatedMemoryForMesh)
00105     {
00106         delete mpMesh;
00107     }
00108 }
00109 
00110 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00111 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::Initialise()
00112 {
00113     HeartEventHandler::BeginEvent(HeartEventHandler::READ_MESH);
00114     if (mpMesh==NULL)
00115     {
00116         // If no mesh has been passed, we get it from the configuration file
00117         try
00118         {
00119             if (HeartConfig::Instance()->GetLoadMesh())
00120             {
00121                 mpMesh = new DistributedTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>(HeartConfig::Instance()->GetMeshPartitioning());
00122                 GenericMeshReader<ELEMENT_DIM, SPACE_DIM> mesh_reader(HeartConfig::Instance()->GetMeshName());
00123                 mpMesh->ConstructFromMeshReader(mesh_reader);
00124             }
00125             else if (HeartConfig::Instance()->GetCreateMesh())
00126             {
00127                 mpMesh = new DistributedTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>(HeartConfig::Instance()->GetMeshPartitioning());
00128                 assert(HeartConfig::Instance()->GetSpaceDimension()==SPACE_DIM);
00129                 double inter_node_space = HeartConfig::Instance()->GetInterNodeSpace();
00130 
00131                 switch(HeartConfig::Instance()->GetSpaceDimension())
00132                 {
00133                     case 1:
00134                     {
00135                         c_vector<double, 1> fibre_length;
00136                         HeartConfig::Instance()->GetFibreLength(fibre_length);
00137                         mpMesh->ConstructRegularSlabMesh(inter_node_space, fibre_length[0]);
00138                         break;
00139                     }
00140                     case 2:
00141                     {
00142                         c_vector<double, 2> sheet_dimensions; //cm
00143                         HeartConfig::Instance()->GetSheetDimensions(sheet_dimensions);
00144                         mpMesh->ConstructRegularSlabMesh(inter_node_space, sheet_dimensions[0], sheet_dimensions[1]);
00145                         break;
00146                     }
00147                     case 3:
00148                     {
00149                         c_vector<double, 3> slab_dimensions; //cm
00150                         HeartConfig::Instance()->GetSlabDimensions(slab_dimensions);
00151                         mpMesh->ConstructRegularSlabMesh(inter_node_space, slab_dimensions[0], slab_dimensions[1], slab_dimensions[2]);
00152                         break;
00153                     }
00154                     default:
00155                         NEVER_REACHED;
00156                 }
00157             }
00158             else
00159             {
00160                 NEVER_REACHED;
00161             }
00162 
00163             mAllocatedMemoryForMesh = true;
00164         }
00165         catch (Exception& e)
00166         {
00167             EXCEPTION(std::string("No mesh given: define it in XML parameters file or call SetMesh()\n") + e.GetShortMessage());
00168         }
00169     }
00170     mpCellFactory->SetMesh( mpMesh );
00171     HeartEventHandler::EndEvent(HeartEventHandler::READ_MESH);
00172 
00173     HeartEventHandler::BeginEvent(HeartEventHandler::INITIALISE);    
00174     
00175     // if the user requested transmural stuff, we fill in the mCellHeterogeneityAreas here.
00176     if (HeartConfig::Instance()->AreCellularTransmuralHeterogeneitiesRequested())
00177     {
00178         mpCellFactory->FillInCellularTransmuralAreas();
00179     }
00180 
00181     delete mpCardiacTissue; // In case we're called twice
00182     mpCardiacTissue = CreateCardiacTissue();
00183 
00184     HeartEventHandler::EndEvent(HeartEventHandler::INITIALISE);
00185 
00186     // Delete any previous solution, so we get a fresh initial condition
00187     if (mSolution)
00188     {
00189         HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00190         VecDestroy(mSolution);
00191         mSolution = NULL;
00192         HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00193     }
00194 
00195     // Always start at time zero
00196     mCurrentTime = 0.0;
00197 
00198     //For Bidomain with bath, this is where we set up the electrodes
00199 
00200     SetElectrodes();
00201 }
00202 
00203 
00204 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00205 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetBoundaryConditionsContainer(boost::shared_ptr<BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> > pBcc)
00206 {
00207     this->mpBoundaryConditionsContainer = pBcc;
00208 }
00209 
00210 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00211 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::PreSolveChecks()
00212 {
00213     if ( mpCardiacTissue == NULL ) // if tissue is NULL, Initialise() probably hasn't been called
00214     {
00215         EXCEPTION("Cardiac tissue is null, Initialise() probably hasn't been called");
00216     }
00217     if ( HeartConfig::Instance()->GetSimulationDuration() <= mCurrentTime)
00218     {
00219         EXCEPTION("End time should be in the future");
00220     }
00221     if (mPrintOutput==true)
00222     {
00223         if( (HeartConfig::Instance()->GetOutputDirectory()=="") || (HeartConfig::Instance()->GetOutputFilenamePrefix()==""))
00224         {
00225             EXCEPTION("Either explicitly specify not to print output (call PrintOutput(false)) or specify the output directory and filename prefix");
00226         }
00227     }
00228 
00229     double end_time = HeartConfig::Instance()->GetSimulationDuration();
00230     double pde_time = HeartConfig::Instance()->GetPdeTimeStep();
00231 
00232     // MatrixIsConstant stuff requires CONSTANT dt - do some checks to make sure the TimeStepper won't find
00233     // non-constant dt.
00234     // Note: printing_time does not have to divide end_time, but dt must divide printing_time and end_time.
00235     // HeartConfig checks pde_dt divides printing dt
00236     if( fabs( end_time - pde_time*round(end_time/pde_time)) > 1e-10 )
00237     {
00238         EXCEPTION("PDE timestep does not seem to divide end time - check parameters");
00239     }
00240 }
00241 
00242 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00243 Vec AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::CreateInitialCondition()
00244 {
00245     DistributedVectorFactory* p_factory = mpMesh->GetDistributedVectorFactory();
00246     Vec initial_condition = p_factory->CreateVec(PROBLEM_DIM);
00247     DistributedVector ic = p_factory->CreateDistributedVector(initial_condition);
00248     std::vector<DistributedVector::Stripe> stripe;
00249     stripe.reserve(PROBLEM_DIM);
00250 
00251     for (unsigned i=0; i<PROBLEM_DIM; i++)
00252     {
00253         stripe.push_back(DistributedVector::Stripe(ic, i));
00254     }
00255 
00256     for (DistributedVector::Iterator index = ic.Begin();
00257          index != ic.End();
00258          ++index)
00259     {
00260         stripe[0][index] = mpCardiacTissue->GetCardiacCell(index.Global)->GetVoltage();
00261         if (PROBLEM_DIM==2)
00262         {
00263             stripe[1][index] = 0;
00264         }
00265     }
00266 
00267     ic.Restore();
00268 
00269     return initial_condition;
00270 }
00271 
00272 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00273 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetMesh(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00274 {
00275     // If this fails the mesh has already been set. We assert rather throw an exception
00276     // to avoid a memory leak when checking it throws correctly
00277     assert(mpMesh==NULL);
00278     assert(pMesh!=NULL);
00279     mAllocatedMemoryForMesh = false;
00280     mpMesh = pMesh;
00281 }
00282 
00283 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00284 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::PrintOutput(bool printOutput)
00285 {
00286     mPrintOutput = printOutput;
00287 }
00288 
00289 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00290 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetWriteInfo(bool writeInfo)
00291 {
00292     mWriteInfo = writeInfo;
00293 }
00294 
00295 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00296 Vec AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetSolution()
00297 {
00298     return mSolution;
00299 }
00300 
00301 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00302 DistributedVector AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetSolutionDistributedVector()
00303 {
00304     return mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(mSolution);
00305 }
00306 
00307 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00308 double AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetCurrentTime()
00309 {
00310     return mCurrentTime;
00311 }
00312 
00313 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00314 AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM> & AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::rGetMesh()
00315 {
00316     assert (mpMesh);
00317     return *mpMesh;
00318 }
00319 
00320 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00321 AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>* AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetTissue()
00322 {
00323     return mpCardiacTissue;
00324 }
00325 
00326 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00327 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetUseTimeAdaptivityController(
00328         bool useAdaptivity, 
00329         AbstractTimeAdaptivityController* pController)
00330 {
00331     if(useAdaptivity)
00332     {
00333         assert(pController);
00334         mpTimeAdaptivityController = pController;
00335     }
00336     else
00337     {
00338         mpTimeAdaptivityController = NULL;
00339     } 
00340 }
00341 
00342 
00343 
00344 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00345 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::Solve()
00346 {
00347     PreSolveChecks();
00348 
00349     if (!mpBoundaryConditionsContainer) // the user didn't supply a bcc
00350     {
00351         // Set up the default bcc
00352         mpDefaultBoundaryConditionsContainer.reset(new BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>);
00353         for (unsigned problem_index=0; problem_index<PROBLEM_DIM; problem_index++)
00354         {
00355             mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(mpMesh, problem_index);
00356         }
00357         mpBoundaryConditionsContainer = mpDefaultBoundaryConditionsContainer;
00358     }
00359 
00360     assert(mpSolver==NULL);
00361     mpSolver = CreateSolver(); // passes mpBoundaryConditionsContainer to solver
00362 
00363     // If we have already run a simulation, use the old solution as initial condition
00364     Vec initial_condition;
00365     if (mSolution)
00366     {
00367         initial_condition = mSolution;
00368     }
00369     else
00370     {
00371         initial_condition = CreateInitialCondition();
00372     }
00373 
00374     std::vector<double> additional_stopping_times;
00375     SetUpAdditionalStoppingTimes(additional_stopping_times);
00376 
00377     TimeStepper stepper(mCurrentTime,
00378                         HeartConfig::Instance()->GetSimulationDuration(),
00379                         HeartConfig::Instance()->GetPrintingTimeStep(),
00380                         false,
00381                         additional_stopping_times);
00382 
00383     std::string progress_reporter_dir;
00384 
00385     if (mPrintOutput)
00386     {
00387         HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
00388         try
00389         {
00390             InitialiseWriter();
00391         }
00392         catch (Exception& e)
00393         {
00394             delete mpWriter;
00395             mpWriter = NULL;
00396             delete mpSolver;
00397             if (mSolution != initial_condition)
00398             {
00399                 // A PETSc Vec is a pointer, so we *don't* need to free the memory if it is free'd somewhere else (e.g.
00400                 // in the destructor). If this is a resumed solution we set initial_condition = mSolution earlier. 
00401                 // mSolution is going to be cleaned up in the constructor. So, only VecDestroy initial_condition when it 
00402                 // is not equal to mSolution (see #1695). 
00403                 VecDestroy(initial_condition);
00404             }
00405             throw e;
00406         }
00407 
00408         // If we are resuming a simulation (i.e. mSolution already exists) there's no need
00409         // of writing the initial timestep,
00410         // since it was already written as the last timestep of the previous run
00411         if (!mSolution)
00412         {
00413             WriteOneStep(stepper.GetTime(), initial_condition);
00414             mpWriter->AdvanceAlongUnlimitedDimension();
00415         }
00416         HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
00417 
00418         progress_reporter_dir = HeartConfig::Instance()->GetOutputDirectory();
00419     }
00420     else
00421     {
00422         progress_reporter_dir = ""; // progress printed to CHASTE_TEST_OUTPUT
00423     }
00424 
00425     // create a progress reporter so users can track how much has gone and
00426     // estimate how much time is left. (Note this has to be done after the
00427     // InitialiseWriter above (if mPrintOutput==true)
00428     ProgressReporter progress_reporter(progress_reporter_dir, mCurrentTime,
00429                                        HeartConfig::Instance()->GetSimulationDuration());
00430     progress_reporter.Update(mCurrentTime);
00431 
00432     mpSolver->SetTimeStep(HeartConfig::Instance()->GetPdeTimeStep());
00433     if(mpTimeAdaptivityController)
00434     {
00435         mpSolver->SetTimeAdaptivityController(mpTimeAdaptivityController);
00436     }
00437 
00438     while ( !stepper.IsTimeAtEnd() )
00439     {
00440         // solve from now up to the next printing time
00441         mpSolver->SetTimes(stepper.GetTime(), stepper.GetNextTime());
00442         mpSolver->SetInitialCondition( initial_condition );
00443 
00444         AtBeginningOfTimestep(stepper.GetTime());
00445 
00446         try
00447         {
00448             mSolution = mpSolver->Solve();
00449         }
00450         catch (Exception &e)
00451         {
00452             // Free memory.
00453             delete mpSolver;
00454             mpSolver=NULL;
00455             if (initial_condition != mSolution) 
00456             {
00457                 // A PETSc Vec is a pointer, so we *don't* need to free the memory if it is free'd somewhere else (e.g.
00458                 // in the destructor). Later, in this while loop we will set initial_condition = mSolution (or, if this is
00459                 // a resumed solution it may also have been done when initial_condition was created). mSolution is going 
00460                 // to be cleaned up in the constructor. So, only VecDestroy initial_condition when it is not equal to 
00461                 // mSolution (see #1695). 
00462                 VecDestroy(initial_condition);
00463             }
00464             
00465 #ifndef NDEBUG
00466             PetscTools::ReplicateException(true);
00467 #endif
00468             // Re-throw
00469             HeartEventHandler::Reset();//EndEvent(HeartEventHandler::EVERYTHING);
00471             CloseFilesAndPostProcess();
00472             throw e;
00473         }
00474 #ifndef NDEBUG
00475         PetscTools::ReplicateException(false);
00476 #endif
00477 
00478         // Free old initial condition
00479         HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00480         VecDestroy(initial_condition);
00481         HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00482 
00483         // Initial condition for next loop is current solution
00484         initial_condition = mSolution;
00485 
00486         // update the current time
00487         stepper.AdvanceOneTimeStep();
00488         mCurrentTime = stepper.GetTime();
00489 
00490         // print out details at current time if asked for
00491         if (mWriteInfo)
00492         {
00493             HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
00494             WriteInfo(stepper.GetTime());
00495             HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
00496         }
00497 
00498         if (mPrintOutput)
00499         {
00500             // Writing data out to the file <FilenamePrefix>.dat
00501             HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
00502             WriteOneStep(stepper.GetTime(), mSolution);
00503             // Just flags that we've finished a time-step; won't actually 'extend' unless new data is written.
00504             mpWriter->AdvanceAlongUnlimitedDimension();
00505 
00506             HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
00507         }
00508 
00509         progress_reporter.Update(stepper.GetTime());
00510 
00511         OnEndOfTimestep(stepper.GetTime());
00512     }
00513 
00514     // Free solver
00515     delete mpSolver;
00516     mpSolver = NULL;
00517 
00518     // close the file that stores voltage values
00519     progress_reporter.PrintFinalising();
00520     CloseFilesAndPostProcess();
00521     HeartEventHandler::EndEvent(HeartEventHandler::EVERYTHING);
00522 }
00523 
00524 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00525 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::CloseFilesAndPostProcess()
00526 {
00527     // Close files
00528     if (!mPrintOutput)
00529     {
00530         //Nothing to do
00531         return;
00532     }
00533     delete mpWriter;
00534     mpWriter = NULL;
00535 
00536     HeartEventHandler::BeginEvent(HeartEventHandler::POST_PROC);
00537     // Only if results files were written and we are outputting all nodes
00538     if (mNodesToOutput.empty())
00539     {
00540         if (HeartConfig::Instance()->GetVisualizeWithMeshalyzer())
00541         {
00542             //Convert simulation data to Meshalyzer format
00543             Hdf5ToMeshalyzerConverter<ELEMENT_DIM,SPACE_DIM> converter(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix(), mpMesh);
00544         }
00545 
00546         if (HeartConfig::Instance()->GetVisualizeWithCmgui())
00547         {
00548             //Convert simulation data to Cmgui format
00549             Hdf5ToCmguiConverter<ELEMENT_DIM,SPACE_DIM> converter(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix(), mpMesh, GetHasBath());
00550         }
00551 
00552         if (HeartConfig::Instance()->GetVisualizeWithVtk())
00553         {
00554             //Convert simulation data to VTK format
00555             Hdf5ToVtkConverter<ELEMENT_DIM,SPACE_DIM> converter(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix(), mpMesh, false);
00556         }
00557 
00558         if (HeartConfig::Instance()->GetVisualizeWithParallelVtk())
00559         {
00560             //Convert simulation data to parallel VTK (pvtu) format
00561             Hdf5ToVtkConverter<ELEMENT_DIM,SPACE_DIM> converter(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix(), mpMesh, true);
00562         }
00563     }
00564 
00565     if(HeartConfig::Instance()->IsPostProcessingRequested())
00566     {
00567         PostProcessingWriter<ELEMENT_DIM, SPACE_DIM> post_writer(*mpMesh, HeartConfig::Instance()->GetOutputDirectory(),
00568                         HeartConfig::Instance()->GetOutputFilenamePrefix(), true);
00569         post_writer.WritePostProcessingFiles();
00570     }
00571 
00572     HeartEventHandler::EndEvent(HeartEventHandler::POST_PROC);
00573 }
00574 
00575 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00576 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::DefineWriterColumns(bool extending)
00577 {
00578     if (!extending)
00579     {
00580         if ( mNodesToOutput.empty() )
00581         {
00582             //Set writer to output all nodes
00583             mpWriter->DefineFixedDimension( mpMesh->GetNumNodes() );
00584         }
00585         else
00586         {
00587             //Output only the nodes indicted
00588             mpWriter->DefineFixedDimension( mNodesToOutput, mpMesh->GetNumNodes() );
00589         }
00590         //mNodeColumnId = mpWriter->DefineVariable("Node", "dimensionless");
00591         mVoltageColumnId = mpWriter->DefineVariable("V","mV");
00592 
00593         // Only used to get an estimate of the # of timesteps below
00594         TimeStepper stepper(mCurrentTime,
00595                             HeartConfig::Instance()->GetSimulationDuration(),
00596                             HeartConfig::Instance()->GetPrintingTimeStep());
00597 
00598         mpWriter->DefineUnlimitedDimension("Time","msecs", stepper.EstimateTimeSteps()+1); // plus one for start and end points
00599     }
00600     else
00601     {
00602         mVoltageColumnId = mpWriter->GetVariableByName("V");
00603     }
00604 }
00605 
00606 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00607 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::DefineExtraVariablesWriterColumns(bool extending)
00608 {
00609     mExtraVariablesId.clear();
00610     // Check if any extra output variables have been requested
00611     if (HeartConfig::Instance()->GetOutputVariablesProvided())
00612     {
00613         // Get their names in a vector
00614         std::vector<std::string> output_variables;
00615         HeartConfig::Instance()->GetOutputVariables(output_variables);
00616         const unsigned num_vars = output_variables.size();
00617         mExtraVariablesId.reserve(num_vars);
00618 
00619         // Loop over them
00620         for (unsigned var_index=0; var_index<num_vars; var_index++)
00621         {
00622             // Get variable name
00623             std::string var_name = output_variables[var_index];
00624 
00625             // Register it (or look it up) in the data writer
00626             unsigned column_id;
00627             if (extending)
00628             {
00629                 column_id = this->mpWriter->GetVariableByName(var_name);
00630             }
00631             else
00632             {
00633                 column_id = this->mpWriter->DefineVariable(var_name, "");
00634             }
00635 
00636             // Store column id
00637             mExtraVariablesId.push_back(column_id);
00638         }
00639     }
00640 }
00641 
00642 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00643 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::WriteExtraVariablesOneStep()
00644 {
00645     // Get the variable names in a vector
00646     std::vector<std::string> output_variables;
00647     unsigned num_vars = mExtraVariablesId.size();
00648     if (num_vars > 0)
00649     {
00650         HeartConfig::Instance()->GetOutputVariables(output_variables);
00651     }
00652     assert(output_variables.size() == num_vars);
00653 
00654     // Loop over the requested variables
00655     for (unsigned var_index=0; var_index<num_vars; var_index++)
00656     {
00657         // Create vector for storing values over the local nodes
00658         Vec variable_data =  this->mpMesh->GetDistributedVectorFactory()->CreateVec();
00659         DistributedVector distributed_var_data = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(variable_data);
00660 
00661         // Loop over the local nodes and gather the data
00662         for (DistributedVector::Iterator index = distributed_var_data.Begin();
00663              index!= distributed_var_data.End();
00664              ++index)
00665         {
00666             // Find the variable in the cell model
00667             AbstractCardiacCell* p_cell = this->mpCardiacTissue->GetCardiacCell(index.Global);
00668             unsigned cell_id = p_cell->GetAnyVariableIndex(output_variables[var_index]);
00669             // Store its value
00670             distributed_var_data[index] = p_cell->GetAnyVariable(cell_id, mCurrentTime);
00671         }
00672         distributed_var_data.Restore();
00673 
00674         // Write it to disc
00675         this->mpWriter->PutVector(mExtraVariablesId[var_index], variable_data);
00676 
00677         VecDestroy(variable_data);
00678     }
00679 }
00680 
00681 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00682 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::InitialiseWriter()
00683 {
00684     bool extend_file = (mSolution != NULL);
00685 
00686     // I think this is impossible to trip; certainly it's very difficult!
00687     assert(!mpWriter);
00688 
00689     try
00690     {
00691         mpWriter = new Hdf5DataWriter(*mpMesh->GetDistributedVectorFactory(),
00692                                       HeartConfig::Instance()->GetOutputDirectory(),
00693                                       HeartConfig::Instance()->GetOutputFilenamePrefix(),
00694                                       !extend_file, // don't clear directory if extension requested
00695                                       extend_file);
00696     }
00697     catch (Exception& e)
00698     {
00699         // The constructor only throws an Exception if we're extending
00700         assert(extend_file);
00701         // Tried to extend and failed, so just create from scratch
00702         extend_file = false;
00703         mpWriter = new Hdf5DataWriter(*mpMesh->GetDistributedVectorFactory(),
00704                                       HeartConfig::Instance()->GetOutputDirectory(),
00705                                       HeartConfig::Instance()->GetOutputFilenamePrefix(),
00706                                       !extend_file,
00707                                       extend_file);
00708     }
00709 
00710     // Define columns, or get the variable IDs from the writer
00711     DefineWriterColumns(extend_file);
00712 
00713     //Possibility of applying a permutation
00714     if (HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering())
00715     {
00716         bool success = mpWriter->ApplyPermutation(mpMesh->rGetNodePermutation());
00717         if (success == false)
00718         {
00719             //It's not really a permutation, so reset
00720             HeartConfig::Instance()->SetOutputUsingOriginalNodeOrdering(false);
00721         }
00722     }
00723 
00724 
00725     if (!extend_file)
00726     {
00727         mpWriter->EndDefineMode();
00728     }
00729 }
00730 
00731 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00732 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetOutputNodes(std::vector<unsigned> &nodesToOutput)
00733 {
00734     mNodesToOutput = nodesToOutput;
00735 }
00736 
00737 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00738 Hdf5DataReader AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetDataReader()
00739 {
00740     if( (HeartConfig::Instance()->GetOutputDirectory()=="") || (HeartConfig::Instance()->GetOutputFilenamePrefix()==""))
00741     {
00742         EXCEPTION("Data reader invalid as data writer cannot be initialised");
00743     }
00744     return Hdf5DataReader(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix());
00745 }
00746 
00747 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00748 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::UseMatrixBasedRhsAssembly(bool usematrix)
00749 {
00750     mUseMatrixBasedRhsAssembly = usematrix;
00751 }
00752 
00753 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00754 bool AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetHasBath()
00755 {
00756     return false;
00757 }
00758 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00759 void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::SetElectrodes()
00760 {
00761 }
00762 
00763 
00765 // Explicit instantiation
00767 
00768 // Monodomain
00769 template class AbstractCardiacProblem<1,1,1>;
00770 template class AbstractCardiacProblem<1,2,1>;
00771 template class AbstractCardiacProblem<1,3,1>;
00772 template class AbstractCardiacProblem<2,2,1>;
00773 template class AbstractCardiacProblem<3,3,1>;
00774 
00775 // Bidomain
00776 template class AbstractCardiacProblem<1,1,2>;
00777 template class AbstractCardiacProblem<2,2,2>;
00778 template class AbstractCardiacProblem<3,3,2>;

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