AbstractCorrectionTermAssembler.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 #include "AbstractCorrectionTermAssembler.hpp"
00030 #include <typeinfo>
00031 
00032 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00033 AbstractCorrectionTermAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AbstractCorrectionTermAssembler(
00034         AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00035         AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>* pTissue,
00036         unsigned numQuadPoints)
00037     : AbstractCardiacFeObjectAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,CARDIAC>(pMesh,pTissue,numQuadPoints)
00038 {
00039     // Work out which elements have the same cell at every node, and hence can have SVI done
00040     mElementsHasIdenticalCellModels.resize(pMesh->GetNumElements(), true);
00041     for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = pMesh->GetElementIteratorBegin();
00042          iter != pMesh->GetElementIteratorEnd();
00043          ++iter)
00044     {
00045         Element<ELEMENT_DIM, SPACE_DIM>& r_element = *iter;
00046         if (r_element.GetOwnership())
00047         {
00048             unsigned node_zero = r_element.GetNodeGlobalIndex(0);
00049             AbstractCardiacCell* p_cell_zero = this->mpCardiacTissue->GetCardiacCellOrHaloCell(node_zero);
00050             const std::type_info& r_zero_info = typeid(*p_cell_zero);
00051             // Check the other nodes match
00052             for (unsigned local_index=1; local_index<r_element.GetNumNodes(); local_index++)
00053             {
00054                 unsigned global_index = r_element.GetNodeGlobalIndex(local_index);
00055                 AbstractCardiacCell* p_cell = this->mpCardiacTissue->GetCardiacCellOrHaloCell(global_index);
00056                 const std::type_info& r_info = typeid(*p_cell);
00057                 if (r_zero_info != r_info)
00058                 {
00059                     mElementsHasIdenticalCellModels[r_element.GetIndex()] = false;
00060                     break;
00061                 }
00062             }
00063         }
00064     }
00065     // Note: the mStateVariables std::vector is resized if correction will be applied to a given element
00066 }
00067 
00068 
00069 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00070 void AbstractCorrectionTermAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ResetInterpolatedQuantities()
00071 {
00072     // reset ionic current, and state variables
00073     mIionicInterp = 0;
00074     for(unsigned i=0; i<mStateVariablesAtQuadPoint.size(); i++)
00075     {
00076         mStateVariablesAtQuadPoint[i] = 0;
00077     }
00078 }
00079 
00080 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00081 void AbstractCorrectionTermAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::IncrementInterpolatedQuantities(
00082             double phiI, const Node<SPACE_DIM>* pNode)
00083 {
00084     // interpolate ionic current, and state variables
00085     
00086     unsigned node_global_index = pNode->GetIndex();
00087 
00088     mIionicInterp  += phiI * this->mpCardiacTissue->rGetIionicCacheReplicated()[ node_global_index ];
00089     for (unsigned i=0; i<mStateVariablesAtQuadPoint.size(); i++)
00090     {
00091         mStateVariablesAtQuadPoint[i] += phiI * this->mpCardiacTissue->GetCardiacCellOrHaloCell(node_global_index)->rGetStateVariables()[i];
00092     }
00093 }
00094 
00095 
00096 
00097 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00098 bool AbstractCorrectionTermAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ElementAssemblyCriterion(Element<ELEMENT_DIM,SPACE_DIM>& rElement)
00099 {
00100     // if element doesn't have identical cell models, can't do SVI.
00101     if (!mElementsHasIdenticalCellModels[rElement.GetIndex()])
00102     {
00103         return false;
00104     }
00105     double DELTA_IIONIC = 1; // tolerance
00106 
00107     //The criterion and the correction both need the ionic cache, so we better make sure that it's up-to-date
00108     assert(this->mpCardiacTissue->GetDoCacheReplication());
00109     ReplicatableVector& r_cache = this->mpCardiacTissue->rGetIionicCacheReplicated();
00110     
00111     double diionic = fabs(r_cache[rElement.GetNodeGlobalIndex(0)] - r_cache[rElement.GetNodeGlobalIndex(1)]);
00112     
00113     if (ELEMENT_DIM > 1)
00114     {
00115         diionic = std::max(diionic, fabs(r_cache[rElement.GetNodeGlobalIndex(0)] - r_cache[rElement.GetNodeGlobalIndex(2)]) );    
00116         diionic = std::max(diionic, fabs(r_cache[rElement.GetNodeGlobalIndex(1)] - r_cache[rElement.GetNodeGlobalIndex(2)]) ); 
00117     }
00118 
00119     if (ELEMENT_DIM > 2)
00120     {
00121         diionic = std::max(diionic, fabs(r_cache[rElement.GetNodeGlobalIndex(0)] - r_cache[rElement.GetNodeGlobalIndex(3)]) ); 
00122         diionic = std::max(diionic, fabs(r_cache[rElement.GetNodeGlobalIndex(1)] - r_cache[rElement.GetNodeGlobalIndex(3)]) );    
00123         diionic = std::max(diionic, fabs(r_cache[rElement.GetNodeGlobalIndex(2)] - r_cache[rElement.GetNodeGlobalIndex(3)]) );
00124     }
00125     
00126     bool will_assemble = (diionic > DELTA_IIONIC);
00127     
00128     if (will_assemble)
00129     {
00130         unsigned any_node = rElement.GetNodeGlobalIndex(0);
00131         mStateVariablesAtQuadPoint.resize(this->mpCardiacTissue->GetCardiacCellOrHaloCell(any_node)->rGetStateVariables().size());
00132     }
00133     
00134     return will_assemble;
00135 }
00136 
00138 // explicit instantiation
00140 
00141 template class AbstractCorrectionTermAssembler<1,1,1>;
00142 template class AbstractCorrectionTermAssembler<1,2,1>;
00143 template class AbstractCorrectionTermAssembler<1,3,1>;
00144 template class AbstractCorrectionTermAssembler<2,2,1>;
00145 template class AbstractCorrectionTermAssembler<3,3,1>;
00146 template class AbstractCorrectionTermAssembler<1,1,2>;
00147 template class AbstractCorrectionTermAssembler<2,2,2>;
00148 template class AbstractCorrectionTermAssembler<3,3,2>;

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