AbstractIsotropicCompressibleMaterialLaw.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "AbstractIsotropicCompressibleMaterialLaw.hpp"
00030 
00031 #include "Debug.hpp"
00032 
00033 template<unsigned DIM>
00034 AbstractIsotropicCompressibleMaterialLaw<DIM>::~AbstractIsotropicCompressibleMaterialLaw()
00035 {
00036 }
00037 
00038 template<unsigned DIM>
00039 void AbstractIsotropicCompressibleMaterialLaw<DIM>::ComputeStressAndStressDerivative(c_matrix<double,DIM,DIM>& rC,
00040                                                                                      c_matrix<double,DIM,DIM>& rInvC,
00041                                                                                      double                    pressure,
00042                                                                                      c_matrix<double,DIM,DIM>& rT,
00043                                                                                      FourthOrderTensor<DIM,DIM,DIM,DIM>&   rDTdE,
00044                                                                                      bool                      computeDTdE)
00045 {
00046     // this is covered, but gcov doesn't see this as being covered
00047     // for some reason, maybe because of optimisations
00048     #define COVERAGE_IGNORE
00049     assert((DIM==2) || (DIM==3));
00050     #undef COVERAGE_IGNORE
00051 
00052     assert(pressure==0.0);
00053 
00054     static c_matrix<double,DIM,DIM> identity = identity_matrix<double>(DIM);
00055 
00056     double I1 = Trace(rC);
00057     double I2 = SecondInvariant(rC);
00058     double I3 = Determinant(rC);
00059 
00060     static c_matrix<double,DIM,DIM> dI2dC;
00061     dI2dC = I1*identity - rC;              // MUST be on separate line to above!
00062 
00063     double w1 = Get_dW_dI1(I1,I2,I3);
00064     double w2 = Get_dW_dI2(I1,I2,I3);
00065     double w3 = Get_dW_dI3(I1,I2,I3);
00066 
00067 
00068     // Compute stress:  **** See FiniteElementImplementations document. ****
00069     //
00070     //  T = dW_dE
00071     //    = 2 dW_dC
00072     //    = 2 (  w1 dI1/dC   +  w2 dI2/dC      +   w3 dI3/dC )
00073     //    = 2 (  w1 I        +  w2 (I1*I - C)  +   w3 I3 inv(C) )
00074     //
00075     //  where w1 = dW/dI1, etc
00076     //
00077     rT = 2*w1*identity + 2*w3*I3*rInvC;
00078     if (DIM==3)
00079     {
00080         rT += 2*w2*dI2dC;
00081     }
00082 
00083     // Compute stress derivative if required:  **** See FiniteElementImplementations document. ****
00084     //
00085     // The stress derivative dT_{MN}/dE_{PQ} is
00086     //
00087     //
00088     //  dT_dE = 2 dT_dC
00089     //        = 4  d/dC ( w1 I  +  w2 (I1*I - C)  +   w3 I3 inv(C) )
00090     //  so (in the following ** represents outer product):
00091     //  (1/4) dT_dE =        w11 I**I          +    w12 I**(I1*I-C)           +     w13 I**inv(C)
00092     //                  +    w21 (I1*I-C)**I   +    w22 (I1*I-C)**(I1*I-C)    +     w23 (I1*I-C)**inv(C)           +   w2 (I**I - dC/dC)
00093     //                  +    w31 I3 inv(C)**I  +    w32 I3 inv(C)**(I1*I-C)   +  (w33 I3 + w3) inv(C)**inv(C)      +   w3 d(invC)/dC
00094     //
00095     //  Here, I**I represents the tensor A[M][N][P][Q] = (M==N)*(P==Q) // ie delta(M,N)delta(P,Q),   etc
00096     //
00097 
00098     if (computeDTdE)
00099     {
00100         double  w11    = Get_d2W_dI1(I1,I2,I3);
00101         double  w22    = Get_d2W_dI2(I1,I2,I3);
00102         double  w33    = Get_d2W_dI3(I1,I2,I3);
00103 
00104         double  w23  = Get_d2W_dI2I3(I1,I2,I3);
00105         double  w13  = Get_d2W_dI1I3(I1,I2,I3);
00106         double  w12  = Get_d2W_dI1I2(I1,I2,I3);
00107 
00108         for (unsigned M=0; M<DIM; M++)
00109         {
00110             for (unsigned N=0; N<DIM; N++)
00111             {
00112                 for (unsigned P=0; P<DIM; P++)
00113                 {
00114                     for (unsigned Q=0; Q<DIM; Q++)
00115                     {
00116                         rDTdE(M,N,P,Q) =   4 * w11  * (M==N) * (P==Q)
00117                                          + 4 * w13  * I3 * ( (M==N) * rInvC(P,Q)  +  rInvC(M,N)*(P==Q) )  // the w13 and w31 terms
00118                                          + 4 * (w33*I3 + w3) * I3 * rInvC(M,N) * rInvC(P,Q)
00119                                          - 4 * w3 * I3 * rInvC(M,P) * rInvC(Q,N);
00120 
00121                         if (DIM==3)
00122                         {
00123                             rDTdE(M,N,P,Q) +=   4 * w22  * dI2dC(M,N) * dI2dC(P,Q)
00124                                               + 4 * w12  * ((M==N)*dI2dC(P,Q) + (P==Q)*dI2dC(M,N))          // the w12 and w21 terms
00125                                               + 4 * w23 * I3 * ( dI2dC(M,N)*rInvC(P,Q) + rInvC(M,N)*dI2dC(P,Q)) // the w23 and w32 terms
00126                                               + 4 * w2   * ((M==N)*(P==Q) - (M==P)*(N==Q));
00127                         }
00128                     }
00129                 }
00130             }
00131         }
00132     }
00133 }
00134 
00135 
00136 
00138 // Explicit instantiation
00140 
00141 
00142 //template class AbstractIsotropicCompressibleMaterialLaw<1>;
00143 template class AbstractIsotropicCompressibleMaterialLaw<2>;
00144 template class AbstractIsotropicCompressibleMaterialLaw<3>;

Generated on Tue May 31 14:31:48 2011 for Chaste by  doxygen 1.5.5