Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractFeSurfaceIntegralAssembler.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifndef ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
37#define ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
38
39#include "AbstractFeAssemblerCommon.hpp"
40#include "GaussianQuadratureRule.hpp"
41#include "BoundaryConditionsContainer.hpp"
42#include "PetscVecTools.hpp"
43#include "PetscMatTools.hpp"
44
45
60template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
61class AbstractFeSurfaceIntegralAssembler : public AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,NORMAL>
62{
63protected:
66
69
72
75
89 // LCOV_EXCL_START
90 virtual c_vector<double, PROBLEM_DIM*ELEMENT_DIM> ComputeVectorSurfaceTerm(
91 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
92 c_vector<double, ELEMENT_DIM>& rPhi,
94 {
95 // If this line is reached this means this method probably hasn't been over-ridden correctly in
96 // the concrete class
98 return zero_vector<double>(ELEMENT_DIM*PROBLEM_DIM);
99 }
100 // LCOV_EXCL_STOP
101
112 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem);
113
114
119
120
121public:
130
135
141 {
142 assert(pBoundaryConditions);
143 this->mpBoundaryConditions = pBoundaryConditions;
144 }
145};
146
147
148template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
152 : AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,false,NORMAL>(),
153 mpMesh(pMesh),
154 mpBoundaryConditions(pBoundaryConditions)
155{
156 assert(pMesh);
157 assert(pBoundaryConditions);
158 // Default to 2nd order quadrature. Our default basis functions are piecewise linear
159 // which means that we are integrating functions which in the worst case (mass matrix)
160 // are quadratic.
161 mpSurfaceQuadRule = new GaussianQuadratureRule<ELEMENT_DIM-1>(2);
162}
163
164template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
169
170
171template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
173{
174 assert(this->mAssembleVector);
175
176 HeartEventHandler::BeginEvent(HeartEventHandler::NEUMANN_BCS);
177
178 // Loop over surface elements with non-zero Neumann boundary conditions
179 if (mpBoundaryConditions->AnyNonZeroNeumannConditions())
180 {
182 neumann_iterator = mpBoundaryConditions->BeginNeumann();
183
184 c_vector<double, PROBLEM_DIM*ELEMENT_DIM> b_surf_elem;
185
186 // Iterate over defined conditions
187 while (neumann_iterator != mpBoundaryConditions->EndNeumann())
188 {
189 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& r_surf_element = *(neumann_iterator->first);
190 AssembleOnSurfaceElement(r_surf_element, b_surf_elem);
191
192 const size_t STENCIL_SIZE=PROBLEM_DIM*ELEMENT_DIM; // problem_dim*num_nodes_on_surface_element
193 unsigned p_indices[STENCIL_SIZE];
194 r_surf_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
195 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_surf_elem);
196 ++neumann_iterator;
197 }
198 }
199
200 HeartEventHandler::EndEvent(HeartEventHandler::NEUMANN_BCS);
201}
202
203
204template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
206 const BoundaryElement<ELEMENT_DIM-1,SPACE_DIM>& rSurfaceElement,
207 c_vector<double, PROBLEM_DIM*ELEMENT_DIM>& rBSurfElem)
208{
209 c_vector<double, SPACE_DIM> weighted_direction;
210 double jacobian_determinant;
211 mpMesh->GetWeightedDirectionForBoundaryElement(rSurfaceElement.GetIndex(), weighted_direction, jacobian_determinant);
212
213 rBSurfElem.clear();
214
215 // Allocate memory for the basis function values
216 c_vector<double, ELEMENT_DIM> phi;
217
218 // Loop over Gauss points
219 for (unsigned quad_index=0; quad_index<mpSurfaceQuadRule->GetNumQuadPoints(); quad_index++)
220 {
221 const ChastePoint<ELEMENT_DIM-1>& quad_point = mpSurfaceQuadRule->rGetQuadPoint(quad_index);
222
223 SurfaceBasisFunction::ComputeBasisFunctions(quad_point, phi);
224
226 // Interpolation: X only
228
229 // The location of the Gauss point in the original element will be stored in x
230 ChastePoint<SPACE_DIM> x(0,0,0);
231
232 this->ResetInterpolatedQuantities();
233 for (unsigned i=0; i<rSurfaceElement.GetNumNodes(); i++)
234 {
235 const c_vector<double, SPACE_DIM> node_loc = rSurfaceElement.GetNode(i)->rGetLocation();
236 x.rGetLocation() += phi(i)*node_loc;
237
238 // Allow the concrete version of the assembler to interpolate any desired quantities
239 this->IncrementInterpolatedQuantities(phi(i), rSurfaceElement.GetNode(i));
240 }
241
242
244 // Create elemental contribution
246
247 double wJ = jacobian_determinant * mpSurfaceQuadRule->GetWeight(quad_index);
249 noalias(rBSurfElem) += ComputeVectorSurfaceTerm(rSurfaceElement, phi, x) * wJ;
250 }
251}
252
253#endif // ABSTRACTFESURFACENTEGRALASSEMBLER_HPP_
#define NEVER_REACHED
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)
GaussianQuadratureRule< ELEMENT_DIM-1 > * mpSurfaceQuadRule
void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned *pIndices) const
std::map< constBoundaryElement< ELEMENT_DIM-1, SPACE_DIM > *, constAbstractBoundaryCondition< SPACE_DIM > * >::const_iterator NeumannMapIterator