Chaste  Release::2018.1
AbstractFeVolumeIntegralAssembler.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 ABSTRACTFEVOLUMEINTEGRALASSEMBLER_HPP_
37 #define ABSTRACTFEVOLUMEINTEGRALASSEMBLER_HPP_
38 
39 #include "AbstractFeAssemblerCommon.hpp"
40 #include "GaussianQuadratureRule.hpp"
41 #include "BoundaryConditionsContainer.hpp"
42 #include "PetscVecTools.hpp"
43 #include "PetscMatTools.hpp"
44 
77 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
79  public AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX,INTERPOLATION_LEVEL>
80 {
81 protected:
84 
87 
90 
111  const c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
112  c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rReturnValue);
113 
119  void DoAssemble();
120 
121 protected:
122 
142  // LCOV_EXCL_START
143  virtual c_matrix<double,PROBLEM_DIM*(ELEMENT_DIM+1),PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeMatrixTerm(
144  c_vector<double, ELEMENT_DIM+1>& rPhi,
145  c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
147  c_vector<double,PROBLEM_DIM>& rU,
148  c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
150  {
151  // If this line is reached this means this method probably hasn't been over-ridden correctly in
152  // the concrete class
154  return zero_matrix<double>(PROBLEM_DIM*(ELEMENT_DIM+1),PROBLEM_DIM*(ELEMENT_DIM+1));
155  }
156  // LCOV_EXCL_STOP
157 
177  // LCOV_EXCL_START
178  virtual c_vector<double,PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeVectorTerm(
179  c_vector<double, ELEMENT_DIM+1>& rPhi,
180  c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
182  c_vector<double,PROBLEM_DIM>& rU,
183  c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
185  {
186  // If this line is reached this means this method probably hasn't been over-ridden correctly in
187  // the concrete class
189  return zero_vector<double>(PROBLEM_DIM*(ELEMENT_DIM+1));
190  }
191  // LCOV_EXCL_STOP
192 
193 
208  virtual void AssembleOnElement(Element<ELEMENT_DIM,SPACE_DIM>& rElement,
209  c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1) >& rAElem,
210  c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)>& rBElem);
211 
220  {
221  return true;
222  }
223 
224 
225 public:
226 
233 
238  {
239  delete mpQuadRule;
240  }
241 };
242 
243 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
246  : AbstractFeAssemblerCommon<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>(),
247  mpMesh(pMesh)
248 {
249  assert(pMesh);
250  // Default to 2nd order quadrature. Our default basis functions are piecewise linear
251  // which means that we are integrating functions which in the worst case (mass matrix)
252  // are quadratic.
254 }
255 
256 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
258 {
259  assert(this->mAssembleMatrix || this->mAssembleVector);
260 
261  HeartEventHandler::EventType assemble_event;
262  if (this->mAssembleMatrix)
263  {
264  assemble_event = HeartEventHandler::ASSEMBLE_SYSTEM;
265  }
266  else
267  {
268  assemble_event = HeartEventHandler::ASSEMBLE_RHS;
269  }
270 
271  if (this->mAssembleMatrix && this->mMatrixToAssemble==nullptr)
272  {
273  EXCEPTION("Matrix to be assembled has not been set");
274  }
275  if (this->mAssembleVector && this->mVectorToAssemble==nullptr)
276  {
277  EXCEPTION("Vector to be assembled has not been set");
278  }
279 
280  HeartEventHandler::BeginEvent(assemble_event);
281 
282  // Zero the matrix/vector if it is to be assembled
283  if (this->mAssembleVector && this->mZeroVectorBeforeAssembly)
284  {
285  PetscVecTools::Zero(this->mVectorToAssemble);
286  }
287  if (this->mAssembleMatrix && this->mZeroMatrixBeforeAssembly)
288  {
289  PetscMatTools::Zero(this->mMatrixToAssemble);
290  }
291 
292  const size_t STENCIL_SIZE=PROBLEM_DIM*(ELEMENT_DIM+1);
293  c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
294  c_vector<double, STENCIL_SIZE> b_elem;
295 
296  // Loop over elements
297  for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = mpMesh->GetElementIteratorBegin();
298  iter != mpMesh->GetElementIteratorEnd();
299  ++iter)
300  {
301  Element<ELEMENT_DIM, SPACE_DIM>& r_element = *iter;
302 
303  // Test for ownership first, since it's pointless to test the criterion on something which we might know nothing about.
304  if (r_element.GetOwnership() == true && ElementAssemblyCriterion(r_element)==true)
305  {
306  AssembleOnElement(r_element, a_elem, b_elem);
307 
308  unsigned p_indices[STENCIL_SIZE];
309  r_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
310 
311  if (this->mAssembleMatrix)
312  {
313  PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mMatrixToAssemble, p_indices, a_elem);
314  }
315 
316  if (this->mAssembleVector)
317  {
318  PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
319  }
320  }
321  }
322 
323  HeartEventHandler::EndEvent(assemble_event);
324 }
325 
326 
328 // Implementation - AssembleOnElement and smaller
330 
331 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
333  const ChastePoint<ELEMENT_DIM>& rPoint,
334  const c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
335  c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rReturnValue)
336 {
337  assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
338  static c_matrix<double, ELEMENT_DIM, ELEMENT_DIM+1> grad_phi;
339 
341  rReturnValue = prod(trans(rInverseJacobian), grad_phi);
342 }
343 
344 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
347  c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1) >& rAElem,
348  c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)>& rBElem)
349 {
355  c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
356  c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
357  double jacobian_determinant;
358 
359  mpMesh->GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
360 
361  if (this->mAssembleMatrix)
362  {
363  rAElem.clear();
364  }
365 
366  if (this->mAssembleVector)
367  {
368  rBElem.clear();
369  }
370 
371  const unsigned num_nodes = rElement.GetNumNodes();
372 
373  // Allocate memory for the basis functions values and derivative values
374  c_vector<double, ELEMENT_DIM+1> phi;
375  c_matrix<double, SPACE_DIM, ELEMENT_DIM+1> grad_phi;
376 
377  // Loop over Gauss points
378  for (unsigned quad_index=0; quad_index < mpQuadRule->GetNumQuadPoints(); quad_index++)
379  {
380  const ChastePoint<ELEMENT_DIM>& quad_point = mpQuadRule->rGetQuadPoint(quad_index);
381 
382  BasisFunction::ComputeBasisFunctions(quad_point, phi);
383 
384  if (this->mAssembleMatrix || INTERPOLATION_LEVEL==NONLINEAR)
385  {
386  ComputeTransformedBasisFunctionDerivatives(quad_point, inverse_jacobian, grad_phi);
387  }
388 
389  // Location of the Gauss point in the original element will be stored in x
390  // Where applicable, u will be set to the value of the current solution at x
391  ChastePoint<SPACE_DIM> x(0,0,0);
392 
393  c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
394  c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
395 
396  // Allow the concrete version of the assembler to interpolate any desired quantities
397  this->ResetInterpolatedQuantities();
398 
399  // Interpolation
400  for (unsigned i=0; i<num_nodes; i++)
401  {
402  const Node<SPACE_DIM>* p_node = rElement.GetNode(i);
403 
404  if (INTERPOLATION_LEVEL != CARDIAC) // don't even interpolate X if cardiac problem
405  {
406  const c_vector<double, SPACE_DIM>& r_node_loc = p_node->rGetLocation();
407  // interpolate x
408  x.rGetLocation() += phi(i)*r_node_loc;
409  }
410 
411  // Interpolate u and grad u if a current solution or guess exists
412  unsigned node_global_index = rElement.GetNodeGlobalIndex(i);
413  if (this->mCurrentSolutionOrGuessReplicated.GetSize() > 0)
414  {
415  for (unsigned index_of_unknown=0; index_of_unknown<(INTERPOLATION_LEVEL!=CARDIAC ? PROBLEM_DIM : 1); index_of_unknown++)
416  {
417  /*
418  * If we have a solution (e.g. this is a dynamic problem) then
419  * interpolate the value at this quadrature point.
420  *
421  * NOTE: the following assumes that if, say, there are two unknowns
422  * u and v, they are stored in the current solution vector as
423  * [U1 V1 U2 V2 ... U_n V_n].
424  */
425  double u_at_node = this->GetCurrentSolutionOrGuessValue(node_global_index, index_of_unknown);
426  u(index_of_unknown) += phi(i)*u_at_node;
427 
428  if (INTERPOLATION_LEVEL==NONLINEAR) // don't need to construct grad_phi or grad_u in other cases
429  {
430  for (unsigned j=0; j<SPACE_DIM; j++)
431  {
432  grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
433  }
434  }
435  }
436  }
437 
438  // Allow the concrete version of the assembler to interpolate any desired quantities
439  this->IncrementInterpolatedQuantities(phi(i), p_node);
440  if (this->mAssembleMatrix || INTERPOLATION_LEVEL==NONLINEAR)
441  {
442  this->IncrementInterpolatedGradientQuantities(grad_phi, i, p_node);
443  }
444  }
445 
446  double wJ = jacobian_determinant * mpQuadRule->GetWeight(quad_index);
447 
448  // Create rAElem and rBElem
449  if (this->mAssembleMatrix)
450  {
451  noalias(rAElem) += ComputeMatrixTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
452  }
453 
454  if (this->mAssembleVector)
455  {
456  noalias(rBElem) += ComputeVectorTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
457  }
458  }
459 }
460 
461 
462 #endif /*ABSTRACTFEVOLUMEINTEGRALASSEMBLER_HPP_*/
Definition: Node.hpp:58
c_vector< double, DIM > & rGetLocation()
Definition: ChastePoint.cpp:76
unsigned GetNodeGlobalIndex(unsigned localIndex) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
LinearBasisFunction< ELEMENT_DIM > BasisFunction
static void ComputeBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, c_matrix< double, ELEMENT_DIM, ELEMENT_DIM+1 > &rReturnValue)
#define NEVER_REACHED
Definition: Exception.hpp:206
virtual c_vector< double, PROBLEM_DIM *(ELEMENT_DIM+1)> ComputeVectorTerm(c_vector< double, ELEMENT_DIM+1 > &rPhi, c_matrix< double, SPACE_DIM, ELEMENT_DIM+1 > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< ELEMENT_DIM, SPACE_DIM > *pElement)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
bool GetOwnership() const
static void Zero(Mat matrix)
GaussianQuadratureRule< ELEMENT_DIM > * mpQuadRule
unsigned GetNumNodes() const
void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned *pIndices) const
virtual void AssembleOnElement(Element< ELEMENT_DIM, SPACE_DIM > &rElement, c_matrix< double, PROBLEM_DIM *(ELEMENT_DIM+1), PROBLEM_DIM *(ELEMENT_DIM+1) > &rAElem, c_vector< double, PROBLEM_DIM *(ELEMENT_DIM+1)> &rBElem)
virtual bool ElementAssemblyCriterion(Element< ELEMENT_DIM, SPACE_DIM > &rElement)
static void Zero(Vec vector)
void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian, c_matrix< double, SPACE_DIM, ELEMENT_DIM+1 > &rReturnValue)
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:139
virtual c_matrix< double, PROBLEM_DIM *(ELEMENT_DIM+1), PROBLEM_DIM *(ELEMENT_DIM+1)> ComputeMatrixTerm(c_vector< double, ELEMENT_DIM+1 > &rPhi, c_matrix< double, SPACE_DIM, ELEMENT_DIM+1 > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< ELEMENT_DIM, SPACE_DIM > *pElement)
unsigned GetIndex() const
AbstractFeVolumeIntegralAssembler(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)