From 0f5161aeb73a9845f952767d5860f860c28309e0 Mon Sep 17 00:00:00 2001 From: James Osborne Date: Tue, 7 Jan 2025 03:31:30 +0000 Subject: [PATCH] Adding Potts PDE --- .../AbstractPottsPdeModifier.cpp | 621 ++++++++++++++++++ .../AbstractPottsPdeModifier.hpp | 299 +++++++++ .../ParabolicPottsPdeModifier.cpp | 413 ++++++++++++ .../ParabolicPottsPdeModifier.hpp | 237 +++++++ src/Test06Angiogenesis/PottsParabolicPde.cpp | 223 +++++++ src/Test06Angiogenesis/PottsParabolicPde.hpp | 282 ++++++++ .../VegfChemotaxisPottsUpdateRule.cpp | 110 ++++ .../VegfChemotaxisPottsUpdateRule.hpp | 111 ++++ test/Test06Angiogenesis.hpp | 111 ++-- 9 files changed, 2363 insertions(+), 44 deletions(-) create mode 100644 src/Test06Angiogenesis/AbstractPottsPdeModifier.cpp create mode 100644 src/Test06Angiogenesis/AbstractPottsPdeModifier.hpp create mode 100644 src/Test06Angiogenesis/ParabolicPottsPdeModifier.cpp create mode 100644 src/Test06Angiogenesis/ParabolicPottsPdeModifier.hpp create mode 100644 src/Test06Angiogenesis/PottsParabolicPde.cpp create mode 100644 src/Test06Angiogenesis/PottsParabolicPde.hpp create mode 100644 src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.cpp create mode 100644 src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.hpp diff --git a/src/Test06Angiogenesis/AbstractPottsPdeModifier.cpp b/src/Test06Angiogenesis/AbstractPottsPdeModifier.cpp new file mode 100644 index 0000000..44ec729 --- /dev/null +++ b/src/Test06Angiogenesis/AbstractPottsPdeModifier.cpp @@ -0,0 +1,621 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "AbstractPottsPdeModifier.hpp" +#include "ReplicatableVector.hpp" +#include "LinearBasisFunction.hpp" +#include "Debug.hpp" +#include "VtkMeshWriter.hpp" + +template +AbstractPottsPdeModifier::AbstractPottsPdeModifier(boost::shared_ptr > pPde, + boost::shared_ptr > pBoundaryCondition, + bool isNeumannBoundaryCondition, + boost::shared_ptr > pMeshCuboid, + double stepSize, + Vec solution) + : AbstractPdeModifier(pPde, + pBoundaryCondition, + isNeumannBoundaryCondition, + solution), + mpMeshCuboid(pMeshCuboid), + mStepSize(stepSize), + mSetBcsOnBoxBoundary(true), + mSetBcsOnBoundingSphere(false), + mSetBcsOnConvexHull(false), + mUseVoronoiCellsForInterpolation(false), + mTypicalCellRadius(0.5) // defaults to 0.5 +{ + if (pMeshCuboid) + { + // We only need to generate mpFeMesh once, as it does not vary with time + this->GenerateFeMesh(mpMeshCuboid, mStepSize); + this->mDeleteFeMesh = true; + //initialise the boundary nodes + this->mIsDirichletBoundaryNode = std::vector(this->mpFeMesh->GetNumNodes(), 0.0); + } +} + +template +AbstractPottsPdeModifier::~AbstractPottsPdeModifier() +{ +} + +template +double AbstractPottsPdeModifier::GetStepSize() +{ + return mStepSize; +} + +template +void AbstractPottsPdeModifier::SetBcsOnBoxBoundary(bool setBcsOnBoxBoundary) +{ + mSetBcsOnBoxBoundary = setBcsOnBoxBoundary; +} + +template +bool AbstractPottsPdeModifier::AreBcsSetOnBoxBoundary() +{ + return mSetBcsOnBoxBoundary; +} + +template +void AbstractPottsPdeModifier::SetBcsOnBoundingSphere(bool setBcsOnBoundingSphere) +{ + mSetBcsOnBoundingSphere = setBcsOnBoundingSphere; +} + +template +bool AbstractPottsPdeModifier::AreBcsSetOnBoundingSphere() +{ + return mSetBcsOnBoundingSphere; +} + +template +void AbstractPottsPdeModifier::SetBcsOnConvexHull(bool setBcsOnConvexHull) +{ + mSetBcsOnConvexHull = setBcsOnConvexHull; +} + +template +bool AbstractPottsPdeModifier::AreBcsSetOnConvexHull() +{ + return mSetBcsOnConvexHull; +} + +template +void AbstractPottsPdeModifier::SetUseVoronoiCellsForInterpolation(bool useVoronoiCellsForInterpolation) +{ + mUseVoronoiCellsForInterpolation = useVoronoiCellsForInterpolation; +} + +template +bool AbstractPottsPdeModifier::GetUseVoronoiCellsForInterpolation() +{ + return mUseVoronoiCellsForInterpolation; +} + +template +void AbstractPottsPdeModifier::SetTypicalCellRadius(double typicalCellRadius) +{ + assert(mTypicalCellRadius>=0.0); + mTypicalCellRadius = typicalCellRadius; +} + +template +double AbstractPottsPdeModifier::GetTypicalCellRadius() +{ + return mTypicalCellRadius; +} + +template +void AbstractPottsPdeModifier::ConstructBoundaryConditionsContainerHelper(AbstractCellPopulation& rCellPopulation, + std::shared_ptr > pBcc) +{ + // Reset mIsDirichletBoundaryNode + for (unsigned i=0; impFeMesh->GetNumNodes(); i++) + { + this->mIsDirichletBoundaryNode[i] = 0.0; + } + + if (!this->mSetBcsOnBoxBoundary) + { + // Here we approximate the cell population by the bounding spehere and apply the boundary conditions outside the sphere. + if (this->mSetBcsOnBoundingSphere) + { + //Cant apply on the convex hull and bounding sphere at same time + if(this->mSetBcsOnConvexHull) + { + NEVER_REACHED; + } + + // First find the centre of the tissue by choosing the mid point of the extrema. + c_vector tissue_maxima = zero_vector(DIM); + c_vector tissue_minima = zero_vector(DIM); + for (unsigned i = 0; i < DIM; i++) + { + tissue_maxima[i] = -DBL_MAX; + tissue_minima[i] = DBL_MAX; + } + + + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + for (unsigned i = 0; i < DIM; i++) + { + if (r_position_of_cell[i] > tissue_maxima[i]) + { + tissue_maxima[i] = r_position_of_cell[i]; + } + if (r_position_of_cell[i] < tissue_minima[i]) + { + tissue_minima[i] = r_position_of_cell[i]; + } + } + } + + c_vector tissue_centre = 0.5*(tissue_maxima + tissue_minima); + + + double tissue_radius = 0.0; + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + double radius = norm_2(tissue_centre - r_position_of_cell.rGetLocation()); + + if (tissue_radius < radius) + { + tissue_radius = radius; + } + } + + // Apply boundary condition to the nodes outside the tissue_radius + if (this->IsNeumannBoundaryCondition()) + { + // Neumann BSC not implemented yet for this geometry + NEVER_REACHED; + } + else + { + // Impose any Dirichlet boundary conditions + for (unsigned i=0; impFeMesh->GetNumNodes(); i++) + { + double radius = norm_2(tissue_centre - this->mpFeMesh->GetNode(i)->rGetLocation()); + + if (radius > tissue_radius) + { + pBcc->AddDirichletBoundaryCondition(this->mpFeMesh->GetNode(i), this->mpBoundaryCondition.get(), 0, false); + this->mIsDirichletBoundaryNode[i] = 1.0; + } + } + } + } + // Here apply the dirichlet conditions to nodes outside the convex hull of the cell centres + // (i.e its applied outside the mesh thats formed by the cetres) + else if (this->mSetBcsOnConvexHull) + { + // Can't apply on the convex hull and bounding sphere at same time + if(this->mSetBcsOnBoundingSphere) + { + NEVER_REACHED; + } + + // Apply boundary condition to the nodes outside the convex hull + if (this->IsNeumannBoundaryCondition()) + { + // Neumann BSC not implemented yet for this geometry + NEVER_REACHED; + } + else + { + TetrahedralMesh* p_convex_mesh = rCellPopulation.GetTetrahedralMeshForPdeModifier(); + + // Impose any Dirichlet boundary conditions + for (unsigned i=0; impFeMesh->GetNumNodes(); i++) + { + std::vector containing_element_indices = p_convex_mesh->GetContainingElementIndices(this->mpFeMesh->GetNode(i)->rGetLocation()); + + if (containing_element_indices.size()==0u) + { + pBcc->AddDirichletBoundaryCondition(this->mpFeMesh->GetNode(i), this->mpBoundaryCondition.get(), 0, false); + this->mIsDirichletBoundaryNode[i] = 1.0; + } + } + } + } + else // Set pde nodes as boundary node if elements dont contain cells or nodes aren't within a specified distance of a cell centre + { + // Get the set of coarse element indices that contain cells + std::set coarse_element_indices_in_map; + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + coarse_element_indices_in_map.insert(this->mCellPdeElementMap[*cell_iter]); + } + + // Find the node indices associated with elements whose indices are NOT in the set coarse_element_indices_in_map + std::set coarse_mesh_boundary_node_indices; + for (unsigned i=0; impFeMesh->GetNumElements(); i++) + { + if (coarse_element_indices_in_map.find(i) == coarse_element_indices_in_map.end()) + { + Element* p_element = this->mpFeMesh->GetElement(i); + for (unsigned j=0; jGetNodeGlobalIndex(j); + coarse_mesh_boundary_node_indices.insert(node_index); + } + } + } + + // Also remove nodes that are within the typical cell radius from the centre of a cell. + std::set nearby_node_indices; + for (std::set::iterator node_iter = coarse_mesh_boundary_node_indices.begin(); + node_iter != coarse_mesh_boundary_node_indices.end(); + ++node_iter) + { + bool remove_node = false; + + c_vector node_location = this->mpFeMesh->GetNode(*node_iter)->rGetLocation(); + + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + double separation = norm_2(node_location - r_position_of_cell.rGetLocation()); + + if (separation <= mTypicalCellRadius) + { + remove_node = true; + break; + } + } + + if (remove_node) + { + // Node near cell so set it to be removed from boundary set + nearby_node_indices.insert(*node_iter); + } + } + + // Remove nodes that are near cells from boundary set + for (std::set::iterator node_iter = nearby_node_indices.begin(); + node_iter != nearby_node_indices.end(); + ++node_iter) + { + coarse_mesh_boundary_node_indices.erase(*node_iter); + } + + + // Apply boundary condition to the nodes in the set coarse_mesh_boundary_node_indices + if (this->IsNeumannBoundaryCondition()) + { + // Neumann BSC not implemented yet for this geometry + NEVER_REACHED; + } + else + { + // Impose any Dirichlet boundary conditions + for (std::set::iterator iter = coarse_mesh_boundary_node_indices.begin(); + iter != coarse_mesh_boundary_node_indices.end(); + ++iter) + { + pBcc->AddDirichletBoundaryCondition(this->mpFeMesh->GetNode(*iter), this->mpBoundaryCondition.get(), 0, false); + this->mIsDirichletBoundaryNode[*iter] = 1.0; + } + } + } + } + else // Apply BC at boundary of box domain FE mesh + { + if (this->IsNeumannBoundaryCondition()) + { + // Impose any Neumann boundary conditions + for (typename TetrahedralMesh::BoundaryElementIterator elem_iter = this->mpFeMesh->GetBoundaryElementIteratorBegin(); + elem_iter != this->mpFeMesh->GetBoundaryElementIteratorEnd(); + ++elem_iter) + { + pBcc->AddNeumannBoundaryCondition(*elem_iter, this->mpBoundaryCondition.get()); + } + } + else + { + // Impose any Dirichlet boundary conditions + for (typename TetrahedralMesh::BoundaryNodeIterator node_iter = this->mpFeMesh->GetBoundaryNodeIteratorBegin(); + node_iter != this->mpFeMesh->GetBoundaryNodeIteratorEnd(); + ++node_iter) + { + pBcc->AddDirichletBoundaryCondition(*node_iter, this->mpBoundaryCondition.get()); + this->mIsDirichletBoundaryNode[(*node_iter)->GetIndex()] = 1.0; + } + } + } +} + + + +template +void AbstractPottsPdeModifier::SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory) +{ + AbstractPdeModifier::SetupSolve(rCellPopulation, outputDirectory); + + InitialiseCellPdeElementMap(rCellPopulation); +} + +template +void AbstractPottsPdeModifier::GenerateFeMesh(boost::shared_ptr > pMeshCuboid, double stepSize) +{ + // Create a regular coarse tetrahedral mesh + this->mpFeMesh = new TetrahedralMesh(); + + GenerateAndReturnFeMesh(pMeshCuboid, stepSize, this->mpFeMesh); +} + +template +void AbstractPottsPdeModifier::GenerateAndReturnFeMesh(boost::shared_ptr > pMeshCuboid, double stepSize, TetrahedralMesh* pMesh) +{ + // Create a regular coarse tetrahedral mesh + switch (DIM) + { + case 1: + pMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0)); + break; + case 2: + pMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0), pMeshCuboid->GetWidth(1)); + break; + case 3: + pMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0), pMeshCuboid->GetWidth(1), pMeshCuboid->GetWidth(2)); + break; + default: + NEVER_REACHED; + } + + // Get centroid of meshCuboid + ChastePoint upper = pMeshCuboid->rGetUpperCorner(); + ChastePoint lower = pMeshCuboid->rGetLowerCorner(); + c_vector centre_of_cuboid = 0.5*(upper.rGetLocation() + lower.rGetLocation()); + + // Find the centre of the PDE mesh + c_vector centre_of_coarse_mesh = zero_vector(DIM); + for (unsigned i=0; iGetNumNodes(); i++) + { + centre_of_coarse_mesh += pMesh->GetNode(i)->rGetLocation(); + } + centre_of_coarse_mesh /= pMesh->GetNumNodes(); + + // Now move the mesh to the correct location + pMesh->Translate(centre_of_cuboid - centre_of_coarse_mesh); +} + +template +void AbstractPottsPdeModifier::UpdateCellData(AbstractCellPopulation& rCellPopulation) +{ + // Store the PDE solution in an accessible form + ReplicatableVector solution_repl(this->mSolution); + + if (mUseVoronoiCellsForInterpolation) + { + NEVER_REACHED; + unsigned num_nodes = rCellPopulation.GetNumNodes(); + + std::vector cell_data(num_nodes, -1); + std::vector num_cells(num_nodes, -1); + + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + unsigned cell_location_index = rCellPopulation.GetLocationIndexUsingCell(*cell_iter); + cell_data[cell_location_index]=0.0; + num_cells[cell_location_index]=0; + } + + // Loop over nodes of the finite element mesh and work out which voronoi region the node is in. + for (typename TetrahedralMesh::NodeIterator node_iter = this->mpFeMesh->GetNodeIteratorBegin(); + node_iter != this->mpFeMesh->GetNodeIteratorEnd(); + ++node_iter) + { + unsigned node_index = node_iter->GetIndex(); + + c_vector node_location = node_iter->rGetLocation(); + + double closest_separation = DBL_MAX; + unsigned nearest_cell = UNSIGNED_UNSET; + + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + c_vector cell_location = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + double separation = norm_2(node_location - cell_location); + + if (separation < closest_separation) + { + closest_separation = separation; + nearest_cell = rCellPopulation.GetLocationIndexUsingCell(*cell_iter); + } + } + assert(closest_separation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + unsigned cell_location_index = rCellPopulation.GetLocationIndexUsingCell(*cell_iter); + + + if (num_cells[cell_location_index]==0) + { + EXCEPTION("One or more of the cells doesnt contain any pde nodes so cant use voroni CellData calculation in the "); + } + + double solution_at_cell = cell_data[cell_location_index]/num_cells[cell_location_index]; + + cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_cell); + } + + if (this->mOutputGradient) + { + // This isnt implemented yet + NEVER_REACHED; + } + } + else // Interpolate solutions + { + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + // The cells are not nodes of the mesh, so we must interpolate + double solution_at_cell = 0.0; + + // Find the element in the FE mesh that contains this cell. CellElementMap has been updated so use this. + unsigned elem_index = mCellPdeElementMap[*cell_iter]; + Element* p_element = this->mpFeMesh->GetElement(elem_index); + + const ChastePoint& node_location = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + c_vector weights = p_element->CalculateInterpolationWeights(node_location); + + for (unsigned i=0; iGetNodeGlobalIndex(i)]; + solution_at_cell += nodal_value * weights(i); + } + + cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_cell); + + if (this->mOutputGradient) + { + // Now calculate the gradient of the solution and store this in CellVecData + c_vector solution_gradient = zero_vector(DIM); + + // Calculate the basis functions at any point (e.g. zero) in the element + c_matrix jacobian, inverse_jacobian; + double jacobian_det; + this->mpFeMesh->GetInverseJacobianForElement(elem_index, jacobian, jacobian_det, inverse_jacobian); + const ChastePoint zero_point; + c_matrix grad_phi; + LinearBasisFunction::ComputeTransformedBasisFunctionDerivatives(zero_point, inverse_jacobian, grad_phi); + + for (unsigned node_index=0; node_indexGetNodeGlobalIndex(node_index)]; + + for (unsigned j=0; jGetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0)); + break; + case 2: + cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0)); + cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1)); + break; + case 3: + cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0)); + cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1)); + cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_z", solution_gradient(2)); + break; + default: + NEVER_REACHED; + } + } + } + } +} + +template +void AbstractPottsPdeModifier::InitialiseCellPdeElementMap(AbstractCellPopulation& rCellPopulation) +{ + mCellPdeElementMap.clear(); + + // Find the element of mpFeMesh that contains each cell and populate mCellPdeElementMap + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + unsigned elem_index = this->mpFeMesh->GetContainingElementIndex(r_position_of_cell); + mCellPdeElementMap[*cell_iter] = elem_index; + } +} + +template +void AbstractPottsPdeModifier::UpdateCellPdeElementMap(AbstractCellPopulation& rCellPopulation) +{ + // Find the element of mpCoarsePdeMesh that contains each cell and populate mCellPdeElementMap + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + unsigned elem_index = this->mpFeMesh->GetContainingElementIndexWithInitialGuess(r_position_of_cell, mCellPdeElementMap[*cell_iter]); + mCellPdeElementMap[*cell_iter] = elem_index; + } +} + +template +void AbstractPottsPdeModifier::OutputSimulationModifierParameters(out_stream& rParamsFile) +{ + // No parameters to output, so just call method on direct parent class + AbstractPdeModifier::OutputSimulationModifierParameters(rParamsFile); +} + +// Explicit instantiation +template class AbstractPottsPdeModifier<1>; +template class AbstractPottsPdeModifier<2>; +template class AbstractPottsPdeModifier<3>; diff --git a/src/Test06Angiogenesis/AbstractPottsPdeModifier.hpp b/src/Test06Angiogenesis/AbstractPottsPdeModifier.hpp new file mode 100644 index 0000000..44e525e --- /dev/null +++ b/src/Test06Angiogenesis/AbstractPottsPdeModifier.hpp @@ -0,0 +1,299 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef ABSTRACTPottsPDEMODIFIER_HPP_ +#define ABSTRACTPottsPDEMODIFIER_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "AbstractPdeModifier.hpp" +#include "BoundaryConditionsContainer.hpp" + +/** + * An abstract modifier class containing functionality common to EllipticPottsPdeModifier + * and ParabolicPottsPdeModifier, which both solve a linear elliptic or parabolic PDE + * coupled to a cell-based simulation on a coarse domain. + */ +template +class AbstractPottsPdeModifier : public AbstractPdeModifier +{ + friend class TestEllipticPottsPdeModifier; + friend class TestParabolicPottsPdeModifier; + friend class TestOffLatticeSimulationWithPdes; + +private: + + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + archive & mpMeshCuboid; + archive & mStepSize; + archive & mSetBcsOnBoxBoundary; + archive & mSetBcsOnBoundingSphere; + archive & mSetBcsOnConvexHull; + archive & mUseVoronoiCellsForInterpolation; + archive & mTypicalCellRadius; + } + +protected: + + /** Map between cells and the elements of the FE mesh containing them. */ + std::map mCellPdeElementMap; + + /** + * Pointer to a ChasteCuboid storing the outer boundary for the FE mesh. + * Stored as a member to facilitate archiving. + */ + boost::shared_ptr > mpMeshCuboid; + + /** + * The step size to be used in the FE mesh. + * Stored as a member to facilitate archiving. + */ + double mStepSize; + + /** + * Whether to set the boundary condition on the edge of the box domain rather than the cell population. + * Default to true. + */ + bool mSetBcsOnBoxBoundary; + + /** + * Whether to set the boundary condition on a sphere which bounds the cell centres of the tissue. + * Only used if mSetBcsOnBoxBoundary is false. + * Default to false. + */ + bool mSetBcsOnBoundingSphere; + + /** + * Whether to set the boundary condition on a the convex hull which bounds the cell centres of the tissue. + * Only used if mSetBcsOnBoxBoundary is false. + * Default to false. + */ + bool mSetBcsOnConvexHull; + + /** + * Whether to use a cell centres voroni region to interpolate the pde solution onto cells. + */ + bool mUseVoronoiCellsForInterpolation; + + /** + * Used to define if a FE node is within a certain radius of a cell centre to help define + * boundary conditions when mSetBcsOnBoxBoundary and mSetBcsOnBoundingSphere are false + * + * defaults to 0.5 CD + */ + double mTypicalCellRadius; + +public: + + /** + * Constructor. + * + * @param pPde A shared pointer to a linear PDE object (defaults to NULL) + * @param pBoundaryCondition A shared pointer to an abstract boundary condition + * (defaults to NULL, corresponding to a constant boundary condition with value zero) + * @param isNeumannBoundaryCondition Whether the boundary condition is Neumann (defaults to true) + * @param pMeshCuboid A shared pointer to a ChasteCuboid specifying the outer boundary for the FE mesh (defaults to NULL) + * @param stepSize step size to be used in the FE mesh (defaults to 1.0, i.e. the default cell size) + * @param solution solution vector (defaults to NULL) + */ + AbstractPottsPdeModifier(boost::shared_ptr > pPde=boost::shared_ptr >(), + boost::shared_ptr > pBoundaryCondition=boost::shared_ptr >(), + bool isNeumannBoundaryCondition=true, + boost::shared_ptr > pMeshCuboid=boost::shared_ptr >(), + double stepSize=1.0, + Vec solution=nullptr); + + /** + * Destructor. + */ + virtual ~AbstractPottsPdeModifier(); + + /** + * @return mStepSize. + */ + double GetStepSize(); + + /** + * Set mSetBcsOnCoarseBoundary. + * + * @param setBcsOnBoxBoundary whether to set the boundary condition on the edge of the box domain rather than the cell population + */ + void SetBcsOnBoxBoundary(bool setBcsOnBoxBoundary); + + /** + * @return mSetBcsOnCoarseBoundary. + */ + bool AreBcsSetOnBoxBoundary(); + + /** + * Set mSetBcsOnBoundingSphere. + * + * @param setBcsOnBoundingSphere whether to set the boundary condition on the bounding spher of the cell population + */ + void SetBcsOnBoundingSphere(bool setBcsOnBoundingSphere); + + /** + * @return mSetBcsOnBoundingSphere. + */ + bool AreBcsSetOnBoundingSphere(); + + /** + * Set mSetBcsOnConvexHull. + * + * @param setBcsOnConvexHull whether to set the boundary condition on the convex hull of than the cell population + */ + void SetBcsOnConvexHull(bool setBcsOnConvexHull); + + /** + * @return mSetBcsOnConvexHull. + */ + bool AreBcsSetOnConvexHull(); + + /** + * Set mUseVoronoiCellsForInterpolation. + * + * @param useVoronoiCellsForInterpolation whether to use the voroni region of cells for interpolation + * of the solution from the FE mesh to the cells. + */ + void SetUseVoronoiCellsForInterpolation(bool useVoronoiCellsForInterpolation); + + /** + * @return mUseVoronoiCellsForInterpolation. + */ + bool GetUseVoronoiCellsForInterpolation(); + + /** + * Set mTypicalCellRadius. + * + * @param typicalCellRadius the radius to use for deining if FE nodes are near cells or not. + */ + void SetTypicalCellRadius(double typicalCellRadius); + + /** + * @return mTypicalCellRadius. + */ + double GetTypicalCellRadius(); + + + /** + * Helper method to construct the boundary conditions container for the PDE. + * + * @param rCellPopulation reference to the cell population + * @param pBcc the boundary conditions container to fill + */ + void ConstructBoundaryConditionsContainerHelper(AbstractCellPopulation& rCellPopulation, + std::shared_ptr > pBcc); + + /** + * Overridden SetupSolve() method. + * + * Specifies what to do in the simulation before the start of the time loop. + * + * Here we just initialize the Cell PDE element map + * + * @param rCellPopulation reference to the cell population + * @param outputDirectory the output directory, relative to where Chaste output is stored + */ + virtual void SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory); + + /** + * Helper method to generate the pde mesh for the first time. + * + * @param pMeshCuboid the outer boundary for the FE mesh. + * @param stepSize the step size to be used in the FE mesh. + */ + void GenerateFeMesh(boost::shared_ptr > pMeshCuboid, double stepSize); + + /** + * Helper method to generate a pde mesh. + * + * @param pMeshCuboid the outer boundary for the FE mesh. + * @param stepSize the step size to be used in the FE mesh. + * @param pMesh a pointer to the mesh to be created + */ + void GenerateAndReturnFeMesh(boost::shared_ptr > pMeshCuboid, double stepSize, TetrahedralMesh* pMesh); + + /** + * Helper method to copy the PDE solution to CellData + * + * Here we need to interpolate from the FE mesh onto the cells. + * + * @param rCellPopulation reference to the cell population + */ + void UpdateCellData(AbstractCellPopulation& rCellPopulation); + + /** + * Initialise mCellPdeElementMap. + * + * @param rCellPopulation reference to the cell population + */ + void InitialiseCellPdeElementMap(AbstractCellPopulation& rCellPopulation); + + /** + * Update the mCellPdeElementMap + * + * This method should be called before sending the element map to a PDE class + * to ensure map is up to date. + * + * @param rCellPopulation reference to the cell population + */ + void UpdateCellPdeElementMap(AbstractCellPopulation& rCellPopulation); + + /** + * Overridden OutputSimulationModifierParameters() method. + * Output any simulation modifier parameters to file. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputSimulationModifierParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +TEMPLATED_CLASS_IS_ABSTRACT_1_UNSIGNED(AbstractPottsPdeModifier) + +#endif /*ABSTRACTPottsPDEMODIFIER_HPP_*/ diff --git a/src/Test06Angiogenesis/ParabolicPottsPdeModifier.cpp b/src/Test06Angiogenesis/ParabolicPottsPdeModifier.cpp new file mode 100644 index 0000000..5c5fa2b --- /dev/null +++ b/src/Test06Angiogenesis/ParabolicPottsPdeModifier.cpp @@ -0,0 +1,413 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "ParabolicPottsPdeModifier.hpp" +#include "SimpleLinearParabolicSolver.hpp" +#include "VtkMeshWriter.hpp" +#include "MutableMesh.hpp" +#include "PottsBasedCellPopulation.hpp" + +template +ParabolicPottsPdeModifier::ParabolicPottsPdeModifier(boost::shared_ptr > pPde, + boost::shared_ptr > pBoundaryCondition, + bool isNeumannBoundaryCondition, + boost::shared_ptr > pMeshCuboid, + double stepSize, + Vec solution) + : AbstractPottsPdeModifier(pPde, + pBoundaryCondition, + isNeumannBoundaryCondition, + pMeshCuboid, + stepSize, + solution), + mMoveSolutionWithCells(false) +{ +} + +template +ParabolicPottsPdeModifier::~ParabolicPottsPdeModifier() +{ +} + +template +void ParabolicPottsPdeModifier::UpdateAtEndOfTimeStep(AbstractCellPopulation& rCellPopulation) +{ + // Set up boundary conditions + std::shared_ptr > p_bcc = ConstructBoundaryConditionsContainer(rCellPopulation); + + this->UpdateCellPdeElementMap(rCellPopulation); + + // When using a PDE mesh which doesn't coincide with the cells, we must set up the source terms before solving the PDE. + // Pass in already updated CellPdeElementMap to speed up finding cells. + this->SetUpSourceTermsForAveragedSourcePde(this->mpFeMesh, &this->mCellPdeElementMap); + + // Use SimpleLinearParabolicSolver as averaged Source PDE + SimpleLinearParabolicSolver solver(this->mpFeMesh, + boost::static_pointer_cast >(this->GetPde()).get(), + p_bcc.get()); + + ///\todo Investigate more than one PDE time step per spatial step + SimulationTime* p_simulation_time = SimulationTime::Instance(); + double current_time = p_simulation_time->GetTime(); + double dt = p_simulation_time->GetTimeStep(); + solver.SetTimes(current_time,current_time + dt); + solver.SetTimeStep(dt); + + + // Use previous solution as the initial condition + Vec previous_solution = this->mSolution; + if (mMoveSolutionWithCells) + { + // interpolate solution from cells movement onto fe mesh + previous_solution = InterpolateSolutionFromCellMovement(rCellPopulation); + } + solver.SetInitialCondition(previous_solution); + + // Note that the linear solver creates a vector, so we have to keep a handle on the old one + // in order to destroy it + this->mSolution = solver.Solve(); + PetscTools::Destroy(previous_solution); + + // Now copy solution to cells + this->UpdateCellData(rCellPopulation); + + // If using a PottsBasedCellPopulation then copy the PDE solution to the lattice sites. + if (dynamic_cast*>(&rCellPopulation)) + { + // Store the PDE solution in an accessible form + ReplicatableVector solution_repl(this->mSolution); + + unsigned num_nodes = rCellPopulation.GetNumNodes(); + + std::vector solution_on_lattice_sites(num_nodes, -1); + + // Loop over the nodes of the PottsMesh and get solution values from the mpFeMesh + for (typename AbstractMesh::NodeIterator node_iter = rCellPopulation.rGetMesh().GetNodeIteratorBegin(); + node_iter != rCellPopulation.rGetMesh().GetNodeIteratorEnd(); + ++node_iter) + { + unsigned node_index = node_iter->GetIndex(); + + const ChastePoint& node_location = node_iter->rGetLocation(); + + // interpolate PDE solution from FEMesh to node_location + // Find the element in the FE mesh that contains this PottsMesh node. + try + { + unsigned fe_elem_index = this->mpFeMesh->GetContainingElementIndex(node_location); + + // Now do the interpolation + Element* p_fe_element = this->mpFeMesh->GetElement(fe_elem_index); + c_vector weights; + + weights = p_fe_element->CalculateInterpolationWeights(node_location); + + double solution_at_node = 0.0; + for (unsigned i=0; iGetNodeGlobalIndex(i)]; + solution_at_node += nodal_value * weights(i); + } + solution_on_lattice_sites[node_index] = solution_at_node; + } + catch (Exception &e) + { + // Potts node is not within the FEMesh + + assert(0); + } + } + // Store the solution in the PottsBasedCellPopulation + static_cast*>(&rCellPopulation)->SetPdeSolution(solution_on_lattice_sites); + } + + + // Finally, if needed store the locations of cells to be used as old loactions in the next timestep + if (mMoveSolutionWithCells) + { + /* + * If required, store the current locations of cell centres. Note that we need to + * use a std::map between cells and locations, rather than (say) a std::vector with + * location indices corresponding to cells, since once we call UpdateCellLocations() + * the location index of each cell may change. This is especially true in the case + * of a CaBasedCellPopulation. + */ + mOldCellLocations.clear(); + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + mOldCellLocations[*cell_iter] = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + } + } +} + +template +void ParabolicPottsPdeModifier::SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory) +{ + AbstractPottsPdeModifier::SetupSolve(rCellPopulation,outputDirectory); + + // Copy the cell data to mSolution (this is the initial condition) + SetupInitialSolutionVector(rCellPopulation); + + // Output the initial conditions on FeMesh + this->UpdateAtEndOfOutputTimeStep(rCellPopulation); + + // If needed store the locations of cells to be used as old locations at the next timestep. + if (mMoveSolutionWithCells) + { + /* + * If required, store the current locations of cell centres. Note that we need to + * use a std::map between cells and locations, rather than (say) a std::vector with + * location indices corresponding to cells, since once we call UpdateCellLocations() + * the location index of each cell may change. This is especially true in the case + * of a CaBasedCellPopulation. + */ + mOldCellLocations.clear(); + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + mOldCellLocations[*cell_iter] = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + } + } +} + +template +std::shared_ptr > ParabolicPottsPdeModifier::ConstructBoundaryConditionsContainer(AbstractCellPopulation& rCellPopulation) +{ + std::shared_ptr > p_bcc(new BoundaryConditionsContainer(false)); + + this->ConstructBoundaryConditionsContainerHelper(rCellPopulation,p_bcc); + + return p_bcc; +} + +template +void ParabolicPottsPdeModifier::SetupInitialSolutionVector(AbstractCellPopulation& rCellPopulation) +{ + // Specify homogeneous initial conditions based upon the values stored in CellData. + // Note need all the CellDataValues to be the same. + + double initial_condition = rCellPopulation.Begin()->GetCellData()->GetItem(this->mDependentVariableName); + + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + double initial_condition_at_cell = cell_iter->GetCellData()->GetItem(this->mDependentVariableName); + UNUSED_OPT(initial_condition_at_cell); + assert(fabs(initial_condition_at_cell - initial_condition)<1e-12); + } + + // Initialise mSolution + this->mSolution = PetscTools::CreateAndSetVec(this->mpFeMesh->GetNumNodes(), initial_condition); +} + +template +Vec ParabolicPottsPdeModifier::InterpolateSolutionFromCellMovement(AbstractCellPopulation& rCellPopulation) +{ + // Store the PDE solution in an accessible form + ReplicatableVector solution_repl(this->mSolution); + + Vec interpolated_solution = PetscTools::CreateAndSetVec(this->mpFeMesh->GetNumNodes(), 0.0); + + // Store max radius so can speed up interpolation. + // TODO replace this with more general exclusion based on dimensions of tissue. + double max_radius = 0.0; + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + const ChastePoint& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + + double radius = norm_2(r_position_of_cell.rGetLocation()); + + if (max_radius < radius) + { + max_radius = radius; + } + } + + // Create mesh from cell centres so can interpolate tissue velocity onto mpFeMesh + std::vector*> temp_nodes; + std::vector> cell_displacements; + unsigned cell_index = 0; + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + c_vector position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter); + cell_index++; + + // Only use cells that were in both this and the previous timestep. + if (mOldCellLocations.find(*cell_iter) != mOldCellLocations.end()) + { + //PRINT_VARIABLE(mOldCellLocations[*cell_iter]); + c_vector displacement = position_of_cell - mOldCellLocations[*cell_iter]; + + cell_displacements.push_back(displacement); + temp_nodes.push_back(new Node(cell_index, position_of_cell)); + } + } + MutableMesh cell_mesh(temp_nodes); + + // Make the deformed mesh. Based on the displacement of cells. + TetrahedralMesh* p_deformed_mesh = new TetrahedralMesh(); + this->GenerateAndReturnFeMesh(this->mpMeshCuboid,this->mStepSize,p_deformed_mesh); + + for (typename TetrahedralMesh::NodeIterator node_iter = p_deformed_mesh->GetNodeIteratorBegin(); + node_iter != p_deformed_mesh->GetNodeIteratorEnd(); + ++node_iter) + { + c_vector node_location = node_iter->rGetLocation(); + + + c_vector new_node_location = node_location; + + if (norm_2(node_location) <= max_radius) + { + // Find the element in the cell mesh that contains this node. + try + { + unsigned elem_index = cell_mesh.GetContainingElementIndex(node_location, false); + + // Now do the interpolation + Element* p_element = cell_mesh.GetElement(elem_index); + c_vector weights = p_element->CalculateInterpolationWeights(node_location); + + c_vector interpolated_cell_displacement = zero_vector(DIM); + for (unsigned i=0; i nodal_value = cell_displacements[p_element->GetNodeGlobalIndex(i)]; + interpolated_cell_displacement += nodal_value * weights(i); + } + new_node_location = node_location + interpolated_cell_displacement; + node_iter->rGetModifiableLocation() = new_node_location; + } + catch (Exception&) // not_in_mesh + { + //Don't do anything as these FE nodes are outside of the Cell mesh. + } + } + } + + // Loop over nodes of the mpFeMesh and get solution values from the deformed mesh + for (typename TetrahedralMesh::NodeIterator node_iter = this->mpFeMesh->GetNodeIteratorBegin(); + node_iter != this->mpFeMesh->GetNodeIteratorEnd(); + ++node_iter) + { + unsigned node_index = node_iter->GetIndex(); + + const ChastePoint& node_location = node_iter->rGetLocation(); + + // Find the element in the deformed mesh that contains this FE node. + std::set elements_to_check = node_iter->rGetContainingElementIndices(); + try + { + unsigned elem_index = p_deformed_mesh->GetContainingElementIndex(node_location, false, elements_to_check); + + // Now do the interpolation + Element* p_element = p_deformed_mesh->GetElement(elem_index); + c_vector weights; + + weights = p_element->CalculateInterpolationWeights(node_location); + + double solution_at_node = 0.0; + for (unsigned i=0; iGetNodeGlobalIndex(i)]; + solution_at_node += nodal_value * weights(i); + } + PetscVecTools::SetElement(interpolated_solution, node_index, solution_at_node); + } + catch (Exception &e) + { + // Handy debug code to work out why the node is not in any element + + // // output the cell mesh + // std::ostringstream time_string1; + // time_string1 << SimulationTime::Instance()->GetTimeStepsElapsed(); + // std::string results_file1 = "cell_mesh_" + this->mDependentVariableName + "_" + time_string1.str(); + // VtkMeshWriter* p_vtk_mesh_writer1 = new VtkMeshWriter(this->mOutputDirectory, results_file1, false); + // p_vtk_mesh_writer1->AddPointData(this->mDependentVariableName, cell_displacements); + // p_vtk_mesh_writer1->WriteFilesUsingMesh(cell_mesh); + // delete p_vtk_mesh_writer1; + + // // output the deformed mesh + // std::ostringstream time_string; + // time_string << SimulationTime::Instance()->GetTimeStepsElapsed(); + // std::string results_file = "deformed_mesh_" + this->mDependentVariableName + "_" + time_string.str(); + // VtkMeshWriter* p_vtk_mesh_writer = new VtkMeshWriter(this->mOutputDirectory, results_file, false); + // p_vtk_mesh_writer->WriteFilesUsingMesh(*p_deformed_mesh); + // delete p_vtk_mesh_writer; + + assert(0); + } + } + + // Tidy Up + delete p_deformed_mesh; + + return interpolated_solution; +} + +template +void ParabolicPottsPdeModifier::SetMoveSolutionWithCells(bool moveSolutionWithCells) +{ + mMoveSolutionWithCells = moveSolutionWithCells; +} + +template +bool ParabolicPottsPdeModifier::GetMoveSolutionWithCells() +{ + return mMoveSolutionWithCells; +} + +template +void ParabolicPottsPdeModifier::OutputSimulationModifierParameters(out_stream& rParamsFile) +{ + // No parameters to output, so just call method on direct parent class + AbstractPottsPdeModifier::OutputSimulationModifierParameters(rParamsFile); +} + +// Explicit instantiation +template class ParabolicPottsPdeModifier<1>; +template class ParabolicPottsPdeModifier<2>; +template class ParabolicPottsPdeModifier<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(ParabolicPottsPdeModifier) diff --git a/src/Test06Angiogenesis/ParabolicPottsPdeModifier.hpp b/src/Test06Angiogenesis/ParabolicPottsPdeModifier.hpp new file mode 100644 index 0000000..3104753 --- /dev/null +++ b/src/Test06Angiogenesis/ParabolicPottsPdeModifier.hpp @@ -0,0 +1,237 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef PARABOLICPottsPDEMODIFIER_HPP_ +#define PARABOLICPottsPDEMODIFIER_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "AbstractPottsPdeModifier.hpp" +#include "BoundaryConditionsContainer.hpp" + +/** + * A modifier class in which a linear parabolic PDE coupled to a cell-based simulation + * is solved on a coarse domain. + * + * The finite element mesh used to solve the PDE numerically is a fixed tessellation of + * a cuboid (box), which must be supplied to the constructor. The value of the dependent + * variable is interpolated between coarse mesh nodes to obtain a value at each cell, + * which is stored and updated in a CellData item. + * + * At each time step the boundary condition supplied to the constructor may be imposed + * either on the boundary of the box domain, or on the boundary of the cell population + * (which is assumed to lie within the box domain). This choice can be made using the + * AbstractPottsPdeModifier method SetBcsOnBoxBoundary(), which is inherited by this + * class. + * + * Examples of PDEs in the source folder that can be solved using this class are + * AveragedSourceParabolicPde and UniformSourceParabolicPde. + */ +template +class ParabolicPottsPdeModifier : public AbstractPottsPdeModifier +{ + friend class TestParabolicPottsPdeModifier; + +private: + + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + archive & mMoveSolutionWithCells; + archive & mOldCellLocations; + } + + /** + * Whether to move the solution along with the cells in the cell population. + * Default to true. + */ + bool mMoveSolutionWithCells; + + /** + * Map used to calculate displacement of cells if moving pde solution with cells + */ + std::map > mOldCellLocations; + + +public: + + /** + * Constructor. + * + * @param pPde A shared pointer to a linear PDE object (defaults to NULL) + * @param pBoundaryCondition A shared pointer to an abstract boundary condition + * (defaults to NULL, corresponding to a constant boundary condition with value zero) + * @param isNeumannBoundaryCondition Whether the boundary condition is Neumann (defaults to true) + * @param pMeshCuboid A shared pointer to a ChasteCuboid specifying the outer boundary for the FE mesh (defaults to NULL) + * @param stepSize step size to be used in the FE mesh (defaults to 1.0, i.e. the default cell size) + * @param solution solution vector (defaults to NULL) + */ + ParabolicPottsPdeModifier(boost::shared_ptr > pPde=boost::shared_ptr >(), + boost::shared_ptr > pBoundaryCondition=boost::shared_ptr >(), + bool isNeumannBoundaryCondition=true, + boost::shared_ptr > pMeshCuboid=boost::shared_ptr >(), + double stepSize=1.0, + Vec solution=nullptr); + + /** + * Destructor. + */ + virtual ~ParabolicPottsPdeModifier(); + + /** + * Overridden UpdateAtEndOfTimeStep() method. + * + * Specifies what to do in the simulation at the end of each time step. + * + * @param rCellPopulation reference to the cell population + */ + virtual void UpdateAtEndOfTimeStep(AbstractCellPopulation& rCellPopulation); + + /** + * Overridden SetupSolve() method. + * + * Specifies what to do in the simulation before the start of the time loop. + * + * @param rCellPopulation reference to the cell population + * @param outputDirectory the output directory, relative to where Chaste output is stored + */ + virtual void SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory); + + /** + * Helper method to construct the boundary conditions container for the PDE. + * + * @param rCellPopulation reference to the cell population + * + * @return the full boundary conditions container + */ + virtual std::shared_ptr > ConstructBoundaryConditionsContainer(AbstractCellPopulation& rCellPopulation); + + /** + * Helper method to initialise the PDE solution using the CellData. + * + * Here we assume a homogeneous initial consition. + * + * TODO use InterpolateSolutionFromCellMovement instead! + * + * @param rCellPopulation reference to the cell population + */ + void SetupInitialSolutionVector(AbstractCellPopulation& rCellPopulation); + + /** + * Helper method to interpolate the PDE solution from cells using the CellData. + * Use the cells voronoi region as all fe nodes in a cell centre voronoi region are given the + * value from cell data. + * + * @param rCellPopulation reference to the cell population + * + * @return the solution interpolated onto the FE Mesh + */ + Vec InterpolateSolutionFromCellMovement(AbstractCellPopulation& rCellPopulation); + + /** + * Set mMoveSolutionWithCells. + * + * @param moveSolutionWithCells whether to move the solution with cells. + */ + void SetMoveSolutionWithCells(bool moveSolutionWithCells); + + /** + * @return mMoveSolutionWithCells. + */ + bool GetMoveSolutionWithCells(); + + /** + * Overridden OutputSimulationModifierParameters() method. + * Output any simulation modifier parameters to file. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputSimulationModifierParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(ParabolicPottsPdeModifier) + +namespace boost +{ +namespace serialization +{ +template +inline void save_construct_data( + Archive & ar, const ParabolicPottsPdeModifier * t, const unsigned int file_version) +{ + if (t->GetSolution()) + { + std::string archive_filename = ArchiveLocationInfo::GetArchiveDirectory() + "solution.vec"; + PetscTools::DumpPetscObject(t->GetSolution(), archive_filename); + } +} + +template +inline void load_construct_data( + Archive & ar, ParabolicPottsPdeModifier * t, const unsigned int file_version) +{ + Vec solution = nullptr; + + std::string archive_filename = ArchiveLocationInfo::GetArchiveDirectory() + "solution.vec"; + FileFinder file_finder(archive_filename, RelativeTo::Absolute); + + if (file_finder.Exists()) + { + PetscTools::ReadPetscObject(solution, archive_filename); + } + + ::new(t)ParabolicPottsPdeModifier(boost::shared_ptr >(), + boost::shared_ptr >(), + true, + boost::shared_ptr >(), + 1.0, + solution); +} +} +} // namespace ... + +#endif /*PARABOLICPottsPDEMODIFIER_HPP_*/ diff --git a/src/Test06Angiogenesis/PottsParabolicPde.cpp b/src/Test06Angiogenesis/PottsParabolicPde.cpp new file mode 100644 index 0000000..1838c38 --- /dev/null +++ b/src/Test06Angiogenesis/PottsParabolicPde.cpp @@ -0,0 +1,223 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "PottsParabolicPde.hpp" +#include "ApoptoticCellProperty.hpp" + +template +PottsParabolicPde::PottsParabolicPde(PottsBasedCellPopulation& rCellPopulation, + double constantCellSourceCoefficient, + double linearCellSourceCoefficient, + double constantSourceCoefficient, + double linearSourceCoefficient, + double diffusionCoefficient, + double duDtCoefficient, + bool scaleByCellVolume) + : mrCellPopulation(rCellPopulation), + mConstantCellSourceCoefficient(constantCellSourceCoefficient), + mLinearCellSourceCoefficient(linearCellSourceCoefficient), + mConstantSourceCoefficient(constantSourceCoefficient), + mLinearSourceCoefficient(linearSourceCoefficient), + mDiffusionCoefficient(diffusionCoefficient), + mDuDtCoefficient(duDtCoefficient), + mScaleByCellVolume(scaleByCellVolume) +{ +} + +template +const PottsBasedCellPopulation& PottsParabolicPde::rGetCellPopulation() const +{ + return mrCellPopulation; +} + +template +double PottsParabolicPde::GetConstantCellCoefficient() const +{ + return mConstantCellSourceCoefficient; +} + +template +double PottsParabolicPde::GetLinearCellCoefficient() const +{ + return mLinearCellSourceCoefficient; +} + +template +double PottsParabolicPde::GetConstantCoefficient() const +{ + return mConstantSourceCoefficient; +} + +template +double PottsParabolicPde::GetLinearCoefficient() const +{ + return mLinearSourceCoefficient; +} + +template +double PottsParabolicPde::GetDiffusionCoefficient() const +{ + return mDiffusionCoefficient; +} + +template +double PottsParabolicPde::GetDuDtCoefficient() const +{ + return mDuDtCoefficient; +} + +template +bool PottsParabolicPde::GetScaleByCellVolume() const +{ + return mScaleByCellVolume; +} + +template +void PottsParabolicPde::SetupSourceTerms(TetrahedralMesh& rCoarseMesh, std::map* pCellPdeElementMap) // must be called before solve +{ + NEVER_REACHED + // Allocate memory + mCellDensityOnCoarseElements.resize(rCoarseMesh.GetNumElements()); + for (unsigned elem_index=0; elem_index::Iterator cell_iter = mrCellPopulation.Begin(); + cell_iter != mrCellPopulation.End(); + ++cell_iter) + { + unsigned elem_index = UNSIGNED_UNSET; + const ChastePoint& r_position_of_cell = mrCellPopulation.GetLocationOfCellCentre(*cell_iter); + + if (pCellPdeElementMap != nullptr) + { + elem_index = (*pCellPdeElementMap)[*cell_iter]; + } + else + { + elem_index = rCoarseMesh.GetContainingElementIndex(r_position_of_cell); + } + + bool cell_is_apoptotic = cell_iter->template HasCellProperty(); + + if (!cell_is_apoptotic) + { + double cell_weight = 1.0; + + if (mScaleByCellVolume) + { + // If scaling by cell volume then use volume here instead of cell count + cell_weight = mrCellPopulation.GetVolumeOfCell(*cell_iter); + + if (cell_weight <1e-6) + { + EXCEPTION("The volume of one of the cells is " << cell_weight << + " and you are scaling by cell volume. Either turn scaling off or use" + " a cell model with non zero areas (i.e. a Bounded Voronoi Tesselation model)."); + } + } + + mCellDensityOnCoarseElements[elem_index] += cell_weight; + } + } + } + + // Then divide each entry of mSourceTermOnCoarseElements by the element's area + c_matrix jacobian; + double det; + for (unsigned elem_index=0; elem_indexCalculateJacobian(jacobian, det); + mCellDensityOnCoarseElements[elem_index] /= rCoarseMesh.GetElement(elem_index)->GetVolume(det); + } +} + +template +double PottsParabolicPde::ComputeDuDtCoefficientFunction(const ChastePoint& ) +{ + return mDuDtCoefficient; +} + +template +double PottsParabolicPde::ComputeSourceTerm(const ChastePoint& rX, double u, Element* pElement) +{ + + assert(!mCellDensityOnCoarseElements.empty()); + + // The source term is (a*density + c)*u + b*density + d + bool is_in_cell = 1; + + double constant_source_term = mConstantCellSourceCoefficient * is_in_cell + mConstantSourceCoefficient; + double linear_source_term_coeficient = mLinearCellSourceCoefficient * is_in_cell + mLinearSourceCoefficient; + + double source_term = constant_source_term + linear_source_term_coeficient * u; + + return source_term; + //return 25*mCellDensityOnCoarseElements[pElement->GetIndex()] - u; +} + +// LCOV_EXCL_START +template +double PottsParabolicPde::ComputeSourceTermAtNode(const Node& rNode, double u) +{ + NEVER_REACHED; + return 0.0; +} +// LCOV_EXCL_STOP + +template +c_matrix PottsParabolicPde::ComputeDiffusionTerm(const ChastePoint& rX, Element* pElement) +{ + return mDiffusionCoefficient*identity_matrix(DIM); +} + +template +double PottsParabolicPde::GetUptakeRateForElement(unsigned elementIndex) +{ + NEVER_REACHED; + return this->mCellDensityOnCoarseElements[elementIndex]; +} + +// Explicit instantiation +template class PottsParabolicPde<1>; +template class PottsParabolicPde<2>; +template class PottsParabolicPde<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(PottsParabolicPde) diff --git a/src/Test06Angiogenesis/PottsParabolicPde.hpp b/src/Test06Angiogenesis/PottsParabolicPde.hpp new file mode 100644 index 0000000..642b126 --- /dev/null +++ b/src/Test06Angiogenesis/PottsParabolicPde.hpp @@ -0,0 +1,282 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef POTTSPARABOLICPDE_HPP_ +#define POTTSPARABOLICPDE_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "PottsBasedCellPopulation.hpp" +#include "TetrahedralMesh.hpp" +#include "AbstractLinearParabolicPde.hpp" + +/** + * A parabolic PDE to be solved numerically using the finite element method, for + * coupling to a cell-based simulation. + * + * The PDE takes the form + * + * c*du/dt = Grad.(D*Grad(u)) + (a*u + b) * rho(x), + * + * where the scalars c, D, a and b are specified by the members mDuDtCoefficient, + * mDiffusionCoefficient, mLinearSourceCoefficient and mConstantSourceCoefficient, + * respectively. Their values must be set in the constructor. + * + * The function rho(x) denotes the local density of non-apoptotic cells. This + * quantity is computed for each element of a 'coarse' finite element mesh that is + * passed to the method SetupSourceTerms() and stored in the member mCellDensityOnCoarseElements. + * For a point x, rho(x) is defined to be the number of non-apoptotic cells whose + * centres lie in each finite element containing that point, scaled by the area of + * that element. + */ +template +class PottsParabolicPde : public AbstractLinearParabolicPde +{ + friend class TestCellBasedParabolicPdes; + +private: + + /** Needed for serialization.*/ + friend class boost::serialization::access; + /** + * Serialize the PDE and its member variables. + * + * @param archive the archive + * @param version the current version of this class + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + archive & mConstantCellSourceCoefficient; + archive & mLinearCellSourceCoefficient; + archive & mConstantSourceCoefficient; + archive & mLinearSourceCoefficient; + archive & mDiffusionCoefficient; + archive & mDuDtCoefficient; + archive & mScaleByCellVolume; + archive & mCellDensityOnCoarseElements; + } + +protected: + + /** The cell population member. */ + PottsBasedCellPopulation& mrCellPopulation; + + /** Coefficient of constant density dependent source term. */ + double mConstantCellSourceCoefficient; + + /** Coefficient of linear density dependent source term. */ + double mLinearCellSourceCoefficient; + + /** Coefficient of constant source term. */ + double mConstantSourceCoefficient; + + /** Coefficient of linear source term. */ + double mLinearSourceCoefficient; + + /** Diffusion coefficient. */ + double mDiffusionCoefficient; + + /** Coefficient of rate of change term. */ + double mDuDtCoefficient; + + /** Whether to scale tems by cell volume*/ + bool mScaleByCellVolume; + + /** Vector of averaged cell densities on elements of the coarse mesh. */ + std::vector mCellDensityOnCoarseElements; + +public: + + /** + * Constructor. + * + * @param rCellPopulation reference to the potts based cell population + * @param constantCellSourceCoefficient the constant density dependent source term coefficient (defaults to 0.0) + * @param linearCellSourceCoefficient the linear density dependent source term coefficient (defaults to 0.0) + * @param constantSourceCoefficient the constant source term coefficient (defaults to 0.0) + * @param linearSourceCoefficient the linear source term coefficient (defaults to 0.0) + * @param diffusionCoefficient the rate of diffusion (defaults to 1.0) + * @param duDtCoefficient rate of reaction (defaults to 1.0) + * @param scaleByCellVolume whether to scale by cell volume (defaults to false) + */ + PottsParabolicPde(PottsBasedCellPopulation& rCellPopulation, + double constantCellSourceCoefficient=0.0, + double linearCellSourceCoefficient=0.0, + double constantSourceCoefficient=0.0, + double linearSourceCoefficient=0.0, + double diffusionCoefficient=1.0, + double duDtCoefficient=1.0, + bool scaleByCellVolume=false); + + /** + * @return const reference to the cell population (used in archiving). + */ + const PottsBasedCellPopulation& rGetCellPopulation() const; + + /** + * @return mConstantCellSourceCoefficient + */ + double GetConstantCellCoefficient() const; + + /** + * @return mLinearCellSourceCoefficient + */ + double GetLinearCellCoefficient() const; + + /** + * @return mConstantSourceCoefficient + */ + double GetConstantCoefficient() const; + + /** + * @return mLinearSourceCoefficient + */ + double GetLinearCoefficient() const; + + /** + * @return mDiffusionCoefficient + */ + double GetDiffusionCoefficient() const; + + /** + * @return mDuDtCoefficient + */ + double GetDuDtCoefficient() const; + + /** + * @return mScaleByCellVolume + */ + bool GetScaleByCellVolume() const; + + /** + * Set up the source terms. + * + * \todo this is identical to the one in PottsEllipticPde so refactor. + * + * @param rCoarseMesh reference to the coarse mesh + * @param pCellPdeElementMap optional pointer to the map from cells to coarse elements + */ + void virtual SetupSourceTerms(TetrahedralMesh& rCoarseMesh, std::map* pCellPdeElementMap=nullptr); + + /** + * Overridden ComputeDuDtCoefficientFunction() method. + * + * @return the function c(x) in "c(x) du/dt = Grad.(DiffusionTerm(x)*Grad(u))+LinearSourceTerm(x)+NonlinearSourceTerm(x, u)" + * + * @param rX the point in space at which the function c is computed + */ + virtual double ComputeDuDtCoefficientFunction(const ChastePoint& rX); + + /** + * Overridden ComputeSourceTerm() method. + * + * @return computed source term. + * + * @param rX the point in space at which the nonlinear source term is computed + * @param u the value of the dependent variable at the point + * @param pElement the mesh element that x is contained in (optional; defaults to NULL). + */ + virtual double ComputeSourceTerm(const ChastePoint& rX, + double u, + Element* pElement=NULL); + + /** + * Overridden ComputeSourceTermAtNode() method. That is never called. + * + * @return computed source term at a node. + * + * @param rNode the node at which the nonlinear source term is computed + * @param u the value of the dependent variable at the node + */ + virtual double ComputeSourceTermAtNode(const Node& rNode, double u); + + /** + * Overridden ComputeDiffusionTerm() method. + * + * @param rX the point in space at which the diffusion term is computed + * @param pElement the mesh element that x is contained in (optional; defaults to NULL). + * + * @return a matrix. + */ + virtual c_matrix ComputeDiffusionTerm(const ChastePoint& rX, Element* pElement=NULL); + + /** + * @return the uptake rate. + * + * @param elementIndex the element we wish to return the uptake rate for + */ + double GetUptakeRateForElement(unsigned elementIndex); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(PottsParabolicPde) + +namespace boost +{ +namespace serialization +{ +/** + * Serialize information required to construct a PottsParabolicPde. + */ +template +inline void save_construct_data( + Archive & ar, const PottsParabolicPde* t, const unsigned int file_version) +{ + // Save data required to construct instance + const PottsBasedCellPopulation* p_cell_population = &(t->rGetCellPopulation()); + ar & p_cell_population; +} + +/** + * De-serialize constructor parameters and initialise a PottsParabolicPde. + */ +template +inline void load_construct_data( + Archive & ar, PottsParabolicPde* t, const unsigned int file_version) +{ + // Retrieve data from archive required to construct new instance + PottsBasedCellPopulation* p_cell_population; + ar >> p_cell_population; + + // Invoke inplace constructor to initialise instance + ::new(t)PottsParabolicPde(*p_cell_population); +} +} +} // namespace ... + +#endif /*AVERAGESOURCEPARABOLICPDE_HPP_*/ diff --git a/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.cpp b/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.cpp new file mode 100644 index 0000000..6ec4d3f --- /dev/null +++ b/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.cpp @@ -0,0 +1,110 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "VegfChemotaxisPottsUpdateRule.hpp" +#include "Debug.hpp" + +template +VegfChemotaxisPottsUpdateRule::VegfChemotaxisPottsUpdateRule() + : AbstractPottsUpdateRule() +{ +} + +template +VegfChemotaxisPottsUpdateRule::~VegfChemotaxisPottsUpdateRule() +{ +} + +template +double VegfChemotaxisPottsUpdateRule::EvaluateHamiltonianContribution(unsigned currentNodeIndex, + unsigned targetNodeIndex, + PottsBasedCellPopulation& rCellPopulation) +{ + double CellBoundaryChemotaxisParameter = 500.0; + //double CellCellChemotaxisParameter = 0.0; + + + std::set containing_elements = rCellPopulation.GetNode(currentNodeIndex)->rGetContainingElementIndices(); + std::set new_location_containing_elements = rCellPopulation.GetNode(targetNodeIndex)->rGetContainingElementIndices(); + + // Every node must each be in at most one element + assert(containing_elements.size() < 2); + assert(new_location_containing_elements.size() < 2); + + + bool current_node_contained = !containing_elements.empty(); + bool target_node_contained = !new_location_containing_elements.empty(); + + if (!current_node_contained && !target_node_contained) + { + EXCEPTION("At least one of the current node or target node must be in an element."); + } + + double delta_H = 0.0; + + if (current_node_contained && target_node_contained) + { + // both cells so return 0 + } + else + { + assert((!current_node_contained && target_node_contained) || (current_node_contained &&!target_node_contained)); + + double current_vegf = rCellPopulation.GetPdeSolutionAtNode(currentNodeIndex); + + double target_vegf = rCellPopulation.GetPdeSolutionAtNode(targetNodeIndex); +//PRINT_2_VARIABLES(current_vegf,target_vegf); + + delta_H = 1000*CellBoundaryChemotaxisParameter * (current_vegf - target_vegf); + + } + return delta_H; +} + +template +void VegfChemotaxisPottsUpdateRule::OutputUpdateRuleParameters(out_stream& rParamsFile) +{ + // Call method on direct parent class + AbstractPottsUpdateRule::OutputUpdateRuleParameters(rParamsFile); +} + +// Explicit instantiation +template class VegfChemotaxisPottsUpdateRule<1>; +template class VegfChemotaxisPottsUpdateRule<2>; +template class VegfChemotaxisPottsUpdateRule<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(VegfChemotaxisPottsUpdateRule) diff --git a/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.hpp b/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.hpp new file mode 100644 index 0000000..e4dd217 --- /dev/null +++ b/src/Test06Angiogenesis/VegfChemotaxisPottsUpdateRule.hpp @@ -0,0 +1,111 @@ +/* + +Copyright (c) 2005-2024, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef VEGFCHEMOTAXISPOTTSUPDATERULE_HPP_ +#define VEGFCHEMOTAXISPOTTSUPDATERULE_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "AbstractPottsUpdateRule.hpp" +#include "PottsBasedCellPopulation.hpp" + +/** + * A simple update rule class to represent simple chemotaxis in the positive x y z direction for use in Potts based simulations. + * + * Note this currently assumes cells don't grow, i.e the target volume is constant + * for each cell over time. + */ +template +class VegfChemotaxisPottsUpdateRule : public AbstractPottsUpdateRule +{ +friend class TestPottsUpdateRules; + +private: + + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + } + +public: + + /** + * Constructor. + */ + VegfChemotaxisPottsUpdateRule(); + + /** + * Destructor. + */ + ~VegfChemotaxisPottsUpdateRule(); + + /** + * Overridden EvaluateHamiltonianContribution() method + * + * Assigns a greater propensity for moving in an increasing x, y (and z) direction. + * + * @param currentNodeIndex The index of the current node/lattice site + * @param targetNodeIndex The index of the target node/lattice site + * @param rCellPopulation The cell population + * + * @return The difference in the Hamiltonian with the configuration of the target node + * having the same spin as the current node with the current configuration. i.e H_1-H_0 + */ + double EvaluateHamiltonianContribution(unsigned currentNodeIndex, + unsigned targetNodeIndex, + PottsBasedCellPopulation& rCellPopulation); + + /** + * Overridden OutputUpdateRuleParameters() method. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputUpdateRuleParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(VegfChemotaxisPottsUpdateRule) + +#endif /*VEGFCHEMOTAXISPOTTSUPDATERULE_HPP_*/ diff --git a/test/Test06Angiogenesis.hpp b/test/Test06Angiogenesis.hpp index fb9b33a..c743b73 100644 --- a/test/Test06Angiogenesis.hpp +++ b/test/Test06Angiogenesis.hpp @@ -49,14 +49,18 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "OnLatticeSimulation.hpp" #include "PottsBasedCellPopulation.hpp" #include "VolumeConstraintPottsUpdateRule.hpp" +#include "SurfaceAreaConstraintPottsUpdateRule.hpp" #include "AdhesionPottsUpdateRule.hpp" +#include "VegfChemotaxisPottsUpdateRule.hpp" #include "DifferentialAdhesionPottsUpdateRule.hpp" #include "TransitCellProliferativeType.hpp" #include "CellLabel.hpp" #include "CellLabelWriter.hpp" -#include "AveragedSourceParabolicPde.hpp" -#include "ParabolicBoxDomainPdeModifier.hpp" +// #include "AveragedSourceEllipticPde.hpp" +// #include "EllipticPottsPdeModifier.hpp" +#include "PottsParabolicPde.hpp" +#include "ParabolicPottsPdeModifier.hpp" class Test06Angiogenesis: public AbstractCellBasedWithTimingsTestSuite { @@ -69,76 +73,95 @@ class Test06Angiogenesis: public AbstractCellBasedWithTimingsTestSuite */ void TestMonolayer() { + double M_DOMAIN_WIDTH = 50.0; - double M_DIFFUSION_COEFICIENT = 1.0; - double M_CONSTANT_SECRETION = 1.0; - double M_LINEAR_CELL_SECRETION = 1.0; + double M_DUDT_COEFICIENT = 1.0; + double M_DIFFUSION_COEFICIENT = 1e2; //1e-3 + double M_CONSTANT_CELL_SECRETION = 1e-2; //1e-3 + double M_LINEAR_SECRETION = 1e-2; //1e-3 double cell_area = 25.0; double box_h = 1.0; + double dt = 0.01; + double end_time = 1.0; + double temperature = 1; EXIT_IF_PARALLEL; - PottsMeshGenerator<2> generator((int) M_DOMAIN_WIDTH, 2, 4, (int) M_DOMAIN_WIDTH, 2, 4); // Parameters are: lattice sites across; num elements across; element width; lattice sites up; num elements up; and element height + PottsMeshGenerator<2> generator((int) M_DOMAIN_WIDTH, 4, 4, (int) M_DOMAIN_WIDTH, 4, 4, 1, 1, 1, 4u); // Parameters are: lattice sites across; num elements across; element width; lattice sites up; num elements up; element height; and element spacing. boost::shared_ptr > p_mesh = generator.GetMesh(); std::vector cells; - MAKE_PTR(TransitCellProliferativeType, p_transit_type); + MAKE_PTR(DifferentiatedCellProliferativeType, p_diff_type); CellsGenerator cells_generator; - cells_generator.GenerateBasicRandom(cells, p_mesh->GetNumElements(), p_transit_type); + cells_generator.GenerateBasicRandom(cells, p_mesh->GetNumElements(), p_diff_type); PottsBasedCellPopulation<2> cell_population(*p_mesh, cells); - cell_population.SetTemperature(0.1); + cell_population.SetTemperature(temperature); cell_population.SetNumSweepsPerTimestep(1); OnLatticeSimulation<2> simulator(cell_population); simulator.SetOutputDirectory("PottsBasedMonolayer"); - simulator.SetEndTime(50.0); - simulator.SetDt(0.1); + simulator.SetEndTime(end_time); + simulator.SetDt(dt); simulator.SetSamplingTimestepMultiple(10); MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule); p_volume_constraint_update_rule->SetMatureCellTargetVolume(16); - p_volume_constraint_update_rule->SetDeformationEnergyParameter(0.2); + p_volume_constraint_update_rule->SetDeformationEnergyParameter(25); simulator.AddUpdateRule(p_volume_constraint_update_rule); + // MAKE_PTR(SurfaceAreaConstraintPottsUpdateRule<2>, p_surface_area_update_rule); + // p_surface_area_update_rule->SetMatureCellTargetSurfaceArea(16); + // p_surface_area_update_rule->SetDeformationEnergyParameter(0.5); + // simulator.AddUpdateRule(p_surface_area_update_rule); + MAKE_PTR(AdhesionPottsUpdateRule<2>, p_adhesion_update_rule); + p_adhesion_update_rule->SetCellCellAdhesionEnergyParameter(40); + p_adhesion_update_rule->SetCellBoundaryAdhesionEnergyParameter(20); simulator.AddUpdateRule(p_adhesion_update_rule); - // Set initial conditions - cell_population.SetDataOnAllCells("VEGF", 0.0); - - // Create a ChasteCuboid on which to base the finite element mesh used to solve the PDE - ChastePoint<2> lower(-0.5,-0.5); - ChastePoint<2> upper(M_DOMAIN_WIDTH-0.5, M_DOMAIN_WIDTH-0.5); - MAKE_PTR_ARGS(ChasteCuboid<2>, p_cuboid, (lower, upper)); - - // Create PDE and BCS - double dudt_coeficient = 0.01; - double diffusion_coeficient = M_DIFFUSION_COEFICIENT; - // double constant_cell_uptake = M_CONSTANT_SECRETION*cell_area; - // double linear_cell_uptake = 0.0; - // double constant_uptake = 0.0; - // double linear_uptake = - M_LINEAR_CELL_SECRETION; - double constant_uptake = M_CONSTANT_SECRETION*cell_area; - double linear_uptake = - M_LINEAR_CELL_SECRETION; - bool is_volume_scaled = false; - MAKE_PTR_ARGS(AveragedSourceParabolicPde<2>, p_pde, (cell_population, constant_uptake, linear_uptake, diffusion_coeficient, dudt_coeficient, is_volume_scaled)); - MAKE_PTR_ARGS(ConstBoundaryCondition<2>, p_bc, (0.0)); + MAKE_PTR(VegfChemotaxisPottsUpdateRule<2>, p_chemotaxis_update_rule); + //p_chemotaxis_update_rule->SetChemotaxisParameter(40); + simulator.AddUpdateRule(p_chemotaxis_update_rule); + + // Set initial conditions + cell_population.SetDataOnAllCells("VEGF", 0.0); + + // Create a ChasteCuboid on which to base the finite element mesh used to solve the PDE + ChastePoint<2> lower(-0.5,-0.5); + ChastePoint<2> upper(M_DOMAIN_WIDTH-0.5, M_DOMAIN_WIDTH-0.5); + MAKE_PTR_ARGS(ChasteCuboid<2>, p_cuboid, (lower, upper)); + + // Create PDE and BCS + double dudt_coeficient = M_DUDT_COEFICIENT; + double diffusion_coeficient = M_DIFFUSION_COEFICIENT; + double constant_cell_uptake = M_CONSTANT_CELL_SECRETION*cell_area; + double linear_cell_uptake = M_LINEAR_SECRETION*cell_area; + double constant_uptake = 0.0; + double linear_uptake = - M_LINEAR_SECRETION; + // double constant_uptake = M_CONSTANT_SECRETION*cell_area; + // double linear_uptake = - M_LINEAR_CELL_SECRETION; + bool is_volume_scaled = false; + MAKE_PTR_ARGS(PottsParabolicPde<2>, p_pde, (cell_population, constant_cell_uptake, linear_cell_uptake, constant_uptake, linear_uptake, diffusion_coeficient, dudt_coeficient, is_volume_scaled)); + //MAKE_PTR_ARGS(AveragedSourceParabolicPde<2>, p_pde, (cell_population, constant_cell_uptake, linear_cell_uptake, constant_uptake, linear_uptake, diffusion_coeficient, dudt_coeficient, is_volume_scaled)); + //MAKE_PTR_ARGS(AveragedSourceEllipticPde<2>, p_pde, (cell_population, constant_cell_uptake, linear_cell_uptake, constant_uptake, linear_uptake, diffusion_coeficient, is_volume_scaled)); + MAKE_PTR_ARGS(ConstBoundaryCondition<2>, p_bc, (0.0)); + + + // Create a PDE modifier and set the name of the dependent variable in the PDE + bool is_neuman_bcs = false; - - // Create a PDE modifier and set the name of the dependent variable in the PDE - bool is_neuman_bcs = false; - - MAKE_PTR_ARGS(ParabolicBoxDomainPdeModifier<2>, p_pde_modifier, (p_pde, p_bc, is_neuman_bcs, p_cuboid, box_h)); - p_pde_modifier->SetDependentVariableName("VEGF"); - p_pde_modifier->SetBcsOnBoxBoundary(true); - p_pde_modifier->SetOutputSolutionAtPdeNodes(true); - //p_pde_modifier->SetSolutionInterval(10); // Speed simulations up a bit - //p_pde_modifier->SetOutputGradient(true); - - simulator.AddSimulationModifier(p_pde_modifier); + MAKE_PTR_ARGS(ParabolicPottsPdeModifier<2>, p_pde_modifier, (p_pde, p_bc, is_neuman_bcs, p_cuboid, box_h)); + //MAKE_PTR_ARGS(EllipticPottsPdeModifier<2>, p_pde_modifier, (p_pde, p_bc, is_neuman_bcs, p_cuboid, box_h)); + p_pde_modifier->SetDependentVariableName("VEGF"); + p_pde_modifier->SetBcsOnBoxBoundary(true); + p_pde_modifier->SetOutputSolutionAtPdeNodes(true); + //p_pde_modifier->SetSolutionInterval(10); // Speed simulations up a bit + //p_pde_modifier->SetOutputGradient(true); + + simulator.AddSimulationModifier(p_pde_modifier); simulator.Solve();