StreeterFibreGenerator.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "UblasCustomFunctions.hpp"
00030 
00031 #include "StreeterFibreGenerator.hpp"
00032 
00033 #include <cmath>
00034 #include <fstream>
00035 #include <sstream>
00036 #include "OutputFileHandler.hpp"
00037 #include "Exception.hpp"
00038 #include "HeartRegionCodes.hpp"
00039 
00040 template<unsigned SPACE_DIM>
00041 double StreeterFibreGenerator<SPACE_DIM>::GetAveragedThicknessLocalNode(
00042         const unsigned nodeIndex, const std::vector<double>& wallThickness) const
00043 {
00044     if (nodeIndex < mrMesh.GetDistributedVectorFactory()->GetLow() ||
00045         nodeIndex >= mrMesh.GetDistributedVectorFactory()->GetHigh() )
00046     {
00047         return 0.0;  //Don't calculate this for nodes which aren't local
00048     }
00049 
00050     // Initialise the average with the value corresponding to the current node
00051     double average = wallThickness[nodeIndex];
00052     unsigned nodes_visited = 1;
00053 
00054     // Use a set to store visited nodes
00055     std::set<unsigned> visited_nodes;
00056     visited_nodes.insert(nodeIndex);
00057 
00058     Node<SPACE_DIM>* p_current_node = mrMesh.GetNode(nodeIndex);
00059 
00060     // Loop over the elements containing the given node
00061     for(typename Node<SPACE_DIM>::ContainingElementIterator element_iterator = p_current_node->ContainingElementsBegin();
00062         element_iterator != p_current_node->ContainingElementsEnd();
00063         ++element_iterator)
00064     {
00065         // Get a pointer to the container element
00066         Element<SPACE_DIM,SPACE_DIM>* p_containing_element = mrMesh.GetElement(*element_iterator);
00067 
00068        // Loop over the nodes of the element
00069        for(unsigned node_local_index=0;
00070            node_local_index<p_containing_element->GetNumNodes();
00071            node_local_index++)
00072        {
00073             Node<SPACE_DIM>* p_neighbour_node = p_containing_element->GetNode(node_local_index);
00074             unsigned neighbour_node_index = p_neighbour_node->GetIndex();
00075 
00076             // Check if the neighbour node has already been visited
00077             if (visited_nodes.find(neighbour_node_index) == visited_nodes.end())
00078             {
00079                 average += wallThickness[neighbour_node_index];
00080                 visited_nodes.insert(neighbour_node_index);
00081                 nodes_visited++;
00082             }
00083        }
00084     }
00085 
00086     return average/nodes_visited;
00087 }
00088 
00089 
00090 
00091 
00092 template<unsigned SPACE_DIM>
00093 double StreeterFibreGenerator<SPACE_DIM>::GetFibreMaxAngle(
00094         const c_vector<HeartRegionType, SPACE_DIM+1>& nodesRegionsForElement) const
00095 {
00096     unsigned lv=0, rv=0;
00097 
00098     for (unsigned index=0; index<SPACE_DIM+1; index++)
00099     {
00100         switch (nodesRegionsForElement[index])
00101         {
00102             case HeartRegionCode::LEFT_VENTRICLE_SURFACE:
00103             case HeartRegionCode::LEFT_VENTRICLE_WALL:
00104             case HeartRegionCode::LEFT_SEPTUM:
00105                 lv++;
00106                 break;
00107 
00108             case HeartRegionCode::RIGHT_VENTRICLE_SURFACE:
00109             case HeartRegionCode::RIGHT_VENTRICLE_WALL:
00110             case HeartRegionCode::RIGHT_SEPTUM:
00111                 rv++;
00112                 break;
00113 
00114             case HeartRegionCode::UNKNOWN:
00115             default:
00116                 NEVER_REACHED;
00117         }
00118     }
00119 
00120     // If most of the nodes are in the right ventricle
00121     if (rv>lv)
00122     {
00123         return M_PI/4.0;
00124     }
00125 
00126     // Anywhere else
00127     return M_PI/3.0;
00128 }
00129 
00130 template<unsigned SPACE_DIM>
00131 StreeterFibreGenerator<SPACE_DIM>::StreeterFibreGenerator(AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>& rMesh)
00132     : mrMesh(rMesh),
00133       mpGeometryInfo(NULL)
00134 {
00135 }
00136 
00137 template<unsigned SPACE_DIM>
00138 StreeterFibreGenerator<SPACE_DIM>::~StreeterFibreGenerator()
00139 {
00140     delete mpGeometryInfo;
00141 }
00142 
00143 template<unsigned SPACE_DIM>
00144 void StreeterFibreGenerator<SPACE_DIM>::SetSurfaceFiles(
00145             const std::string &epicardiumFile,
00146             const std::string &rightVentricleFile,
00147             const std::string &leftVentricleFile)
00148 {
00149     // Compute the distance map of each surface
00150      mpGeometryInfo = new HeartGeometryInformation<SPACE_DIM>(mrMesh, epicardiumFile, leftVentricleFile, rightVentricleFile, false);
00151 }
00152 
00153 template<unsigned SPACE_DIM>
00154 void StreeterFibreGenerator<SPACE_DIM>::GenerateOrthotropicFibreOrientation(
00155             std::string outputDirectory,
00156             std::string fibreOrientationFile,
00157             bool logInfo)
00158 {
00159     if (mpGeometryInfo == NULL)
00160     {
00161         EXCEPTION("Files defining the heart surfaces not set");
00162     }
00163     unsigned num_elements = mrMesh.GetNumElements();
00164 
00165     // Open files
00166     OutputFileHandler handler(outputDirectory, false);
00167     out_stream p_fibre_file, p_regions_file, p_thickness_file, p_ave_thickness_file;
00168 
00169     //Make sure that only the master process makes an empty output file
00170     //and writes the log files
00171     if (PetscTools::AmMaster())
00172     {
00173         //Open file and close it, but don't write to it yet
00174         p_fibre_file=handler.OpenOutputFile(fibreOrientationFile);
00175         // First line of the fibre file: number of elements of the mesh
00176         *p_fibre_file << num_elements << std::endl;
00177         p_fibre_file->close();
00178     }
00179     else
00180     {
00181         logInfo=false;
00182     }
00183 
00184     if (logInfo)
00185     {
00186         p_regions_file  = handler.OpenOutputFile("node_regions.data");
00187         p_thickness_file = handler.OpenOutputFile("wall_thickness.data");
00188         p_ave_thickness_file = handler.OpenOutputFile("averaged_thickness.data");
00189     }
00190 
00191     // Compute the distance map of each surface
00192 
00193     CheckVentricleAlignment();
00194 
00195     // Compute wall thickness parameter
00196     unsigned num_nodes = mrMesh.GetNumNodes();
00197     std::vector<double> wall_thickness(num_nodes);
00198     for (unsigned node_index=0; node_index<num_nodes; node_index++)
00199     {
00200         double dist_epi, dist_endo;
00201 
00202         HeartRegionType node_region = mpGeometryInfo->GetHeartRegion(node_index);
00203 
00204         switch(node_region)
00205         {
00206             case HeartRegionCode::LEFT_VENTRICLE_SURFACE:
00207             case HeartRegionCode::LEFT_VENTRICLE_WALL:
00208                 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
00209                 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00210                 break;
00211 
00212             case HeartRegionCode::RIGHT_VENTRICLE_SURFACE:
00213             case HeartRegionCode::RIGHT_VENTRICLE_WALL:
00214                 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
00215                 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00216                 break;
00217 
00218             case HeartRegionCode::LEFT_SEPTUM:
00219                 dist_epi = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00220                 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00221                 break;
00222 
00223             case HeartRegionCode::RIGHT_SEPTUM:
00224                 dist_epi = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00225                 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00226                 break;
00227 
00228             case HeartRegionCode::UNKNOWN:
00229                 #define COVERAGE_IGNORE
00230                 std::cerr << "Wrong distances node: " << node_index << "\t"
00231                           << "Epi " << mpGeometryInfo->rGetDistanceMapEpicardium()[node_index] << "\t"
00232                           << "RV " << mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index] << "\t"
00233                           << "LV " << mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index]
00234                           << std::endl;
00235 
00236                 // Make wall_thickness=0 as in Martin's code
00237                 dist_epi = 1;
00238                 dist_endo = 0;
00239                 break;
00240                 #undef COVERAGE_IGNORE
00241 
00242             default:
00243                 NEVER_REACHED;
00244         }
00245 
00246         wall_thickness[node_index] = dist_endo / (dist_endo + dist_epi);
00247 
00248         if (std::isnan(wall_thickness[node_index]))
00249         {
00250             #define COVERAGE_IGNORE
00251             /*
00252              *  A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
00253              *  By setting its value to 0 we consider it contained only on the lv (or rv) surface.
00254              */
00255             wall_thickness[node_index] = 0;
00256             #undef COVERAGE_IGNORE
00257         }
00258 
00259         if (logInfo)
00260         {
00261             *p_regions_file << node_region*100 << "\n";
00262             *p_thickness_file << wall_thickness[node_index] << "\n";
00263         }
00264     }
00265 
00266     /*
00267      *  For each node, average its value of e with the values of all the neighbours
00268      */
00269     std::vector<double> my_averaged_wall_thickness(num_nodes);
00270     std::vector<double> averaged_wall_thickness(num_nodes);
00271     for (unsigned node_index=0; node_index<num_nodes; node_index++)
00272     {
00273         my_averaged_wall_thickness[node_index] = GetAveragedThicknessLocalNode(node_index, wall_thickness);
00274     }
00275 
00276     //Non-local information appear as zeros in the vector
00277     MPI_Allreduce(&my_averaged_wall_thickness[0], &averaged_wall_thickness[0], num_nodes,
00278                   MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00279     if (logInfo)
00280     {
00281         for (unsigned node_index=0; node_index<num_nodes; node_index++)
00282         {
00283              *p_ave_thickness_file << averaged_wall_thickness[node_index] << "\n";
00284         }
00285 
00286     }
00287 
00288     /*
00289      *  Compute the gradient of the averaged wall thickness at the centroid of each tetrahedral
00290      */
00291     c_vector<double,SPACE_DIM> grad_ave_wall_thickness;
00292 
00293     bool fibre_file_is_open= false;
00294     for (unsigned element_index=0; element_index<num_elements; element_index++)
00295     {
00296         //PRINT_VARIABLE(element_index);
00297         if (mrMesh.CalculateDesignatedOwnershipOfElement(element_index))
00298         {
00299             //PRINT_VARIABLE(element_index);
00300             //Locally own this element
00301             //Make sure that writing will proceed
00302             //The same file is written to by multiple processes, therefore we need
00303             //this barrier - Barrier happens BEFORE opening and AFTER closing
00304 
00305             PetscTools::Barrier("GenerateOrthotropicFibreOrientation");
00306             if (fibre_file_is_open != true)
00307             {
00308                 //Open and append
00309                 p_fibre_file = handler.OpenOutputFile(fibreOrientationFile, std::ios::app);
00310                 fibre_file_is_open=true;
00311             }
00312 
00313 
00314             Element<SPACE_DIM,SPACE_DIM>* p_element = mrMesh.GetElement(element_index);
00315 
00316             /*
00317              *  The gradient of the averaged thickness at the element is:
00318              *
00319              *     grad_ave_wall_thickness[element_index] = ave' * BF * inv(J)
00320              *
00321              *  being : ave, averaged thickness values of the nodes defining the element
00322              *          J,   the Jacobian of the element as defined in class Element.
00323              *                               (-1 -1 -1)
00324              *          BF,  basis functions ( 1  0  0)
00325              *                               ( 0  1  0)
00326              *                               ( 0  0  1)
00327              *
00328              *  Defined as u in Streeter paper
00329              */
00330 
00331             c_vector<double, SPACE_DIM+1> elem_nodes_ave_thickness;
00332             double element_averaged_thickness = 0.0;
00333             c_vector<HeartRegionType, SPACE_DIM+1> elem_nodes_region;
00334 
00335             for (unsigned local_node_index=0; local_node_index<SPACE_DIM+1; local_node_index++)
00336             {
00337                 // Get node's global index
00338                 unsigned global_node_index = p_element->GetNode(local_node_index)->GetIndex();
00339 
00340                 elem_nodes_ave_thickness[local_node_index] = averaged_wall_thickness[global_node_index];
00341                 elem_nodes_region[local_node_index] = mpGeometryInfo->GetHeartRegion(global_node_index);
00342 
00343                 // Calculate wall thickness averaged value for the element
00344                 element_averaged_thickness +=  wall_thickness[global_node_index];
00345             }
00346 
00347             element_averaged_thickness /= SPACE_DIM+1;
00348 
00349             c_matrix<double, SPACE_DIM+1, SPACE_DIM> basis_functions( zero_matrix<double>(4u,3u) );
00350             basis_functions(0,0) = basis_functions(0,1) = basis_functions(0,2) = -1.0;
00351             basis_functions(1,0) = basis_functions(2,1) = basis_functions(3,2) =  1.0;
00352 
00353             c_matrix<double, SPACE_DIM+1, SPACE_DIM> temp;
00354             c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian, inverse_jacobian;
00355             double jacobian_det;
00356             mrMesh.GetInverseJacobianForElement(element_index, jacobian, jacobian_det, inverse_jacobian);
00357             noalias(temp) = prod (basis_functions, inverse_jacobian);
00358             noalias(grad_ave_wall_thickness) = prod(elem_nodes_ave_thickness, temp);
00359 
00360             grad_ave_wall_thickness /= norm_2(grad_ave_wall_thickness);
00361 
00362             /*
00363              *   Normal to the gradient (v in Streeter paper) which is then the circumferential direction
00364              * (it will be the fibre direction after rotation)
00365              *
00366              *  Computed as the cross product with the x-axis (assuming base-apex axis is x). The output vector is not normal,
00367              * since the angle between them may be != 90, normalise it.
00368              */
00369             c_vector<double, SPACE_DIM> fibre_direction = VectorProduct(grad_ave_wall_thickness, Create_c_vector(1.0, 0.0, 0.0));
00370             fibre_direction /= norm_2(fibre_direction);
00371 
00372             /*
00373              *  Longitude direction (w in Streeter paper)
00374              */
00375             c_vector<double, SPACE_DIM> longitude_direction = VectorProduct(grad_ave_wall_thickness, fibre_direction);
00376 
00377             /*
00378              *  Compute fibre to v angle: alpha = R*(1-2e)^3
00379              *
00380              *    R is the maximum angle between the fibre and the v axis (heart region dependant)
00381              *    (1 - 2e)^3 scales it by a value in [-1, 1] defining the rotation of the fibre based
00382              *       on the position in the wall
00383              */
00384             double alpha = GetFibreMaxAngle(elem_nodes_region) * SmallPow( (1 - 2*element_averaged_thickness), 3 );
00385 
00386             /*
00387              *  Apply alpha rotation about the u axis to the orthonormal basis
00388              *
00389              *               ( u(1) v(1) w(1) )
00390              *   (u, v, w) = ( u(2) v(2) w(2) )
00391              *               ( u(3) v(3) w(3) )
00392              *
00393              *  The following matrix defines a rotation about the u axis
00394              *
00395              *                 ( 1        0           0      ) (u')
00396              *   R = (u, v, w) ( 0    cos(alpha) -sin(alpha) ) (v')
00397              *                 ( 0    sin(alpha)  cos(alpha) ) (w')
00398              *
00399              *  The rotated basis is computed like:
00400              *
00401              *                                             ( 1        0           0      )
00402              *  (u, v_r, w_r ) = R * (u, v, w) = (u, v, w) ( 0    cos(alpha) -sin(alpha) )
00403              *                                             ( 0    sin(alpha)  cos(alpha) )
00404              *
00405              *  Which simplifies to:
00406              *
00407              *   v_r =  v*cos(alpha) + w*sin(alpha)
00408              *   w_r = -v*sin(alpha) + w*cos(alpha)
00409              */
00410             c_vector<double, SPACE_DIM> rotated_fibre_direction = fibre_direction*cos(alpha) + longitude_direction*sin(alpha);
00411             c_vector<double, SPACE_DIM> rotated_longitude_direction = -fibre_direction*sin(alpha) + longitude_direction*cos(alpha);
00412 
00413 
00414             /*
00415              * Test the orthonormality of the basis
00416              */
00417             assert( fabs(norm_2(rotated_fibre_direction) - 1) < 100*DBL_EPSILON );
00418             assert( fabs(norm_2(grad_ave_wall_thickness) - 1) < 100*DBL_EPSILON );
00419             assert( fabs(norm_2(rotated_longitude_direction) - 1) < 100*DBL_EPSILON );
00420 
00421             assert( fabs(inner_prod(rotated_fibre_direction, grad_ave_wall_thickness)) < 100*DBL_EPSILON );
00422             assert( fabs(inner_prod(rotated_fibre_direction, rotated_longitude_direction)) < 100*DBL_EPSILON);
00423             assert( fabs(inner_prod(grad_ave_wall_thickness, rotated_longitude_direction)) < 100*DBL_EPSILON);
00424 
00425             /*
00426              *  Output the direction of the myofibre, the transverse to it in the plane of the myocite laminae and the normal to this laminae (in that order)
00427              *
00428              *  See Fig. 1 "Laminar Structure of the Heart: a mathematical model" LeGrice et al. 97
00429              *
00430              */
00431             *p_fibre_file << rotated_fibre_direction[0]     << " " << rotated_fibre_direction[1]     << " "  << rotated_fibre_direction[2]     << " "
00432                           << grad_ave_wall_thickness[0]     << " " << grad_ave_wall_thickness[1]     << " "  << grad_ave_wall_thickness[2]     << " "
00433                           << rotated_longitude_direction[0] << " " << rotated_longitude_direction[1] << " "  << rotated_longitude_direction[2] << std::endl;
00434         }
00435         else
00436         {
00437             //We don't locally own this element
00438             //Give up the file pointer
00439             if (fibre_file_is_open)
00440             {
00441                 p_fibre_file->close();
00442                 fibre_file_is_open=false;
00443             }
00444             PetscTools::Barrier("GenerateOrthotropicFibreOrientation");
00445 
00446         }
00447     }
00448 
00449     if (fibre_file_is_open)
00450     {
00451         p_fibre_file->close();
00452     }
00453 
00454     if (logInfo)
00455     {
00456         p_regions_file->close();
00457         p_thickness_file->close();
00458         p_ave_thickness_file->close();
00459     }
00460 
00461 }
00462 
00463 template<unsigned SPACE_DIM>
00464 void StreeterFibreGenerator<SPACE_DIM>::CheckVentricleAlignment()
00465 {
00466     assert(SPACE_DIM == 3);
00467 
00468     ChasteCuboid<SPACE_DIM> lv_bounds=mpGeometryInfo->CalculateBoundingBoxOfLV();
00469     ChasteCuboid<SPACE_DIM> rv_bounds=mpGeometryInfo->CalculateBoundingBoxOfRV();
00470 
00471     //Check that LV midway point is not inside the RV interval
00472     double lv_y_midway=lv_bounds.rGetUpperCorner()[1] + lv_bounds.rGetLowerCorner()[1];
00473     lv_y_midway /= 2.0;
00474     assert(lv_y_midway < rv_bounds.rGetLowerCorner()[1]  || lv_y_midway > rv_bounds.rGetUpperCorner()[1]);
00475 
00476     //Check that RV midway point is not inside the LV interval
00477     double rv_y_midway=rv_bounds.rGetUpperCorner()[1] + rv_bounds.rGetLowerCorner()[1];
00478     rv_y_midway /= 2.0;
00479     assert(rv_y_midway < lv_bounds.rGetLowerCorner()[1]  || rv_y_midway > lv_bounds.rGetUpperCorner()[1]);
00480 }
00481 
00483 // Explicit instantiation
00485 #define COVERAGE_IGNORE
00486 template class StreeterFibreGenerator<3>;
00487 #undef COVERAGE_IGNORE

Generated by  doxygen 1.6.2