Chaste  Release::2018.1
PapillaryFibreCalculator.hpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 #ifndef PAPILLARYFIBRECALCULATOR_HPP_
36 #define PAPILLARYFIBRECALCULATOR_HPP_
37 
38 #include "UblasCustomFunctions.hpp"
39 #include "TetrahedralMesh.hpp"
40 
47 {
48 // Allow the test class to use the private functions.
49 friend class TestPapillaryFibreCalculator;
50 
51 private:
55  std::vector< c_vector<double, 3> > mRadiusVectors;
57  std::vector< c_matrix<double,3,3> > mStructureTensors;
59  std::vector< c_matrix<double,3,3> > mSmoothedStructureTensors;
60 
68  c_vector<double,3> GetRadiusVectorForOneElement(unsigned elementIndex);
69 
74  void GetRadiusVectors();
75 
76 
83 
94 
95 public:
102 
108  std::vector<c_vector<double,3> > CalculateFibreOrientations();
109 };
110 
111 // PUBLIC METHODS
113  : mrMesh(rMesh)
114 {
118 }
119 
121 {
123 
125 
127 
128  // Calculate eigenvalues
129  std::vector<c_vector<double,3> > fibre_orientations(mrMesh.GetNumElements());
130  for (unsigned i=0; i<fibre_orientations.size(); i++)
131  {
133  }
134 
135  return fibre_orientations;
136 }
137 
138 // PRIVATE METHODS
139 c_vector<double,3> PapillaryFibreCalculator::GetRadiusVectorForOneElement(unsigned elementIndex)
140 {
141  c_vector<double, 3> centroid = (mrMesh.GetElement(elementIndex))->CalculateCentroid();
142  // Loops over all papillary face nodes
143  c_vector<double,3> coordinates;
144 
145  double nearest_r_squared=DBL_MAX;
146  unsigned nearest_face_node = 0;
147 
149  while (bound_node_iter != mrMesh.GetBoundaryNodeIteratorEnd())
150  {
151  unsigned bound_node_index = (*bound_node_iter)->GetIndex();
152  coordinates=mrMesh.GetNode(bound_node_index)->rGetLocation();
153 
154  // Calculates the distance between the papillary face node and the centroid
155  double r_squared = norm_2(centroid-coordinates);
156  // Checks to see if it is the smallest so far - if it is, update the current smallest distance
157  if (r_squared < nearest_r_squared)
158  {
159  nearest_r_squared = r_squared;
160  nearest_face_node = bound_node_index;
161  }
162  ++bound_node_iter;
163  }
164 
165  coordinates = mrMesh.GetNode(nearest_face_node)->rGetLocation();
166  c_vector<double,3> radial_vector = coordinates-centroid;
167  return radial_vector;
168 }
169 
171 {
172  // Loops over all elements finding radius vector
174  iter != mrMesh.GetElementIteratorEnd();
175  ++iter)
176  {
177  unsigned element_index = iter->GetIndex();
178  mRadiusVectors[element_index] = GetRadiusVectorForOneElement(element_index);
179  }
180 }
181 
183 {
184  for (unsigned i=0;i<mRadiusVectors.size();i++)
185  {
186  mStructureTensors[i] = outer_prod(mRadiusVectors[i],mRadiusVectors[i]);
187  }
188 }
189 
191 {
192  const double sigma = 0.05; //cm
193  const double r_max = 0.1; //cm
194  double g_factor_sum = 0;
195  double g_factor = 0;
196 
198  elem_iter != mrMesh.GetElementIteratorEnd();
199  ++elem_iter)
200  {
201  mSmoothedStructureTensors[ elem_iter->GetIndex()] = zero_matrix<double>(3,3);
202 
203  c_vector<double, 3> centroid = elem_iter->CalculateCentroid();
204  g_factor_sum = 0;
205 
207  iter_2 != mrMesh.GetElementIteratorEnd();
208  ++iter_2)
209  {
210  c_vector<double, 3> centroid_2 = iter_2->CalculateCentroid();
211  double r = norm_2(centroid-centroid_2);
212  if (r < r_max)
213  {
214  g_factor = exp(-r/(2*sigma*sigma));
215 
216  g_factor_sum += g_factor;
217 
218  mSmoothedStructureTensors[elem_iter->GetIndex()] += g_factor*mStructureTensors[iter_2->GetIndex()];
219  }
220  }
221 
222  mSmoothedStructureTensors[elem_iter->GetIndex()] /= g_factor_sum;
223  }
224 }
225 
226 #endif /*PAPILLARYFIBRECALCULATOR_HPP_*/
227 
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
c_vector< double, 3 > CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix< double, 3, 3 > &rA)
Node< SPACE_DIM > * GetNode(unsigned index) const
virtual unsigned GetNumElements() const
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
BoundaryNodeIterator GetBoundaryNodeIteratorBegin() const
c_vector< double, 3 > GetRadiusVectorForOneElement(unsigned elementIndex)
std::vector< c_matrix< double, 3, 3 > > mSmoothedStructureTensors
std::vector< c_vector< double, 3 > > CalculateFibreOrientations()
std::vector< c_matrix< double, 3, 3 > > mStructureTensors
PapillaryFibreCalculator(TetrahedralMesh< 3, 3 > &rMesh)
BoundaryNodeIterator GetBoundaryNodeIteratorEnd() const
TetrahedralMesh< 3, 3 > & mrMesh
std::vector< c_vector< double, 3 > > mRadiusVectors