Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractBoxDomainPdeModifier.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 "AbstractBoxDomainPdeModifier.hpp"
37#include "ReplicatableVector.hpp"
38#include "LinearBasisFunction.hpp"
39
40template<unsigned DIM>
42 boost::shared_ptr<AbstractBoundaryCondition<DIM> > pBoundaryCondition,
43 bool isNeumannBoundaryCondition,
44 boost::shared_ptr<ChasteCuboid<DIM> > pMeshCuboid,
45 double stepSize,
46 Vec solution)
47 : AbstractPdeModifier<DIM>(pPde,
48 pBoundaryCondition,
49 isNeumannBoundaryCondition,
50 solution),
51 mpMeshCuboid(pMeshCuboid),
52 mStepSize(stepSize),
53 mSetBcsOnBoxBoundary(true)
54{
55 if (pMeshCuboid)
56 {
57 // We only need to generate mpFeMesh once, as it does not vary with time
59 this->mDeleteFeMesh = true;
60 }
61}
62
63template<unsigned DIM>
67
68template<unsigned DIM>
70{
71 return mStepSize;
72}
73
74template<unsigned DIM>
76{
77 mSetBcsOnBoxBoundary = setBcsOnBoxBoundary;
78}
79
80template<unsigned DIM>
82{
83 return mSetBcsOnBoxBoundary;
84}
85
86template<unsigned DIM>
87void AbstractBoxDomainPdeModifier<DIM>::SetupSolve(AbstractCellPopulation<DIM,DIM>& rCellPopulation, std::string outputDirectory)
88{
89 AbstractPdeModifier<DIM>::SetupSolve(rCellPopulation, outputDirectory);
90
91 InitialiseCellPdeElementMap(rCellPopulation);
92}
93
94template<unsigned DIM>
95void AbstractBoxDomainPdeModifier<DIM>::GenerateFeMesh(boost::shared_ptr<ChasteCuboid<DIM> > pMeshCuboid, double stepSize)
96{
97 // Create a regular coarse tetrahedral mesh
98 this->mpFeMesh = new TetrahedralMesh<DIM,DIM>();
99 switch (DIM)
100 {
101 case 1:
102 this->mpFeMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0));
103 break;
104 case 2:
105 this->mpFeMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0), pMeshCuboid->GetWidth(1));
106 break;
107 case 3:
108 this->mpFeMesh->ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0), pMeshCuboid->GetWidth(1), pMeshCuboid->GetWidth(2));
109 break;
110 default:
112 }
113
114 // Get centroid of meshCuboid
115 ChastePoint<DIM> upper = pMeshCuboid->rGetUpperCorner();
116 ChastePoint<DIM> lower = pMeshCuboid->rGetLowerCorner();
117 c_vector<double,DIM> centre_of_cuboid = 0.5*(upper.rGetLocation() + lower.rGetLocation());
118
119 // Find the centre of the PDE mesh
120 c_vector<double,DIM> centre_of_coarse_mesh = zero_vector<double>(DIM);
121 for (unsigned i=0; i<this->mpFeMesh->GetNumNodes(); i++)
122 {
123 centre_of_coarse_mesh += this->mpFeMesh->GetNode(i)->rGetLocation();
124 }
125 centre_of_coarse_mesh /= this->mpFeMesh->GetNumNodes();
126
127 // Now move the mesh to the correct location
128 this->mpFeMesh->Translate(centre_of_cuboid - centre_of_coarse_mesh);
129}
130
131template<unsigned DIM>
133{
134 // Store the PDE solution in an accessible form
135 ReplicatableVector solution_repl(this->mSolution);
136
137 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
138 cell_iter != rCellPopulation.End();
139 ++cell_iter)
140 {
141 // The cells are not nodes of the mesh, so we must interpolate
142 double solution_at_cell = 0.0;
143
144 // Find the element in the FE mesh that contains this cell. CellElementMap has been updated so use this.
145 unsigned elem_index = mCellPdeElementMap[*cell_iter];
146 Element<DIM,DIM>* p_element = this->mpFeMesh->GetElement(elem_index);
147
148 const ChastePoint<DIM>& node_location = rCellPopulation.GetLocationOfCellCentre(*cell_iter);
149
150 c_vector<double,DIM+1> weights = p_element->CalculateInterpolationWeights(node_location);
151
152 for (unsigned i=0; i<DIM+1; i++)
153 {
154 double nodal_value = solution_repl[p_element->GetNodeGlobalIndex(i)];
155 solution_at_cell += nodal_value * weights(i);
156 }
157
158 cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_cell);
159
160 if (this->mOutputGradient)
161 {
162 // Now calculate the gradient of the solution and store this in CellVecData
163 c_vector<double, DIM> solution_gradient = zero_vector<double>(DIM);
164
165 // Calculate the basis functions at any point (e.g. zero) in the element
166 c_matrix<double, DIM, DIM> jacobian, inverse_jacobian;
167 double jacobian_det;
168 this->mpFeMesh->GetInverseJacobianForElement(elem_index, jacobian, jacobian_det, inverse_jacobian);
169 const ChastePoint<DIM> zero_point;
170 c_matrix<double, DIM, DIM+1> grad_phi;
171 LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(zero_point, inverse_jacobian, grad_phi);
172
173 for (unsigned node_index=0; node_index<DIM+1; node_index++)
174 {
175 double nodal_value = solution_repl[p_element->GetNodeGlobalIndex(node_index)];
176
177 for (unsigned j=0; j<DIM; j++)
178 {
179 solution_gradient(j) += nodal_value* grad_phi(j, node_index);
180 }
181 }
182
183 switch (DIM)
184 {
185 case 1:
186 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
187 break;
188 case 2:
189 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
190 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1));
191 break;
192 case 3:
193 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
194 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1));
195 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_z", solution_gradient(2));
196 break;
197 default:
199 }
200 }
201 }
202}
203
204template<unsigned DIM>
206{
207 mCellPdeElementMap.clear();
208
209 // Find the element of mpFeMesh that contains each cell and populate mCellPdeElementMap
210 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
211 cell_iter != rCellPopulation.End();
212 ++cell_iter)
213 {
214 const ChastePoint<DIM>& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter);
215 unsigned elem_index = this->mpFeMesh->GetContainingElementIndex(r_position_of_cell);
216 mCellPdeElementMap[*cell_iter] = elem_index;
217 }
218}
219
220template<unsigned DIM>
222{
223 // Find the element of mpCoarsePdeMesh that contains each cell and populate mCellPdeElementMap
224 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
225 cell_iter != rCellPopulation.End();
226 ++cell_iter)
227 {
228 const ChastePoint<DIM>& r_position_of_cell = rCellPopulation.GetLocationOfCellCentre(*cell_iter);
229 unsigned elem_index = this->mpFeMesh->GetContainingElementIndexWithInitialGuess(r_position_of_cell, mCellPdeElementMap[*cell_iter]);
230 mCellPdeElementMap[*cell_iter] = elem_index;
231 }
232}
233
234template<unsigned DIM>
236{
237 // No parameters to output, so just call method on direct parent class
239}
240
241// Explicit instantiation
#define NEVER_REACHED
void OutputSimulationModifierParameters(out_stream &rParamsFile)
AbstractBoxDomainPdeModifier(boost::shared_ptr< AbstractLinearPde< DIM, DIM > > pPde=boost::shared_ptr< AbstractLinearPde< DIM, DIM > >(), boost::shared_ptr< AbstractBoundaryCondition< DIM > > pBoundaryCondition=boost::shared_ptr< AbstractBoundaryCondition< DIM > >(), bool isNeumannBoundaryCondition=true, boost::shared_ptr< ChasteCuboid< DIM > > pMeshCuboid=boost::shared_ptr< ChasteCuboid< DIM > >(), double stepSize=1.0, Vec solution=nullptr)
virtual void SetupSolve(AbstractCellPopulation< DIM, DIM > &rCellPopulation, std::string outputDirectory)
void GenerateFeMesh(boost::shared_ptr< ChasteCuboid< DIM > > pMeshCuboid, double stepSize)
boost::shared_ptr< ChasteCuboid< DIM > > mpMeshCuboid
void UpdateCellData(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void SetBcsOnBoxBoundary(bool setBcsOnBoxBoundary)
void UpdateCellPdeElementMap(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void InitialiseCellPdeElementMap(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
virtual c_vector< double, SPACE_DIM > GetLocationOfCellCentre(CellPtr pCell)=0
unsigned GetNodeGlobalIndex(unsigned localIndex) const
virtual void SetupSolve(AbstractCellPopulation< DIM, DIM > &rCellPopulation, std::string outputDirectory)
void OutputSimulationModifierParameters(out_stream &rParamsFile)
c_vector< double, DIM > & rGetLocation()
c_vector< double, SPACE_DIM+1 > CalculateInterpolationWeights(const ChastePoint< SPACE_DIM > &rTestPoint)
Definition Element.cpp:224
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM, ELEMENT_DIM+1 > &rReturnValue)