Chaste Commit::ca8ccdedf819b6e02855bc0e8e6f50bdecbc5208
IncompressibleNonlinearElasticitySolver.cpp
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/*
37 * NOTE ON COMPILATION ERRORS:
38 *
39 * This file won't compile with Intel icpc version 9.1.039, with error message:
40 * "Terminate with:
41 (0): internal error: backend signals"
42 *
43 * Try recompiling with icpc version 10.0.025.
44 */
45
46#include "IncompressibleNonlinearElasticitySolver.hpp"
47#include "LinearBasisFunction.hpp"
48#include "QuadraticBasisFunction.hpp"
49#include <algorithm>
50
51template<size_t DIM>
53 bool assembleJacobian)
54{
55 // Check we've actually been asked to do something!
56 assert(assembleResidual || assembleJacobian);
57 assert(this->mCurrentSolution.size()==this->mNumDofs);
58
59 // Zero the matrix/vector if it is to be assembled
60 if (assembleResidual)
61 {
62 PetscVecTools::Finalise(this->mResidualVector);
63 PetscVecTools::Zero(this->mResidualVector);
64 }
65 if (assembleJacobian)
66 {
67 PetscMatTools::Zero(this->mrJacobianMatrix);
68 PetscMatTools::Zero(this->mPreconditionMatrix);
69 }
70
71 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
72 // The (element) preconditioner matrix: this is the same as the jacobian, but
73 // with the mass matrix (ie \intgl phi_i phi_j) in the pressure-pressure block.
74 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem_precond;
75 c_vector<double, STENCIL_SIZE> b_elem;
76
77 // Loop over elements
78 for (typename AbstractTetrahedralMesh<DIM, DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
79 iter != this->mrQuadMesh.GetElementIteratorEnd();
80 ++iter)
81 {
82 // LCOV_EXCL_START
83 // Note: if assembleJacobian only
84 if (CommandLineArguments::Instance()->OptionExists("-mech_very_verbose") && assembleJacobian)
85 {
86 std::cout << "\r[" << PetscTools::GetMyRank() << "]: Element " << (*iter).GetIndex() << " of " << this->mrQuadMesh.GetNumElements() << std::flush;
87 }
88 // LCOV_EXCL_STOP
89
90 Element<DIM, DIM>& element = *iter;
91
92 if (element.GetOwnership() == true)
93 {
94 AssembleOnElement(element, a_elem, a_elem_precond, b_elem, assembleResidual, assembleJacobian);
95
99 //for (unsigned i=0; i<STENCIL_SIZE; i++)
100 //{
101 // for (unsigned j=0; j<STENCIL_SIZE; j++)
102 // {
103 // a_elem(i,j)=1.0;
104 // }
105 //}
106
107
109 // See comments about ordering at the elemental level vs ordering of the global mat/vec
110 // in eg AbstractContinuumMechanicsAssembler
112
113 unsigned p_indices[STENCIL_SIZE];
114 for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
115 {
116 for (unsigned j=0; j<DIM; j++)
117 {
118 // note: DIM+1 is the problem dimension (= this->mProblemDimension)
119 p_indices[DIM*i+j] = (DIM+1)*element.GetNodeGlobalIndex(i) + j;
120 }
121 }
122
123 for (unsigned i=0; i<NUM_VERTICES_PER_ELEMENT; i++)
124 {
125 // We assume the vertices are the first num_vertices nodes in the list of nodes
126 // in the element. Hence:
127 unsigned vertex_index = element.GetNodeGlobalIndex(i);
128 // note: DIM+1 is the problem dimension (= this->mProblemDimension)
129 p_indices[DIM*NUM_NODES_PER_ELEMENT + i] = (DIM+1)*vertex_index + DIM;
130 }
131
132 if (assembleJacobian)
133 {
134 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_elem);
135 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_elem_precond);
136 }
137
138 if (assembleResidual)
139 {
140 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mResidualVector, p_indices, b_elem);
141 }
142 }
143 }
144
145 // Loop over specified boundary elements and compute surface traction terms
146 c_vector<double, BOUNDARY_STENCIL_SIZE> b_boundary_elem; // note BOUNDARY_STENCIL_SIZE = DIM*NUM_BOUNDARY_NODES, as all pressure block is zero
147 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE> a_boundary_elem;
148
149 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() != NO_TRACTIONS)
150 {
151 for (unsigned bc_index=0; bc_index<this->mrProblemDefinition.rGetTractionBoundaryElements().size(); bc_index++)
152 {
153 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(this->mrProblemDefinition.rGetTractionBoundaryElements()[bc_index]);
154
155 // If the BCs are tractions applied on a given surface, the boundary integral is independent of u,
156 // so a_boundary_elem will be zero (no contribution to jacobian).
157 // If the BCs are normal pressure applied to the deformed body, the boundary depends on the deformation,
158 // so there is a contribution to the jacobian, and a_boundary_elem is non-zero. Note however that
159 // the AssembleOnBoundaryElement() method might decide not to include this, as it can actually
160 // cause divergence if the current guess is not close to the true solution
161 this->AssembleOnBoundaryElement(r_boundary_element, a_boundary_elem, b_boundary_elem, assembleResidual, assembleJacobian, bc_index);
162
163 unsigned p_indices[BOUNDARY_STENCIL_SIZE];
164 for (unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
165 {
166 for (unsigned j=0; j<DIM; j++)
167 {
168 // note: DIM+1, on the right hand side of the below, is the problem dimension (= this->mProblemDimension)
169 p_indices[DIM*i+j] = (DIM+1)*r_boundary_element.GetNodeGlobalIndex(i) + j;
170 }
171 }
172
173 if (assembleJacobian)
174 {
175 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_boundary_elem);
176 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_boundary_elem);
177 }
178
179 if (assembleResidual)
180 {
181 PetscVecTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mResidualVector, p_indices, b_boundary_elem);
182 }
183 }
184 }
185
186
187 if (assembleResidual)
188 {
189 PetscVecTools::Finalise(this->mResidualVector);
190 }
191 if (assembleJacobian)
192 {
193 PetscMatTools::SwitchWriteMode(this->mrJacobianMatrix);
194 PetscMatTools::SwitchWriteMode(this->mPreconditionMatrix);
195 }
196
197 if (assembleJacobian)
198 {
199 this->AddIdentityBlockForDummyPressureVariables(NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING);
200 }
201 else if (assembleResidual)
202 {
203 this->AddIdentityBlockForDummyPressureVariables(NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY);
204 }
205
206 this->FinishAssembleSystem(assembleResidual, assembleJacobian);
207}
208
209template<size_t DIM>
211 Element<DIM, DIM>& rElement,
212 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
213 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElemPrecond,
214 c_vector<double, STENCIL_SIZE>& rBElem,
215 bool assembleResidual,
216 bool assembleJacobian)
217{
218 static c_matrix<double,DIM,DIM> jacobian;
219 static c_matrix<double,DIM,DIM> inverse_jacobian;
220 double jacobian_determinant;
221
222 this->mrQuadMesh.GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
223
224 if (assembleJacobian)
225 {
226 rAElem.clear();
227 rAElemPrecond.clear();
228 }
229
230 if (assembleResidual)
231 {
232 rBElem.clear();
233 }
234
235 // Get the current displacement at the nodes
236 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
237 static c_vector<double,NUM_VERTICES_PER_ELEMENT> element_current_pressures;
238 for (unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
239 {
240 for (unsigned JJ=0; JJ<DIM; JJ++)
241 {
242 // note: DIM+1, on the right hand side of the below, is the problem dimension (= this->mProblemDimension)
243 element_current_displacements(JJ,II) = this->mCurrentSolution[(DIM+1)*rElement.GetNodeGlobalIndex(II) + JJ];
244 }
245 }
246
247 // Get the current pressure at the vertices
248 for (unsigned II=0; II<NUM_VERTICES_PER_ELEMENT; II++)
249 {
250 // At the moment we assume the vertices are the first num_vertices nodes in the list of nodes
251 // in the mesh. Hence:
252 unsigned vertex_index = rElement.GetNodeGlobalIndex(II);
253
254 // note: DIM+1, on the right hand side of the below, is the problem dimension (= this->mProblemDimension)
255 element_current_pressures(II) = this->mCurrentSolution[(DIM+1)*vertex_index + DIM];
256 }
257
258 // Allocate memory for the basis functions values and derivative values
259 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
260 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
261 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
262 static c_matrix<double, NUM_NODES_PER_ELEMENT, DIM> trans_grad_quad_phi;
263
264 // Get the material law
266 = this->mrProblemDefinition.GetIncompressibleMaterialLaw(rElement.GetIndex());
267
268 static c_matrix<double,DIM,DIM> grad_u; // grad_u = (du_i/dX_M)
269
270 static c_matrix<double,DIM,DIM> F; // the deformation gradient, F = dx/dX, F_{iM} = dx_i/dX_M
271 static c_matrix<double,DIM,DIM> C; // Green deformation tensor, C = F^T F
272 static c_matrix<double,DIM,DIM> inv_C; // inverse(C)
273 static c_matrix<double,DIM,DIM> inv_F; // inverse(F)
274 static c_matrix<double,DIM,DIM> T; // Second Piola-Kirchoff stress tensor (= dW/dE = 2dW/dC)
275
276 static c_matrix<double,DIM,DIM> F_T; // F*T
277 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> F_T_grad_quad_phi; // F*T*grad_quad_phi
278
279 c_vector<double,DIM> body_force;
280
281 static FourthOrderTensor<DIM,DIM,DIM,DIM> dTdE; // dTdE(M,N,P,Q) = dT_{MN}/dE_{PQ}
282 static FourthOrderTensor<DIM,DIM,DIM,DIM> dSdF; // dSdF(M,i,N,j) = dS_{Mi}/dF_{jN}
283
286
287 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> temp_matrix;
288 static c_matrix<double,NUM_NODES_PER_ELEMENT,DIM> grad_quad_phi_times_invF;
289
290
291 if (this->mSetComputeAverageStressPerElement)
292 {
293 this->mAverageStressesPerElement[rElement.GetIndex()] = zero_vector<double>(DIM*(DIM+1)/2);
294 }
295
296 // Loop over Gauss points
297 for (unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
298 {
299 // This is needed by the cardiac mechanics solver
300 unsigned current_quad_point_global_index = rElement.GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
301 + quadrature_index;
302
303 double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
304
305 const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
306
307 // Set up basis function information
308 LinearBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, linear_phi);
309 QuadraticBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, quad_phi);
310 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
311 trans_grad_quad_phi = trans(grad_quad_phi);
312
313 // Get the body force, interpolating X if necessary
314 if (assembleResidual)
315 {
316 switch (this->mrProblemDefinition.GetBodyForceType())
317 {
318 case FUNCTIONAL_BODY_FORCE:
319 {
320 c_vector<double,DIM> X = zero_vector<double>(DIM);
321 // interpolate X (using the vertices and the /linear/ bases, as no curvilinear elements
322 for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
323 {
324 X += linear_phi(node_index)*this->mrQuadMesh.GetNode( rElement.GetNodeGlobalIndex(node_index) )->rGetLocation();
325 }
326 body_force = this->mrProblemDefinition.EvaluateBodyForceFunction(X, this->mCurrentTime);
327 break;
328 }
329 case CONSTANT_BODY_FORCE:
330 {
331 body_force = this->mrProblemDefinition.GetConstantBodyForce();
332 break;
333 }
334 default:
336 }
337 }
338
339 // Interpolate grad_u and p
340 grad_u = zero_matrix<double>(DIM,DIM);
341
342 for (unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
343 {
344 for (unsigned i=0; i<DIM; i++)
345 {
346 for (unsigned M=0; M<DIM; M++)
347 {
348 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
349 }
350 }
351 }
352
353 double pressure = 0;
354 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
355 {
356 pressure += linear_phi(vertex_index)*element_current_pressures(vertex_index);
357 }
358
359 // Calculate C, inv(C) and T
360 for (unsigned i=0; i<DIM; i++)
361 {
362 for (unsigned M=0; M<DIM; M++)
363 {
364 F(i,M) = (i==M?1:0) + grad_u(i,M);
365 }
366 }
367
368 C = prod(trans(F),F);
369 inv_C = Inverse(C);
370 inv_F = Inverse(F);
371
372 double detF = Determinant(F);
373
374 // Compute the passive stress, and dTdE corresponding to passive stress
375 this->SetupChangeOfBasisMatrix(rElement.GetIndex(), current_quad_point_global_index);
376 p_material_law->SetChangeOfBasisMatrix(this->mChangeOfBasisMatrix);
377 p_material_law->ComputeStressAndStressDerivative(C, inv_C, pressure, T, dTdE, assembleJacobian);
378
379 if (this->mIncludeActiveTension)
380 {
381 // Add any active stresses, if there are any. Requires subclasses to overload this method,
382 // see for example the cardiac mechanics assemblers.
383 this->AddActiveStressAndStressDerivative(C, rElement.GetIndex(), current_quad_point_global_index,
384 T, dTdE, assembleJacobian);
385 }
386
387 if (this->mSetComputeAverageStressPerElement)
388 {
389 this->AddStressToAverageStressPerElement(T,rElement.GetIndex());
390 }
391
392 // Residual vector
393 if (assembleResidual)
394 {
395 F_T = prod(F,T);
396 F_T_grad_quad_phi = prod(F_T, grad_quad_phi);
397
398 for (unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
399 {
400 unsigned spatial_dim = index%DIM;
401 unsigned node_index = (index-spatial_dim)/DIM;
402
403 rBElem(index) += - this->mrProblemDefinition.GetDensity()
404 * body_force(spatial_dim)
405 * quad_phi(node_index)
406 * wJ;
407
408 // The T(M,N)*F(spatial_dim,M)*grad_quad_phi(N,node_index) term
409 rBElem(index) += F_T_grad_quad_phi(spatial_dim,node_index)
410 * wJ;
411 }
412
413 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
414 {
415 rBElem( NUM_NODES_PER_ELEMENT*DIM + vertex_index ) += linear_phi(vertex_index)
416 * (detF - 1)
417 * wJ;
418 }
419 }
420
421 // Jacobian matrix
422 if (assembleJacobian)
423 {
424 // Save trans(grad_quad_phi) * invF
425 grad_quad_phi_times_invF = prod(trans_grad_quad_phi, inv_F);
426
428 // Set up the tensor dSdF
429 //
430 // dSdF as a function of T and dTdE (which is what the material law returns) is given by:
431 //
432 // dS_{Mi}/dF_{jN} = (dT_{MN}/dC_{PQ}+dT_{MN}/dC_{PQ}) F{iP} F_{jQ} + T_{MN} delta_{ij}
433 //
434 // todo1: this should probably move into the material law (but need to make sure
435 // memory is handled efficiently
436 // todo2: get material law to return this immediately, not dTdE
438
439 // Set up the tensor 0.5(dTdE(M,N,P,Q) + dTdE(M,N,Q,P))
440 for (unsigned M=0; M<DIM; M++)
441 {
442 for (unsigned N=0; N<DIM; N++)
443 {
444 for (unsigned P=0; P<DIM; P++)
445 {
446 for (unsigned Q=0; Q<DIM; Q++)
447 {
448 // This is NOT dSdF, just using this as storage space
449 dSdF(M,N,P,Q) = 0.5*(dTdE(M,N,P,Q) + dTdE(M,N,Q,P));
450 }
451 }
452 }
453 }
454
455 // This is NOT dTdE, just reusing memory. A^{MdPQ} = F^d_N * dTdE_sym^{MNPQ}
456 dTdE.template SetAsContractionOnSecondDimension<DIM>(F, dSdF);
457
458 // dSdF{MdPe} := F^d_N * F^e_Q * dTdE_sym^{MNPQ}
459 dSdF.template SetAsContractionOnFourthDimension<DIM>(F, dTdE);
460
461 // Now add the T_{MN} delta_{ij} term
462 for (unsigned M=0; M<DIM; M++)
463 {
464 for (unsigned N=0; N<DIM; N++)
465 {
466 for (unsigned i=0; i<DIM; i++)
467 {
468 dSdF(M,i,N,i) += T(M,N);
469 }
470 }
471 }
472
474 // Set up the tensor
475 // dSdF_quad_quad(node_index1, spatial_dim1, node_index2, spatial_dim2)
476 // = dS_{M,spatial_dim1}/d_F{spatial_dim2,N}
477 // * grad_quad_phi(M,node_index1)
478 // * grad_quad_phi(P,node_index2)
479 //
480 // = dSdF(M,spatial_index1,N,spatial_index2)
481 // * grad_quad_phi(M,node_index1)
482 // * grad_quad_phi(P,node_index2)
483 //
485 temp_tensor.template SetAsContractionOnFirstDimension<DIM>( trans_grad_quad_phi, dSdF );
486 dSdF_quad_quad.template SetAsContractionOnThirdDimension<DIM>( trans_grad_quad_phi, temp_tensor );
487
488 for (unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
489 {
490 unsigned spatial_dim1 = index1%DIM;
491 unsigned node_index1 = (index1-spatial_dim1)/DIM;
492
493
494 for (unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
495 {
496 unsigned spatial_dim2 = index2%DIM;
497 unsigned node_index2 = (index2-spatial_dim2)/DIM;
498
499 // The dSdF*grad_quad_phi*grad_quad_phi term
500 rAElem(index1,index2) += dSdF_quad_quad(node_index1,spatial_dim1,node_index2,spatial_dim2)
501 * wJ;
502 }
503
504 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
505 {
506 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
507
508 // The -invF(M,spatial_dim1)*grad_quad_phi(M,node_index1)*linear_phi(vertex_index) term
509 rAElem(index1,index2) += - grad_quad_phi_times_invF(node_index1,spatial_dim1)
510 * linear_phi(vertex_index)
511 * wJ;
512 }
513 }
514
515 for (unsigned vertex_index=0; vertex_index<NUM_VERTICES_PER_ELEMENT; vertex_index++)
516 {
517 unsigned index1 = NUM_NODES_PER_ELEMENT*DIM + vertex_index;
518
519 for (unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
520 {
521 unsigned spatial_dim2 = index2%DIM;
522 unsigned node_index2 = (index2-spatial_dim2)/DIM;
523
524 // Same as (negative of) the opposite block (ie a few lines up), except for detF
525 rAElem(index1,index2) += detF
526 * grad_quad_phi_times_invF(node_index2,spatial_dim2)
527 * linear_phi(vertex_index)
528 * wJ;
529 }
530
532 // Preconditioner matrix
533 // Fill the mass matrix (ie \intgl phi_i phi_j) in the
534 // pressure-pressure block. Note, the rest of the
535 // entries are filled in at the end
537 for (unsigned vertex_index2=0; vertex_index2<NUM_VERTICES_PER_ELEMENT; vertex_index2++)
538 {
539 unsigned index2 = NUM_NODES_PER_ELEMENT*DIM + vertex_index2;
540 rAElemPrecond(index1,index2) += linear_phi(vertex_index)
541 * linear_phi(vertex_index2)
542 * wJ;
543 }
544 }
545 }
546 }
547
548 if (assembleJacobian)
549 {
550 if (this->mPetscDirectSolve)
551 {
552 // Petsc will do an LU factorisation of the preconditioner, which we
553 // set equal to [ A B1^T ]
554 // [ B2 M ]
555 // The reason for the mass matrix is to avoid zeros on the diagonal
556 rAElemPrecond = rAElemPrecond + rAElem;
557 }
558 else
559 {
560 // Fill in the other blocks of the preconditioner matrix, by adding
561 // the Jacobian matrix (this doesn't effect the pressure-pressure block
562 // of rAElemPrecond as the pressure-pressure block of rAElem is zero),
563 // and the zero a block.
564 //
565 // The following altogether gives the preconditioner [ A B1^T ]
566 // [ 0 M ]
567 rAElemPrecond = rAElemPrecond + rAElem;
568
569 for (unsigned i=NUM_NODES_PER_ELEMENT*DIM; i<STENCIL_SIZE; i++)
570 {
571 for (unsigned j=0; j<NUM_NODES_PER_ELEMENT*DIM; j++)
572 {
573 rAElemPrecond(i,j) = 0.0;
574 }
575 }
576 }
577 }
578
579 if (this->mSetComputeAverageStressPerElement)
580 {
581 for (unsigned i=0; i<DIM*(DIM+1)/2; i++)
582 {
583 this->mAverageStressesPerElement[rElement.GetIndex()](i) /= this->mpQuadratureRule->GetNumQuadPoints();
584 }
585 }
586}
587
588template<size_t DIM>
590{
591 this->mCurrentSolution.resize(this->mNumDofs, 0.0);
592
593 for (typename AbstractTetrahedralMesh<DIM, DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
594 iter != this->mrQuadMesh.GetElementIteratorEnd();
595 ++iter)
596 {
598 double zero_strain_pressure
599 = this->mrProblemDefinition.GetIncompressibleMaterialLaw(iter->GetIndex())->GetZeroStrainPressure();
600
601
602 // Loop over vertices and set pressure solution to be zero-strain-pressure
603 for (unsigned j=0; j<NUM_VERTICES_PER_ELEMENT; j++)
604 {
605 // We assume the vertices are the first num_vertices nodes in the list of nodes
606 // in the element. Hence:
607 unsigned vertex_index = iter->GetNodeGlobalIndex(j);
608 // note: DIM+1 is the problem dimension (= this->mProblemDimension)
609 this->mCurrentSolution[ (DIM+1)*vertex_index + DIM ] = zero_strain_pressure;
610 }
611 }
612}
613
614template<size_t DIM>
617 SolidMechanicsProblemDefinition<DIM>& rProblemDefinition,
618 std::string outputDirectory)
619 : AbstractNonlinearElasticitySolver<DIM>(rQuadMesh,
620 rProblemDefinition,
621 outputDirectory,
622 INCOMPRESSIBLE)
623{
624 if (rProblemDefinition.GetCompressibilityType() != INCOMPRESSIBLE)
625 {
626 EXCEPTION("SolidMechanicsProblemDefinition object contains compressible material laws");
627 }
628
630}
631
632// Explicit instantiation
#define EXCEPTION(message)
#define NEVER_REACHED
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
bool GetOwnership() const
unsigned GetIndex() const
virtual void ComputeStressAndStressDerivative(c_matrix< double, DIM, DIM > &rC, c_matrix< double, DIM, DIM > &rInvC, double pressure, c_matrix< double, DIM, DIM > &rT, FourthOrderTensor< DIM, DIM, DIM, DIM > &rDTdE, bool computeDTdE)=0
void SetChangeOfBasisMatrix(c_matrix< double, DIM, DIM > &rChangeOfBasisMatrix)
bool OptionExists(const std::string &rOption)
static CommandLineArguments * Instance()
IncompressibleNonlinearElasticitySolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)
void AssembleSystem(bool assembleResidual, bool assembleJacobian)
virtual void AssembleOnElement(Element< DIM, DIM > &rElement, c_matrix< double, STENCIL_SIZE, STENCIL_SIZE > &rAElem, c_matrix< double, STENCIL_SIZE, STENCIL_SIZE > &rAElemPrecond, c_vector< double, STENCIL_SIZE > &rBElem, bool assembleResidual, bool assembleJacobian)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double, ELEMENT_DIM+1 > &rReturnValue)
static void Zero(Mat matrix)
static void SwitchWriteMode(Mat matrix)
static unsigned GetMyRank()
static void Finalise(Vec vector)
static void Zero(Vec vector)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)