BidomainWithBathMatrixBasedAssembler.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 "BidomainWithBathMatrixBasedAssembler.hpp"
00030 
00031 #include <boost/numeric/ublas/vector_proxy.hpp>
00032 
00033 #include "ConstBoundaryCondition.hpp"
00034 #include "HeartRegionCodes.hpp"
00035 
00036 
00038 // BidomainWithBathRhsMatrixAssembler
00040 
00041 template<unsigned DIM>
00042 c_matrix<double,2*(DIM+1),2*(DIM+1)> BidomainWithBathRhsMatrixAssembler<DIM>::ComputeMatrixTerm(
00043             c_vector<double, DIM+1> &rPhi,
00044             c_matrix<double, DIM, DIM+1> &rGradPhi,
00045             ChastePoint<DIM> &rX,
00046             c_vector<double,2> &u,
00047             c_matrix<double,2,DIM> &rGradU /* not used */,
00048             Element<DIM,DIM>* pElement)
00049 {
00050     c_matrix<double,2*(DIM+1),2*(DIM+1)> ret = zero_matrix<double>(2*(DIM+1), 2*(DIM+1));
00051 
00052     if (pElement->GetRegion() == HeartRegionCode::TISSUE)
00053     {          
00054         c_matrix<double, DIM+1, DIM+1> basis_outer_prod = outer_prod(rPhi, rPhi);
00055  
00056         // even rows, even columns
00057         matrix_slice<c_matrix<double, 2*DIM+2, 2*DIM+2> >
00058         slice00(ret, slice (0, 2, DIM+1), slice (0, 2, DIM+1));
00059         slice00 =  basis_outer_prod;
00060 
00061         // odd rows, even columns: are zero
00062         // even rows, odd columns: are zero
00063         // odd rows, odd columns: are zero
00064     }
00065 
00066     return ret;
00067 }
00068 
00069 template<unsigned DIM>
00070 c_vector<double,2*(DIM+1)> BidomainWithBathRhsMatrixAssembler<DIM>::ComputeVectorTerm(
00071     c_vector<double, DIM+1> &rPhi,
00072     c_matrix<double, DIM, DIM+1> &rGradPhi,
00073     ChastePoint<DIM> &rX,
00074     c_vector<double,2> &u,
00075     c_matrix<double, 2, DIM> &rGradU /* not used */,
00076     Element<DIM,DIM>* pElement)
00077 
00078 {
00079     #define COVERAGE_IGNORE
00080     NEVER_REACHED;
00081     return zero_vector<double>(2*(DIM+1));
00082     #undef COVERAGE_IGNORE
00083 }
00084 
00085 
00086 template<unsigned DIM>
00087 c_vector<double, 2*DIM> BidomainWithBathRhsMatrixAssembler<DIM>::ComputeVectorSurfaceTerm(
00088     const BoundaryElement<DIM-1,DIM> &rSurfaceElement,
00089     c_vector<double, DIM> &rPhi,
00090     ChastePoint<DIM> &rX )
00091 {
00092     #define COVERAGE_IGNORE
00093     NEVER_REACHED;
00094     return zero_vector<double>(2*DIM);
00095     #undef COVERAGE_IGNORE
00096 }
00097 
00098 
00099 template<unsigned DIM>
00100 BidomainWithBathRhsMatrixAssembler<DIM>::BidomainWithBathRhsMatrixAssembler(AbstractMesh<DIM,DIM>* pMesh)
00101     : AbstractLinearAssembler<DIM,DIM,2,false,BidomainWithBathRhsMatrixAssembler<DIM> >()
00102 {
00103     this->mpMesh = pMesh;
00104 
00105     // this needs to be set up, though no boundary condition values are used in the matrix
00106     this->mpBoundaryConditions = new BoundaryConditionsContainer<DIM,DIM,2>;
00107     this->mpBoundaryConditions->DefineZeroNeumannOnMeshBoundary(pMesh);
00108 
00109     //DistributedVector::SetProblemSize(this->mpMesh->GetNumNodes()); WOULD BE WRONG -- we need the maintain an uneven distribution, if given
00110     Vec template_vec = DistributedVector::CreateVec(2);
00111     this->mpLinearSystem = new LinearSystem(template_vec);
00112     VecDestroy(template_vec);
00113 
00114 
00115     this->AssembleSystem(false,true);
00116 }
00117 
00118 template<unsigned DIM>
00119 BidomainWithBathRhsMatrixAssembler<DIM>::~BidomainWithBathRhsMatrixAssembler()
00120 {
00121     delete this->mpBoundaryConditions;
00122 }
00123 
00124 template<unsigned DIM>
00125 Mat* BidomainWithBathRhsMatrixAssembler<DIM>::GetMatrix()
00126 {
00127     return &(this->mpLinearSystem->rGetLhsMatrix());
00128 }
00129 
00131 // BidomainMatrixBasedAssembler
00133 
00134 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00135 BidomainWithBathMatrixBasedAssembler<ELEMENT_DIM,SPACE_DIM>::BidomainWithBathMatrixBasedAssembler(
00136             AbstractMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00137             BidomainPde<SPACE_DIM>* pPde,
00138             BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, 2>* pBcc,
00139             unsigned numQuadPoints)
00140     : BidomainDg0Assembler<ELEMENT_DIM,SPACE_DIM>(pMesh, pPde, pBcc, numQuadPoints), 
00141       BidomainMatrixBasedAssembler<ELEMENT_DIM,SPACE_DIM>(pMesh, pPde, pBcc, numQuadPoints),
00142       BidomainWithBathAssembler<ELEMENT_DIM,SPACE_DIM>(pMesh, pPde, pBcc, numQuadPoints)
00143       
00144     
00145 {
00146     // construct matrix using the helper class
00147     mpBidomainWithBathRhsMatrixAssembler = new BidomainWithBathRhsMatrixAssembler<SPACE_DIM>(pMesh);
00148     this->mpMatrixForMatrixBasedRhsAssembly = mpBidomainWithBathRhsMatrixAssembler->GetMatrix();
00149     
00152 }
00153 
00154 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00155 BidomainWithBathMatrixBasedAssembler<ELEMENT_DIM,SPACE_DIM>::~BidomainWithBathMatrixBasedAssembler()
00156 {
00157     delete mpBidomainWithBathRhsMatrixAssembler;    
00158 }
00159 
00160 
00161 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00162 void BidomainWithBathMatrixBasedAssembler<ELEMENT_DIM,SPACE_DIM>::ConstructVectorForMatrixBasedRhsAssembly(
00163         Vec currentSolution)
00164 {
00165     // dist stripe for the current Voltage
00166     DistributedVector distributed_current_solution(currentSolution);
00167     DistributedVector::Stripe distributed_current_solution_vm(distributed_current_solution, 0); 
00168          
00169     // dist stripe for z
00170     DistributedVector dist_vec_matrix_based(this->mVectorForMatrixBasedRhsAssembly);     
00171     DistributedVector::Stripe dist_vec_matrix_based_vm(dist_vec_matrix_based, 0);
00172     DistributedVector::Stripe dist_vec_matrix_based_phie(dist_vec_matrix_based, 1);
00173 
00174     double Am = HeartConfig::Instance()->GetSurfaceAreaToVolumeRatio();
00175     double Cm  = HeartConfig::Instance()->GetCapacitance();
00176     
00177     for (DistributedVector::Iterator index = DistributedVector::Begin();
00178          index!= DistributedVector::End();
00179          ++index)
00180     {
00181         if(this->mpMesh->GetNode(index.Global)->GetRegion() == HeartRegionCode::TISSUE)
00182         {
00183             double V = distributed_current_solution_vm[index];
00184             double F = - Am*this->mpBidomainPde->rGetIionicCacheReplicated()[index.Global] 
00185                        - this->mpBidomainPde->rGetIntracellularStimulusCacheReplicated()[index.Global]; 
00186             
00187             dist_vec_matrix_based_vm[index] = Am*Cm*V*this->mDtInverse + F;
00188         }
00189         else
00190         {
00191             dist_vec_matrix_based_vm[index] = 0.0;
00192         } 
00193 
00194         dist_vec_matrix_based_phie[index] = 0.0;
00195     }
00196 
00197     dist_vec_matrix_based.Restore();
00198     
00199     VecAssemblyBegin(this->mVectorForMatrixBasedRhsAssembly);
00200     VecAssemblyEnd(this->mVectorForMatrixBasedRhsAssembly); 
00201 }
00202 
00203 
00205 // Explicit instantiation
00207 
00208 template class BidomainWithBathMatrixBasedAssembler<1,1>;
00209 template class BidomainWithBathMatrixBasedAssembler<2,2>;
00210 template class BidomainWithBathMatrixBasedAssembler<3,3>;
00211 
00212 template class BidomainWithBathRhsMatrixAssembler<1>;
00213 template class BidomainWithBathRhsMatrixAssembler<2>;
00214 template class BidomainWithBathRhsMatrixAssembler<3>;

Generated on Wed Mar 18 12:51:52 2009 for Chaste by  doxygen 1.5.5