00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef ABSTRACTFUNCTIONALCALCULATOR_HPP_
00029 #define ABSTRACTFUNCTIONALCALCULATOR_HPP_
00030
00031 #include "LinearBasisFunction.hpp"
00032 #include "GaussianQuadratureRule.hpp"
00033 #include "AbstractTetrahedralMesh.hpp"
00034 #include "GaussianQuadratureRule.hpp"
00035 #include "ReplicatableVector.hpp"
00036
00037
00050 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00051 class AbstractFunctionalCalculator
00052 {
00053 private:
00054
00056 ReplicatableVector mSolutionReplicated;
00057
00065 virtual double GetIntegrand(ChastePoint<SPACE_DIM>& rX,
00066 c_vector<double,PROBLEM_DIM>& rU,
00067 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU)=0;
00068
00069 public:
00070
00074 virtual ~AbstractFunctionalCalculator()
00075 {
00076 }
00077
00088 double Calculate(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh, Vec solution);
00089
00095 double CalculateOnElement(Element<ELEMENT_DIM,SPACE_DIM>& rElement);
00096 };
00097
00098
00100
00102
00103
00104 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00105 double AbstractFunctionalCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::CalculateOnElement(Element<ELEMENT_DIM, SPACE_DIM>& rElement)
00106 {
00107 double result_on_element = 0;
00108
00109 GaussianQuadratureRule<ELEMENT_DIM> quad_rule(2);
00110
00113 double jacobian_determinant;
00114 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00115 c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
00116 rElement.CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
00117
00118 const unsigned num_nodes = rElement.GetNumNodes();
00119
00120
00121 for (unsigned quad_index=0; quad_index < quad_rule.GetNumQuadPoints(); quad_index++)
00122 {
00123 const ChastePoint<ELEMENT_DIM>& quad_point = quad_rule.rGetQuadPoint(quad_index);
00124
00125 c_vector<double, ELEMENT_DIM+1> phi;
00126 LinearBasisFunction<ELEMENT_DIM>::ComputeBasisFunctions(quad_point, phi);
00127 c_matrix<double, ELEMENT_DIM, ELEMENT_DIM+1> grad_phi;
00128 LinearBasisFunction<ELEMENT_DIM>::ComputeTransformedBasisFunctionDerivatives(quad_point, inverse_jacobian, grad_phi);
00129
00130
00131 ChastePoint<SPACE_DIM> x(0,0,0);
00132 c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
00133 c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
00134
00135 for (unsigned i=0; i<num_nodes; i++)
00136 {
00137 const c_vector<double, SPACE_DIM>& r_node_loc = rElement.GetNode(i)->rGetLocation();
00138
00139
00140 x.rGetLocation() += phi(i)*r_node_loc;
00141
00142
00143 unsigned node_global_index = rElement.GetNodeGlobalIndex(i);
00144 for (unsigned index_of_unknown=0; index_of_unknown<PROBLEM_DIM; index_of_unknown++)
00145 {
00146
00147
00148
00149 unsigned index_into_vec = PROBLEM_DIM*node_global_index + index_of_unknown;
00150
00151 double u_at_node = mSolutionReplicated[index_into_vec];
00152 u(index_of_unknown) += phi(i)*u_at_node;
00153 for (unsigned j=0; j<SPACE_DIM; j++)
00154 {
00155 grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
00156 }
00157 }
00158 }
00159
00160 double wJ = jacobian_determinant * quad_rule.GetWeight(quad_index);
00161 result_on_element += GetIntegrand(x, u, grad_u) * wJ;
00162 }
00163
00164 return result_on_element;
00165 }
00166
00167 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00168 double AbstractFunctionalCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::Calculate(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh, Vec solution)
00169 {
00170 assert(solution);
00171 mSolutionReplicated.ReplicatePetscVector(solution);
00172 if (mSolutionReplicated.GetSize() != rMesh.GetNumNodes() * PROBLEM_DIM)
00173 {
00174 EXCEPTION("The solution size does not match the mesh");
00175 }
00176
00177 double local_result = 0;
00178
00179 try
00180 {
00181 for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator iter = rMesh.GetElementIteratorBegin();
00182 iter != rMesh.GetElementIteratorEnd();
00183 ++iter)
00184 {
00185 if (rMesh.CalculateDesignatedOwnershipOfElement((*iter).GetIndex()) == true)
00186 {
00187 local_result += CalculateOnElement(*iter);
00188 }
00189 }
00190 }
00191 catch (Exception &exception_in_integral)
00192 {
00193 PetscTools::ReplicateException(true);
00194 throw exception_in_integral;
00195 }
00196 PetscTools::ReplicateException(false);
00197
00198 double final_result;
00199 MPI_Allreduce(&local_result, &final_result, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00200 return final_result;
00201 }
00202
00203 #endif