00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
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
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
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
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
00194
00195 #include "HeartConfigDefaults.hpp"
00196
00197
00198
00199
00200 std::auto_ptr<HeartConfig> HeartConfig::mpInstance;
00201
00202
00203
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
00227
00228
00229 mEpiFraction = -1.0;
00230 mEndoFraction = -1.0;
00231 mMidFraction = -1.0;
00232 mUserAskedForCellularTransmuralHeterogeneities = false;
00233
00234
00235
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
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
00267
00268
00269
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
00283
00284 ::xml_schema::namespace_infomap map;
00285
00286 map[""].schema = "ChasteParameters_1_1.xsd";
00287
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
00301
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
00308 if (useArchiveLocationInfo)
00309 {
00310 CopySchema(output_dirname);
00311 }
00312
00313 }
00314
00315 void HeartConfig::LoadFromCheckpoint()
00316 {
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326 assert(mpParameters.use_count() > 0);
00327 boost::shared_ptr<cp::chaste_parameters_type> p_new_parameters = mpParameters;
00328
00329
00330
00331
00332
00333
00334
00335
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
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
00352 if (p_new_parameters->ResumeSimulation().present())
00353 {
00354 UpdateParametersFromResumeSimulation(p_new_parameters);
00355 }
00356
00357 CheckTimeSteps();
00358 }
00359
00360 void HeartConfig::CopySchema(const std::string& rToDirectory)
00361 {
00362
00363
00364
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
00371 schema_location.SetPath(schema_name, RelativeTo::CWD);
00372 if (!schema_location.Exists())
00373 {
00374
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
00391 std::string root_dir = std::string(ChasteBuildInfo::GetRootDir()) + "/heart/src/io/";
00392
00393 mSchemaLocations[""] = root_dir + "ChasteParameters_1_1.xsd";
00394
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
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
00471
00472 try
00473 {
00474
00475 XmlTools::Finalizer finalizer(false);
00476
00477 xsd::cxx::xml::dom::auto_ptr<xercesc::DOMDocument> p_doc = XmlTools::ReadXmlFile(rFileName, props);
00478
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)
00483 {
00484 XmlTransforms::TransformIonicModelDefinitions(p_doc.get(), p_root_elt);
00485 }
00486 if (version < 2001)
00487 {
00488 XmlTransforms::TransformArchiveDirectory(p_doc.get(), p_root_elt);
00489 XmlTransforms::CheckForIluPreconditioner(p_doc.get(), p_root_elt);
00490 }
00491 if (version < 3001)
00492 {
00493 XmlTransforms::MoveConductivityHeterogeneities(p_doc.get(), p_root_elt);
00494 }
00495 if (version < 3003)
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
00501 std::auto_ptr<cp::chaste_parameters_type> p_params(cp::ChasteParameters(*p_doc, ::xml_schema::flags::dont_initialize, props));
00502
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
00511 mpParameters.reset();
00512 EXCEPTION("XML parsing error in configuration file: " + rFileName);
00513 }
00514 catch (...)
00515 {
00516
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();
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
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
00549 HeartConfig::Instance()->SetSimulationDuration(pResumeParameters->ResumeSimulation()->SimulationDuration());
00550
00551
00552 if (pResumeParameters->ResumeSimulation()->Stimuli().present())
00553 {
00554 mpParameters->Simulation()->Stimuli().set(pResumeParameters->ResumeSimulation()->Stimuli().get());
00555 }
00556
00557
00558
00559 if (pResumeParameters->ResumeSimulation()->CellHeterogeneities().present())
00560 {
00561 if (!mpParameters->Simulation()->CellHeterogeneities().present())
00562 {
00563
00564 mpParameters->Simulation()->CellHeterogeneities().set(pResumeParameters->ResumeSimulation()->CellHeterogeneities().get());
00565 }
00566 else
00567 {
00568
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
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
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
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
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
00654 mpInstance.reset();
00655
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
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;
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
00980
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
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
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
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
01219 mUserAskedForCellularTransmuralHeterogeneities = user_asking_for_transmural_layers;
01220
01221 if ( mUserAskedForCellularTransmuralHeterogeneities )
01222 {
01223
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
01230
01231
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
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
01525
01526
01527
01528
01529
01530
01531
01532 if (bathRegion == UINT_MAX)
01533 {
01534
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
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
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
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
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
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
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
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
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
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
02313
02314
02315
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
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
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
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
02442 WARNING("Code has changed the KSP solver type from "<<GetKSPSolver()<<" to "<< kspSolver);
02443 }
02444
02445
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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);
03019 XmlTools::WrapContentInElement(pDocument, p_elt_list[0], X("Hardcoded"));
03020
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);
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);
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);
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
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
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
03118 #include "SerializationExportWrapperForCpp.hpp"
03119 CHASTE_CLASS_EXPORT(HeartConfig)