Chaste  Release::2018.1
BoundaryConditionsContainer.hpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #ifndef _BOUNDARYCONDITIONSCONTAINER_HPP_
37 #define _BOUNDARYCONDITIONSCONTAINER_HPP_
38 
39 #include <map>
40 #include "ChasteSerialization.hpp"
41 #include <boost/serialization/split_member.hpp>
42 #include <boost/serialization/map.hpp>
43 
44 #include "AbstractBoundaryConditionsContainer.hpp"
45 #include "AbstractBoundaryCondition.hpp"
46 #include "AbstractTetrahedralMesh.hpp"
47 #include "BoundaryElement.hpp"
48 #include "Node.hpp"
49 #include "LinearSystem.hpp"
50 #include "PetscException.hpp"
51 #include "ChastePoint.hpp"
52 #include "ConstBoundaryCondition.hpp"
53 #include "DistributedVectorFactory.hpp"
54 
66 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
67 class BoundaryConditionsContainer : public AbstractBoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>
68 {
69 public:
70 
72  typedef typename std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>*, const AbstractBoundaryCondition<SPACE_DIM>* >::const_iterator
74 
77 
78 private:
79 
80  std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *, const AbstractBoundaryCondition<SPACE_DIM>* >*
81  mpNeumannMap[PROBLEM_DIM];
84  std::map< const Node<SPACE_DIM> *, const Node<SPACE_DIM> * >* mpPeriodicBcMap[PROBLEM_DIM];
85 
90 
95 
98 
101 
102 public:
103 
110  BoundaryConditionsContainer(bool deleteConditions=true);
111 
117 
130  void AddDirichletBoundaryCondition(const Node<SPACE_DIM>* pBoundaryNode,
131  const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
132  unsigned indexOfUnknown = 0,
133  bool checkIfBoundaryNode = true);
134 
154  const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
155  unsigned indexOfUnknown = 0);
156 
157 
166  const Node<SPACE_DIM>* pNode2);
167 
168 
177  unsigned indexOfUnknown = 0);
178 
188  double value,
189  unsigned indexOfUnknown = 0);
190 
199  unsigned indexOfUnknown = 0);
200 
215  void ApplyDirichletToLinearProblem(LinearSystem& rLinearSystem,
216  bool applyToMatrix = true,
217  bool applyToRhsVector = true);
218 
238  void ApplyPeriodicBcsToLinearProblem(LinearSystem& rLinearSystem,
239  bool applyToMatrix = true,
240  bool applyToRhsVector = true);
241 
254  void ApplyDirichletToNonlinearResidual(const Vec currentSolution, Vec residual,
255  DistributedVectorFactory& rFactory);
256 
267  void ApplyDirichletToNonlinearJacobian(Mat jacobian);
268 
280 
290  double GetNeumannBCValue(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>* pSurfaceElement,
291  const ChastePoint<SPACE_DIM>& rX,
292  unsigned indexOfUnknown = 0);
293 
304  unsigned indexOfUnknown = 0);
305 
310 
315 
320 
334  template <class Archive>
336  {
337  if (mLoadedFromArchive)
338  {
339  return;
340  }
341 
342  MergeFromArchive(archive, pMesh);
343  }
344 
354  template <class Archive>
356 
357 private:
358 
360  friend class boost::serialization::access;
361 
368  template<class Archive>
369  void save(Archive & archive, const unsigned int version) const;
370 
388  template<class Archive>
389  void load(Archive & archive, const unsigned int version)
390  {
391  }
392 
393  BOOST_SERIALIZATION_SPLIT_MEMBER()
394 };
395 
396 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
397 template<class Archive>
399  Archive & archive, const unsigned int version) const
400 {
401  typedef typename std::map<unsigned, const AbstractBoundaryCondition<SPACE_DIM> *> archive_map_type;
402 
403  // Save Dirichlet conditions
404  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
405  {
406  archive_map_type bc_map;
407  typename BaseClassType::DirichletIteratorType it = this->mpDirichletMap[index_of_unknown]->begin();
408  while (it != this->mpDirichletMap[index_of_unknown]->end() )
409  {
410  unsigned node_index = it->first->GetIndex();
411  const AbstractBoundaryCondition<SPACE_DIM> * p_cond = it->second;
412  bc_map[node_index] = p_cond;
413 
414  it++;
415  }
416  archive & bc_map;
417  }
418 
419  // Save Neumann conditions
420  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
421  {
422  archive_map_type bc_map;
423  for (NeumannMapIterator it = mpNeumannMap[index_of_unknown]->begin();
424  it != mpNeumannMap[index_of_unknown]->end();
425  ++it)
426  {
427  unsigned elem_index = it->first->GetIndex();
428  const AbstractBoundaryCondition<SPACE_DIM>* p_cond = it->second;
429  bc_map[elem_index] = p_cond;
430  }
431  archive & bc_map;
432  }
433 }
434 
435 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
436 template<class Archive>
438  Archive & archive, AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
439 {
440  mLoadedFromArchive = true;
441 
442  typedef typename std::map<unsigned, AbstractBoundaryCondition<SPACE_DIM>*> archive_map_type;
443 
444  // Keep track of conditions that might need deleting
445  std::set<const AbstractBoundaryCondition<SPACE_DIM>*> maybe_unused_bcs;
446  std::set<const AbstractBoundaryCondition<SPACE_DIM>*> used_bcs;
447 
448  // Load Dirichlet conditions
449  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
450  {
451  archive_map_type bc_map;
452  archive & bc_map;
453  for (typename archive_map_type::iterator it = bc_map.begin();
454  it != bc_map.end();
455  ++it)
456  {
457  unsigned node_index = it->first;
458  this->mHasDirichletBCs=true; //We know that a Dirichlet is being added, even if not by this process
459  Node<SPACE_DIM>* p_node;
460  try
461  {
462  p_node = pMesh->GetNodeFromPrePermutationIndex(node_index);
463  }
464  catch (Exception&)
465  {
466  // It's a distributed mesh and we don't own this node - skip to the next BC
467  maybe_unused_bcs.insert(it->second);
468  continue;
469  }
470  AddDirichletBoundaryCondition(p_node, it->second, index_of_unknown, false);
471  used_bcs.insert(it->second);
472  }
473  }
474  this->mCheckedAndCommunicatedIfDirichletBcs=true; // Whether the Dirichlet BCC was empty or not, all processes know the status.
475 
476  // Load Neumann conditions
477  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
478  {
479  archive_map_type bc_map;
480  archive & bc_map;
481  for (typename archive_map_type::iterator it = bc_map.begin();
482  it != bc_map.end();
483  ++it)
484  {
485  unsigned boundary_element_index = it->first;
486  BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element;
487  try
488  {
489  p_boundary_element = pMesh->GetBoundaryElement(boundary_element_index);
490  }
491  catch (Exception&)
492  {
493  // It's a distributed mesh and we don't own this element - skip to the next BC
494  maybe_unused_bcs.insert(it->second);
495  continue;
496  }
497  AddNeumannBoundaryCondition(p_boundary_element, it->second, index_of_unknown);
498  used_bcs.insert(it->second);
499  }
500  }
501 
502  // Free any unused BCs
503  for (typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator it=maybe_unused_bcs.begin();
504  it != maybe_unused_bcs.end();
505  ++it)
506  {
507  typename std::set<const AbstractBoundaryCondition<SPACE_DIM>*>::iterator used = used_bcs.find(*it);
508  if (used == used_bcs.end())
509  {
510  delete (*it);
511  }
512  }
513 }
514 
515 #endif //_BOUNDARYCONDITIONSCONTAINER_HPP_
void load(Archive &archive, const unsigned int version)
void AddPeriodicBoundaryCondition(const Node< SPACE_DIM > *pNode1, const Node< SPACE_DIM > *pNode2)
Definition: Node.hpp:58
void ApplyDirichletToNonlinearResidual(const Vec currentSolution, Vec residual, DistributedVectorFactory &rFactory)
NeumannMapIterator mLastNeumannCondition[PROBLEM_DIM]
void DefineZeroDirichletOnMeshBoundary(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, unsigned indexOfUnknown=0)
void AddNeumannBoundaryCondition(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *pBoundaryElement, const AbstractBoundaryCondition< SPACE_DIM > *pBoundaryCondition, unsigned indexOfUnknown=0)
void LoadFromArchive(Archive &archive, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
std::map< const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *, const AbstractBoundaryCondition< SPACE_DIM > * >::const_iterator NeumannMapIterator
void ApplyDirichletToLinearProblem(LinearSystem &rLinearSystem, bool applyToMatrix=true, bool applyToRhsVector=true)
double GetNeumannBCValue(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *pSurfaceElement, const ChastePoint< SPACE_DIM > &rX, unsigned indexOfUnknown=0)
void ApplyPeriodicBcsToLinearProblem(LinearSystem &rLinearSystem, bool applyToMatrix=true, bool applyToRhsVector=true)
void AddDirichletBoundaryCondition(const Node< SPACE_DIM > *pBoundaryNode, const AbstractBoundaryCondition< SPACE_DIM > *pBoundaryCondition, unsigned indexOfUnknown=0, bool checkIfBoundaryNode=true)
std::map< const Node< SPACE_DIM > *, const Node< SPACE_DIM > * > * mpPeriodicBcMap[PROBLEM_DIM]
void MergeFromArchive(Archive &archive, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
bool Validate(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
ConstBoundaryCondition< SPACE_DIM > * mpZeroBoundaryCondition
BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * GetBoundaryElement(unsigned index) const
void save(Archive &archive, const unsigned int version) const
std::map< const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *, const AbstractBoundaryCondition< SPACE_DIM > * > * mpNeumannMap[PROBLEM_DIM]
void DefineConstantDirichletOnMeshBoundary(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, double value, unsigned indexOfUnknown=0)
std::map< const Node< SPACE_DIM > *, const AbstractBoundaryCondition< SPACE_DIM > *, LessThanNode< SPACE_DIM > >::const_iterator DirichletIteratorType
AbstractBoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > BaseClassType
void DefineZeroNeumannOnMeshBoundary(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, unsigned indexOfUnknown=0)
bool HasNeumannBoundaryCondition(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *pSurfaceElement, unsigned indexOfUnknown=0)
bool mAnyNonZeroNeumannConditionsForUnknown[PROBLEM_DIM]
Node< SPACE_DIM > * GetNodeFromPrePermutationIndex(unsigned index) const