Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
StreeterFibreGenerator.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
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"
49static PetscBool StreeterCite = PETSC_FALSE;
50const 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
61template<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
110template<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
148template<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
162template<unsigned SPACE_DIM>
167
168template<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
179template<unsigned SPACE_DIM>
181{
182 *(this->mpMasterFile) << this->mpMesh->GetNumElements();
183 *(this->mpMasterFile) << std::setprecision(16);
184}
185
186template<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
246 // 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
313template<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
446template<unsigned SPACE_DIM>
447void 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
457template<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
468template<unsigned SPACE_DIM>
470{
471 mLogInfo = logInfo;
472}
473
474// Explicit instantiation
475// LCOV_EXCL_START
476template class StreeterFibreGenerator<3>;
477// LCOV_EXCL_STOP
#define EXCEPTION(message)
#define NEVER_REACHED
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetNumNodes() const
unsigned GetIndex() const
virtual unsigned GetNumNodes() const
static void Register(const char *pCitation, PetscBool *pSet)
Definition Citations.cpp:43
std::vector< double > & rGetDistanceMapEpicardium()
Definition Node.hpp:59
ContainingElementIterator ContainingElementsEnd() const
Definition Node.hpp:493
ContainingElementIterator ContainingElementsBegin() const
Definition Node.hpp:485
unsigned GetIndex() const
Definition Node.cpp:158
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static bool AmMaster()
double GetAveragedThicknessLocalNode(const unsigned nodeIndex, const std::vector< double > &wallThickness) const
void Visit(Element< SPACE_DIM, SPACE_DIM > *pElement, unsigned localElementIndex, c_vector< double, SPACE_DIM *SPACE_DIM > &rData)
std::vector< double > mAveragedWallThickness
void SetApexToBase(const c_vector< double, SPACE_DIM > &apexToBase)
void PreWriteCalculations(OutputFileHandler &rOutputDirectory)
void SetLogInfo(bool logInfo=true)
std::vector< double > mWallThickness
double GetFibreMaxAngle(const c_vector< HeartRegionType, SPACE_DIM+1 > &nodesRegionsForElement) const
StreeterFibreGenerator(AbstractTetrahedralMesh< SPACE_DIM, SPACE_DIM > &rMesh)
void SetSurfaceFiles(const std::string &rEpicardiumFile, const std::string &rRightVentricleFile, const std::string &rLeftVentricleFile, bool indexFromZero)