Chaste  Release::2018.1
AbstractCardiacMechanicsSolver.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "AbstractCardiacMechanicsSolver.hpp"
37 #include "AbstractContractionCellFactory.hpp"
38 #include "FakeBathContractionModel.hpp"
39 
40 template<class ELASTICITY_SOLVER,unsigned DIM>
42  ElectroMechanicsProblemDefinition<DIM>& rProblemDefinition,
43  std::string outputDirectory)
44  : ELASTICITY_SOLVER(rQuadMesh,
45  rProblemDefinition,
46  outputDirectory),
47  mpMeshPair(NULL),
48  mCurrentTime(DBL_MAX),
49  mNextTime(DBL_MAX),
50  mOdeTimestep(DBL_MAX),
51  mrElectroMechanicsProblemDefinition(rProblemDefinition)
52 {
53 }
54 
55 template<class ELASTICITY_SOLVER,unsigned DIM>
57 {
58  // compute total num quad points
59  unsigned num_quad_pts_per_element = this->mpQuadratureRule->GetNumQuadPoints();
60  mTotalQuadPoints = this->mrQuadMesh.GetNumElements()*num_quad_pts_per_element;
61 
62  std::vector<ElementAndWeights<DIM> > fine_elements = mpMeshPair->rGetElementsAndWeights();
63  assert(fine_elements.size()==mTotalQuadPoints);
64  assert(mpMeshPair!=NULL);
65 
66  AbstractContractionCellFactory<DIM>* p_factory = mrElectroMechanicsProblemDefinition.GetContractionCellFactory();
67 
69  iter != this->mrQuadMesh.GetElementIteratorEnd();
70  ++iter)
71  {
72  Element<DIM, DIM>& element = *iter;
73 
74  if (element.GetOwnership() == true)
75  {
76  for (unsigned j=0; j<num_quad_pts_per_element; j++)
77  {
78  unsigned quad_pt_global_index = element.GetIndex()*num_quad_pts_per_element + j;
79 
80  // We construct a set of data to be assigned to each quadrature point
81  // this includes a contraction cell model set as bath or by the contraction
82  // cell factory.
83  DataAtQuadraturePoint data_at_quad_point;
84  data_at_quad_point.Stretch = 1.0;
85  data_at_quad_point.StretchLastTimeStep = 1.0;
86 
87  if (mpMeshPair->GetFineMesh().GetElement(fine_elements[quad_pt_global_index].ElementNum)
88  ->GetUnsignedAttribute() == HeartRegionCode::GetValidBathId() )
89  {
90  // Bath
91  data_at_quad_point.ContractionModel = new FakeBathContractionModel;
92  }
93  else
94  {
95  // Tissue
96  data_at_quad_point.ContractionModel = p_factory->CreateContractionCellForElement( &element );
97  }
98  mQuadPointToDataAtQuadPointMap[quad_pt_global_index] = data_at_quad_point;
99  }
100  }
101  }
102 
103  // initialise the iterator to point at the beginning
104  mMapIterator = mQuadPointToDataAtQuadPointMap.begin();
105 
106  // initialise fibre/sheet direction matrix to be the identity, fibres in X-direction, and sheet in XY-plane
107  mConstantFibreSheetDirections = zero_matrix<double>(DIM,DIM);
108  for (unsigned i=0; i<DIM; i++)
109  {
110  mConstantFibreSheetDirections(i,i) = 1.0;
111  }
112 
113  mpVariableFibreSheetDirections = NULL;
114 
115  // Check that we are using the right kind of solver.
116  for (std::map<unsigned,DataAtQuadraturePoint>::iterator iter = this->mQuadPointToDataAtQuadPointMap.begin();
117  iter != this->mQuadPointToDataAtQuadPointMap.end();
118  iter++)
119  {
120  if (!IsImplicitSolver() && (*iter).second.ContractionModel->IsStretchRateDependent())
121  {
122  EXCEPTION("stretch-rate-dependent contraction model requires an IMPLICIT cardiac mechanics solver.");
123  }
124 
125  if (!IsImplicitSolver() && (*iter).second.ContractionModel->IsStretchDependent())
126  {
127  WARN_ONCE_ONLY("stretch-dependent contraction model may require an IMPLICIT cardiac mechanics solver.");
128  }
129  }
130 }
131 
132 
133 template<class ELASTICITY_SOLVER,unsigned DIM>
135 {
136  assert(pMeshPair!=NULL);
137  if (pMeshPair->GetCoarseMesh().GetNumElements() != this->mrQuadMesh.GetNumElements())
138  {
139  EXCEPTION("When setting a mesh pair into the solver, the coarse mesh of the mesh pair must be the same as the quadratic mesh");
140  }
141  mpMeshPair = pMeshPair;
142 }
143 
144 template<class ELASTICITY_SOLVER,unsigned DIM>
146 {
147  for (mMapIterator = mQuadPointToDataAtQuadPointMap.begin();
148  mMapIterator != mQuadPointToDataAtQuadPointMap.end();
149  ++mMapIterator)
150  {
151  AbstractContractionModel* p_model = mMapIterator->second.ContractionModel;
152  if (p_model)
153  {
154  delete p_model;
155  }
156  }
157 
158  if (mpVariableFibreSheetDirections)
159  {
160  delete mpVariableFibreSheetDirections;
161  }
162 }
163 
164 template<class ELASTICITY_SOLVER,unsigned DIM>
166  std::vector<double>& rVoltages)
167 
168 {
169  assert(rCalciumConcentrations.size() == mTotalQuadPoints);
170  assert(rVoltages.size() == mTotalQuadPoints);
171 
172  ContractionModelInputParameters input_parameters;
173 
174  for (unsigned i=0; i<rCalciumConcentrations.size(); i++)
175  {
176  input_parameters.intracellularCalciumConcentration = rCalciumConcentrations[i];
177  input_parameters.voltage = rVoltages[i];
178 
180  std::map<unsigned,DataAtQuadraturePoint>::iterator iter = mQuadPointToDataAtQuadPointMap.find(i);
181  if (iter != mQuadPointToDataAtQuadPointMap.end())
182  {
183  iter->second.ContractionModel->SetInputParameters(input_parameters);
184  }
185  }
186 }
187 
188 
189 template<class ELASTICITY_SOLVER,unsigned DIM>
191  unsigned currentQuadPointGlobalIndex)
192 {
193  if (!mpVariableFibreSheetDirections) // constant fibre directions
194  {
195  this->mChangeOfBasisMatrix = mConstantFibreSheetDirections;
196  }
197  else if (!mFibreSheetDirectionsDefinedByQuadraturePoint) // fibre directions defined for each mechanics mesh element
198  {
199  this->mChangeOfBasisMatrix = (*mpVariableFibreSheetDirections)[elementIndex];
200  }
201  else // fibre directions defined for each mechanics mesh quadrature point
202  {
203  this->mChangeOfBasisMatrix = (*mpVariableFibreSheetDirections)[currentQuadPointGlobalIndex];
204  }
205 }
206 
207 template<class ELASTICITY_SOLVER,unsigned DIM>
209  unsigned elementIndex,
210  unsigned currentQuadPointGlobalIndex,
211  c_matrix<double,DIM,DIM>& rT,
213  bool addToDTdE)
214 {
215  for (unsigned i=0; i<DIM; i++)
216  {
217  mCurrentElementFibreDirection(i) = this->mChangeOfBasisMatrix(i,0);
218  }
219 
220  //Compute the active tension and add to the stress and stress-derivative
221  double I4_fibre = inner_prod(mCurrentElementFibreDirection, prod(rC, mCurrentElementFibreDirection));
222  double lambda_fibre = sqrt(I4_fibre);
223 
224  double active_tension = 0;
225  double d_act_tension_dlam = 0.0; // Set and used if assembleJacobian==true
226  double d_act_tension_d_dlamdt = 0.0; // Set and used if assembleJacobian==true
227 
228  GetActiveTensionAndTensionDerivs(lambda_fibre, currentQuadPointGlobalIndex, addToDTdE,
229  active_tension, d_act_tension_dlam, d_act_tension_d_dlamdt);
230 
231 
232  double detF = sqrt(Determinant(rC));
233  rT += (active_tension*detF/I4_fibre)*outer_prod(mCurrentElementFibreDirection,mCurrentElementFibreDirection);
234 
235  // amend the stress and dTdE using the active tension
236  double dTdE_coeff1 = -2*active_tension*detF/(I4_fibre*I4_fibre); // note: I4_fibre*I4_fibre = lam^4
237  double dTdE_coeff2 = active_tension*detF/I4_fibre;
238  double dTdE_coeff_s1 = 0.0; // only set non-zero if we apply cross fibre tension (in 2/3D)
239  double dTdE_coeff_s2 = 0.0; // only set non-zero if we apply cross fibre tension (in 2/3D)
240  double dTdE_coeff_s3 = 0.0; // only set non-zero if we apply cross fibre tension and implicit (in 2/3D)
241  double dTdE_coeff_n1 = 0.0; // only set non-zero if we apply cross fibre tension in 3D
242  double dTdE_coeff_n2 = 0.0; // only set non-zero if we apply cross fibre tension in 3D
243  double dTdE_coeff_n3 = 0.0; // only set non-zero if we apply cross fibre tension in 3D and implicit
244 
245  if (IsImplicitSolver())
246  {
247  double dt = mNextTime-mCurrentTime;
248  //std::cout << "d sigma / d lamda = " << d_act_tension_dlam << ", d sigma / d lamdat = " << d_act_tension_d_dlamdt << "\n" << std::flush;
249  dTdE_coeff1 += (d_act_tension_dlam + d_act_tension_d_dlamdt/dt)*detF/(lambda_fibre*I4_fibre); // note: I4_fibre*lam = lam^3
250  }
251 
252  bool apply_cross_fibre_tension = (this->mrElectroMechanicsProblemDefinition.GetApplyCrossFibreTension()) && (DIM > 1);
253  if (apply_cross_fibre_tension)
254  {
255  double sheet_cross_fraction = mrElectroMechanicsProblemDefinition.GetSheetTensionFraction();
256 
257  for (unsigned i=0; i<DIM; i++)
258  {
259  mCurrentElementSheetDirection(i) = this->mChangeOfBasisMatrix(i,1);
260  }
261 
262  double I4_sheet = inner_prod(mCurrentElementSheetDirection, prod(rC, mCurrentElementSheetDirection));
263 
264  // amend the stress and dTdE using the active tension
265  dTdE_coeff_s1 = -2*sheet_cross_fraction*detF*active_tension/(I4_sheet*I4_sheet); // note: I4*I4 = lam^4
266 
267  if (IsImplicitSolver())
268  {
269  double dt = mNextTime-mCurrentTime;
270  dTdE_coeff_s3 = sheet_cross_fraction*(d_act_tension_dlam + d_act_tension_d_dlamdt/dt)*detF/(lambda_fibre*I4_sheet); // note: I4*lam = lam^3
271  }
272 
273  rT += sheet_cross_fraction*(active_tension*detF/I4_sheet)*outer_prod(mCurrentElementSheetDirection,mCurrentElementSheetDirection);
274 
275  dTdE_coeff_s2 = active_tension*sheet_cross_fraction*detF/I4_sheet;
276 
277  if (DIM>2)
278  {
279  double sheet_normal_cross_fraction = mrElectroMechanicsProblemDefinition.GetSheetNormalTensionFraction();
280  for (unsigned i=0; i<DIM; i++)
281  {
282  mCurrentElementSheetNormalDirection(i) = this->mChangeOfBasisMatrix(i,2);
283  }
284 
285  double I4_sheet_normal = inner_prod(mCurrentElementSheetNormalDirection, prod(rC, mCurrentElementSheetNormalDirection));
286 
287  dTdE_coeff_n1 =-2*sheet_normal_cross_fraction*detF*active_tension/(I4_sheet_normal*I4_sheet_normal); // note: I4*I4 = lam^4
288 
289  rT += sheet_normal_cross_fraction*(active_tension*detF/I4_sheet_normal)*outer_prod(mCurrentElementSheetNormalDirection,mCurrentElementSheetNormalDirection);
290 
291  dTdE_coeff_n2 = active_tension*sheet_normal_cross_fraction*detF/I4_sheet_normal;
292  if (IsImplicitSolver())
293  {
294  double dt = mNextTime-mCurrentTime;
295  dTdE_coeff_n3 = sheet_normal_cross_fraction*(d_act_tension_dlam + d_act_tension_d_dlamdt/dt)*detF/(lambda_fibre*I4_sheet_normal); // note: I4*lam = lam^3
296  }
297  }
298  }
299 
300 
301  if (addToDTdE)
302  {
303  c_matrix<double,DIM,DIM> invC = Inverse(rC);
304 
305  for (unsigned M=0; M<DIM; M++)
306  {
307  for (unsigned N=0; N<DIM; N++)
308  {
309  for (unsigned P=0; P<DIM; P++)
310  {
311  for (unsigned Q=0; Q<DIM; Q++)
312  {
313  rDTdE(M,N,P,Q) += dTdE_coeff1 * mCurrentElementFibreDirection(M)
314  * mCurrentElementFibreDirection(N)
315  * mCurrentElementFibreDirection(P)
316  * mCurrentElementFibreDirection(Q)
317 
318  + dTdE_coeff2 * mCurrentElementFibreDirection(M)
319  * mCurrentElementFibreDirection(N)
320  * invC(P,Q);
321  if (apply_cross_fibre_tension)
322  {
323  rDTdE(M,N,P,Q) += dTdE_coeff_s1 * mCurrentElementSheetDirection(M)
324  * mCurrentElementSheetDirection(N)
325  * mCurrentElementSheetDirection(P)
326  * mCurrentElementSheetDirection(Q)
327 
328  + dTdE_coeff_s2 * mCurrentElementSheetDirection(M)
329  * mCurrentElementSheetDirection(N)
330  * invC(P,Q)
331 
332  + dTdE_coeff_s3 * mCurrentElementSheetDirection(M)
333  * mCurrentElementSheetDirection(N)
334  * mCurrentElementFibreDirection(P)
335  * mCurrentElementFibreDirection(Q);
336  if (DIM>2)
337  {
338  rDTdE(M,N,P,Q) += dTdE_coeff_n1 * mCurrentElementSheetNormalDirection(M)
339  * mCurrentElementSheetNormalDirection(N)
340  * mCurrentElementSheetNormalDirection(P)
341  * mCurrentElementSheetNormalDirection(Q)
342 
343  + dTdE_coeff_n2 * mCurrentElementSheetNormalDirection(M)
344  * mCurrentElementSheetNormalDirection(N)
345  * invC(P,Q)
346 
347  + dTdE_coeff_n3 * mCurrentElementSheetNormalDirection(M)
348  * mCurrentElementSheetNormalDirection(N)
349  * mCurrentElementFibreDirection(P)
350  * mCurrentElementFibreDirection(Q);
351  }
352  }
353  }
354  }
355  }
356  }
357  }
358 
359 // ///\todo #2180 The code below applies a cross fibre tension in the 2D case. Things that need doing:
360 // // * Refactor the common code between the block below and the block above to avoid duplication.
361 // // * Handle the 3D case.
362 // if (this->mrElectroMechanicsProblemDefinition.GetApplyCrossFibreTension() && DIM > 1)
363 // {
364 // double sheet_cross_fraction = mrElectroMechanicsProblemDefinition.GetSheetTensionFraction();
365 //
366 // for (unsigned i=0; i<DIM; i++)
367 // {
368 // mCurrentElementSheetDirection(i) = this->mChangeOfBasisMatrix(i,1);
369 // }
370 //
371 // double I4_sheet = inner_prod(mCurrentElementSheetDirection, prod(rC, mCurrentElementSheetDirection));
372 //
373 // // amend the stress and dTdE using the active tension
374 // double dTdE_coeff_s1 = -2*sheet_cross_fraction*detF*active_tension/(I4_sheet*I4_sheet); // note: I4*I4 = lam^4
375 //
376 // ///\todo #2180 The code below is specific to the implicit cardiac mechanics solver. Currently
377 // // the cross-fibre code is only tested using the explicit solver so the code below fails coverage.
378 // // This will need to be added back in once an implicit test is in place.
379 // double lambda_sheet = sqrt(I4_sheet);
380 // if (IsImplicitSolver())
381 // {
382 // double dt = mNextTime-mCurrentTime;
383 // dTdE_coeff_s1 += (d_act_tension_dlam + d_act_tension_d_dlamdt/dt)/(lambda_sheet*I4_sheet); // note: I4*lam = lam^3
384 // }
385 //
386 // rT += sheet_cross_fraction*(active_tension*detF/I4_sheet)*outer_prod(mCurrentElementSheetDirection,mCurrentElementSheetDirection);
387 //
388 // double dTdE_coeff_s2 = active_tension*detF/I4_sheet;
389 // if (addToDTdE)
390 // {
391 // for (unsigned M=0; M<DIM; M++)
392 // {
393 // for (unsigned N=0; N<DIM; N++)
394 // {
395 // for (unsigned P=0; P<DIM; P++)
396 // {
397 // for (unsigned Q=0; Q<DIM; Q++)
398 // {
399 // rDTdE(M,N,P,Q) += dTdE_coeff_s1 * mCurrentElementSheetDirection(M)
400 // * mCurrentElementSheetDirection(N)
401 // * mCurrentElementSheetDirection(P)
402 // * mCurrentElementSheetDirection(Q)
403 //
404 // + dTdE_coeff_s2 * mCurrentElementFibreDirection(M)
405 // * mCurrentElementFibreDirection(N)
406 // * invC(P,Q);
407 //
408 // }
409 // }
410 // }
411 // }
412 // }
413 // }
414 }
415 
416 
417 template<class ELASTICITY_SOLVER,unsigned DIM>
419  std::vector<c_matrix<double,DIM,DIM> >& rDeformationGradients,
420  std::vector<double>& rStretches)
421 {
422  assert(rStretches.size()==this->mrQuadMesh.GetNumElements());
423 
424  // this will only work currently if the coarse mesh fibre info is defined per element, not per quad point
425  assert(!mpVariableFibreSheetDirections || !mFibreSheetDirectionsDefinedByQuadraturePoint);
426 
427  static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> element_current_displacements;
428  static c_matrix<double,DIM,NUM_VERTICES_PER_ELEMENT> grad_lin_phi;
429  static c_matrix<double,DIM,DIM> F; // the deformation gradient, F = dx/dX, F_{iM} = dx_i/dX_M
430 
431  static c_matrix<double,DIM,DIM> jacobian;
432  static c_matrix<double,DIM,DIM> inverse_jacobian;
433  double jacobian_determinant;
434  ChastePoint<DIM> quadrature_point; // not needed, but has to be passed in
435 
436  // loop over all the elements
437  for (unsigned elem_index=0; elem_index<this->mrQuadMesh.GetNumElements(); elem_index++)
438  {
439  Element<DIM,DIM>* p_elem = this->mrQuadMesh.GetElement(elem_index);
440 
441  // get the fibre direction for this element
442  SetupChangeOfBasisMatrix(elem_index, 0); // 0 is quad index, and doesn't matter as checked that fibres not defined by quad pt above.
443  for (unsigned i=0; i<DIM; i++)
444  {
445  mCurrentElementFibreDirection(i) = this->mChangeOfBasisMatrix(i,0);
446  }
447 
448  // get the displacement at this element
449  for (unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
450  {
451  for (unsigned JJ=0; JJ<DIM; JJ++)
452  {
453  // mProblemDimension = DIM for compressible elasticity and DIM+1 for incompressible elasticity
454  element_current_displacements(JJ,II) = this->mCurrentSolution[this->mProblemDimension*p_elem->GetNodeGlobalIndex(II) + JJ];
455  }
456  }
457 
458  // set up the linear basis functions
459  this->mrQuadMesh.GetInverseJacobianForElement(elem_index, jacobian, jacobian_determinant, inverse_jacobian);
460  LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_lin_phi);
461 
462  F = identity_matrix<double>(DIM,DIM);
463 
464  // loop over the vertices and interpolate F, the deformation gradient
465  // (note: could use matrix-mult if this becomes inefficient
466  for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
467  {
468  for (unsigned i=0; i<DIM; i++)
469  {
470  for (unsigned M=0; M<DIM; M++)
471  {
472  F(i,M) += grad_lin_phi(M,node_index)*element_current_displacements(i,node_index);
473  }
474  }
475  }
476 
477  rDeformationGradients[elem_index] = F;
478 
479  // compute and save the stretch: m^T C m = ||Fm||
480  c_vector<double,DIM> deformed_fibre = prod(F, mCurrentElementFibreDirection);
481  rStretches[elem_index] = norm_2(deformed_fibre);
482  }
483 }
484 
485 
486 template<class ELASTICITY_SOLVER,unsigned DIM>
488 {
489  mFibreSheetDirectionsDefinedByQuadraturePoint = definedPerQuadraturePoint;
490 
491  FibreReader<DIM> reader(rOrthoFile, ORTHO);
492 
493  unsigned num_entries = reader.GetNumLinesOfData();
494 
495  if (!mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=this->mrQuadMesh.GetNumElements()) )
496  {
497  EXCEPTION("Number of entries defined at top of file " << rOrthoFile.GetAbsolutePath() <<
498  " does not match number of elements in the mesh, " << "found " << num_entries <<
499  ", expected " << this->mrQuadMesh.GetNumElements());
500  }
501 
502  if (mFibreSheetDirectionsDefinedByQuadraturePoint && (num_entries!=mTotalQuadPoints) )
503  {
504  EXCEPTION("Number of entries defined at top of file " << rOrthoFile.GetAbsolutePath() <<
505  " does not match number of quadrature points defined, " << "found " << num_entries <<
506  ", expected " << mTotalQuadPoints);
507  }
508 
509  mpVariableFibreSheetDirections = new std::vector<c_matrix<double,DIM,DIM> >(num_entries, zero_matrix<double>(DIM,DIM));
510  for (unsigned index=0; index<num_entries; index++)
511  {
512  reader.GetFibreSheetAndNormalMatrix(index, (*mpVariableFibreSheetDirections)[index] );
513  }
514 }
515 
516 template<class ELASTICITY_SOLVER,unsigned DIM>
518 {
519  mConstantFibreSheetDirections = rFibreSheetMatrix;
520  // check orthogonality
521  c_matrix<double,DIM,DIM> temp = prod(trans(rFibreSheetMatrix),rFibreSheetMatrix);
522  for (unsigned i=0; i<DIM; i++)
523  {
524  for (unsigned j=0; j<DIM; j++)
525  {
526  double val = (i==j ? 1.0 : 0.0);
527  if (fabs(temp(i,j)-val)>1e-4)
528  {
529  EXCEPTION("The given fibre-sheet matrix, " << rFibreSheetMatrix << ", is not orthogonal");
530  }
531  }
532  }
533 }
534 
539 
540 
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
void ComputeDeformationGradientAndStretchInEachElement(std::vector< c_matrix< double, DIM, DIM > > &rDeformationGradients, std::vector< double > &rStretches)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetNumLinesOfData()
void SetConstantFibreSheetDirections(const c_matrix< double, DIM, DIM > &rFibreSheetMatrix)
std::string GetAbsolutePath() const
Definition: FileFinder.cpp:224
AbstractCardiacMechanicsSolver(QuadraticMesh< DIM > &rQuadMesh, ElectroMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)
virtual unsigned GetNumElements() const
#define EXCEPTION(message)
Definition: Exception.hpp:143
static HeartRegionType GetValidBathId()
const AbstractTetrahedralMesh< DIM, DIM > & GetCoarseMesh() const
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM, ELEMENT_DIM+1 > &rReturnValue)
bool GetOwnership() const
void AddActiveStressAndStressDerivative(c_matrix< double, DIM, DIM > &rC, unsigned elementIndex, unsigned currentQuadPointGlobalIndex, c_matrix< double, DIM, DIM > &rT, FourthOrderTensor< DIM, DIM, DIM, DIM > &rDTdE, bool addToDTdE)
void SetFineCoarseMeshPair(FineCoarseMeshPair< DIM > *pMeshPair)
AbstractContractionModel * ContractionModel
void SetCalciumAndVoltage(std::vector< double > &rCalciumConcentrations, std::vector< double > &rVoltages)
unsigned GetIndex() const
void SetVariableFibreSheetDirections(const FileFinder &rOrthoFile, bool definedPerQuadraturePoint)
virtual AbstractContractionModel * CreateContractionCellForElement(Element< DIM, DIM > *pElement)=0
void GetFibreSheetAndNormalMatrix(unsigned fibreIndex, c_matrix< double, DIM, DIM > &rFibreMatrix, bool checkOrthogonality=true)
void SetupChangeOfBasisMatrix(unsigned elementIndex, unsigned currentQuadPointGlobalIndex)