Chaste Release::3.1
|
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)