00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00048 mTotalQuadPoints = pQuadMesh->GetNumElements()*this->mpQuadratureRule->GetNumQuadPoints();
00049
00050
00051 mLambda.resize(mTotalQuadPoints, 1.0);
00052 mLambdaLastTimeStep.resize(mTotalQuadPoints, 1.0);
00053 mCellMechSystems.resize(mTotalQuadPoints);
00054
00055
00056
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);
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 time, double nextTime, double odeTimestep)
00101 {
00102
00103 assert(time < nextTime);
00104 mCurrentTime = time;
00105 mNextTime = nextTime;
00106 mOdeTimestep = odeTimestep;
00107
00108
00109 NonlinearElasticityAssembler<DIM>::Solve();
00110
00111
00112
00113
00114 this->AssembleSystem(true,false);
00115
00116
00117
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
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
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
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
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
00182 AbstractIncompressibleMaterialLaw<DIM>* p_material_law;
00183 if(this->mMaterialLaws.size()==1)
00184 {
00185
00186 p_material_law = this->mMaterialLaws[0];
00187 }
00188 else
00189 {
00190
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
00198
00199
00200
00201
00202
00203
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
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
00230
00232
00234 static c_matrix<double,DIM,DIM> grad_u;
00235 grad_u = zero_matrix<double>(DIM,DIM);
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
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
00283
00284
00285
00286
00287
00289
00290
00291
00292
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
00301
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
00321 double d_act_tension_dlam = 0.0;
00322 double d_act_tension_d_dlamdt = 0.0;
00323
00324 if(assembleJacobian)
00325 {
00326
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
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
00343
00344
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 {
00355 #define COVERAGE_IGNORE
00356 LOG(2, "WARNING in ImplicitCardiacMechanicsAssembler!\n");
00357 active_tension = 1e10;
00358
00359 #undef COVERAGE_IGNORE
00360 }
00361
00362
00363
00364
00365
00366 p_material_law->ComputeStressAndStressDerivative(C,inv_C,pressure,T,this->dTdE,assembleJacobian);
00367
00368
00369 T(0,0) += active_tension/C(0,0);
00370 this->dTdE(0,0,0,0) -= 2*active_tension/(C(0,0)*C(0,0));
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402 static FourthOrderTensor<DIM> dTdE_F;
00403 static FourthOrderTensor<DIM> dTdE_FF1;
00404 static FourthOrderTensor<DIM> dTdE_FF2;
00405
00406 dTdE_F.SetAsProduct(this->dTdE, F, 0);
00407 dTdE_FF1.SetAsProduct(dTdE_F, F, 3);
00408 dTdE_FF2.SetAsProduct(dTdE_F, F, 2);
00409
00410
00411
00412
00413
00414
00416
00418 if (assembleResidual)
00419 {
00420 for(unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
00421 {
00422 unsigned spatial_dim = index%DIM;
00423 unsigned node_index = (index-spatial_dim)/DIM;
00424
00425 assert(node_index < NUM_NODES_PER_ELEMENT);
00426
00427
00428
00429 for (unsigned M=0; M<DIM; M++)
00430 {
00431 for (unsigned N=0; N<DIM; N++)
00432 {
00433 rBElem(index) += T(M,N)
00434 * F(spatial_dim,M)
00435 * grad_quad_phi(N,node_index)
00436 * wJ;
00437 }
00438 }
00439 }
00440
00441 for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00442 {
00443 rBElem( NUM_NODES_PER_ELEMENT*DIM + vertex_index ) += linear_phi(vertex_index)
00444 * (detF - 1)
00445 * wJ;
00446 }
00447 }
00448
00450
00452 if(assembleJacobian)
00453 {
00454 for(unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
00455 {
00456 unsigned spatial_dim1 = index1%DIM;
00457 unsigned node_index1 = (index1-spatial_dim1)/DIM;
00458
00459
00460 for(unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00461 {
00462 unsigned spatial_dim2 = index2%DIM;
00463 unsigned node_index2 = (index2-spatial_dim2)/DIM;
00464
00465 for (unsigned M=0; M<DIM; M++)
00466 {
00467 for (unsigned N=0; N<DIM; N++)
00468 {
00469 rAElem(index1,index2) += T(M,N)
00470 * grad_quad_phi(N,node_index1)
00471 * grad_quad_phi(M,node_index2)
00472 * (spatial_dim1==spatial_dim2?1:0)
00473 * wJ;
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493 }
00494 }
00495
00496 for (unsigned N=0; N<DIM; N++)
00497 {
00498 for (unsigned P=0; P<DIM; P++)
00499 {
00500 rAElem(index1,index2) += 0.5
00501 * dTdE_FF1(spatial_dim1,N,P,spatial_dim2)
00502 * grad_quad_phi(P,node_index2)
00503 * grad_quad_phi(N,node_index1)
00504 * wJ;
00505 }
00506
00507 for (unsigned Q=0; Q<DIM; Q++)
00508 {
00509 rAElem(index1,index2) += 0.5
00510 * dTdE_FF2(spatial_dim1,N,spatial_dim2,Q)
00511 * grad_quad_phi(Q,node_index2)
00512 * grad_quad_phi(N,node_index1)
00513 * wJ;
00514 }
00515 }
00516
00517
00518
00519
00520
00521
00522
00523
00524 rAElem(index1,index2) += (
00525 d_act_tension_dlam
00526 +
00527 d_act_tension_d_dlamdt/(mNextTime-mCurrentTime)
00528 )
00529 * (F(spatial_dim2,0)/lam)
00530 * grad_quad_phi(0,node_index2)
00531 * F(spatial_dim1,0)
00532 * grad_quad_phi(0,node_index1)
00533 * wJ;
00534
00535
00536
00537 }
00538
00539
00540 for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00541 {
00542 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00543
00544 for (unsigned M=0; M<DIM; M++)
00545 {
00546 for (unsigned N=0; N<DIM; N++)
00547 {
00548 rAElem(index1,index2) += - F(spatial_dim1,M)
00549 * inv_C(M,N)
00550 * grad_quad_phi(N,node_index1)
00551 * linear_phi(vertex_index)
00552 * wJ;
00553 }
00554 }
00555 }
00556 }
00557
00558 for(unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
00559 {
00560 unsigned index1 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
00561
00562 for(unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
00563 {
00564 unsigned spatial_dim2 = index2%DIM;
00565 unsigned node_index2 = (index2-spatial_dim2)/DIM;
00566
00567 for (unsigned M=0; M<DIM; M++)
00568 {
00569 rAElem(index1,index2) += linear_phi(vertex_index)
00570 * detF
00571 * inv_F(M,spatial_dim2)
00572 * grad_quad_phi(M,node_index2)
00573 * wJ;
00574 }
00575 }
00576
00578
00579
00580
00581
00583 for(unsigned vertex_index2=0; vertex_index2< NUM_VERTICES_PER_ELEMENT; vertex_index2++)
00584 {
00585 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index2;
00586 rAElemPrecond(index1,index2) += linear_phi(vertex_index)
00587 * linear_phi(vertex_index2)
00588 * wJ;
00589 }
00590 }
00591 }
00592 }
00593
00594
00595 if (assembleJacobian)
00596 {
00597
00598
00599
00600 rAElemPrecond = rAElemPrecond + rAElem;
00601 }
00602 }
00603
00604 template class ImplicitCardiacMechanicsAssembler<2>;
00605 template class ImplicitCardiacMechanicsAssembler<3>;
00606
00607