Electrodes.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "Electrodes.hpp"
00030 #include "DistributedTetrahedralMesh.hpp"
00031 #include "IsNan.hpp"
00032 #include "HeartConfig.hpp"
00033 
00034 template<unsigned DIM>
00035 Electrodes<DIM>::Electrodes(AbstractTetrahedralMesh<DIM,DIM>& rMesh)
00036     : mpMesh(&rMesh),
00037       mLeftElectrodeArea(0.0),
00038       mRightElectrodeArea(0.0)
00039 {
00040     unsigned axis_index;
00041     double magnitude, duration;
00042     
00043     HeartConfig::Instance()->GetElectrodeParameters(mGroundSecondElectrode, axis_index, magnitude, mStartTime, duration);
00044     
00045     assert(axis_index < DIM);
00046     assert(duration > 0);
00047     mEndTime = mStartTime + duration;
00048     mAreActive = false; // consider electrodes initially switched off!
00049 
00050     ChasteCuboid<DIM> bounding_box = mpMesh->CalculateBoundingBox();
00051     double global_min = bounding_box.rGetLowerCorner()[axis_index];
00052     double global_max = bounding_box.rGetUpperCorner()[axis_index];
00053 
00054     mpBoundaryConditionsContainer.reset(new BoundaryConditionsContainer<DIM,DIM,2>);
00055 
00056 
00057     double input_flux;
00058     double output_flux;
00059 
00060     try
00061     {
00062         ComputeElectrodesAreasAndCheckEquality(axis_index, global_min, global_max);
00063         input_flux = magnitude;
00064         output_flux = -input_flux;
00065 
00066     }
00067     catch (Exception& e)
00068     {
00069         // magnitude of second electrode scaled so that left_area*magnitude_left = -right_area*magnitude_right
00070         input_flux = magnitude;
00071         output_flux = -input_flux*mLeftElectrodeArea/mRightElectrodeArea;
00072 
00073         // Paranoia. In case one of the areas is 0
00074         assert( ! std::isnan(output_flux));
00075         assert( output_flux != 0.0);
00076     }
00077 
00078     ConstBoundaryCondition<DIM>* p_bc_flux_in = new ConstBoundaryCondition<DIM>(input_flux);
00079     ConstBoundaryCondition<DIM>* p_bc_flux_out = new ConstBoundaryCondition<DIM>(output_flux);
00080 
00081     // loop over boundary elements and add a non-zero phi_e boundary condition (ie extracellular
00082     // stimulus) if (assuming axis_index=0, etc) x=lowerValue (where x is the x-value of the centroid)
00083     for (typename AbstractTetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00084             = mpMesh->GetBoundaryElementIteratorBegin();
00085        iter != mpMesh->GetBoundaryElementIteratorEnd();
00086        iter++)
00087     {
00088         if (fabs((*iter)->CalculateCentroid()[axis_index] - global_min) < 1e-6)
00089         {
00090             mpBoundaryConditionsContainer->AddNeumannBoundaryCondition(*iter, p_bc_flux_in,  1);
00091         }
00092 
00093         if (!mGroundSecondElectrode)
00094         {
00095             if (fabs((*iter)->CalculateCentroid()[axis_index] - global_max) < 1e-6)
00096             {
00097                 mpBoundaryConditionsContainer->AddNeumannBoundaryCondition(*iter, p_bc_flux_out, 1);
00098             }
00099         }
00100     }
00101 
00102     // set up mGroundedNodes using opposite surface ie second electrode is
00103     // grounded
00104     if (mGroundSecondElectrode)
00105     {
00106         ConstBoundaryCondition<DIM>* p_zero_bc = new ConstBoundaryCondition<DIM>(0.0);
00107         // We will need to recalculate this when HasDirichletBoundaryConditions() is called.
00108         this->mpBoundaryConditionsContainer->ResetDirichletCommunication();
00109 
00110         for (typename AbstractTetrahedralMesh<DIM,DIM>::NodeIterator iter=mpMesh->GetNodeIteratorBegin();
00111              iter != mpMesh->GetNodeIteratorEnd();
00112              ++iter)
00113         {
00114             if (fabs((*iter).rGetLocation()[axis_index]-global_max) < 1e-6)
00115             {
00116                 mpBoundaryConditionsContainer->AddDirichletBoundaryCondition(&(*iter), p_zero_bc, 1);
00117             }
00118         }
00119 
00120         //Unused boundary conditions will not be deleted by the b.c. container
00121         delete p_bc_flux_out;
00122     }
00123 }
00124 
00125 
00126 
00127 template<unsigned DIM>
00128 boost::shared_ptr<BoundaryConditionsContainer<DIM,DIM,2> > Electrodes<DIM>::GetBoundaryConditionsContainer()
00129 {
00130     //assert(mAreActive);
00131     return mpBoundaryConditionsContainer;
00132 }
00133 
00134 
00135 template<unsigned DIM>
00136 bool Electrodes<DIM>::SwitchOff(double time)
00137 {
00138     // This smidge has to be the same as the one below.
00139     double smidge = 1e-10;
00140     if (mAreActive && time>mEndTime-smidge)
00141     {
00142         mAreActive = false;
00143         return true;
00144     }
00145 
00146     return false;
00147 }
00148 
00149 template<unsigned DIM>
00150 bool Electrodes<DIM>::SwitchOn(double time)
00151 {
00152     // This smidge has to be the same as the one above.
00153     double smidge = 1e-10;
00154     if (!mAreActive && time>=mStartTime && time<=mEndTime - smidge)
00155     {
00156         mAreActive = true;
00157         return true;
00158     }
00159 
00160     return false;
00161 }
00162 
00163 template<unsigned DIM>
00164 void Electrodes<DIM>::ComputeElectrodesAreasAndCheckEquality(unsigned dimensionIndex, double lowerValue, double upperValue)
00165 {
00166     //
00167     // loop over boundary elements and determine areas of the electrode faces
00168     //
00169     double local_left_area = 0.0;
00170     double local_right_area = 0.0;
00171 
00172     c_vector<double,DIM> weighted_direction;
00173     double jacobian_determinant;
00174 
00175 
00176     for (typename AbstractTetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00177              = mpMesh->GetBoundaryElementIteratorBegin();
00178          iter != mpMesh->GetBoundaryElementIteratorEnd();
00179          iter++)
00180     {
00181         if ( mpMesh->CalculateDesignatedOwnershipOfBoundaryElement( (*iter)->GetIndex() ))
00182         {
00183             if (fabs((*iter)->CalculateCentroid()[dimensionIndex] - lowerValue) < 1e-6)
00184             {
00185                 mpMesh->GetWeightedDirectionForBoundaryElement((*iter)->GetIndex(), weighted_direction, jacobian_determinant);
00186                 local_left_area += jacobian_determinant;
00187             }
00188 
00189             if (fabs((*iter)->CalculateCentroid()[dimensionIndex] - upperValue) < 1e-6)
00190             {
00191                 mpMesh->GetWeightedDirectionForBoundaryElement((*iter)->GetIndex(), weighted_direction, jacobian_determinant);
00192                 local_right_area += jacobian_determinant;
00193             }
00194         }
00195     }
00196 
00197     if(DIM==3)
00198     {
00199         // if the dimension of the face is 1, the mapping is from the face to the canonical element [0,1], so the
00200         // jacobian_determinants used above will be exactly the areas of the faces
00201         // if the dimension of the face is 1, the mapping is from the face to the canonical element, the triangle
00202         // with nodes (0,0), (0,1), (1,0), which has area 0.5, so scale the jacobian_determinants by 0.5 to
00203         // get the face areas, ie scale the total area by 0.5:
00204         // (Not technically needed as it is only the ratio of these that is used but since we are calling the variables
00205         // 'area's we ought to be correct).
00206         local_left_area /= 2.0;
00207         local_right_area /= 2.0;
00208     }
00209 
00210     int mpi_ret = MPI_Allreduce(&local_left_area, &mLeftElectrodeArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00211     assert(mpi_ret == MPI_SUCCESS);
00212 
00213     mpi_ret = MPI_Allreduce(&local_right_area, &mRightElectrodeArea, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00214     assert(mpi_ret == MPI_SUCCESS);
00215 
00216     if (mLeftElectrodeArea != mRightElectrodeArea)
00217     {
00218         EXCEPTION("Electrodes have different area");
00219     }
00220 }
00221 
00222 
00224 // Explicit instantiation
00226 
00227 template class Electrodes<1>;
00228 template class Electrodes<2>;
00229 template class Electrodes<3>;
00230 
00231 
00232 // Serialization for Boost >= 1.36
00233 #include "SerializationExportWrapperForCpp.hpp"
00234 EXPORT_TEMPLATE_CLASS_SAME_DIMS(Electrodes)

Generated on Mon Nov 1 12:35:17 2010 for Chaste by  doxygen 1.5.5