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

Generated on Tue Aug 4 16:10:22 2009 for Chaste by  doxygen 1.5.5