QuadraticBasisFunction.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 <boost/numeric/ublas/matrix.hpp>
00030 #include <boost/numeric/ublas/matrix_proxy.hpp>
00031 
00032 #include "QuadraticBasisFunction.hpp"
00033 #include "PetscTools.hpp"
00034 
00035 
00044 double QuadraticBasisFunction<0>::ComputeBasisFunction(const ChastePoint<0>& rPoint, unsigned basisIndex)
00045 {
00046     assert(basisIndex == 0);
00047     return 1.0;
00048 }
00049 
00056 void QuadraticBasisFunction<0>::ComputeBasisFunctions(const ChastePoint<0>& rPoint,
00057                                                       c_vector<double, 1>& rReturnValue)
00058 {
00059     rReturnValue(0) = ComputeBasisFunction(rPoint, 0);
00060 }
00061 
00071 template <unsigned ELEMENT_DIM>
00072 double QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunction(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex)
00073 {
00074     assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0);
00075     double x, y, z;
00076     switch (ELEMENT_DIM)
00077     {
00078     case 0:
00079         assert(basisIndex == 0);
00080         return 1.0;
00081         break;
00082 
00083     case 1:
00084         x = rPoint[0];
00085         switch (basisIndex)
00086         {
00087             case 0:
00088                 return 2.0*(x-1.0)*(x-0.5);
00089                 break;
00090             case 1:
00091                 return 2.0*x*(x-0.5);
00092                 break;
00093             case 2:
00094                 return 4.0*x*(1.0-x);
00095                 break;
00096             default:
00097                 NEVER_REACHED;
00098         }
00099         break;
00100 
00101     case 2:
00102         x = rPoint[0];
00103         y = rPoint[1];
00104         switch (basisIndex)
00105         {
00106             case 0: // the node at (0,0)
00107                 return 2.0 * (1.0 - x - y) * (0.5 - x - y);
00108                 break;
00109             case 1: // the node at (1,0)
00110                 return 2.0*x*(x-0.5);
00111                 break;
00112             case 2: // the node at (0,1)
00113                 return 2.0*y*(y-0.5);
00114                 break;
00115             case 3: // the node opposite 0, which is (1/2,1/2)
00116                 return 4.0 * y * x;
00117                 break;
00118             case 4: // the node opposite 1, which is (0,1/2)
00119                 return 4.0 * (1.0 - x - y) * y;
00120                 break;
00121             case 5: // the node opposite 2, which is (1/2,0)
00122                 return 4.0 * (1.0 - x - y) * x;
00123                 break;
00124             default:
00125                 NEVER_REACHED;
00126         }
00127         break;
00128 
00129     case 3:
00130         x = rPoint[0];
00131         y = rPoint[1];
00132         z = rPoint[2];
00133         switch (basisIndex)
00134         {
00135             case 0: // the node at (0,0,0)
00136                 return 2.0 * (1.0 - x - y - z) * (0.5 - x - y - z);
00137                 break;
00138             case 1: // the node at (1,0,0)
00139                 return 2.0*x*(x-0.5);
00140                 break;
00141             case 2: // the node at (0,1,0)
00142                 return 2.0*y*(y-0.5);
00143                 break;
00144             case 3: // the node at (0,0,1)
00145                 return 2.0*z*(z-0.5);
00146                 break;
00147             case 4: // our (tetgen convention), node4 is between nodes 0 and 1, (1/2,0,0)
00148                 return 4.0 * (1.0 - x - y - z) * x;
00149                 break;
00150             case 5: // our (tetgen convention), node5 is between nodes 1 and 2, (1/2,1/2,0)
00151                 return 4 * x * y;
00152                 break;
00153             case 6: // our (tetgen convention), node6 is between nodes 0 and 2, (0,1/2,0)
00154                 return 4.0 * (1.0 - x - y - z) * y;
00155                 break;
00156             case 7: // our (tetgen convention), node7 is between nodes 0 and 3, (0,0,1/2)
00157                 return 4.0 * (1.0 - x - y - z) * z;
00158                 break;
00159             case 8: // our (tetgen convention), node8 is between nodes 1 and 3, (1/2,0,1/2)
00160                 return 4.0 * x * z;
00161                 break;
00162             case 9: // our (tetgen convention), node9 is between nodes 2 and 3, (0,1/2,1/2)
00163                 return 4.0 * y * z;
00164                 break;
00165             default:
00166                 NEVER_REACHED;
00167         }
00168         break;
00169     }
00170     return 0.0; // Avoid compiler warning
00171 }
00172 
00183 template <unsigned ELEMENT_DIM>
00184 c_vector<double, ELEMENT_DIM> QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivative(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex)
00185 {
00186     c_vector<double, ELEMENT_DIM> gradN;
00187     assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
00188 
00189     double x, y, z;
00190     switch (ELEMENT_DIM)
00191     {
00192     case 1:
00193         x = rPoint[0];
00194         switch (basisIndex)
00195         {
00196             case 0:
00197                 gradN(0) = 4.0*x-3.0;
00198                 break;
00199             case 1:
00200                 gradN(0) = 4.0*x-1.0;
00201                 break;
00202             case 2:
00203                 gradN(0) = 4.0-8.0*x;
00204                 break;
00205             default:
00206                 NEVER_REACHED;
00207         }
00208         break;
00209 
00210     case 2:
00211         x = rPoint[0];
00212         y = rPoint[1];
00213         switch (basisIndex)
00214         {
00215             case 0:
00216                 gradN(0) = -3.0 + 4.0*x + 4.0*y;
00217                 gradN(1) = -3.0 + 4.0*x + 4.0*y;
00218                 break;
00219             case 1:
00220                 gradN(0) = 4.0*x - 1.0;
00221                 gradN(1) = 0.0;
00222                 break;
00223             case 2:
00224                 gradN(0) = 0.0;
00225                 gradN(1) = 4.0*y - 1.0;
00226                 break;
00227             case 3:
00228                 gradN(0) = 4.0*y;
00229                 gradN(1) = 4.0*x;
00230                 break;
00231             case 4:
00232                 gradN(0) = -4.0*y;
00233                 gradN(1) = 4.0-4.0*x-8.0*y;
00234                 break;
00235             case 5:
00236                 gradN(0) = 4.0-8.0*x-4.0*y;
00237                 gradN(1) = -4.0*x;
00238                 break;
00239             default:
00240                 NEVER_REACHED;
00241         }
00242         break;
00243 
00244     case 3:
00245         x = rPoint[0];
00246         y = rPoint[1];
00247         z = rPoint[2];
00248         switch (basisIndex)
00249         {
00250             case 0:
00251                 gradN(0) = -3.0 + 4.0*(x+y+z);
00252                 gradN(1) = -3.0 + 4.0*(x+y+z);
00253                 gradN(2) = -3.0 + 4.0*(x+y+z);
00254                 break;
00255             case 1:
00256                 gradN(0) =  4.0*x-1.0;
00257                 gradN(1) =  0;
00258                 gradN(2) =  0;
00259                 break;
00260             case 2:
00261                 gradN(0) =  0;
00262                 gradN(1) =  4.0*y-1.0;
00263                 gradN(2) =  0;
00264                 break;
00265             case 3:
00266                 gradN(0) =  0;
00267                 gradN(1) =  0;
00268                 gradN(2) =  4.0*z-1.0;
00269                 break;
00270             case 4:
00271                 gradN(0) =  4.0-8.0*x-4.0*y-4.0*z;
00272                 gradN(1) =  -4.0*x;
00273                 gradN(2) =  -4.0*x;
00274                 break;
00275             case 5:
00276                 gradN(0) =  4.0*y;
00277                 gradN(1) =  4.0*x;
00278                 gradN(2) =  0.0;
00279                 break;
00280             case 6:
00281                 gradN(0) =  -4.0*y;
00282                 gradN(1) =  4.0-4.0*x-8.0*y-4.0*z;
00283                 gradN(2) =  -4.0*y;
00284                 break;
00285             case 7:
00286                 gradN(0) =  -4.0*z;
00287                 gradN(1) =  -4.0*z;
00288                 gradN(2) =  4.0-4.0*x-4.0*y-8.0*z;
00289                 break;
00290             case 8:
00291                 gradN(0) =  4.0*z;
00292                 gradN(1) =  0;
00293                 gradN(2) =  4.0*x;
00294                 break;
00295             case 9:
00296                 gradN(0) =  0;
00297                 gradN(1) =  4.0*z;
00298                 gradN(2) =  4.0*y;
00299                 break;
00300             default:
00301                 NEVER_REACHED;
00302         }
00303         break;
00304     }
00305     return gradN;
00306 }
00307 
00315 template <unsigned ELEMENT_DIM>
00316 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctions(const ChastePoint<ELEMENT_DIM>& rPoint,
00317                                                              c_vector<double, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
00318 {
00319     assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0);
00320 
00321     for (unsigned i=0; i<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; i++)
00322     {
00323         rReturnValue(i) = ComputeBasisFunction(rPoint, i);
00324     }
00325 }
00326 
00336 template <unsigned ELEMENT_DIM>
00337 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
00338                                                                     c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
00339 {
00340     assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
00341 
00342     for (unsigned j=0; j<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; j++)
00343     {
00344         matrix_column<c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2> > column(rReturnValue, j);
00345         column = ComputeBasisFunctionDerivative(rPoint, j);
00346     }
00347 }
00348 
00362 template <unsigned ELEMENT_DIM>
00363 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeTransformedBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
00364                                                                                   const c_matrix<double, ELEMENT_DIM, ELEMENT_DIM>& rInverseJacobian,
00365                                                                                   c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
00366 {
00367     assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
00368 
00369     ComputeBasisFunctionDerivatives(rPoint, rReturnValue);
00370     rReturnValue = prod(trans(rInverseJacobian), rReturnValue);
00371 }
00372 
00373 
00375 // Explicit instantiation
00377 
00378 template class QuadraticBasisFunction<1>;
00379 template class QuadraticBasisFunction<2>;
00380 template class QuadraticBasisFunction<3>;

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