Chaste  Build::
AbstractUnstructuredGridDiscreteContinuumSolver.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 <vtkPolyData.h>
40 #include <vtkProbeFilter.h>
41 #include <vtkImageData.h>
42 #include <vtkSmartPointer.h>
43 #include "AbstractUnstructuredGridDiscreteContinuumSolver.hpp"
44 
45 template<unsigned DIM>
48  mpVtkSolution(),
49  mpMesh()
50 {
51  this->mHasUnstructuredGrid = true;
52 }
53 
54 template<unsigned DIM>
56 {
57 
58 }
59 
60 template<unsigned DIM>
61 boost::shared_ptr<DiscreteContinuumMesh<DIM> > AbstractUnstructuredGridDiscreteContinuumSolver<DIM>::GetMesh()
62 {
63  if(!this->mpMesh)
64  {
65  EXCEPTION("A mesh has not been set.");
66  }
67  return this->mpMesh;
68 }
69 
70 template<unsigned DIM>
72 {
73  if(!this->mpVtkSolution)
74  {
75  this->Setup();
76  }
77 
78  std::vector<c_vector<double, DIM> > centroids = mpMesh->GetElementCentroids();
79  std::vector<units::quantity<unit::concentration> > sampled_solution(centroids.size(), 0.0*this->mReferenceConcentration);
80 
81  // Sample the field at these locations
82  vtkSmartPointer<vtkPolyData> p_polydata = vtkSmartPointer<vtkPolyData>::New();
83  vtkSmartPointer<vtkPoints> p_points = vtkSmartPointer<vtkPoints>::New();
84  p_points->SetNumberOfPoints(centroids.size());
85  for(unsigned idx=0; idx< centroids.size(); idx++)
86  {
87  if(DIM==3)
88  {
89  p_points->SetPoint(idx, centroids[idx][0], centroids[idx][1], centroids[idx][2]);
90  }
91  else
92  {
93  p_points->SetPoint(idx, centroids[idx][0], centroids[idx][1], 0.0);
94  }
95  }
96  p_polydata->SetPoints(p_points);
97 
98  vtkSmartPointer<vtkProbeFilter> p_probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
99  #if VTK_MAJOR_VERSION <= 5
100  p_probe_filter->SetInput(p_polydata);
101  p_probe_filter->SetSource(this->mpVtkSolution);
102  #else
103  p_probe_filter->SetInputData(p_polydata);
104  p_probe_filter->SetSourceData(this->mpVtkSolution);
105  #endif
106  p_probe_filter->Update();
107  vtkSmartPointer<vtkPointData> p_point_data = p_probe_filter->GetOutput()->GetPointData();
108 
109  unsigned num_points = p_point_data->GetArray(this->mLabel.c_str())->GetNumberOfTuples();
110  for(unsigned idx=0; idx<num_points; idx++)
111  {
112  sampled_solution[idx] = p_point_data->GetArray(this->mLabel.c_str())->GetTuple1(idx)*this->mReferenceConcentration;
113  }
114  return sampled_solution;
115 }
116 
117 template<unsigned DIM>
118 std::vector<units::quantity<unit::concentration> > AbstractUnstructuredGridDiscreteContinuumSolver<DIM>::GetConcentrations(const std::vector<DimensionalChastePoint<DIM> >& samplePoints)
119 {
120  if(!this->mpVtkSolution)
121  {
122  this->Setup();
123  }
124 
125  std::vector<units::quantity<unit::concentration> > sampled_solution(samplePoints.size(), 0.0*this->mReferenceConcentration);
126 
127  // Sample the field at these locations
128  vtkSmartPointer<vtkPolyData> p_polydata = vtkSmartPointer<vtkPolyData>::New();
129  vtkSmartPointer<vtkPoints> p_points = vtkSmartPointer<vtkPoints>::New();
130 
131  p_points->SetNumberOfPoints(samplePoints.size());
132  for(unsigned idx=0; idx< samplePoints.size(); idx++)
133  {
134  c_vector<double, DIM> location = samplePoints[idx].GetLocation(this->mpMesh->GetReferenceLengthScale());
135  if(DIM==3)
136  {
137  p_points->SetPoint(idx, location[0], location[1], location[2]);
138  }
139  else
140  {
141  p_points->SetPoint(idx, location[0], location[1], 0.0);
142  }
143  }
144  p_polydata->SetPoints(p_points);
145 
146  vtkSmartPointer<vtkProbeFilter> p_probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
147  #if VTK_MAJOR_VERSION <= 5
148  p_probe_filter->SetInput(p_polydata);
149  p_probe_filter->SetSource(this->mpVtkSolution);
150  #else
151  p_probe_filter->SetInputData(p_polydata);
152  p_probe_filter->SetSourceData(this->mpVtkSolution);
153  #endif
154  p_probe_filter->Update();
155  vtkSmartPointer<vtkPointData> p_point_data = p_probe_filter->GetOutput()->GetPointData();
156 
157  unsigned num_points = p_point_data->GetArray(this->mLabel.c_str())->GetNumberOfTuples();
158  for(unsigned idx=0; idx<num_points; idx++)
159  {
160  sampled_solution[idx] = p_point_data->GetArray(this->mLabel.c_str())->GetTuple1(idx)*this->mReferenceConcentration;
161  }
162  return sampled_solution;
163 }
164 
165 template<unsigned DIM>
166 std::vector<units::quantity<unit::concentration> > AbstractUnstructuredGridDiscreteContinuumSolver<DIM>::GetConcentrations(boost::shared_ptr<RegularGrid<DIM> > pGrid)
167 {
168  return this->GetConcentrations(pGrid->GetLocations());
169 }
170 
171 template<unsigned DIM>
172 std::vector<units::quantity<unit::concentration> > AbstractUnstructuredGridDiscreteContinuumSolver<DIM>::GetConcentrations(boost::shared_ptr<DiscreteContinuumMesh<DIM> > pMesh)
173 {
174  if(this->mpMesh == pMesh)
175  {
176  return this->GetConcentrations();
177  }
178  else
179  {
180  return this->GetConcentrations(pMesh->GetNodeLocationsAsPoints());
181  }
182 }
183 
184 template<unsigned DIM>
186 {
187  if(!this->mpVtkSolution)
188  {
189  this->Setup();
190  }
191 
192  std::vector<double> sampled_solution(rSamplePoints.size(), 0.0);
193 
194  // Sample the field at these locations
195  vtkSmartPointer<vtkPolyData> p_polydata = vtkSmartPointer<vtkPolyData>::New();
196  vtkSmartPointer<vtkPoints> p_points = vtkSmartPointer<vtkPoints>::New();
197  p_points->SetNumberOfPoints(rSamplePoints.size());
198  for(unsigned idx=0; idx< rSamplePoints.size(); idx++)
199  {
200  c_vector<double, DIM> location = rSamplePoints[idx].GetLocation(this->mpMesh->GetReferenceLengthScale());
201  if(DIM==3)
202  {
203  p_points->SetPoint(idx, location[0], location[1], location[2]);
204  }
205  else
206  {
207  p_points->SetPoint(idx, location[0], location[1], 0.0);
208  }
209  }
210  p_polydata->SetPoints(p_points);
211 
212  vtkSmartPointer<vtkProbeFilter> p_probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
213  #if VTK_MAJOR_VERSION <= 5
214  p_probe_filter->SetInput(p_polydata);
215  p_probe_filter->SetSource(this->mpVtkSolution);
216  #else
217  p_probe_filter->SetInputData(p_polydata);
218  p_probe_filter->SetSourceData(this->mpVtkSolution);
219  #endif
220  p_probe_filter->Update();
221  vtkSmartPointer<vtkPointData> p_point_data = p_probe_filter->GetOutput()->GetPointData();
222 
223  unsigned num_points = p_point_data->GetArray(this->mLabel.c_str())->GetNumberOfTuples();
224  for(unsigned idx=0; idx<num_points; idx++)
225  {
226  sampled_solution[idx] = p_point_data->GetArray(this->mLabel.c_str())->GetTuple1(idx);
227  }
228  return sampled_solution;
229 }
230 
231 template<unsigned DIM>
233 {
234  return this->GetSolution(pGrid->GetLocations());
235 }
236 
237 template<unsigned DIM>
239 {
240  if(this->mpMesh == pMesh)
241  {
242  return this->GetSolution();
243  }
244  else
245  {
246  return this->GetSolution(pMesh->GetNodeLocationsAsPoints());
247  }
248 }
249 
250 template<unsigned DIM>
252 {
253  if(!this->mpVtkSolution)
254  {
255  this->Setup();
256  }
257  return this->mpVtkSolution;
258 }
259 
260 template<unsigned DIM>
262 {
263  this->mpMesh = pMesh;
264 }
265 
266 template<unsigned DIM>
268 {
269  if(!this->mpMesh)
270  {
271  EXCEPTION("Mesh needed before Setup can be called.");
272  }
273 
274  // Set up the VTK solution
275  this->mpVtkSolution = mpMesh->GetAsVtkUnstructuredGrid();
276 
277  unsigned num_nodes = this->mpMesh->GetNodeLocationsAsPoints().size();
278  this->mSolution = std::vector<double>(0.0, num_nodes);
279 
280  if(this->mpNetwork)
281  {
282  this->mpMesh->SetVesselNetwork(this->mpNetwork) ;
283  }
284 }
285 
286 template<unsigned DIM>
288 {
289  if(!this->mpVtkSolution)
290  {
291  this->Setup();
292  }
293 
294  vtkSmartPointer<vtkDoubleArray> pPointData = vtkSmartPointer<vtkDoubleArray>::New();
295  pPointData->SetNumberOfComponents(1);
296  pPointData->SetNumberOfTuples(data.size());
297  pPointData->SetName(this->GetLabel().c_str());
298  for (unsigned i = 0; i < data.size(); i++)
299  {
300  pPointData->SetValue(i, data[i]);
301  }
302  this->mpVtkSolution->GetPointData()->AddArray(pPointData);
303 
304  // Note, if the data vector being passed in is mPointSolution, then it will be overwritten with zeros.
305  this->mSolution = std::vector<double>(data.size(), 0.0);
306  for (unsigned i = 0; i < data.size(); i++)
307  {
308  this->mSolution[i] = data[i];
309  }
310 }
311 
312 template<unsigned DIM>
314 {
315  if(!this->mpVtkSolution)
316  {
317  this->Setup();
318  }
319 
320  vtkSmartPointer<vtkDoubleArray> pPointData = vtkSmartPointer<vtkDoubleArray>::New();
321  pPointData->SetNumberOfComponents(1);
322  pPointData->SetNumberOfTuples(data.size());
323  pPointData->SetName(this->GetLabel().c_str());
324  for (unsigned i = 0; i < data.size(); i++)
325  {
326  pPointData->SetValue(i, data[i]);
327  }
328  this->mpVtkSolution->GetCellData()->AddArray(pPointData);
329 }
330 
331 template<unsigned DIM>
332 void AbstractUnstructuredGridDiscreteContinuumSolver<DIM>::UpdateSolution(const std::vector<units::quantity<unit::concentration> >& data)
333 {
334  if(!this->mpVtkSolution)
335  {
336  this->Setup();
337  }
338 
339  vtkSmartPointer<vtkDoubleArray> pPointData = vtkSmartPointer<vtkDoubleArray>::New();
340  pPointData->SetNumberOfComponents(1);
341  pPointData->SetNumberOfTuples(data.size());
342  pPointData->SetName(this->GetLabel().c_str());
343  for (unsigned i = 0; i < data.size(); i++)
344  {
345  pPointData->SetValue(i, data[i]/this->mReferenceConcentration);
346  }
347  this->mpVtkSolution->GetPointData()->AddArray(pPointData);
348 
349  // Note, if the data vector being passed in is mPointSolution, then it will be overwritten with zeros.
350  this->mConcentrations = std::vector<units::quantity<unit::concentration> >(data.size(), 0.0*this->mReferenceConcentration);
351  for (unsigned i = 0; i < data.size(); i++)
352  {
353  this->mConcentrations[i] = data[i];
354  }
355 }
356 
357 template<unsigned DIM>
359 {
360  if(!this->mpVtkSolution)
361  {
362  this->Setup();
363  }
364 
365  if(!this->CellPopulationIsSet())
366  {
367  EXCEPTION("The DiscreteContinuum solver needs a cell population for this operation.");
368  }
369 
370  // Update for FEM
371 // this->mpRegularGrid->SetCellPopulation(*(this->mpCellPopulation));
372 // std::vector<std::vector<CellPtr> > point_cell_map = this->mpRegularGrid->GetPointCellMap();
373 // for(unsigned idx=0; idx<point_cell_map.size(); idx++)
374 // {
375 // for(unsigned jdx=0; jdx<point_cell_map[idx].size(); jdx++)
376 // {
377 // point_cell_map[idx][jdx]->GetCellData()->SetItem(this->mLabel, this->mSolution[idx]);
378 // }
379 // }
380 }
381 
382 template<unsigned DIM>
384 {
385 
386 }
387 
388 template<unsigned DIM>
390 {
391  if(!this->mpVtkSolution)
392  {
393  this->Setup();
394  }
395 
396  if(!this->mpOutputFileHandler)
397  {
398  EXCEPTION("An output file handler has not been set for the DiscreteContinuum solver.");
399  }
400 
401  if(PetscTools::AmMaster())
402  {
403  vtkSmartPointer<vtkXMLUnstructuredGridWriter> p_writer1 = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
404  if(!this->mFilename.empty())
405  {
406  p_writer1->SetFileName((this->mpOutputFileHandler->GetOutputDirectoryFullPath() + "/" + this->mFilename+".vtu").c_str());
407  }
408  else
409  {
410  p_writer1->SetFileName((this->mpOutputFileHandler->GetOutputDirectoryFullPath() + "/solution.vtu").c_str());
411  }
412  #if VTK_MAJOR_VERSION <= 5
413  p_writer1->SetInput(this->mpVtkSolution);
414  #else
415  p_writer1->SetInputData(this->mpVtkSolution);
416  #endif
417  p_writer1->Update();
418  p_writer1->Write();
419  }
420 }
421 
422 // Explicit instantiation
virtual std::vector< units::quantity< unit::concentration > > GetConcentrations()
void SetMesh(boost::shared_ptr< DiscreteContinuumMesh< DIM, DIM > > pMesh)
virtual std::vector< units::quantity< unit::concentration > > GetConcentrationsAtCentroids()
boost::shared_ptr< OutputFileHandler > mpOutputFileHandler
std::vector< units::quantity< unit::concentration > > mConcentrations
boost::shared_ptr< VesselNetwork< DIM > > mpNetwork
units::quantity< unit::concentration > mReferenceConcentration
boost::shared_ptr< DiscreteContinuumMesh< DIM, DIM > > mpMesh