63 const unsigned nodeIndex,
const std::vector<double>& wallThickness)
const
65 if (nodeIndex < this->mpMesh->GetDistributedVectorFactory()->GetLow() ||
66 nodeIndex >= this->mpMesh->GetDistributedVectorFactory()->GetHigh() )
72 double average = wallThickness[nodeIndex];
73 unsigned nodes_visited = 1;
76 std::set<unsigned> visited_nodes;
77 visited_nodes.insert(nodeIndex);
90 for (
unsigned node_local_index=0;
91 node_local_index<p_containing_element->
GetNumNodes();
95 unsigned neighbour_node_index = p_neighbour_node->
GetIndex();
98 if (visited_nodes.find(neighbour_node_index) == visited_nodes.end())
100 average += wallThickness[neighbour_node_index];
101 visited_nodes.insert(neighbour_node_index);
107 return average/nodes_visited;
189 assert(SPACE_DIM == 3);
190 if (mpGeometryInfo == NULL)
192 EXCEPTION(
"Files defining the heart surfaces not set");
196 out_stream p_regions_file, p_thickness_file, p_ave_thickness_file;
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");
209 if (fabs(norm_2(mApexToBase)) < DBL_EPSILON)
211 EXCEPTION(
"Apex to base vector has not been set");
215 unsigned num_nodes = this->mpMesh->GetNumNodes();
216 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
218 double dist_epi, dist_endo;
220 HeartRegionType node_region = mpGeometryInfo->GetHeartRegion(node_index);
226 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
227 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
232 dist_epi = mpGeometryInfo->rGetDistanceMapEpicardium()[node_index];
233 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
237 dist_epi = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
238 dist_endo = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
242 dist_epi = mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index];
243 dist_endo = mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index];
248 std::cerr <<
"Wrong distances node: " << node_index <<
"\t"
250 <<
"RV " << mpGeometryInfo->rGetDistanceMapRightVentricle()[node_index] <<
"\t"
251 <<
"LV " << mpGeometryInfo->rGetDistanceMapLeftVentricle()[node_index]
264 mWallThickness[node_index] = dist_endo / (dist_endo + dist_epi);
266 if (std::isnan(mWallThickness[node_index]))
273 mWallThickness[node_index] = 0;
279 *p_regions_file << node_region*100 <<
"\n";
280 *p_thickness_file << mWallThickness[node_index] <<
"\n";
287 std::vector<double> my_averaged_wall_thickness(num_nodes);
288 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
290 my_averaged_wall_thickness[node_index] = GetAveragedThicknessLocalNode(node_index, mWallThickness);
294 MPI_Allreduce(&my_averaged_wall_thickness[0], &mAveragedWallThickness[0], num_nodes,
295 MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
299 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
301 *p_ave_thickness_file << mAveragedWallThickness[node_index] <<
"\n";
307 p_regions_file->close();
308 p_thickness_file->close();
309 p_ave_thickness_file->close();
315 unsigned localElementIndex,
316 c_vector<double, SPACE_DIM*SPACE_DIM>& rData)
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;
337 for (
unsigned local_node_index=0; local_node_index<SPACE_DIM+1; local_node_index++)
340 unsigned global_node_index = pElement->
GetNode(local_node_index)->GetIndex();
342 elem_nodes_ave_thickness[local_node_index] = mAveragedWallThickness[global_node_index];
343 elem_nodes_region[local_node_index] = mpGeometryInfo->GetHeartRegion(global_node_index);
346 element_averaged_thickness += mWallThickness[global_node_index];
349 element_averaged_thickness /= SPACE_DIM+1;
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;
355 c_matrix<double, SPACE_DIM+1, SPACE_DIM> temp;
356 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian, inverse_jacobian;
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);
363 grad_ave_wall_thickness /= norm_2(grad_ave_wall_thickness);
372 c_vector<double, SPACE_DIM> fibre_direction =
VectorProduct(grad_ave_wall_thickness, mApexToBase);
373 fibre_direction /= norm_2(fibre_direction);
378 c_vector<double, SPACE_DIM> longitude_direction =
VectorProduct(grad_ave_wall_thickness, fibre_direction);
387 double alpha = GetFibreMaxAngle(elem_nodes_region) * SmallPow( (1 - 2*element_averaged_thickness), 3 );
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);
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 );
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);
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];