Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractFeCableIntegralAssembler.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 ABSTRACTFECABLEINTEGRALASSEMBLER_HPP_
37#define ABSTRACTFECABLEINTEGRALASSEMBLER_HPP_
38
39#include "AbstractFeAssemblerCommon.hpp"
40#include "GaussianQuadratureRule.hpp"
41#include "PetscVecTools.hpp"
42#include "PetscMatTools.hpp"
43
50template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
51class AbstractFeCableIntegralAssembler : public AbstractFeAssemblerCommon<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,CAN_ASSEMBLE_VECTOR,CAN_ASSEMBLE_MATRIX,INTERPOLATION_LEVEL>
52{
53protected:
55 static const unsigned CABLE_ELEMENT_DIM = 1;
56
58 static const unsigned NUM_CABLE_ELEMENT_NODES = 2;
59
62
65
68
84 const c_matrix<double, CABLE_ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
85 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES>& rReturnValue);
86
92 void DoAssemble();
93
116 // LCOV_EXCL_START
117 virtual c_matrix<double,PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES,PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES> ComputeCableMatrixTerm(
118 c_vector<double, NUM_CABLE_ELEMENT_NODES>& rPhi,
119 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES>& rGradPhi,
121 c_vector<double,PROBLEM_DIM>& rU,
122 c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
124 {
125 // If this line is reached this means this method probably hasn't been over-ridden correctly in
126 // the concrete class
128 return zero_matrix<double>(PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES,PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES);
129 }
130 // LCOV_EXCL_STOP
131
152 // LCOV_EXCL_START
153 virtual c_vector<double,PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES> ComputeCableVectorTerm(
154 c_vector<double, NUM_CABLE_ELEMENT_NODES>& rPhi,
155 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES>& rGradPhi,
157 c_vector<double,PROBLEM_DIM>& rU,
158 c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
160 {
161 // If this line is reached this means this method probably hasn't been over-ridden correctly in
162 // the concrete class
164 return zero_vector<double>(PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES);
165 }
166 // LCOV_EXCL_STOP
167
183 c_matrix<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES >& rAElem,
184 c_vector<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES>& rBElem);
185
186
195 {
196 return true;
197 }
198
199public:
200
207
212 {
213 delete mpCableQuadRule;
214 }
215};
216
217template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
220 : AbstractFeAssemblerCommon<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX, INTERPOLATION_LEVEL>(),
221 mpMesh(pMesh)
222{
223 assert(pMesh);
224 assert(CAN_ASSEMBLE_VECTOR || CAN_ASSEMBLE_MATRIX);
225 // Default to 2nd order quadrature. Our default basis functions are piecewise linear
226 // which means that we are integrating functions which in the worst case (mass matrix)
227 // are quadratic.
229
230 // Not supporting this yet - if a nonlinear assembler on cable elements is required, uncomment code
231 // in AssembleOnCableElement below (search for NONLINEAR)
232 assert(INTERPOLATION_LEVEL!=NONLINEAR);
233}
234
235
236template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
238{
239 HeartEventHandler::EventType assemble_event;
240 if (this->mAssembleMatrix)
241 {
242 assemble_event = HeartEventHandler::ASSEMBLE_SYSTEM;
243 }
244 else
245 {
246 assemble_event = HeartEventHandler::ASSEMBLE_RHS;
247 }
248
249 if (this->mAssembleMatrix && this->mMatrixToAssemble==NULL)
250 {
251 EXCEPTION("Matrix to be assembled has not been set");
252 }
253 if (this->mAssembleVector && this->mVectorToAssemble==NULL)
254 {
255 EXCEPTION("Vector to be assembled has not been set");
256 }
257
258 // This has to be below PrepareForAssembleSystem as in that method the ODEs are solved in cardiac problems
259 HeartEventHandler::BeginEvent(assemble_event);
260
261 // Zero the matrix/vector if it is to be assembled
262 if (this->mAssembleVector && this->mZeroVectorBeforeAssembly)
263 {
264 PetscVecTools::Zero(this->mVectorToAssemble);
265 }
266 if (this->mAssembleMatrix && this->mZeroMatrixBeforeAssembly)
267 {
268 PetscMatTools::Zero(this->mMatrixToAssemble);
269 }
270
271 const size_t STENCIL_SIZE=PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES;
272 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
273 c_vector<double, STENCIL_SIZE> b_elem;
274
275 // Loop over elements
276 if (this->mAssembleMatrix || this->mAssembleVector)
277 {
278 for (typename MixedDimensionMesh<CABLE_ELEMENT_DIM, SPACE_DIM>::CableElementIterator iter = mpMesh->GetCableElementIteratorBegin();
279 iter != mpMesh->GetCableElementIteratorEnd();
280 ++iter)
281 {
282 Element<CABLE_ELEMENT_DIM, SPACE_DIM>& r_element = *(*iter);
283
284 // Test for ownership first, since it's pointless to test the criterion on something which we might know nothing about.
285 if (r_element.GetOwnership() == true && ElementAssemblyCriterion(r_element)==true)
286 {
287 AssembleOnCableElement(r_element, a_elem, b_elem);
288
289 unsigned p_indices[STENCIL_SIZE];
290 r_element.GetStiffnessMatrixGlobalIndices(PROBLEM_DIM, p_indices);
291
292 if (this->mAssembleMatrix)
293 {
294 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mMatrixToAssemble, p_indices, a_elem);
295 }
296
297 if (this->mAssembleVector)
298 {
299 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mVectorToAssemble, p_indices, b_elem);
300 }
301 }
302 }
303 }
304
305 HeartEventHandler::EndEvent(assemble_event);
306}
307
309// Implementation - AssembleOnCableElement and smaller
311
312template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
314 const ChastePoint<CABLE_ELEMENT_DIM>& rPoint,
315 const c_matrix<double, CABLE_ELEMENT_DIM, SPACE_DIM>& rInverseJacobian,
316 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES>& rReturnValue)
317{
318 static c_matrix<double, CABLE_ELEMENT_DIM, NUM_CABLE_ELEMENT_NODES> grad_phi;
319
321 rReturnValue = prod(trans(rInverseJacobian), grad_phi);
322}
323
324template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
327 c_matrix<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES >& rAElem,
328 c_vector<double, PROBLEM_DIM*NUM_CABLE_ELEMENT_NODES>& rBElem)
329{
335 c_matrix<double, SPACE_DIM, CABLE_ELEMENT_DIM> jacobian;
336 c_matrix<double, CABLE_ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
337 double jacobian_determinant;
338
341 // mpMesh->GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
343 rElement.CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
344
345 if (this->mAssembleMatrix)
346 {
347 rAElem.clear();
348 }
349
350 if (this->mAssembleVector)
351 {
352 rBElem.clear();
353 }
354
355 const unsigned num_nodes = rElement.GetNumNodes();
356
357 // Allocate memory for the basis functions values and derivative values
358 c_vector<double, NUM_CABLE_ELEMENT_NODES> phi;
359 c_matrix<double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES> grad_phi;
360
361 // Loop over Gauss points
362 for (unsigned quad_index=0; quad_index < mpCableQuadRule->GetNumQuadPoints(); quad_index++)
363 {
364 const ChastePoint<CABLE_ELEMENT_DIM>& quad_point = mpCableQuadRule->rGetQuadPoint(quad_index);
365
366 CableBasisFunction::ComputeBasisFunctions(quad_point, phi);
367
368 if (this->mAssembleMatrix || INTERPOLATION_LEVEL==NONLINEAR)
369 {
370 ComputeTransformedBasisFunctionDerivatives(quad_point, inverse_jacobian, grad_phi);
371 }
372
373 // Location of the Gauss point in the original element will be stored in x
374 // Where applicable, u will be set to the value of the current solution at x
375 ChastePoint<SPACE_DIM> x(0,0,0);
376
377 c_vector<double,PROBLEM_DIM> u = zero_vector<double>(PROBLEM_DIM);
378 c_matrix<double,PROBLEM_DIM,SPACE_DIM> grad_u = zero_matrix<double>(PROBLEM_DIM,SPACE_DIM);
379
380 // Allow the concrete version of the assembler to interpolate any desired quantities
381 this->ResetInterpolatedQuantities();
382
383 // Interpolation
384 for (unsigned i=0; i<num_nodes; i++)
385 {
386 const Node<SPACE_DIM>* p_node = rElement.GetNode(i);
387
388 if (INTERPOLATION_LEVEL != CARDIAC) // don't interpolate X if cardiac problem
389 {
390 const c_vector<double, SPACE_DIM>& r_node_loc = p_node->rGetLocation();
391 // interpolate x
392 x.rGetLocation() += phi(i)*r_node_loc;
393 }
394
395 // Interpolate u and grad u if a current solution or guess exists
396 unsigned node_global_index = rElement.GetNodeGlobalIndex(i);
397 if (this->mCurrentSolutionOrGuessReplicated.GetSize() > 0)
398 {
399 for (unsigned index_of_unknown=0; index_of_unknown<(INTERPOLATION_LEVEL!=CARDIAC ? PROBLEM_DIM : 1); index_of_unknown++)
400 {
401 /*
402 * If we have a solution (e.g. this is a dynamic problem) then
403 * interpolate the value at this quadrature point.
404 *
405 * NOTE: the following assumes that if, say, there are two unknowns
406 * u and v, they are stored in the current solution vector as
407 * [U1 V1 U2 V2 ... U_n V_n].
408 */
409 double u_at_node = this->GetCurrentSolutionOrGuessValue(node_global_index, index_of_unknown);
410 u(index_of_unknown) += phi(i)*u_at_node;
411
413// if (INTERPOLATION_LEVEL==NONLINEAR) // don't need to construct grad_phi or grad_u in other cases
414// {
415// for (unsigned j=0; j<SPACE_DIM; j++)
416// {
417// grad_u(index_of_unknown,j) += grad_phi(j,i)*u_at_node;
418// }
419// }
420 }
421 }
422
423 // Allow the concrete version of the assembler to interpolate any desired quantities
424 this->IncrementInterpolatedQuantities(phi(i), p_node);
425 }
426
427 double wJ = jacobian_determinant * mpCableQuadRule->GetWeight(quad_index);
428
429 // Create rAElem and rBElem
430 if (this->mAssembleMatrix)
431 {
432 noalias(rAElem) += ComputeCableMatrixTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
433 }
434
435 if (this->mAssembleVector)
436 {
437 noalias(rBElem) += ComputeCableVectorTerm(phi, grad_phi, x, u, grad_u, &rElement) * wJ;
438 }
439 }
440}
441
442#endif /*ABSTRACTFECABLEINTEGRALASSEMBLER_HPP_*/
#define EXCEPTION(message)
#define NEVER_REACHED
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
bool GetOwnership() const
virtual c_vector< double, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES > ComputeCableVectorTerm(c_vector< double, NUM_CABLE_ELEMENT_NODES > &rPhi, c_matrix< double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< CABLE_ELEMENT_DIM, SPACE_DIM > *pElement)
virtual c_matrix< double, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES > ComputeCableMatrixTerm(c_vector< double, NUM_CABLE_ELEMENT_NODES > &rPhi, c_matrix< double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< CABLE_ELEMENT_DIM, SPACE_DIM > *pElement)
virtual bool ElementAssemblyCriterion(Element< CABLE_ELEMENT_DIM, SPACE_DIM > &rElement)
AbstractFeCableIntegralAssembler(MixedDimensionMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< CABLE_ELEMENT_DIM > &rPoint, const c_matrix< double, CABLE_ELEMENT_DIM, SPACE_DIM > &rInverseJacobian, c_matrix< double, SPACE_DIM, NUM_CABLE_ELEMENT_NODES > &rReturnValue)
MixedDimensionMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
virtual void AssembleOnCableElement(Element< CABLE_ELEMENT_DIM, SPACE_DIM > &rElement, c_matrix< double, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES > &rAElem, c_vector< double, PROBLEM_DIM *NUM_CABLE_ELEMENT_NODES > &rBElem)
void CalculateInverseJacobian(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian)
void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned *pIndices) const
static void ComputeBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, c_matrix< double, ELEMENT_DIM, ELEMENT_DIM+1 > &rReturnValue)
std::vector< Element< 1, SPACE_DIM > * >::const_iterator CableElementIterator
Definition Node.hpp:59
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
static void Zero(Mat matrix)
static void Zero(Vec vector)