Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractGrowingDomainPdeModifier.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 "AbstractGrowingDomainPdeModifier.hpp"
37#include "VertexBasedCellPopulation.hpp"
38#include "MeshBasedCellPopulation.hpp"
39#include "CaBasedCellPopulation.hpp"
40#include "NodeBasedCellPopulation.hpp"
41#include "ReplicatableVector.hpp"
42#include "LinearBasisFunction.hpp"
43
44template <unsigned DIM>
46 boost::shared_ptr<AbstractBoundaryCondition<DIM> > pBoundaryCondition,
47 bool isNeumannBoundaryCondition,
48 Vec solution)
49 : AbstractPdeModifier<DIM>(pPde,
50 pBoundaryCondition,
51 isNeumannBoundaryCondition,
52 solution)
53{
54}
55
56template<unsigned DIM>
60
61template<unsigned DIM>
63{
64 if (this->mDeleteFeMesh)
65 {
66 // If a mesh has been created on a previous time step then we need to tidy it up
67 assert(this->mpFeMesh != nullptr);
68 delete this->mpFeMesh;
69 }
70 else
71 {
73 // This placement assumes that if this->mDeleteFeMesh is false it is uninitialised and needs to
74 // be checked. If true, it has been checked elsewhere.
75 this->mDeleteFeMesh = (dynamic_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation) == nullptr);
76 }
77
78 // Get the finite element mesh via the cell population. Set to NULL first in case mesh generation fails.
79 this->mpFeMesh = nullptr;
80 this->mpFeMesh = rCellPopulation.GetTetrahedralMeshForPdeModifier();
81}
82
83template<unsigned DIM>
85{
86 // Store the PDE solution in an accessible form
87 ReplicatableVector solution_repl(this->mSolution);
88
89 // Local cell index used by the CA simulation
90 unsigned cell_index = 0;
91
92 unsigned index_in_solution_repl = 0;
93 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
94 cell_iter != rCellPopulation.End();
95 ++cell_iter)
96 {
97 unsigned tet_node_index = rCellPopulation.GetLocationIndexUsingCell(*cell_iter);
98
100 if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) != nullptr)
101 {
102 // Offset to relate elements in vertex mesh to nodes in tetrahedral mesh
103 tet_node_index += rCellPopulation.GetNumNodes();
104 }
105 else if (dynamic_cast<CaBasedCellPopulation<DIM>*>(&rCellPopulation) != nullptr)
106 {
107 // Here local cell index corresponds to tet node
108 tet_node_index = cell_index;
109 cell_index++;
110 }
111 else if (dynamic_cast<NodeBasedCellPopulation<DIM>*>(&rCellPopulation) != nullptr)
112 {
113 tet_node_index = index_in_solution_repl;
114 index_in_solution_repl++;
115 }
116
117 double solution_at_node = solution_repl[tet_node_index];
118
119 cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_node);
120
121 if (this->mOutputGradient)
122 {
123 // Now calculate the gradient of the solution and store this in CellVecData
124 c_vector<double, DIM> solution_gradient = zero_vector<double>(DIM);
125
126 Node<DIM>* p_tet_node = this->mpFeMesh->GetNode(tet_node_index);
127
128 // Get the containing elements and average the contribution from each one
129 for (typename Node<DIM>::ContainingElementIterator element_iter = p_tet_node->ContainingElementsBegin();
130 element_iter != p_tet_node->ContainingElementsEnd();
131 ++element_iter)
132 {
133 // Calculate the basis functions at any point (eg zero) in the element
134 c_matrix<double, DIM, DIM> jacobian, inverse_jacobian;
135 double jacobian_det;
136 this->mpFeMesh->GetInverseJacobianForElement(*element_iter, jacobian, jacobian_det, inverse_jacobian);
137 const ChastePoint<DIM> zero_point;
138 c_matrix<double, DIM, DIM+1> grad_phi;
139 LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(zero_point, inverse_jacobian, grad_phi);
140
141 // Add the contribution from this element
142 for (unsigned node_index=0; node_index<DIM+1; node_index++)
143 {
144 double nodal_value = solution_repl[this->mpFeMesh->GetElement(*element_iter)->GetNodeGlobalIndex(node_index)];
145
146 for (unsigned j=0; j<DIM; j++)
147 {
148 solution_gradient(j) += nodal_value* grad_phi(j, node_index);
149 }
150 }
151 }
152
153 // Divide by number of containing elements
154 solution_gradient /= p_tet_node->GetNumContainingElements();
155
156 switch (DIM)
157 {
158 case 1:
159 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
160 break;
161 case 2:
162 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
163 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1));
164 break;
165 case 3:
166 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_x", solution_gradient(0));
167 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_y", solution_gradient(1));
168 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+"_grad_z", solution_gradient(2));
169 break;
170 default:
172 }
173 }
174 }
175}
176
177template<unsigned DIM>
179{
180 // No parameters to output, so just call method on direct parent class
182}
183
184// Explicit instantiation
#define NEVER_REACHED
unsigned GetLocationIndexUsingCell(CellPtr pCell)
virtual TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * GetTetrahedralMeshForPdeModifier()=0
virtual unsigned GetNumNodes()=0
void GenerateFeMesh(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void OutputSimulationModifierParameters(out_stream &rParamsFile)
AbstractGrowingDomainPdeModifier(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, Vec solution=nullptr)
void UpdateCellData(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void OutputSimulationModifierParameters(out_stream &rParamsFile)
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)
Definition Node.hpp:59
ContainingElementIterator ContainingElementsEnd() const
Definition Node.hpp:493
unsigned GetNumContainingElements() const
Definition Node.cpp:312
ContainingElementIterator ContainingElementsBegin() const
Definition Node.hpp:485