BoundaryConditionsContainer.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 #ifndef _BOUNDARYCONDITIONSCONTAINER_HPP_
00029 #define _BOUNDARYCONDITIONSCONTAINER_HPP_
00030 
00031 #include <map>
00032 #include "ChasteSerialization.hpp"
00033 #include <boost/serialization/split_member.hpp>
00034 #include <boost/serialization/map.hpp>
00035 
00036 #include "AbstractBoundaryConditionsContainer.hpp"
00037 #include "AbstractBoundaryCondition.hpp"
00038 #include "AbstractTetrahedralMesh.hpp"
00039 #include "BoundaryElement.hpp"
00040 #include "Node.hpp"
00041 #include "LinearSystem.hpp"
00042 #include "PetscException.hpp"
00043 #include "ChastePoint.hpp"
00044 #include "ConstBoundaryCondition.hpp"
00045 #include "DistributedVectorFactory.hpp"
00046 
00047 
00059 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00060 class BoundaryConditionsContainer : public AbstractBoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>
00061 {
00062 public:
00063 
00065     typedef typename std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>*, const AbstractBoundaryCondition<SPACE_DIM>* >::const_iterator
00066         NeumannMapIterator;
00067 
00069     typedef AbstractBoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM> BaseClassType;
00070 
00071 private:
00072 
00073     std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *, const AbstractBoundaryCondition<SPACE_DIM>* >*
00074         mpNeumannMap[PROBLEM_DIM]; 
00079     NeumannMapIterator mLastNeumannCondition[PROBLEM_DIM];
00080 
00084     bool mAnyNonZeroNeumannConditionsForUnknown[PROBLEM_DIM];
00085 
00087     ConstBoundaryCondition<SPACE_DIM>* mpZeroBoundaryCondition;
00088 
00090     bool mLoadedFromArchive;
00091 
00092 public:
00093 
00100     BoundaryConditionsContainer(bool deleteConditions=true);
00101 
00106     ~BoundaryConditionsContainer();
00107 
00120     void AddDirichletBoundaryCondition(const Node<SPACE_DIM>* pBoundaryNode,
00121                                        const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
00122                                        unsigned indexOfUnknown = 0,
00123                                        bool checkIfBoundaryNode = true);
00124 
00143     void AddNeumannBoundaryCondition(const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* pBoundaryElement,
00144                                      const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
00145                                      unsigned indexOfUnknown = 0);
00146 
00154     void DefineZeroDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00155                                            unsigned indexOfUnknown = 0);
00156 
00165     void DefineConstantDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00166                                                double value,
00167                                                unsigned indexOfUnknown = 0);
00168 
00176     void DefineZeroNeumannOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00177                                          unsigned indexOfUnknown = 0);
00178 
00193     void ApplyDirichletToLinearProblem(LinearSystem& rLinearSystem,
00194                                        bool applyToMatrix = true,
00195                                        bool applyToRhsVector = true);
00196 
00209     void ApplyDirichletToNonlinearResidual(const Vec currentSolution, Vec residual,
00210                                            DistributedVectorFactory& rFactory);
00211 
00222     void ApplyDirichletToNonlinearJacobian(Mat jacobian);
00223 
00234     bool Validate(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00235 
00245     double GetNeumannBCValue(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00246                              const ChastePoint<SPACE_DIM>& rX,
00247                              unsigned indexOfUnknown = 0);
00248 
00258     bool HasNeumannBoundaryCondition(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00259                                      unsigned indexOfUnknown = 0);
00260 
00264     bool AnyNonZeroNeumannConditions();
00265 
00269     NeumannMapIterator BeginNeumann();
00270 
00274     NeumannMapIterator EndNeumann();
00275 
00289     template <class Archive>
00290     void LoadFromArchive(Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00291     {
00292         if (mLoadedFromArchive)
00293         {
00294             return;
00295         }
00296 
00297         MergeFromArchive(archive, pMesh);
00298     }
00299 
00309     template <class Archive>
00310     void MergeFromArchive(Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00311 
00312 private:
00314     friend class boost::serialization::access;
00315 
00322     template<class Archive>
00323     void save(Archive & archive, const unsigned int version) const;
00324 
00342     template<class Archive>
00343     void load(Archive & archive, const unsigned int version)
00344     {
00345     }
00346 
00347     BOOST_SERIALIZATION_SPLIT_MEMBER()
00348 
00349 };
00350 
00351 
00352 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00353 template<class Archive>
00354 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::save(
00355         Archive & archive, const unsigned int version) const
00356 {
00357     typedef typename std::map<unsigned, const AbstractBoundaryCondition<SPACE_DIM> *> archive_map_type;
00358     // Save Dirichlet conditions
00359     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00360     {
00361         archive_map_type bc_map;
00362         typename BaseClassType::DirichletIteratorType it = this->mpDirichletMap[index_of_unknown]->begin();
00363         while (it != this->mpDirichletMap[index_of_unknown]->end() )
00364         {
00365             unsigned node_index = it->first->GetIndex();
00366             const AbstractBoundaryCondition<SPACE_DIM> * p_cond = it->second;
00367             bc_map[node_index] = p_cond;
00368 
00369             it++;
00370         }
00371         archive & bc_map;
00372     }
00373 
00374     // Save Neumann conditions
00375     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00376     {
00377         archive_map_type bc_map;
00378         for (NeumannMapIterator it = mpNeumannMap[index_of_unknown]->begin();
00379              it != mpNeumannMap[index_of_unknown]->end();
00380              ++it)
00381         {
00382             unsigned elem_index = it->first->GetIndex();
00383             const AbstractBoundaryCondition<SPACE_DIM>* p_cond = it->second;
00384             bc_map[elem_index] = p_cond;
00385         }
00386         archive & bc_map;
00387     }
00388 }
00389 
00390 
00391 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00392 template<class Archive>
00393 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::MergeFromArchive(
00394         Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00395 {
00396     mLoadedFromArchive = true;
00397 
00398     typedef typename std::map<unsigned, AbstractBoundaryCondition<SPACE_DIM>*> archive_map_type;
00399     // Keep track of conditions that might need deleting
00400     std::set<const AbstractBoundaryCondition<SPACE_DIM>*> maybe_unused_bcs;
00401     std::set<const AbstractBoundaryCondition<SPACE_DIM>*> used_bcs;
00402     
00403     // Load Dirichlet conditions
00404     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00405     {
00406         archive_map_type bc_map;
00407         archive & bc_map;
00408         for (typename archive_map_type::iterator it = bc_map.begin();
00409              it != bc_map.end();
00410              ++it)
00411         {
00412             unsigned node_index = it->first;
00413             this->mHasDirichletBCs=true;  //We know that a Dirichlet is being added, even if not by this process
00414             Node<SPACE_DIM>* p_node;
00415             try
00416             {
00417                 p_node = pMesh->GetNodeFromPrePermutationIndex(node_index);
00418             }
00419             catch (Exception& e)
00420             {
00421                 // It's a distributed mesh and we don't own this node - skip to the next BC
00422                 maybe_unused_bcs.insert(it->second);
00423                 continue;
00424             }
00425             AddDirichletBoundaryCondition(p_node, it->second, index_of_unknown, false);
00426             used_bcs.insert(it->second);
00427         }
00428     }
00429     this->mCheckedAndCommunicatedIfDirichletBcs=true; // Whether the Dirichlet BCC was empty or not, all processes know the status.
00430 
00431     // Load Neumann conditions
00432     for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00433     {
00434         archive_map_type bc_map;
00435         archive & bc_map;
00436         for (typename archive_map_type::iterator it = bc_map.begin();
00437              it != bc_map.end();
00438              ++it)
00439         {
00440             unsigned boundary_element_index = it->first;
00441             BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element;
00442             try
00443             {
00444                 p_boundary_element = pMesh->GetBoundaryElement(boundary_element_index);
00445             }
00446             catch (Exception& e)
00447             {
00448                 // It's a distributed mesh and we don't own this element - skip to the next BC
00449                 maybe_unused_bcs.insert(it->second);
00450                 continue;
00451             }
00452             AddNeumannBoundaryCondition(p_boundary_element, it->second, index_of_unknown);
00453             used_bcs.insert(it->second);
00454         }
00455     }
00456     
00457     // Free any unused BCs
00458     for (typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator it=maybe_unused_bcs.begin();
00459          it != maybe_unused_bcs.end();
00460          ++it)
00461     {
00462         typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator used = used_bcs.find(*it);
00463         if (used == used_bcs.end())
00464         {
00465             delete (*it);
00466         }
00467     }
00468 }
00469 
00470 #endif //_BOUNDARYCONDITIONSCONTAINER_HPP_

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