ExtendedBidomainAssembler.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 "ExtendedBidomainAssembler.hpp"
00030 #include <boost/numeric/ublas/vector_proxy.hpp>
00031 
00032 #include "Exception.hpp"
00033 #include "DistributedVector.hpp"
00034 #include "PdeSimulationTime.hpp"
00035 #include "ConstBoundaryCondition.hpp"
00036 
00037 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00038 c_matrix<double,3*(ELEMENT_DIM+1),3*(ELEMENT_DIM+1)>
00039     ExtendedBidomainAssembler<ELEMENT_DIM,SPACE_DIM>::ComputeMatrixTerm(
00040             c_vector<double, ELEMENT_DIM+1> &rPhi,
00041             c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> &rGradPhi,
00042             ChastePoint<SPACE_DIM> &rX,
00043             c_vector<double,3> &rU,
00044             c_matrix<double, 3, SPACE_DIM> &rGradU /* not used */,
00045             Element<ELEMENT_DIM,SPACE_DIM>* pElement)
00046 {
00047     // get bidomain parameters
00048     double Am1 = mpExtendedBidomainTissue->GetAmFirstCell();
00049     double Am2 = mpExtendedBidomainTissue->GetAmSecondCell();
00050     double Cm1 = mpExtendedBidomainTissue->GetCmFirstCell();
00051     double Cm2 = mpExtendedBidomainTissue->GetCmSecondCell();
00052 
00053     const c_matrix<double, SPACE_DIM, SPACE_DIM>& sigma_i_first_cell = mpExtendedBidomainTissue->rGetIntracellularConductivityTensor(pElement->GetIndex());
00054     const c_matrix<double, SPACE_DIM, SPACE_DIM>& sigma_i_second_cell = mpExtendedBidomainTissue->rGetIntracellularConductivityTensorSecondCell(pElement->GetIndex());
00055     const c_matrix<double, SPACE_DIM, SPACE_DIM>& sigma_e = mpExtendedBidomainTissue->rGetExtracellularConductivityTensor(pElement->GetIndex());
00056 
00057     double delta_t = PdeSimulationTime::GetPdeTimeStep();
00058 
00059     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> temp_1 = prod(sigma_i_first_cell, rGradPhi);
00060     c_matrix<double, ELEMENT_DIM+1, ELEMENT_DIM+1> grad_phi_sigma_i_first_cell_grad_phi =
00061         prod(trans(rGradPhi), temp_1);
00062 
00063     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> temp_2 = prod(sigma_i_second_cell, rGradPhi);
00064     c_matrix<double, ELEMENT_DIM+1, ELEMENT_DIM+1> grad_phi_sigma_i_second_cell_grad_phi =
00065         prod(trans(rGradPhi), temp_2);
00066 
00067     c_matrix<double, ELEMENT_DIM+1, ELEMENT_DIM+1> basis_outer_prod =
00068         outer_prod(rPhi, rPhi);
00069 
00070     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> temp_ext = prod(sigma_e, rGradPhi);
00071     c_matrix<double, ELEMENT_DIM+1, ELEMENT_DIM+1> grad_phi_sigma_e_grad_phi =
00072         prod(trans(rGradPhi), temp_ext);
00073 
00074 
00075     c_matrix<double,3*(ELEMENT_DIM+1),3*(ELEMENT_DIM+1)> ret;
00076 
00077     // first equation, first unknown
00078     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00079     slice100(ret, slice (0, 3, ELEMENT_DIM+1), slice (0, 3, ELEMENT_DIM+1));
00080     slice100 = (Am1*Cm1/delta_t)*basis_outer_prod + grad_phi_sigma_i_first_cell_grad_phi;
00081 
00082     // first equation, second unknown
00083     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00084     slice200(ret, slice (0, 3, ELEMENT_DIM+1), slice (1, 3, ELEMENT_DIM+1));
00085     slice200 = zero_matrix<double>(ELEMENT_DIM+1, ELEMENT_DIM+1);
00086 
00087     // first equation, third unknown
00088     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00089     slice300(ret, slice (0, 3, ELEMENT_DIM+1), slice (2, 3, ELEMENT_DIM+1));
00090     slice300 = - (Am1*Cm1/delta_t)*basis_outer_prod;
00091 
00092     // second equation, first unknown
00093     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00094     slice010(ret, slice (1, 3, ELEMENT_DIM+1), slice (0, 3, ELEMENT_DIM+1));
00095     slice010 = zero_matrix<double>(ELEMENT_DIM+1, ELEMENT_DIM+1);
00096 
00097     // second equation, second unknown
00098     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00099     slice020(ret, slice (1, 3, ELEMENT_DIM+1), slice (1, 3, ELEMENT_DIM+1));
00100     slice020 = (Am2*Cm2/delta_t)*basis_outer_prod + grad_phi_sigma_i_second_cell_grad_phi;
00101 
00102     // second equation, third unknown
00103     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00104     slice030(ret, slice (1, 3, ELEMENT_DIM+1), slice (2, 3, ELEMENT_DIM+1));
00105     slice030 = - (Am2*Cm2/delta_t)*basis_outer_prod;
00106 
00107     // third equation, first unknown
00108     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00109     slice001(ret, slice (2, 3, ELEMENT_DIM+1), slice (0, 3, ELEMENT_DIM+1));
00110     slice001 =   - grad_phi_sigma_i_first_cell_grad_phi;
00111 
00112     // third equation, second unknown
00113     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00114     slice002(ret, slice (2, 3, ELEMENT_DIM+1), slice (1, 3, ELEMENT_DIM+1));
00115     slice002 =  - grad_phi_sigma_i_second_cell_grad_phi;
00116 
00117     // third equation, third unknown
00118     matrix_slice<c_matrix<double, 3*ELEMENT_DIM+3, 3*ELEMENT_DIM+3> >
00119     slice003(ret, slice (2, 3, ELEMENT_DIM+1), slice (2, 3, ELEMENT_DIM+1));
00120     slice003 =  - grad_phi_sigma_e_grad_phi;
00121 
00122     return ret;
00123 }
00124 
00125 
00126 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00127 ExtendedBidomainAssembler<ELEMENT_DIM,SPACE_DIM>::ExtendedBidomainAssembler(
00128                                 AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00129                                 ExtendedBidomainTissue<SPACE_DIM>* pTissue,
00130                                 unsigned numQuadPoints)
00131     : AbstractCardiacFeVolumeIntegralAssembler<ELEMENT_DIM,SPACE_DIM,3,true,true,NORMAL>(pMesh,pTissue,numQuadPoints),
00132               mpExtendedBidomainTissue(pTissue)
00133 {
00134     assert(pTissue != NULL);
00135     mpConfig = HeartConfig::Instance();
00136 }
00137 
00138 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00139 ExtendedBidomainAssembler<ELEMENT_DIM,SPACE_DIM>::~ExtendedBidomainAssembler()
00140 {
00141 }
00142 
00144 // Explicit instantiation
00146 
00147 template class ExtendedBidomainAssembler<1,1>;
00148 template class ExtendedBidomainAssembler<2,2>;
00149 template class ExtendedBidomainAssembler<3,3>;