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
00029 #ifndef ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
00030 #define ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
00031
00032 #include "AbstractFeAssemblerCommon.hpp"
00033 #include "GaussianQuadratureRule.hpp"
00034 #include "BoundaryConditionsContainer.hpp"
00035 #include "PetscVecTools.hpp"
00036 #include "PetscMatTools.hpp"
00037
00038
00053 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00054 class AbstractFeSurfaceIntegralAssembler : public AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,NORMAL>
00055 {
00056 protected:
00058 AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* mpMesh;
00059
00061 BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* mpBoundaryConditions;
00062
00064 GaussianQuadratureRule<ELEMENT_DIM-1>* mpSurfaceQuadRule;
00065
00067 typedef LinearBasisFunction<ELEMENT_DIM-1> SurfaceBasisFunction;
00068
00082 virtual c_vector<double, PROBLEM_DIM*ELEMENT_DIM> ComputeVectorSurfaceTerm(
00083 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00084 c_vector<double, ELEMENT_DIM>& rPhi,
00085 ChastePoint<SPACE_DIM>& rX)
00086 {
00087
00088
00089 NEVER_REACHED;
00090 return zero_vector<double>(ELEMENT_DIM*PROBLEM_DIM);
00091 }
00092
00102 virtual void AssembleOnSurfaceElement(const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00103 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem);
00104
00105
00109 void DoAssemble();
00110
00111
00112 public:
00120 AbstractFeSurfaceIntegralAssembler(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00121 BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* pBoundaryConditions,
00122 unsigned numQuadPoints = 2);
00123
00127 ~AbstractFeSurfaceIntegralAssembler();
00128
00133 void ResetBoundaryConditionsContainer(BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* pBoundaryConditions)
00134 {
00135 assert(pBoundaryConditions);
00136 this->mpBoundaryConditions = pBoundaryConditions;
00137 }
00138 };
00139
00140
00141 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00142 AbstractFeSurfaceIntegralAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AbstractFeSurfaceIntegralAssembler(
00143 AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00144 BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>* pBoundaryConditions,
00145 unsigned numQuadPoints)
00146 : AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,NORMAL>(),
00147 mpMesh(pMesh),
00148 mpBoundaryConditions(pBoundaryConditions)
00149 {
00150 assert(pMesh);
00151 assert(pBoundaryConditions);
00152 assert(numQuadPoints > 0);
00153
00154 mpSurfaceQuadRule = new GaussianQuadratureRule<ELEMENT_DIM-1>(numQuadPoints);
00155 }
00156
00157 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00158 AbstractFeSurfaceIntegralAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::~AbstractFeSurfaceIntegralAssembler()
00159 {
00160 delete mpSurfaceQuadRule;
00161 }
00162
00163
00164 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00165 void AbstractFeSurfaceIntegralAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::DoAssemble()
00166 {
00167 assert(this->mAssembleVector);
00168
00169 HeartEventHandler::BeginEvent(HeartEventHandler::NEUMANN_BCS);
00170
00171
00172 if (mpBoundaryConditions->AnyNonZeroNeumannConditions())
00173 {
00174 typename BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::NeumannMapIterator
00175 neumann_iterator = mpBoundaryConditions->BeginNeumann();
00176
00177 c_vector<double, PROBLEM_DIM*ELEMENT_DIM> b_surf_elem;
00178
00179
00180 while (neumann_iterator != mpBoundaryConditions->EndNeumann())
00181 {
00182 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& r_surf_element = *(neumann_iterator->first);
00183 AssembleOnSurfaceElement(r_surf_element, b_surf_elem);
00184
00185 const size_t STENCIL_SIZE=PROBLEM_DIM*ELEMENT_DIM;
00186 unsigned p_indices[STENCIL_SIZE];
00187 r_surf_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
00188 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_surf_elem);
00189 ++neumann_iterator;
00190 }
00191 }
00192
00193 HeartEventHandler::EndEvent(HeartEventHandler::NEUMANN_BCS);
00194 }
00195
00196
00197 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00198 void AbstractFeSurfaceIntegralAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AssembleOnSurfaceElement(
00199 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
00200 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem)
00201 {
00202 c_vector<double, SPACE_DIM> weighted_direction;
00203 double jacobian_determinant;
00204 mpMesh->GetWeightedDirectionForBoundaryElement(rSurfaceElement.GetIndex(), weighted_direction, jacobian_determinant);
00205
00206 rBSurfElem.clear();
00207
00208
00209 c_vector<double, ELEMENT_DIM> phi;
00210
00211
00212 for (unsigned quad_index=0; quad_index<mpSurfaceQuadRule->GetNumQuadPoints(); quad_index++)
00213 {
00214 const ChastePoint<ELEMENT_DIM-1>& quad_point = mpSurfaceQuadRule->rGetQuadPoint(quad_index);
00215
00216 SurfaceBasisFunction::ComputeBasisFunctions(quad_point, phi);
00217
00219
00221
00222
00223 ChastePoint<SPACE_DIM> x(0,0,0);
00224
00225 this->ResetInterpolatedQuantities();
00226 for (unsigned i=0; i<rSurfaceElement.GetNumNodes(); i++)
00227 {
00228 const c_vector<double, SPACE_DIM> node_loc = rSurfaceElement.GetNode(i)->rGetLocation();
00229 x.rGetLocation() += phi(i)*node_loc;
00230
00231
00232 this->IncrementInterpolatedQuantities(phi(i), rSurfaceElement.GetNode(i));
00233 }
00234
00235
00237
00239
00240 double wJ = jacobian_determinant * mpSurfaceQuadRule->GetWeight(quad_index);
00242 noalias(rBSurfElem) += ComputeVectorSurfaceTerm(rSurfaceElement, phi, x) * wJ;
00243 }
00244 };
00245
00246 #endif // ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_