Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryLinearMembraneForce.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 "ImmersedBoundaryLinearMembraneForce.hpp"
37
38template <unsigned DIM>
41 mElementSpringConst(1e6),
42 mElementRestLength(0.5),
43 mLaminaSpringConst(1e6),
44 mLaminaRestLength(0.5)
45{
46}
47
48template <unsigned DIM>
52
53template <unsigned DIM>
55 std::vector<std::pair<Node<DIM>*, Node<DIM>*> >& rNodePairs,
57{
58 // Data common across the entire cell population
59 double intrinsic_spacing_squared = rCellPopulation.GetIntrinsicSpacing() * rCellPopulation.GetIntrinsicSpacing();
60
61 // Loop over all elements ( <DIM, DIM> )
62 for (auto elem_it = rCellPopulation.rGetMesh().GetElementIteratorBegin();
63 elem_it != rCellPopulation.rGetMesh().GetElementIteratorEnd();
64 ++elem_it)
65 {
66 CalculateForcesOnElement(*elem_it, rCellPopulation, intrinsic_spacing_squared);
67 }
68
69 // Loop over all laminas ( <DIM-1, DIM> )
70 for (auto lam_it = rCellPopulation.rGetMesh().GetLaminaIteratorBegin();
71 lam_it != rCellPopulation.rGetMesh().GetLaminaIteratorEnd();
72 ++lam_it)
73 {
74 CalculateForcesOnElement(*lam_it, rCellPopulation, intrinsic_spacing_squared);
75 }
76
77 if (this->mAdditiveNormalNoise)
78 {
79 this->AddNormalNoiseToNodes(rCellPopulation);
80 }
81}
82
83template <unsigned DIM>
84template <unsigned ELEMENT_DIM>
88 double intrinsicSpacingSquared)
89{
90 // Get index and number of nodes of current element
91 unsigned elem_idx = rElement.GetIndex();
92 unsigned num_nodes = rElement.GetNumNodes();
93
94 // Make a vector to store the force on node i+1 from node i
95 std::vector<c_vector<double, DIM> > elastic_force_to_next_node(num_nodes);
96
97 /*
98 * Get the node spacing ratio for this element. The rest length and spring
99 * constant are derived from this characteristic length.
100 *
101 * The spring constant is derived with reference to the intrinsic spacing,
102 * so that with different node spacings the user-defined parameters do not
103 * have to be updated.
104 *
105 * The correct factor to increase the spring constant by is (intrinsic
106 * spacing / node_spacing)^2. One factor takes into account the energy
107 * considerations of the elastic springs, and the other takes account of the
108 * factor of node_spacing used in discretising the force relation.
109 */
110 double node_spacing = 0.0;
111 double spring_constant = 0.0;
112 double rest_length = 0.0;
113
114 // Determine if we're in a lamina or not
115 if (ELEMENT_DIM < DIM)
116 {
117 node_spacing = rCellPopulation.rGetMesh().GetAverageNodeSpacingOfLamina(elem_idx, false);
118
119 spring_constant = mLaminaSpringConst * intrinsicSpacingSquared / (node_spacing * node_spacing);
120 rest_length = mLaminaRestLength * node_spacing;
121 }
122 else // regular element
123 {
124 node_spacing = rCellPopulation.rGetMesh().GetAverageNodeSpacingOfElement(elem_idx, false);
125
126 spring_constant = mElementSpringConst * intrinsicSpacingSquared / (node_spacing * node_spacing);
127 rest_length = mElementRestLength * node_spacing;
128 }
129
130 // Loop over nodes and calculate the force exerted on node i+1 by node i
131 for (unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
132 {
133 // Index of the next node, calculated modulo number of nodes in this element
134 unsigned next_idx = (node_idx + 1) % num_nodes;
135
136 double modified_spring_constant = spring_constant;
137 double modified_rest_length = rest_length;
138
139 // Hooke's law linear spring force
140 elastic_force_to_next_node[node_idx] = rCellPopulation.rGetMesh().GetVectorFromAtoB(rElement.GetNodeLocation(node_idx), rElement.GetNodeLocation(next_idx));
141 double normed_dist = norm_2(elastic_force_to_next_node[node_idx]);
142 elastic_force_to_next_node[node_idx] *= modified_spring_constant * (normed_dist - modified_rest_length) / normed_dist;
143 }
144
145 // Add the contributions of springs adjacent to each node
146 for (unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
147 {
148 // Get index of previous node
149 unsigned prev_idx = (node_idx + num_nodes - 1) % num_nodes;
150
151 c_vector<double, DIM> aggregate_force = elastic_force_to_next_node[node_idx] - elastic_force_to_next_node[prev_idx];
152
153 // Add the aggregate force contribution to the node
154 rElement.GetNode(node_idx)->AddAppliedForceContribution(aggregate_force);
155 }
156}
157
158template <unsigned DIM>
160 out_stream& rParamsFile)
161{
162 *rParamsFile << "\t\t\t<ElementSpringConstant>" << mElementSpringConst << "</ElementSpringConstant>\n";
163 *rParamsFile << "\t\t\t<ElementRestLength>" << mElementRestLength << "</ElementRestLength>\n";
164 *rParamsFile << "\t\t\t<LaminaSpringConstant>" << mLaminaSpringConst << "</LaminaSpringConstant>\n";
165 *rParamsFile << "\t\t\t<LaminaRestLength>" << mLaminaRestLength << "</LaminaRestLength>\n";
166
167 // Call method on direct parent class
169}
170
171template <unsigned DIM>
173{
174 return mElementSpringConst;
175}
176
177template <unsigned DIM>
179 double elementSpringConst)
180{
181 mElementSpringConst = elementSpringConst;
182}
183
184template <unsigned DIM>
186{
187 return mElementRestLength;
188}
189
190template <unsigned DIM>
192 double elementRestLength)
193{
194 mElementRestLength = elementRestLength;
195}
196
197template <unsigned DIM>
199{
200 return mLaminaSpringConst;
201}
202
203template <unsigned DIM>
205 double laminaSpringConst)
206{
207 mLaminaSpringConst = laminaSpringConst;
208}
209
210template <unsigned DIM>
212{
213 return mLaminaRestLength;
214}
215
216template <unsigned DIM>
218 double laminaRestLength)
219{
220 mLaminaRestLength = laminaRestLength;
221}
222
223// Explicit instantiation
227
228// 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()
void OutputImmersedBoundaryForceParameters(out_stream &rParamsFile)
void CalculateForcesOnElement(ImmersedBoundaryElement< ELEMENT_DIM, DIM > &rElement, ImmersedBoundaryCellPopulation< DIM > &rCellPopulation, double intrinsicSpacingSquared)
void AddImmersedBoundaryForceContribution(std::vector< std::pair< Node< DIM > *, Node< DIM > * > > &rNodePairs, ImmersedBoundaryCellPopulation< DIM > &rCellPopulation)
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)
Definition Node.hpp:59