00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00098 BoundaryConditionsContainer();
00099
00104 ~BoundaryConditionsContainer();
00105
00118 void AddDirichletBoundaryCondition(const Node<SPACE_DIM>* pBoundaryNode,
00119 const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
00120 unsigned indexOfUnknown = 0,
00121 bool checkIfBoundaryNode = true);
00122
00141 void AddNeumannBoundaryCondition(const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* pBoundaryElement,
00142 const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
00143 unsigned indexOfUnknown = 0);
00144
00152 void DefineZeroDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00153 unsigned indexOfUnknown = 0);
00154
00163 void DefineConstantDirichletOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00164 double value,
00165 unsigned indexOfUnknown = 0);
00166
00174 void DefineZeroNeumannOnMeshBoundary(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00175 unsigned indexOfUnknown = 0);
00176
00191 void ApplyDirichletToLinearProblem(LinearSystem& rLinearSystem,
00192 bool applyToMatrix = true,
00193 bool applyToRhsVector = true);
00194
00207 void ApplyDirichletToNonlinearResidual(const Vec currentSolution, Vec residual,
00208 DistributedVectorFactory& rFactory);
00209
00220 void ApplyDirichletToNonlinearJacobian(Mat jacobian);
00221
00232 bool Validate(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00233
00243 double GetNeumannBCValue(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00244 const ChastePoint<SPACE_DIM>& rX,
00245 unsigned indexOfUnknown = 0);
00246
00257 bool HasNeumannBoundaryCondition(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
00258 unsigned indexOfUnknown = 0);
00259
00263 bool AnyNonZeroNeumannConditions();
00264
00268 NeumannMapIterator BeginNeumann();
00269
00273 NeumannMapIterator EndNeumann();
00274
00288 template <class Archive>
00289 void LoadFromArchive(Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00290 {
00291 if (mLoadedFromArchive)
00292 {
00293 return;
00294 }
00295
00296 MergeFromArchive(archive, pMesh);
00297 }
00298
00308 template <class Archive>
00309 void MergeFromArchive(Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00310
00311 private:
00313 friend class boost::serialization::access;
00314
00321 template<class Archive>
00322 void save(Archive & archive, const unsigned int version) const;
00323
00341 template<class Archive>
00342 void load(Archive & archive, const unsigned int version)
00343 {
00344 }
00345
00346 BOOST_SERIALIZATION_SPLIT_MEMBER()
00347
00348 };
00349
00350
00351 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00352 template<class Archive>
00353 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::save(
00354 Archive & archive, const unsigned int version) const
00355 {
00356 typedef typename std::map<unsigned, const AbstractBoundaryCondition<SPACE_DIM> *> archive_map_type;
00357
00358 for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00359 {
00360 archive_map_type bc_map;
00361 typename BaseClassType::DirichletIteratorType it = this->mpDirichletMap[index_of_unknown]->begin();
00362 while (it != this->mpDirichletMap[index_of_unknown]->end() )
00363 {
00364 unsigned node_index = it->first->GetIndex();
00365 const AbstractBoundaryCondition<SPACE_DIM> * p_cond = it->second;
00366 bc_map[node_index] = p_cond;
00367
00368 it++;
00369 }
00370 archive & bc_map;
00371 }
00372
00373
00374 for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00375 {
00376 archive_map_type bc_map;
00377 for (NeumannMapIterator it = mpNeumannMap[index_of_unknown]->begin();
00378 it != mpNeumannMap[index_of_unknown]->end();
00379 ++it)
00380 {
00381 unsigned elem_index = it->first->GetIndex();
00382 const AbstractBoundaryCondition<SPACE_DIM>* p_cond = it->second;
00383 bc_map[elem_index] = p_cond;
00384 }
00385 archive & bc_map;
00386 }
00387 }
00388
00389
00390 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00391 template<class Archive>
00392 void BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::MergeFromArchive(
00393 Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00394 {
00395 mLoadedFromArchive = true;
00396
00397 typedef typename std::map<unsigned, AbstractBoundaryCondition<SPACE_DIM>*> archive_map_type;
00398
00399 std::set<const AbstractBoundaryCondition<SPACE_DIM>*> maybe_unused_bcs;
00400 std::set<const AbstractBoundaryCondition<SPACE_DIM>*> used_bcs;
00401
00402
00403 for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00404 {
00405 archive_map_type bc_map;
00406 archive & bc_map;
00407 for (typename archive_map_type::iterator it = bc_map.begin();
00408 it != bc_map.end();
00409 ++it)
00410 {
00411 unsigned node_index = it->first;
00412 this->mHasDirichletBCs=true;
00413 Node<SPACE_DIM>* p_node;
00414 try
00415 {
00416 p_node = pMesh->GetNodeFromPrePermutationIndex(node_index);
00417 }
00418 catch (Exception& e)
00419 {
00420
00421 maybe_unused_bcs.insert(it->second);
00422 continue;
00423 }
00424 AddDirichletBoundaryCondition(p_node, it->second, index_of_unknown, false);
00425 used_bcs.insert(it->second);
00426 }
00427 }
00428 this->mCheckedAndCommunicatedIfDirichletBcs=true;
00429
00430
00431 for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00432 {
00433 archive_map_type bc_map;
00434 archive & bc_map;
00435 for (typename archive_map_type::iterator it = bc_map.begin();
00436 it != bc_map.end();
00437 ++it)
00438 {
00439 unsigned boundary_element_index = it->first;
00440 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element;
00441 try
00442 {
00443 p_boundary_element = pMesh->GetBoundaryElement(boundary_element_index);
00444 }
00445 catch (Exception& e)
00446 {
00447
00448 maybe_unused_bcs.insert(it->second);
00449 continue;
00450 }
00451 AddNeumannBoundaryCondition(p_boundary_element, it->second, index_of_unknown);
00452 used_bcs.insert(it->second);
00453 }
00454 }
00455
00456
00457 for (typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator it=maybe_unused_bcs.begin();
00458 it != maybe_unused_bcs.end();
00459 ++it)
00460 {
00461 typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator used = used_bcs.find(*it);
00462 if (used == used_bcs.end())
00463 {
00464 delete (*it);
00465 }
00466 }
00467 }
00468
00469 #endif //_BOUNDARYCONDITIONSCONTAINER_HPP_