SimpleBathProblemSetup.hpp

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #ifndef SIMPLEBATHPROBLEMSETUP_HPP_
00029 #define SIMPLEBATHPROBLEMSETUP_HPP_
00030 
00037 #include "AbstractCardiacCellFactory.hpp"
00038 #include "SimpleStimulus.hpp"
00039 #include "LuoRudy1991.hpp"
00040 #include "HeartRegionCodes.hpp"
00041 
00046 template<unsigned DIM, class CELLTYPE=CellLuoRudy1991FromCellML>
00047 class BathCellFactory : public AbstractCardiacCellFactory<DIM>
00048 {
00049 private:
00051     boost::shared_ptr<SimpleStimulus> mpStimulus;
00053     c_vector<double,DIM> mStimulatedPoint;
00054 
00055 public:
00062     BathCellFactory(double stimulusMagnitude, c_vector<double,DIM> stimulatedPoint)
00063         : AbstractCardiacCellFactory<DIM>(),
00064           mpStimulus(new SimpleStimulus(stimulusMagnitude, 0.5)),
00065           mStimulatedPoint(stimulatedPoint)
00066     {
00067     }
00068 
00073     AbstractCardiacCell* CreateCardiacCellForTissueNode(unsigned node)
00074     {
00075         // paranoia - check this is really a tissue node
00076         assert(HeartRegionCode::IsRegionTissue( this->GetMesh()->GetNode(node)->GetRegion() ));
00077 
00078         // stimulate centre node normally..
00079         bool is_centre;
00080 
00081         if (DIM==1)
00082         {
00083             is_centre = (fabs(this->GetMesh()->GetNode(node)->GetPoint()[0]-mStimulatedPoint(0)) < 1e-6);
00084         }
00085         else if (DIM==2)
00086         {
00087             is_centre = (    (fabs(this->GetMesh()->GetNode(node)->GetPoint()[0]-mStimulatedPoint(0)) < 1e-6)
00088                           && (fabs(this->GetMesh()->GetNode(node)->GetPoint()[1]-mStimulatedPoint(1)) < 1e-6) );
00089         }
00090         else
00091         {
00092             is_centre = (    (fabs(this->GetMesh()->GetNode(node)->GetPoint()[0]-mStimulatedPoint(0)) < 1e-6)
00093                           && (fabs(this->GetMesh()->GetNode(node)->GetPoint()[1]-mStimulatedPoint(1)) < 1e-6)
00094                           && (fabs(this->GetMesh()->GetNode(node)->GetPoint()[2]-mStimulatedPoint(2)) < 1e-6) );
00095         }
00096 
00097         if (is_centre)
00098         {
00099             return new CELLTYPE(this->mpSolver, mpStimulus);
00100         }
00101         else
00102         {
00103             return new CELLTYPE(this->mpSolver, this->mpZeroStimulus);
00104         }
00105     }
00106 };
00107 
00116 template<class MeshType>
00117 void SetCircularTissueIn2dMesh(MeshType* pMesh,
00118                                double centreX, double centreY, double radius)
00119 {
00120     for (typename MeshType::ElementIterator it = pMesh->GetElementIteratorBegin();
00121          it != pMesh->GetElementIteratorEnd();
00122          ++it)
00123     {
00124         double x = it->CalculateCentroid()[0];
00125         double y = it->CalculateCentroid()[1];
00126         if ( (x-centreX)*(x-centreX) + (y-centreY)*(y-centreY) > radius*radius )
00127         {
00128             it->SetRegion(HeartRegionCode::GetValidBathId());
00129         }
00130     }
00131     pMesh->SetMeshHasChangedSinceLoading();
00132 }
00133 
00142 template<class MeshType>
00143 MeshType* Load2dMeshAndSetCircularTissue(const std::string& rMeshPath,
00144                                          double centreX, double centreY, double radius)
00145 {
00146     TrianglesMeshReader<2,2> reader(rMeshPath);
00147     MeshType* p_mesh = new MeshType;
00148     p_mesh->ConstructFromMeshReader(reader);
00149 
00150     SetCircularTissueIn2dMesh(p_mesh, centreX, centreY, radius);
00151 
00152     return p_mesh;
00153 }
00154 
00163 template<>
00164 DistributedTetrahedralMesh<2,2>* Load2dMeshAndSetCircularTissue(const std::string& rMeshPath,
00165                                                                 double centreX, double centreY, double radius)
00166 {
00167     TrianglesMeshReader<2,2> reader(rMeshPath);
00168     // Force dumb partitioning so migration tests pass!
00169     DistributedTetrahedralMesh<2,2>* p_mesh = new DistributedTetrahedralMesh<2,2>(DistributedTetrahedralMeshPartitionType::DUMB);
00170     p_mesh->ConstructFromMeshReader(reader);
00171 
00172     SetCircularTissueIn2dMesh(p_mesh, centreX, centreY, radius);
00173 
00174     return p_mesh;
00175 }
00176 
00177 #endif /*SIMPLEBATHPROBLEMSETUP_HPP_*/

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