Merge branch 'main' of https://github.com/Chaste/OpenVT

Author: jmosborne

src/Test02MonolayerGrowth/GeneralisedLinearSpringForceWithMinDistanceItem.cpp

@@ -35,6 +35,10 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "GeneralisedLinearSpringForceWithMinDistanceItem.hpp"
 
+#include "AbstractCentreBasedCellPopulation.hpp"
+#include "MeshBasedCellPopulation.hpp"
+#include "NodeBasedCellPopulation.hpp"
+
 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
 GeneralisedLinearSpringForceWithMinDistanceItem<ELEMENT_DIM,SPACE_DIM>::GeneralisedLinearSpringForceWithMinDistanceItem()
    : AbstractTwoBodyInteractionForce<ELEMENT_DIM,SPACE_DIM>(),

src/Test02MonolayerGrowth/QuadraticRepulsionForce.cpp

@@ -35,6 +35,10 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include "QuadraticRepulsionForce.hpp"
 
+#include "AbstractCentreBasedCellPopulation.hpp"
+#include "MeshBasedCellPopulation.hpp"
+#include "NodeBasedCellPopulation.hpp"
+
 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
 QuadraticRepulsionForce<ELEMENT_DIM,SPACE_DIM>::QuadraticRepulsionForce()
    : AbstractTwoBodyInteractionForce<ELEMENT_DIM,SPACE_DIM>(),

src/Test02MonolayerGrowth/SimplifiedNagaiHondaForce.cpp

@@ -0,0 +1,290 @@
+/*
+
+Copyright (c) 2005-2025, 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 "SimplifiedNagaiHondaForce.hpp"
+
+template<unsigned DIM>
+SimplifiedNagaiHondaForce<DIM>::SimplifiedNagaiHondaForce()
+   : AbstractForce<DIM>(),
+     mNagaiHondaDeformationEnergyParameter(100.0), // This is 1.0 in the Nagai & Honda paper.
+     mNagaiHondaMembraneSurfaceEnergyParameter(10.0), // This is 0.1 in the Nagai & Honda paper.
+     mNagaiHondaCellCellAdhesionEnergyParameter(0.5), // This corresponds to a value of 1.0 for
+                                                      // the sigma parameter in the Nagai & Honda
+                                                      // paper. In the paper, the sigma value is
+                                                      // set to 0.01.
+     mNagaiHondaCellBoundaryAdhesionEnergyParameter(1.0), // This is 0.01 in the Nagai & Honda paper
+     mNagaiHondaTargetAreaParameter(0.5*sqrt(3.0)), // This corresponds to a regular hexagon with edge length 0.5
+     mNagaiHondaTargetPerimeterParameter(3.0)
+{
+}
+
+template<unsigned DIM>
+SimplifiedNagaiHondaForce<DIM>::~SimplifiedNagaiHondaForce()
+{
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::AddForceContribution(AbstractCellPopulation<DIM>& rCellPopulation)
+{
+    // Throw an exception message if not using a VertexBasedCellPopulation
+    if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) == nullptr)
+    {
+        EXCEPTION("SimplifiedNagaiHondaForce is to be used with a VertexBasedCellPopulation only");
+    }
+
+    // Define some helper variables
+    VertexBasedCellPopulation<DIM>* p_cell_population = static_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation);
+    unsigned num_nodes = p_cell_population->GetNumNodes();
+    unsigned num_elements = p_cell_population->GetNumElements();
+
+    /*
+     * Check if a subclass of AbstractTargetAreaModifier is being used by
+     * interrogating the first cell (assuming either all cells, or no cells, in
+     * the population have the CellData item "target area").
+     */
+    bool using_target_area_modifier = p_cell_population->Begin()->GetCellData()->HasItem("target area");
+
+    if (using_target_area_modifier)
+    {
+        // This simple model does not support a target area modifier
+        NEVER_REACHED;
+    }
+
+    // Begin by computing the area and perimeter of each element in the mesh, to avoid having to do this multiple times
+    std::vector<double> element_areas(num_elements);
+    std::vector<double> element_perimeters(num_elements);
+    std::vector<double> target_areas(num_elements);
+    std::vector<double> target_perimeters(num_elements);
+
+    for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = p_cell_population->rGetMesh().GetElementIteratorBegin();
+         elem_iter != p_cell_population->rGetMesh().GetElementIteratorEnd();
+         ++elem_iter)
+    {
+        unsigned elem_index = elem_iter->GetIndex();
+        element_areas[elem_index] = p_cell_population->rGetMesh().GetVolumeOfElement(elem_index);
+        element_perimeters[elem_index] = p_cell_population->rGetMesh().GetSurfaceAreaOfElement(elem_index);
+        
+        target_areas[elem_index] = mNagaiHondaTargetAreaParameter;
+        target_perimeters[elem_index] = mNagaiHondaTargetPerimeterParameter;
+    }
+
+    // Iterate over vertices in the cell population
+    for (unsigned node_index=0; node_index<num_nodes; node_index++)
+    {
+        Node<DIM>* p_this_node = p_cell_population->GetNode(node_index);
+
+        /*
+         * The force on this Node is given by the gradient of the total free
+         * energy of the CellPopulation, evaluated at the position of the vertex. This
+         * free energy is the sum of the free energies of all CellPtrs in
+         * the cell population. The free energy of each CellPtr is comprised of three
+         * parts - a cell deformation energy, a membrane surface tension energy
+         * and an adhesion energy.
+         *
+         * Note that since the movement of this Node only affects the free energy
+         * of the CellPtrs containing it, we can just consider the contributions
+         * to the free energy gradient from each of these CellPtrs.
+         */
+        c_vector<double, DIM> deformation_contribution = zero_vector<double>(DIM);
+        c_vector<double, DIM> membrane_surface_tension_contribution = zero_vector<double>(DIM);
+        c_vector<double, DIM> adhesion_contribution = zero_vector<double>(DIM);
+
+        // Find the indices of the elements owned by this node
+        std::set<unsigned> containing_elem_indices = p_cell_population->GetNode(node_index)->rGetContainingElementIndices();
+
+        // Iterate over these elements
+        for (std::set<unsigned>::iterator iter = containing_elem_indices.begin();
+             iter != containing_elem_indices.end();
+             ++iter)
+        {
+            // Get this element, its index and its number of nodes
+            VertexElement<DIM, DIM>* p_element = p_cell_population->GetElement(*iter);
+            unsigned elem_index = p_element->GetIndex();
+            unsigned num_nodes_elem = p_element->GetNumNodes();
+
+            // Find the local index of this node in this element
+            unsigned local_index = p_element->GetNodeLocalIndex(node_index);
+
+            // Add the force contribution from this cell's deformation energy (note the minus sign)
+            c_vector<double, DIM> element_area_gradient = p_cell_population->rGetMesh().GetAreaGradientOfElementAtNode(p_element, local_index);
+            deformation_contribution -= 2*GetNagaiHondaDeformationEnergyParameter()*(element_areas[elem_index] - target_areas[elem_index])*element_area_gradient;
+
+            // Get the previous and next nodes in this element
+            unsigned previous_node_local_index = (num_nodes_elem+local_index-1)%num_nodes_elem;
+            Node<DIM>* p_previous_node = p_element->GetNode(previous_node_local_index);
+
+            unsigned next_node_local_index = (local_index+1)%num_nodes_elem;
+            Node<DIM>* p_next_node = p_element->GetNode(next_node_local_index);
+
+            // Compute the adhesion parameter for each of these edges
+            double previous_edge_adhesion_parameter = GetAdhesionParameter(p_previous_node, p_this_node, *p_cell_population);
+            double next_edge_adhesion_parameter = GetAdhesionParameter(p_this_node, p_next_node, *p_cell_population);
+
+            // Compute the gradient of each these edges, computed at the present node
+            c_vector<double, DIM> previous_edge_gradient = -p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, previous_node_local_index);
+            c_vector<double, DIM> next_edge_gradient = p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, local_index);
+
+            // Add the force contribution from cell-cell and cell-boundary adhesion (note the minus sign)
+            adhesion_contribution -= previous_edge_adhesion_parameter*previous_edge_gradient + next_edge_adhesion_parameter*next_edge_gradient;
+
+            // Add the force contribution from this cell's membrane surface tension (note the minus sign)
+            c_vector<double, DIM> element_perimeter_gradient;
+            element_perimeter_gradient = previous_edge_gradient + next_edge_gradient;
+            //double cell_target_perimeter = 2*sqrt(M_PI*target_areas[elem_index]);
+            membrane_surface_tension_contribution -= 2*GetNagaiHondaMembraneSurfaceEnergyParameter()*(element_perimeters[elem_index] - target_perimeters[elem_index])*element_perimeter_gradient;
+        }
+
+        c_vector<double, DIM> force_on_node = deformation_contribution + membrane_surface_tension_contribution + adhesion_contribution;
+        p_cell_population->GetNode(node_index)->AddAppliedForceContribution(force_on_node);
+    }
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetAdhesionParameter(Node<DIM>* pNodeA, Node<DIM>* pNodeB, VertexBasedCellPopulation<DIM>& rVertexCellPopulation)
+{
+    // Find the indices of the elements owned by each node
+    std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
+    std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
+
+    // Find common elements
+    std::set<unsigned> shared_elements;
+    std::set_intersection(elements_containing_nodeA.begin(),
+                          elements_containing_nodeA.end(),
+                          elements_containing_nodeB.begin(),
+                          elements_containing_nodeB.end(),
+                          std::inserter(shared_elements, shared_elements.begin()));
+
+    // Check that the nodes have a common edge
+    assert(!shared_elements.empty());
+
+    double adhesion_parameter = GetNagaiHondaCellCellAdhesionEnergyParameter();
+
+    // If the edge corresponds to a single element, then the cell is on the boundary
+    if (shared_elements.size() == 1)
+    {
+        adhesion_parameter = GetNagaiHondaCellBoundaryAdhesionEnergyParameter();
+    }
+
+    return adhesion_parameter;
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetNagaiHondaDeformationEnergyParameter()
+{
+    return mNagaiHondaDeformationEnergyParameter;
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetNagaiHondaMembraneSurfaceEnergyParameter()
+{
+    return mNagaiHondaMembraneSurfaceEnergyParameter;
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetNagaiHondaCellCellAdhesionEnergyParameter()
+{
+    return mNagaiHondaCellCellAdhesionEnergyParameter;
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetNagaiHondaCellBoundaryAdhesionEnergyParameter()
+{
+    return mNagaiHondaCellBoundaryAdhesionEnergyParameter;
+}
+
+template<unsigned DIM>
+double SimplifiedNagaiHondaForce<DIM>::GetNagaiHondaTargetAreaParameter()
+{
+    return mNagaiHondaTargetAreaParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaDeformationEnergyParameter(double deformationEnergyParameter)
+{
+    mNagaiHondaDeformationEnergyParameter = deformationEnergyParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaMembraneSurfaceEnergyParameter(double membraneSurfaceEnergyParameter)
+{
+    mNagaiHondaMembraneSurfaceEnergyParameter = membraneSurfaceEnergyParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaCellCellAdhesionEnergyParameter(double cellCellAdhesionEnergyParameter)
+{
+    mNagaiHondaCellCellAdhesionEnergyParameter = cellCellAdhesionEnergyParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaCellBoundaryAdhesionEnergyParameter(double cellBoundaryAdhesionEnergyParameter)
+{
+    mNagaiHondaCellBoundaryAdhesionEnergyParameter = cellBoundaryAdhesionEnergyParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaTargetAreaParameter(double nagaiHondaTargetAreaParameter)
+{
+    mNagaiHondaTargetAreaParameter = nagaiHondaTargetAreaParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::SetNagaiHondaTargetPerimeterParameter(double nagaiHondaTargetPerimeterParameter)
+{
+    mNagaiHondaTargetPerimeterParameter = nagaiHondaTargetPerimeterParameter;
+}
+
+template<unsigned DIM>
+void SimplifiedNagaiHondaForce<DIM>::OutputForceParameters(out_stream& rParamsFile)
+{
+    *rParamsFile << "\t\t\t<NagaiHondaDeformationEnergyParameter>" << mNagaiHondaDeformationEnergyParameter << "</NagaiHondaDeformationEnergyParameter>\n";
+    *rParamsFile << "\t\t\t<NagaiHondaMembraneSurfaceEnergyParameter>" << mNagaiHondaMembraneSurfaceEnergyParameter << "</NagaiHondaMembraneSurfaceEnergyParameter>\n";
+    *rParamsFile << "\t\t\t<NagaiHondaCellCellAdhesionEnergyParameter>" << mNagaiHondaCellCellAdhesionEnergyParameter << "</NagaiHondaCellCellAdhesionEnergyParameter>\n";
+    *rParamsFile << "\t\t\t<NagaiHondaCellBoundaryAdhesionEnergyParameter>" << mNagaiHondaCellBoundaryAdhesionEnergyParameter << "</NagaiHondaCellBoundaryAdhesionEnergyParameter>\n";
+    *rParamsFile << "\t\t\t<NagaiHondaTargetAreaParameter>" << mNagaiHondaTargetAreaParameter << "</NagaiHondaTargetAreaParameter>\n";
+    *rParamsFile << "\t\t\t<NagaiHondaTargetPerimeterParameter>" << mNagaiHondaTargetPerimeterParameter << "</NagaiHondaTargetPerimeterParameter>\n";
+
+    // Call method on direct parent class
+    AbstractForce<DIM>::OutputForceParameters(rParamsFile);
+}
+
+// Explicit instantiation
+template class SimplifiedNagaiHondaForce<1>;
+template class SimplifiedNagaiHondaForce<2>;
+template class SimplifiedNagaiHondaForce<3>;
+
+// Serialization for Boost >= 1.36
+#include "SerializationExportWrapperForCpp.hpp"
+EXPORT_TEMPLATE_CLASS_SAME_DIMS(SimplifiedNagaiHondaForce)