Chaste  Release::3.4
AbstractFunctionalCalculator.hpp
1 /*
2 
3 Copyright (c) 2005-2016, 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  for (unsigned j=0; j<SPACE_DIM; j++)
168  {
169  grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
170  }
171  }
172  }
173 
174  double wJ = jacobian_determinant * quad_rule.GetWeight(quad_index);
175  result_on_element += GetIntegrand(x, u, grad_u) * wJ;
176  }
177 
178  return result_on_element;
179 }
180 
181 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
183 {
184  assert(solution);
185  mSolutionReplicated.ReplicatePetscVector(solution);
186  if (mSolutionReplicated.GetSize() != rMesh.GetNumNodes() * PROBLEM_DIM)
187  {
188  EXCEPTION("The solution size does not match the mesh");
189  }
190 
191  double local_result = 0;
192 
193  try
194  {
196  iter != rMesh.GetElementIteratorEnd();
197  ++iter)
198  {
199  if (rMesh.CalculateDesignatedOwnershipOfElement((*iter).GetIndex()) == true && !ShouldSkipThisElement(*iter))
200  {
201  local_result += CalculateOnElement(*iter);
202  }
203  }
204  }
205  catch (Exception &exception_in_integral)
206  {
208  throw exception_in_integral;
209  }
211 
212  double final_result;
213  MPI_Allreduce(&local_result, &final_result, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
214  return final_result;
215 }
216 
217 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
219 {
220  return false;
221 }
222 
223 #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