ImplicitCardiacMechanicsAssembler.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 "ImplicitCardiacMechanicsAssembler.hpp"
00030 
00031 template<unsigned DIM>
00032 ImplicitCardiacMechanicsAssembler<DIM>::ImplicitCardiacMechanicsAssembler(
00033                                   QuadraticMesh<DIM>* pQuadMesh,
00034                                   std::string outputDirectory,
00035                                   std::vector<unsigned>& rFixedNodes,
00036                                   AbstractIncompressibleMaterialLaw<DIM>* pMaterialLaw)
00037     : NonlinearElasticityAssembler<DIM>(pQuadMesh, 
00038                                         pMaterialLaw!=NULL ? pMaterialLaw : new NashHunterPoleZeroLaw<DIM>,
00039                                         zero_vector<double>(DIM),
00040                                         DOUBLE_UNSET,
00041                                         outputDirectory, 
00042                                         rFixedNodes),
00043       mCurrentTime(DBL_MAX),
00044       mNextTime(DBL_MAX),
00045       mOdeTimestep(DBL_MAX)
00046 {
00047     // compute total num quad points
00048     mTotalQuadPoints = pQuadMesh->GetNumElements()*this->mpQuadratureRule->GetNumQuadPoints();
00049 
00050     // initialise stores
00051     mLambda.resize(mTotalQuadPoints, 1.0);
00052     mLambdaLastTimeStep.resize(mTotalQuadPoints, 1.0);
00053     mCellMechSystems.resize(mTotalQuadPoints);
00054 
00055     // note that if pMaterialLaw is NULL a new NashHunter law was sent to the 
00056     // NonlinElas constuctor (see above)
00057     mAllocatedMaterialLawMemory = (pMaterialLaw==NULL);
00058 }
00059 
00060 template<unsigned DIM>
00061 ImplicitCardiacMechanicsAssembler<DIM>::~ImplicitCardiacMechanicsAssembler()
00062 {
00063     if(mAllocatedMaterialLawMemory)
00064     {
00065         assert(this->mMaterialLaws.size()==1); // haven't implemented heterogeniety yet
00066         delete this->mMaterialLaws[0];
00067     }
00068 }
00069 
00070 template<unsigned DIM>
00071 unsigned ImplicitCardiacMechanicsAssembler<DIM>::GetTotalNumQuadPoints()
00072 {
00073     return mTotalQuadPoints;
00074 }
00075 
00076 template<unsigned DIM>
00077 GaussianQuadratureRule<DIM>* ImplicitCardiacMechanicsAssembler<DIM>::GetQuadratureRule()
00078 {
00079     return this->mpQuadratureRule;
00080 }
00081 
00082 template<unsigned DIM>
00083 void ImplicitCardiacMechanicsAssembler<DIM>::SetIntracellularCalciumConcentrations(std::vector<double>& caI)
00084 {
00085     assert(caI.size() == mTotalQuadPoints);
00086     for(unsigned i=0; i<caI.size(); i++)
00087     {
00088         mCellMechSystems[i].SetIntracellularCalciumConcentration(caI[i]);
00089     }
00090 }
00091 
00092 template<unsigned DIM>
00093 std::vector<double>& ImplicitCardiacMechanicsAssembler<DIM>::rGetLambda()
00094 {
00095     return mLambda;
00096 }
00097 
00098 
00099 template<unsigned DIM>
00100 void ImplicitCardiacMechanicsAssembler<DIM>::Solve(double currentTime, double nextTime, double odeTimestep)
00101 {
00102     // set the times, which are used in AssembleOnElement
00103     assert(currentTime < nextTime);
00104     mCurrentTime = currentTime;
00105     mNextTime = nextTime;
00106     mOdeTimestep = odeTimestep;
00107 
00108     // solve
00109     NonlinearElasticityAssembler<DIM>::Solve();
00110 
00111     // assemble residual again (to solve the cell models implicitly again 
00112     // using the correct value of the deformation x (in case this wasn't the
00113     // last thing that was done 
00114     this->AssembleSystem(true,false);
00115 
00116     // now update state variables, and set lambda at last timestep. Note
00117     // lambda was set in AssembleOnElement
00118     for(unsigned i=0; i<mCellMechSystems.size(); i++)
00119     {
00120          mCellMechSystems[i].UpdateStateVariables();
00121          mLambdaLastTimeStep[i] = mCellMechSystems[i].GetLambda();
00122     }
00123 }
00124 
00125 
00126 
00127 template<unsigned DIM>
00128 void ImplicitCardiacMechanicsAssembler<DIM>::AssembleOnElement(Element<DIM, DIM>& rElement,
00129                        c_matrix<double,STENCIL_SIZE,STENCIL_SIZE>& rAElem,
00130                        c_matrix<double,STENCIL_SIZE,STENCIL_SIZE>& rAElemPrecond,
00131                        c_vector<double,STENCIL_SIZE>& rBElem,
00132                        bool assembleResidual,
00133                        bool assembleJacobian)
00134 {
00135     // check these have been set
00136     assert(mCurrentTime != DBL_MAX);
00137     assert(mNextTime != DBL_MAX);
00138     assert(mOdeTimestep != DBL_MAX);
00139     
00140     c_matrix<double, DIM, DIM> jacobian, inverse_jacobian;
00141     double jacobian_determinant;
00142     this->mpQuadMesh->GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
00143 
00144     if (assembleJacobian)
00145     {
00146         rAElem.clear();
00147         rAElemPrecond.clear();
00148     }
00149 
00150     if (assembleResidual)
00151     {
00152         rBElem.clear();
00153     }
00154 
00156     // Get the current displacement at the nodes
00158     static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
00159     static c_vector<double,NUM_VERTICES_PER_ELEMENT> element_current_pressures;
00160     for(unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
00161     {
00162         for(unsigned JJ=0; JJ<DIM; JJ++)
00163         {
00164             element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*rElement.GetNodeGlobalIndex(II) + JJ];
00165         }
00166     }
00167 
00169     // Get the current pressure at the vertices
00171     for(unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
00172     {
00173         element_current_pressures(II) = this->mCurrentSolution[DIM*this->mpQuadMesh->GetNumNodes() + rElement.GetNodeGlobalIndex(II)];
00174     }
00175 
00176     // allocate memory for the basis functions values and derivative values
00177     c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
00178     c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
00179     c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
00180 
00181     // get the material law
00182     AbstractIncompressibleMaterialLaw<DIM>* p_material_law;
00183     if(this->mMaterialLaws.size()==1)
00184     {
00185         // homogeneous
00186         p_material_law = this->mMaterialLaws[0];
00187     }
00188     else
00189     {
00190         // heterogeneous
00191         #define COVERAGE_IGNORE // not going to have tests in cts for everything
00192         p_material_law = this->mMaterialLaws[rElement.GetIndex()];
00193         #undef COVERAGE_IGNORE
00194     }
00195     
00197 //    assert(DIM==2);
00198 //    double   theta = 0.785398163/5 * elementIter->vertex(0)[0]; //0->pi/20
00199 //    this->mFibreSheetMat[0][0] =  cos(theta);
00200 //    this->mFibreSheetMat[0][1] =  sin(theta);
00201 //    this->mFibreSheetMat[1][0] = -sin(theta);
00202 //    this->mFibreSheetMat[1][1] =  cos(theta);
00203 //    this->mTransFibreSheetMat = transpose(this->mFibreSheetMat);
00204 
00205 
00211     for (unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
00212     {
00213         unsigned current_quad_point_global_index =   rElement.GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
00214                                                    + quadrature_index;
00215 
00216         double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
00217 
00218         const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
00219 
00221         // set up basis function info
00223         LinearBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, linear_phi);
00224         QuadraticBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, quad_phi);
00225         QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
00226         
00228         // (dont get the body force)
00230 
00232         // interpolate grad_u and p
00234         static c_matrix<double,DIM,DIM> grad_u; // grad_u = (du_i/dX_M)
00235         grad_u = zero_matrix<double>(DIM,DIM);  // must be on new line!!
00236 
00237         for(unsigned node_index=0; 
00238             node_index<NUM_NODES_PER_ELEMENT; 
00239             node_index++)
00240         {
00241             for (unsigned i=0; i<DIM; i++)
00242             {
00243                 for(unsigned M=0; M<DIM; M++)
00244                 {
00245                     grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
00246                 }
00247             }
00248         }
00249 
00250         double pressure = 0;
00251         for(unsigned vertex_index=0;
00252             vertex_index<NUM_VERTICES_PER_ELEMENT;
00253             vertex_index++)
00254         {
00255             pressure += linear_phi(vertex_index)*element_current_pressures(vertex_index);
00256         }
00257 
00259         // calculate C, inv(C) and T
00261         static c_matrix<double,DIM,DIM> F;
00262         static c_matrix<double,DIM,DIM> C;
00263         static c_matrix<double,DIM,DIM> inv_C;
00264         static c_matrix<double,DIM,DIM> inv_F;
00265         static c_matrix<double,DIM,DIM> T;
00266 
00267         for (unsigned i=0; i<DIM; i++)
00268         {
00269             for (unsigned M=0; M<DIM; M++)
00270             {
00271                 F(i,M) = (i==M?1:0) + grad_u(i,M);
00272             }
00273         }
00274 
00275         C = prod(trans(F),F);
00276         inv_C = Inverse(C);
00277         inv_F = Inverse(F);
00278 
00279         double detF = Determinant(F);
00280 
00281         /************************************
00282          *  The cardiac-specific code PART 1
00283          ************************************/
00284         //static c_matrix<double,DIM,DIM> C_fibre;          // C when transformed to fibre-sheet axes
00285         //static c_matrix<double,DIM,DIM> inv_C_fibre;      // C^{-1} transformed to fibre-sheet axes
00286         //static c_matrix<double,DIM,DIM> T_fibre; // T when transformed to fibre-sheet axes
00287         //
00289         //C_fibre = this->mTransFibreSheetMat * C * this->mFibreSheetMat;
00290         //inv_C_fibre = this->mTransFibreSheetMat * inv_C * this->mFibreSheetMat;
00291 
00292         // store the stretch in the fibre direction
00293         this->mLambda[current_quad_point_global_index] = sqrt(C(0,0));
00294 
00295         double lam = sqrt(C(0,0));
00296         double dlam_dt = (lam-mLambdaLastTimeStep[current_quad_point_global_index])/(mNextTime-mCurrentTime);
00297 
00298         NhsSystemWithImplicitSolver& system = mCellMechSystems[current_quad_point_global_index];
00299 
00300         // get proper active tension
00301         // see NOTE below
00302         system.SetLambdaAndDerivative(lam, dlam_dt);
00303 
00304         double active_tension=0.0;
00305         try
00306         {
00307             system.SolveDoNotUpdate(mCurrentTime,mNextTime,mOdeTimestep);
00308             active_tension = system.GetActiveTensionAtNextTime(); 
00309         }
00310         catch (Exception& e)
00311         {
00312             #define COVERAGE_IGNORE
00313             if(assembleJacobian) 
00314             {
00315                 EXCEPTION("Failed in solving NHS systems when assembling Jacobian");
00316             }
00317             #undef COVERAGE_IGNORE
00318         }
00319 
00320         // compute the derivative of the active tension wrt lam and dlam_dt
00321         double d_act_tension_dlam = 0.0; //Set and used if assembleJacobian==true
00322         double d_act_tension_d_dlamdt = 0.0; //Set and used if assembleJacobian==true
00323 
00324         if(assembleJacobian)
00325         {
00326             // get active tension for (lam+h,dlamdt)
00327             double h1 = std::max(1e-6, lam/100);
00328             system.SetLambdaAndDerivative(lam+h1, dlam_dt);
00329             system.SolveDoNotUpdate(mCurrentTime,mNextTime,mOdeTimestep);
00330             double active_tension_at_lam_plus_h = system.GetActiveTensionAtNextTime();
00331 
00332             // get active tension for (lam,dlamdt+h)
00333             double h2 = std::max(1e-6, dlam_dt/100);
00334             system.SetLambdaAndDerivative(lam, dlam_dt+h2);
00335             system.SolveDoNotUpdate(mCurrentTime,mNextTime,mOdeTimestep);
00336             double active_tension_at_dlamdt_plus_h = system.GetActiveTensionAtNextTime();
00337 
00338             d_act_tension_dlam = (active_tension_at_lam_plus_h - active_tension)/h1;
00339             d_act_tension_d_dlamdt = (active_tension_at_dlamdt_plus_h - active_tension)/h2;
00340         }
00341 
00342         // NOTE - have to get the active tension again, this must be done last!!
00343         // As if this turns out to be the correct solution, the state vars will be updated!
00344         // TODO: sort out this inefficiency
00345         system.SetLambdaAndDerivative(lam, dlam_dt);
00346         system.SetActiveTensionInitialGuess(active_tension);
00347 
00348         try
00349         {
00350             system.SolveDoNotUpdate(mCurrentTime,mNextTime,mOdeTimestep);
00351             assert( fabs(system.GetActiveTensionAtNextTime()-active_tension)<1e-8);
00352         }
00353         catch (Exception& e)
00354         {
00356             #define COVERAGE_IGNORE
00357             LOG(2, "WARNING in ImplicitCardiacMechanicsAssembler!\n");
00358             active_tension = 1e10;
00359             // should have done something above..
00360             #undef COVERAGE_IGNORE
00361         }
00362 
00363         //this->mDTdE_fibre.Zero();
00364 
00365         // compute the transformed tension. The material law should law be a cardiac-
00366         // specific law which assumes the x-axes in the fibre, the z-axes the sheet normal
00367         p_material_law->ComputeStressAndStressDerivative(C,inv_C,pressure,T,this->dTdE,assembleJacobian);
00368 
00369         // amend the stress and dTdE using the active tension
00370         T(0,0) += active_tension/C(0,0);
00371         this->dTdE(0,0,0,0) -= 2*active_tension/(C(0,0)*C(0,0));
00372 
00373 //            //// could do this without the loop now
00374 //            // transform T back to real coordinates
00375 //            for(unsigned M=0; M<DIM; M++)
00376 //            {
00377 //                for(unsigned N=0; N<DIM; N++)
00378 //                {
00379 //                    T[M][N] = 0;
00380 //                    for(unsigned al=0; al<DIM; al++)
00381 //                    {
00382 //                        for(unsigned be=0; be<DIM; be++)
00383 //                        {
00384 //                            T[M][N] +=                      T_fibre [al][be]
00385 //                                        *      this->mFibreSheetMat [M] [al]
00386 //                                        * this->mTransFibreSheetMat [be][N];
00387 //                        }
00388 //                    }
00389 //                }
00390 //            }
00391 //            static FourthOrderTensor<DIM> temp1;
00392 //            static FourthOrderTensor<DIM> temp2;
00393 //            static FourthOrderTensor<DIM> temp3;
00394 //
00395 //            temp1.SetAsProduct(this->mDTdE_fibre, this->mFibreSheetMat, 0);
00396 //            temp2.SetAsProduct(temp1,             this->mFibreSheetMat, 1);
00397 //            temp3.SetAsProduct(temp2,             this->mFibreSheetMat, 2);
00398 //
00399 //            this->dTdE.SetAsProduct(temp3, this->mFibreSheetMat, 3);
00400 
00401 
00402         /*************************************
00403          * end of cardiac specific code PART 1
00404          *************************************/
00405 
00407         // residual vector
00409         if (assembleResidual)
00410         {
00411             for(unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
00412             {
00413                 unsigned spatial_dim = index%DIM;
00414                 unsigned node_index = (index-spatial_dim)/DIM;
00415 
00416                 assert(node_index < NUM_NODES_PER_ELEMENT);
00417 
00418                 // no body force bit as body force = 0
00419 
00420                 for (unsigned M=0; M<DIM; M++)
00421                 {
00422                     for (unsigned N=0; N<DIM; N++)
00423                     {
00424                         rBElem(index) +=   T(M,N)
00425                                          * F(spatial_dim,M)
00426                                          * grad_quad_phi(N,node_index)
00427                                          * wJ;
00428                     }
00429                 }
00430             }
00431             
00432             for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00433             {
00434                 rBElem( NUM_NODES_PER_ELEMENT*DIM + vertex_index ) +=   linear_phi(vertex_index)
00435                                                                       * (detF - 1)
00436                                                                       * wJ;
00437             }
00438         }
00439 
00441         // Jacobian matrix
00443         if(assembleJacobian)
00444         {
00445             for(unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
00446             {
00447                 unsigned spatial_dim1 = index1%DIM;
00448                 unsigned node_index1 = (index1-spatial_dim1)/DIM;
00449                 
00450                 
00451                 for(unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00452                 {
00453                     unsigned spatial_dim2 = index2%DIM;
00454                     unsigned node_index2 = (index2-spatial_dim2)/DIM;
00455 
00456                     for (unsigned M=0; M<DIM; M++)
00457                     {
00458                         for (unsigned N=0; N<DIM; N++)
00459                         {
00460                             rAElem(index1,index2) +=   T(M,N)
00461                                                      * grad_quad_phi(N,node_index1)
00462                                                      * grad_quad_phi(M,node_index2)
00463                                                      * (spatial_dim1==spatial_dim2?1:0)
00464                                                      * wJ;
00465 
00466                             for (unsigned P=0; P<DIM; P++)
00467                             {
00468                                 for (unsigned Q=0; Q<DIM; Q++)
00469                                 {
00470                                     rAElem(index1,index2)  +=   0.5
00471                                                               * this->dTdE(M,N,P,Q)
00472                                                               * (
00473                                                                   grad_quad_phi(P,node_index2)
00474                                                                 * F(spatial_dim2,Q)
00475                                                                    +
00476                                                                   grad_quad_phi(Q,node_index2)
00477                                                                 * F(spatial_dim2,P)
00478                                                                  )
00479                                                               * F(spatial_dim1,M)
00480                                                               * grad_quad_phi(N,node_index1)
00481                                                               * wJ;
00482                                 }
00483                             }
00484                         }
00485                     }
00486                     
00487                     /************************************
00488                      *  The cardiac-specific code PART 2
00489                      ************************************/
00490                     rAElem(index1,index2) +=  (
00491                                                    d_act_tension_dlam
00492                                                  +
00493                                                    d_act_tension_d_dlamdt/(mNextTime-mCurrentTime)
00494                                                 )
00495                                                 * (F(spatial_dim2,0)/lam)
00496                                                 * grad_quad_phi(0,node_index2)
00497                                                 //* (fe_values.shape_grad(j,q_point)[0]/C[0][0])
00498                                                 * F(spatial_dim1,0)
00499                                                 * grad_quad_phi(0,node_index1)
00500                                                 //* fe_values.shape_grad(i,q_point)[0]
00501                                                 * wJ;                        
00502                    /************************************
00503                     *  End cardiac-specific code PART 2
00504                     ************************************/
00505                 }
00506                 
00507                 
00508                 for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00509                 {
00510                     unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00511                     
00512                     for (unsigned M=0; M<DIM; M++)
00513                     {
00514                         for (unsigned N=0; N<DIM; N++)
00515                         {
00516                             rAElem(index1,index2)  +=  - F(spatial_dim1,M)
00517                                                        * inv_C(M,N)
00518                                                        * grad_quad_phi(N,node_index1)
00519                                                        * linear_phi(vertex_index)
00520                                                        * wJ;
00521                         }
00522                     }                 
00523                 }
00524             }
00525 
00526             for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00527             {
00528                 unsigned index1 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00529 
00530                 for(unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00531                 {
00532                     unsigned spatial_dim2 = index2%DIM;
00533                     unsigned node_index2 = (index2-spatial_dim2)/DIM;
00534 
00535                     for (unsigned M=0; M<DIM; M++)
00536                     {
00537                         rAElem(index1,index2) +=   linear_phi(vertex_index)
00538                                                  * detF
00539                                                  * inv_F(M,spatial_dim2)
00540                                                  * grad_quad_phi(M,node_index2)
00541                                                  * wJ;
00542                     }
00543                 }
00544 
00546                 // Preconditioner matrix
00547                 // Fill the mass matrix (ie \intgl phi_i phi_j) in the 
00548                 // pressure-pressure block. Note, the rest of the 
00549                 // entries are filled in at the end
00551                 for(unsigned vertex_index2=0; vertex_index2< NUM_VERTICES_PER_ELEMENT; vertex_index2++) 
00552                 { 
00553                     unsigned index2 =  NUM_NODES_PER_ELEMENT*DIM + vertex_index2;
00554                     rAElemPrecond(index1,index2) +=   linear_phi(vertex_index)
00555                                                     * linear_phi(vertex_index2) 
00556                                                     * wJ; 
00557                 } 
00558             }
00559         }
00560     }
00561 
00562 
00563     if (assembleJacobian) 
00564     { 
00565         // Fill in the other blocks of the preconditioner matrix. (This doesn't
00566         // effect the pressure-pressure block of the rAElemPrecond but the 
00567         // pressure-pressure block of rAElem is zero
00568         rAElemPrecond = rAElemPrecond + rAElem; 
00569     }
00570 }
00571 
00572 template class ImplicitCardiacMechanicsAssembler<2>;
00573 template class ImplicitCardiacMechanicsAssembler<3>;
00574 
00575 

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