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