Chaste  Release::2018.1
CompressibleNonlinearElasticitySolver.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
53 template<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
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 
174 template<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 
457 template<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 
472 template<size_t DIM>
474 {
475 }
476 
477 // Explicit instantiation
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
void AssembleSystem(bool assembleResidual, bool assembleJacobian)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
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)
#define NEVER_REACHED
Definition: Exception.hpp:206
bool OptionExists(const std::string &rOption)
bool GetOwnership() const
static void Zero(Mat matrix)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2 > &rReturnValue)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double, ELEMENT_DIM+1 > &rReturnValue)
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 Zero(Vec vector)
unsigned GetIndex() const
static CommandLineArguments * Instance()
static unsigned GetMyRank()
Definition: PetscTools.cpp:114
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)
CompressibleNonlinearElasticitySolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory)
static void Finalise(Vec vector)