BoundaryConditionsContainerImplementation.hpp

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

Generated by  doxygen 1.6.2