Chaste  Release::2018.1
Electrodes.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "Electrodes.hpp"
37 #include "DistributedTetrahedralMesh.hpp"
38 #include "IsNan.hpp"
39 #include "HeartConfig.hpp"
40 
41 template<unsigned DIM>
43  : mpMesh(&rMesh),
44  mLeftElectrodeArea(0.0),
45  mRightElectrodeArea(0.0)
46 {
47  unsigned axis_index;
48  double magnitude, duration;
49 
51 
52  assert(axis_index < DIM);
53  assert(duration > 0);
54  mEndTime = mStartTime + duration;
55  mAreActive = false; // consider electrodes initially switched off!
56 
58  double global_min = bounding_box.rGetLowerCorner()[axis_index];
59  double global_max = bounding_box.rGetUpperCorner()[axis_index];
60 
62 
63 
64  double input_flux;
65  double output_flux;
66 
67  try
68  {
69  ComputeElectrodesAreasAndCheckEquality(axis_index, global_min, global_max);
70  input_flux = magnitude;
71  output_flux = -input_flux;
72 
73  }
74  catch (Exception&)
75  {
76  // magnitude of second electrode scaled so that left_area*magnitude_left = -right_area*magnitude_right
77  input_flux = magnitude;
78  output_flux = -input_flux*mLeftElectrodeArea/mRightElectrodeArea;
79 
80  // Paranoia. In case one of the areas is 0
81  assert( ! std::isnan(output_flux));
82  assert( output_flux != 0.0);
83  }
84 
85  ConstBoundaryCondition<DIM>* p_bc_flux_in = new ConstBoundaryCondition<DIM>(input_flux);
86  ConstBoundaryCondition<DIM>* p_bc_flux_out = new ConstBoundaryCondition<DIM>(output_flux);
87 
88  // loop over boundary elements and add a non-zero phi_e boundary condition (ie extracellular
89  // stimulus) if (assuming axis_index=0, etc) x=lowerValue (where x is the x-value of the centroid)
93  iter++)
94  {
95  if (fabs((*iter)->CalculateCentroid()[axis_index] - global_min) < 1e-6)
96  {
97  mpBoundaryConditionsContainer->AddNeumannBoundaryCondition(*iter, p_bc_flux_in, 1);
98  }
99 
101  {
102  if (fabs((*iter)->CalculateCentroid()[axis_index] - global_max) < 1e-6)
103  {
104  mpBoundaryConditionsContainer->AddNeumannBoundaryCondition(*iter, p_bc_flux_out, 1);
105  }
106  }
107  }
108 
109  // set up mGroundedNodes using opposite surface ie second electrode is
110  // grounded
112  {
114  // We will need to recalculate this when HasDirichletBoundaryConditions() is called.
115  this->mpBoundaryConditionsContainer->ResetDirichletCommunication();
116 
118  iter != mpMesh->GetNodeIteratorEnd();
119  ++iter)
120  {
121  if (fabs((*iter).rGetLocation()[axis_index]-global_max) < 1e-6)
122  {
123  mpBoundaryConditionsContainer->AddDirichletBoundaryCondition(&(*iter), p_zero_bc, 1);
124  }
125  }
126 
127  //Unused boundary conditions will not be deleted by the b.c. container
128  delete p_bc_flux_out;
129  }
130 }
131 
132 template<unsigned DIM>
133 boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,2> > Electrodes<DIM>::GetBoundaryConditionsContainer()
134 {
135  //assert(mAreActive);
136  return mpBoundaryConditionsContainer;
137 }
138 
139 template<unsigned DIM>
140 bool Electrodes<DIM>::SwitchOff(double time)
141 {
142  // This smidge has to be the same as the one below
144  double smidge = 1e-10;
145  if (mAreActive && time>mEndTime-smidge)
146  {
147  mAreActive = false;
148  return true;
149  }
150 
151  return false;
152 }
153 
154 template<unsigned DIM>
155 bool Electrodes<DIM>::SwitchOn(double time)
156 {
157  // This smidge has to be the same as the one above.
159  double smidge = 1e-10;
160  if (!mAreActive && time>=mStartTime && time<=mEndTime - smidge)
161  {
162  mAreActive = true;
163  return true;
164  }
165 
166  return false;
167 }
168 
169 template<unsigned DIM>
170 void Electrodes<DIM>::ComputeElectrodesAreasAndCheckEquality(unsigned dimensionIndex, double lowerValue, double upperValue)
171 {
172  //
173  // loop over boundary elements and determine areas of the electrode faces
174  //
175  double local_left_area = 0.0;
176  double local_right_area = 0.0;
177 
178  c_vector<double,DIM> weighted_direction;
179  double jacobian_determinant;
180 
183  iter != mpMesh->GetBoundaryElementIteratorEnd();
184  iter++)
185  {
186  if (mpMesh->CalculateDesignatedOwnershipOfBoundaryElement((*iter)->GetIndex()))
187  {
188  if (fabs((*iter)->CalculateCentroid()[dimensionIndex] - lowerValue) < 1e-6)
189  {
190  mpMesh->GetWeightedDirectionForBoundaryElement((*iter)->GetIndex(), weighted_direction, jacobian_determinant);
191  local_left_area += jacobian_determinant;
192  }
193 
194  if (fabs((*iter)->CalculateCentroid()[dimensionIndex] - upperValue) < 1e-6)
195  {
196  mpMesh->GetWeightedDirectionForBoundaryElement((*iter)->GetIndex(), weighted_direction, jacobian_determinant);
197  local_right_area += jacobian_determinant;
198  }
199  }
200  }
201 
202  if (DIM == 3)
203  {
204  // if the dimension of the face is 1, the mapping is from the face to the canonical element [0,1], so the
205  // jacobian_determinants used above will be exactly the areas of the faces
206  // if the dimension of the face is 1, the mapping is from the face to the canonical element, the triangle
207  // with nodes (0,0), (0,1), (1,0), which has area 0.5, so scale the jacobian_determinants by 0.5 to
208  // get the face areas, ie scale the total area by 0.5:
209  // (Not technically needed as it is only the ratio of these that is used but since we are calling the variables
210  // 'area's we ought to be correct).
211  local_left_area /= 2.0;
212  local_right_area /= 2.0;
213  }
214 
215  int mpi_ret = MPI_Allreduce(&local_left_area, &mLeftElectrodeArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
216  UNUSED_OPT(mpi_ret);
217  assert(mpi_ret == MPI_SUCCESS);
218 
219  mpi_ret = MPI_Allreduce(&local_right_area, &mRightElectrodeArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
220  assert(mpi_ret == MPI_SUCCESS);
221 
222  if (mLeftElectrodeArea != mRightElectrodeArea)
223  {
224  EXCEPTION("Electrodes have different area");
225  }
226 }
227 
228 // Explicit instantiation
229 template class Electrodes<1>;
230 template class Electrodes<2>;
231 template class Electrodes<3>;
232 
233 // Serialization for Boost >= 1.36
bool mGroundSecondElectrode
Definition: Electrodes.hpp:72
BoundaryElementIterator GetBoundaryElementIteratorBegin() const
double mEndTime
Definition: Electrodes.hpp:78
#define EXCEPTION(message)
Definition: Exception.hpp:143
double mLeftElectrodeArea
Definition: Electrodes.hpp:89
NodeIterator GetNodeIteratorEnd()
boost::shared_ptr< BoundaryConditionsContainer< DIM, DIM, 2 > > mpBoundaryConditionsContainer
Definition: Electrodes.hpp:74
AbstractTetrahedralMesh< DIM, DIM > * mpMesh
Definition: Electrodes.hpp:86
const ChastePoint< SPACE_DIM > & rGetUpperCorner() const
boost::shared_ptr< BoundaryConditionsContainer< DIM, DIM, 2 > > GetBoundaryConditionsContainer()
Definition: Electrodes.cpp:133
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
bool mAreActive
Definition: Electrodes.hpp:80
void GetElectrodeParameters(bool &rGroundSecondElectrode, unsigned &rIndex, double &rMagnitude, double &rStartTime, double &rDuration)
NodeIterator GetNodeIteratorBegin(bool skipDeletedNodes=true)
#define UNUSED_OPT(var)
Definition: Exception.hpp:216
ChasteCuboid< SPACE_DIM > CalculateBoundingBox(const std::vector< Node< SPACE_DIM > * > &rNodes) const
double mStartTime
Definition: Electrodes.hpp:76
double mRightElectrodeArea
Definition: Electrodes.hpp:92
bool SwitchOn(double time)
Definition: Electrodes.cpp:155
static HeartConfig * Instance()
void ComputeElectrodesAreasAndCheckEquality(unsigned index, double lowerValue, double upperValue)
Definition: Electrodes.cpp:170
bool SwitchOff(double time)
Definition: Electrodes.cpp:140
BoundaryElementIterator GetBoundaryElementIteratorEnd() const
const ChastePoint< SPACE_DIM > & rGetLowerCorner() const