BidomainProblem.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #include "BidomainProblem.hpp"
00031 #include "BasicBidomainSolver.hpp"
00032 #include "MatrixBasedBidomainSolver.hpp"
00033 #include "HeartConfig.hpp"
00034 #include "Exception.hpp"
00035 #include "DistributedVector.hpp"
00036 #include "ReplicatableVector.hpp"
00037 
00038 template <unsigned DIM>
00039 void BidomainProblem<DIM>::AnalyseMeshForBath()
00040 {
00041     // Annotate bath notes with correct region code
00042     if (mHasBath)
00043     {
00044         // Initialize all nodes to be bath nodes
00045         for (typename AbstractTetrahedralMesh<DIM,DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
00046              iter != this->mpMesh->GetNodeIteratorEnd();
00047             ++iter)
00048         {
00049             (*iter).SetRegion(HeartRegionCode::GetValidBathId());
00050         }
00051 
00052         bool any_bath_element_found = false;
00053 
00054         // Set nodes that are part of a heart element to be heart nodes
00055         //for (unsigned i=0; i<this->mpMesh->GetNumElements(); i++)
00056         for (typename AbstractTetrahedralMesh<DIM,DIM>::ElementIterator it = this->mpMesh->GetElementIteratorBegin();
00057              it != this->mpMesh->GetElementIteratorEnd();
00058              ++it)
00059         {
00060             Element<DIM, DIM>& r_element = *it;
00061 
00062             if (HeartRegionCode::IsRegionTissue( r_element.GetRegion() ))
00063             {
00064                 for (unsigned j=0; j<r_element.GetNumNodes(); j++)
00065                 {
00066                     r_element.GetNode(j)->SetRegion(HeartRegionCode::GetValidTissueId());
00067                 }
00068             }
00069             else
00070             {
00071                 assert(HeartRegionCode::IsRegionBath( r_element.GetRegion() ));
00072                 any_bath_element_found = true;
00073             }
00074         }
00075 
00076         if (!PetscTools::ReplicateBool(any_bath_element_found))
00077         {
00078            EXCEPTION("No bath element found");
00079         }
00080     }
00081 }
00082 
00083 
00084 template<unsigned DIM>
00085 Vec BidomainProblem<DIM>::CreateInitialCondition()
00086 {
00087     Vec init_cond = AbstractCardiacProblem<DIM,DIM,2>::CreateInitialCondition();
00088     if (mHasBath)
00089     {
00090         // get the voltage stripe
00091         DistributedVector ic = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(init_cond);
00092         DistributedVector::Stripe voltage_stripe = DistributedVector::Stripe(ic,0);
00093 
00094         for (DistributedVector::Iterator index = ic.Begin();
00095              index!= ic.End();
00096              ++index)
00097         {
00098             if (HeartRegionCode::IsRegionBath( this->mpMesh->GetNode( index.Global )->GetRegion() ))
00099             {
00100                 voltage_stripe[index] = 0.0;
00101             }
00102         }
00103         ic.Restore();
00104     }
00105 
00106     return init_cond;
00107 }
00108 
00109 template<unsigned DIM>
00110 AbstractCardiacTissue<DIM> * BidomainProblem<DIM>::CreateCardiacTissue()
00111 {
00112     AnalyseMeshForBath();
00113     mpBidomainTissue = new BidomainTissue<DIM>(this->mpCellFactory, HeartConfig::Instance()->GetUseStateVariableInterpolation());
00114     return mpBidomainTissue;
00115 }
00116 
00117 template<unsigned DIM>
00118 AbstractDynamicLinearPdeSolver<DIM, DIM, 2>* BidomainProblem<DIM>::CreateSolver()
00119 {
00120     /*
00121      * NOTE: The this->mpBoundaryConditionsContainer.get() lines below convert a
00122      * boost::shared_ptr to a normal pointer, as this is what the solvers are
00123      * expecting. We have to be a bit careful though as boost could decide to delete
00124      * them whenever it feels like as it won't count the assembers as using them.
00125      *
00126      * As long as they are kept as member variables here for as long as they are
00127      * required in the solvers it should all work OK.
00128      */
00129     if (!this->mUseMatrixBasedRhsAssembly)
00130     {
00131         mpSolver = new BasicBidomainSolver<DIM,DIM>(mHasBath, 
00132                                                     this->mpMesh,
00133                                                     mpBidomainTissue,
00134                                                     this->mpBoundaryConditionsContainer.get(),
00135                                                     2);
00136     }
00137     else
00138     {
00139         mpSolver = new MatrixBasedBidomainSolver<DIM,DIM>(mHasBath, 
00140                                                           this->mpMesh,
00141                                                           mpBidomainTissue,
00142                                                           this->mpBoundaryConditionsContainer.get(),
00143                                                           2);
00144     }
00145 
00146     try
00147     {
00148         mpSolver->SetFixedExtracellularPotentialNodes(mFixedExtracellularPotentialNodes);
00149         mpSolver->SetRowForAverageOfPhiZeroed(mRowForAverageOfPhiZeroed);
00150     }
00151     catch (const Exception& e)
00152     {
00153         delete mpSolver;
00154         throw e;
00155     }
00156 
00157     return mpSolver;
00158 }
00159 
00160 template<unsigned DIM>
00161 BidomainProblem<DIM>::BidomainProblem(
00162             AbstractCardiacCellFactory<DIM>* pCellFactory, bool hasBath)
00163     : AbstractCardiacProblem<DIM,DIM, 2>(pCellFactory),
00164       mpBidomainTissue(NULL),
00165       mRowForAverageOfPhiZeroed(INT_MAX),
00166       mHasBath(hasBath)
00167 {
00168     mFixedExtracellularPotentialNodes.resize(0);
00169 }
00170 
00171 template<unsigned DIM>
00172 BidomainProblem<DIM>::BidomainProblem()
00173     : AbstractCardiacProblem<DIM, DIM, 2>(),
00174       mpBidomainTissue(NULL),
00175       mRowForAverageOfPhiZeroed(INT_MAX)
00176 {
00177     mFixedExtracellularPotentialNodes.resize(0);
00178 }
00179 
00180 
00181 
00182 template<unsigned DIM>
00183 void BidomainProblem<DIM>::SetFixedExtracellularPotentialNodes(std::vector<unsigned> nodes)
00184 {
00185     mFixedExtracellularPotentialNodes.resize(nodes.size());
00186     for (unsigned i=0; i<nodes.size(); i++)
00187     {
00188         // the solver checks that the nodes[i] is less than
00189         // the number of nodes in the mesh so this is not done here
00190         mFixedExtracellularPotentialNodes[i] = nodes[i];
00191     }
00192 }
00193 
00194 template<unsigned DIM>
00195 void BidomainProblem<DIM>::SetNodeForAverageOfPhiZeroed(unsigned node)
00196 {
00197     mRowForAverageOfPhiZeroed = 2*node+1;
00198 }
00199 
00200 template<unsigned DIM>
00201 BidomainTissue<DIM>* BidomainProblem<DIM>::GetBidomainTissue()
00202 {
00203     assert(mpBidomainTissue!=NULL);
00204     return mpBidomainTissue;
00205 }
00206 
00207 template<unsigned DIM>
00208 void BidomainProblem<DIM>::WriteInfo(double time)
00209 {
00210     if (PetscTools::AmMaster())
00211     {
00212         std::cout << "Solved to time " << time << "\n" << std::flush;
00213     }
00214 
00215     double v_max, v_min, phi_max, phi_min;
00216 
00217     VecSetBlockSize( this->mSolution, 2 );
00218 
00219     VecStrideMax( this->mSolution, 0, PETSC_NULL, &v_max );
00220     VecStrideMin( this->mSolution, 0, PETSC_NULL, &v_min );
00221 
00222     VecStrideMax( this->mSolution, 1, PETSC_NULL, &phi_max );
00223     VecStrideMin( this->mSolution, 1, PETSC_NULL, &phi_min );
00224 
00225     if (PetscTools::AmMaster())
00226     {
00227         std::cout << " V; phi_e = " << "[" <<v_min << ", " << v_max << "]" << ";\t"
00228                   << "[" <<phi_min << ", " << phi_max << "]" << "\n"
00229                   << std::flush;
00230     }
00231 }
00232 
00233 template<unsigned DIM>
00234 void BidomainProblem<DIM>::DefineWriterColumns(bool extending)
00235 {
00236     AbstractCardiacProblem<DIM,DIM,2>::DefineWriterColumns(extending);
00237     if (extending)
00238     {
00239         mExtracelluarColumnId = this->mpWriter->GetVariableByName("Phi_e");
00240     }
00241     else
00242     {
00243         mExtracelluarColumnId = this->mpWriter->DefineVariable("Phi_e","mV");
00244     }
00245     AbstractCardiacProblem<DIM,DIM,2>::DefineExtraVariablesWriterColumns(extending);
00246 }
00247 
00248 template<unsigned DIM>
00249 void BidomainProblem<DIM>::WriteOneStep(double time, Vec voltageVec)
00250 {
00251     this->mpWriter->PutUnlimitedVariable(time);
00252     std::vector<int> variable_ids;
00253     variable_ids.push_back(this->mVoltageColumnId);
00254     variable_ids.push_back(mExtracelluarColumnId);
00255     this->mpWriter->PutStripedVector(variable_ids, voltageVec);
00256     AbstractCardiacProblem<DIM,DIM,2>::WriteExtraVariablesOneStep();
00257 }
00258 
00259 template<unsigned DIM>
00260 void BidomainProblem<DIM>::PreSolveChecks()
00261 {
00262     AbstractCardiacProblem<DIM,DIM, 2>::PreSolveChecks();
00263     if (mFixedExtracellularPotentialNodes.empty())
00264     {
00265         // We're not pinning any nodes.
00266         if (mRowForAverageOfPhiZeroed==INT_MAX)
00267         {
00268             // We're not using the constrain Average phi_e to 0 method, hence use a null space
00269             // Check that the KSP solver isn't going to do anything stupid:
00270             // phi_e is not bounded, so it'd be wrong to use a relative tolerance
00271             if (HeartConfig::Instance()->GetUseRelativeTolerance())
00272             {
00273                 EXCEPTION("Bidomain external voltage is not bounded in this simulation - use KSP *absolute* tolerance");
00274             }
00275         }
00276     }
00277 }
00278 
00279 template<unsigned DIM>
00280 void BidomainProblem<DIM>::SetElectrodes()
00281 {
00282     if (!mHasBath)
00283     {
00284         //Cannot set electrodes when problem has been defined to not have a bath
00285         return;
00286     }
00287     
00288     assert(this->mpMesh!=NULL);
00289     
00290     if (HeartConfig::Instance()->IsElectrodesPresent())
00291     { 
00292         mpElectrodes = (boost::shared_ptr<Electrodes<DIM> >) new Electrodes<DIM>(*(this->mpMesh));
00293     }
00294 }
00295 
00296 template<unsigned DIM>
00297 void BidomainProblem<DIM>::AtBeginningOfTimestep(double time)
00298 {
00299     if ( mpElectrodes && mpElectrodes->SwitchOn(time) )
00300     {
00301         // At the moment mpBcc and mpDefaultBcc point to a set default BC
00302         assert(this->mpBoundaryConditionsContainer);
00303         //assert(this->mpDefaultBoundaryConditionsContainer);
00304 
00305         // Note, no point calling this->SetBoundaryConditionsContainer() as the
00306         // solver has already been created..
00307         mpSolver->ResetBoundaryConditionsContainer(mpElectrodes->GetBoundaryConditionsContainer().get());
00308 
00309         // ..but we set mpBcc anyway, so the local mpBcc is
00310         // the same as the one being used in the solver...
00311         this->mpBoundaryConditionsContainer = mpElectrodes->GetBoundaryConditionsContainer();
00312 
00314         this->mpDefaultBoundaryConditionsContainer = this->mpBoundaryConditionsContainer;
00315 
00316         // At t==0 or after checkpointing we won't have a system assembled at this stage: BCs will be applied once the matrix
00317         // is assembled. Dirichlet BCs will be present at the time of assembly and no null space will be created either.
00318         if ( mpSolver->GetLinearSystem() != NULL )
00319         {
00320             // System matrix is assembled once at the beginning of the simulation. After that, nobody will take care
00321             // of applying new BC to the system matrix. Must be triggered explicitly.
00322             if (mpElectrodes->HasGroundedElectrode())
00323             {
00324                 this->mpBoundaryConditionsContainer->ApplyDirichletToLinearProblem( *(mpSolver->GetLinearSystem()),
00325                                                                                    true, false);
00326             }
00327 
00328             // If a grounded electrode is switched on, the linear system is not singular anymore. Remove the null space.
00329             if (mpElectrodes->HasGroundedElectrode())
00330             {
00331                 mpSolver->GetLinearSystem()->RemoveNullSpace();
00332             }
00333         }
00334     }
00335 }
00336 
00337 template<unsigned DIM>
00338 void BidomainProblem<DIM>::OnEndOfTimestep(double time)
00339 {
00340     if ( mpElectrodes && mpElectrodes->SwitchOff(time) )
00341     {
00342         // At the moment mpBcc should exist and therefore
00343         // mpDefaultBcc should be empty (not if electrodes switched on after 0ms)
00344         assert(this->mpBoundaryConditionsContainer);
00345         //assert(! this->mpDefaultBoundaryConditionsContainer);
00346 
00347         // Set up default boundary conditions container - no Neumann fluxes
00348         // or Dirichlet fixed nodes
00349         this->mpDefaultBoundaryConditionsContainer.reset(new BoundaryConditionsContainer<DIM,DIM,2>);
00350         for (unsigned problem_index=0; problem_index<2; problem_index++)
00351         {
00352             this->mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(this->mpMesh, problem_index);
00353         }
00354 
00355         // If there's a grounded electrode, we must remove BC from linear system. At the moment, we don't
00356         // have a sensible way of doing this, therefore we reassemble the system.
00357         if (mpElectrodes->HasGroundedElectrode())
00358         {
00359             delete mpSolver;
00360             AbstractCardiacProblem<DIM,DIM,2>::mpSolver = CreateSolver();
00361             mpSolver->SetTimeStep(HeartConfig::Instance()->GetPdeTimeStep());
00362         }
00363 
00364         // Note, no point calling this->SetBoundaryConditionsContainer() as the
00365         // solver has already been created..
00366         mpSolver->ResetBoundaryConditionsContainer(this->mpDefaultBoundaryConditionsContainer.get());
00367         // ..but we set mpBcc to be mpDefaultBcc anyway, so the local mpBcc is
00368         // the same as the one being used in the solver...
00369         this->mpBoundaryConditionsContainer = this->mpDefaultBoundaryConditionsContainer;
00370     }
00371 }
00372 
00373 
00374 
00375 template<unsigned DIM>
00376 void BidomainProblem<DIM>::SetUpAdditionalStoppingTimes(std::vector<double>& rAdditionalStoppingTimes)
00377 {
00378     if ( mpElectrodes )
00379     {
00380         rAdditionalStoppingTimes.push_back( mpElectrodes->GetSwitchOnTime() );
00381         rAdditionalStoppingTimes.push_back( mpElectrodes->GetSwitchOffTime() );
00382     }
00383 }
00384 
00385 template<unsigned DIM>
00386 bool BidomainProblem<DIM>::GetHasBath()
00387 {
00388     return mHasBath;
00389 }
00390 
00392 // Explicit instantiation
00394 
00395 template class BidomainProblem<1>;
00396 template class BidomainProblem<2>;
00397 template class BidomainProblem<3>;
00398 
00399 
00400 // Serialization for Boost >= 1.36
00401 #include "SerializationExportWrapperForCpp.hpp"
00402 EXPORT_TEMPLATE_CLASS_SAME_DIMS(BidomainProblem)

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