Chaste  Release::2018.1
StreeterFibreGenerator.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 "UblasCustomFunctions.hpp"
37 
38 #include "StreeterFibreGenerator.hpp"
39 
40 #include <cmath>
41 #include <fstream>
42 #include <sstream>
43 #include "OutputFileHandler.hpp"
44 #include "Exception.hpp"
45 //#include "HeartRegionCodes.hpp"
46 
47 // Add the citation for original Streeter paper.
48 #include "Citations.hpp"
49 static PetscBool StreeterCite = PETSC_FALSE;
50 const char StreeterCitation[] = "@article{streeter1969fiber,\n"
51 " title={Fiber orientation in the canine left ventricle during diastole and systole},\n"
52 " author={Streeter, Daniel D and Spotnitz, Henry M and Patel, Dali P and Ross, John and Sonnenblick, Edmund H},\n"
53 " journal={Circulation research},\n"
54 " volume={24},\n"
55 " number={3},\n"
56 " pages={339--347},\n"
57 " year={1969},\n"
58 " publisher={Am Heart Assoc}\n"
59 "}\n";
60 
61 template<unsigned SPACE_DIM>
63  const unsigned nodeIndex, const std::vector<double>& wallThickness) const
64 {
65  if (nodeIndex < this->mpMesh->GetDistributedVectorFactory()->GetLow() ||
66  nodeIndex >= this->mpMesh->GetDistributedVectorFactory()->GetHigh() )
67  {
68  return 0.0; //Don't calculate this for nodes which aren't local
69  }
70 
71  // Initialise the average with the value corresponding to the current node
72  double average = wallThickness[nodeIndex];
73  unsigned nodes_visited = 1;
74 
75  // Use a set to store visited nodes
76  std::set<unsigned> visited_nodes;
77  visited_nodes.insert(nodeIndex);
78 
79  Node<SPACE_DIM>* p_current_node = this->mpMesh->GetNode(nodeIndex);
80 
81  // Loop over the elements containing the given node
82  for (typename Node<SPACE_DIM>::ContainingElementIterator element_iterator = p_current_node->ContainingElementsBegin();
83  element_iterator != p_current_node->ContainingElementsEnd();
84  ++element_iterator)
85  {
86  // Get a pointer to the container element
87  Element<SPACE_DIM,SPACE_DIM>* p_containing_element = this->mpMesh->GetElement(*element_iterator);
88 
89  // Loop over the nodes of the element
90  for (unsigned node_local_index=0;
91  node_local_index<p_containing_element->GetNumNodes();
92  node_local_index++)
93  {
94  Node<SPACE_DIM>* p_neighbour_node = p_containing_element->GetNode(node_local_index);
95  unsigned neighbour_node_index = p_neighbour_node->GetIndex();
96 
97  // Check if the neighbour node has already been visited
98  if (visited_nodes.find(neighbour_node_index) == visited_nodes.end())
99  {
100  average += wallThickness[neighbour_node_index];
101  visited_nodes.insert(neighbour_node_index);
102  nodes_visited++;
103  }
104  }
105  }
106 
107  return average/nodes_visited;
108 }
109 
110 template<unsigned SPACE_DIM>
112  const c_vector<HeartRegionType, SPACE_DIM+1>& nodesRegionsForElement) const
113 {
114  unsigned lv=0, rv=0;
115 
116  for (unsigned index=0; index<SPACE_DIM+1; index++)
117  {
118  switch (nodesRegionsForElement[index])
119  {
123  lv++;
124  break;
125 
129  rv++;
130  break;
131 
133  default:
135  }
136  }
137 
138  // If most of the nodes are in the right ventricle
139  if (rv>lv)
140  {
141  return M_PI/4.0;
142  }
143 
144  // Anywhere else
145  return M_PI/3.0;
146 }
147 
148 template<unsigned SPACE_DIM>
150  : AbstractPerElementWriter<SPACE_DIM,SPACE_DIM,SPACE_DIM*SPACE_DIM>(&rMesh),
151  mpGeometryInfo(NULL),
152  mApexToBase(zero_vector<double>(SPACE_DIM)),
153  mLogInfo(false)
154 {
155  // Record a reference for the calculations performed here, can be extracted with the '-citations' flag.
156  Citations::Register(StreeterCitation, &StreeterCite);
157 
158  mWallThickness.resize(rMesh.GetNumNodes());
159  mAveragedWallThickness.resize(rMesh.GetNumNodes());
160 }
161 
162 template<unsigned SPACE_DIM>
164 {
165  delete mpGeometryInfo;
166 }
167 
168 template<unsigned SPACE_DIM>
170  const std::string& rEpicardiumFile,
171  const std::string& rRightVentricleFile,
172  const std::string& rLeftVentricleFile,
173  bool indexFromZero)
174 {
175  // Compute the distance map of each surface
176  mpGeometryInfo = new HeartGeometryInformation<SPACE_DIM>(*(this->mpMesh), rEpicardiumFile, rLeftVentricleFile, rRightVentricleFile, indexFromZero);
177 }
178 
179 template<unsigned SPACE_DIM>
181 {
182  *(this->mpMasterFile) << this->mpMesh->GetNumElements();
183  *(this->mpMasterFile) << std::setprecision(16);
184 }
185 
186 template<unsigned SPACE_DIM>
188 {
189  assert(SPACE_DIM == 3);
190  if (mpGeometryInfo == NULL)
191  {
192  EXCEPTION("Files defining the heart surfaces not set");
193  }
194 
195  // Open files
196  out_stream p_regions_file, p_thickness_file, p_ave_thickness_file;
197 
198  //Make sure that only the master process writes the log files if requested.
199  bool logInfo = PetscTools::AmMaster() && mLogInfo;
200 
201  if (logInfo)
202  {
203  p_regions_file = rOutputDirectory.OpenOutputFile("node_regions.data");
204  p_thickness_file = rOutputDirectory.OpenOutputFile("wall_thickness.data");
205  p_ave_thickness_file = rOutputDirectory.OpenOutputFile("averaged_thickness.data");
206  }
207 
208  //We expect that the apex to base has been set
209  if (fabs(norm_2(mApexToBase)) < DBL_EPSILON)
210  {
211  EXCEPTION("Apex to base vector has not been set");
212  }
213 
214  // Compute wall thickness parameter
215  unsigned num_nodes = this->mpMesh->GetNumNodes();
216  for (unsigned node_index=0; node_index<num_nodes; node_index++)
217  {
218  double dist_epi, dist_endo;
219 
220  HeartRegionType node_region = mpGeometryInfo->GetHeartRegion(node_index);
221 
222  switch(node_region)
223  {
226  dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
227  dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
228  break;
229 
232  dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
233  dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
234  break;
235 
237  dist_epi = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
238  dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
239  break;
240 
242  dist_epi = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
243  dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
244  break;
245 
247  // LCOV_EXCL_START
248  std::cerr << "Wrong distances node: " << node_index << "\t"
249  << "Epi " << mpGeometryInfo->rGetDistanceMapEpicardium()[node_index] << "\t"
250  << "RV " << mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index] << "\t"
251  << "LV " << mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index]
252  << std::endl;
253 
254  // Make wall_thickness=0 as in Martin's code
255  dist_epi = 1;
256  dist_endo = 0;
257  break;
258  // LCOV_EXCL_STOP
259 
260  default:
262  }
263 
264  mWallThickness[node_index] = dist_endo / (dist_endo + dist_epi);
265 
266  if (std::isnan(mWallThickness[node_index]))
267  {
268  // LCOV_EXCL_START
269  /*
270  * A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
271  * By setting its value to 0 we consider it contained only on the lv (or rv) surface.
272  */
273  mWallThickness[node_index] = 0;
274  // LCOV_EXCL_STOP
275  }
276 
277  if (logInfo)
278  {
279  *p_regions_file << node_region*100 << "\n";
280  *p_thickness_file << mWallThickness[node_index] << "\n";
281  }
282  }
283 
284  /*
285  * For each node, average its value of e with the values of all the neighbours
286  */
287  std::vector<double> my_averaged_wall_thickness(num_nodes);
288  for (unsigned node_index=0; node_index<num_nodes; node_index++)
289  {
290  my_averaged_wall_thickness[node_index] = GetAveragedThicknessLocalNode(node_index, mWallThickness);
291  }
292 
293  // Non-local information appear as zeros in the vector
294  MPI_Allreduce(&my_averaged_wall_thickness[0], &mAveragedWallThickness[0], num_nodes,
295  MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
296 
297  if (logInfo)
298  {
299  for (unsigned node_index=0; node_index<num_nodes; node_index++)
300  {
301  *p_ave_thickness_file << mAveragedWallThickness[node_index] << "\n";
302  }
303  }
304 
305  if (logInfo)
306  {
307  p_regions_file->close();
308  p_thickness_file->close();
309  p_ave_thickness_file->close();
310  }
311 }
312 
313 template<unsigned SPACE_DIM>
315  unsigned localElementIndex,
316  c_vector<double, SPACE_DIM*SPACE_DIM>& rData)
317 {
318  /*
319  * The gradient of the averaged thickness at the element is:
320  *
321  * grad_ave_wall_thickness[element_index] = ave' * BF * inv(J)
322  *
323  * being : ave, averaged thickness values of the nodes defining the element
324  * J, the Jacobian of the element as defined in class Element.
325  * (-1 -1 -1)
326  * BF, basis functions ( 1 0 0)
327  * ( 0 1 0)
328  * ( 0 0 1)
329  *
330  * Defined as u in Streeter paper.
331  */
332  c_vector<double, SPACE_DIM> grad_ave_wall_thickness;
333  c_vector<double, SPACE_DIM+1> elem_nodes_ave_thickness;
334  double element_averaged_thickness = 0.0;
335  c_vector<HeartRegionType, SPACE_DIM+1> elem_nodes_region;
336 
337  for (unsigned local_node_index=0; local_node_index<SPACE_DIM+1; local_node_index++)
338  {
339  // Get node's global index
340  unsigned global_node_index = pElement->GetNode(local_node_index)->GetIndex();
341 
342  elem_nodes_ave_thickness[local_node_index] = mAveragedWallThickness[global_node_index];
343  elem_nodes_region[local_node_index] = mpGeometryInfo->GetHeartRegion(global_node_index);
344 
345  // Calculate wall thickness averaged value for the element
346  element_averaged_thickness += mWallThickness[global_node_index];
347  }
348 
349  element_averaged_thickness /= SPACE_DIM+1;
350 
351  c_matrix<double, SPACE_DIM+1, SPACE_DIM> basis_functions( zero_matrix<double>(4u,3u) );
352  basis_functions(0,0) = basis_functions(0,1) = basis_functions(0,2) = -1.0;
353  basis_functions(1,0) = basis_functions(2,1) = basis_functions(3,2) = 1.0;
354 
355  c_matrix<double, SPACE_DIM+1, SPACE_DIM> temp;
356  c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian, inverse_jacobian;
357  double jacobian_det;
358  unsigned element_index = pElement->GetIndex();
359  this->mpMesh->GetInverseJacobianForElement(element_index, jacobian, jacobian_det, inverse_jacobian);
360  noalias(temp) = prod (basis_functions, inverse_jacobian);
361  noalias(grad_ave_wall_thickness) = prod(elem_nodes_ave_thickness, temp);
362 
363  grad_ave_wall_thickness /= norm_2(grad_ave_wall_thickness);
364 
365  /*
366  * Normal to the gradient (v in Streeter paper) which is then the circumferential direction
367  * (it will be the fibre direction after rotation)
368  *
369  * Computed as the cross product with the base-apex direction (originally assumed base-apex axis is x). The output vector is not normal,
370  * since the angle between them may be != 90, normalise it.
371  */
372  c_vector<double, SPACE_DIM> fibre_direction = VectorProduct(grad_ave_wall_thickness, mApexToBase);
373  fibre_direction /= norm_2(fibre_direction);
374 
375  /*
376  * Longitude direction (w in Streeter paper)
377  */
378  c_vector<double, SPACE_DIM> longitude_direction = VectorProduct(grad_ave_wall_thickness, fibre_direction);
379 
380  /*
381  * Compute fibre to v angle: alpha = R*(1-2e)^3
382  *
383  * R is the maximum angle between the fibre and the v axis (heart region dependant)
384  * (1 - 2e)^3 scales it by a value in [-1, 1] defining the rotation of the fibre based
385  * on the position in the wall
386  */
387  double alpha = GetFibreMaxAngle(elem_nodes_region) * SmallPow( (1 - 2*element_averaged_thickness), 3 );
388 
389  /*
390  * Apply alpha rotation about the u axis to the orthonormal basis
391  *
392  * ( u(1) v(1) w(1) )
393  * (u, v, w) = ( u(2) v(2) w(2) )
394  * ( u(3) v(3) w(3) )
395  *
396  * The following matrix defines a rotation about the u axis
397  *
398  * ( 1 0 0 ) (u')
399  * R = (u, v, w) ( 0 cos(alpha) -sin(alpha) ) (v')
400  * ( 0 sin(alpha) cos(alpha) ) (w')
401  *
402  * The rotated basis is computed like:
403  *
404  * ( 1 0 0 )
405  * (u, v_r, w_r ) = R * (u, v, w) = (u, v, w) ( 0 cos(alpha) -sin(alpha) )
406  * ( 0 sin(alpha) cos(alpha) )
407  *
408  * Which simplifies to:
409  *
410  * v_r = v*cos(alpha) + w*sin(alpha)
411  * w_r = -v*sin(alpha) + w*cos(alpha)
412  */
413  c_vector<double, SPACE_DIM> rotated_fibre_direction = fibre_direction*cos(alpha) + longitude_direction*sin(alpha);
414  c_vector<double, SPACE_DIM> rotated_longitude_direction = -fibre_direction*sin(alpha) + longitude_direction*cos(alpha);
415 
416 
417  /*
418  * Test the orthonormality of the basis
419  */
420  assert( fabs(norm_2(rotated_fibre_direction) - 1) < 100*DBL_EPSILON );
421  assert( fabs(norm_2(grad_ave_wall_thickness) - 1) < 100*DBL_EPSILON );
422  assert( fabs(norm_2(rotated_longitude_direction) - 1) < 100*DBL_EPSILON );
423 
424  assert( fabs(inner_prod(rotated_fibre_direction, grad_ave_wall_thickness)) < 100*DBL_EPSILON );
425  assert( fabs(inner_prod(rotated_fibre_direction, rotated_longitude_direction)) < 100*DBL_EPSILON);
426  assert( fabs(inner_prod(grad_ave_wall_thickness, rotated_longitude_direction)) < 100*DBL_EPSILON);
427 
428  /*
429  * Output the direction of the myofibre, the transverse to it in the plane
430  * of the myocite laminae and the normal to this laminae (in that order)
431  *
432  * See Fig. 1 "Laminar Structure of the Heart: a mathematical model" LeGrice et al. 97
433  *
434  */
435  rData[0] = rotated_fibre_direction[0];
436  rData[1] = rotated_fibre_direction[1];
437  rData[2] = rotated_fibre_direction[2];
438  rData[3] = grad_ave_wall_thickness[0];
439  rData[4] = grad_ave_wall_thickness[1];
440  rData[5] = grad_ave_wall_thickness[2];
441  rData[6] = rotated_longitude_direction[0];
442  rData[7] = rotated_longitude_direction[1];
443  rData[8] = rotated_longitude_direction[2];
444 }
445 
446 template<unsigned SPACE_DIM>
447 void StreeterFibreGenerator<SPACE_DIM>::SetApexToBase(const c_vector<double, SPACE_DIM>& apexToBase)
448 {
449  double norm = norm_2(apexToBase);
450  if (norm < DBL_EPSILON)
451  {
452  EXCEPTION("Apex to base vector should be non-zero");
453  }
454  mApexToBase = apexToBase / norm;
455 }
456 
457 template<unsigned SPACE_DIM>
459 {
460  if (axis >= SPACE_DIM)
461  {
462  EXCEPTION("Apex to base coordinate axis was out of range");
463  }
464  mApexToBase = zero_vector<double>(SPACE_DIM);
465  mApexToBase[axis] = 1.0;
466 }
467 
468 template<unsigned SPACE_DIM>
470 {
471  mLogInfo = logInfo;
472 }
473 
474 // Explicit instantiation
475 // LCOV_EXCL_START
476 template class StreeterFibreGenerator<3>;
477 // LCOV_EXCL_STOP
void SetLogInfo(bool logInfo=true)
Definition: Node.hpp:58
void SetApexToBase(const c_vector< double, SPACE_DIM > &apexToBase)
double SmallPow(double x, unsigned exponent)
#define EXCEPTION(message)
Definition: Exception.hpp:143
static bool AmMaster()
Definition: PetscTools.cpp:120
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
virtual unsigned GetNumNodes() const
std::vector< double > mAveragedWallThickness
#define NEVER_REACHED
Definition: Exception.hpp:206
std::vector< double > & rGetDistanceMapEpicardium()
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
StreeterFibreGenerator(AbstractTetrahedralMesh< SPACE_DIM, SPACE_DIM > &rMesh)
ContainingElementIterator ContainingElementsBegin() const
Definition: Node.hpp:485
double GetAveragedThicknessLocalNode(const unsigned nodeIndex, const std::vector< double > &wallThickness) const
unsigned GetNumNodes() const
void SetSurfaceFiles(const std::string &rEpicardiumFile, const std::string &rRightVentricleFile, const std::string &rLeftVentricleFile, bool indexFromZero)
double GetFibreMaxAngle(const c_vector< HeartRegionType, SPACE_DIM+1 > &nodesRegionsForElement) const
void PreWriteCalculations(OutputFileHandler &rOutputDirectory)
static void Register(const char *pCitation, PetscBool *pSet)
Definition: Citations.cpp:43
std::vector< double > mWallThickness
unsigned GetIndex() const
unsigned GetIndex() const
Definition: Node.cpp:158
ContainingElementIterator ContainingElementsEnd() const
Definition: Node.hpp:493
void Visit(Element< SPACE_DIM, SPACE_DIM > *pElement, unsigned localElementIndex, c_vector< double, SPACE_DIM *SPACE_DIM > &rData)