Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryMorseMembraneForce.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "ImmersedBoundaryMorseMembraneForce.hpp"
37
38template <unsigned DIM>
41 mElementWellDepth(1e6),
42 mElementRestLength(0.5),
43 mLaminaWellDepth(1e6),
44 mLaminaRestLength(0.5),
45 mWellWidth(0.25)
46{
47}
48
49template <unsigned DIM>
53
54template <unsigned DIM>
56 std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs,
58{
59 // Data common across the entire cell population
60 double intrinsic_spacing_squared = rCellPopulation.GetIntrinsicSpacing() * rCellPopulation.GetIntrinsicSpacing();
61
62 // Loop over all elements ( <DIM, DIM> )
63 for (auto elem_it = rCellPopulation.rGetMesh().GetElementIteratorBegin();
64 elem_it != rCellPopulation.rGetMesh().GetElementIteratorEnd();
65 ++elem_it)
66 {
67 CalculateForcesOnElement(*elem_it, rCellPopulation, intrinsic_spacing_squared);
68 }
69
70 // Loop over all laminas ( <DIM-1, DIM> )
71 for (auto lam_it = rCellPopulation.rGetMesh().GetLaminaIteratorBegin();
72 lam_it != rCellPopulation.rGetMesh().GetLaminaIteratorEnd();
73 ++lam_it)
74 {
75 CalculateForcesOnElement(*lam_it, rCellPopulation, intrinsic_spacing_squared);
76 }
77
78 if (this->mAdditiveNormalNoise)
79 {
80 this->AddNormalNoiseToNodes(rCellPopulation);
81 }
82}
83
84template <unsigned DIM>
85template <unsigned ELEMENT_DIM>
89 double intrinsicSpacingSquared)
90{
91 // Get index and number of nodes of current element
92 unsigned elem_idx = rElement.GetIndex();
93 unsigned num_nodes = rElement.GetNumNodes();
94
95 // Make a vector to store the force on node i+1 from node i
96 std::vector<c_vector<double, DIM> > force_to_next(num_nodes);
97
98 /*
99 * Get the node spacing ratio for this element. The rest length and spring
100 * constant are derived from this characteristic length.
101 *
102 * The spring constant is derived with reference to the intrinsic spacing,
103 * so that with different node spacings the user-defined parameters do not
104 * have to be updated.
105 *
106 * The correct factor to increase the spring constant by is (intrinsic
107 * spacing / node_spacing)^2. One factor takes into account the energy
108 * considerations of the elastic springs, and the other takes account of the
109 * factor of node_spacing used in discretising the force relation.
110 */
111 double node_spacing = 0.0;
112 double well_depth = 0.0;
113 double rest_length = 0.0;
114 double well_width = mWellWidth;
115
116 // Determine if we're in a lamina or not
117 if (ELEMENT_DIM < DIM)
118 {
119 node_spacing = rCellPopulation.rGetMesh().GetAverageNodeSpacingOfLamina(elem_idx, false);
120
121 well_depth = mLaminaWellDepth * intrinsicSpacingSquared / (node_spacing * node_spacing);
122 rest_length = mLaminaRestLength * node_spacing;
123 well_width *= node_spacing;
124 }
125 else // regular element
126 {
127 node_spacing = rCellPopulation.rGetMesh().GetAverageNodeSpacingOfElement(elem_idx, false);
128
129 well_depth = mElementWellDepth * intrinsicSpacingSquared / (node_spacing * node_spacing);
130 rest_length = mElementRestLength * node_spacing;
131 well_width *= node_spacing;
132 }
133
134 // Loop over nodes and calculate the force exerted on node i+1 by node i
135 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
136 {
137 // Index of the next node, calculated modulo number of nodes in this element
138 unsigned next_idx = (node_idx + 1) % num_nodes;
139
140 // Morse force (derivative of Morse potential wrt distance between nodes
141 force_to_next[node_idx] = rCellPopulation.rGetMesh().GetVectorFromAtoB(rElement.GetNodeLocation(node_idx),
142 rElement.GetNodeLocation(next_idx));
143 double normed_dist = norm_2(force_to_next[node_idx]);
144
145 double morse_exp = exp((rest_length - normed_dist) / well_width);
146 force_to_next[node_idx] *= 2.0 * well_width * well_depth * morse_exp * (1.0 - morse_exp) / normed_dist;
147 }
148
149 // Add the contributions of springs adjacent to each node
150 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
151 {
152 // Get index of previous node
153 unsigned prev_idx = (node_idx + num_nodes - 1) % num_nodes;
154
155 c_vector<double, DIM> aggregate_force = force_to_next[node_idx] - force_to_next[prev_idx];
156
157 // Add the aggregate force contribution to the node
158 rElement.GetNode(node_idx)->AddAppliedForceContribution(aggregate_force);
159 }
160}
161
162template <unsigned DIM>
164 out_stream& rParamsFile)
165{
166 *rParamsFile << "\t\t\t<ElementWellDepth>" << mElementWellDepth << "</ElementWellDepth>\n";
167 *rParamsFile << "\t\t\t<ElementRestLength>" << mElementRestLength << "</ElementRestLength>\n";
168 *rParamsFile << "\t\t\t<LaminaWellDepth>" << mLaminaWellDepth << "</LaminaWellDepth>\n";
169 *rParamsFile << "\t\t\t<LaminaRestLength>" << mLaminaRestLength << "</LaminaRestLength>\n";
170 *rParamsFile << "\t\t\t<WellWidth>" << mWellWidth << "</WellWidth>\n";
171
172 // Call method on direct parent class
174}
175
176template <unsigned DIM>
178{
179 return mElementWellDepth;
180}
181
182template <unsigned DIM>
184 double elementWellDepth)
185{
186 mElementWellDepth = elementWellDepth;
187}
188
189template <unsigned DIM>
191{
192 return mElementRestLength;
193}
194
195template <unsigned DIM>
197 double elementRestLength)
198{
199 mElementRestLength = elementRestLength;
200}
201
202template <unsigned DIM>
204{
205 return mLaminaWellDepth;
206}
207
208template <unsigned DIM>
210 double laminaWellDepth)
211{
212 mLaminaWellDepth = laminaWellDepth;
213}
214
215template <unsigned DIM>
217{
218 return mLaminaRestLength;
219}
220
221template <unsigned DIM>
223 double laminaRestLength)
224{
225 mLaminaRestLength = laminaRestLength;
226}
227
228template <unsigned DIM>
230{
231 return mWellWidth;
232}
233
234template <unsigned DIM>
236{
237 mWellWidth = wellWidth;
238}
239
240// Explicit instantiation
244
245// Serialization for Boost >= 1.36
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
unsigned GetNumNodes() const
unsigned GetIndex() const
virtual void OutputImmersedBoundaryForceParameters(out_stream &rParamsFile)=0
ImmersedBoundaryMesh< DIM, DIM > & rGetMesh()
ImmersedBoundaryElementIterator GetElementIteratorEnd()
ImmersedBoundaryLaminaIterator GetLaminaIteratorBegin(bool skipDeletedLaminas=true)
ImmersedBoundaryElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
ImmersedBoundaryLaminaIterator GetLaminaIteratorEnd()
double GetAverageNodeSpacingOfElement(unsigned index, bool recalculate=true)
double GetAverageNodeSpacingOfLamina(unsigned index, bool recalculate=true)
c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocation1, const c_vector< double, SPACE_DIM > &rLocation2)
void CalculateForcesOnElement(ImmersedBoundaryElement< ELEMENT_DIM, DIM > &rElement, ImmersedBoundaryCellPopulation< DIM > &rCellPopulation, double intrinsicSpacingSquared)
void OutputImmersedBoundaryForceParameters(out_stream &rParamsFile)
void AddImmersedBoundaryForceContribution(std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs, ImmersedBoundaryCellPopulation< DIM > &rCellPopulation)
Definition Node.hpp:59