Chaste  Release::2018.1
BoundaryConditionsContainerImplementation.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 _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_
37 #define _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_
38 
39 #include "BoundaryConditionsContainer.hpp"
40 #include "ConstBoundaryCondition.hpp"
41 #include "DistributedVector.hpp"
42 #include "Exception.hpp"
43 #include "HeartEventHandler.hpp"
44 #include "PetscMatTools.hpp"
45 #include "PetscVecTools.hpp"
46 #include "ReplicatableVector.hpp"
47 
48 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
50  : AbstractBoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>(deleteConditions)
51 {
52  mLoadedFromArchive = false;
53 
54  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
55  {
56  mpNeumannMap[index_of_unknown] = new std::map< const BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *, const AbstractBoundaryCondition<SPACE_DIM>*>;
57 
58  mAnyNonZeroNeumannConditionsForUnknown[index_of_unknown] = false;
59  mLastNeumannCondition[index_of_unknown] = mpNeumannMap[index_of_unknown]->begin();
60 
61  mpPeriodicBcMap[index_of_unknown] = new std::map< const Node<SPACE_DIM> *, const Node<SPACE_DIM> * >;
62  }
63 
64  // This zero boundary condition is only used in AddNeumannBoundaryCondition
66 }
67 
68 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
70 {
71  // Keep track of what boundary condition objects we've deleted
72  std::set<const AbstractBoundaryCondition<SPACE_DIM>*> deleted_conditions;
73  for (unsigned i=0; i<PROBLEM_DIM; i++)
74  {
75  NeumannMapIterator neumann_iterator = mpNeumannMap[i]->begin();
76  while (neumann_iterator != mpNeumannMap[i]->end() )
77  {
78 
79  if (deleted_conditions.count(neumann_iterator->second) == 0)
80  {
81  deleted_conditions.insert(neumann_iterator->second);
82  //Leave the zero boundary condition until last
83  if (neumann_iterator->second != mpZeroBoundaryCondition)
84  {
85  if (this->mDeleteConditions)
86  {
87  delete neumann_iterator->second;
88  }
89  }
90  }
91  neumann_iterator++;
92  }
93  delete(mpNeumannMap[i]);
94  delete(mpPeriodicBcMap[i]);
95  }
96 
97  delete mpZeroBoundaryCondition;
98 
99  if (this->mDeleteConditions)
100  {
101  this->DeleteDirichletBoundaryConditions(deleted_conditions);
102  }
103 }
104 
105 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
107  const AbstractBoundaryCondition<SPACE_DIM>* pBoundaryCondition,
108  unsigned indexOfUnknown,
109  bool checkIfBoundaryNode)
110 {
111  assert(indexOfUnknown < PROBLEM_DIM);
112  if (checkIfBoundaryNode)
113  {
114  assert(pBoundaryNode->IsBoundaryNode());
115  }
116 
117  (*(this->mpDirichletMap[indexOfUnknown]))[pBoundaryNode] = pBoundaryCondition;
118 }
119 
120 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
122  const Node<SPACE_DIM>* pNode2)
123 {
124  assert(pNode1->IsBoundaryNode());
125  assert(pNode2->IsBoundaryNode());
126 
127  // will assume the periodic BC is to be applied to ALL unknowns, can't really imagine a
128  // situation where this isn't going to be true. If necessary can easily change this method
129  // to take in the index of the unknown
130  for (unsigned i=0; i<PROBLEM_DIM; i++)
131  {
132  (*(this->mpPeriodicBcMap[i]))[pNode1] = pNode2;
133  }
134 }
135 
136 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
138  const AbstractBoundaryCondition<SPACE_DIM> * pBoundaryCondition,
139  unsigned indexOfUnknown)
140 {
141  assert(indexOfUnknown < PROBLEM_DIM);
142 
143  /*
144  * If this condition is constant, we can test whether it is zero.
145  * Otherwise we assume that this could be a non-zero boundary condition.
146  */
147  const ConstBoundaryCondition<SPACE_DIM>* p_const_cond = dynamic_cast<const ConstBoundaryCondition<SPACE_DIM>*>(pBoundaryCondition);
148  if (p_const_cond)
149  {
150  if (p_const_cond->GetValue(pBoundaryElement->GetNode(0)->GetPoint()) != 0.0)
151  {
152  mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = true;
153  }
154  }
155  else
156  {
157  mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = true;
158  }
159 
160  for (unsigned unknown=0; unknown<PROBLEM_DIM; unknown++)
161  {
162  if (unknown == indexOfUnknown)
163  {
164  (*(mpNeumannMap[indexOfUnknown]))[pBoundaryElement] = pBoundaryCondition;
165  }
166  else
167  {
168  // If can't find pBoundaryElement in map[unknown]
169  if (mpNeumannMap[unknown]->find(pBoundaryElement)==mpNeumannMap[unknown]->end())
170  {
171  // Add zero bc to other unknowns (so all maps are in sync)
172  (*(mpNeumannMap[unknown]))[pBoundaryElement] = mpZeroBoundaryCondition;
173  }
174  }
175  }
176 }
177 
178 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
180  unsigned indexOfUnknown)
181 {
182  this->DefineConstantDirichletOnMeshBoundary(pMesh, 0.0, indexOfUnknown);
183 }
184 
185 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
187  double value,
188  unsigned indexOfUnknown)
189 {
190  assert(indexOfUnknown < PROBLEM_DIM);
191  // In applying a condition to the boundary, we need to be sure that the boundary exists
192  assert(PetscTools::ReplicateBool( pMesh->GetNumBoundaryNodes() > 0 ) );
193 
194  ConstBoundaryCondition<SPACE_DIM>* p_boundary_condition = new ConstBoundaryCondition<SPACE_DIM>(value);
195 
197  iter = pMesh->GetBoundaryNodeIteratorBegin();
198  while (iter != pMesh->GetBoundaryNodeIteratorEnd())
199  {
200  AddDirichletBoundaryCondition(*iter, p_boundary_condition, indexOfUnknown);
201  iter++;
202  }
203 }
204 
205 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
207  unsigned indexOfUnknown)
208 {
209  assert(indexOfUnknown < PROBLEM_DIM);
210 
211  // In applying a condition to the boundary, we need to be sure that the boundary exists
212  assert(pMesh->GetNumBoundaryElements() > 0);
213  ConstBoundaryCondition<SPACE_DIM>* p_zero_boundary_condition = new ConstBoundaryCondition<SPACE_DIM>( 0.0 );
214 
216  iter = pMesh->GetBoundaryElementIteratorBegin();
217  while (iter != pMesh->GetBoundaryElementIteratorEnd())
218  {
219  AddNeumannBoundaryCondition(*iter, p_zero_boundary_condition, indexOfUnknown);
220  iter++;
221  }
222 
223  mAnyNonZeroNeumannConditionsForUnknown[indexOfUnknown] = false;
224 }
225 
256 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
258  LinearSystem& rLinearSystem,
259  bool applyToMatrix,
260  bool applyToRhsVector)
261 {
262  HeartEventHandler::BeginEvent(HeartEventHandler::DIRICHLET_BCS);
263 
264  if (applyToMatrix)
265  {
266  if (!this->HasDirichletBoundaryConditions())
267  {
268  // Short-circuit the replication if there are no conditions
269  HeartEventHandler::EndEvent(HeartEventHandler::DIRICHLET_BCS);
270  return;
271  }
272 
273  bool matrix_is_symmetric = rLinearSystem.IsMatrixSymmetric();
274 
275  if (matrix_is_symmetric)
276  {
277  /*
278  * Modifications to the RHS are stored in the Dirichlet boundary
279  * conditions vector. This is done so that they can be reapplied
280  * at each time step.
281  * Make a new vector to store the Dirichlet offsets in.
282  */
283  Vec& r_bcs_vec = rLinearSystem.rGetDirichletBoundaryConditionsVector();
284  if (!r_bcs_vec)
285  {
286  VecDuplicate(rLinearSystem.rGetRhsVector(), &r_bcs_vec);
287  }
288  PetscVecTools::Zero(r_bcs_vec);
289  /*
290  * If the matrix is symmetric, calls to GetMatrixRowDistributed()
291  * require the matrix to be in assembled state. Otherwise we can
292  * defer it.
293  */
294  rLinearSystem.AssembleFinalLinearSystem();
295  }
296 
297  // Work out where we're setting Dirichlet boundary conditions *everywhere*, not just those locally known
298  ReplicatableVector dirichlet_conditions(rLinearSystem.GetSize());
299  unsigned lo, hi;
300  {
301  PetscInt lo_s, hi_s;
302  rLinearSystem.GetOwnershipRange(lo_s, hi_s);
303  lo = lo_s; hi = hi_s;
304  }
305  // Initialise all local entries to DBL_MAX, i.e. don't know if there's a condition
306  for (unsigned i=lo; i<hi; i++)
307  {
308  dirichlet_conditions[i] = DBL_MAX;
309  }
310  // Now fill in the ones we know
311  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
312  {
313  this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
314 
315  while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
316  {
317  unsigned node_index = this->mDirichIterator->first->GetIndex();
318  double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
319  assert(value != DBL_MAX);
320 
321  unsigned row = PROBLEM_DIM*node_index + index_of_unknown;
322  dirichlet_conditions[row] = value;
323 
324  this->mDirichIterator++;
325  }
326  }
327 
328  // And replicate
329  dirichlet_conditions.Replicate(lo, hi);
330 
331  // Which rows have conditions?
332  std::vector<unsigned> rows_to_zero;
333  for (unsigned i=0; i<dirichlet_conditions.GetSize(); i++)
334  {
335  if (dirichlet_conditions[i] != DBL_MAX)
336  {
337  rows_to_zero.push_back(i);
338  }
339  }
340 
341  if (matrix_is_symmetric)
342  {
343  // Modify the matrix columns
344  for (unsigned i=0; i<rows_to_zero.size(); i++)
345  {
346  unsigned col = rows_to_zero[i];
347  double minus_value = -dirichlet_conditions[col];
348 
349  /*
350  * Get a vector which will store the column of the matrix (column d,
351  * where d is the index of the row (and column) to be altered for the
352  * boundary condition. Since the matrix is symmetric when get row
353  * number "col" and treat it as a column. PETSc uses compressed row
354  * format and therefore getting rows is far more efficient than getting
355  * columns.
356  */
357  Vec matrix_col = rLinearSystem.GetMatrixRowDistributed(col);
358 
359  // Zero the correct entry of the column
360  PetscVecTools::SetElement(matrix_col, col, 0.0);
361 
362  /*
363  * Set up the RHS Dirichlet boundary conditions vector.
364  * Assuming one boundary at the zeroth node (x_0 = value), this is equal to
365  * -value*[0 a_21 a_31 .. a_N1]
366  * and will be added to the RHS.
367  */
368  PetscVecTools::AddScaledVector(rLinearSystem.rGetDirichletBoundaryConditionsVector(), matrix_col, minus_value);
369  PetscTools::Destroy(matrix_col);
370  }
371  }
372 
373  /*
374  * Now zero the appropriate rows and columns of the matrix. If the matrix
375  * is symmetric we apply the boundary conditions in a way the symmetry isn't
376  * lost (rows and columns). If not only the row is zeroed.
377  */
378  if (matrix_is_symmetric)
379  {
380  rLinearSystem.ZeroMatrixRowsAndColumnsWithValueOnDiagonal(rows_to_zero, 1.0);
381  }
382  else
383  {
384  rLinearSystem.ZeroMatrixRowsWithValueOnDiagonal(rows_to_zero, 1.0);
385  }
386  }
387 
388  if (applyToRhsVector)
389  {
390  // Apply the RHS boundary conditions modification if required.
391  if (rLinearSystem.rGetDirichletBoundaryConditionsVector())
392  {
394  }
395 
396  /*
397  * Apply the actual boundary condition to the RHS, note this must be
398  * done after the modification to the RHS vector.
399  */
400  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
401  {
402  this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
403 
404  while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
405  {
406  unsigned node_index = this->mDirichIterator->first->GetIndex();
407  double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
408 
409  unsigned row = PROBLEM_DIM*node_index + index_of_unknown;
410 
411  rLinearSystem.SetRhsVectorElement(row, value);
412 
413  this->mDirichIterator++;
414  }
415  }
416  }
417 
418  HeartEventHandler::EndEvent(HeartEventHandler::DIRICHLET_BCS);
419 }
420 
421 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
423  bool applyToMatrix,
424  bool applyToRhsVector)
425 {
426  bool has_periodic_bcs = false;
427  for (unsigned i=0; i<PROBLEM_DIM; i++)
428  {
429  if (!mpPeriodicBcMap[i]->empty())
430  {
431  has_periodic_bcs = true;
432  break;
433  }
434  }
435 
436  EXCEPT_IF_NOT(has_periodic_bcs);
437 
438  if (applyToMatrix)
439  {
440  std::vector<unsigned> rows_to_zero;
441  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
442  {
443  for (typename std::map< const Node<SPACE_DIM> *, const Node<SPACE_DIM> * >::const_iterator iter = mpPeriodicBcMap[index_of_unknown]->begin();
444  iter != mpPeriodicBcMap[index_of_unknown]->end();
445  ++iter)
446  {
447  unsigned node_index_1 = iter->first->GetIndex();
448  unsigned row_index_1 = PROBLEM_DIM*node_index_1 + index_of_unknown;
449  rows_to_zero.push_back(row_index_1);
450  }
451  }
452 
453  rLinearSystem.ZeroMatrixRowsWithValueOnDiagonal(rows_to_zero, 1.0);
454 
455  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
456  {
457  for (typename std::map< const Node<SPACE_DIM> *, const Node<SPACE_DIM> * >::const_iterator iter = mpPeriodicBcMap[index_of_unknown]->begin();
458  iter != mpPeriodicBcMap[index_of_unknown]->end();
459  ++iter)
460  {
461  unsigned node_index_1 = iter->first->GetIndex();
462  unsigned node_index_2 = iter->second->GetIndex();
463 
464  unsigned mat_index1 = PROBLEM_DIM*node_index_1 + index_of_unknown;
465  unsigned mat_index2 = PROBLEM_DIM*node_index_2 + index_of_unknown;
466  PetscMatTools::SetElement(rLinearSystem.rGetLhsMatrix(), mat_index1, mat_index2, -1.0);
467  }
468  }
469  }
470 
471  if (applyToRhsVector)
472  {
473  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
474  {
475  for (typename std::map< const Node<SPACE_DIM> *, const Node<SPACE_DIM> * >::const_iterator iter = mpPeriodicBcMap[index_of_unknown]->begin();
476  iter != mpPeriodicBcMap[index_of_unknown]->end();
477  ++iter)
478  {
479  unsigned node_index = iter->first->GetIndex();
480  unsigned row_index = PROBLEM_DIM*node_index + index_of_unknown;
481  rLinearSystem.SetRhsVectorElement(row_index, 0.0);
482  }
483  }
484  }
485 }
486 
487 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
489  const Vec currentSolution,
490  Vec residual,
491  DistributedVectorFactory& rFactory)
492 {
493  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
494  {
495  this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
496 
497  DistributedVector solution_distributed = rFactory.CreateDistributedVector(currentSolution, true /*Read-only*/);
498  DistributedVector residual_distributed = rFactory.CreateDistributedVector(residual);
499 
500 
501  while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
502  {
503  DistributedVector::Stripe solution_stripe(solution_distributed, index_of_unknown);
504  DistributedVector::Stripe residual_stripe(residual_distributed, index_of_unknown);
505 
506  unsigned node_index = this->mDirichIterator->first->GetIndex();
507 
508  double value = this->mDirichIterator->second->GetValue(this->mDirichIterator->first->GetPoint());
509 
510  if (solution_distributed.IsGlobalIndexLocal(node_index))
511  {
512  residual_stripe[node_index]=solution_stripe[node_index] - value;
513  }
514  this->mDirichIterator++;
515  }
516  // Don't restore the read-only one: solution_distributed.Restore();
517  residual_distributed.Restore();
518  }
519 }
520 
521 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
523 {
524  unsigned num_boundary_conditions = 0;
525  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
526  {
527  num_boundary_conditions += this->mpDirichletMap[index_of_unknown]->size();
528  }
529 
530  std::vector<unsigned> rows_to_zero(num_boundary_conditions);
531 
532  unsigned counter=0;
533  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
534  {
535  this->mDirichIterator = this->mpDirichletMap[index_of_unknown]->begin();
536 
537  while (this->mDirichIterator != this->mpDirichletMap[index_of_unknown]->end() )
538  {
539  unsigned node_index = this->mDirichIterator->first->GetIndex();
540  rows_to_zero[counter++] = PROBLEM_DIM*node_index + index_of_unknown;
541  this->mDirichIterator++;
542  }
543  }
544  PetscMatTools::Finalise(jacobian);
545  PetscMatTools::ZeroRowsWithValueOnDiagonal(jacobian, rows_to_zero, 1.0);
546 }
547 
548 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
550 {
551  bool valid = true;
552 
553  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
554  {
555  // Iterate over surface elements
558  while (valid && elt_iter != pMesh->GetBoundaryElementIteratorEnd())
559  {
560  if (!HasNeumannBoundaryCondition(*elt_iter, index_of_unknown))
561  {
562  // Check for Dirichlet conditions on this element's nodes
563  for (unsigned i=0; i<(*elt_iter)->GetNumNodes(); i++)
564  {
565  if (!this->HasDirichletBoundaryCondition((*elt_iter)->GetNode(i)))
566  {
567  valid = false;
568  }
569  }
570  }
571  elt_iter++;
572  }
573  }
574  return valid;
575 }
576 
577 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
579  const ChastePoint<SPACE_DIM>& rX,
580  unsigned indexOfUnknown)
581 {
582  assert(indexOfUnknown < PROBLEM_DIM);
583 
584  // Did we see this condition on the last search we did?
585  if (mLastNeumannCondition[indexOfUnknown] == mpNeumannMap[indexOfUnknown]->end() ||
586  mLastNeumannCondition[indexOfUnknown]->first != pSurfaceElement)
587  {
588  mLastNeumannCondition[indexOfUnknown] = mpNeumannMap[indexOfUnknown]->find(pSurfaceElement);
589  }
590  if (mLastNeumannCondition[indexOfUnknown] == mpNeumannMap[indexOfUnknown]->end())
591  {
592  // No Neumann condition is equivalent to a zero Neumann condition
593  return 0.0;
594  }
595  else
596  {
597  return mLastNeumannCondition[indexOfUnknown]->second->GetValue(rX);
598  }
599 }
600 
601 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
603  unsigned indexOfUnknown)
604 {
605  assert(indexOfUnknown < PROBLEM_DIM);
606 
607  mLastNeumannCondition[indexOfUnknown] = mpNeumannMap[indexOfUnknown]->find(pSurfaceElement);
608 
609  return (mLastNeumannCondition[indexOfUnknown] != mpNeumannMap[indexOfUnknown]->end());
610 }
611 
612 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
614 {
615  bool ret = false;
616  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
617  {
618  if (mAnyNonZeroNeumannConditionsForUnknown[index_of_unknown] == true)
619  {
620  ret = true;
621  }
622  }
623  return ret;
624 }
625 
626 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
628 {
629  // [0] is ok as all maps will be in sync due to the way ApplyNeumannBoundaryCondition works
630  return mpNeumannMap[0]->begin();
631 }
632 
633 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
635 {
636  // [0] is ok as all maps will be in sync due to the way ApplyNeumannBoundaryCondition works
637  return mpNeumannMap[0]->end();
638 }
639 
640 #endif // _BOUNDARYCONDITIONSCONTAINERIMPLEMENTATION_HPP_
void AssembleFinalLinearSystem()
virtual unsigned GetNumBoundaryElements() const
void AddPeriodicBoundaryCondition(const Node< SPACE_DIM > *pNode1, const Node< SPACE_DIM > *pNode2)
void GetOwnershipRange(PetscInt &lo, PetscInt &hi)
Definition: Node.hpp:58
static bool ReplicateBool(bool flag)
Definition: PetscTools.cpp:186
BoundaryElementIterator GetBoundaryElementIteratorBegin() const
void ApplyDirichletToNonlinearResidual(const Vec currentSolution, Vec residual, DistributedVectorFactory &rFactory)
void ZeroMatrixRowsWithValueOnDiagonal(std::vector< unsigned > &rRows, double diagonalValue)
unsigned GetSize() const
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
NeumannMapIterator mLastNeumannCondition[PROBLEM_DIM]
void DefineZeroDirichletOnMeshBoundary(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, unsigned indexOfUnknown=0)
bool IsBoundaryNode() const
Definition: Node.cpp:164
static void SetElement(Vec vector, PetscInt row, double value)
void AddNeumannBoundaryCondition(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *pBoundaryElement, const AbstractBoundaryCondition< SPACE_DIM > *pBoundaryCondition, unsigned indexOfUnknown=0)
Vec & rGetDirichletBoundaryConditionsVector()
std::map< const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *, const AbstractBoundaryCondition< SPACE_DIM > * >::const_iterator NeumannMapIterator
void ZeroMatrixRowsAndColumnsWithValueOnDiagonal(std::vector< unsigned > &rRowColIndices, double diagonalValue)
bool IsGlobalIndexLocal(unsigned globalIndex)
Mat & rGetLhsMatrix()
BoundaryNodeIterator GetBoundaryNodeIteratorBegin() const
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)
Vec & rGetRhsVector()
void AddDirichletBoundaryCondition(const Node< SPACE_DIM > *pBoundaryNode, const AbstractBoundaryCondition< SPACE_DIM > *pBoundaryCondition, unsigned indexOfUnknown=0, bool checkIfBoundaryNode=true)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
std::map< const Node< SPACE_DIM > *, const Node< SPACE_DIM > * > * mpPeriodicBcMap[PROBLEM_DIM]
unsigned GetNumBoundaryNodes() const
static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
static void ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector< unsigned > &rRows, double diagonalValue)
bool Validate(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
std::vector< Node< SPACE_DIM > * >::const_iterator BoundaryNodeIterator
static void AddScaledVector(Vec y, Vec x, double scaleFactor)
#define EXCEPT_IF_NOT(test)
Definition: Exception.hpp:158
ConstBoundaryCondition< SPACE_DIM > * mpZeroBoundaryCondition
static void Zero(Vec vector)
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
void SetRhsVectorElement(PetscInt row, double value)
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)
bool IsMatrixSymmetric()
static void Finalise(Mat matrix)
void DefineZeroNeumannOnMeshBoundary(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, unsigned indexOfUnknown=0)
BoundaryNodeIterator GetBoundaryNodeIteratorEnd() const
std::vector< BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > * >::const_iterator BoundaryElementIterator
double GetValue(const ChastePoint< SPACE_DIM > &rX) const
bool HasNeumannBoundaryCondition(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *pSurfaceElement, unsigned indexOfUnknown=0)
bool mAnyNonZeroNeumannConditionsForUnknown[PROBLEM_DIM]
Vec GetMatrixRowDistributed(unsigned rowIndex)
BoundaryElementIterator GetBoundaryElementIteratorEnd() const