StreeterFibreGenerator.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 HeartGeometryInformation<SPACE_DIM>::LEFT_VENTRICLE_SURFACE:
00103             case HeartGeometryInformation<SPACE_DIM>::LEFT_VENTRICLE_WALL:
00104             case HeartGeometryInformation<SPACE_DIM>::LEFT_SEPTUM:
00105                 lv++;
00106                 break;
00107 
00108             case HeartGeometryInformation<SPACE_DIM>::RIGHT_VENTRICLE_SURFACE:
00109             case HeartGeometryInformation<SPACE_DIM>::RIGHT_VENTRICLE_WALL:
00110             case HeartGeometryInformation<SPACE_DIM>::RIGHT_SEPTUM:
00111                 rv++;
00112                 break;
00113 
00114             case HeartGeometryInformation<SPACE_DIM>::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       mApexToBase(zero_vector<double>(SPACE_DIM))
00135 {
00136 }
00137 
00138 template<unsigned SPACE_DIM>
00139 StreeterFibreGenerator<SPACE_DIM>::~StreeterFibreGenerator()
00140 {
00141     delete mpGeometryInfo;
00142 }
00143 
00144 template<unsigned SPACE_DIM>
00145 void StreeterFibreGenerator<SPACE_DIM>::SetSurfaceFiles(
00146             const std::string &epicardiumFile,
00147             const std::string &rightVentricleFile,
00148             const std::string &leftVentricleFile,
00149             bool indexFromZero)
00150 {
00151     // Compute the distance map of each surface
00152      mpGeometryInfo = new HeartGeometryInformation<SPACE_DIM>(mrMesh, epicardiumFile, leftVentricleFile, rightVentricleFile, indexFromZero);
00153 }
00154 
00155 template<unsigned SPACE_DIM>
00156 void StreeterFibreGenerator<SPACE_DIM>::GenerateOrthotropicFibreOrientation(
00157             std::string outputDirectory,
00158             std::string fibreOrientationFile,
00159             bool logInfo)
00160 {
00161     assert(SPACE_DIM == 3);
00162     if (mpGeometryInfo == NULL)
00163     {
00164         EXCEPTION("Files defining the heart surfaces not set");
00165     }
00166     unsigned num_elements = mrMesh.GetNumElements();
00167 
00168     // Open files
00169     OutputFileHandler handler(outputDirectory, false);
00170     out_stream p_fibre_file, p_regions_file, p_thickness_file, p_ave_thickness_file;
00171 
00172     //Make sure that only the master process makes an empty output file
00173     //and writes the log files
00174     if (PetscTools::AmMaster())
00175     {
00176         //Open file and close it, but don't write to it yet
00177         p_fibre_file=handler.OpenOutputFile(fibreOrientationFile);
00178         // First line of the fibre file: number of elements of the mesh
00179         *p_fibre_file << num_elements << std::endl;
00180         p_fibre_file->close();
00181     }
00182     else
00183     {
00184         logInfo=false;
00185     }
00186 
00187     if (logInfo)
00188     {
00189         p_regions_file  = handler.OpenOutputFile("node_regions.data");
00190         p_thickness_file = handler.OpenOutputFile("wall_thickness.data");
00191         p_ave_thickness_file = handler.OpenOutputFile("averaged_thickness.data");
00192     }
00193 
00194     
00195     //We expect that the apex to base has been set
00196     if (fabs(norm_2(mApexToBase)) < DBL_EPSILON)
00197     {
00198         EXCEPTION("Apex to base vector has not been set");
00199     }
00200 
00201     // Compute wall thickness parameter
00202     unsigned num_nodes = mrMesh.GetNumNodes();
00203     std::vector<double> wall_thickness(num_nodes);
00204     for (unsigned node_index=0; node_index<num_nodes; node_index++)
00205     {
00206         double dist_epi, dist_endo;
00207 
00208         HeartRegionType node_region = mpGeometryInfo->GetHeartRegion(node_index);
00209 
00210         switch(node_region)
00211         {
00212             case HeartGeometryInformation<SPACE_DIM>::LEFT_VENTRICLE_SURFACE:
00213             case HeartGeometryInformation<SPACE_DIM>::LEFT_VENTRICLE_WALL:
00214                 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
00215                 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00216                 break;
00217 
00218             case HeartGeometryInformation<SPACE_DIM>::RIGHT_VENTRICLE_SURFACE:
00219             case HeartGeometryInformation<SPACE_DIM>::RIGHT_VENTRICLE_WALL:
00220                 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
00221                 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00222                 break;
00223 
00224             case HeartGeometryInformation<SPACE_DIM>::LEFT_SEPTUM:
00225                 dist_epi = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00226                 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00227                 break;
00228 
00229             case HeartGeometryInformation<SPACE_DIM>::RIGHT_SEPTUM:
00230                 dist_epi = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
00231                 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
00232                 break;
00233 
00234             case HeartGeometryInformation<SPACE_DIM>::UNKNOWN:
00235                 #define COVERAGE_IGNORE
00236                 std::cerr << "Wrong distances node: " << node_index << "\t"
00237                           << "Epi " << mpGeometryInfo->rGetDistanceMapEpicardium()[node_index] << "\t"
00238                           << "RV " << mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index] << "\t"
00239                           << "LV " << mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index]
00240                           << std::endl;
00241 
00242                 // Make wall_thickness=0 as in Martin's code
00243                 dist_epi = 1;
00244                 dist_endo = 0;
00245                 break;
00246                 #undef COVERAGE_IGNORE
00247 
00248             default:
00249                 NEVER_REACHED;
00250         }
00251 
00252         wall_thickness[node_index] = dist_endo / (dist_endo + dist_epi);
00253 
00254         if (std::isnan(wall_thickness[node_index]))
00255         {
00256             #define COVERAGE_IGNORE
00257             /*
00258              *  A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
00259              *  By setting its value to 0 we consider it contained only on the lv (or rv) surface.
00260              */
00261             wall_thickness[node_index] = 0;
00262             #undef COVERAGE_IGNORE
00263         }
00264 
00265         if (logInfo)
00266         {
00267             *p_regions_file << node_region*100 << "\n";
00268             *p_thickness_file << wall_thickness[node_index] << "\n";
00269         }
00270     }
00271 
00272     /*
00273      *  For each node, average its value of e with the values of all the neighbours
00274      */
00275     std::vector<double> my_averaged_wall_thickness(num_nodes);
00276     std::vector<double> averaged_wall_thickness(num_nodes);
00277     for (unsigned node_index=0; node_index<num_nodes; node_index++)
00278     {
00279         my_averaged_wall_thickness[node_index] = GetAveragedThicknessLocalNode(node_index, wall_thickness);
00280     }
00281 
00282     //Non-local information appear as zeros in the vector
00283     MPI_Allreduce(&my_averaged_wall_thickness[0], &averaged_wall_thickness[0], num_nodes,
00284                   MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00285     if (logInfo)
00286     {
00287         for (unsigned node_index=0; node_index<num_nodes; node_index++)
00288         {
00289              *p_ave_thickness_file << averaged_wall_thickness[node_index] << "\n";
00290         }
00291 
00292     }
00293 
00294     /*
00295      *  Compute the gradient of the averaged wall thickness at the centroid of each tetrahedron
00296      */
00297     c_vector<double,SPACE_DIM> grad_ave_wall_thickness;
00298 
00299     bool fibre_file_is_open= false;
00300     for (unsigned element_index=0; element_index<num_elements; element_index++)
00301     {
00302         //PRINT_VARIABLE(element_index);
00303         if (mrMesh.CalculateDesignatedOwnershipOfElement(element_index))
00304         {
00305             //PRINT_VARIABLE(element_index);
00306             //Locally own this element
00307             //Make sure that writing will proceed
00308             //The same file is written to by multiple processes, therefore we need
00309             //this barrier - Barrier happens BEFORE opening and AFTER closing
00310 
00311             PetscTools::Barrier("GenerateOrthotropicFibreOrientation");
00312             if (fibre_file_is_open != true)
00313             {
00314                 //Open and append
00315                 p_fibre_file = handler.OpenOutputFile(fibreOrientationFile, std::ios::app);
00316                 fibre_file_is_open=true;
00317             }
00318 
00319 
00320             Element<SPACE_DIM,SPACE_DIM>* p_element = mrMesh.GetElement(element_index);
00321 
00322             /*
00323              *  The gradient of the averaged thickness at the element is:
00324              *
00325              *     grad_ave_wall_thickness[element_index] = ave' * BF * inv(J)
00326              *
00327              *  being : ave, averaged thickness values of the nodes defining the element
00328              *          J,   the Jacobian of the element as defined in class Element.
00329              *                               (-1 -1 -1)
00330              *          BF,  basis functions ( 1  0  0)
00331              *                               ( 0  1  0)
00332              *                               ( 0  0  1)
00333              *
00334              *  Defined as u in Streeter paper
00335              */
00336 
00337             c_vector<double, SPACE_DIM+1> elem_nodes_ave_thickness;
00338             double element_averaged_thickness = 0.0;
00339             c_vector<HeartRegionType, SPACE_DIM+1> elem_nodes_region;
00340 
00341             for (unsigned local_node_index=0; local_node_index<SPACE_DIM+1; local_node_index++)
00342             {
00343                 // Get node's global index
00344                 unsigned global_node_index = p_element->GetNode(local_node_index)->GetIndex();
00345 
00346                 elem_nodes_ave_thickness[local_node_index] = averaged_wall_thickness[global_node_index];
00347                 elem_nodes_region[local_node_index] = mpGeometryInfo->GetHeartRegion(global_node_index);
00348 
00349                 // Calculate wall thickness averaged value for the element
00350                 element_averaged_thickness +=  wall_thickness[global_node_index];
00351             }
00352 
00353             element_averaged_thickness /= SPACE_DIM+1;
00354 
00355             c_matrix<double, SPACE_DIM+1, SPACE_DIM> basis_functions( zero_matrix<double>(4u,3u) );
00356             basis_functions(0,0) = basis_functions(0,1) = basis_functions(0,2) = -1.0;
00357             basis_functions(1,0) = basis_functions(2,1) = basis_functions(3,2) =  1.0;
00358 
00359             c_matrix<double, SPACE_DIM+1, SPACE_DIM> temp;
00360             c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian, inverse_jacobian;
00361             double jacobian_det;
00362             mrMesh.GetInverseJacobianForElement(element_index, jacobian, jacobian_det, inverse_jacobian);
00363             noalias(temp) = prod (basis_functions, inverse_jacobian);
00364             noalias(grad_ave_wall_thickness) = prod(elem_nodes_ave_thickness, temp);
00365 
00366             grad_ave_wall_thickness /= norm_2(grad_ave_wall_thickness);
00367 
00368             /*
00369              *   Normal to the gradient (v in Streeter paper) which is then the circumferential direction
00370              * (it will be the fibre direction after rotation)
00371              *
00372              *  Computed as the cross product with the base-apex direction (originally assumed base-apex axis is x). The output vector is not normal,
00373              * since the angle between them may be != 90, normalise it.
00374              */
00375             c_vector<double, SPACE_DIM> fibre_direction = VectorProduct(grad_ave_wall_thickness, mApexToBase);
00376             fibre_direction /= norm_2(fibre_direction);
00377 
00378             /*
00379              *  Longitude direction (w in Streeter paper)
00380              */
00381             c_vector<double, SPACE_DIM> longitude_direction = VectorProduct(grad_ave_wall_thickness, fibre_direction);
00382 
00383             /*
00384              *  Compute fibre to v angle: alpha = R*(1-2e)^3
00385              *
00386              *    R is the maximum angle between the fibre and the v axis (heart region dependant)
00387              *    (1 - 2e)^3 scales it by a value in [-1, 1] defining the rotation of the fibre based
00388              *       on the position in the wall
00389              */
00390             double alpha = GetFibreMaxAngle(elem_nodes_region) * SmallPow( (1 - 2*element_averaged_thickness), 3 );
00391 
00392             /*
00393              *  Apply alpha rotation about the u axis to the orthonormal basis
00394              *
00395              *               ( u(1) v(1) w(1) )
00396              *   (u, v, w) = ( u(2) v(2) w(2) )
00397              *               ( u(3) v(3) w(3) )
00398              *
00399              *  The following matrix defines a rotation about the u axis
00400              *
00401              *                 ( 1        0           0      ) (u')
00402              *   R = (u, v, w) ( 0    cos(alpha) -sin(alpha) ) (v')
00403              *                 ( 0    sin(alpha)  cos(alpha) ) (w')
00404              *
00405              *  The rotated basis is computed like:
00406              *
00407              *                                             ( 1        0           0      )
00408              *  (u, v_r, w_r ) = R * (u, v, w) = (u, v, w) ( 0    cos(alpha) -sin(alpha) )
00409              *                                             ( 0    sin(alpha)  cos(alpha) )
00410              *
00411              *  Which simplifies to:
00412              *
00413              *   v_r =  v*cos(alpha) + w*sin(alpha)
00414              *   w_r = -v*sin(alpha) + w*cos(alpha)
00415              */
00416             c_vector<double, SPACE_DIM> rotated_fibre_direction = fibre_direction*cos(alpha) + longitude_direction*sin(alpha);
00417             c_vector<double, SPACE_DIM> rotated_longitude_direction = -fibre_direction*sin(alpha) + longitude_direction*cos(alpha);
00418 
00419 
00420             /*
00421              * Test the orthonormality of the basis
00422              */
00423             assert( fabs(norm_2(rotated_fibre_direction) - 1) < 100*DBL_EPSILON );
00424             assert( fabs(norm_2(grad_ave_wall_thickness) - 1) < 100*DBL_EPSILON );
00425             assert( fabs(norm_2(rotated_longitude_direction) - 1) < 100*DBL_EPSILON );
00426 
00427             assert( fabs(inner_prod(rotated_fibre_direction, grad_ave_wall_thickness)) < 100*DBL_EPSILON );
00428             assert( fabs(inner_prod(rotated_fibre_direction, rotated_longitude_direction)) < 100*DBL_EPSILON);
00429             assert( fabs(inner_prod(grad_ave_wall_thickness, rotated_longitude_direction)) < 100*DBL_EPSILON);
00430 
00431             /*
00432              *  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)
00433              *
00434              *  See Fig. 1 "Laminar Structure of the Heart: a mathematical model" LeGrice et al. 97
00435              *
00436              */
00437             *p_fibre_file << rotated_fibre_direction[0]     << " " << rotated_fibre_direction[1]     << " "  << rotated_fibre_direction[2]     << " "
00438                           << grad_ave_wall_thickness[0]     << " " << grad_ave_wall_thickness[1]     << " "  << grad_ave_wall_thickness[2]     << " "
00439                           << rotated_longitude_direction[0] << " " << rotated_longitude_direction[1] << " "  << rotated_longitude_direction[2] << std::endl;
00440         }
00441         else
00442         {
00443             //We don't locally own this element
00444             //Give up the file pointer
00445             if (fibre_file_is_open)
00446             {
00447                 p_fibre_file->close();
00448                 fibre_file_is_open=false;
00449             }
00450             PetscTools::Barrier("GenerateOrthotropicFibreOrientation");
00451 
00452         }
00453     }
00454 
00455     if (fibre_file_is_open)
00456     {
00457         p_fibre_file->close();
00458     }
00459 
00460     if (logInfo)
00461     {
00462         p_regions_file->close();
00463         p_thickness_file->close();
00464         p_ave_thickness_file->close();
00465     }
00466 
00467 }
00468 
00469 
00470 template<unsigned SPACE_DIM>
00471 void StreeterFibreGenerator<SPACE_DIM>::SetApexToBase(const c_vector<double, SPACE_DIM>& apexToBase)
00472 {
00473     double norm = norm_2(apexToBase);
00474     if (norm < DBL_EPSILON)
00475     {
00476         EXCEPTION("Apex to base vector should be non-zero");
00477     }
00478     mApexToBase = apexToBase / norm;
00479 }
00480 template<unsigned SPACE_DIM>
00481 void StreeterFibreGenerator<SPACE_DIM>::SetApexToBase(unsigned axis)
00482 {
00483     if (axis >= SPACE_DIM)
00484     {
00485         EXCEPTION("Apex to base coordinate axis was out of range");
00486     }
00487     mApexToBase = zero_vector<double>(SPACE_DIM);
00488     mApexToBase[axis] = 1.0;
00489 }
00490 
00492 // Explicit instantiation
00494 #define COVERAGE_IGNORE
00495 template class StreeterFibreGenerator<3>;
00496 #undef COVERAGE_IGNORE

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5