QuadraticBasisFunction.cpp

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

Generated by  doxygen 1.6.2