Chaste  Build::
GreensFunctionSolver.cpp
1 /*
2 
3 Copyright (c) 2005-2016, 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 
36 #define _BACKWARD_BACKWARD_WARNING_H 1 //Cut out the vtk deprecated warning for now (gcc4.3)
37 #include <vtkDoubleArray.h>
38 #include <vtkPointData.h>
39 #include <vtkSmartPointer.h>
40 #include "Vessel.hpp"
41 #include "VesselSegment.hpp"
42 #include "ChastePoint.hpp"
43 #include "LinearSystem.hpp"
44 #include "ReplicatableVector.hpp"
45 #include "UblasMatrixInclude.hpp"
46 #include "UnitCollection.hpp"
47 #include "RegularGridWriter.hpp"
48 #include "GeometryWriter.hpp"
49 #include "GreensFunctionSolver.hpp"
50 #include "BaseUnits.hpp"
51 
52 template<unsigned DIM>
55  mpDomain(),
56  mSinkCoordinates(),
57  mSinkPointMap(),
58  mSubSegmentCoordinates(),
59  mSubSegmentLengths(),
60  mSinkRates(),
61  mSourceRates(),
62  mSegmentConcentration(),
63  mSegmentPointMap(),
64  mGtt(),
65  mGvv(),
66  mGvt(),
67  mGtv(),
68  mSubsegmentCutoff(1.0*unit::microns)
69 {
70 
71 }
72 
73 template<unsigned DIM>
75 {
76 
77 }
78 
79 template<unsigned DIM>
80 void GreensFunctionSolver<DIM>::SetSubSegmentCutoff(units::quantity<unit::length> value)
81 {
82  mSubsegmentCutoff = value;
83 }
84 
85 template<unsigned DIM>
87 {
88  // Set up the sub-segment and tissue point co-ordinates
91 
92  // Generate the greens function matrices
93  UpdateGreensFunctionMatrices(1, 1, 1, 1);
94 
95  // Get the sink rates
96  unsigned number_of_sinks = mSinkCoordinates.size();
97  units::quantity<unit::concentration_flow_rate> sink_rate = this->mpPde->ComputeConstantInUSourceTerm();
98  units::quantity<unit::volume> sink_volume = units::pow<3>(this->mpRegularGrid->GetSpacing());
99  mSinkRates = std::vector<units::quantity<unit::molar_flow_rate> >(number_of_sinks, sink_rate * sink_volume);
100  units::quantity<unit::molar_flow_rate> total_sink_rate = std::accumulate(mSinkRates.begin(), mSinkRates.end(), 0.0*unit::mole_per_second);
101 
102  // Get the sink substance demand on each vessel subsegment
103  unsigned number_of_subsegments = mSubSegmentCoordinates.size();
104  units::quantity<unit::diffusivity> diffusivity = this->mpPde->ComputeIsotropicDiffusionTerm();
105  std::vector<units::quantity<unit::concentration> > sink_demand_per_subsegment(number_of_subsegments, 0.0*this->mReferenceConcentration);
106  for (unsigned idx = 0; idx < number_of_subsegments; idx++)
107  {
108  for (unsigned jdx = 0; jdx < number_of_sinks; jdx++)
109  {
110  sink_demand_per_subsegment[idx] += ((*mGvt)[idx][jdx] / diffusivity) * mSinkRates[jdx];
111  }
112  }
113 
114  mSegmentConcentration = std::vector<units::quantity<unit::concentration> >(number_of_subsegments, 1.0*this->mReferenceConcentration);
115  this->mConcentrations = std::vector<units::quantity<unit::concentration> >(number_of_sinks, 0.0*this->mReferenceConcentration);
116 
117  // Solve for the subsegment source rates required to meet the sink substance demand
118  double tolerance = 1.e-10;
119  units::quantity<unit::concentration> g0 = 0.0 * this->mReferenceConcentration;
120  mSourceRates = std::vector<units::quantity<unit::molar_flow_rate> >(number_of_subsegments, 0.0*unit::mole_per_second);
121 
122  units::quantity<unit::time> reference_time = BaseUnits::Instance()->GetReferenceTimeScale();
123  units::quantity<unit::concentration> reference_concentration = this->mReferenceConcentration;
124  units::quantity<unit::amount> reference_amount(1.0*unit::moles);
125  LinearSystem linear_system(number_of_subsegments + 1, number_of_subsegments + 1);
126  linear_system.SetKspType("bcgs");
127 
128  for (unsigned iteration = 0; iteration < 10; iteration++)
129  {
130  linear_system.AssembleIntermediateLinearSystem();
131  for (unsigned i = 0; i < number_of_subsegments; i++)
132  {
133  linear_system.SetRhsVectorElement(i, (mSegmentConcentration[i] - sink_demand_per_subsegment[i])/reference_concentration);
134  }
135 
136  linear_system.SetRhsVectorElement(number_of_subsegments, -total_sink_rate*(reference_time/reference_amount));
137 
138  // Set up Linear system matrix
139  for (unsigned iter = 0; iter < number_of_subsegments; iter++)
140  {
141  for (unsigned jter = 0; jter < number_of_subsegments; jter++)
142  {
143  linear_system.SetMatrixElement(iter, jter, ((*mGvv)[iter][jter] / diffusivity)*(reference_amount/(reference_time*reference_concentration)));
144  }
145  linear_system.SetMatrixElement(number_of_subsegments, iter, 1.0);
146  linear_system.SetMatrixElement(iter, number_of_subsegments, 1.0);
147  }
148  linear_system.SetMatrixElement(number_of_subsegments, number_of_subsegments, 0.0);
149 
150  // Solve the linear system
151  linear_system.AssembleFinalLinearSystem();
152  ReplicatableVector soln_repl(linear_system.Solve());
153 
154  // Populate the solution vector
155  std::vector<units::quantity<unit::molar_flow_rate> > solution_vector(number_of_subsegments + 1);
156  for (unsigned row = 0; row < number_of_subsegments + 1; row++)
157  {
158  (solution_vector)[row] = soln_repl[row]*(reference_amount/reference_time);
159  }
160 
161  // Check convergence
162  bool all_in_tolerance = true;
163  for (unsigned i = 0; i < number_of_subsegments; i++)
164  {
165  units::quantity<unit::molar_flow_rate> diff = units::abs(mSourceRates[i] - solution_vector[i]);
166  if (diff > tolerance*unit::mole_per_second)
167  {
168  all_in_tolerance = false;
169  break;
170  }
171  }
172  // Retrieve the solution
173  for (unsigned i = 0; i < number_of_subsegments; i++)
174  {
175  mSourceRates[i] = solution_vector[i];
176  }
177  g0 = solution_vector[number_of_subsegments]*reference_concentration*(reference_time/reference_amount);
178 
179  if (all_in_tolerance)
180  {
181  break;
182  }
183  else
184  {
185  if (iteration == 9)
186  {
187  std::cout << "Did not converge\n";
188  }
189  }
190  }
191 
192  // Get the tissue concentration and write the solution
193  for (unsigned i = 0; i < number_of_sinks; i++)
194  {
195  this->mConcentrations[i] = 0.0 * this->mReferenceConcentration;
196  for (unsigned j = 0; j < number_of_sinks; j++)
197  {
198  this->mConcentrations[i] += (*mGtt)[i][j] * mSinkRates[j] / diffusivity;
199  }
200 
201  for (unsigned j = 0; j < number_of_subsegments; j++)
202  {
203  this->mConcentrations[i] += (*mGtv)[i][j] * mSourceRates[j] / diffusivity;
204  }
205  this->mConcentrations[i] += g0;
206  }
207 
208  std::map<std::string, std::vector<units::quantity<unit::concentration> > > segmentPointData;
209  segmentPointData[this->mLabel] = mSegmentConcentration;
210 
211  this->UpdateSolution(this->mConcentrations);
212  if(this->mWriteSolution)
213  {
214  this->WriteSolution(segmentPointData);
215  }
216 }
217 
218 template<unsigned DIM>
220 {
221  // Set up the sub-segment points and map to original segments
222  units::quantity<unit::length> max_subsegment_length = mSubsegmentCutoff;
223 
224  std::vector<boost::shared_ptr<Vessel<DIM> > > vessels = this->mpNetwork->GetVessels();
225  typename std::vector<boost::shared_ptr<Vessel<DIM> > >::iterator vessel_iter;
226  typename std::vector<boost::shared_ptr<VesselSegment<DIM> > >::iterator segment_iter;
227 
228  // Iterate over all segments and store midpoints and lengths of subsegment regions for
229  // the greens functions calculation. Create a map of subsegment index to the parent segment
230  // for later use.
231  for (vessel_iter = vessels.begin(); vessel_iter != vessels.end(); vessel_iter++)
232  {
233  std::vector<boost::shared_ptr<VesselSegment<DIM> > > segments = (*vessel_iter)->GetSegments();
234  for (segment_iter = segments.begin(); segment_iter != segments.end(); segment_iter++)
235  {
236  units::quantity<unit::length> segment_length = (*segment_iter)->GetLength();
237 
238  // If the segment is shorter than the max length just use its mid-point
239  if (segment_length < 1.01 * max_subsegment_length)
240  {
241  mSubSegmentCoordinates.push_back((*segment_iter)->GetMidPoint());
242  mSubSegmentLengths.push_back(segment_length);
243  mSegmentPointMap[mSubSegmentCoordinates.size() - 1] = (*segment_iter);
244  }
245  // Otherwise generate subsegment points along its length
246  else
247  {
248  DimensionalChastePoint<DIM> start_point = (*segment_iter)->GetNode(0)->rGetLocation();
249  DimensionalChastePoint<DIM> end_point = (*segment_iter)->GetNode(1)->rGetLocation();
250 
251  double subsegment_length_factor = segment_length / max_subsegment_length;
252  unsigned num_subsegments = std::floor(subsegment_length_factor) + 1;
253  units::quantity<unit::length> subsegment_length = segment_length / double(num_subsegments);
254  DimensionalChastePoint<DIM> increment = (end_point - start_point);
255  for (unsigned i = 0; i < num_subsegments; i++)
256  {
257  mSubSegmentCoordinates.push_back(start_point + (increment*(double(i) + 0.5)));
258  mSubSegmentLengths.push_back(subsegment_length);
259  mSegmentPointMap[mSubSegmentCoordinates.size() - 1] = (*segment_iter);
260  }
261  }
262  }
263  }
264 }
265 
266 template<unsigned DIM>
268 {
269  unsigned num_points = this->mpRegularGrid->GetNumberOfPoints();
270  mSinkCoordinates = std::vector<DimensionalChastePoint<DIM> >(num_points);
271  mSinkPointMap = std::vector<unsigned>(num_points);
272  for(unsigned idx=0; idx<num_points; idx++)
273  {
274  mSinkCoordinates[idx] = this->mpRegularGrid->GetLocationOf1dIndex(idx);
275  mSinkPointMap[idx] = idx;
276  }
277 }
278 
279 template<unsigned DIM>
280 void GreensFunctionSolver<DIM>::UpdateGreensFunctionMatrices(bool updateGtt, bool updateGvv, bool updateGtv,
281  bool updateGvt)
282 {
283  // Get the Greens Function coefficient matrices
284  if (updateGtt)
285  {
287  }
288  if (updateGvv)
289  {
291  }
292  if (updateGtv)
293  {
295  }
296  if (updateGvt)
297  {
299  }
300 }
301 
302 template<unsigned DIM>
303 boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > GreensFunctionSolver<DIM>::GetVesselVesselInteractionMatrix()
304 {
305  typedef boost::multi_array<units::quantity<unit::per_length>, 2>::index index;
306  unsigned num_sub_segments = mSubSegmentCoordinates.size();
307  double coefficient = 1.0 / (4.0 * M_PI);
308 
309  boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > p_interaction_matrix(new boost::multi_array<units::quantity<unit::per_length>, 2>(boost::extents[num_sub_segments][num_sub_segments]));
310  for (index iter = 0; iter < num_sub_segments; iter++)
311  {
312  for (index iter2 = 0; iter2 < num_sub_segments; iter2++)
313  {
314  if (iter <= iter2)
315  {
316  units::quantity<unit::length> distance = mSubSegmentCoordinates[iter2].GetDistance(mSubSegmentCoordinates[iter]);
317  units::quantity<unit::per_length> term;
318  if (distance < mSegmentPointMap[iter]->GetRadius())
319  {
320  units::quantity<unit::length> radius = mSegmentPointMap[iter]->GetRadius();
321  units::quantity<unit::length> max_segment_length = std::max(mSubSegmentLengths[iter], mSubSegmentLengths[iter2]);
322  double green_correction = 0.6 * std::exp(-0.45 * max_segment_length /radius);
323  if (iter != iter2)
324  {
325  distance = radius;
326  }
327  term = (1.298 / (1.0 + 0.297 * pow(max_segment_length /radius, 0.838))- green_correction * pow(distance / radius, 2))*coefficient/radius;
328  }
329  else
330  {
331  term = coefficient / distance;
332  }
333  (*p_interaction_matrix)[iter][iter2] = term;
334  (*p_interaction_matrix)[iter2][iter] = term;
335  }
336  }
337  }
338  return p_interaction_matrix;
339 }
340 
341 template<unsigned DIM>
342 boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > GreensFunctionSolver<DIM>::GetTissueTissueInteractionMatrix()
343 {
344  typedef boost::multi_array<units::quantity<unit::per_length>, 2>::index index;
345  unsigned num_points = mSinkCoordinates.size();
346  double coefficient = 1.0 / (4.0 * M_PI);
347  units::quantity<unit::volume> tissue_point_volume = units::pow<3>(this->mpRegularGrid->GetSpacing());
348  units::quantity<unit::length> equivalent_tissue_point_radius = units::root<3>(tissue_point_volume * 0.75 / M_PI);
349  boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2 > > p_interaction_matrix(new boost::multi_array<units::quantity<unit::per_length>, 2>(boost::extents[num_points][num_points]));
350  for (index iter = 0; iter < num_points; iter++)
351  {
352  for (index iter2 = 0; iter2 < num_points; iter2++)
353  {
354  if (iter < iter2)
355  {
356  units::quantity<unit::length> distance = mSinkCoordinates[iter2].GetDistance(mSinkCoordinates[iter]);
357  (*p_interaction_matrix)[iter][iter2] = coefficient / distance;
358  (*p_interaction_matrix)[iter2][iter] = coefficient / distance;
359  }
360  else if (iter == iter2)
361  {
362  (*p_interaction_matrix)[iter][iter2] = 1.2 * coefficient / equivalent_tissue_point_radius;
363  }
364  }
365  }
366  return p_interaction_matrix;
367 }
368 
369 template<unsigned DIM>
370 boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > GreensFunctionSolver<DIM>::GetTissueVesselInteractionMatrix()
371 {
372  typedef boost::multi_array<units::quantity<unit::per_length>, 2>::index index;
373  unsigned num_sinks = mSinkCoordinates.size();
374  unsigned num_subsegments = mSubSegmentCoordinates.size();
375 
376  units::quantity<unit::volume> tissue_point_volume = units::pow<3>(this->mpRegularGrid->GetSpacing());
377  units::quantity<unit::length> equivalent_tissue_point_radius = units::root<3>(tissue_point_volume * 0.75 / M_PI);
378  double coefficient = 1.0 / (4.0 * M_PI);
379 
380  boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > p_interaction_matrix(new boost::multi_array<units::quantity<unit::per_length>, 2>(boost::extents[num_sinks][num_subsegments]));
381  for (index iter = 0; iter < num_sinks; iter++)
382  {
383  for (index iter2 = 0; iter2 < num_subsegments; iter2++)
384  {
385  units::quantity<unit::length> distance = mSinkCoordinates[iter2].GetDistance(mSinkCoordinates[iter]);
386  units::quantity<unit::per_length> term;
387  if (distance <= equivalent_tissue_point_radius)
388  {
389  term = coefficient * (1.5 - 0.5 * (pow(distance / equivalent_tissue_point_radius, 2))) / equivalent_tissue_point_radius;
390  }
391  else
392  {
393  term = coefficient / distance;
394  }
395  (*p_interaction_matrix)[iter][iter2] = term;
396  }
397  }
398  return p_interaction_matrix;
399 }
400 
401 template<unsigned DIM>
402 boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > GreensFunctionSolver<DIM>::GetVesselTissueInteractionMatrix()
403 {
404  typedef boost::multi_array<units::quantity<unit::per_length>, 2>::index index;
405  unsigned num_subsegments = mSubSegmentCoordinates.size();
406  unsigned num_sinks = mSinkCoordinates.size();
407  double coefficient = 1.0 / (4.0 * M_PI);
408 
409  units::quantity<unit::volume> tissue_point_volume = units::pow<3>(this->mpRegularGrid->GetSpacing());
410  units::quantity<unit::length> equivalent_tissue_point_radius = units::root<3>(tissue_point_volume * 0.75 / M_PI);
411 
412  boost::shared_ptr<boost::multi_array<units::quantity<unit::per_length>, 2> > p_interaction_matrix(new boost::multi_array<units::quantity<unit::per_length>, 2>(boost::extents[num_subsegments][num_sinks]));
413  for (index iter = 0; iter < num_subsegments; iter++)
414  {
415  for (index iter2 = 0; iter2 < num_sinks; iter2++)
416  {
417  units::quantity<unit::length> distance = mSinkCoordinates[iter2].GetDistance(mSinkCoordinates[iter]);
418  units::quantity<unit::per_length> term;
419  if (distance <= equivalent_tissue_point_radius)
420  {
421  term = coefficient * (1.5 - 0.5 * (pow(distance / equivalent_tissue_point_radius, 2)))/ equivalent_tissue_point_radius;
422  }
423  else
424  {
425  term = coefficient / distance;
426  }
427  (*p_interaction_matrix)[iter][iter2] = term;
428  }
429  }
430  return p_interaction_matrix;
431 }
432 
433 template<unsigned DIM>
434 void GreensFunctionSolver<DIM>::WriteSolution(std::map<std::string, std::vector<units::quantity<unit::concentration> > >& segmentPointData)
435 {
436  // Write the tissue point data
437  RegularGridWriter writer;
438  writer.SetFilename(this->mpOutputFileHandler->GetOutputDirectoryFullPath() + "/solution.vti");
439  writer.SetImage(this->mpVtkSolution);
440  writer.Write();
441 
442  // Add the segment points
443  vtkSmartPointer<vtkPolyData> pPolyData = vtkSmartPointer<vtkPolyData>::New();
444  vtkSmartPointer<vtkPoints> pPoints = vtkSmartPointer<vtkPoints>::New();
445  for (unsigned i = 0; i < mSubSegmentCoordinates.size(); i++)
446  {
447  c_vector<double, DIM> location = mSubSegmentCoordinates[i].GetLocation(this->mpRegularGrid->GetReferenceLengthScale());
448  pPoints->InsertNextPoint(location[0], location[1], location[2]);
449  }
450  pPolyData->SetPoints(pPoints);
451 
452  // Add the segment point data
453  std::map<std::string, std::vector<units::quantity<unit::concentration> > >::iterator segment_iter;
454  for (segment_iter = segmentPointData.begin(); segment_iter != segmentPointData.end(); ++segment_iter)
455  {
456  vtkSmartPointer<vtkDoubleArray> pInfo = vtkSmartPointer<vtkDoubleArray>::New();
457  pInfo->SetNumberOfComponents(1);
458  pInfo->SetNumberOfTuples(mSubSegmentCoordinates.size());
459  pInfo->SetName(segment_iter->first.c_str());
460 
461  for (unsigned i = 0; i < mSubSegmentCoordinates.size(); i++)
462  {
463  pInfo->SetValue(i, segmentPointData[segment_iter->first][i]/this->mReferenceConcentration);
464  }
465  pPolyData->GetPointData()->AddArray(pInfo);
466  }
467 
468  GeometryWriter geometry_writer;
469  geometry_writer.SetFileName(this->mpOutputFileHandler->GetOutputDirectoryFullPath() + "/segments.vtp");
470  geometry_writer.SetInput(pPolyData);
471  geometry_writer.Write();
472 }
473 
474 // Explicit instantiation
475 template class GreensFunctionSolver<2>;
476 template class GreensFunctionSolver<3>;
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > mGtt
std::vector< units::quantity< unit::length > > mSubSegmentLengths
boost::shared_ptr< AbstractDiscreteContinuumLinearEllipticPde< DIM, DIM > > mpPde
std::vector< units::quantity< unit::molar_flow_rate > > mSourceRates
units::quantity< unit::length > mSubsegmentCutoff
void SetSubSegmentCutoff(units::quantity< unit::length > value)
std::vector< DimensionalChastePoint< DIM > > mSubSegmentCoordinates
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > GetTissueVesselInteractionMatrix()
void SetFilename(const std::string &rFilename)
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > mGvt
std::vector< units::quantity< unit::concentration > > mSegmentConcentration
units::quantity< unit::time > GetReferenceTimeScale()
Definition: BaseUnits.cpp:72
boost::shared_ptr< OutputFileHandler > mpOutputFileHandler
std::vector< units::quantity< unit::molar_flow_rate > > mSinkRates
std::map< unsigned, boost::shared_ptr< VesselSegment< DIM > > > mSegmentPointMap
void SetInput(vtkSmartPointer< vtkPolyData > pSurface)
std::vector< DimensionalChastePoint< DIM > > mSinkCoordinates
std::vector< units::quantity< unit::concentration > > mConcentrations
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > mGtv
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > mGvv
static BaseUnits * Instance()
Definition: BaseUnits.cpp:53
void UpdateGreensFunctionMatrices(bool updateGtt=0, bool updateGvv=0, bool updateGtv=0, bool updateGvt=0)
boost::shared_ptr< VesselNetwork< DIM > > mpNetwork
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > GetVesselVesselInteractionMatrix()
units::quantity< unit::concentration > mReferenceConcentration
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > GetVesselTissueInteractionMatrix()
void SetImage(vtkSmartPointer< vtkImageData > pImage)
boost::shared_ptr< boost::multi_array< units::quantity< unit::per_length >, 2 > > GetTissueTissueInteractionMatrix()
void SetFileName(const std::string &rFileName)
void WriteSolution(std::map< std::string, std::vector< units::quantity< unit::concentration > > > &segmentPointData)
std::vector< unsigned > mSinkPointMap