36#ifndef ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
37#define ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
39#include "AbstractFeAssemblerCommon.hpp"
40#include "GaussianQuadratureRule.hpp"
41#include "BoundaryConditionsContainer.hpp"
42#include "PetscVecTools.hpp"
43#include "PetscMatTools.hpp"
60template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
92 c_vector<double, ELEMENT_DIM>& rPhi,
98 return zero_vector<double>(ELEMENT_DIM*PROBLEM_DIM);
112 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem);
142 assert(pBoundaryConditions);
143 this->mpBoundaryConditions = pBoundaryConditions;
148template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
154 mpBoundaryConditions(pBoundaryConditions)
157 assert(pBoundaryConditions);
164template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
167 delete mpSurfaceQuadRule;
171template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
174 assert(this->mAssembleVector);
179 if (mpBoundaryConditions->AnyNonZeroNeumannConditions())
182 neumann_iterator = mpBoundaryConditions->BeginNeumann();
184 c_vector<double, PROBLEM_DIM*ELEMENT_DIM> b_surf_elem;
187 while (neumann_iterator != mpBoundaryConditions->
EndNeumann())
189 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& r_surf_element = *(neumann_iterator->first);
190 AssembleOnSurfaceElement(r_surf_element, b_surf_elem);
192 const size_t STENCIL_SIZE=PROBLEM_DIM*ELEMENT_DIM;
193 unsigned p_indices[STENCIL_SIZE];
195 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_surf_elem);
204template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
207 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem)
209 c_vector<double, SPACE_DIM> weighted_direction;
210 double jacobian_determinant;
211 mpMesh->GetWeightedDirectionForBoundaryElement(rSurfaceElement.
GetIndex(), weighted_direction, jacobian_determinant);
216 c_vector<double, ELEMENT_DIM> phi;
219 for (
unsigned quad_index=0; quad_index<mpSurfaceQuadRule->GetNumQuadPoints(); quad_index++)
221 const ChastePoint<ELEMENT_DIM-1>& quad_point = mpSurfaceQuadRule->rGetQuadPoint(quad_index);
223 SurfaceBasisFunction::ComputeBasisFunctions(quad_point, phi);
232 this->ResetInterpolatedQuantities();
233 for (
unsigned i=0; i<rSurfaceElement.
GetNumNodes(); i++)
235 const c_vector<double, SPACE_DIM> node_loc = rSurfaceElement.
GetNode(i)->rGetLocation();
236 x.rGetLocation() += phi(i)*node_loc;
239 this->IncrementInterpolatedQuantities(phi(i), rSurfaceElement.
GetNode(i));
247 double wJ = jacobian_determinant * mpSurfaceQuadRule->GetWeight(quad_index);
249 noalias(rBSurfElem) += ComputeVectorSurfaceTerm(rSurfaceElement, phi, x) * wJ;
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetNumNodes() const
unsigned GetIndex() const
LinearBasisFunction< ELEMENT_DIM-1 > SurfaceBasisFunction
BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpBoundaryConditions
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
virtual c_vector< double, PROBLEM_DIM *ELEMENT_DIM > ComputeVectorSurfaceTerm(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > &rSurfaceElement, c_vector< double, ELEMENT_DIM > &rPhi, ChastePoint< SPACE_DIM > &rX)
AbstractFeSurfaceIntegralAssembler(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > *pBoundaryConditions)
void ResetBoundaryConditionsContainer(BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > *pBoundaryConditions)
virtual void AssembleOnSurfaceElement(const BoundaryElement< ELEMENT_DIM-1, SPACE_DIM > &rSurfaceElement, c_vector< double, PROBLEM_DIM *ELEMENT_DIM > &rBSurfElem)
virtual ~AbstractFeSurfaceIntegralAssembler()
GaussianQuadratureRule< ELEMENT_DIM-1 > * mpSurfaceQuadRule
void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned *pIndices) const
NeumannMapIterator EndNeumann()
std::map< constBoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *, constAbstractBoundaryCondition< SPACE_DIM > * >::const_iterator NeumannMapIterator
static void BeginEvent(unsigned event)
static void EndEvent(unsigned event)