Chaste Commit::30a3e656d4b131f8c595cc6eb2becd297337570f
AveragedSourceParabolicPde.cpp
1/*
2
3Copyright (c) 2005-2025, 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 "AveragedSourceParabolicPde.hpp"
37#include "ApoptoticCellProperty.hpp"
38
39template<unsigned DIM>
41 double constantSourceCoefficient,
42 double linearSourceCoefficient,
43 double diffusionCoefficient,
44 double duDtCoefficient,
45 bool scaleByCellVolume)
46 : mrCellPopulation(rCellPopulation),
47 mConstantSourceCoefficient(constantSourceCoefficient),
48 mLinearSourceCoefficient(linearSourceCoefficient),
49 mDiffusionCoefficient(diffusionCoefficient),
50 mDuDtCoefficient(duDtCoefficient),
51 mScaleByCellVolume(scaleByCellVolume)
52{
53}
54
55template<unsigned DIM>
60
61template<unsigned DIM>
63{
64 return mConstantSourceCoefficient;
65}
66
67template<unsigned DIM>
69{
70 return mLinearSourceCoefficient;
71}
72
73template<unsigned DIM>
75{
76 return mDiffusionCoefficient;
77}
78
79template<unsigned DIM>
81{
82 return mDuDtCoefficient;
83}
84
85template<unsigned DIM>
87{
88 return mScaleByCellVolume;
89}
90
91template<unsigned DIM>
92void AveragedSourceParabolicPde<DIM>::SetupSourceTerms(TetrahedralMesh<DIM,DIM>& rCoarseMesh, std::map<CellPtr, unsigned>* pCellPdeElementMap) // must be called before solve
93{
94 // Allocate memory
95 mCellDensityOnCoarseElements.resize(rCoarseMesh.GetNumElements());
96 for (double & mCellDensityOnCoarseElement : mCellDensityOnCoarseElements)
97 {
98 mCellDensityOnCoarseElement = 0.0;
99 }
100
101 // Loop over cells, find which coarse element it is in, and add 1 to mSourceTermOnCoarseElements[elem_index]
102 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mrCellPopulation.Begin();
103 cell_iter != mrCellPopulation.End();
104 ++cell_iter)
105 {
106 unsigned elem_index = 0;
107 const ChastePoint<DIM>& r_position_of_cell = mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
108
109 if (pCellPdeElementMap != nullptr)
110 {
111 elem_index = (*pCellPdeElementMap)[*cell_iter];
112 }
113 else
114 {
115 elem_index = rCoarseMesh.GetContainingElementIndex(r_position_of_cell);
116 }
117
118 // Update element map if cell has moved
119 bool cell_is_apoptotic = cell_iter->template HasCellProperty<ApoptoticCellProperty>();
120
121 if (!cell_is_apoptotic)
122 {
123 double cell_weight = 1.0;
124
125 if (mScaleByCellVolume)
126 {
127 // If scaling by cell volume then use volume her instead of cell count
128 cell_weight = mrCellPopulation.GetVolumeOfCell(*cell_iter);
129
130 if (cell_weight <1e-6)
131 {
132 EXCEPTION("The volume of one of the cells is " << cell_weight <<
133 " and you are scaling by cell volume. Either turn scaling off or use"
134 " a cell model with non zero areas (i.e. a Bounded Voronoi Tesselation model).");
135 }
136 }
137
138 mCellDensityOnCoarseElements[elem_index] += cell_weight;
139 }
140 }
141
142 // Then divide each entry of mSourceTermOnCoarseElements by the element's area
143 c_matrix<double, DIM, DIM> jacobian;
144 double det;
145 for (unsigned elem_index=0; elem_index<mCellDensityOnCoarseElements.size(); elem_index++)
146 {
147 rCoarseMesh.GetElement(elem_index)->CalculateJacobian(jacobian, det);
148 mCellDensityOnCoarseElements[elem_index] /= rCoarseMesh.GetElement(elem_index)->GetVolume(det);
149 }
150}
151
152template<unsigned DIM>
154{
155 return mDuDtCoefficient;
156}
157
158template<unsigned DIM>
160{
161 assert(!mCellDensityOnCoarseElements.empty());
162
163 // The source term is (a*u + b)*density
164 double source_term = (mLinearSourceCoefficient * u + mConstantSourceCoefficient) * mCellDensityOnCoarseElements[pElement->GetIndex()];
165
166 return source_term;
167}
168
169// LCOV_EXCL_START
170template<unsigned DIM>
175// LCOV_EXCL_STOP
176
177template<unsigned DIM>
179{
180 return mDiffusionCoefficient*identity_matrix<double>(DIM);
181}
182
183template<unsigned DIM>
185{
186 return this->mCellDensityOnCoarseElements[elementIndex];
187}
188
189// Explicit instantiation
190template class AveragedSourceParabolicPde<1>;
191template class AveragedSourceParabolicPde<2>;
192template class AveragedSourceParabolicPde<3>;
193
194// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
unsigned GetIndex() const
unsigned GetContainingElementIndex(const ChastePoint< SPACE_DIM > &rTestPoint, bool strict=false, std::set< unsigned > testElements=std::set< unsigned >(), bool onlyTryWithTestElements=false)
Element< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
virtual unsigned GetNumElements() const
double ComputeSourceTerm(const ChastePoint< DIM > &rX, double u, Element< DIM, DIM > *pElement) override
double ComputeDuDtCoefficientFunction(const ChastePoint< DIM > &rX) override
c_matrix< double, DIM, DIM > ComputeDiffusionTerm(const ChastePoint< DIM > &rX, Element< DIM, DIM > *pElement=nullptr) override
AveragedSourceParabolicPde(AbstractCellPopulation< DIM, DIM > &rCellPopulation, double constantSourceCoefficient=0.0, double linearSourceCoefficient=0.0, double diffusionCoefficient=1.0, double duDtCoefficient=1.0, bool scaleByCellVolume=false)
double GetUptakeRateForElement(unsigned elementIndex)
double ComputeSourceTermAtNode(const Node< DIM > &rNode, double u) override
const AbstractCellPopulation< DIM > & rGetCellPopulation() const
virtual void SetupSourceTerms(TetrahedralMesh< DIM, DIM > &rCoarseMesh, std::map< CellPtr, unsigned > *pCellPdeElementMap=nullptr)
Definition Node.hpp:59