00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "CellBasedPdeHandler.hpp"
00030 #include "MeshBasedCellPopulationWithGhostNodes.hpp"
00031 #include "NodeBasedCellPopulation.hpp"
00032 #include "BoundaryConditionsContainer.hpp"
00033 #include "SimpleLinearEllipticSolver.hpp"
00034 #include "CellBasedPdeSolver.hpp"
00035 #include "CellwiseData.hpp"
00036 #include "Exception.hpp"
00037
00038 template<unsigned DIM>
00039 CellBasedPdeHandler<DIM>::CellBasedPdeHandler(AbstractCellPopulation<DIM>* pCellPopulation,
00040 bool deleteMemberPointersInDestructor)
00041 : mpCellPopulation(pCellPopulation),
00042 mWriteAverageRadialPdeSolution(false),
00043 mWriteDailyAverageRadialPdeSolution(false),
00044 mSetBcsOnCoarseBoundary(true),
00045 mNumRadialIntervals(UNSIGNED_UNSET),
00046 mpCoarsePdeMesh(NULL),
00047 mDeleteMemberPointersInDestructor(deleteMemberPointersInDestructor)
00048 {
00049
00051 assert((dynamic_cast<NodeBasedCellPopulation<DIM>*>(mpCellPopulation) != NULL) || (dynamic_cast<MeshBasedCellPopulation<DIM>*>(mpCellPopulation) != NULL));
00052 assert(dynamic_cast<MeshBasedCellPopulationWithGhostNodes<DIM>*>(mpCellPopulation) == NULL);
00053 assert(mpCellPopulation->GetNumRealCells() != 0);
00054 }
00055
00056 template<unsigned DIM>
00057 CellBasedPdeHandler<DIM>::~CellBasedPdeHandler()
00058 {
00059
00060
00061
00062
00063
00064
00065 if (mDeleteMemberPointersInDestructor)
00066 {
00067 for (unsigned i=0; i<mPdeAndBcCollection.size(); i++)
00068 {
00069 delete mPdeAndBcCollection[i];
00070 }
00071 }
00072 if (mpCoarsePdeMesh)
00073 {
00074 delete mpCoarsePdeMesh;
00075 }
00076 }
00077
00078 template<unsigned DIM>
00079 const AbstractCellPopulation<DIM>* CellBasedPdeHandler<DIM>::GetCellPopulation() const
00080 {
00081 return mpCellPopulation;
00082 }
00083
00084 template<unsigned DIM>
00085 TetrahedralMesh<DIM,DIM>* CellBasedPdeHandler<DIM>::GetCoarsePdeMesh()
00086 {
00087 return mpCoarsePdeMesh;
00088 }
00089
00090 template<unsigned DIM>
00091 void CellBasedPdeHandler<DIM>::AddPdeAndBc(PdeAndBoundaryConditions<DIM>* pPdeAndBc)
00092 {
00093 mPdeAndBcCollection.push_back(pPdeAndBc);
00094 }
00095
00096 template<unsigned DIM>
00097 Vec CellBasedPdeHandler<DIM>::GetPdeSolution(unsigned pdeIndex)
00098 {
00099 assert(pdeIndex<mPdeAndBcCollection.size());
00100 return mPdeAndBcCollection[pdeIndex]->GetSolution();
00101 }
00102
00103 template<unsigned DIM>
00104 void CellBasedPdeHandler<DIM>::InitialiseCellPdeElementMap()
00105 {
00106 if (mpCoarsePdeMesh == NULL)
00107 {
00108 EXCEPTION("InitialiseCellPdeElementMap() should only be called if mpCoarsePdeMesh is set up.");
00109 }
00110
00111 mCellPdeElementMap.clear();
00112
00113
00114 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00115 cell_iter != mpCellPopulation->End();
00116 ++cell_iter)
00117 {
00118 const ChastePoint<DIM>& r_position_of_cell = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
00119 unsigned elem_index = mpCoarsePdeMesh->GetContainingElementIndex(r_position_of_cell);
00120 mCellPdeElementMap[*cell_iter] = elem_index;
00121 }
00122 }
00123
00124 template<unsigned DIM>
00125 void CellBasedPdeHandler<DIM>::UpdateCellPdeElementMap()
00126 {
00127
00128 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00129 cell_iter != mpCellPopulation->End();
00130 ++cell_iter)
00131 {
00132 const ChastePoint<DIM>& r_position_of_cell = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
00133 unsigned elem_index = mpCoarsePdeMesh->GetContainingElementIndexWithInitialGuess(r_position_of_cell, mCellPdeElementMap[*cell_iter]);
00134 mCellPdeElementMap[*cell_iter] = elem_index;
00135 }
00136 }
00137
00138 template<unsigned DIM>
00139 void CellBasedPdeHandler<DIM>::OpenResultsFiles(std::string outputDirectory)
00140 {
00141
00142 if ((dynamic_cast<NodeBasedCellPopulation<DIM>*>(mpCellPopulation) != NULL) && mpCoarsePdeMesh==NULL)
00143 {
00144 EXCEPTION("Trying to solve a PDE on a NodeBasedCellPopulation without setting up a coarse mesh. Try calling UseCoarsePdeMesh().");
00145 }
00146
00147 if (mpCoarsePdeMesh != NULL)
00148 {
00149 InitialiseCellPdeElementMap();
00150
00151
00152 TrianglesMeshWriter<DIM,DIM> mesh_writer(outputDirectory+"/coarse_mesh_output", "coarse_mesh",false);
00153 mesh_writer.WriteFilesUsingMesh(*mpCoarsePdeMesh);
00154 }
00155
00156 if (PetscTools::AmMaster())
00157 {
00158 OutputFileHandler output_file_handler(outputDirectory+"/", false);
00159
00160 if (mpCoarsePdeMesh != NULL)
00161 {
00162 mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizcoarsepdesolution");
00163 }
00164 else
00165 {
00166 mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizpdesolution");
00167 }
00168
00169 if (mWriteAverageRadialPdeSolution)
00170 {
00171 mpAverageRadialPdeSolutionResultsFile = output_file_handler.OpenOutputFile("radial_dist.dat");
00172 }
00173 }
00174
00175 double current_time = SimulationTime::Instance()->GetTime();
00176 WritePdeSolution(current_time);
00177 }
00178
00179 template<unsigned DIM>
00180 void CellBasedPdeHandler<DIM>::CloseResultsFiles()
00181 {
00182
00183 if (PetscTools::AmMaster())
00184 {
00185 mpVizPdeSolutionResultsFile->close();
00186 if (mWriteAverageRadialPdeSolution)
00187 {
00188 WriteAverageRadialPdeSolution(SimulationTime::Instance()->GetTime());
00189 mpAverageRadialPdeSolutionResultsFile->close();
00190 }
00191 }
00192 }
00193
00194 template<unsigned DIM>
00195 void CellBasedPdeHandler<DIM>::UseCoarsePdeMesh(double stepSize, double meshWidth)
00196 {
00197
00198 if (mPdeAndBcCollection.empty())
00199 {
00200 EXCEPTION("mPdeAndBcCollection should be populated prior to calling UseCoarsePdeMesh().");
00201 }
00202 for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00203 {
00204 if (mPdeAndBcCollection[pde_index]->HasAveragedSourcePde() == false)
00205 {
00206 EXCEPTION("UseCoarsePdeMesh() should only be called if averaged-source PDEs are specified.");
00207 }
00208 }
00209
00210
00211 mpCoarsePdeMesh = new TetrahedralMesh<DIM,DIM>;
00212 switch (DIM)
00213 {
00214 case 1:
00215 mpCoarsePdeMesh->ConstructRegularSlabMesh(stepSize, meshWidth);
00216 break;
00217 case 2:
00218 mpCoarsePdeMesh->ConstructRegularSlabMesh(stepSize, meshWidth, meshWidth);
00219 break;
00220 case 3:
00221 mpCoarsePdeMesh->ConstructRegularSlabMesh(stepSize, meshWidth, meshWidth, meshWidth);
00222 break;
00223 default:
00224 NEVER_REACHED;
00225 }
00226
00227
00228 c_vector<double,DIM> centre_of_coarse_mesh = zero_vector<double>(DIM);
00229 for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
00230 {
00231 centre_of_coarse_mesh += mpCoarsePdeMesh->GetNode(i)->rGetLocation();
00232 }
00233 centre_of_coarse_mesh /= mpCoarsePdeMesh->GetNumNodes();
00234
00235
00236 c_vector<double,DIM> centre_of_cell_population = mpCellPopulation->GetCentroidOfCellPopulation();
00237 mpCoarsePdeMesh->Translate(centre_of_cell_population - centre_of_coarse_mesh);
00238 }
00239
00240 template<unsigned DIM>
00241 void CellBasedPdeHandler<DIM>::SolvePdeAndWriteResultsToFile(unsigned samplingTimestepMultiple)
00242 {
00243
00244 bool using_coarse_pde_mesh = (mpCoarsePdeMesh != NULL);
00245
00246
00247 assert(!mPdeAndBcCollection.empty());
00248 for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00249 {
00250 assert(mPdeAndBcCollection[pde_index]);
00251 assert(mPdeAndBcCollection[pde_index]->HasAveragedSourcePde() == using_coarse_pde_mesh);
00252 }
00253
00254
00255 CellwiseData<DIM>::Instance()->ReallocateMemory();
00256
00257
00258 mpCellPopulation->Update();
00259
00260
00261 TetrahedralMesh<DIM,DIM>* p_mesh;
00262 if (using_coarse_pde_mesh)
00263 {
00264 p_mesh = mpCoarsePdeMesh;
00265 }
00266 else
00267 {
00268
00269 p_mesh = &(static_cast<MeshBasedCellPopulation<DIM>*>(mpCellPopulation)->rGetMesh());
00270 }
00271
00272
00273 for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00274 {
00275
00276 PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];
00277
00278
00279 AbstractBoundaryCondition<DIM>* p_bc = p_pde_and_bc->GetBoundaryCondition();
00280 BoundaryConditionsContainer<DIM,DIM,1> bcc(false);
00281
00282 if (p_pde_and_bc->IsNeumannBoundaryCondition())
00283 {
00284 if (using_coarse_pde_mesh)
00285 {
00287 EXCEPTION("Neumann BCs not yet implemented when using a coarse PDE mesh");
00288 }
00289 else
00290 {
00291 for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator elem_iter = p_mesh->GetBoundaryElementIteratorBegin();
00292 elem_iter != p_mesh->GetBoundaryElementIteratorEnd();
00293 ++elem_iter)
00294 {
00295 bcc.AddNeumannBoundaryCondition(*elem_iter, p_bc);
00296 }
00297 }
00298 }
00299 else
00300 {
00301 if (using_coarse_pde_mesh && !mSetBcsOnCoarseBoundary)
00302 {
00303
00304 std::set<unsigned> coarse_element_indices_in_map;
00305 for (typename AbstractCentreBasedCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00306 cell_iter != mpCellPopulation->End();
00307 ++cell_iter)
00308 {
00309 coarse_element_indices_in_map.insert(mCellPdeElementMap[*cell_iter]);
00310 }
00311
00312
00313 std::set<unsigned> coarse_mesh_boundary_node_indices;
00314 for (unsigned i=0; i<p_mesh->GetNumElements(); i++)
00315 {
00316 if (coarse_element_indices_in_map.find(i) == coarse_element_indices_in_map.end())
00317 {
00318 Element<DIM,DIM>* p_element = p_mesh->GetElement(i);
00319 for (unsigned j=0; j<DIM+1; j++)
00320 {
00321 unsigned node_index = p_element->GetNodeGlobalIndex(j);
00322 coarse_mesh_boundary_node_indices.insert(node_index);
00323 }
00324 }
00325 }
00326
00327
00328 for (std::set<unsigned>::iterator iter = coarse_mesh_boundary_node_indices.begin();
00329 iter != coarse_mesh_boundary_node_indices.end();
00330 ++iter)
00331 {
00332 bcc.AddDirichletBoundaryCondition(p_mesh->GetNode(*iter), p_bc, 0, false);
00333 }
00334 }
00335 else
00336 {
00337 for (typename TetrahedralMesh<DIM,DIM>::BoundaryNodeIterator node_iter = p_mesh->GetBoundaryNodeIteratorBegin();
00338 node_iter != p_mesh->GetBoundaryNodeIteratorEnd();
00339 ++node_iter)
00340 {
00341 bcc.AddDirichletBoundaryCondition(*node_iter, p_bc);
00342 }
00343 }
00344 }
00345
00346
00347 PetscInt previous_solution_size = 0;
00348 if (p_pde_and_bc->GetSolution())
00349 {
00350 VecGetSize(p_pde_and_bc->GetSolution(), &previous_solution_size);
00351 }
00352
00353
00354 bool is_previous_solution_size_correct = (previous_solution_size == (int)p_mesh->GetNumNodes());
00355
00356
00357 Vec initial_guess;
00358 if (is_previous_solution_size_correct)
00359 {
00360
00361 VecDuplicate(p_pde_and_bc->GetSolution(), &initial_guess);
00362 VecCopy(p_pde_and_bc->GetSolution(), initial_guess);
00363 p_pde_and_bc->DestroySolution();
00364 }
00365 else
00366 {
00368 if (!using_coarse_pde_mesh && p_pde_and_bc->GetSolution())
00369 {
00370 assert(previous_solution_size != 0);
00371 p_pde_and_bc->DestroySolution();
00372 }
00373 }
00374
00375
00376 if (using_coarse_pde_mesh)
00377 {
00378
00379
00380 this->UpdateCellPdeElementMap();
00381 p_pde_and_bc->SetUpSourceTermsForAveragedSourcePde(p_mesh, &mCellPdeElementMap);
00382
00383 SimpleLinearEllipticSolver<DIM,DIM> solver(p_mesh, p_pde_and_bc->GetPde(), &bcc);
00384
00385
00386 if (is_previous_solution_size_correct)
00387 {
00388 p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
00389 VecDestroy(initial_guess);
00390 }
00391 else
00392 {
00393 p_pde_and_bc->SetSolution(solver.Solve());
00394 }
00395 }
00396 else
00397 {
00398 CellBasedPdeSolver<DIM> solver(p_mesh, p_pde_and_bc->GetPde(), &bcc);
00399
00400
00401 if (is_previous_solution_size_correct)
00402 {
00403 p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
00404 VecDestroy(initial_guess);
00405 }
00406 else
00407 {
00408 p_pde_and_bc->SetSolution(solver.Solve());
00409 }
00410 }
00411
00412
00413 ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
00414
00415
00416 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00417 cell_iter != mpCellPopulation->End();
00418 ++cell_iter)
00419 {
00420 unsigned node_index = mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
00421 double solution_at_node = 0.0;
00422
00423 if (using_coarse_pde_mesh)
00424 {
00425
00426
00427
00428 unsigned elem_index = mCellPdeElementMap[*cell_iter];
00429 Element<DIM,DIM>* p_element = mpCoarsePdeMesh->GetElement(elem_index);
00430
00431 const ChastePoint<DIM>& node_location = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
00432
00433 c_vector<double,DIM+1> weights = p_element->CalculateInterpolationWeights(node_location);
00434 for (unsigned i=0; i<DIM+1; i++)
00435 {
00436 double nodal_value = solution_repl[p_element->GetNodeGlobalIndex(i)];
00437 solution_at_node += nodal_value * weights(i);
00438 }
00439 }
00440 else
00441 {
00442 solution_at_node = solution_repl[node_index];
00443 }
00444
00445 CellwiseData<DIM>::Instance()->SetValue(solution_at_node, node_index, pde_index);
00446 }
00447 }
00448
00449
00450 SimulationTime* p_time = SimulationTime::Instance();
00451 double time_next_step = p_time->GetTime() + p_time->GetTimeStep();
00452 if ((p_time->GetTimeStepsElapsed()+1)%samplingTimestepMultiple == 0)
00453 {
00454 WritePdeSolution(time_next_step);
00455 }
00456
00457 #define COVERAGE_IGNORE
00458
00459 if (!using_coarse_pde_mesh)
00460 {
00461 if (mWriteDailyAverageRadialPdeSolution)
00462 {
00464 SimulationTime* p_time = SimulationTime::Instance();
00465 double time_next_step = p_time->GetTime() + p_time->GetTimeStep();
00466 unsigned num_timesteps_per_day = (unsigned) (DBL_EPSILON + 24/SimulationTime::Instance()->GetTimeStep());
00467 if ((p_time->GetTimeStepsElapsed()+1) % num_timesteps_per_day == 0)
00468 {
00469 WriteAverageRadialPdeSolution(time_next_step);
00470 }
00471 }
00472 }
00473 #undef COVERAGE_IGNORE
00474 }
00475
00476 template<unsigned DIM>
00477 unsigned CellBasedPdeHandler<DIM>::FindCoarseElementContainingCell(CellPtr pCell)
00478 {
00479
00480 unsigned old_element_index = mCellPdeElementMap[pCell];
00481
00482
00483 std::set<unsigned> test_elements;
00484 test_elements.insert(old_element_index);
00485
00486 Element<DIM,DIM>* p_element = mpCoarsePdeMesh->GetElement(old_element_index);
00487 for (unsigned local_index=0; local_index<DIM+1; local_index++)
00488 {
00489 std::set<unsigned> element_indices = p_element->GetNode(local_index)->rGetContainingElementIndices();
00490 for (std::set<unsigned>::iterator iter = element_indices.begin();
00491 iter != element_indices.end();
00492 ++iter)
00493 {
00494 test_elements.insert(*iter);
00495 }
00496 }
00497
00498
00499 const ChastePoint<DIM>& r_cell_position = mpCellPopulation->GetLocationOfCellCentre(pCell);
00500 unsigned new_element_index = mpCoarsePdeMesh->GetContainingElementIndex(r_cell_position, false, test_elements);
00501
00502
00503 mCellPdeElementMap[pCell] = new_element_index;
00504
00505 return new_element_index;
00506 }
00507
00508 template<unsigned DIM>
00509 void CellBasedPdeHandler<DIM>::WritePdeSolution(double time)
00510 {
00511 if (PetscTools::AmMaster())
00512 {
00513 (*mpVizPdeSolutionResultsFile) << time << "\t";
00514
00515 for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00516 {
00517 if (mpCoarsePdeMesh != NULL)
00518 {
00519 PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];
00520
00521 for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
00522 {
00523 (*mpVizPdeSolutionResultsFile) << i << " ";
00524 c_vector<double,DIM> location = mpCoarsePdeMesh->GetNode(i)->rGetLocation();
00525 for (unsigned k=0; k<DIM; k++)
00526 {
00527 (*mpVizPdeSolutionResultsFile) << location[k] << " ";
00528 }
00529
00530 if (p_pde_and_bc->GetSolution())
00531 {
00532 ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
00533 (*mpVizPdeSolutionResultsFile) << solution_repl[i] << " ";
00534 }
00535 else
00536 {
00538
00539
00540 unsigned nearest_node_index = 0;
00541 double nearest_node_distance = DBL_MAX;
00542 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00543 cell_iter != mpCellPopulation->End();
00544 ++cell_iter)
00545 {
00546 c_vector<double, DIM> node_location = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
00547 if (norm_2(node_location - location) < nearest_node_distance)
00548 {
00549 nearest_node_index = mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
00550 }
00551 }
00552
00553 CellPtr p_cell = mpCellPopulation->GetCellUsingLocationIndex(nearest_node_index);
00554 double solution = CellwiseData<DIM>::Instance()->GetValue(p_cell, pde_index);
00555 (*mpVizPdeSolutionResultsFile) << solution << " ";
00556 }
00557 }
00558 }
00559 else
00560 {
00561 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00562 cell_iter != mpCellPopulation->End();
00563 ++cell_iter)
00564 {
00565 unsigned node_index = mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
00566 (*mpVizPdeSolutionResultsFile) << node_index << " ";
00567 const c_vector<double,DIM>& position = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
00568 for (unsigned i=0; i<DIM; i++)
00569 {
00570 (*mpVizPdeSolutionResultsFile) << position[i] << " ";
00571 }
00572 double solution = CellwiseData<DIM>::Instance()->GetValue(*cell_iter, pde_index);
00573 (*mpVizPdeSolutionResultsFile) << solution << " ";
00574 }
00575 }
00576 }
00577 (*mpVizPdeSolutionResultsFile) << "\n";
00578 }
00579 }
00580
00581 template<unsigned DIM>
00582 void CellBasedPdeHandler<DIM>::SetWriteAverageRadialPdeSolution(unsigned numRadialIntervals, bool writeDailyResults)
00583 {
00584 mWriteAverageRadialPdeSolution = true;
00585 mNumRadialIntervals = numRadialIntervals;
00586 mWriteDailyAverageRadialPdeSolution = writeDailyResults;
00587 }
00588
00589 template<unsigned DIM>
00590 void CellBasedPdeHandler<DIM>::SetImposeBcsOnCoarseBoundary(bool setBcsOnCoarseBoundary)
00591 {
00592 mSetBcsOnCoarseBoundary = setBcsOnCoarseBoundary;
00593 }
00594
00595 template<unsigned DIM>
00596 void CellBasedPdeHandler<DIM>::WriteAverageRadialPdeSolution(double time)
00597 {
00598 (*mpAverageRadialPdeSolutionResultsFile) << time << " ";
00599
00600
00601 c_vector<double,DIM> centre = mpCellPopulation->GetCentroidOfCellPopulation();
00602
00603
00604 std::map<double, CellPtr> radius_cell_map;
00605 double max_distance_from_centre = 0.0;
00606 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
00607 cell_iter != mpCellPopulation->End();
00608 ++cell_iter)
00609 {
00610 double distance = norm_2(mpCellPopulation->GetLocationOfCellCentre(*cell_iter) - centre);
00611 radius_cell_map[distance] = *cell_iter;
00612
00613 if (distance > max_distance_from_centre)
00614 {
00615 max_distance_from_centre = distance;
00616 }
00617 }
00618
00619
00620 std::vector<double> radius_intervals;
00621 for (unsigned i=0; i<mNumRadialIntervals; i++)
00622 {
00623 double upper_radius = max_distance_from_centre*((double) i+1)/((double) mNumRadialIntervals);
00624 radius_intervals.push_back(upper_radius);
00625 }
00626
00627
00628 double lower_radius = 0.0;
00629 for (unsigned i=0; i<mNumRadialIntervals; i++)
00630 {
00631 unsigned counter = 0;
00632 double average_solution = 0.0;
00633
00634 for (std::map<double, CellPtr>::iterator iter = radius_cell_map.begin(); iter != radius_cell_map.end(); ++iter)
00635 {
00636 if (iter->first > lower_radius && iter->first <= radius_intervals[i])
00637 {
00638 average_solution += CellwiseData<DIM>::Instance()->GetValue(iter->second);
00639 counter++;
00640 }
00641 }
00642 if (counter != 0)
00643 {
00644 average_solution /= (double) counter;
00645 }
00646
00647
00648 (*mpAverageRadialPdeSolutionResultsFile) << radius_intervals[i] << " " << average_solution << " ";
00649 lower_radius = radius_intervals[i];
00650 }
00651 (*mpAverageRadialPdeSolutionResultsFile) << "\n";
00652 }
00653
00654 template<unsigned DIM>
00655 bool CellBasedPdeHandler<DIM>::GetWriteAverageRadialPdeSolution()
00656 {
00657 return mWriteAverageRadialPdeSolution;
00658 }
00659
00660 template<unsigned DIM>
00661 bool CellBasedPdeHandler<DIM>::GetWriteDailyAverageRadialPdeSolution()
00662 {
00663 return mWriteDailyAverageRadialPdeSolution;
00664 }
00665
00666 template<unsigned DIM>
00667 bool CellBasedPdeHandler<DIM>::GetImposeBcsOnCoarseBoundary()
00668 {
00669 return mSetBcsOnCoarseBoundary;
00670 }
00671
00672 template<unsigned DIM>
00673 unsigned CellBasedPdeHandler<DIM>::GetNumRadialIntervals()
00674 {
00675 return mNumRadialIntervals;
00676 }
00677
00678 template<unsigned DIM>
00679 void CellBasedPdeHandler<DIM>::OutputParameters(out_stream& rParamsFile)
00680 {
00681 std::string type = GetIdentifier();
00682
00683 *rParamsFile << "\t\t<" << type << ">\n";
00684 *rParamsFile << "\t\t<WriteAverageRadialPdeSolution>" << mWriteAverageRadialPdeSolution << "</WriteAverageRadialPdeSolution>\n";
00685 *rParamsFile << "\t\t<WriteDailyAverageRadialPdeSolution>" << mWriteDailyAverageRadialPdeSolution << "</WriteDailyAverageRadialPdeSolution>\n";
00686 *rParamsFile << "\t\t<SetBcsOnCoarseBoundary>" << mSetBcsOnCoarseBoundary << "</SetBcsOnCoarseBoundary>\n";
00687 *rParamsFile << "\t\t<NumRadialIntervals>" << mNumRadialIntervals << "</NumRadialIntervals>\n";
00688 *rParamsFile << "\t\t</" << type << ">\n";
00689 }
00690
00691
00692 #include "SerializationExportWrapperForCpp.hpp"
00693 EXPORT_TEMPLATE_CLASS_SAME_DIMS(CellBasedPdeHandler)
00694
00695
00696
00698
00699 template class CellBasedPdeHandler<1>;
00700 template class CellBasedPdeHandler<2>;
00701 template class CellBasedPdeHandler<3>;