36#include "HeartGeometryInformation.hpp"
41#include "OutputFileHandler.hpp"
46template<
unsigned SPACE_DIM>
49template<
unsigned SPACE_DIM>
52template<
unsigned SPACE_DIM>
54 const std::string& rEpiFile,
55 const std::string& rEndoFile,
71template<
unsigned SPACE_DIM>
73 const std::string& rEpiFile,
74 const std::string& rLVFile,
75 const std::string& rRVFile,
93 EXCEPTION(
"At least one of left ventricle or right ventricle files is required");
107template<
unsigned SPACE_DIM>
111 std::ifstream heterogeneity_file;
112 heterogeneity_file.open(nodeHeterogeneityFileName.c_str());
114 if (!(heterogeneity_file.is_open()))
116 heterogeneity_file.close();
117 EXCEPTION(
"Could not open heterogeneities file (" << nodeHeterogeneityFileName <<
")");
120 while (!heterogeneity_file.eof())
124 heterogeneity_file >> value;
127 if ((heterogeneity_file.fail() && !heterogeneity_file.eof()) || value < 0 || value > 2)
129 heterogeneity_file.close();
130 EXCEPTION(
"A value in the heterogeneities file (" << nodeHeterogeneityFileName
131 <<
") is out of range (or not an integer). It should be epi = 0, mid = 1, endo = 2");
134 if (!heterogeneity_file.eof())
138 mLayerForEachNode.push_back(EPI);
142 mLayerForEachNode.push_back(MID);
147 mLayerForEachNode.push_back(ENDO);
152 heterogeneity_file.close();
155template<
unsigned SPACE_DIM>
157 const std::string& rLineFromFile,
158 std::set<unsigned>& rSurfaceNodeIndexSet,
159 unsigned offset)
const
161 std::stringstream line_stream(rLineFromFile);
162 while (!line_stream.eof())
167 if (item == 0 && offset != 0)
169 EXCEPTION(
"Error when reading surface file. It was assumed not to be indexed from zero, but zero appeared in the list.");
171 rSurfaceNodeIndexSet.insert(item-offset);
175template<
unsigned SPACE_DIM>
177 const std::string& rSurfaceFileName,
178 std::vector<unsigned>& rSurfaceNodes,
179 bool indexFromZero)
const
182 std::ifstream file_stream;
184 if (indexFromZero ==
false)
189 file_stream.open(rSurfaceFileName.c_str());
190 if (!file_stream.is_open())
192 EXCEPTION(
"Wrong surface definition file name " + rSurfaceFileName);
196 std::set<unsigned> surface_original_node_index_set;
200 getline(file_stream, line);
203 ProcessLine(line, surface_original_node_index_set, offset);
205 getline(file_stream, line);
207 while (!file_stream.eof());
211 rSurfaceNodes.reserve(surface_original_node_index_set.size());
213 if (mpMesh->rGetNodePermutation().empty())
216 for (std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
217 node_index_it != surface_original_node_index_set.end();
220 rSurfaceNodes.push_back(*node_index_it);
226 for (std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
227 node_index_it != surface_original_node_index_set.end();
230 rSurfaceNodes.push_back(mpMesh->rGetNodePermutation()[*node_index_it]);
235template<
unsigned SPACE_DIM>
239 if (mDistMapRightVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
240 mDistMapRightVentricle[nodeIndex] >= mDistMapLeftVentricle[nodeIndex])
242 return LEFT_VENTRICLE_WALL;
245 if (mDistMapLeftVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
246 mDistMapLeftVentricle[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
248 return RIGHT_VENTRICLE_WALL;
251 if (mDistMapEpicardium[nodeIndex] >= mDistMapLeftVentricle[nodeIndex] &&
252 mDistMapEpicardium[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
254 if (mDistMapLeftVentricle[nodeIndex]
255 < LEFT_SEPTUM_SIZE*(mDistMapLeftVentricle[nodeIndex] + mDistMapRightVentricle[nodeIndex]))
269template<
unsigned SPACE_DIM>
273 if (mNumberOfSurfacesProvided == 3)
278 case LEFT_VENTRICLE_WALL:
279 case LEFT_VENTRICLE_SURFACE:
280 return mDistMapLeftVentricle[nodeIndex];
283 case RIGHT_VENTRICLE_WALL:
284 case RIGHT_VENTRICLE_SURFACE:
285 return mDistMapRightVentricle[nodeIndex];
289 return mDistMapLeftVentricle[nodeIndex];
293 return mDistMapRightVentricle[nodeIndex] ;
298 std::cerr <<
"Wrong distances node: " << nodeIndex <<
"\t"
299 <<
"Epi " << mDistMapEpicardium[nodeIndex] <<
"\t"
300 <<
"RV " << mDistMapRightVentricle[nodeIndex] <<
"\t"
301 <<
"LV " << mDistMapLeftVentricle[nodeIndex]
316 return mDistMapEndocardium[nodeIndex];
324template<
unsigned SPACE_DIM>
327 return mDistMapEpicardium[nodeIndex];
330template<
unsigned SPACE_DIM>
334 double dist_endo = GetDistanceToEndo(nodeIndex);
335 double dist_epi = GetDistanceToEpi(nodeIndex);
336 assert( (dist_endo + dist_epi) != 0 );
341 double relative_position = dist_endo / (dist_endo + dist_epi);
342 return relative_position;
345template<
unsigned SPACE_DIM>
348 if (epiFraction+endoFraction>1)
350 EXCEPTION(
"The sum of fractions of epicardial and endocardial layers must be lesser than 1");
353 if ((endoFraction<0) || (epiFraction<0))
355 EXCEPTION(
"A fraction of a layer must be positive");
358 mLayerForEachNode.resize(mpMesh->GetNumNodes());
359 for (
unsigned i=0; i<mpMesh->GetNumNodes(); i++)
361 double position = CalculateRelativeWallPosition(i);
362 if (position<endoFraction)
364 mLayerForEachNode[i] = ENDO;
366 else if (position<(1-epiFraction))
368 mLayerForEachNode[i] = MID;
372 mLayerForEachNode[i] = EPI;
377template<
unsigned SPACE_DIM>
385 assert(mLayerForEachNode.size()>0);
386 for (
unsigned i=0; i<mpMesh->GetNumNodes(); i++)
388 if (mLayerForEachNode[i]==EPI)
392 else if (mLayerForEachNode[i]==MID)
407template<
unsigned SPACE_DIM>
409 const std::vector<unsigned>& rSurfaceNodes)
412 assert(rSurfaceNodes.size()>0);
414 c_vector<double, SPACE_DIM> my_minimum_point = scalar_vector<double>(SPACE_DIM, DBL_MAX);
416 c_vector<double, SPACE_DIM> my_maximum_point=-my_minimum_point;
419 for (
unsigned surface_index=0; surface_index<rSurfaceNodes.size(); surface_index++)
421 unsigned global_index=rSurfaceNodes[surface_index];
422 if (mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index) )
424 const c_vector<double, SPACE_DIM>& r_position = mpMesh->GetNode(global_index)->rGetLocation();
426 for (
unsigned i=0; i<SPACE_DIM; i++)
428 if (r_position[i] < my_minimum_point[i])
430 my_minimum_point[i] = r_position[i];
432 if (r_position[i] > my_maximum_point[i])
434 my_maximum_point[i] = r_position[i];
441 c_vector<double, SPACE_DIM> global_minimum_point;
442 c_vector<double, SPACE_DIM> global_maximum_point;
443 MPI_Allreduce(&my_minimum_point[0], &global_minimum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
444 MPI_Allreduce(&my_maximum_point[0], &global_maximum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
#define EXCEPTION(message)
void ComputeDistanceMap(const std::vector< unsigned > &rSourceNodeIndices, std::vector< double > &rNodeDistances)
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const