Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
CompressibleNonlinearElasticitySolver.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 * (The following applies to IncompressibleNonlinearElasticityAssembler; possibly/probably holds for this class too).
40 *
41 * This file won't compile with Intel icpc version 9.1.039, with error message:
42 * "Terminate with:
43 (0): internal error: backend signals"
44 *
45 * Try recompiling with icpc version 10.0.025.
46 */
47
48#include "CompressibleNonlinearElasticitySolver.hpp"
49#include "LinearBasisFunction.hpp"
50#include "QuadraticBasisFunction.hpp"
51#include <algorithm>
52
53template<size_t DIM>
55 bool assembleJacobian)
56{
57 // Check we've actually been asked to do something!
58 assert(assembleResidual || assembleJacobian);
59 assert(this->mCurrentSolution.size()==this->mNumDofs);
60
61 // Zero the matrix/vector if it is to be assembled
62 if (assembleResidual)
63 {
64 PetscVecTools::Finalise(this->mResidualVector);
65 PetscVecTools::Zero(this->mResidualVector);
66 }
67 if (assembleJacobian)
68 {
69 PetscMatTools::Zero(this->mrJacobianMatrix);
70 PetscMatTools::Zero(this->mPreconditionMatrix);
71 }
72
73 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem;
74 // The (element) preconditioner matrix: this is the same as the jacobian, but
75 // with the mass matrix (ie \intgl phi_i phi_j) in the pressure-pressure block.
76 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE> a_elem_precond;
77 c_vector<double, STENCIL_SIZE> b_elem;
78
79 // Loop over elements
80 for (typename AbstractTetrahedralMesh<DIM, DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
81 iter != this->mrQuadMesh.GetElementIteratorEnd();
82 ++iter)
83 {
84 Element<DIM, DIM>& element = *iter;
85
86 if (element.GetOwnership() == true)
87 {
88 // LCOV_EXCL_START
89 // note: if assembleJacobian only
90 if (CommandLineArguments::Instance()->OptionExists("-mech_very_verbose") && assembleJacobian)
91 {
92 std::cout << "\r[" << PetscTools::GetMyRank() << "]: Element " << (*iter).GetIndex() << " of " << this->mrQuadMesh.GetNumElements() << std::flush;
93 }
94 // LCOV_EXCL_STOP
95
96 AssembleOnElement(element, a_elem, a_elem_precond, b_elem, assembleResidual, assembleJacobian);
97
101 //for (unsigned i=0; i<STENCIL_SIZE; i++)
102 //{
103 // for (unsigned j=0; j<STENCIL_SIZE; j++)
104 // {
105 // a_elem(i,j)=1.0;
106 // }
107 //}
108
109 unsigned p_indices[STENCIL_SIZE];
110 for (unsigned i=0; i<NUM_NODES_PER_ELEMENT; i++)
111 {
112 for (unsigned j=0; j<DIM; j++)
113 {
114 p_indices[DIM*i+j] = DIM*element.GetNodeGlobalIndex(i) + j;
115 }
116 }
117
118 if (assembleJacobian)
119 {
120 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_elem);
121 PetscMatTools::AddMultipleValues<STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_elem_precond);
122 }
123
124 if (assembleResidual)
125 {
126 PetscVecTools::AddMultipleValues<STENCIL_SIZE>(this->mResidualVector, p_indices, b_elem);
127 }
128 }
129 }
130
131 // Loop over specified boundary elements and compute surface traction terms
132 c_vector<double, BOUNDARY_STENCIL_SIZE> b_boundary_elem;
133 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE> a_boundary_elem;
134 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() != NO_TRACTIONS)
135 {
136 for (unsigned bc_index=0; bc_index<this->mrProblemDefinition.rGetTractionBoundaryElements().size(); bc_index++)
137 {
138 BoundaryElement<DIM-1,DIM>& r_boundary_element = *(this->mrProblemDefinition.rGetTractionBoundaryElements()[bc_index]);
139
140 // If the BCs are tractions applied on a given surface, the boundary integral is independent of u,
141 // so a_boundary_elem will be zero (no contribution to jacobian).
142 // If the BCs are normal pressure applied to the deformed body, the boundary depends on the deformation,
143 // so there is a contribution to the jacobian, and a_boundary_elem is non-zero. Note however that
144 // the AssembleOnBoundaryElement() method might decide not to include this, as it can actually
145 // cause divergence if the current guess is not close to the true solution
146
147 this->AssembleOnBoundaryElement(r_boundary_element, a_boundary_elem, b_boundary_elem, assembleResidual, assembleJacobian, bc_index);
148
149 unsigned p_indices[BOUNDARY_STENCIL_SIZE];
150 for (unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
151 {
152 for (unsigned j=0; j<DIM; j++)
153 {
154 p_indices[DIM*i+j] = DIM*r_boundary_element.GetNodeGlobalIndex(i) + j;
155 }
156 }
157
158 if (assembleJacobian)
159 {
160 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mrJacobianMatrix, p_indices, a_boundary_elem);
161 PetscMatTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mPreconditionMatrix, p_indices, a_boundary_elem);
162 }
163
164 if (assembleResidual)
165 {
166 PetscVecTools::AddMultipleValues<BOUNDARY_STENCIL_SIZE>(this->mResidualVector, p_indices, b_boundary_elem);
167 }
168 }
169 }
170
171 this->FinishAssembleSystem(assembleResidual, assembleJacobian);
172}
173
174template<size_t DIM>
176 Element<DIM, DIM>& rElement,
177 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElem,
178 c_matrix<double, STENCIL_SIZE, STENCIL_SIZE >& rAElemPrecond,
179 c_vector<double, STENCIL_SIZE>& rBElem,
180 bool assembleResidual,
181 bool assembleJacobian)
182{
183 static c_matrix<double,DIM,DIM> jacobian;
184 static c_matrix<double,DIM,DIM> inverse_jacobian;
185 double jacobian_determinant;
186
187 this->mrQuadMesh.GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
188
189 if (assembleJacobian)
190 {
191 rAElem.clear();
192 rAElemPrecond.clear();
193 }
194
195 if (assembleResidual)
196 {
197 rBElem.clear();
198 }
199
200 // Get the current displacement at the nodes
201 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
202 for (unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
203 {
204 for (unsigned JJ=0; JJ<DIM; JJ++)
205 {
206 element_current_displacements(JJ,II) = this->mCurrentSolution[DIM*rElement.GetNodeGlobalIndex(II) + JJ];
207 }
208 }
209
210 // Allocate memory for the basis functions values and derivative values
211 static c_vector<double, NUM_VERTICES_PER_ELEMENT> linear_phi;
212 static c_vector<double, NUM_NODES_PER_ELEMENT> quad_phi;
213 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
214 static c_matrix<double, NUM_NODES_PER_ELEMENT, DIM> trans_grad_quad_phi;
215
216 // Get the material law
218 = this->mrProblemDefinition.GetCompressibleMaterialLaw(rElement.GetIndex());
219
220
221 static c_matrix<double,DIM,DIM> grad_u; // grad_u = (du_i/dX_M)
222
223 static c_matrix<double,DIM,DIM> F; // the deformation gradient, F = dx/dX, F_{iM} = dx_i/dX_M
224 static c_matrix<double,DIM,DIM> C; // Green deformation tensor, C = F^T F
225 static c_matrix<double,DIM,DIM> inv_C; // inverse(C)
226 static c_matrix<double,DIM,DIM> inv_F; // inverse(F)
227 static c_matrix<double,DIM,DIM> T; // Second Piola-Kirchoff stress tensor (= dW/dE = 2dW/dC)
228
229 static c_matrix<double,DIM,DIM> F_T; // F*T
230 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> F_T_grad_quad_phi; // F*T*grad_quad_phi
231
232 c_vector<double,DIM> body_force;
233
234 static FourthOrderTensor<DIM,DIM,DIM,DIM> dTdE; // dTdE(M,N,P,Q) = dT_{MN}/dE_{PQ}
235 static FourthOrderTensor<DIM,DIM,DIM,DIM> dSdF; // dSdF(M,i,N,j) = dS_{Mi}/dF_{jN}
236
239
240 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> temp_matrix;
241 static c_matrix<double,NUM_NODES_PER_ELEMENT,DIM> grad_quad_phi_times_invF;
242
243 if (this->mSetComputeAverageStressPerElement)
244 {
245 this->mAverageStressesPerElement[rElement.GetIndex()] = zero_vector<double>(DIM*(DIM+1)/2);
246 }
247
248 // Loop over Gauss points
249 for (unsigned quadrature_index=0; quadrature_index < this->mpQuadratureRule->GetNumQuadPoints(); quadrature_index++)
250 {
251 // This is needed by the cardiac mechanics solver
252 unsigned current_quad_point_global_index = rElement.GetIndex()*this->mpQuadratureRule->GetNumQuadPoints()
253 + quadrature_index;
254
255 double wJ = jacobian_determinant * this->mpQuadratureRule->GetWeight(quadrature_index);
256
257 const ChastePoint<DIM>& quadrature_point = this->mpQuadratureRule->rGetQuadPoint(quadrature_index);
258
259 // Set up basis function information
260 LinearBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, linear_phi);
261 QuadraticBasisFunction<DIM>::ComputeBasisFunctions(quadrature_point, quad_phi);
262 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
263 trans_grad_quad_phi = trans(grad_quad_phi);
264
265 // Get the body force, interpolating X if necessary
266 if (assembleResidual)
267 {
268 switch (this->mrProblemDefinition.GetBodyForceType())
269 {
270 case FUNCTIONAL_BODY_FORCE:
271 {
272 c_vector<double,DIM> X = zero_vector<double>(DIM);
273 // interpolate X (using the vertices and the /linear/ bases, as no curvilinear elements
274 for (unsigned node_index=0; node_index<NUM_VERTICES_PER_ELEMENT; node_index++)
275 {
276 X += linear_phi(node_index)*this->mrQuadMesh.GetNode( rElement.GetNodeGlobalIndex(node_index) )->rGetLocation();
277 }
278 body_force = this->mrProblemDefinition.EvaluateBodyForceFunction(X, this->mCurrentTime);
279 break;
280 }
281 case CONSTANT_BODY_FORCE:
282 {
283 body_force = this->mrProblemDefinition.GetConstantBodyForce();
284 break;
285 }
286 default:
288 }
289 }
290
291 // Interpolate grad_u
292 grad_u = zero_matrix<double>(DIM,DIM);
293 for (unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
294 {
295 for (unsigned i=0; i<DIM; i++)
296 {
297 for (unsigned M=0; M<DIM; M++)
298 {
299 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
300 }
301 }
302 }
303
304 // Calculate C, inv(C) and T
305 for (unsigned i=0; i<DIM; i++)
306 {
307 for (unsigned M=0; M<DIM; M++)
308 {
309 F(i,M) = (i==M?1:0) + grad_u(i,M);
310 }
311 }
312
313 C = prod(trans(F),F);
314 inv_C = Inverse(C);
315 inv_F = Inverse(F);
316
317 // Compute the passive stress, and dTdE corresponding to passive stress
318 this->SetupChangeOfBasisMatrix(rElement.GetIndex(), current_quad_point_global_index);
319 p_material_law->SetChangeOfBasisMatrix(this->mChangeOfBasisMatrix);
320 p_material_law->ComputeStressAndStressDerivative(C, inv_C, 0.0, T, dTdE, assembleJacobian);
321
322 if (this->mIncludeActiveTension)
323 {
324 // Add any active stresses, if there are any. Requires subclasses to overload this method,
325 // see for example the cardiac mechanics assemblers.
326 this->AddActiveStressAndStressDerivative(C, rElement.GetIndex(), current_quad_point_global_index,
327 T, dTdE, assembleJacobian);
328 }
329
330 if (this->mSetComputeAverageStressPerElement)
331 {
332 this->AddStressToAverageStressPerElement(T,rElement.GetIndex());
333 }
334
335 // Residual vector
336 if (assembleResidual)
337 {
338 F_T = prod(F,T);
339 F_T_grad_quad_phi = prod(F_T, grad_quad_phi);
340
341 for (unsigned index=0; index<NUM_NODES_PER_ELEMENT*DIM; index++)
342 {
343 unsigned spatial_dim = index%DIM;
344 unsigned node_index = (index-spatial_dim)/DIM;
345
346 rBElem(index) += - this->mrProblemDefinition.GetDensity()
347 * body_force(spatial_dim)
348 * quad_phi(node_index)
349 * wJ;
350
351 // The T(M,N)*F(spatial_dim,M)*grad_quad_phi(N,node_index) term
352 rBElem(index) += F_T_grad_quad_phi(spatial_dim,node_index)
353 * wJ;
354 }
355 }
356
357 // Jacobian matrix
358 if (assembleJacobian)
359 {
360 // Save trans(grad_quad_phi) * invF
361 grad_quad_phi_times_invF = prod(trans_grad_quad_phi, inv_F);
362
364 // Set up the tensor dSdF
365 //
366 // dSdF as a function of T and dTdE (which is what the material law returns) is given by:
367 //
368 // dS_{Mi}/dF_{jN} = (dT_{MN}/dC_{PQ}+dT_{MN}/dC_{PQ}) F{iP} F_{jQ} + T_{MN} delta_{ij}
369 //
370 // todo1: this should probably move into the material law (but need to make sure
371 // memory is handled efficiently
372 // todo2: get material law to return this immediately, not dTdE
374
375 // Set up the tensor 0.5(dTdE(M,N,P,Q) + dTdE(M,N,Q,P))
376 for (unsigned M=0; M<DIM; M++)
377 {
378 for (unsigned N=0; N<DIM; N++)
379 {
380 for (unsigned P=0; P<DIM; P++)
381 {
382 for (unsigned Q=0; Q<DIM; Q++)
383 {
384 // this is NOT dSdF, just using this as storage space
385 dSdF(M,N,P,Q) = 0.5*(dTdE(M,N,P,Q) + dTdE(M,N,Q,P));
386 }
387 }
388 }
389 }
390
391 // This is NOT dTdE, just reusing memory. A^{MdPQ} = F^d_N * dTdE_sym^{MNPQ}
392 dTdE.template SetAsContractionOnSecondDimension<DIM>(F, dSdF);
393
394 // dSdF{MdPe} := F^d_N * F^e_Q * dTdE_sym^{MNPQ}
395 dSdF.template SetAsContractionOnFourthDimension<DIM>(F, dTdE);
396
397 // Now add the T_{MN} delta_{ij} term
398 for (unsigned M=0; M<DIM; M++)
399 {
400 for (unsigned N=0; N<DIM; N++)
401 {
402 for (unsigned i=0; i<DIM; i++)
403 {
404 dSdF(M,i,N,i) += T(M,N);
405 }
406 }
407 }
408
410 // Set up the tensor
411 // dSdF_quad_quad(node_index1, spatial_dim1, node_index2, spatial_dim2)
412 // = dS_{M,spatial_dim1}/d_F{spatial_dim2,N}
413 // * grad_quad_phi(M,node_index1)
414 // * grad_quad_phi(P,node_index2)
415 //
416 // = dSdF(M,spatial_index1,N,spatial_index2)
417 // * grad_quad_phi(M,node_index1)
418 // * grad_quad_phi(P,node_index2)
419 //
421 temp_tensor.template SetAsContractionOnFirstDimension<DIM>(trans_grad_quad_phi, dSdF);
422 dSdF_quad_quad.template SetAsContractionOnThirdDimension<DIM>(trans_grad_quad_phi, temp_tensor);
423
424 for (unsigned index1=0; index1<NUM_NODES_PER_ELEMENT*DIM; index1++)
425 {
426 unsigned spatial_dim1 = index1%DIM;
427 unsigned node_index1 = (index1-spatial_dim1)/DIM;
428
429 for (unsigned index2=0; index2<NUM_NODES_PER_ELEMENT*DIM; index2++)
430 {
431 unsigned spatial_dim2 = index2%DIM;
432 unsigned node_index2 = (index2-spatial_dim2)/DIM;
433
434 // The dSdF*grad_quad_phi*grad_quad_phi term
435 rAElem(index1,index2) += dSdF_quad_quad(node_index1,spatial_dim1,node_index2,spatial_dim2)
436 * wJ;
437 }
438 }
439 }
440 }
441
442 rAElemPrecond.clear();
443 if (assembleJacobian)
444 {
445 rAElemPrecond = rAElem;
446 }
447
448 if (this->mSetComputeAverageStressPerElement)
449 {
450 for (unsigned i=0; i<DIM*(DIM+1)/2; i++)
451 {
452 this->mAverageStressesPerElement[rElement.GetIndex()](i) /= this->mpQuadratureRule->GetNumQuadPoints();
453 }
454 }
455}
456
457template<size_t DIM>
459 SolidMechanicsProblemDefinition<DIM>& rProblemDefinition,
460 std::string outputDirectory)
461 : AbstractNonlinearElasticitySolver<DIM>(rQuadMesh,
462 rProblemDefinition,
463 outputDirectory,
464 COMPRESSIBLE)
465{
466 if (rProblemDefinition.GetCompressibilityType() != COMPRESSIBLE)
467 {
468 EXCEPTION("SolidMechanicsProblemDefinition object contains incompressible material laws");
469 }
470}
471
472template<size_t DIM>
476
477// Explicit instantiation
#define EXCEPTION(message)
#define NEVER_REACHED
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()
CompressibleNonlinearElasticitySolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)
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)
void AssembleSystem(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 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)