- Generated on Mon Nov 5 2012 12:47:33 for Chaste by
1.7.4
|
Chaste Release::3.1
|
00001 /* 00002 00003 Copyright (c) 2005-2012, University of Oxford. 00004 All rights reserved. 00005 00006 University of Oxford means the Chancellor, Masters and Scholars of the 00007 University of Oxford, having an administrative office at Wellington 00008 Square, Oxford OX1 2JD, UK. 00009 00010 This file is part of Chaste. 00011 00012 Redistribution and use in source and binary forms, with or without 00013 modification, are permitted provided that the following conditions are met: 00014 * Redistributions of source code must retain the above copyright notice, 00015 this list of conditions and the following disclaimer. 00016 * Redistributions in binary form must reproduce the above copyright notice, 00017 this list of conditions and the following disclaimer in the documentation 00018 and/or other materials provided with the distribution. 00019 * Neither the name of the University of Oxford nor the names of its 00020 contributors may be used to endorse or promote products derived from this 00021 software without specific prior written permission. 00022 00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE 00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00033 00034 */ 00035 00036 #include <boost/numeric/ublas/matrix.hpp> 00037 #include <boost/numeric/ublas/matrix_proxy.hpp> 00038 00039 #include "QuadraticBasisFunction.hpp" 00040 #include "Exception.hpp" 00041 00051 double QuadraticBasisFunction<0>::ComputeBasisFunction(const ChastePoint<0>& rPoint, unsigned basisIndex) 00052 { 00053 assert(basisIndex == 0); 00054 return 1.0; 00055 } 00056 00064 void QuadraticBasisFunction<0>::ComputeBasisFunctions(const ChastePoint<0>& rPoint, 00065 c_vector<double, 1>& rReturnValue) 00066 { 00067 rReturnValue(0) = ComputeBasisFunction(rPoint, 0); 00068 } 00069 00079 template <unsigned ELEMENT_DIM> 00080 double QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunction(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex) 00081 { 00082 assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0); 00083 double x, y, z; 00084 switch (ELEMENT_DIM) 00085 { 00086 case 0: 00087 assert(basisIndex == 0); 00088 return 1.0; 00089 break; 00090 00091 case 1: 00092 x = rPoint[0]; 00093 switch (basisIndex) 00094 { 00095 case 0: 00096 return 2.0*(x-1.0)*(x-0.5); 00097 break; 00098 case 1: 00099 return 2.0*x*(x-0.5); 00100 break; 00101 case 2: 00102 return 4.0*x*(1.0-x); 00103 break; 00104 default: 00105 NEVER_REACHED; 00106 } 00107 break; 00108 00109 case 2: 00110 x = rPoint[0]; 00111 y = rPoint[1]; 00112 switch (basisIndex) 00113 { 00114 case 0: // the node at (0,0) 00115 return 2.0 * (1.0 - x - y) * (0.5 - x - y); 00116 break; 00117 case 1: // the node at (1,0) 00118 return 2.0*x*(x-0.5); 00119 break; 00120 case 2: // the node at (0,1) 00121 return 2.0*y*(y-0.5); 00122 break; 00123 case 3: // the node opposite 0, which is (1/2,1/2) 00124 return 4.0 * y * x; 00125 break; 00126 case 4: // the node opposite 1, which is (0,1/2) 00127 return 4.0 * (1.0 - x - y) * y; 00128 break; 00129 case 5: // the node opposite 2, which is (1/2,0) 00130 return 4.0 * (1.0 - x - y) * x; 00131 break; 00132 default: 00133 NEVER_REACHED; 00134 } 00135 break; 00136 00137 case 3: 00138 x = rPoint[0]; 00139 y = rPoint[1]; 00140 z = rPoint[2]; 00141 switch (basisIndex) 00142 { 00143 case 0: // the node at (0,0,0) 00144 return 2.0 * (1.0 - x - y - z) * (0.5 - x - y - z); 00145 break; 00146 case 1: // the node at (1,0,0) 00147 return 2.0*x*(x-0.5); 00148 break; 00149 case 2: // the node at (0,1,0) 00150 return 2.0*y*(y-0.5); 00151 break; 00152 case 3: // the node at (0,0,1) 00153 return 2.0*z*(z-0.5); 00154 break; 00155 case 4: // our (tetgen convention), node4 is between nodes 0 and 1, (1/2,0,0) 00156 return 4.0 * (1.0 - x - y - z) * x; 00157 break; 00158 case 5: // our (tetgen convention), node5 is between nodes 1 and 2, (1/2,1/2,0) 00159 return 4 * x * y; 00160 break; 00161 case 6: // our (tetgen convention), node6 is between nodes 0 and 2, (0,1/2,0) 00162 return 4.0 * (1.0 - x - y - z) * y; 00163 break; 00164 case 7: // our (tetgen convention), node7 is between nodes 0 and 3, (0,0,1/2) 00165 return 4.0 * (1.0 - x - y - z) * z; 00166 break; 00167 case 8: // our (tetgen convention), node8 is between nodes 1 and 3, (1/2,0,1/2) 00168 return 4.0 * x * z; 00169 break; 00170 case 9: // our (tetgen convention), node9 is between nodes 2 and 3, (0,1/2,1/2) 00171 return 4.0 * y * z; 00172 break; 00173 default: 00174 NEVER_REACHED; 00175 } 00176 break; 00177 } 00178 return 0.0; // Avoid compiler warning 00179 } 00180 00191 template <unsigned ELEMENT_DIM> 00192 c_vector<double, ELEMENT_DIM> QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivative(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex) 00193 { 00194 c_vector<double, ELEMENT_DIM> gradN; 00195 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0); 00196 00197 double x, y, z; 00198 switch (ELEMENT_DIM) 00199 { 00200 case 1: 00201 x = rPoint[0]; 00202 switch (basisIndex) 00203 { 00204 case 0: 00205 gradN(0) = 4.0*x-3.0; 00206 break; 00207 case 1: 00208 gradN(0) = 4.0*x-1.0; 00209 break; 00210 case 2: 00211 gradN(0) = 4.0-8.0*x; 00212 break; 00213 default: 00214 NEVER_REACHED; 00215 } 00216 break; 00217 00218 case 2: 00219 x = rPoint[0]; 00220 y = rPoint[1]; 00221 switch (basisIndex) 00222 { 00223 case 0: 00224 gradN(0) = -3.0 + 4.0*x + 4.0*y; 00225 gradN(1) = -3.0 + 4.0*x + 4.0*y; 00226 break; 00227 case 1: 00228 gradN(0) = 4.0*x - 1.0; 00229 gradN(1) = 0.0; 00230 break; 00231 case 2: 00232 gradN(0) = 0.0; 00233 gradN(1) = 4.0*y - 1.0; 00234 break; 00235 case 3: 00236 gradN(0) = 4.0*y; 00237 gradN(1) = 4.0*x; 00238 break; 00239 case 4: 00240 gradN(0) = -4.0*y; 00241 gradN(1) = 4.0-4.0*x-8.0*y; 00242 break; 00243 case 5: 00244 gradN(0) = 4.0-8.0*x-4.0*y; 00245 gradN(1) = -4.0*x; 00246 break; 00247 default: 00248 NEVER_REACHED; 00249 } 00250 break; 00251 00252 case 3: 00253 x = rPoint[0]; 00254 y = rPoint[1]; 00255 z = rPoint[2]; 00256 switch (basisIndex) 00257 { 00258 case 0: 00259 gradN(0) = -3.0 + 4.0*(x+y+z); 00260 gradN(1) = -3.0 + 4.0*(x+y+z); 00261 gradN(2) = -3.0 + 4.0*(x+y+z); 00262 break; 00263 case 1: 00264 gradN(0) = 4.0*x-1.0; 00265 gradN(1) = 0; 00266 gradN(2) = 0; 00267 break; 00268 case 2: 00269 gradN(0) = 0; 00270 gradN(1) = 4.0*y-1.0; 00271 gradN(2) = 0; 00272 break; 00273 case 3: 00274 gradN(0) = 0; 00275 gradN(1) = 0; 00276 gradN(2) = 4.0*z-1.0; 00277 break; 00278 case 4: 00279 gradN(0) = 4.0-8.0*x-4.0*y-4.0*z; 00280 gradN(1) = -4.0*x; 00281 gradN(2) = -4.0*x; 00282 break; 00283 case 5: 00284 gradN(0) = 4.0*y; 00285 gradN(1) = 4.0*x; 00286 gradN(2) = 0.0; 00287 break; 00288 case 6: 00289 gradN(0) = -4.0*y; 00290 gradN(1) = 4.0-4.0*x-8.0*y-4.0*z; 00291 gradN(2) = -4.0*y; 00292 break; 00293 case 7: 00294 gradN(0) = -4.0*z; 00295 gradN(1) = -4.0*z; 00296 gradN(2) = 4.0-4.0*x-4.0*y-8.0*z; 00297 break; 00298 case 8: 00299 gradN(0) = 4.0*z; 00300 gradN(1) = 0; 00301 gradN(2) = 4.0*x; 00302 break; 00303 case 9: 00304 gradN(0) = 0; 00305 gradN(1) = 4.0*z; 00306 gradN(2) = 4.0*y; 00307 break; 00308 default: 00309 NEVER_REACHED; 00310 } 00311 break; 00312 } 00313 return gradN; 00314 } 00315 00323 template <unsigned ELEMENT_DIM> 00324 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctions(const ChastePoint<ELEMENT_DIM>& rPoint, 00325 c_vector<double, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue) 00326 { 00327 assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0); 00328 00329 for (unsigned i=0; i<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; i++) 00330 { 00331 rReturnValue(i) = ComputeBasisFunction(rPoint, i); 00332 } 00333 } 00334 00344 template <unsigned ELEMENT_DIM> 00345 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint, 00346 c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue) 00347 { 00348 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0); 00349 00350 for (unsigned j=0; j<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; j++) 00351 { 00352 matrix_column<c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2> > column(rReturnValue, j); 00353 column = ComputeBasisFunctionDerivative(rPoint, j); 00354 } 00355 } 00356 00370 template <unsigned ELEMENT_DIM> 00371 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeTransformedBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint, 00372 const c_matrix<double, ELEMENT_DIM, ELEMENT_DIM>& rInverseJacobian, 00373 c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue) 00374 { 00375 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0); 00376 00377 ComputeBasisFunctionDerivatives(rPoint, rReturnValue); 00378 rReturnValue = prod(trans(rInverseJacobian), rReturnValue); 00379 } 00380 00382 // Explicit instantiation 00384 00385 template class QuadraticBasisFunction<1>; 00386 template class QuadraticBasisFunction<2>; 00387 template class QuadraticBasisFunction<3>;
1.7.4