BoundaryConditionsContainerImplementation.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 #ifndef _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_
00030 #define _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_
00031 
00032 #include "PetscVecTools.hpp"
00033 #include "PetscMatTools.hpp"
00034 #include "BoundaryConditionsContainer.hpp"
00035 #include "DistributedVector.hpp"
00036 #include "ReplicatableVector.hpp"
00037 #include "ConstBoundaryCondition.hpp"
00038 
00039 #include "HeartEventHandler.hpp"
00040 
00041 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00042 BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::BoundaryConditionsContainer(bool deleteConditions)
00043             : AbstractBoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>(deleteConditions)
00044 {
00045     mLoadedFromArchive = false;
00046 
00047     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00048     {
00049         mpNeumannMap[index_of_unknown] = new std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *, const AbstractBoundaryCondition<SPACE_DIM>*>;
00050 
00051         mAnyNonZeroNeumannConditionsForUnknown[index_of_unknown] = false;
00052         mLastNeumannCondition[index_of_unknown] = mpNeumannMap[index_of_unknown]->begin();
00053     }
00054 
00055     // This zero boundary condition is only used in AddNeumannBoundaryCondition
00056     mpZeroBoundaryCondition = new ConstBoundaryCondition<SPACE_DIM>(0.0);
00057 }
00058 
00059 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00060 BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::~BoundaryConditionsContainer()
00061 {
00062 
00063     // Keep track of what boundary condition objects we've deleted
00064     std::set<const AbstractBoundaryCondition<SPACE_DIM>*> deleted_conditions;
00065     for (unsigned i=0; i<PROBLEM_DIM; i++)
00066     {
00067         NeumannMapIterator neumann_iterator = mpNeumannMap[i]->begin();
00068         while (neumann_iterator != mpNeumannMap[i]->end() )
00069         {
00070 
00071             if (deleted_conditions.count(neumann_iterator->second) == 0)
00072             {
00073                 deleted_conditions.insert(neumann_iterator->second);
00074                 //Leave the zero boundary condition until last
00075                 if (neumann_iterator->second != mpZeroBoundaryCondition)
00076                 {
00077                     if (this->mDeleteConditions)
00078                     {
00079                         delete neumann_iterator->second;
00080                     }
00081                 }
00082             }
00083             neumann_iterator++;
00084         }
00085         delete(mpNeumannMap[i]);
00086     }
00087 
00088     delete mpZeroBoundaryCondition;
00089 
00090     if (this->mDeleteConditions)
00091     {
00092         this->DeleteDirichletBoundaryConditions(deleted_conditions);
00093     }
00094 }
00095 
00096 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00097 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AddDirichletBoundaryCondition( const Node<SPACE_DIM> *  pBoundaryNode,
00098                                         const AbstractBoundaryCondition<SPACE_DIM> * pBoundaryCondition,
00099                                         unsigned indexOfUnknown,
00100                                         bool checkIfBoundaryNode)
00101 {
00102     assert(indexOfUnknown < PROBLEM_DIM);
00103     if (checkIfBoundaryNode)
00104     {
00105         assert( pBoundaryNode->IsBoundaryNode());
00106     }
00107 
00108     (*(this->mpDirichletMap[indexOfUnknown]))[pBoundaryNode] = pBoundaryCondition;
00109 }
00110 
00111 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00112 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AddNeumannBoundaryCondition( const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> * pBoundaryElement,
00113                                       const AbstractBoundaryCondition<SPACE_DIM> * pBoundaryCondition,
00114                                       unsigned indexOfUnknown)
00115 {
00116     assert(indexOfUnknown < PROBLEM_DIM);
00117 
00118     // If this condition is constant, we can test whether it is zero.
00119     // Otherwise we assume that this could be a non-zero boundary condition.
00120     const ConstBoundaryCondition<SPACE_DIM>* p_const_cond = dynamic_cast<const ConstBoundaryCondition<SPACE_DIM>*>(pBoundaryCondition);
00121     if (p_const_cond)
00122     {
00123         if (p_const_cond->GetValue(pBoundaryElement->GetNode(0)->GetPoint()) != 0.0)
00124         {
00125             mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = true;
00126         }
00127     }
00128     else
00129     {
00130         mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = true;
00131     }
00132 
00133     for (unsigned unknown=0; unknown<PROBLEM_DIM; unknown++)
00134     {
00135         if (unknown==indexOfUnknown)
00136         {
00137             (*(mpNeumannMap[indexOfUnknown]))[pBoundaryElement] = pBoundaryCondition;
00138         }
00139         else
00140         {
00141             // if can't find pBoundaryElement in map[unknown]
00142             if ( mpNeumannMap[unknown]->find(pBoundaryElement)==mpNeumannMap[unknown]->end() )
00143             {
00144                 // add zero bc to other unknowns (so all maps are in sync)
00145                 (*(mpNeumannMap[unknown]))[pBoundaryElement] = mpZeroBoundaryCondition;
00146             }
00147         }
00148     }
00149 }
00150 
00151 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00152 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::DefineZeroDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00153                                            unsigned indexOfUnknown)
00154 {
00155     this->DefineConstantDirichletOnMeshBoundary(pMesh, 0.0, indexOfUnknown);
00156 }
00157 
00158 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00159 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::DefineConstantDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00160                                                double value,
00161                                                unsigned indexOfUnknown)
00162 {
00163     assert(indexOfUnknown < PROBLEM_DIM);
00164     //In applying a condition to the boundary, we need to be sure that the boundary exists
00165     assert(pMesh->GetNumBoundaryNodes() > 0);
00166 
00167     ConstBoundaryCondition<SPACE_DIM>* p_boundary_condition =
00168         new ConstBoundaryCondition<SPACE_DIM>( value );
00169 
00170     typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryNodeIterator iter;
00171     iter = pMesh->GetBoundaryNodeIteratorBegin();
00172     while (iter != pMesh->GetBoundaryNodeIteratorEnd())
00173     {
00174         AddDirichletBoundaryCondition(*iter, p_boundary_condition, indexOfUnknown);
00175         iter++;
00176     }
00177 }
00178 
00179 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00180 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::DefineZeroNeumannOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00181                                          unsigned indexOfUnknown)
00182 {
00183     assert(indexOfUnknown < PROBLEM_DIM);
00184     //In applying a condition to the boundary, we need to be sure that the boundary exists
00185     assert(pMesh->GetNumBoundaryElements() > 0);
00186     ConstBoundaryCondition<SPACE_DIM>* p_zero_boundary_condition =
00187         new ConstBoundaryCondition<SPACE_DIM>( 0.0 );
00188 
00189     typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryElementIterator iter;
00190     iter = pMesh->GetBoundaryElementIteratorBegin();
00191     while (iter != pMesh->GetBoundaryElementIteratorEnd())
00192     {
00193         AddNeumannBoundaryCondition(*iter, p_zero_boundary_condition, indexOfUnknown);
00194         iter++;
00195     }
00196 
00197     mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = false;
00198 }
00229 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00230 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ApplyDirichletToLinearProblem(
00231         LinearSystem& rLinearSystem,
00232         bool applyToMatrix,
00233         bool applyToRhsVector)
00234 {
00235     HeartEventHandler::BeginEvent(HeartEventHandler::DIRICHLET_BCS);
00236 
00237     if (applyToMatrix)
00238     {
00239         if (!this->HasDirichletBoundaryConditions())
00240         {
00241             // Short-circuit the replication if there are no conditions
00242             HeartEventHandler::EndEvent(HeartEventHandler::DIRICHLET_BCS);
00243             return;
00244         }
00245 
00246         bool matrix_is_symmetric = rLinearSystem.IsMatrixSymmetric();
00247 
00248         if (matrix_is_symmetric)
00249         {
00250             //Modifications to the RHS are stored in the Dirichlet boundary conditions vector. This is done so
00251             //that they can be reapplied at each time step.
00252             //Make a new vector to store the Dirichlet offsets in
00253             VecDuplicate(rLinearSystem.rGetRhsVector(), &(rLinearSystem.rGetDirichletBoundaryConditionsVector()));
00254             PetscVecTools::Zero(rLinearSystem.rGetDirichletBoundaryConditionsVector());
00255             // If the matrix is symmetric, calls to GetMatrixRowDistributed() require the matrix to be
00256             // in assembled state. Otherwise we can defer it.
00257             rLinearSystem.AssembleFinalLinearSystem();
00258         }
00259 
00260         // Work out where we're setting dirichlet boundary conditions *everywhere*, not just those locally known
00261         ReplicatableVector dirichlet_conditions(rLinearSystem.GetSize());
00262         unsigned lo, hi;
00263         {
00264             PetscInt lo_s, hi_s;
00265             rLinearSystem.GetOwnershipRange(lo_s, hi_s);
00266             lo = lo_s; hi = hi_s;
00267         }
00268         // Initialise all local entries to DBL_MAX, i.e. don't know if there's a condition
00269         for (unsigned i=lo; i<hi; i++)
00270         {
00271             dirichlet_conditions[i] = DBL_MAX;
00272         }
00273         // Now fill in the ones we know
00274         for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00275         {
00276             this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
00277 
00278             while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
00279             {
00280                 unsigned node_index = this->mDirichIterator->first->GetIndex();
00281                 double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
00282                 assert(value != DBL_MAX);
00283 
00284                 unsigned row = PROBLEM_DIM*node_index + index_of_unknown; // assumes vm and phie equations are interleaved
00285                 dirichlet_conditions[row] = value;
00286 
00287                 this->mDirichIterator++;
00288             }
00289         }
00290         // And replicate
00291         dirichlet_conditions.Replicate(lo, hi);
00292 
00293         // Which rows have conditions?
00294         std::vector<unsigned> rows_to_zero;
00295         for (unsigned i=0; i<dirichlet_conditions.GetSize(); i++)
00296         {
00297             if (dirichlet_conditions[i] != DBL_MAX)
00298             {
00299                 rows_to_zero.push_back(i);
00300             }
00301         }
00302 
00303         if (matrix_is_symmetric)
00304         {
00305             // Modify the matrix columns
00306             for (unsigned i=0; i<rows_to_zero.size(); i++)
00307             {
00308                 unsigned col = rows_to_zero[i];
00309                 double minus_value = -dirichlet_conditions[col];
00310 
00311                 // Get a vector which will store the column of the matrix (column d, where d is
00312                 // the index of the row (and column) to be altered for the boundary condition.
00313                 // Since the matrix is symmetric when get row number "col" and treat it as a column.
00314                 // PETSc uses compressed row format and therefore getting rows is far more efficient
00315                 // than getting columns.
00316                 Vec matrix_col = rLinearSystem.GetMatrixRowDistributed(col);
00317 
00318                 // Zero the correct entry of the column
00319                 PetscVecTools::SetElement(matrix_col, col, 0.0);
00320 
00321                 // Set up the RHS Dirichlet boundary conditions vector
00322                 // Assuming one boundary at the zeroth node (x_0 = value), this is equal to
00323                 //   -value*[0 a_21 a_31 .. a_N1]
00324                 // and will be added to the RHS.
00325                 PetscVecTools::AddScaledVector(rLinearSystem.rGetDirichletBoundaryConditionsVector(), matrix_col, minus_value);
00326                 VecDestroy(matrix_col);
00327             }
00328         }
00329 
00330         // Now zero the appropriate rows and columns of the matrix. If the matrix is symmetric we apply the
00331         // boundary conditions in a way the symmetry isn't lost (rows and columns). If not only the row is
00332         // zeroed.
00333         if (matrix_is_symmetric)
00334         {
00335             rLinearSystem.ZeroMatrixRowsAndColumnsWithValueOnDiagonal(rows_to_zero, 1.0);
00336         }
00337         else
00338         {
00339             rLinearSystem.ZeroMatrixRowsWithValueOnDiagonal(rows_to_zero, 1.0);
00340         }
00341     }
00342 
00343     if(applyToRhsVector)
00344     {
00345         //Apply the RHS boundary conditions modification if required.
00346         if (rLinearSystem.rGetDirichletBoundaryConditionsVector())
00347         {
00348             PetscVecTools::AddScaledVector(rLinearSystem.rGetRhsVector(), rLinearSystem.rGetDirichletBoundaryConditionsVector(), 1.0);
00349         }
00350     
00351         //Apply the actual boundary condition to the RHS, note this must be done after the modification to the
00352         //RHS vector.
00353         for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00354         {
00355             this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
00356     
00357             while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
00358             {
00359                 unsigned node_index = this->mDirichIterator->first->GetIndex();
00360                 double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
00361     
00362                 unsigned row = PROBLEM_DIM*node_index + index_of_unknown;
00363     
00364                 rLinearSystem.SetRhsVectorElement(row, value);
00365     
00366                 this->mDirichIterator++;
00367             }
00368         }
00369     }
00370 
00371     HeartEventHandler::EndEvent(HeartEventHandler::DIRICHLET_BCS);
00372 }
00373 
00374 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00375 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ApplyDirichletToNonlinearResidual(
00376         const Vec currentSolution,
00377         Vec residual,
00378         DistributedVectorFactory& rFactory)
00379 {
00380     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00381     {
00382         this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
00383 
00384         DistributedVector solution_distributed = rFactory.CreateDistributedVector(currentSolution);
00385         DistributedVector residual_distributed = rFactory.CreateDistributedVector(residual);
00386 
00387 
00388         while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
00389         {
00390             DistributedVector::Stripe solution_stripe(solution_distributed, index_of_unknown);
00391             DistributedVector::Stripe residual_stripe(residual_distributed, index_of_unknown);
00392 
00393             unsigned node_index = this->mDirichIterator->first->GetIndex();
00394 
00395             double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
00396 
00397             if (solution_distributed.IsGlobalIndexLocal(node_index))
00398             {
00399                 residual_stripe[node_index]=solution_stripe[node_index] - value;
00400             }
00401             this->mDirichIterator++;
00402         }
00403         solution_distributed.Restore();
00404         residual_distributed.Restore();
00405     }
00406 }
00407 
00408 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00409 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ApplyDirichletToNonlinearJacobian(Mat jacobian)
00410 {
00411     unsigned num_boundary_conditions = 0;
00412     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00413     {
00414         num_boundary_conditions += this->mpDirichletMap[index_of_unknown]->size();
00415     }
00416 
00417     std::vector<unsigned> rows_to_zero(num_boundary_conditions);
00418 
00419     unsigned counter=0;
00420     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00421     {
00422         this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
00423 
00424         while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
00425         {
00426             unsigned node_index = this->mDirichIterator->first->GetIndex();
00427             rows_to_zero[counter++] = PROBLEM_DIM*node_index + index_of_unknown;
00428             this->mDirichIterator++;
00429         }
00430     }
00431     PetscMatTools::AssembleFinal(jacobian);
00432     PetscMatTools::ZeroRowsWithValueOnDiagonal(jacobian, rows_to_zero, 1.0);
00433 }
00434 
00435 
00436 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00437 bool BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::Validate(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00438 {
00439     bool valid = true;
00440 
00441     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00442     {
00443         // Iterate over surface elements
00444         typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::BoundaryElementIterator elt_iter
00445         = pMesh->GetBoundaryElementIteratorBegin();
00446         while (valid && elt_iter != pMesh->GetBoundaryElementIteratorEnd())
00447         {
00448             if (!HasNeumannBoundaryCondition(*elt_iter, index_of_unknown))
00449             {
00450                 // Check for Dirichlet conditions on this element's nodes
00451                 for (unsigned i=0; i<(*elt_iter)->GetNumNodes(); i++)
00452                 {
00453                     if (!HasDirichletBoundaryCondition((*elt_iter)->GetNode(i)))
00454                     {
00455                         valid = false;
00456                     }
00457                 }
00458             }
00459             elt_iter++;
00460         }
00461     }
00462     return valid;
00463 }
00464 
00465 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00466 double BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::GetNeumannBCValue(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00467                              const ChastePoint<SPACE_DIM>& rX,
00468                              unsigned indexOfUnknown)
00469 {
00470     assert(indexOfUnknown < PROBLEM_DIM);
00471 
00472     // Did we see this condition on the last search we did?
00473     if (mLastNeumannCondition[indexOfUnknown] == mpNeumannMap[indexOfUnknown]->end() ||
00474         mLastNeumannCondition[indexOfUnknown]->first != pSurfaceElement)
00475     {
00476         mLastNeumannCondition[indexOfUnknown] = mpNeumannMap[indexOfUnknown]->find(pSurfaceElement);
00477     }
00478     if (mLastNeumannCondition[indexOfUnknown] == mpNeumannMap[indexOfUnknown]->end())
00479     {
00480         // No Neumann condition is equivalent to a zero Neumann condition
00481         return 0.0;
00482     }
00483     else
00484     {
00485         return mLastNeumannCondition[indexOfUnknown]->second->GetValue(rX);
00486     }
00487 }
00488 
00489 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00490 bool BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::HasNeumannBoundaryCondition(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00491                                      unsigned indexOfUnknown)
00492 {
00493     assert(indexOfUnknown < PROBLEM_DIM);
00494 
00495     mLastNeumannCondition[indexOfUnknown] = mpNeumannMap[indexOfUnknown]->find(pSurfaceElement);
00496 
00497     return (mLastNeumannCondition[indexOfUnknown] != mpNeumannMap[indexOfUnknown]->end());
00498 }
00499 
00500 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00501 bool BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::AnyNonZeroNeumannConditions()
00502 {
00503     bool ret = false;
00504     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00505     {
00506         if (mAnyNonZeroNeumannConditionsForUnknown[index_of_unknown] == true)
00507         {
00508             ret = true;
00509         }
00510     }
00511     return ret;
00512 }
00513 
00514 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00515 typename BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::NeumannMapIterator BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::BeginNeumann()
00516 {
00517     // [0] is ok as all maps will be in sync due to the way ApplyNeumannBoundaryCondition works
00518     return mpNeumannMap[0]->begin();
00519 }
00520 
00521 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00522 typename BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::NeumannMapIterator BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::EndNeumann()
00523 {
00524     // [0] is ok as all maps will be in sync due to the way ApplyNeumannBoundaryCondition works
00525     return mpNeumannMap[0]->end();
00526 }
00527 
00528 #endif // _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_

Generated on Tue May 31 14:31:48 2011 for Chaste by  doxygen 1.5.5