00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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 : AbstractCardiacFeVolumeIntegralAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,CARDIAC>(pMesh,pTissue,numQuadPoints)
00038 {
00039
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
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
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
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
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
00101 if (!mElementsHasIdenticalCellModels[rElement.GetIndex()])
00102 {
00103 return false;
00104 }
00105 double DELTA_IIONIC = 1;
00106
00107
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
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>;