HeartConfig.cpp

00001 /*
00002 
00003 Copyright (c) 2005-2015, University of Oxford.
00004 All rights reserved.
00005 
00006 University of Oxford means the Chancellor, Masters and Scholars of the
00007 University of Oxford, having an administrative office at Wellington
00008 Square, Oxford OX1 2JD, UK.
00009 
00010 This file is part of Chaste.
00011 
00012 Redistribution and use in source and binary forms, with or without
00013 modification, are permitted provided that the following conditions are met:
00014  * Redistributions of source code must retain the above copyright notice,
00015    this list of conditions and the following disclaimer.
00016  * Redistributions in binary form must reproduce the above copyright notice,
00017    this list of conditions and the following disclaimer in the documentation
00018    and/or other materials provided with the distribution.
00019  * Neither the name of the University of Oxford nor the names of its
00020    contributors may be used to endorse or promote products derived from this
00021    software without specific prior written permission.
00022 
00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00033 
00034 */
00035 
00036 // Work-around for newer Boost versions
00037 #include "CheckpointArchiveTypesIfNeeded.hpp"
00038 
00039 #include "UblasCustomFunctions.hpp"
00040 
00041 #include "HeartConfig.hpp"
00042 #include "OutputFileHandler.hpp"
00043 #include "Exception.hpp"
00044 #include "ChastePoint.hpp"
00045 #include "Version.hpp"
00046 #include "AbstractChasteRegion.hpp"
00047 #include "HeartFileFinder.hpp"
00048 #include "Warnings.hpp"
00049 
00050 #include "HeartRegionCodes.hpp"
00051 
00052 #include "SimpleStimulus.hpp"
00053 #include "RegularStimulus.hpp"
00054 
00055 #include <string>
00056 #include <istream>
00057 #include <fstream>
00058 #include <cassert>
00059 #include <map>
00060 
00061 #include "XmlTools.hpp"
00062 #include <xsd/cxx/tree/exceptions.hxx>
00063 using namespace xsd::cxx::tree;
00064 
00065 // Coping with changes to XSD interface
00066 #if (XSD_INT_VERSION >= 3000000L)
00067 #define XSD_SEQUENCE_TYPE(base) base##_sequence
00068 #define XSD_ITERATOR_TYPE(base) base##_iterator
00069 #define XSD_NESTED_TYPE(t) t##_type
00070 #define XSD_ANON_TYPE(t1, t2) \
00071     t1::t2##_type
00072 #else
00073 #define XSD_SEQUENCE_TYPE(base) base::container
00074 #define XSD_ITERATOR_TYPE(base) base::iterator
00075 #define XSD_NESTED_TYPE(t) t::type
00076 #define XSD_ANON_TYPE(t1, t2) \
00077     t1::t2::_xsd_##t2##_::t2
00078 #endif
00079 
00080 // These are for convenience
00081 #define XSD_ANON_SEQUENCE_TYPE(t1, t2, t3) \
00082     XSD_SEQUENCE_TYPE(XSD_ANON_TYPE(t1, t2)::t3)
00083 #define XSD_ANON_ITERATOR_TYPE(t1, t2, t3) \
00084     XSD_ITERATOR_TYPE(XSD_ANON_TYPE(t1, t2)::t3)
00085 
00086 // Newer versions don't allow you to set fixed attributes
00087 #if (XSD_INT_VERSION >= 3020000L)
00088 #define XSD_CREATE_WITH_FIXED_ATTR(type, name, attr) \
00089     type name
00090 #define XSD_CREATE_WITH_FIXED_ATTR1(type, name, arg1, attr) \
00091     type name(arg1)
00092 #define XSD_CREATE_WITH_FIXED_ATTR2(type, name, arg1, arg2, attr) \
00093     type name(arg1, arg2)
00094 #define XSD_CREATE_WITH_FIXED_ATTR3(type, name, arg1, arg2, arg3, attr) \
00095     type name(arg1, arg2, arg3)
00096 #else
00097 #define XSD_CREATE_WITH_FIXED_ATTR(type, name, attr) \
00098     type name(attr)
00099 #define XSD_CREATE_WITH_FIXED_ATTR1(type, name, arg1, attr) \
00100     type name(arg1, attr)
00101 #define XSD_CREATE_WITH_FIXED_ATTR2(type, name, arg1, arg2, attr) \
00102     type name(arg1, arg2, attr)
00103 #define XSD_CREATE_WITH_FIXED_ATTR3(type, name, arg1, arg2, arg3, attr) \
00104     type name(arg1, arg2, arg3, attr)
00105 #endif
00106 
00110 #define ENSURE_SECTION_PRESENT(location, type) \
00111     if (!location.present())                   \
00112     {                                          \
00113         type empty_item;                       \
00114         location.set(empty_item);              \
00115     }
00116 
00117 #include <boost/current_function.hpp>
00125 #define CHECK_EXISTS(test, path)                                                         \
00126     do {                                                                                 \
00127         if (!test) {                                                                     \
00128             EXCEPTION("No XML element " << path << " found in parameters when calling '" \
00129                       << BOOST_CURRENT_FUNCTION << "'");                                 \
00130     }} while (false)
00131 
00132 
00137 class XmlTransforms
00138 {
00139 public:
00147     static void TransformIonicModelDefinitions(xercesc::DOMDocument* pDocument,
00148                                                xercesc::DOMElement* pRootElement);
00149 
00158     static void TransformArchiveDirectory(xercesc::DOMDocument* pDocument,
00159                                           xercesc::DOMElement* pRootElement);
00160 
00168     static void CheckForIluPreconditioner(xercesc::DOMDocument* pDocument,
00169                                           xercesc::DOMElement* pRootElement);
00170 
00178     static void MoveConductivityHeterogeneities(xercesc::DOMDocument* pDocument,
00179                                                 xercesc::DOMElement* pRootElement);
00180 
00188     static void SetDefaultVisualizer(xercesc::DOMDocument* pDocument,
00189                                      xercesc::DOMElement* pRootElement);
00190 };
00191 
00192 //
00193 // Default settings
00194 //
00195 #include "HeartConfigDefaults.hpp"
00196 
00197 //
00198 // Definition of static member variables
00199 //
00200 std::auto_ptr<HeartConfig> HeartConfig::mpInstance;
00201 
00202 //
00203 // Methods
00204 //
00205 
00206 HeartConfig* HeartConfig::Instance()
00207 {
00208     if (mpInstance.get() == NULL)
00209     {
00210         mpInstance.reset(new HeartConfig);
00211     }
00212     return mpInstance.get();
00213 }
00214 
00215 HeartConfig::HeartConfig()
00216     : mUseMassLumping(false),
00217       mUseMassLumpingForPrecond(false),
00218       mUseFixedNumberIterations(false),
00219       mEvaluateNumItsEveryNSolves(UINT_MAX)
00220 {
00221     assert(mpInstance.get() == NULL);
00222     mUseFixedSchemaLocation = true;
00223     SetDefaultSchemaLocations();
00224 
00225     mpParameters = CreateDefaultParameters();
00226     //CheckTimeSteps(); // necessity of this line of code is not tested -- remove with caution!
00227 
00228     //initialise the member variable of the layers
00229     mEpiFraction = -1.0;
00230     mEndoFraction =  -1.0;
00231     mMidFraction = -1.0;
00232     mUserAskedForCellularTransmuralHeterogeneities = false;
00233     // initialise to senseless values (these should be only 0, 1 and 2)
00234     // note: the 'minus 3' is for checking purposes as we need to add 0, 1 or 2 to this initial value
00235     // and UINT_MAX+1 seems to be 0
00236     mIndexMid = UINT_MAX-3u;
00237     mIndexEpi = UINT_MAX-3u;
00238     mIndexEndo = UINT_MAX-3u;
00239 
00240     mUseReactionDiffusionOperatorSplitting = false;
00241 
00243     mTissueIdentifiers.insert(0);
00244     mBathIdentifiers.insert(1);
00245 }
00246 
00247 HeartConfig::~HeartConfig()
00248 {
00249 }
00250 
00251 
00252 void HeartConfig::Write(bool useArchiveLocationInfo, std::string subfolderName)
00253 {
00254     //Output file
00255     std::string output_dirname;
00256     if (useArchiveLocationInfo)
00257     {
00258         output_dirname = ArchiveLocationInfo::GetArchiveDirectory();
00259     }
00260     else
00261     {
00262         OutputFileHandler handler(GetOutputDirectory() + "/" + subfolderName, false);
00263         output_dirname = handler.GetOutputDirectoryFullPath();
00264     }
00265 
00266     // Sometimes this method is called collectively, sometimes not,
00267     // in any case the below file operations only want to be performed by
00268     // the master - so exit. Caller takes responsibility for nice
00269     // exception handling.
00270     if (!PetscTools::AmMaster())
00271     {
00272         return;
00273     }
00274 
00275     out_stream p_parameters_file( new std::ofstream( (output_dirname+"ChasteParameters.xml").c_str() ) );
00276 
00277     if (!p_parameters_file->is_open())
00278     {
00279         EXCEPTION("Could not open XML file in HeartConfig");
00280     }
00281 
00282     //Schema map
00283     //Note - this location is relative to where we are storing the xml
00284     ::xml_schema::namespace_infomap map;
00285     // Release 1.1 (and earlier) didn't use a namespace
00286     map[""].schema = "ChasteParameters_1_1.xsd";
00287     // Later releases use namespaces of the form https://chaste.comlab.ox.ac.uk/nss/parameters/N_M
00288     map["cp20"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/2_0";
00289     map["cp20"].schema = "ChasteParameters_2_0.xsd";
00290     map["cp21"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/2_1";
00291     map["cp21"].schema = "ChasteParameters_2_1.xsd";
00292     map["cp22"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/2_2";
00293     map["cp22"].schema = "ChasteParameters_2_2.xsd";
00294     map["cp23"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/2_3";
00295     map["cp23"].schema = "ChasteParameters_2_3.xsd";
00296     map["cp30"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/3_0";
00297     map["cp30"].schema = "ChasteParameters_3_0.xsd";
00298     map["cp31"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/3_1";
00299     map["cp31"].schema = "ChasteParameters_3_1.xsd";
00300     // We use 'cp' as prefix for the latest version to avoid having to change saved
00301     // versions for comparison at every release.
00302     map["cp"].name = "https://chaste.comlab.ox.ac.uk/nss/parameters/3_3";
00303     map["cp"].schema = "ChasteParameters_3_3.xsd";
00304 
00305     cp::ChasteParameters(*p_parameters_file, *mpParameters, map);
00306 
00307     // If we're archiving, try to save a copy of the latest schema too
00308     if (useArchiveLocationInfo)
00309     {
00310         CopySchema(output_dirname);
00311     }
00312 
00313 }
00314 
00315 void HeartConfig::LoadFromCheckpoint()
00316 {
00317     /*
00318      *  This method implements the logic required by HeartConfig to be able to handle resuming a simulation via the executable.
00319      *
00320      *  When the control reaches the method mpParameters points to the file specified as resuming parameters.
00321      *  However SetParametersFile() will set this variable to point to the archived parameters.
00322      *
00323      *  We make a temporary copy of mpParameters so we don't lose its content.
00324      *  At the end of the method we update the new mpParameters with the resuming parameters.
00325      */
00326     assert(mpParameters.use_count() > 0);
00327     boost::shared_ptr<cp::chaste_parameters_type> p_new_parameters = mpParameters;
00328 
00329     /*
00330      *  When we unarchive a simulation, we load the old parameters file in order to inherit things such
00331      *  as default cell model, stimuli, heterogeneities, ... This has the side effect of inheriting the
00332      *  <CheckpointSimulation> element (if defined).
00333      *
00334      *  We disable checkpointing definition coming from the unarchived config file. We will enable it again
00335      *  if defined in the resume config file.
00336      */
00337     std::string parameters_filename_xml = ArchiveLocationInfo::GetArchiveDirectory() + "ChasteParameters.xml";
00338     mpParameters = ReadFile(parameters_filename_xml);
00339     mParametersFilePath.SetPath(parameters_filename_xml, RelativeTo::AbsoluteOrCwd);
00340 
00341     // Release 3.0 and earlier wrote a separate defaults file in the checkpoint
00342     std::string defaults_filename_xml = ArchiveLocationInfo::GetArchiveDirectory() + "ChasteDefaults.xml";
00343     if (FileFinder(defaults_filename_xml).Exists())
00344     {
00345         boost::shared_ptr<cp::chaste_parameters_type> p_defaults = ReadFile(defaults_filename_xml);
00346         MergeDefaults(mpParameters, p_defaults);
00347     }
00348 
00349     HeartConfig::Instance()->SetCheckpointSimulation(false);
00350 
00351     // If we are resuming a simulation, some parameters can be altered at this point.
00352     if (p_new_parameters->ResumeSimulation().present())
00353     {
00354         UpdateParametersFromResumeSimulation(p_new_parameters);
00355     }
00356 
00357     CheckTimeSteps(); // For consistency with SetParametersFile
00358 }
00359 
00360 void HeartConfig::CopySchema(const std::string& rToDirectory)
00361 {
00362     // N.B. This method should only be called by the master process,
00363     // in a situation where it can handle EXCEPTION()s nicely, e.g.
00364     // TRY_IF_MASTER(CopySchema(...));
00365 
00366     std::string schema_name("ChasteParameters_3_3.xsd");
00367     FileFinder schema_location("heart/src/io/" + schema_name, RelativeTo::ChasteSourceRoot);
00368     if (!schema_location.Exists())
00369     {
00370         // Try a relative path instead
00371         schema_location.SetPath(schema_name, RelativeTo::CWD);
00372         if (!schema_location.Exists())
00373         {
00374             // Warn the user
00375             std::string message("Unable to locate schema file " + schema_name +
00376                                 ". You will need to ensure it is available when resuming from the checkpoint.");
00377             WARN_ONCE_ONLY(message);
00378         }
00379     }
00380     if (schema_location.Exists())
00381     {
00382         FileFinder output_directory(rToDirectory, RelativeTo::Absolute);
00383         schema_location.CopyTo(output_directory);
00384     }
00385 }
00386 
00387 void HeartConfig::SetDefaultSchemaLocations()
00388 {
00389     mSchemaLocations.clear();
00390     // Location of schemas in the source tree
00391     std::string root_dir = std::string(ChasteBuildInfo::GetRootDir()) + "/heart/src/io/";
00392     // Release 1.1 (and earlier) didn't use a namespace
00393     mSchemaLocations[""] = root_dir + "ChasteParameters_1_1.xsd";
00394     // Later releases use namespaces of the form https://chaste.comlab.ox.ac.uk/nss/parameters/N_M
00395     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/2_0"] = root_dir + "ChasteParameters_2_0.xsd";
00396     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/2_1"] = root_dir + "ChasteParameters_2_1.xsd";
00397     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/2_2"] = root_dir + "ChasteParameters_2_2.xsd";
00398     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/2_3"] = root_dir + "ChasteParameters_2_3.xsd";
00399     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/3_0"] = root_dir + "ChasteParameters_3_0.xsd";
00400     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/3_1"] = root_dir + "ChasteParameters_3_1.xsd";
00401     mSchemaLocations["https://chaste.comlab.ox.ac.uk/nss/parameters/3_3"] = root_dir + "ChasteParameters_3_3.xsd";
00402 }
00403 
00404 unsigned HeartConfig::GetVersionFromNamespace(const std::string& rNamespaceUri)
00405 {
00406     unsigned version_major = 0;
00407     unsigned version_minor = 0;
00408     if (rNamespaceUri == "")
00409     {
00410         version_major = 1;
00411         version_minor = 1;
00412     }
00413     else
00414     {
00415         std::string uri_base("https://chaste.comlab.ox.ac.uk/nss/parameters/");
00416         if (rNamespaceUri.substr(0, uri_base.length()) == uri_base)
00417         {
00418             std::istringstream version_string(rNamespaceUri.substr(uri_base.length()));
00419             version_string >> version_major;
00420             version_string.ignore(1);
00421             version_string >> version_minor;
00422             if (version_string.fail())
00423             {
00424                 version_major = 0;
00425                 version_minor = 0;
00426             }
00427         }
00428     }
00429 
00430     unsigned version = version_major * 1000 + version_minor;
00431     if (version == 0)
00432     {
00433         EXCEPTION(rNamespaceUri + " is not a recognised Chaste parameters namespace.");
00434     }
00435     return version;
00436 }
00437 
00438 void HeartConfig::SetFixedSchemaLocations(const SchemaLocationsMap& rSchemaLocations)
00439 {
00440     mSchemaLocations = rSchemaLocations;
00441     SetUseFixedSchemaLocation(true);
00442 }
00443 
00444 void HeartConfig::SetUseFixedSchemaLocation(bool useFixedSchemaLocation)
00445 {
00446     mUseFixedSchemaLocation = useFixedSchemaLocation;
00447 }
00448 
00449 boost::shared_ptr<cp::chaste_parameters_type> HeartConfig::ReadFile(const std::string& rFileName)
00450 {
00451     // Determine whether to use the schema path given in the input XML, or our own schema
00452     ::xml_schema::properties props;
00453     if (mUseFixedSchemaLocation)
00454     {
00455         for (SchemaLocationsMap::iterator it = mSchemaLocations.begin();
00456              it != mSchemaLocations.end();
00457              ++it)
00458         {
00459             if (it->first == "")
00460             {
00461                 props.no_namespace_schema_location(XmlTools::EscapeSpaces(it->second));
00462             }
00463             else
00464             {
00465                 props.schema_location(it->first, XmlTools::EscapeSpaces(it->second));
00466             }
00467         }
00468     }
00469 
00470     // Get the parameters using the method 'ChasteParameters(rFileName)',
00471     // which returns a std::auto_ptr. We convert to a shared_ptr for easier semantics.
00472     try
00473     {
00474         // Make sure Xerces finalization happens
00475         XmlTools::Finalizer finalizer(false);
00476         // Parse XML to DOM
00477         xsd::cxx::xml::dom::auto_ptr<xercesc::DOMDocument> p_doc = XmlTools::ReadXmlFile(rFileName, props);
00478         // Test the namespace on the root element
00479         xercesc::DOMElement* p_root_elt = p_doc->getDocumentElement();
00480         std::string namespace_uri(X2C(p_root_elt->getNamespaceURI()));
00481         const unsigned version = GetVersionFromNamespace(namespace_uri);
00482         if (version < 2000) // Changes made in release 2.0
00483         {
00484             XmlTransforms::TransformIonicModelDefinitions(p_doc.get(), p_root_elt);
00485         }
00486         if (version < 2001) // Changes made in release 2.1
00487         {
00488             XmlTransforms::TransformArchiveDirectory(p_doc.get(), p_root_elt);
00489             XmlTransforms::CheckForIluPreconditioner(p_doc.get(), p_root_elt);
00490         }
00491         if (version < 3001) // Changes made in release 3.1
00492         {
00493             XmlTransforms::MoveConductivityHeterogeneities(p_doc.get(), p_root_elt);
00494         }
00495         if (version < 3003) // Not the latest release
00496         {
00497             XmlTransforms::SetDefaultVisualizer(p_doc.get(), p_root_elt);
00498             XmlTools::SetNamespace(p_doc.get(), p_root_elt, "https://chaste.comlab.ox.ac.uk/nss/parameters/3_3");
00499         }
00500         // Parse DOM to object model
00501         std::auto_ptr<cp::chaste_parameters_type> p_params(cp::ChasteParameters(*p_doc, ::xml_schema::flags::dont_initialize, props));
00502         // Get rid of the DOM stuff
00503         p_doc.reset();
00504 
00505         return boost::shared_ptr<cp::chaste_parameters_type>(p_params);
00506     }
00507     catch (const xml_schema::exception& e)
00508     {
00509         std::cerr << e << std::endl;
00510         // Make sure we don't store invalid parameters
00511         mpParameters.reset();
00512         EXCEPTION("XML parsing error in configuration file: " + rFileName);
00513     }
00514     catch (...)
00515     {
00516         // Make sure we don't store invalid parameters
00517         mpParameters.reset();
00518         throw;
00519     }
00520 }
00521 
00522 void HeartConfig::SetParametersFile(const std::string& rFileName)
00523 {
00524     mpParameters = ReadFile(rFileName);
00525     MergeDefaults(mpParameters, CreateDefaultParameters());
00526     mParametersFilePath.SetPath(rFileName, RelativeTo::AbsoluteOrCwd);
00527 
00528     if (IsSimulationDefined())
00529     {
00530         CheckTimeSteps(); // Resume files might not have time steps defined
00531     }
00532 }
00533 
00534 FileFinder HeartConfig::GetParametersFilePath()
00535 {
00536     return mParametersFilePath;
00537 }
00538 
00539 void HeartConfig::UpdateParametersFromResumeSimulation(boost::shared_ptr<cp::chaste_parameters_type> pResumeParameters)
00540 {
00541     // Check for user foolishness
00542     if ( (pResumeParameters->ResumeSimulation()->SpaceDimension() != HeartConfig::Instance()->GetSpaceDimension())
00543          ||(pResumeParameters->ResumeSimulation()->Domain() != HeartConfig::Instance()->GetDomain()))
00544     {
00545         EXCEPTION("Problem type and space dimension should match when restarting a simulation.");
00546     }
00547 
00548     // New simulation duration
00549     HeartConfig::Instance()->SetSimulationDuration(pResumeParameters->ResumeSimulation()->SimulationDuration());
00550 
00551     // Stimulus definition.  For these we always replace any previous definitions (at least for now...)
00552     if (pResumeParameters->ResumeSimulation()->Stimuli().present())
00553     {
00554         mpParameters->Simulation()->Stimuli().set(pResumeParameters->ResumeSimulation()->Stimuli().get());
00555     }
00556 
00557     // Cell heterogeneities.  Note that while we copy the elements here, other code in CardiacSimulation actually updates
00558     // the loaded simulation to take account of the new settings.
00559     if (pResumeParameters->ResumeSimulation()->CellHeterogeneities().present())
00560     {
00561         if (!mpParameters->Simulation()->CellHeterogeneities().present())
00562         {
00563             // Original parameters had no heterogeneities, so just copy the whole element
00564             mpParameters->Simulation()->CellHeterogeneities().set(pResumeParameters->ResumeSimulation()->CellHeterogeneities().get());
00565         }
00566         else
00567         {
00568             // Need to append the new heterogeneity defitions to the original sequence
00569             XSD_SEQUENCE_TYPE(cp::cell_heterogeneities_type::CellHeterogeneity)&
00570                 new_seq = pResumeParameters->ResumeSimulation()->CellHeterogeneities()->CellHeterogeneity();
00571             XSD_SEQUENCE_TYPE(cp::cell_heterogeneities_type::CellHeterogeneity)&
00572                 orig_seq = mpParameters->Simulation()->CellHeterogeneities()->CellHeterogeneity();
00573             for (XSD_ITERATOR_TYPE(cp::cell_heterogeneities_type::CellHeterogeneity) i = new_seq.begin();
00574                  i != new_seq.end();
00575                  ++i)
00576             {
00577                 orig_seq.push_back(*i);
00578             }
00579         }
00580     }
00581 
00582     // Whether to checkpoint the resumed simulation
00583     if (pResumeParameters->ResumeSimulation()->CheckpointSimulation().present())
00584     {
00585         HeartConfig::Instance()->SetCheckpointSimulation(true,
00586                                                          pResumeParameters->ResumeSimulation()->CheckpointSimulation()->timestep(),
00587                                                          pResumeParameters->ResumeSimulation()->CheckpointSimulation()->max_checkpoints_on_disk());
00588     }
00589 
00590     //Visualization parameters are no longer compulsory
00591     if (pResumeParameters->ResumeSimulation()->OutputVisualizer().present())
00592     {
00593         HeartConfig::Instance()->SetVisualizeWithParallelVtk(pResumeParameters->ResumeSimulation()->OutputVisualizer()->parallel_vtk() == cp::yesno_type::yes);
00594         HeartConfig::Instance()->SetVisualizeWithVtk(pResumeParameters->ResumeSimulation()->OutputVisualizer()->vtk() == cp::yesno_type::yes);
00595         HeartConfig::Instance()->SetVisualizeWithCmgui(pResumeParameters->ResumeSimulation()->OutputVisualizer()->cmgui() == cp::yesno_type::yes);
00596         HeartConfig::Instance()->SetVisualizeWithMeshalyzer(pResumeParameters->ResumeSimulation()->OutputVisualizer()->meshalyzer() == cp::yesno_type::yes);
00597     }
00598 
00599     // Numerical parameters may be overridden
00600     {
00601         cp::numerical_type& r_resume = pResumeParameters->Numerical();
00602         cp::numerical_type& r_user = mpParameters->Numerical();
00603         if (r_resume.TimeSteps().present())
00604         {
00605             r_user.TimeSteps().set(r_resume.TimeSteps().get());
00606         }
00607         if (r_resume.KSPTolerances().present())
00608         {
00609             r_user.KSPTolerances().set(r_resume.KSPTolerances().get());
00610         }
00611         if (r_resume.KSPSolver().present())
00612         {
00613             r_user.KSPSolver().set(r_resume.KSPSolver().get());
00614         }
00615         if (r_resume.KSPPreconditioner().present())
00616         {
00617             r_user.KSPPreconditioner().set(r_resume.KSPPreconditioner().get());
00618         }
00619         if (r_resume.AdaptivityParameters().present())
00620         {
00621             r_user.AdaptivityParameters().set(r_resume.AdaptivityParameters().get());
00622         }
00623     }
00624 
00625     // Post-processing parameters may be overridden
00626     if (pResumeParameters->PostProcessing().present())
00627     {
00628         EnsurePostProcessingSectionPresent();
00629         cp::postprocessing_type& r_resume = pResumeParameters->PostProcessing().get();
00630         cp::postprocessing_type& r_user = mpParameters->PostProcessing().get();
00631         if (!r_resume.ActionPotentialDurationMap().empty())
00632         {
00633             r_user.ActionPotentialDurationMap() = r_resume.ActionPotentialDurationMap();
00634         }
00635         if (!r_resume.UpstrokeTimeMap().empty())
00636         {
00637             r_user.UpstrokeTimeMap() = r_resume.UpstrokeTimeMap();
00638         }
00639         if (!r_resume.MaxUpstrokeVelocityMap().empty())
00640         {
00641             r_user.MaxUpstrokeVelocityMap() = r_resume.MaxUpstrokeVelocityMap();
00642         }
00643         if (!r_resume.ConductionVelocityMap().empty())
00644         {
00645             r_user.ConductionVelocityMap() = r_resume.ConductionVelocityMap();
00646         }
00647     }
00648 }
00649 
00650 
00651 void HeartConfig::Reset()
00652 {
00653     // Throw it away first, so that mpInstance is NULL when we...
00654     mpInstance.reset();
00655     // ...make a new one
00656     mpInstance.reset(new HeartConfig);
00657 }
00658 
00659 bool HeartConfig::IsSimulationDefined() const
00660 {
00661     return mpParameters->Simulation().present();
00662 }
00663 
00664 bool HeartConfig::IsSimulationResumed() const
00665 {
00666     return mpParameters->ResumeSimulation().present();
00667 }
00668 
00669 
00670 void HeartConfig::CheckSimulationIsDefined(std::string callingMethod) const
00671 {
00672     if (IsSimulationResumed())
00673     {
00674         EXCEPTION(callingMethod + " information is not available in a resumed simulation.");
00675     }
00676 }
00677 
00678 void HeartConfig::CheckResumeSimulationIsDefined(std::string callingMethod) const
00679 {
00680     if (IsSimulationDefined())
00681     {
00682         EXCEPTION(callingMethod + " information is not available in a standard (non-resumed) simulation.");
00683     }
00684 }
00685 
00686 unsigned HeartConfig::GetSpaceDimension() const
00687 {
00688     if (IsSimulationDefined())
00689     {
00690         CHECK_EXISTS(mpParameters->Simulation()->SpaceDimension().present(), "Simulation/SpaceDimension");
00691         return mpParameters->Simulation()->SpaceDimension().get();
00692     }
00693     else
00694     {
00695         return mpParameters->ResumeSimulation()->SpaceDimension();
00696     }
00697 }
00698 
00699 double HeartConfig::GetSimulationDuration() const
00700 {
00701     if (IsSimulationDefined())
00702     {
00703         CHECK_EXISTS(mpParameters->Simulation()->SimulationDuration().present(), "Simulation/SimulationDuration");
00704         return mpParameters->Simulation()->SimulationDuration().get();
00705     }
00706     else // IsSimulationResumed
00707     {
00708         return mpParameters->ResumeSimulation()->SimulationDuration();
00709     }
00710 }
00711 
00712 cp::domain_type HeartConfig::GetDomain() const
00713 {
00714     if (IsSimulationDefined())
00715     {
00716         CHECK_EXISTS(mpParameters->Simulation()->Domain().present(), "Simulation/Domain");
00717         return mpParameters->Simulation()->Domain().get();
00718     }
00719     else
00720     {
00721         return mpParameters->ResumeSimulation()->Domain();
00722     }
00723 }
00724 
00725 cp::ionic_model_selection_type HeartConfig::GetDefaultIonicModel() const
00726 {
00727     CheckSimulationIsDefined("DefaultIonicModel");
00728 
00729     return mpParameters->Simulation()->IonicModels()->Default();
00730 }
00731 
00732 template<unsigned DIM>
00733 void HeartConfig::GetIonicModelRegions(std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& definedRegions,
00734                                        std::vector<cp::ionic_model_selection_type>& ionicModels) const
00735 {
00736     CheckSimulationIsDefined("IonicModelRegions");
00737     definedRegions.clear();
00738     ionicModels.clear();
00739 
00740     XSD_SEQUENCE_TYPE(cp::ionic_models_type::Region)&
00741          regions = mpParameters->Simulation()->IonicModels()->Region();
00742 
00743     for (XSD_ITERATOR_TYPE(cp::ionic_models_type::Region) i = regions.begin();
00744          i != regions.end();
00745          ++i)
00746     {
00747         cp::ionic_model_region_type ionic_model_region(*i);
00748 
00749         if (ionic_model_region.Location().Cuboid().present() || ionic_model_region.Location().Ellipsoid().present())
00750         {
00751             if (ionic_model_region.Location().Cuboid().present())
00752             {
00753                 cp::point_type point_a = ionic_model_region.Location().Cuboid()->LowerCoordinates();
00754                 cp::point_type point_b = ionic_model_region.Location().Cuboid()->UpperCoordinates();
00755 
00756                 switch (DIM)
00757                 {
00758                     case 1:
00759                     {
00760                         ChastePoint<DIM> chaste_point_a ( point_a.x() );
00761                         ChastePoint<DIM> chaste_point_b ( point_b.x() );
00762                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b )));
00763                         break;
00764                     }
00765                     case 2:
00766                     {
00767                         ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y() );
00768                         ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y() );
00769                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b )));
00770                         break;
00771                     }
00772                     case 3:
00773                     {
00774                         ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y(), point_a.z() );
00775                         ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y(), point_b.z() );
00776                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b )) );
00777                         break;
00778                     }
00779                     default:
00780                         NEVER_REACHED;
00781                         break;
00782                 }
00783             }
00784             else if (ionic_model_region.Location().Ellipsoid().present())
00785             {
00786                 cp::point_type centre = ionic_model_region.Location().Ellipsoid()->Centre();
00787                 cp::point_type radii  = ionic_model_region.Location().Ellipsoid()->Radii();
00788                 switch (DIM)
00789                 {
00790                     case 1:
00791                     {
00792                         ChastePoint<DIM> chaste_point_a ( centre.x() );
00793                         ChastePoint<DIM> chaste_point_b ( radii.x() );
00794                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b )) );
00795                         break;
00796                     }
00797                     case 2:
00798                     {
00799                         ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y() );
00800                         ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y() );
00801                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b )) );
00802                         break;
00803                     }
00804                     case 3:
00805                     {
00806                         ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y(), centre.z() );
00807                         ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y(), radii.z() );
00808                         definedRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b )) );
00809                         break;
00810                     }
00811                     default:
00812                     {
00813                         NEVER_REACHED;
00814                         break;
00815                     }
00816                 }
00817             }
00818             else
00819             {
00820                 NEVER_REACHED;
00821             }
00822 
00823             ionicModels.push_back(ionic_model_region.IonicModel());
00824         }
00825         else if(ionic_model_region.Location().EpiLayer().present() || ionic_model_region.Location().MidLayer().present() || ionic_model_region.Location().EndoLayer().present() )
00826         {
00828             EXCEPTION("Definition of transmural layers is not yet supported for defining different ionic models, please use cuboids instead");
00829         }
00830         else
00831         {
00832             EXCEPTION("Invalid region type for ionic model definition");
00833         }
00834     }
00835 }
00836 
00837 
00838 bool HeartConfig::IsMeshProvided() const
00839 {
00840     CheckSimulationIsDefined("Mesh");
00841     return mpParameters->Simulation()->Mesh().present();
00842 }
00843 
00844 bool HeartConfig::GetCreateMesh() const
00845 {
00846     CheckSimulationIsDefined("Mesh");
00847     CHECK_EXISTS(IsMeshProvided(), "Simulation/Mesh");
00848     cp::mesh_type mesh = mpParameters->Simulation()->Mesh().get();
00849     return (mesh.Slab().present() || mesh.Sheet().present() || mesh.Fibre().present());
00850 }
00851 
00852 bool HeartConfig::GetCreateSlab() const
00853 {
00854     CheckSimulationIsDefined("Mesh");
00855     CHECK_EXISTS(IsMeshProvided(), "Simulation/Mesh");
00856     cp::mesh_type mesh = mpParameters->Simulation()->Mesh().get();
00857     return (mesh.Slab().present());
00858 }
00859 
00860 bool HeartConfig::GetCreateSheet() const
00861 {
00862     CheckSimulationIsDefined("Mesh");
00863     CHECK_EXISTS(IsMeshProvided(), "Simulation/Mesh");
00864     cp::mesh_type mesh = mpParameters->Simulation()->Mesh().get();
00865     return (mesh.Sheet().present());
00866 }
00867 
00868 bool HeartConfig::GetCreateFibre() const
00869 {
00870     CheckSimulationIsDefined("Mesh");
00871     CHECK_EXISTS(IsMeshProvided(), "Simulation/Mesh");
00872     cp::mesh_type mesh = mpParameters->Simulation()->Mesh().get();
00873     return (mesh.Fibre().present());
00874 }
00875 
00876 
00877 bool HeartConfig::GetLoadMesh() const
00878 {
00879     CheckSimulationIsDefined("Mesh");
00880     CHECK_EXISTS(IsMeshProvided(), "Simulation/Mesh");
00881     return (mpParameters->Simulation()->Mesh()->LoadMesh().present());
00882 }
00883 
00884 void HeartConfig::GetSlabDimensions(c_vector<double, 3>& slabDimensions) const
00885 {
00886     CheckSimulationIsDefined("Slab");
00887 
00888     if (GetSpaceDimension()!=3 || !GetCreateSlab())
00889     {
00890         EXCEPTION("Tissue slabs can only be defined in 3D");
00891     }
00892 
00893     optional<cp::slab_type, false> slab_dimensions = mpParameters->Simulation()->Mesh()->Slab();
00894 
00895     slabDimensions[0] = slab_dimensions->x();
00896     slabDimensions[1] = slab_dimensions->y();
00897     slabDimensions[2] = slab_dimensions->z();
00898 }
00899 
00900 void HeartConfig::GetSheetDimensions(c_vector<double, 2>& sheetDimensions) const
00901 {
00902     CheckSimulationIsDefined("Sheet");
00903 
00904     if (GetSpaceDimension()!=2 || !GetCreateSheet())
00905     {
00906         EXCEPTION("Tissue sheets can only be defined in 2D");
00907     }
00908 
00909     optional<cp::sheet_type, false> sheet_dimensions = mpParameters->Simulation()->Mesh()->Sheet();
00910 
00911     sheetDimensions[0] = sheet_dimensions->x();
00912     sheetDimensions[1] = sheet_dimensions->y();
00913 }
00914 
00915 void HeartConfig::GetFibreLength(c_vector<double, 1>& fibreLength) const
00916 {
00917     CheckSimulationIsDefined("Fibre");
00918 
00919     if (GetSpaceDimension()!=1 || !GetCreateFibre())
00920     {
00921         EXCEPTION("Tissue fibres can only be defined in 1D");
00922     }
00923 
00924     optional<cp::fibre_type, false> fibre_length = mpParameters->Simulation()->Mesh()->Fibre();
00925 
00926     fibreLength[0] = fibre_length->x();
00927 }
00928 
00929 double HeartConfig::GetInterNodeSpace() const
00930 {
00931     CheckSimulationIsDefined("InterNodeSpace");
00932 
00933     switch (GetSpaceDimension())
00934     {
00935         case 3:
00936             CHECK_EXISTS(GetCreateSlab(), "Simulation/Mesh/Slab");
00937             return mpParameters->Simulation()->Mesh()->Slab()->inter_node_space();
00938             break;
00939         case 2:
00940             CHECK_EXISTS(GetCreateSheet(), "Simulation/Mesh/Sheet");
00941             return mpParameters->Simulation()->Mesh()->Sheet()->inter_node_space();
00942             break;
00943         case 1:
00944             CHECK_EXISTS(GetCreateFibre(), "Simulation/Mesh/Fibre");
00945             return mpParameters->Simulation()->Mesh()->Fibre()->inter_node_space();
00946             break;
00947         default:
00948             NEVER_REACHED;
00949 #define COVERAGE_IGNORE
00950             return 0.0; //To fool the compiler
00951 #undef COVERAGE_IGNORE
00952     }
00953 }
00954 
00955 std::string HeartConfig::GetMeshName() const
00956 {
00957     CheckSimulationIsDefined("LoadMesh");
00958     CHECK_EXISTS(GetLoadMesh(), "Mesh/LoadMesh");
00959 
00960     return mpParameters->Simulation()->Mesh()->LoadMesh()->name();
00961 }
00962 
00963 cp::media_type HeartConfig::GetConductivityMedia() const
00964 {
00965     CheckSimulationIsDefined("LoadMesh");
00966     CHECK_EXISTS(GetLoadMesh(), "Mesh/LoadMesh");
00967 
00968     return mpParameters->Simulation()->Mesh()->LoadMesh()->conductivity_media();
00969 }
00970 
00971 template <unsigned DIM>
00972 void HeartConfig::GetStimuli(std::vector<boost::shared_ptr<AbstractStimulusFunction> >& rStimuliApplied,
00973                              std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rStimulatedAreas) const
00974 {
00975     CheckSimulationIsDefined("Stimuli");
00976 
00977     if (!mpParameters->Simulation()->Stimuli().present())
00978     {
00979         // Finding no stimuli defined is allowed (although HeartConfigRelatedFactory does
00980         // throw an exception is no stimuli and no electrodes)
00981         return;
00982     }
00983 
00984     XSD_SEQUENCE_TYPE(cp::stimuli_type::Stimulus) stimuli = mpParameters->Simulation()->Stimuli()->Stimulus();
00985 
00986     for (XSD_ITERATOR_TYPE(cp::stimuli_type::Stimulus) i = stimuli.begin();
00987          i != stimuli.end();
00988          ++i)
00989     {
00990         cp::stimulus_type stimulus(*i);
00991         if (stimulus.Location().Cuboid().present() || stimulus.Location().Ellipsoid().present())
00992         {
00993             boost::shared_ptr<AbstractChasteRegion<DIM> > area_ptr;
00994             if (stimulus.Location().Cuboid().present() )
00995             {
00996                 cp::point_type point_a = stimulus.Location().Cuboid()->LowerCoordinates();
00997                 cp::point_type point_b = stimulus.Location().Cuboid()->UpperCoordinates();
00998                 switch (DIM)
00999                 {
01000                     case 1:
01001                     {
01002                         ChastePoint<DIM> chaste_point_a ( point_a.x() );
01003                         ChastePoint<DIM> chaste_point_b ( point_b.x() );
01004                         area_ptr.reset(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b ) );
01005                         break;
01006                     }
01007                     case 2:
01008                     {
01009                         ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y() );
01010                         ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y() );
01011                         area_ptr.reset(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b ) );
01012                         break;
01013                     }
01014                     case 3:
01015                     {
01016                         ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y(), point_a.z() );
01017                         ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y(), point_b.z() );
01018                         area_ptr.reset(new ChasteCuboid<DIM>( chaste_point_a, chaste_point_b ) );
01019                         break;
01020                     }
01021                     default:
01022                         NEVER_REACHED;
01023                         break;
01024                 }
01025             }
01026             else if (stimulus.Location().Ellipsoid().present())
01027             {
01028                 cp::point_type centre = stimulus.Location().Ellipsoid()->Centre();
01029                 cp::point_type radii  = stimulus.Location().Ellipsoid()->Radii();
01030                 switch (DIM)
01031                 {
01032                     case 1:
01033                     {
01034                         ChastePoint<DIM> chaste_point_a ( centre.x() );
01035                         ChastePoint<DIM> chaste_point_b ( radii.x() );
01036                         area_ptr.reset( new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b ) );
01037                         break;
01038                     }
01039                     case 2:
01040                     {
01041                         ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y() );
01042                         ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y() );
01043                         area_ptr.reset( new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b ) );
01044                         break;
01045                     }
01046                     case 3:
01047                     {
01048                         ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y(), centre.z() );
01049                         ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y(), radii.z() );
01050                         area_ptr.reset( new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b ) );
01051                         break;
01052                     }
01053                     default:
01054                     {
01055                         NEVER_REACHED;
01056                         break;
01057                     }
01058                 }
01059             }
01060             rStimulatedAreas.push_back(area_ptr);
01061 
01062             boost::shared_ptr<AbstractStimulusFunction> stim;
01063 
01064             if (stimulus.Period().present())
01065             {
01066                 if (stimulus.StopTime().present())
01067                 {
01068                     stim.reset(new RegularStimulus(stimulus.Strength(),
01069                                                    stimulus.Duration(),
01070                                                    stimulus.Period().get(),
01071                                                    stimulus.Delay(),
01072                                                    stimulus.StopTime().get()));
01073                 }
01074                 else
01075                 {
01076                     stim.reset(new RegularStimulus(stimulus.Strength(),
01077                                                    stimulus.Duration(),
01078                                                    stimulus.Period().get(),
01079                                                    stimulus.Delay()));
01080                 }
01081 
01082             }
01083             else
01084             {
01085                 if (stimulus.StopTime().present())
01086                 {
01087                     EXCEPTION("Stop time can not be defined for SimpleStimulus. Use Duration instead.");
01088                 }
01089 
01090                 stim.reset(new SimpleStimulus(stimulus.Strength(),
01091                                               stimulus.Duration(),
01092                                               stimulus.Delay()));
01093             }
01094             rStimuliApplied.push_back( stim );
01095         }
01096         else if(stimulus.Location().EpiLayer().present() || stimulus.Location().MidLayer().present() || stimulus.Location().EndoLayer().present() )
01097         {
01098             EXCEPTION("Definition of transmural layers is not yet supported for specifying stimulated areas, please use cuboids instead");
01099         }
01100         else
01101         {
01102             EXCEPTION("Invalid region type for stimulus definition");
01103         }
01104     }
01105 }
01106 
01107 
01108 template<unsigned DIM>
01109 void HeartConfig::GetCellHeterogeneities(std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rCellHeterogeneityRegions,
01110                                          std::vector<double>& rScaleFactorGks,
01111                                          std::vector<double>& rScaleFactorIto,
01112                                          std::vector<double>& rScaleFactorGkr,
01113                                          std::vector<std::map<std::string, double> >* pParameterSettings)
01114 {
01115     CheckSimulationIsDefined("CellHeterogeneities");
01116 
01117     if (!mpParameters->Simulation()->CellHeterogeneities().present())
01118     {
01119         // finding no heterogeneities defined is allowed
01120         return;
01121     }
01122     XSD_SEQUENCE_TYPE(cp::cell_heterogeneities_type::CellHeterogeneity) cell_heterogeneity
01123             = mpParameters->Simulation()->CellHeterogeneities()->CellHeterogeneity();
01124 
01125     bool user_supplied_negative_value = false;
01126     bool user_asking_for_transmural_layers = false;
01127     bool user_asked_for_cuboids_or_ellipsoids = false;
01128     unsigned counter_of_heterogeneities = 0;
01129 
01130     for (XSD_ITERATOR_TYPE(cp::cell_heterogeneities_type::CellHeterogeneity) i = cell_heterogeneity.begin();
01131          i != cell_heterogeneity.end();
01132          ++i)
01133     {
01134         cp::cell_heterogeneity_type ht(*i);
01135 
01136         if (ht.Location().Cuboid().present())
01137         {
01138             user_asked_for_cuboids_or_ellipsoids = true;
01139             cp::point_type point_a = ht.Location().Cuboid()->LowerCoordinates();
01140             cp::point_type point_b = ht.Location().Cuboid()->UpperCoordinates();
01141 
01142             ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y(), point_a.z() );
01143             ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y(), point_b.z() );
01144 
01145             rCellHeterogeneityRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteCuboid<DIM> ( chaste_point_a, chaste_point_b )) );
01146         }
01147         else if (ht.Location().Ellipsoid().present())
01148         {
01149             user_asked_for_cuboids_or_ellipsoids = true;
01150             cp::point_type centre = ht.Location().Ellipsoid()->Centre();
01151             cp::point_type radii  = ht.Location().Ellipsoid()->Radii();
01152 
01153             ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y(), centre.z() );
01154             ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y(), radii.z() );
01155             rCellHeterogeneityRegions.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >(new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b )) );
01156         }
01157         else if (ht.Location().EpiLayer().present())
01158         {
01159             mEpiFraction  =  ht.Location().EpiLayer().get();
01160 
01161             user_asking_for_transmural_layers = true;
01162             if (mEpiFraction <0)
01163             {
01164                 user_supplied_negative_value=true;
01165             }
01166             mIndexEpi = counter_of_heterogeneities;
01167         }
01168         else if (ht.Location().EndoLayer().present())
01169         {
01170             mEndoFraction  =  ht.Location().EndoLayer().get();
01171 
01172             user_asking_for_transmural_layers = true;
01173             if (mEndoFraction <0)
01174             {
01175                 user_supplied_negative_value=true;
01176             }
01177             mIndexEndo = counter_of_heterogeneities;
01178         }
01179         else if (ht.Location().MidLayer().present())
01180         {
01181             mMidFraction  =  ht.Location().MidLayer().get();
01182 
01183             user_asking_for_transmural_layers = true;
01184             if (mMidFraction <0)
01185             {
01186                 user_supplied_negative_value=true;
01187             }
01188             mIndexMid =  counter_of_heterogeneities;
01189         }
01190         else
01191         {
01192             EXCEPTION("Invalid region type for cell heterogeneity definition");
01193         }
01194 
01195         // Old scale factors
01196         rScaleFactorGks.push_back(ht.ScaleFactorGks().present() ? (double)ht.ScaleFactorGks().get() : 1.0);
01197         rScaleFactorIto.push_back(ht.ScaleFactorIto().present() ? (double)ht.ScaleFactorIto().get() : 1.0);
01198         rScaleFactorGkr.push_back(ht.ScaleFactorGkr().present() ? (double)ht.ScaleFactorGkr().get() : 1.0);
01199 
01200         // Named parameters
01201         if (pParameterSettings)
01202         {
01203             std::map<std::string, double> param_settings;
01204             XSD_SEQUENCE_TYPE(cp::cell_heterogeneity_type::SetParameter)& params = ht.SetParameter();
01205             for (XSD_ITERATOR_TYPE(cp::cell_heterogeneity_type::SetParameter) param_it = params.begin();
01206                  param_it != params.end();
01207                  ++param_it)
01208             {
01209                 cp::set_parameter_type param(*param_it);
01210                 param_settings[param.name()] = param.value();
01211             }
01212             pParameterSettings->push_back(param_settings);
01213         }
01214 
01215         counter_of_heterogeneities++;
01216     }
01217 
01218     //set the flag for request of transmural layers
01219     mUserAskedForCellularTransmuralHeterogeneities = user_asking_for_transmural_layers;
01220 
01221     if ( mUserAskedForCellularTransmuralHeterogeneities )
01222     {
01223         // cuboids/ellipsoids and layers at the same time are not yet supported
01224         if (user_asked_for_cuboids_or_ellipsoids )
01225         {
01226             EXCEPTION ("Specification of cellular heterogeneities by cuboids/ellipsoids and layers at the same time is not yet supported");
01227         }
01228 
01229         //check that the user supplied all three layers, the indexes should be 0, 1 and 2.
01230         // As they are initialised to a higher value, if their summation is higher than 3,
01231         // one (or more) is missing
01232         if ((mIndexMid+mIndexEndo+mIndexEpi) > 3)
01233         {
01234             EXCEPTION ("Three specifications of layers must be supplied");
01235         }
01236         if (fabs((mEndoFraction+mMidFraction+mEpiFraction)-1)>1e-2)
01237         {
01238             EXCEPTION ("Summation of epicardial, midmyocardial and endocardial fractions should be 1");
01239         }
01240         if (user_supplied_negative_value)
01241         {
01242            EXCEPTION ("Fractions must be positive");
01243         }
01244     }
01245 }
01246 
01247 bool HeartConfig::AreCellularTransmuralHeterogeneitiesRequested()
01248 {
01249     return mUserAskedForCellularTransmuralHeterogeneities;
01250 }
01251 
01252 double HeartConfig::GetEpiLayerFraction()
01253 {
01254     return mEpiFraction;
01255 }
01256 
01257 double HeartConfig::GetEndoLayerFraction()
01258 {
01259     return mEndoFraction;
01260 }
01261 
01262 double HeartConfig::GetMidLayerFraction()
01263 {
01264     return mMidFraction;
01265 }
01266 
01267 unsigned HeartConfig::GetEpiLayerIndex()
01268 {
01269     return mIndexEpi;
01270 }
01271 
01272 unsigned HeartConfig::GetEndoLayerIndex()
01273 {
01274     return mIndexEndo;
01275 }
01276 
01277 unsigned HeartConfig::GetMidLayerIndex()
01278 {
01279     return mIndexMid;
01280 }
01281 
01282 bool HeartConfig::GetConductivityHeterogeneitiesProvided() const
01283 {
01284     CheckSimulationIsDefined("ConductivityHeterogeneities");
01285     return mpParameters->Physiological().ConductivityHeterogeneities().present();
01286 }
01287 
01288 template<unsigned DIM>
01289 void HeartConfig::GetConductivityHeterogeneities(
01290         std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& conductivitiesHeterogeneityAreas,
01291         std::vector< c_vector<double,3> >& intraConductivities,
01292         std::vector< c_vector<double,3> >& extraConductivities) const
01293 {
01294     CheckSimulationIsDefined("ConductivityHeterogeneities");
01295     CHECK_EXISTS(GetConductivityHeterogeneitiesProvided(), "Physiological/ConductivityHeterogeneities");
01296     XSD_ANON_SEQUENCE_TYPE(cp::physiological_type, ConductivityHeterogeneities, ConductivityHeterogeneity)&
01297          conductivity_heterogeneity = mpParameters->Physiological().ConductivityHeterogeneities()->ConductivityHeterogeneity();
01298 
01299     for (XSD_ANON_ITERATOR_TYPE(cp::physiological_type, ConductivityHeterogeneities, ConductivityHeterogeneity) i = conductivity_heterogeneity.begin();
01300          i != conductivity_heterogeneity.end();
01301          ++i)
01302     {
01303         cp::conductivity_heterogeneity_type ht(*i);
01304 
01305         if (ht.Location().Cuboid().present())
01306         {
01307             cp::point_type point_a = ht.Location().Cuboid()->LowerCoordinates();
01308             cp::point_type point_b = ht.Location().Cuboid()->UpperCoordinates();
01309             ChastePoint<DIM> chaste_point_a ( point_a.x(), point_a.y(), point_a.z() );
01310             ChastePoint<DIM> chaste_point_b ( point_b.x(), point_b.y(), point_b.z() );
01311             conductivitiesHeterogeneityAreas.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >   (  new ChasteCuboid<DIM> ( chaste_point_a, chaste_point_b )) );
01312         }
01313         else if (ht.Location().Ellipsoid().present())
01314         {
01315             cp::point_type centre = ht.Location().Ellipsoid()->Centre();
01316             cp::point_type radii  = ht.Location().Ellipsoid()->Radii();
01317             ChastePoint<DIM> chaste_point_a ( centre.x(), centre.y(), centre.z() );
01318             ChastePoint<DIM> chaste_point_b ( radii.x(), radii.y(), radii.z() );
01319             conductivitiesHeterogeneityAreas.push_back(boost::shared_ptr<AbstractChasteRegion<DIM> >   (  new ChasteEllipsoid<DIM> ( chaste_point_a, chaste_point_b )) );
01320         }
01321         else if (ht.Location().EpiLayer().present() || ht.Location().MidLayer().present() || ht.Location().EndoLayer().present() )
01322         {
01324             EXCEPTION("Definition of transmural layers is not allowed for conductivities heterogeneities, you may use fibre orientation support instead");
01325         }
01326         else
01327         {
01328             EXCEPTION("Invalid region type for conductivity definition");
01329         }
01330 
01331         if (ht.IntracellularConductivities().present())
01332         {
01333             double intra_x = ht.IntracellularConductivities()->longi();
01334             double intra_y = ht.IntracellularConductivities()->trans();
01335             double intra_z = ht.IntracellularConductivities()->normal();
01336 
01337             intraConductivities.push_back( Create_c_vector(intra_x, intra_y, intra_z) );
01338         }
01339         else
01340         {
01341             c_vector<double, 3> intra_conductivities;
01342             GetIntracellularConductivities(intra_conductivities);
01343             intraConductivities.push_back(intra_conductivities);
01344         }
01345 
01346         if (ht.ExtracellularConductivities().present())
01347         {
01348             double extra_x = ht.ExtracellularConductivities()->longi();
01349             double extra_y = ht.ExtracellularConductivities()->trans();
01350             double extra_z = ht.ExtracellularConductivities()->normal();
01351 
01352             extraConductivities.push_back( Create_c_vector(extra_x, extra_y, extra_z) );
01353         }
01354         else
01355         {
01356             c_vector<double, 3> extra_conductivities;
01357             GetExtracellularConductivities(extra_conductivities);
01358             extraConductivities.push_back(extra_conductivities);
01359         }
01360 
01361     }
01362 }
01363 
01364 std::string HeartConfig::GetOutputDirectory() const
01365 {
01366     CheckSimulationIsDefined("Simulation/OutputDirectory");
01367     CHECK_EXISTS(mpParameters->Simulation()->OutputDirectory().present(), "Simulation/OutputDirectory");
01368     return mpParameters->Simulation()->OutputDirectory().get();
01369 }
01370 
01371 std::string HeartConfig::GetOutputFilenamePrefix() const
01372 {
01373     CheckSimulationIsDefined("Simulation/OutputFilenamePrefix");
01374     CHECK_EXISTS(mpParameters->Simulation()->OutputFilenamePrefix().present(), "Simulation/OutputFilenamePrefix");
01375     return mpParameters->Simulation()->OutputFilenamePrefix().get();
01376 }
01377 
01378 bool HeartConfig::GetOutputVariablesProvided() const
01379 {
01380     CheckSimulationIsDefined("OutputVariables");
01381     return mpParameters->Simulation()->OutputVariables().present();
01382 }
01383 
01384 void HeartConfig::GetOutputVariables(std::vector<std::string>& rOutputVariables) const
01385 {
01386     CHECK_EXISTS(GetOutputVariablesProvided(), "Simulation/OutputVariables");
01387     XSD_SEQUENCE_TYPE(cp::output_variables_type::Var)&
01388          output_variables = mpParameters->Simulation()->OutputVariables()->Var();
01389     rOutputVariables.clear();
01390 
01391     for (XSD_ITERATOR_TYPE(cp::output_variables_type::Var) i = output_variables.begin();
01392          i != output_variables.end();
01393          ++i)
01394     {
01395         cp::var_type& r_var(*i);
01396 
01397         // Add to outputVariables the string returned by var.name()
01398         rOutputVariables.push_back(r_var.name());
01399     }
01400 }
01401 
01402 bool HeartConfig::GetOutputUsingOriginalNodeOrdering()
01403 {
01404     CheckSimulationIsDefined("OutputUsingOriginalNodeOrdering");
01405     bool result = false;
01406     if (mpParameters->Simulation()->OutputUsingOriginalNodeOrdering().present())
01407     {
01408         result = (mpParameters->Simulation()->OutputUsingOriginalNodeOrdering().get() == cp::yesno_type::yes);
01409     }
01410     return result;
01411 }
01412 
01413 bool HeartConfig::GetCheckpointSimulation() const
01414 {
01415     return IsSimulationDefined() && mpParameters->Simulation()->CheckpointSimulation().present();
01416 }
01417 
01418 double HeartConfig::GetCheckpointTimestep() const
01419 {
01420     CHECK_EXISTS(GetCheckpointSimulation(), "Simulation/CheckpointSimulation");
01421     return mpParameters->Simulation()->CheckpointSimulation()->timestep();
01422 }
01423 
01424 unsigned HeartConfig::GetMaxCheckpointsOnDisk() const
01425 {
01426     CHECK_EXISTS(GetCheckpointSimulation(), "Simulation/CheckpointSimulation");
01427     return mpParameters->Simulation()->CheckpointSimulation()->max_checkpoints_on_disk();
01428 }
01429 
01430 
01431 HeartFileFinder HeartConfig::GetArchivedSimulationDir() const
01432 {
01433     CheckResumeSimulationIsDefined("GetArchivedSimulationDir");
01434 
01435     return HeartFileFinder(mpParameters->ResumeSimulation()->ArchiveDirectory());
01436 }
01437 
01438 
01439 void HeartConfig::GetIntracellularConductivities(c_vector<double, 3>& intraConductivities) const
01440 {
01441     CHECK_EXISTS(mpParameters->Physiological().IntracellularConductivities().present(), "Physiological/IntracellularConductivities");
01442     cp::conductivities_type intra_conductivities
01443         = mpParameters->Physiological().IntracellularConductivities().get();
01444     double intra_x_cond = intra_conductivities.longi();
01445     double intra_y_cond = intra_conductivities.trans();
01446     double intra_z_cond = intra_conductivities.normal();;
01447 
01448     assert(intra_y_cond != DBL_MAX);
01449     assert(intra_z_cond != DBL_MAX);
01450 
01451     intraConductivities[0] = intra_x_cond;
01452     intraConductivities[1] = intra_y_cond;
01453     intraConductivities[2] = intra_z_cond;
01454 }
01455 
01456 void HeartConfig::GetIntracellularConductivities(c_vector<double, 2>& intraConductivities) const
01457 {
01458     CHECK_EXISTS(mpParameters->Physiological().IntracellularConductivities().present(), "Physiological/IntracellularConductivities");
01459     cp::conductivities_type intra_conductivities
01460         = mpParameters->Physiological().IntracellularConductivities().get();
01461     double intra_x_cond = intra_conductivities.longi();
01462     double intra_y_cond = intra_conductivities.trans();
01463 
01464     assert(intra_y_cond != DBL_MAX);
01465 
01466     intraConductivities[0] = intra_x_cond;
01467     intraConductivities[1] = intra_y_cond;
01468 }
01469 
01470 void HeartConfig::GetIntracellularConductivities(c_vector<double, 1>& intraConductivities) const
01471 {
01472     CHECK_EXISTS(mpParameters->Physiological().IntracellularConductivities().present(), "Physiological/IntracellularConductivities");
01473     cp::conductivities_type intra_conductivities
01474         = mpParameters->Physiological().IntracellularConductivities().get();
01475     double intra_x_cond = intra_conductivities.longi();
01476 
01477     intraConductivities[0] = intra_x_cond;
01478 }
01479 
01480 void HeartConfig::GetExtracellularConductivities(c_vector<double, 3>& extraConductivities) const
01481 {
01482     CHECK_EXISTS(mpParameters->Physiological().ExtracellularConductivities().present(), "Physiological/ExtracellularConductivities");
01483     cp::conductivities_type extra_conductivities
01484         = mpParameters->Physiological().ExtracellularConductivities().get();
01485     double extra_x_cond = extra_conductivities.longi();
01486     double extra_y_cond = extra_conductivities.trans();
01487     double extra_z_cond = extra_conductivities.normal();;
01488 
01489     assert(extra_y_cond != DBL_MAX);
01490     assert(extra_z_cond != DBL_MAX);
01491 
01492     extraConductivities[0] = extra_x_cond;
01493     extraConductivities[1] = extra_y_cond;
01494     extraConductivities[2] = extra_z_cond;
01495 }
01496 
01497 void HeartConfig::GetExtracellularConductivities(c_vector<double, 2>& extraConductivities) const
01498 {
01499     CHECK_EXISTS(mpParameters->Physiological().ExtracellularConductivities().present(), "Physiological/ExtracellularConductivities");
01500     cp::conductivities_type extra_conductivities
01501         = mpParameters->Physiological().ExtracellularConductivities().get();
01502     double extra_x_cond = extra_conductivities.longi();
01503     double extra_y_cond = extra_conductivities.trans();
01504 
01505     assert(extra_y_cond != DBL_MAX);
01506 
01507     extraConductivities[0] = extra_x_cond;
01508     extraConductivities[1] = extra_y_cond;
01509 }
01510 
01511 void HeartConfig::GetExtracellularConductivities(c_vector<double, 1>& extraConductivities) const
01512 {
01513     CHECK_EXISTS(mpParameters->Physiological().ExtracellularConductivities().present(), "Physiological/ExtracellularConductivities");
01514     cp::conductivities_type extra_conductivities
01515         = mpParameters->Physiological().ExtracellularConductivities().get();
01516     double extra_x_cond = extra_conductivities.longi();
01517 
01518     extraConductivities[0] = extra_x_cond;
01519 }
01520 
01521 double HeartConfig::GetBathConductivity(unsigned bathRegion) const
01522 {
01523     /*
01524      *  We have to consider three cases: The user asks for ...
01525      *    a) ... the default conductivity (bathRegion=UINT_MAX)
01526      *    b) ... the conductivity of region defined to be heterogeneous
01527      *    c) ... the conductivity of region NOT defined to be heterogeneous
01528      *
01529      *  a) and c) should return the same
01530      */
01531 
01532     if (bathRegion == UINT_MAX)
01533     {
01534         /*bath conductivity mS/cm*/
01535         CHECK_EXISTS(mpParameters->Physiological().BathConductivity().present(), "Physiological/BathConductivity");
01536         return mpParameters->Physiological().BathConductivity().get();
01537     }
01538     else
01539     {
01540         assert(HeartRegionCode::IsRegionBath(bathRegion));
01541 
01542         std::map<unsigned, double>::const_iterator map_entry = mBathConductivities.find(bathRegion);
01543 
01544         if (map_entry != mBathConductivities.end())
01545         {
01546             return map_entry->second;
01547         }
01548         else
01549         {
01550             /*bath conductivity mS/cm*/
01551             CHECK_EXISTS(mpParameters->Physiological().BathConductivity().present(), "Physiological/BathConductivity");
01552             return mpParameters->Physiological().BathConductivity().get();
01553         }
01554     }
01555 }
01556 const std::set<unsigned>&  HeartConfig::rGetTissueIdentifiers()
01557 {
01558     return mTissueIdentifiers;
01559 }
01560 
01561 const std::set<unsigned>&  HeartConfig::rGetBathIdentifiers()
01562 {
01563     return mBathIdentifiers;
01564 }
01565 
01566 double HeartConfig::GetSurfaceAreaToVolumeRatio() const
01567 {
01568     CHECK_EXISTS(mpParameters->Physiological().SurfaceAreaToVolumeRatio().present(), "Physiological/SurfaceAreaToVolumeRatio");
01569     return mpParameters->Physiological().SurfaceAreaToVolumeRatio().get();
01570 }
01571 
01572 double HeartConfig::GetCapacitance() const
01573 {
01574     CHECK_EXISTS(mpParameters->Physiological().Capacitance().present(), "Physiological/Capacitance");
01575     return mpParameters->Physiological().Capacitance().get();
01576 }
01577 
01578 double HeartConfig::GetOdeTimeStep() const
01579 {
01580     CHECK_EXISTS(mpParameters->Numerical().TimeSteps().present(), "Numerical/TimeSteps");
01581     return mpParameters->Numerical().TimeSteps()->ode();
01582 }
01583 
01584 double HeartConfig::GetPdeTimeStep() const
01585 {
01586     CHECK_EXISTS(mpParameters->Numerical().TimeSteps().present(), "Numerical/TimeSteps");
01587     return mpParameters->Numerical().TimeSteps()->pde();
01588 }
01589 
01590 double HeartConfig::GetPrintingTimeStep() const
01591 {
01592     CHECK_EXISTS(mpParameters->Numerical().TimeSteps().present(), "Numerical/TimeSteps");
01593     return mpParameters->Numerical().TimeSteps()->printing();
01594 }
01595 
01596 bool HeartConfig::GetUseAbsoluteTolerance() const
01597 {
01598     CHECK_EXISTS(mpParameters->Numerical().KSPTolerances().present(), "Numerical/KSPTolerances");
01599     return mpParameters->Numerical().KSPTolerances()->KSPAbsolute().present();
01600 }
01601 
01602 double HeartConfig::GetAbsoluteTolerance() const
01603 {
01604     CHECK_EXISTS(mpParameters->Numerical().KSPTolerances().present(), "Numerical/KSPTolerances");
01605     if (!GetUseAbsoluteTolerance())
01606     {
01607         EXCEPTION("Absolute tolerance is not set in Chaste parameters");
01608     }
01609     return mpParameters->Numerical().KSPTolerances()->KSPAbsolute().get();
01610 }
01611 
01612 bool HeartConfig::GetUseRelativeTolerance() const
01613 {
01614     CHECK_EXISTS(mpParameters->Numerical().KSPTolerances().present(), "Numerical/KSPTolerances");
01615     return mpParameters->Numerical().KSPTolerances()->KSPRelative().present();
01616 }
01617 
01618 double HeartConfig::GetRelativeTolerance() const
01619 {
01620     CHECK_EXISTS(mpParameters->Numerical().KSPTolerances().present(), "Numerical/KSPTolerances");
01621     if (!GetUseRelativeTolerance())
01622     {
01623         EXCEPTION("Relative tolerance is not set in Chaste parameters");
01624     }
01625     return mpParameters->Numerical().KSPTolerances()->KSPRelative().get();
01626 }
01627 
01628 const char* HeartConfig::GetKSPSolver() const
01629 {
01630     CHECK_EXISTS(mpParameters->Numerical().KSPSolver().present(), "Numerical/KSPSolver");
01631     switch (mpParameters->Numerical().KSPSolver().get())
01632     {
01633         case cp::ksp_solver_type::gmres :
01634             return "gmres";
01635         case cp::ksp_solver_type::cg :
01636             return "cg";
01637         case cp::ksp_solver_type::symmlq :
01638             return "symmlq";
01639         case cp::ksp_solver_type::chebychev :
01640             return "chebychev";
01641     }
01642 #define COVERAGE_IGNORE
01643     EXCEPTION("Unknown ksp solver");
01644 #undef COVERAGE_IGNORE
01645 }
01646 
01647 const char* HeartConfig::GetKSPPreconditioner() const
01648 {
01649     CHECK_EXISTS(mpParameters->Numerical().KSPPreconditioner().present(), "Numerical/KSPPreconditioner");
01650     switch ( mpParameters->Numerical().KSPPreconditioner().get() )
01651     {
01652         case cp::ksp_preconditioner_type::jacobi :
01653             return "jacobi";
01654         case cp::ksp_preconditioner_type::bjacobi :
01655             return "bjacobi";
01656         case cp::ksp_preconditioner_type::hypre :
01657             return "hypre";
01658         case cp::ksp_preconditioner_type::ml :
01659             return "ml";
01660         case cp::ksp_preconditioner_type::spai :
01661             return "spai";
01662         case cp::ksp_preconditioner_type::blockdiagonal :
01663             return "blockdiagonal";
01664         case cp::ksp_preconditioner_type::ldufactorisation :
01665             return "ldufactorisation";
01666         case cp::ksp_preconditioner_type::twolevelsblockdiagonal :
01667             return "twolevelsblockdiagonal";
01668         case cp::ksp_preconditioner_type::none :
01669             return "none";
01670 
01671     }
01672 #define COVERAGE_IGNORE
01673     EXCEPTION("Unknown ksp preconditioner");
01674 #undef COVERAGE_IGNORE
01675 }
01676 
01677 DistributedTetrahedralMeshPartitionType::type HeartConfig::GetMeshPartitioning() const
01678 {
01679     CHECK_EXISTS(mpParameters->Numerical().MeshPartitioning().present(), "Numerical/MeshPartitioning");
01680     switch ( mpParameters->Numerical().MeshPartitioning().get() )
01681     {
01682         case cp::mesh_partitioning_type::dumb :
01683             return DistributedTetrahedralMeshPartitionType::DUMB;
01684         case cp::mesh_partitioning_type::metis :
01685             return DistributedTetrahedralMeshPartitionType::METIS_LIBRARY;
01686         case cp::mesh_partitioning_type::parmetis :
01687             return DistributedTetrahedralMeshPartitionType::PARMETIS_LIBRARY;
01688         case cp::mesh_partitioning_type::petsc :
01689             return DistributedTetrahedralMeshPartitionType::PETSC_MAT_PARTITION;
01690     }
01691 #define COVERAGE_IGNORE
01692     EXCEPTION("Unknown mesh partitioning type");
01693 #undef COVERAGE_IGNORE
01694 }
01695 
01696 bool HeartConfig::IsAdaptivityParametersPresent() const
01697 {
01698     bool IsAdaptivityParametersPresent = mpParameters->Numerical().AdaptivityParameters().present();
01699     if (IsAdaptivityParametersPresent)
01700     {
01701         WARNING("Use of the Adaptivity library is deprecated");
01702     }
01703     return IsAdaptivityParametersPresent;
01704 }
01705 
01706 /*
01707  * PostProcessing
01708  */
01709 
01710 bool HeartConfig::IsPostProcessingSectionPresent() const
01711 {
01712     return mpParameters->PostProcessing().present();
01713 }
01714 
01715 void HeartConfig::EnsurePostProcessingSectionPresent()
01716 {
01717     ENSURE_SECTION_PRESENT(mpParameters->PostProcessing(), cp::postprocessing_type);
01718 }
01719 
01720 bool HeartConfig::IsPostProcessingRequested() const
01721 {
01722     if (IsPostProcessingSectionPresent() == false)
01723     {
01724         return false;
01725     }
01726     else
01727     {
01728         return(IsApdMapsRequested() ||
01729                IsUpstrokeTimeMapsRequested() ||
01730                IsMaxUpstrokeVelocityMapRequested() ||
01731                IsConductionVelocityMapsRequested() ||
01732                IsAnyNodalTimeTraceRequested()||
01733                IsPseudoEcgCalculationRequested());
01734     }
01735 }
01736 bool HeartConfig::IsApdMapsRequested() const
01737 {
01738     bool result = false;
01739     if (IsPostProcessingSectionPresent())
01740     {
01741         XSD_SEQUENCE_TYPE(cp::postprocessing_type::ActionPotentialDurationMap)&
01742                 apd_maps = mpParameters->PostProcessing()->ActionPotentialDurationMap();
01743         result = (apd_maps.begin() != apd_maps.end());
01744     }
01745     return result;
01746 }
01747 
01748 void HeartConfig::GetApdMaps(std::vector<std::pair<double,double> >& apd_maps) const
01749 {
01750     CHECK_EXISTS(IsApdMapsRequested(), "PostProcessing/ActionPotentialDurationMap");
01751     apd_maps.clear();
01752 
01753     XSD_SEQUENCE_TYPE(cp::postprocessing_type::ActionPotentialDurationMap)&
01754         apd_maps_sequence = mpParameters->PostProcessing()->ActionPotentialDurationMap();
01755 
01756     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::ActionPotentialDurationMap) i = apd_maps_sequence.begin();
01757          i != apd_maps_sequence.end();
01758          ++i)
01759     {
01760         std::pair<double,double> map(i->repolarisation_percentage(),i->threshold());
01761 
01762         apd_maps.push_back(map);
01763     }
01764 }
01765 
01766 bool HeartConfig::IsUpstrokeTimeMapsRequested() const
01767 {
01768     bool result = false;
01769     if (IsPostProcessingSectionPresent())
01770     {
01771         XSD_SEQUENCE_TYPE(cp::postprocessing_type::UpstrokeTimeMap)&
01772             upstroke_map = mpParameters->PostProcessing()->UpstrokeTimeMap();
01773         result = (upstroke_map.begin() != upstroke_map.end());
01774     }
01775     return result;
01776 }
01777 void HeartConfig::GetUpstrokeTimeMaps (std::vector<double>& upstroke_time_maps) const
01778 {
01779     CHECK_EXISTS(IsUpstrokeTimeMapsRequested(), "PostProcessing/UpstrokeTimeMap");
01780     assert(upstroke_time_maps.size() == 0);
01781 
01782     XSD_SEQUENCE_TYPE(cp::postprocessing_type::UpstrokeTimeMap)&
01783         upstroke_maps_sequence = mpParameters->PostProcessing()->UpstrokeTimeMap();
01784 
01785     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::UpstrokeTimeMap) i = upstroke_maps_sequence.begin();
01786          i != upstroke_maps_sequence.end();
01787          ++i)
01788     {
01789         upstroke_time_maps.push_back(i->threshold());
01790     }
01791 }
01792 
01793 bool HeartConfig::IsMaxUpstrokeVelocityMapRequested() const
01794 {
01795     bool result = false;
01796     if (IsPostProcessingSectionPresent())
01797     {
01798         XSD_SEQUENCE_TYPE(cp::postprocessing_type::MaxUpstrokeVelocityMap)&
01799             max_upstroke_velocity_map = mpParameters->PostProcessing()->MaxUpstrokeVelocityMap();
01800         result = (max_upstroke_velocity_map.begin() != max_upstroke_velocity_map.end());
01801     }
01802     return result;
01803 }
01804 
01805 void HeartConfig::GetMaxUpstrokeVelocityMaps(std::vector<double>& upstroke_velocity_maps) const
01806 {
01807     CHECK_EXISTS(IsMaxUpstrokeVelocityMapRequested(), "PostProcessing/MaxUpstrokeVelocityMap");
01808     assert(upstroke_velocity_maps.size() == 0);
01809 
01810     XSD_SEQUENCE_TYPE(cp::postprocessing_type::MaxUpstrokeVelocityMap)&
01811         max_upstroke_velocity_maps_sequence = mpParameters->PostProcessing()->MaxUpstrokeVelocityMap();
01812 
01813     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::MaxUpstrokeVelocityMap) i = max_upstroke_velocity_maps_sequence.begin();
01814          i != max_upstroke_velocity_maps_sequence.end();
01815          ++i)
01816     {
01817         upstroke_velocity_maps.push_back(i->threshold());
01818     }
01819 }
01820 
01821 bool HeartConfig::IsConductionVelocityMapsRequested() const
01822 {
01823     bool result = false;
01824     if (IsPostProcessingSectionPresent())
01825     {
01826         XSD_SEQUENCE_TYPE(cp::postprocessing_type::ConductionVelocityMap)&
01827             cond_vel_maps = mpParameters->PostProcessing()->ConductionVelocityMap();
01828         result = (cond_vel_maps.begin() != cond_vel_maps.end());
01829     }
01830     return result;
01831 }
01832 
01833 void HeartConfig::GetConductionVelocityMaps(std::vector<unsigned>& conduction_velocity_maps) const
01834 {
01835     CHECK_EXISTS(IsConductionVelocityMapsRequested(), "PostProcessing/ConductionVelocityMap");
01836     assert(conduction_velocity_maps.size() == 0);
01837 
01838     XSD_SEQUENCE_TYPE(cp::postprocessing_type::ConductionVelocityMap)&
01839         cond_vel_maps_sequence = mpParameters->PostProcessing()->ConductionVelocityMap();
01840 
01841     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::ConductionVelocityMap) i = cond_vel_maps_sequence.begin();
01842          i != cond_vel_maps_sequence.end();
01843          ++i)
01844     {
01845         conduction_velocity_maps.push_back(i->origin_node());
01846     }
01847 }
01848 
01849 bool HeartConfig::IsAnyNodalTimeTraceRequested() const
01850 {
01851     bool result = false;
01852     if (IsPostProcessingSectionPresent())
01853     {
01854         XSD_SEQUENCE_TYPE(cp::postprocessing_type::TimeTraceAtNode)&
01855             requested_nodes = mpParameters->PostProcessing()->TimeTraceAtNode();
01856         result = (requested_nodes.begin() != requested_nodes.end());
01857     }
01858     return result;
01859 }
01860 
01861 void HeartConfig::GetNodalTimeTraceRequested(std::vector<unsigned>& rRequestedNodes) const
01862 {
01863     CHECK_EXISTS(IsAnyNodalTimeTraceRequested(), "PostProcessing/TimeTraceAtNode");
01864     assert(rRequestedNodes.size() == 0);
01865 
01866     XSD_SEQUENCE_TYPE(cp::postprocessing_type::TimeTraceAtNode)&
01867         req_nodes = mpParameters->PostProcessing()->TimeTraceAtNode();
01868 
01869     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::TimeTraceAtNode) i = req_nodes.begin();
01870          i != req_nodes.end();
01871          ++i)
01872     {
01873         rRequestedNodes.push_back(i->node_number());
01874     }
01875 }
01876 
01877 
01878 bool HeartConfig::IsPseudoEcgCalculationRequested() const
01879 {
01880     bool result = false;
01881     if (IsPostProcessingSectionPresent())
01882     {
01883         XSD_SEQUENCE_TYPE(cp::postprocessing_type::PseudoEcgElectrodePosition)&
01884             electrodes = mpParameters->PostProcessing()->PseudoEcgElectrodePosition();
01885         result = (electrodes.begin() != electrodes.end());
01886     }
01887     return result;
01888 }
01889 
01890 template<unsigned SPACE_DIM>
01891 void HeartConfig::GetPseudoEcgElectrodePositions(std::vector<ChastePoint<SPACE_DIM> >& rPseudoEcgElectrodePositions) const
01892 {
01893     rPseudoEcgElectrodePositions.clear();
01894     XSD_SEQUENCE_TYPE(cp::postprocessing_type::PseudoEcgElectrodePosition)&
01895         electrodes = mpParameters->PostProcessing()->PseudoEcgElectrodePosition();
01896     for (XSD_ITERATOR_TYPE(cp::postprocessing_type::PseudoEcgElectrodePosition) i = electrodes.begin();
01897          i != electrodes.end();
01898          ++i)
01899     {
01900         rPseudoEcgElectrodePositions.push_back(ChastePoint<SPACE_DIM>(i->x(), i->y(), i->z()));
01901     }
01902 }
01903 
01904 
01905 /*
01906  * Output visualization
01907  */
01908 
01909 bool HeartConfig::IsOutputVisualizerPresent() const
01910 {
01911     CheckSimulationIsDefined("OutputVisualizer");
01912 
01913     return mpParameters->Simulation()->OutputVisualizer().present();
01914 }
01915 
01916 bool HeartConfig::GetVisualizeWithMeshalyzer() const
01917 {
01918     if (!IsOutputVisualizerPresent())
01919     {
01920         return false;
01921     }
01922     else
01923     {
01924         return mpParameters->Simulation()->OutputVisualizer()->meshalyzer() == cp::yesno_type::yes;
01925     }
01926 }
01927 
01928 bool HeartConfig::GetVisualizeWithCmgui() const
01929 {
01930     if (!IsOutputVisualizerPresent())
01931     {
01932         return false;
01933     }
01934     else
01935     {
01936         return mpParameters->Simulation()->OutputVisualizer()->cmgui() == cp::yesno_type::yes;
01937     }
01938 }
01939 
01940 bool HeartConfig::GetVisualizeWithParallelVtk() const
01941 {
01942     if (!IsOutputVisualizerPresent())
01943     {
01944         return false;
01945     }
01946     else
01947     {
01948         return mpParameters->Simulation()->OutputVisualizer()->parallel_vtk() == cp::yesno_type::yes;
01949     }
01950 }
01951 
01952 bool HeartConfig::GetVisualizeWithVtk() const
01953 {
01954     if (!IsOutputVisualizerPresent())
01955     {
01956         return false;
01957     }
01958     else
01959     {
01960         return mpParameters->Simulation()->OutputVisualizer()->vtk() == cp::yesno_type::yes;
01961     }
01962 }
01963 
01964 unsigned HeartConfig::GetVisualizerOutputPrecision()
01965 {
01966     if (!IsOutputVisualizerPresent())
01967     {
01968         return 0u;
01969     }
01970     else
01971     {
01972         return mpParameters->Simulation()->OutputVisualizer()->precision();
01973     }
01974 }
01975 
01976 
01977 bool HeartConfig::IsElectrodesPresent() const
01978 {
01979     return mpParameters->Simulation()->Electrodes().present();
01980 }
01981 
01982 /*
01983  *  Set methods
01984  */
01985 void HeartConfig::SetSpaceDimension(unsigned spaceDimension)
01986 {
01987     mpParameters->Simulation()->SpaceDimension().set(spaceDimension);
01988 }
01989 
01990 void HeartConfig::SetSimulationDuration(double simulationDuration)
01991 {
01992     XSD_CREATE_WITH_FIXED_ATTR1(cp::time_type, time, simulationDuration, "ms");
01993     mpParameters->Simulation()->SimulationDuration().set(time);
01994 }
01995 
01996 void HeartConfig::SetDomain(const cp::domain_type& rDomain)
01997 {
01998     mpParameters->Simulation()->Domain().set(rDomain);
01999 }
02000 
02001 void HeartConfig::SetDefaultIonicModel(const cp::ionic_models_available_type& rIonicModel)
02002 {
02003     cp::ionic_model_selection_type ionic_model;
02004     ionic_model.Hardcoded(rIonicModel);
02005     cp::ionic_models_type container(ionic_model);
02006     mpParameters->Simulation()->IonicModels().set(container);
02007 }
02008 
02009 void HeartConfig::SetSlabDimensions(double x, double y, double z, double inter_node_space)
02010 {
02011     if ( ! mpParameters->Simulation()->Mesh().present())
02012     {
02013         XSD_CREATE_WITH_FIXED_ATTR(cp::mesh_type, mesh_to_load, "cm");
02014         mpParameters->Simulation()->Mesh().set(mesh_to_load);
02015     }
02016 
02017     cp::slab_type slab_definition(x, y, z, inter_node_space);
02018     mpParameters->Simulation()->Mesh()->Slab().set(slab_definition);
02019 }
02020 
02021 void HeartConfig::SetSheetDimensions(double x, double y, double inter_node_space)
02022 {
02023     if ( ! mpParameters->Simulation()->Mesh().present())
02024     {
02025         XSD_CREATE_WITH_FIXED_ATTR(cp::mesh_type, mesh_to_load, "cm");
02026         mpParameters->Simulation()->Mesh().set(mesh_to_load);
02027     }
02028 
02029     cp::sheet_type sheet_definition(x, y, inter_node_space);
02030     mpParameters->Simulation()->Mesh()->Sheet().set(sheet_definition);
02031 }
02032 
02033 void HeartConfig::SetFibreLength(double x, double inter_node_space)
02034 {
02035     if ( ! mpParameters->Simulation()->Mesh().present())
02036     {
02037         XSD_CREATE_WITH_FIXED_ATTR(cp::mesh_type, mesh_to_load, "cm");
02038         mpParameters->Simulation()->Mesh().set(mesh_to_load);
02039     }
02040 
02041     cp::fibre_type fibre_definition(x, inter_node_space);
02042     mpParameters->Simulation()->Mesh()->Fibre().set(fibre_definition);
02043 }
02044 
02045 void HeartConfig::SetMeshFileName(std::string meshPrefix, cp::media_type fibreDefinition)
02046 {
02047     if ( ! mpParameters->Simulation()->Mesh().present())
02048     {
02049         XSD_CREATE_WITH_FIXED_ATTR(cp::mesh_type, mesh_to_load, "cm");
02050         mpParameters->Simulation()->Mesh().set(mesh_to_load);
02051     }
02052 
02053     XSD_NESTED_TYPE(cp::mesh_type::LoadMesh) mesh_prefix(meshPrefix, fibreDefinition);
02054     mpParameters->Simulation()->Mesh()->LoadMesh().set(mesh_prefix);
02055 }
02056 
02057 void HeartConfig::SetIonicModelRegions(std::vector<ChasteCuboid<3> >& rDefinedRegions,
02058                                        std::vector<cp::ionic_model_selection_type>& rIonicModels) const
02059 {
02060     assert(rDefinedRegions.size() == rIonicModels.size());
02061     // You need to have defined a default model first...
02062     assert(mpParameters->Simulation()->IonicModels().present());
02063     XSD_SEQUENCE_TYPE(cp::ionic_models_type::Region)&
02064         regions = mpParameters->Simulation()->IonicModels()->Region();
02065     regions.clear();
02066     for (unsigned region_index=0; region_index<rDefinedRegions.size(); region_index++)
02067     {
02068         cp::point_type point_a(rDefinedRegions[region_index].rGetLowerCorner()[0],
02069                                rDefinedRegions[region_index].rGetLowerCorner()[1],
02070                                rDefinedRegions[region_index].rGetLowerCorner()[2]);
02071 
02072         cp::point_type point_b(rDefinedRegions[region_index].rGetUpperCorner()[0],
02073                                rDefinedRegions[region_index].rGetUpperCorner()[1],
02074                                rDefinedRegions[region_index].rGetUpperCorner()[2]);
02075 
02076         XSD_CREATE_WITH_FIXED_ATTR(cp::location_type, locn, "cm");
02077         locn.Cuboid().set(cp::box_type(point_a, point_b));
02078 
02079         cp::ionic_model_region_type region(rIonicModels[region_index], locn);
02080         regions.push_back(region);
02081     }
02082 }
02083 
02084 void HeartConfig::SetConductivityHeterogeneities(std::vector<ChasteCuboid<3> >& conductivityAreas,
02085         std::vector< c_vector<double,3> >& intraConductivities,
02086         std::vector< c_vector<double,3> >& extraConductivities)
02087 {
02088     assert ( conductivityAreas.size() == intraConductivities.size() );
02089     assert ( intraConductivities.size() == extraConductivities.size());
02090 
02091     XSD_ANON_SEQUENCE_TYPE(cp::physiological_type, ConductivityHeterogeneities, ConductivityHeterogeneity) heterogeneities_container;
02092 
02093     for (unsigned region_index=0; region_index<conductivityAreas.size(); region_index++)
02094     {
02095         cp::point_type point_a(conductivityAreas[region_index].rGetLowerCorner()[0],
02096                            conductivityAreas[region_index].rGetLowerCorner()[1],
02097                            conductivityAreas[region_index].rGetLowerCorner()[2]);
02098 
02099         cp::point_type point_b(conductivityAreas[region_index].rGetUpperCorner()[0],
02100                            conductivityAreas[region_index].rGetUpperCorner()[1],
02101                            conductivityAreas[region_index].rGetUpperCorner()[2]);
02102 
02103         XSD_CREATE_WITH_FIXED_ATTR(cp::location_type, locn, "cm");
02104         locn.Cuboid().set(cp::box_type(point_a, point_b));
02105         cp::conductivity_heterogeneity_type ht(locn);
02106 
02107         XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, intra,
02108                                     intraConductivities[region_index][0],
02109                                     intraConductivities[region_index][1],
02110                                     intraConductivities[region_index][2],
02111                                     "mS/cm");
02112 
02113         ht.IntracellularConductivities(intra);
02114 
02115         XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, extra,
02116                                     extraConductivities[region_index][0],
02117                                     extraConductivities[region_index][1],
02118                                     extraConductivities[region_index][2],
02119                                     "mS/cm");
02120 
02121         ht.ExtracellularConductivities(extra);
02122 
02123         heterogeneities_container.push_back(ht);
02124     }
02125 
02126     XSD_ANON_TYPE(cp::physiological_type, ConductivityHeterogeneities) heterogeneities_object;
02127     heterogeneities_object.ConductivityHeterogeneity(heterogeneities_container);
02128 
02129     mpParameters->Physiological().ConductivityHeterogeneities().set(heterogeneities_object);
02130 }
02131 
02132 void HeartConfig::SetConductivityHeterogeneitiesEllipsoid(std::vector<ChasteEllipsoid<3> >& conductivityAreas,
02133         std::vector< c_vector<double,3> >& intraConductivities,
02134         std::vector< c_vector<double,3> >& extraConductivities)
02135 {
02136     assert ( conductivityAreas.size() == intraConductivities.size() );
02137     assert ( intraConductivities.size() == extraConductivities.size());
02138 
02139     XSD_ANON_SEQUENCE_TYPE(cp::physiological_type, ConductivityHeterogeneities, ConductivityHeterogeneity) heterogeneities_container;
02140 
02141     for (unsigned region_index=0; region_index<conductivityAreas.size(); region_index++)
02142     {
02143         cp::point_type centre(conductivityAreas[region_index].rGetCentre()[0],
02144                               conductivityAreas[region_index].rGetCentre()[1],
02145                               conductivityAreas[region_index].rGetCentre()[2]);
02146 
02147         cp::point_type radii(conductivityAreas[region_index].rGetRadii()[0],
02148                              conductivityAreas[region_index].rGetRadii()[1],
02149                              conductivityAreas[region_index].rGetRadii()[2]);
02150 
02151         XSD_CREATE_WITH_FIXED_ATTR(cp::location_type, locn, "cm");
02152         locn.Ellipsoid().set(cp::ellipsoid_type(centre, radii));
02153         cp::conductivity_heterogeneity_type ht(locn);
02154 
02155         XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, intra,
02156                                     intraConductivities[region_index][0],
02157                                     intraConductivities[region_index][1],
02158                                     intraConductivities[region_index][2],
02159                                     "mS/cm");
02160 
02161         ht.IntracellularConductivities(intra);
02162 
02163         XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, extra,
02164                                     extraConductivities[region_index][0],
02165                                     extraConductivities[region_index][1],
02166                                     extraConductivities[region_index][2],
02167                                     "mS/cm");
02168 
02169         ht.ExtracellularConductivities(extra);
02170 
02171         heterogeneities_container.push_back(ht);
02172     }
02173 
02174     XSD_ANON_TYPE(cp::physiological_type, ConductivityHeterogeneities) heterogeneities_object;
02175     heterogeneities_object.ConductivityHeterogeneity(heterogeneities_container);
02176 
02177     mpParameters->Physiological().ConductivityHeterogeneities().set(heterogeneities_object);
02178 }
02179 
02180 void HeartConfig::SetOutputDirectory(const std::string& rOutputDirectory)
02181 {
02182     mpParameters->Simulation()->OutputDirectory().set(rOutputDirectory);
02183 }
02184 
02185 void HeartConfig::SetOutputFilenamePrefix(const std::string& rOutputFilenamePrefix)
02186 {
02187     mpParameters->Simulation()->OutputFilenamePrefix().set(rOutputFilenamePrefix);
02188 }
02189 
02190 void HeartConfig::SetOutputVariables(const std::vector<std::string>& rOutputVariables)
02191 {
02192     if ( ! mpParameters->Simulation()->OutputVariables().present())
02193     {
02194         cp::output_variables_type variables_requested;
02195         mpParameters->Simulation()->OutputVariables().set(variables_requested);
02196     }
02197 
02198     XSD_SEQUENCE_TYPE(cp::output_variables_type::Var)&
02199         var_type_sequence = mpParameters->Simulation()->OutputVariables()->Var();
02200     //Erase or create a sequence
02201     var_type_sequence.clear();
02202 
02203     for (unsigned i=0; i<rOutputVariables.size(); i++)
02204     {
02205         cp::var_type temp(rOutputVariables[i]);
02206         var_type_sequence.push_back(temp);
02207     }
02208 }
02209 
02210 void  HeartConfig::SetOutputUsingOriginalNodeOrdering(bool useOriginal)
02211 {
02212     //What if it doesn't exist?
02213     mpParameters->Simulation()->OutputUsingOriginalNodeOrdering().set(useOriginal? cp::yesno_type::yes : cp::yesno_type::no);
02214 }
02215 
02216 void HeartConfig::SetCheckpointSimulation(bool saveSimulation, double checkpointTimestep, unsigned maxCheckpointsOnDisk)
02217 {
02218     if (saveSimulation)
02219     {
02220         // Make sure values for the optional parameters have been provided
02221         assert(checkpointTimestep!=-1.0 && maxCheckpointsOnDisk!=UINT_MAX);
02222 
02223         XSD_CREATE_WITH_FIXED_ATTR2(cp::simulation_type::XSD_NESTED_TYPE(CheckpointSimulation),
02224                                     cs,
02225                                     checkpointTimestep,
02226                                     maxCheckpointsOnDisk,
02227                                     "ms");
02228         mpParameters->Simulation()->CheckpointSimulation().set(cs);
02229     }
02230     else
02231     {
02232         mpParameters->Simulation()->CheckpointSimulation().reset();
02233     }
02234 
02235     CheckTimeSteps();
02236 }
02237 
02238 // Physiological
02239 
02240 void HeartConfig::SetIntracellularConductivities(const c_vector<double, 3>& intraConductivities)
02241 {
02242     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, intra,
02243                                 intraConductivities[0],
02244                                 intraConductivities[1],
02245                                 intraConductivities[2],
02246                                 "mS/cm");
02247 
02248     mpParameters->Physiological().IntracellularConductivities().set(intra);
02249 }
02250 
02251 void HeartConfig::SetIntracellularConductivities(const c_vector<double, 2>& intraConductivities)
02252 {
02253     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, intra,
02254                                 intraConductivities[0],
02255                                 intraConductivities[1],
02256                                 0.0, "mS/cm");
02257 
02258     mpParameters->Physiological().IntracellularConductivities().set(intra);
02259 }
02260 
02261 void HeartConfig::SetIntracellularConductivities(const c_vector<double, 1>& intraConductivities)
02262 {
02263     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, intra,
02264                                 intraConductivities[0],
02265                                 0.0, 0.0, "mS/cm");
02266 
02267     mpParameters->Physiological().IntracellularConductivities().set(intra);
02268 }
02269 
02270 void HeartConfig::SetExtracellularConductivities(const c_vector<double, 3>& extraConductivities)
02271 {
02272     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, extra,
02273                                 extraConductivities[0],
02274                                 extraConductivities[1],
02275                                 extraConductivities[2],
02276                                 "mS/cm");
02277 
02278     mpParameters->Physiological().ExtracellularConductivities().set(extra);
02279 }
02280 
02281 void HeartConfig::SetExtracellularConductivities(const c_vector<double, 2>& extraConductivities)
02282 {
02283     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, extra,
02284                                 extraConductivities[0],
02285                                 extraConductivities[1],
02286                                 0.0, "mS/cm");
02287 
02288     mpParameters->Physiological().ExtracellularConductivities().set(extra);
02289 }
02290 
02291 void HeartConfig::SetExtracellularConductivities(const c_vector<double, 1>& extraConductivities)
02292 {
02293     XSD_CREATE_WITH_FIXED_ATTR3(cp::conductivities_type, extra,
02294                                 extraConductivities[0],
02295                                 0.0, 0.0, "mS/cm");
02296 
02297     mpParameters->Physiological().ExtracellularConductivities().set(extra);
02298 }
02299 
02300 void HeartConfig::SetBathConductivity(double bathConductivity)
02301 {
02302     XSD_CREATE_WITH_FIXED_ATTR1(cp::conductivity_type, cond, bathConductivity, "mS/cm");
02303     mpParameters->Physiological().BathConductivity().set(cond);
02304 }
02305 
02306 void HeartConfig::SetBathMultipleConductivities(std::map<unsigned, double> bathConductivities)
02307 {
02309     mBathConductivities = bathConductivities;
02310 }
02311 
02312 //void HeartConfig::SetTissueIdentifiers(const std::set<unsigned>& tissueIds)
02313 //{
02314 //    std::set<unsigned> empty_bath_identifiers;  //Too dangerous (see GetValidBathId)
02315 //    SetTissueAndBathIdentifiers(tissueIds, mBathIdentifiers);
02316 //}
02317 
02318 void HeartConfig::SetTissueAndBathIdentifiers(const std::set<unsigned>& tissueIds, const std::set<unsigned>& bathIds)
02319 {
02320     if (tissueIds.empty() || bathIds.empty() )
02321     {
02322         EXCEPTION("Identifying set must be non-empty");
02323     }
02324     std::set<unsigned> shared_identifiers;
02325     std::set_intersection(tissueIds.begin(),
02326                           tissueIds.end(),
02327                           bathIds.begin(),
02328                           bathIds.end(),
02329                           std::inserter(shared_identifiers, shared_identifiers.begin()));
02330 
02331     if (!shared_identifiers.empty())
02332     {
02333         EXCEPTION("Tissue identifiers and bath identifiers overlap");
02334     }
02335     mTissueIdentifiers=tissueIds;
02336     mBathIdentifiers=bathIds;
02337 }
02338 
02339 void HeartConfig::SetSurfaceAreaToVolumeRatio(double ratio)
02340 {
02341     XSD_CREATE_WITH_FIXED_ATTR1(cp::inverse_length_type, ratio_object, ratio, "1/cm");
02342     mpParameters->Physiological().SurfaceAreaToVolumeRatio().set(ratio_object);
02343 }
02344 
02345 void HeartConfig::SetCapacitance(double capacitance)
02346 {
02347     XSD_CREATE_WITH_FIXED_ATTR1(cp::capacitance_type, capacitance_object, capacitance, "uF/cm^2");
02348     mpParameters->Physiological().Capacitance().set(capacitance_object);
02349 }
02350 
02351 
02352 // Numerical
02353 void HeartConfig::SetOdePdeAndPrintingTimeSteps(double odeTimeStep, double pdeTimeStep, double printingTimeStep)
02354 {
02355     XSD_CREATE_WITH_FIXED_ATTR3(cp::time_steps_type, time_steps,
02356                                 odeTimeStep, pdeTimeStep, printingTimeStep, "ms");
02357     mpParameters->Numerical().TimeSteps().set(time_steps);
02358     CheckTimeSteps();
02359 }
02360 
02361 void HeartConfig::SetOdeTimeStep(double odeTimeStep)
02362 {
02363     SetOdePdeAndPrintingTimeSteps(odeTimeStep, GetPdeTimeStep(), GetPrintingTimeStep());
02364 }
02365 
02366 void HeartConfig::SetPdeTimeStep(double pdeTimeStep)
02367 {
02368     SetOdePdeAndPrintingTimeSteps(GetOdeTimeStep(), pdeTimeStep, GetPrintingTimeStep());
02369 }
02370 
02371 void HeartConfig::SetPrintingTimeStep(double printingTimeStep)
02372 {
02373     SetOdePdeAndPrintingTimeSteps(GetOdeTimeStep(), GetPdeTimeStep(), printingTimeStep);
02374 }
02375 
02376 void HeartConfig::CheckTimeSteps() const
02377 {
02378     if (GetOdeTimeStep() <= 0)
02379     {
02380         EXCEPTION("Ode time-step should be positive");
02381     }
02382     if (GetPdeTimeStep() <= 0)
02383     {
02384         EXCEPTION("Pde time-step should be positive");
02385     }
02386     if (GetPrintingTimeStep() <= 0.0)
02387     {
02388         EXCEPTION("Printing time-step should be positive");
02389     }
02390 
02391     if (GetPdeTimeStep()>GetPrintingTimeStep())
02392     {
02393         EXCEPTION("Printing time-step should not be smaller than PDE time-step");
02394     }
02395 
02396     if ( !Divides(GetPdeTimeStep(), GetPrintingTimeStep()) )
02397     {
02398         EXCEPTION("Printing time-step should be a multiple of PDE time step");
02399     }
02400 
02401     if ( GetOdeTimeStep() > GetPdeTimeStep() )
02402     {
02403         EXCEPTION("Ode time-step should not be greater than PDE time-step");
02404     }
02405 
02406     if (GetCheckpointSimulation())
02407     {
02408         if (GetCheckpointTimestep() <= 0.0)
02409         {
02410             EXCEPTION("Checkpoint time-step should be positive");
02411         }
02412 
02413         if ( !Divides(GetPrintingTimeStep(), GetCheckpointTimestep()) )
02414         {
02415             EXCEPTION("Checkpoint time-step should be a multiple of printing time-step");
02416         }
02417     }
02418 }
02419 
02420 
02421 void HeartConfig::SetUseRelativeTolerance(double relativeTolerance)
02422 {
02423     ENSURE_SECTION_PRESENT(mpParameters->Numerical().KSPTolerances(), cp::ksp_tolerances_type);
02424     //Remove any reference to tolerances is user parameters
02425     mpParameters->Numerical().KSPTolerances()->KSPAbsolute().reset();
02426     mpParameters->Numerical().KSPTolerances()->KSPRelative().set(relativeTolerance);
02427 }
02428 
02429 void HeartConfig::SetUseAbsoluteTolerance(double absoluteTolerance)
02430 {
02431     ENSURE_SECTION_PRESENT(mpParameters->Numerical().KSPTolerances(), cp::ksp_tolerances_type);
02432     //Remove any reference to tolerances is user parameters
02433     mpParameters->Numerical().KSPTolerances()->KSPRelative().reset();
02434     mpParameters->Numerical().KSPTolerances()->KSPAbsolute().set(absoluteTolerance);
02435 }
02436 
02437 void HeartConfig::SetKSPSolver(const char* kspSolver, bool warnOfChange)
02438 {
02439     if (warnOfChange && strcmp(GetKSPSolver(), kspSolver) != 0)
02440     {
02441         //Warn
02442         WARNING("Code has changed the KSP solver type from "<<GetKSPSolver()<<" to "<< kspSolver);
02443     }
02444 
02445     /* Note that changes in these conditions need to be reflected in the Doxygen*/
02446     if ( strcmp(kspSolver, "gmres") == 0)
02447     {
02448         mpParameters->Numerical().KSPSolver().set(cp::ksp_solver_type::gmres);
02449         return;
02450     }
02451     if ( strcmp(kspSolver, "cg") == 0)
02452     {
02453         mpParameters->Numerical().KSPSolver().set(cp::ksp_solver_type::cg);
02454         return;
02455     }
02456     if ( strcmp(kspSolver, "symmlq") == 0)
02457     {
02458         mpParameters->Numerical().KSPSolver().set(cp::ksp_solver_type::symmlq);
02459         return;
02460     }
02461     if ( strcmp(kspSolver, "chebychev") == 0)
02462     {
02463         mpParameters->Numerical().KSPSolver().set(cp::ksp_solver_type::chebychev);
02464         return;
02465     }
02466 
02467     EXCEPTION("Unknown solver type provided");
02468 }
02469 
02470 void HeartConfig::SetKSPPreconditioner(const char* kspPreconditioner)
02471 {
02472     /* Note that changes in these conditions need to be reflected in the Doxygen*/
02473     if ( strcmp(kspPreconditioner, "jacobi") == 0)
02474     {
02475         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::jacobi);
02476         return;
02477     }
02478     if ( strcmp(kspPreconditioner, "bjacobi") == 0)
02479     {
02480         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::bjacobi);
02481         return;
02482     }
02483     if ( strcmp(kspPreconditioner, "hypre") == 0)
02484     {
02485         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::hypre);
02486         return;
02487     }
02488     if ( strcmp(kspPreconditioner, "ml") == 0)
02489     {
02490         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::ml);
02491         return;
02492     }
02493     if ( strcmp(kspPreconditioner, "spai") == 0)
02494     {
02495         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::spai);
02496         return;
02497     }
02498     if ( strcmp(kspPreconditioner, "twolevelsblockdiagonal") == 0)
02499     {
02500         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::twolevelsblockdiagonal);
02501         return;
02502     }
02503     if ( strcmp(kspPreconditioner, "blockdiagonal") == 0)
02504     {
02505         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::blockdiagonal);
02506         return;
02507     }
02508     if ( strcmp(kspPreconditioner, "ldufactorisation") == 0)
02509     {
02510         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::ldufactorisation);
02511         return;
02512     }
02513     if ( strcmp(kspPreconditioner, "none") == 0)
02514     {
02515         mpParameters->Numerical().KSPPreconditioner().set(cp::ksp_preconditioner_type::none);
02516         return;
02517     }
02518 
02519     EXCEPTION("Unknown preconditioner type provided");
02520 }
02521 
02522 void HeartConfig::SetMeshPartitioning(const char* meshPartioningMethod)
02523 {
02524     /* Note that changes in these conditions need to be reflected in the Doxygen*/
02525     if ( strcmp(meshPartioningMethod, "dumb") == 0)
02526     {
02527         mpParameters->Numerical().MeshPartitioning().set(cp::mesh_partitioning_type::dumb);
02528         return;
02529     }
02530     if ( strcmp(meshPartioningMethod, "metis") == 0)
02531     {
02532         mpParameters->Numerical().MeshPartitioning().set(cp::mesh_partitioning_type::metis);
02533         return;
02534     }
02535     if ( strcmp(meshPartioningMethod, "parmetis") == 0)
02536     {
02537         mpParameters->Numerical().MeshPartitioning().set(cp::mesh_partitioning_type::parmetis);
02538         return;
02539     }
02540     if ( strcmp(meshPartioningMethod, "petsc") == 0)
02541     {
02542         mpParameters->Numerical().MeshPartitioning().set(cp::mesh_partitioning_type::petsc);
02543         return;
02544     }
02545 
02546     EXCEPTION("Unknown mesh partitioning method provided");
02547 }
02548 
02549 
02550 void HeartConfig::SetApdMaps(const std::vector<std::pair<double,double> >& apdMaps)
02551 {
02552     EnsurePostProcessingSectionPresent();
02553     XSD_SEQUENCE_TYPE(cp::postprocessing_type::ActionPotentialDurationMap)& apd_maps_sequence
02554         = mpParameters->PostProcessing()->ActionPotentialDurationMap();
02555     //Erase or create a sequence
02556     apd_maps_sequence.clear();
02557 
02558     for (unsigned i=0; i<apdMaps.size(); i++)
02559     {
02560         XSD_CREATE_WITH_FIXED_ATTR2(cp::apd_map_type, temp,
02561                                     apdMaps[i].first, apdMaps[i].second,
02562                                     "mV");
02563         apd_maps_sequence.push_back( temp);
02564     }
02565 }
02566 
02567 
02568 void HeartConfig::SetUpstrokeTimeMaps (std::vector<double>& upstrokeTimeMaps)
02569 {
02570     EnsurePostProcessingSectionPresent();
02571     XSD_SEQUENCE_TYPE(cp::postprocessing_type::UpstrokeTimeMap)& var_type_sequence
02572         = mpParameters->PostProcessing()->UpstrokeTimeMap();
02573 
02574     //Erase or create a sequence
02575     var_type_sequence.clear();
02576 
02577     for (unsigned i=0; i<upstrokeTimeMaps.size(); i++)
02578     {
02579         XSD_CREATE_WITH_FIXED_ATTR1(cp::upstrokes_map_type, temp,
02580                                     upstrokeTimeMaps[i],
02581                                     "mV");
02582         var_type_sequence.push_back(temp);
02583     }
02584 }
02585 
02586 void HeartConfig::SetMaxUpstrokeVelocityMaps (std::vector<double>& maxUpstrokeVelocityMaps)
02587 {
02588     EnsurePostProcessingSectionPresent();
02589     XSD_SEQUENCE_TYPE(cp::postprocessing_type::MaxUpstrokeVelocityMap)& max_upstroke_velocity_maps_sequence
02590         = mpParameters->PostProcessing()->MaxUpstrokeVelocityMap();
02591 
02592     //Erase or create a sequence
02593     max_upstroke_velocity_maps_sequence.clear();
02594 
02595     for (unsigned i=0; i<maxUpstrokeVelocityMaps.size(); i++)
02596     {
02597         XSD_CREATE_WITH_FIXED_ATTR1(cp::max_upstrokes_velocity_map_type, temp,
02598                                     maxUpstrokeVelocityMaps[i],
02599                                     "mV");
02600 
02601 
02602         max_upstroke_velocity_maps_sequence.push_back(temp);
02603     }
02604 
02605 }
02606 
02607 void HeartConfig::SetConductionVelocityMaps (std::vector<unsigned>& conductionVelocityMaps)
02608 {
02609     EnsurePostProcessingSectionPresent();
02610     XSD_SEQUENCE_TYPE(cp::postprocessing_type::ConductionVelocityMap)& conduction_velocity_maps_sequence
02611         = mpParameters->PostProcessing()->ConductionVelocityMap();
02612 
02613     //Erase or create a sequence
02614     conduction_velocity_maps_sequence.clear();
02615 
02616     for (unsigned i=0; i<conductionVelocityMaps.size(); i++)
02617     {
02618         cp::conduction_velocity_map_type temp(conductionVelocityMaps[i]);
02619         conduction_velocity_maps_sequence.push_back(temp);
02620     }
02621 }
02622 
02623 void HeartConfig::SetRequestedNodalTimeTraces (std::vector<unsigned>& requestedNodes)
02624 {
02625     EnsurePostProcessingSectionPresent();
02626     XSD_SEQUENCE_TYPE(cp::postprocessing_type::TimeTraceAtNode)& requested_nodes_sequence
02627         = mpParameters->PostProcessing()->TimeTraceAtNode();
02628 
02629     //Erase or create a sequence
02630     requested_nodes_sequence.clear();
02631 
02632     for (unsigned i=0; i<requestedNodes.size(); i++)
02633     {
02634         cp::node_number_type temp(requestedNodes[i]);
02635         requested_nodes_sequence.push_back(temp);
02636     }
02637 }
02638 
02639 template<unsigned SPACE_DIM>
02640 void HeartConfig::SetPseudoEcgElectrodePositions(const std::vector<ChastePoint<SPACE_DIM> >& rPseudoEcgElectrodePositions)
02641 {
02642     EnsurePostProcessingSectionPresent();
02643     XSD_SEQUENCE_TYPE(cp::postprocessing_type::PseudoEcgElectrodePosition)& electrodes_sequence
02644         = mpParameters->PostProcessing()->PseudoEcgElectrodePosition();
02645 
02646     //Erase or create a sequence
02647     electrodes_sequence.clear();
02648 
02649     for (unsigned i=0; i<rPseudoEcgElectrodePositions.size(); i++)
02650     {
02651         cp::point_type temp(rPseudoEcgElectrodePositions[i].GetWithDefault(0),
02652                             rPseudoEcgElectrodePositions[i].GetWithDefault(1),
02653                             rPseudoEcgElectrodePositions[i].GetWithDefault(2));
02654         electrodes_sequence.push_back(temp);
02655     }
02656 }
02657 
02658 
02659 /*
02660  * Output visualizer
02661  */
02662 
02663 void HeartConfig::EnsureOutputVisualizerExists()
02664 {
02665     ENSURE_SECTION_PRESENT(mpParameters->Simulation()->OutputVisualizer(), cp::output_visualizer_type);
02666 }
02667 
02668 void HeartConfig::SetVisualizeWithMeshalyzer(bool useMeshalyzer)
02669 {
02670     EnsureOutputVisualizerExists();
02671 
02672     mpParameters->Simulation()->OutputVisualizer()->meshalyzer(
02673         useMeshalyzer ? cp::yesno_type::yes : cp::yesno_type::no);
02674 }
02675 
02676 void HeartConfig::SetVisualizeWithCmgui(bool useCmgui)
02677 {
02678     EnsureOutputVisualizerExists();
02679 
02680     mpParameters->Simulation()->OutputVisualizer()->cmgui(
02681         useCmgui ? cp::yesno_type::yes : cp::yesno_type::no);
02682 }
02683 
02684 void HeartConfig::SetVisualizeWithVtk(bool useVtk)
02685 {
02686     EnsureOutputVisualizerExists();
02687 
02688     mpParameters->Simulation()->OutputVisualizer()->vtk(
02689         useVtk ? cp::yesno_type::yes : cp::yesno_type::no);
02690 }
02691 
02692 void HeartConfig::SetVisualizeWithParallelVtk(bool useParallelVtk)
02693 {
02694     EnsureOutputVisualizerExists();
02695 
02696     mpParameters->Simulation()->OutputVisualizer()->parallel_vtk(
02697         useParallelVtk ? cp::yesno_type::yes : cp::yesno_type::no);
02698 }
02699 
02700 void HeartConfig::SetVisualizerOutputPrecision(unsigned numberOfDigits)
02701 {
02702     EnsureOutputVisualizerExists();
02703 
02704     mpParameters->Simulation()->OutputVisualizer()->precision(numberOfDigits);
02705 }
02706 
02707 
02708 void HeartConfig::SetElectrodeParameters(bool groundSecondElectrode,
02709                                          unsigned index, double magnitude,
02710                                          double startTime, double duration )
02711 {
02712     assert(index < 3);
02713 
02714     cp::axis_type axis = cp::axis_type::x;
02715     if (index==1)
02716     {
02717         axis = cp::axis_type::y;
02718     }
02719     else if (index==2)
02720     {
02721         axis = cp::axis_type::z;
02722     }
02723 
02724     XSD_CREATE_WITH_FIXED_ATTR1(cp::surface_stimulus_strength_type, strength, magnitude, "uA/cm^2");
02725     XSD_CREATE_WITH_FIXED_ATTR1(cp::time_type, start_time, startTime, "ms");
02726     XSD_CREATE_WITH_FIXED_ATTR1(cp::time_type, duration_time, duration, "ms");
02727 
02728     if (!IsElectrodesPresent())
02729     {
02730         cp::electrodes_type element( groundSecondElectrode ? cp::yesno_type::yes : cp::yesno_type::no,
02731                                      axis,
02732                                      strength,
02733                                      start_time,
02734                                      duration_time );
02735         mpParameters->Simulation()->Electrodes().set(element);
02736     }
02737     else
02738     {
02739         mpParameters->Simulation()->Electrodes()->GroundSecondElectrode(groundSecondElectrode ? cp::yesno_type::yes : cp::yesno_type::no);
02740         mpParameters->Simulation()->Electrodes()->PerpendicularToAxis(axis);
02741         mpParameters->Simulation()->Electrodes()->Strength(strength);
02742         mpParameters->Simulation()->Electrodes()->StartTime(start_time);
02743         mpParameters->Simulation()->Electrodes()->Duration(duration_time);
02744     }
02745 }
02746 
02747 void HeartConfig::GetElectrodeParameters(bool& rGroundSecondElectrode,
02748                                          unsigned& rIndex, double& rMagnitude,
02749                                          double& rStartTime, double& rDuration )
02750 {
02751     if (!IsElectrodesPresent())
02752     {
02753         EXCEPTION("Attempted to get electrodes that have not been defined.");
02754     }
02755     else
02756     {
02757         rGroundSecondElectrode = (mpParameters->Simulation()->Electrodes()->GroundSecondElectrode()==cp::yesno_type::yes);
02758 
02759         cp::axis_type axis = mpParameters->Simulation()->Electrodes()->PerpendicularToAxis();
02760         if (axis==cp::axis_type::x)
02761         {
02762             rIndex = 0;
02763         }
02764         else if (axis==cp::axis_type::y)
02765         {
02766             rIndex = 1;
02767         }
02768         else
02769         {
02770             rIndex = 2;
02771         }
02772 
02773         rMagnitude = mpParameters->Simulation()->Electrodes()->Strength();
02774         rStartTime = mpParameters->Simulation()->Electrodes()->StartTime();
02775         rDuration = mpParameters->Simulation()->Electrodes()->Duration();
02776     }
02777 
02778 }
02779 
02780 bool HeartConfig::GetUseStateVariableInterpolation() const
02781 {
02782     // If it's an older version parameters & defaults (we're loading a checkpoint) say 'no'
02783     bool result = false;
02784     if (mpParameters->Numerical().UseStateVariableInterpolation().present())
02785     {
02786         result = mpParameters->Numerical().UseStateVariableInterpolation().get() == cp::yesno_type::yes;
02787     }
02788     return result;
02789 }
02790 
02791 void HeartConfig::SetUseStateVariableInterpolation(bool useStateVariableInterpolation)
02792 {
02793     if (useStateVariableInterpolation)
02794     {
02795         mpParameters->Numerical().UseStateVariableInterpolation().set(cp::yesno_type::yes);
02796     }
02797     else
02798     {
02799         mpParameters->Numerical().UseStateVariableInterpolation().set(cp::yesno_type::no);
02800     }
02801 }
02802 
02803 
02804 
02805 bool HeartConfig::HasDrugDose() const
02806 {
02807     return mpParameters->Physiological().ApplyDrug().present();
02808 }
02809 
02810 double HeartConfig::GetDrugDose() const
02811 {
02812     CHECK_EXISTS(HasDrugDose(), "Physiological/ApplyDrug");
02813     return mpParameters->Physiological().ApplyDrug()->concentration();
02814 }
02815 
02816 void HeartConfig::SetDrugDose(double drugDose)
02817 {
02818     if (!mpParameters->Physiological().ApplyDrug().present())
02819     {
02820         cp::apply_drug_type drug(drugDose);
02821         mpParameters->Physiological().ApplyDrug().set(drug);
02822     }
02823     else
02824     {
02825         mpParameters->Physiological().ApplyDrug()->concentration(drugDose);
02826     }
02827 }
02828 
02829 std::map<std::string, std::pair<double, double> > HeartConfig::GetIc50Values()
02830 {
02831     CHECK_EXISTS(HasDrugDose(), "Physiological/ApplyDrug");
02832     std::map<std::string, std::pair<double, double> > ic50s;
02833 
02834     XSD_SEQUENCE_TYPE(cp::apply_drug_type::IC50)&
02835         ic50_seq = mpParameters->Physiological().ApplyDrug()->IC50();
02836 
02837     for (XSD_ITERATOR_TYPE(cp::apply_drug_type::IC50) i = ic50_seq.begin();
02838          i != ic50_seq.end();
02839          ++i)
02840     {
02841         std::pair<double, double> ic50_hill(*i, i->hill());
02842         std::string current = i->current();
02843         ic50s[current] = ic50_hill;
02844     }
02845 
02846     return ic50s;
02847 }
02848 
02849 void HeartConfig::SetIc50Value(const std::string& rCurrentName, double ic50, double hill)
02850 {
02851     if (!mpParameters->Physiological().ApplyDrug().present())
02852     {
02853         SetDrugDose(0.0);
02854     }
02855     XSD_SEQUENCE_TYPE(cp::apply_drug_type::IC50)& ic50_seq = mpParameters->Physiological().ApplyDrug()->IC50();
02856     if (ic50_seq.empty())
02857     {
02858         // Erase or create a sequence
02859         ic50_seq.clear();
02860     }
02861     bool entry_exists = false;
02862     cp::ic50_type ic50_elt(ic50, rCurrentName);
02863     ic50_elt.hill(hill);
02864     for (XSD_ITERATOR_TYPE(cp::apply_drug_type::IC50) i = ic50_seq.begin();
02865          i != ic50_seq.end();
02866          ++i)
02867     {
02868         if (i->current() == rCurrentName)
02869         {
02870             entry_exists = true;
02871             *i = ic50_elt;
02872             break;
02873         }
02874     }
02875     if (!entry_exists)
02876     {
02877         ic50_seq.push_back(ic50_elt);
02878     }
02879 }
02880 
02881 
02882 
02883 void HeartConfig::SetUseMassLumping(bool useMassLumping)
02884 {
02885     mUseMassLumping = useMassLumping;
02886 }
02887 
02888 bool HeartConfig::GetUseMassLumping()
02889 {
02890     return mUseMassLumping;
02891 }
02892 
02893 void HeartConfig::SetUseMassLumpingForPrecond(bool useMassLumping)
02894 {
02895     mUseMassLumpingForPrecond = useMassLumping;
02896 }
02897 
02898 bool HeartConfig::GetUseMassLumpingForPrecond()
02899 {
02900     return mUseMassLumpingForPrecond;
02901 }
02902 
02903 void HeartConfig::SetUseReactionDiffusionOperatorSplitting(bool useOperatorSplitting)
02904 {
02905     mUseReactionDiffusionOperatorSplitting = useOperatorSplitting;
02906 }
02907 
02908 bool HeartConfig::GetUseReactionDiffusionOperatorSplitting()
02909 {
02910     return mUseReactionDiffusionOperatorSplitting;
02911 }
02912 
02913 void HeartConfig::SetUseFixedNumberIterationsLinearSolver(bool useFixedNumberIterations, unsigned evaluateNumItsEveryNSolves)
02914 {
02915     mUseFixedNumberIterations = useFixedNumberIterations;
02916     mEvaluateNumItsEveryNSolves = evaluateNumItsEveryNSolves;
02917 }
02918 
02919 bool HeartConfig::GetUseFixedNumberIterationsLinearSolver()
02920 {
02921     return mUseFixedNumberIterations;
02922 }
02923 
02924 unsigned HeartConfig::GetEvaluateNumItsEveryNSolves()
02925 {
02926     return mEvaluateNumItsEveryNSolves;
02927 }
02928 
02929 //
02930 // Purkinje methods
02931 //
02932 
02933 bool HeartConfig::HasPurkinje()
02934 {
02935     CheckSimulationIsDefined("Purkinje");
02936     return mpParameters->Simulation()->Purkinje().present();
02937 }
02938 
02939 double HeartConfig::GetPurkinjeCapacitance()
02940 {
02941     CHECK_EXISTS(mpParameters->Physiological().Purkinje().present(), "Physiological/Purkinje");
02942     CHECK_EXISTS(mpParameters->Physiological().Purkinje()->Capacitance().present(),
02943                  "Physiological/Purkinje/Capacitance");
02944     return mpParameters->Physiological().Purkinje()->Capacitance().get();
02945 }
02946 
02947 void HeartConfig::SetPurkinjeCapacitance(double capacitance)
02948 {
02949     ENSURE_SECTION_PRESENT(mpParameters->Physiological().Purkinje(), cp::purkinje_physiological_type);
02950     XSD_CREATE_WITH_FIXED_ATTR1(cp::capacitance_type, purk_Cm, capacitance, "uF/cm^2");
02951     mpParameters->Physiological().Purkinje()->Capacitance().set(purk_Cm);
02952 }
02953 
02954 
02955 double HeartConfig::GetPurkinjeSurfaceAreaToVolumeRatio()
02956 {
02957     CHECK_EXISTS(mpParameters->Physiological().Purkinje().present(), "Physiological/Purkinje");
02958     CHECK_EXISTS(mpParameters->Physiological().Purkinje()->SurfaceAreaToVolumeRatio().present(),
02959                  "Physiological/Purkinje/SurfaceAreaToVolumeRatio");
02960     return mpParameters->Physiological().Purkinje()->SurfaceAreaToVolumeRatio().get();
02961 }
02962 
02963 void HeartConfig::SetPurkinjeSurfaceAreaToVolumeRatio(double ratio)
02964 {
02965     ENSURE_SECTION_PRESENT(mpParameters->Physiological().Purkinje(), cp::purkinje_physiological_type);
02966     XSD_CREATE_WITH_FIXED_ATTR1(cp::inverse_length_type, purk_Am, ratio, "1/cm");
02967     mpParameters->Physiological().Purkinje()->SurfaceAreaToVolumeRatio().set(purk_Am);
02968 }
02969 
02970 double HeartConfig::GetPurkinjeConductivity()
02971 {
02972     CHECK_EXISTS(mpParameters->Physiological().Purkinje().present(), "Physiological/Purkinje");
02973     CHECK_EXISTS(mpParameters->Physiological().Purkinje()->Conductivity().present(),
02974             "Physiological/Purkinje/Conductivity");
02975     return mpParameters->Physiological().Purkinje()->Conductivity().get();
02976 }
02977 
02978 void HeartConfig::SetPurkinjeConductivity(double conductivity)
02979 {
02980     ENSURE_SECTION_PRESENT(mpParameters->Physiological().Purkinje(), cp::purkinje_physiological_type);
02981     XSD_CREATE_WITH_FIXED_ATTR1(cp::conductivity_type, purkinje_conductivity, conductivity, "mS/cm");
02982     mpParameters->Physiological().Purkinje()->Conductivity().set(purkinje_conductivity);
02983 }
02984 
02985 /**********************************************************************
02986  *                                                                    *
02987  *                                                                    *
02988  *            Utility methods for reading/transforming XML            *
02989  *                                                                    *
02990  *                                                                    *
02991  **********************************************************************/
02992 
02993 
02994 void XmlTransforms::TransformArchiveDirectory(xercesc::DOMDocument* pDocument,
02995                                               xercesc::DOMElement* pRootElement)
02996 {
02997     using namespace xercesc;
02998     std::vector<xercesc::DOMElement*> elts = XmlTools::FindElements(
02999         pRootElement,
03000         "ResumeSimulation/ArchiveDirectory");
03001     if (elts.size() > 0)
03002     {
03003         // We have an ArchiveDirectory element, so add the relative_to='chaste_test_output' attribute
03004         DOMElement* p_dir_elt = elts[0];
03005         p_dir_elt->setAttribute(X("relative_to"), X("chaste_test_output"));
03006     }
03007 }
03008 
03009 void XmlTransforms::TransformIonicModelDefinitions(xercesc::DOMDocument* pDocument,
03010                                                    xercesc::DOMElement* pRootElement)
03011 {
03012     // Default ionic model
03013     std::vector<xercesc::DOMElement*> p_elt_list = XmlTools::FindElements(
03014         pRootElement,
03015         "Simulation/IonicModels/Default");
03016     if (p_elt_list.size() > 0)
03017     {
03018         assert(p_elt_list.size() == 1); // Asserted by schema
03019         XmlTools::WrapContentInElement(pDocument, p_elt_list[0], X("Hardcoded"));
03020         // Now do any region-specific definitions
03021         p_elt_list = XmlTools::FindElements(pRootElement, "Simulation/IonicModels/Region/IonicModel");
03022         for (unsigned i=0; i<p_elt_list.size(); i++)
03023         {
03024             XmlTools::WrapContentInElement(pDocument, p_elt_list[i], X("Hardcoded"));
03025         }
03026     }
03027 }
03028 
03029 void XmlTransforms::CheckForIluPreconditioner(xercesc::DOMDocument* pDocument,
03030                                               xercesc::DOMElement* pRootElement)
03031 {
03032     std::vector<xercesc::DOMElement*> p_elt_list = XmlTools::FindElements(
03033         pRootElement,
03034         "Numerical/KSPPreconditioner");
03035     if (p_elt_list.size() > 0)
03036     {
03037         assert(p_elt_list.size() == 1); // Asserted by schema
03038         std::string text_value = X2C(p_elt_list[0]->getTextContent());
03039         if (text_value == "ilu")
03040         {
03041             EXCEPTION("PETSc does not have a parallel implementation of ilu, so we no longer allow it as an option.  Use bjacobi instead.");
03042         }
03043     }
03044 }
03045 
03046 void XmlTransforms::MoveConductivityHeterogeneities(xercesc::DOMDocument* pDocument,
03047                                                     xercesc::DOMElement* pRootElement)
03048 {
03049     std::vector<xercesc::DOMElement*> p_elt_list = XmlTools::FindElements(
03050             pRootElement,
03051             "Simulation/ConductivityHeterogeneities");
03052     if (p_elt_list.size() > 0)
03053     {
03054         assert(p_elt_list.size() == 1); // Asserted by schema
03055         xercesc::DOMNode* p_parent = p_elt_list[0]->getParentNode();
03056         xercesc::DOMNode* p_child = p_parent->removeChild(p_elt_list[0]);
03057         std::vector<xercesc::DOMElement*> p_phys_list = XmlTools::FindElements(pRootElement, "Physiological");
03058         assert(p_phys_list.size() == 1); // Asserted by schema
03059         p_phys_list[0]->appendChild(p_child);
03060     }
03061 }
03062 
03063 void XmlTransforms::SetDefaultVisualizer(xercesc::DOMDocument* pDocument,
03064                                          xercesc::DOMElement* pRootElement)
03065 {
03066     std::vector<xercesc::DOMElement*> p_sim_list = XmlTools::FindElements(pRootElement, "Simulation");
03067     if (p_sim_list.size() > 0)
03068     {
03069         std::vector<xercesc::DOMElement*> p_viz_list = XmlTools::FindElements(p_sim_list[0], "OutputVisualizer");
03070         if (p_viz_list.empty())
03071         {
03072             // Create the element and set meshalyzer (only) to on
03073             xercesc::DOMElement* p_viz_elt = pDocument->createElementNS(X("https://chaste.comlab.ox.ac.uk/nss/parameters/3_3"), X("OutputVisualizer"));
03074             p_sim_list[0]->appendChild(p_viz_elt);
03075             p_viz_elt->setAttribute(X("meshalyzer"), X("yes"));
03076         }
03077     }
03078 }
03079 
03081 // Explicit instantiation of the templated functions
03083 #define COVERAGE_IGNORE //These methods are covered above with DIM=1,2,3 but the instantiations may fail spuriously
03084 
03088 template void HeartConfig::GetIonicModelRegions<3u>(std::vector<boost::shared_ptr<AbstractChasteRegion<3u> > >& , std::vector<cp::ionic_model_selection_type>&) const;
03089 template void HeartConfig::GetStimuli<3u>(std::vector<boost::shared_ptr<AbstractStimulusFunction> >& , std::vector<boost::shared_ptr<AbstractChasteRegion<3u> > >& ) const;
03090 template void HeartConfig::GetCellHeterogeneities<3u>(std::vector<boost::shared_ptr<AbstractChasteRegion<3u> > >& ,std::vector<double>& ,std::vector<double>& ,std::vector<double>& ,std::vector<std::map<std::string, double> >*) ;
03091 template void HeartConfig::GetConductivityHeterogeneities<3u>(std::vector<boost::shared_ptr<AbstractChasteRegion<3u> > >& ,std::vector< c_vector<double,3> >& ,std::vector< c_vector<double,3> >& ) const;
03092 
03093 template void HeartConfig::GetIonicModelRegions<2u>(std::vector<boost::shared_ptr<AbstractChasteRegion<2u> > >& , std::vector<cp::ionic_model_selection_type>&) const;
03094 template void HeartConfig::GetStimuli<2u>(std::vector<boost::shared_ptr<AbstractStimulusFunction> >& , std::vector<boost::shared_ptr<AbstractChasteRegion<2u> > >& ) const;
03095 template void HeartConfig::GetCellHeterogeneities<2u>(std::vector<boost::shared_ptr<AbstractChasteRegion<2u> > >& ,std::vector<double>& ,std::vector<double>& ,std::vector<double>& ,std::vector<std::map<std::string, double> >*) ;
03096 template void HeartConfig::GetConductivityHeterogeneities<2u>(std::vector<boost::shared_ptr<AbstractChasteRegion<2u> > >& ,std::vector< c_vector<double,3> >& ,std::vector< c_vector<double,3> >& ) const;
03097 
03098 template void HeartConfig::GetIonicModelRegions<1u>(std::vector<boost::shared_ptr<AbstractChasteRegion<1u> > >& , std::vector<cp::ionic_model_selection_type>&) const;
03099 template void HeartConfig::GetStimuli<1u>(std::vector<boost::shared_ptr<AbstractStimulusFunction> >& , std::vector<boost::shared_ptr<AbstractChasteRegion<1u> > >& ) const;
03100 template void HeartConfig::GetCellHeterogeneities<1u>(std::vector<boost::shared_ptr<AbstractChasteRegion<1u> > >& ,std::vector<double>& ,std::vector<double>& ,std::vector<double>& ,std::vector<std::map<std::string, double> >*);
03101 template void HeartConfig::GetConductivityHeterogeneities<1u>(std::vector<boost::shared_ptr<AbstractChasteRegion<1u> > >& ,std::vector< c_vector<double,3> >& ,std::vector< c_vector<double,3> >& ) const;
03102 
03103 template void HeartConfig::GetPseudoEcgElectrodePositions(std::vector<ChastePoint<1u> >& rPseudoEcgElectrodePositions) const;
03104 template void HeartConfig::GetPseudoEcgElectrodePositions(std::vector<ChastePoint<2u> >& rPseudoEcgElectrodePositions) const;
03105 template void HeartConfig::GetPseudoEcgElectrodePositions(std::vector<ChastePoint<3u> >& rPseudoEcgElectrodePositions) const;
03106 
03107 template void HeartConfig::SetPseudoEcgElectrodePositions(const std::vector<ChastePoint<1u> >& rPseudoEcgElectrodePositions);
03108 template void HeartConfig::SetPseudoEcgElectrodePositions(const std::vector<ChastePoint<2u> >& rPseudoEcgElectrodePositions);
03109 template void HeartConfig::SetPseudoEcgElectrodePositions(const std::vector<ChastePoint<3u> >& rPseudoEcgElectrodePositions);
03114 #undef COVERAGE_IGNORE //These methods are covered above with DIM=1,2,3 but the instantiations may fail spuriously
03115 
03116 
03117 // Serialization for Boost >= 1.36
03118 #include "SerializationExportWrapperForCpp.hpp"
03119 CHASTE_CLASS_EXPORT(HeartConfig)

Generated by  doxygen 1.6.2