ExtendedBidomainAssembler.cpp

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