Chaste  Release::2018.1
AbstractFunctionalCalculator.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 ABSTRACTFUNCTIONALCALCULATOR_HPP_
37 #define ABSTRACTFUNCTIONALCALCULATOR_HPP_
38 
39 #include "LinearBasisFunction.hpp"
40 #include "GaussianQuadratureRule.hpp"
41 #include "AbstractTetrahedralMesh.hpp"
42 #include "ReplicatableVector.hpp"
43 
57 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
59 {
60 private:
61 
64 
72  virtual double GetIntegrand(ChastePoint<SPACE_DIM>& rX,
73  c_vector<double,PROBLEM_DIM>& rU,
74  c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU)=0;
75 
83 
84 public:
85 
90  {
91  }
92 
104 
111 };
112 
114 // Implementation
116 
117 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
119 {
120  double result_on_element = 0;
121 
122  // Third order quadrature. Note that the functional may be non-polynomial (see documentation of class).
124 
127  double jacobian_determinant;
128  c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
129  c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
130  rElement.CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
131 
132  const unsigned num_nodes = rElement.GetNumNodes();
133 
134  // Loop over Gauss points
135  for (unsigned quad_index=0; quad_index < quad_rule.GetNumQuadPoints(); quad_index++)
136  {
137  const ChastePoint<ELEMENT_DIM>& quad_point = quad_rule.rGetQuadPoint(quad_index);
138 
139  c_vector<double, ELEMENT_DIM+1> phi;
141  c_matrix<double, ELEMENT_DIM, ELEMENT_DIM+1> grad_phi;
143 
144  // Location of the Gauss point in the original element will be stored in x
145  ChastePoint<SPACE_DIM> x(0,0,0);
146  c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
147  c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
148 
149  for (unsigned i=0; i<num_nodes; i++)
150  {
151  const c_vector<double, SPACE_DIM>& r_node_loc = rElement.GetNode(i)->rGetLocation();
152 
153  // Interpolate x
154  x.rGetLocation() += phi(i)*r_node_loc;
155 
156  // Interpolate u and grad u
157  unsigned node_global_index = rElement.GetNodeGlobalIndex(i);
158  for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
159  {
160  // NOTE - following assumes that, if say there are two unknowns u and v, they
161  // are stored in the current solution vector as
162  // [U1 V1 U2 V2 ... U_n V_n]
163  unsigned index_into_vec = PROBLEM_DIM*node_global_index + index_of_unknown;
164 
165  double u_at_node = mSolutionReplicated[index_into_vec];
166  u(index_of_unknown) += phi(i)*u_at_node;
167  // NB. grad_u is PROBLEM_DIM x SPACE_DIM but grad_phi is ELEMENT_DIM x (ELEMENT_DIM+1)
168  // Assume here that SPACE_DIM == ELEMENT_DIM and assert it in calling function
169  for (unsigned j=0; j<ELEMENT_DIM; j++)
170  {
171  grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
172  }
173  }
174  }
175 
176  double wJ = jacobian_determinant * quad_rule.GetWeight(quad_index);
177  result_on_element += GetIntegrand(x, u, grad_u) * wJ;
178  }
179 
180  return result_on_element;
181 }
182 
183 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
185 {
186  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
187  assert(solution);
188  mSolutionReplicated.ReplicatePetscVector(solution);
189  if (mSolutionReplicated.GetSize() != rMesh.GetNumNodes() * PROBLEM_DIM)
190  {
191  EXCEPTION("The solution size does not match the mesh");
192  }
193 
194  double local_result = 0;
195 
196  try
197  {
199  iter != rMesh.GetElementIteratorEnd();
200  ++iter)
201  {
202  if (rMesh.CalculateDesignatedOwnershipOfElement((*iter).GetIndex()) == true && !ShouldSkipThisElement(*iter))
203  {
204  local_result += CalculateOnElement(*iter);
205  }
206  }
207  }
208  catch (Exception &exception_in_integral)
209  {
211  throw exception_in_integral;
212  }
214 
215  double final_result;
216  MPI_Allreduce(&local_result, &final_result, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
217  return final_result;
218 }
219 
220 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
222 {
223  return false;
224 }
225 
226 #endif /*ABSTRACTFUNCTIONALCALCULATOR_HPP_*/
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
virtual bool ShouldSkipThisElement(Element< ELEMENT_DIM, SPACE_DIM > &rElement)
c_vector< double, DIM > & rGetLocation()
Definition: ChastePoint.cpp:76
unsigned GetNodeGlobalIndex(unsigned localIndex) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
double CalculateOnElement(Element< ELEMENT_DIM, SPACE_DIM > &rElement)
virtual unsigned GetNumNodes() const
unsigned GetNumQuadPoints() const
void CalculateInverseJacobian(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM, ELEMENT_DIM+1 > &rReturnValue)
double Calculate(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, Vec solution)
unsigned GetNumNodes() const
static void ReplicateException(bool flag)
Definition: PetscTools.cpp:198
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double, ELEMENT_DIM+1 > &rReturnValue)
double GetWeight(unsigned index) const
virtual bool CalculateDesignatedOwnershipOfElement(unsigned elementIndex)
virtual double GetIntegrand(ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU)=0
const ChastePoint< ELEMENT_DIM > & rGetQuadPoint(unsigned index) const