Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
QuadraticBasisFunction.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "QuadraticBasisFunction.hpp"
37#include "Exception.hpp"
38#include "UblasIncludes.hpp"
39
49double QuadraticBasisFunction<0>::ComputeBasisFunction(const ChastePoint<0>& rPoint, unsigned basisIndex)
50{
51 assert(basisIndex == 0);
52 return 1.0;
53}
54
62void QuadraticBasisFunction<0>::ComputeBasisFunctions(const ChastePoint<0>& rPoint,
63 c_vector<double, 1>& rReturnValue)
64{
65 rReturnValue(0) = ComputeBasisFunction(rPoint, 0);
66}
67
77template <unsigned ELEMENT_DIM>
78double QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunction(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex)
79{
80 assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0);
81 double x, y, z;
82 switch (ELEMENT_DIM)
83 {
84 case 0:
85 assert(basisIndex == 0);
86 return 1.0;
87 break;
88
89 case 1:
90 x = rPoint[0];
91 switch (basisIndex)
92 {
93 case 0:
94 return 2.0*(x-1.0)*(x-0.5);
95 break;
96 case 1:
97 return 2.0*x*(x-0.5);
98 break;
99 case 2:
100 return 4.0*x*(1.0-x);
101 break;
102 default:
103 NEVER_REACHED;
104 }
105 break;
106
107 case 2:
108 x = rPoint[0];
109 y = rPoint[1];
110 switch (basisIndex)
111 {
112 case 0: // the node at (0,0)
113 return 2.0 * (1.0 - x - y) * (0.5 - x - y);
114 break;
115 case 1: // the node at (1,0)
116 return 2.0*x*(x-0.5);
117 break;
118 case 2: // the node at (0,1)
119 return 2.0*y*(y-0.5);
120 break;
121 case 3: // the node opposite 0, which is (1/2,1/2)
122 return 4.0 * y * x;
123 break;
124 case 4: // the node opposite 1, which is (0,1/2)
125 return 4.0 * (1.0 - x - y) * y;
126 break;
127 case 5: // the node opposite 2, which is (1/2,0)
128 return 4.0 * (1.0 - x - y) * x;
129 break;
130 default:
131 NEVER_REACHED;
132 }
133 break;
134
135 case 3:
136 x = rPoint[0];
137 y = rPoint[1];
138 z = rPoint[2];
139 switch (basisIndex)
140 {
141 case 0: // the node at (0,0,0)
142 return 2.0 * (1.0 - x - y - z) * (0.5 - x - y - z);
143 break;
144 case 1: // the node at (1,0,0)
145 return 2.0*x*(x-0.5);
146 break;
147 case 2: // the node at (0,1,0)
148 return 2.0*y*(y-0.5);
149 break;
150 case 3: // the node at (0,0,1)
151 return 2.0*z*(z-0.5);
152 break;
153 case 4: // our (tetgen convention), node4 is between nodes 0 and 1, (1/2,0,0)
154 return 4.0 * (1.0 - x - y - z) * x;
155 break;
156 case 5: // our (tetgen convention), node5 is between nodes 1 and 2, (1/2,1/2,0)
157 return 4 * x * y;
158 break;
159 case 6: // our (tetgen convention), node6 is between nodes 0 and 2, (0,1/2,0)
160 return 4.0 * (1.0 - x - y - z) * y;
161 break;
162 case 7: // our (tetgen convention), node7 is between nodes 0 and 3, (0,0,1/2)
163 return 4.0 * (1.0 - x - y - z) * z;
164 break;
165 case 8: // our (tetgen convention), node8 is between nodes 1 and 3, (1/2,0,1/2)
166 return 4.0 * x * z;
167 break;
168 case 9: // our (tetgen convention), node9 is between nodes 2 and 3, (0,1/2,1/2)
169 return 4.0 * y * z;
170 break;
171 default:
172 NEVER_REACHED;
173 }
174 break;
175 }
176 return 0.0; // Avoid compiler warning
177}
178
189template <unsigned ELEMENT_DIM>
190c_vector<double, ELEMENT_DIM> QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivative(const ChastePoint<ELEMENT_DIM>& rPoint, unsigned basisIndex)
191{
192 c_vector<double, ELEMENT_DIM> gradN;
193 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
194
195 double x, y, z;
196 switch (ELEMENT_DIM)
197 {
198 case 1:
199 x = rPoint[0];
200 switch (basisIndex)
201 {
202 case 0:
203 gradN(0) = 4.0*x-3.0;
204 break;
205 case 1:
206 gradN(0) = 4.0*x-1.0;
207 break;
208 case 2:
209 gradN(0) = 4.0-8.0*x;
210 break;
211 default:
212 NEVER_REACHED;
213 }
214 break;
215
216 case 2:
217 x = rPoint[0];
218 y = rPoint[1];
219 switch (basisIndex)
220 {
221 case 0:
222 gradN(0) = -3.0 + 4.0*x + 4.0*y;
223 gradN(1) = -3.0 + 4.0*x + 4.0*y;
224 break;
225 case 1:
226 gradN(0) = 4.0*x - 1.0;
227 gradN(1) = 0.0;
228 break;
229 case 2:
230 gradN(0) = 0.0;
231 gradN(1) = 4.0*y - 1.0;
232 break;
233 case 3:
234 gradN(0) = 4.0*y;
235 gradN(1) = 4.0*x;
236 break;
237 case 4:
238 gradN(0) = -4.0*y;
239 gradN(1) = 4.0-4.0*x-8.0*y;
240 break;
241 case 5:
242 gradN(0) = 4.0-8.0*x-4.0*y;
243 gradN(1) = -4.0*x;
244 break;
245 default:
246 NEVER_REACHED;
247 }
248 break;
249
250 case 3:
251 x = rPoint[0];
252 y = rPoint[1];
253 z = rPoint[2];
254 switch (basisIndex)
255 {
256 case 0:
257 gradN(0) = -3.0 + 4.0*(x+y+z);
258 gradN(1) = -3.0 + 4.0*(x+y+z);
259 gradN(2) = -3.0 + 4.0*(x+y+z);
260 break;
261 case 1:
262 gradN(0) = 4.0*x-1.0;
263 gradN(1) = 0;
264 gradN(2) = 0;
265 break;
266 case 2:
267 gradN(0) = 0;
268 gradN(1) = 4.0*y-1.0;
269 gradN(2) = 0;
270 break;
271 case 3:
272 gradN(0) = 0;
273 gradN(1) = 0;
274 gradN(2) = 4.0*z-1.0;
275 break;
276 case 4:
277 gradN(0) = 4.0-8.0*x-4.0*y-4.0*z;
278 gradN(1) = -4.0*x;
279 gradN(2) = -4.0*x;
280 break;
281 case 5:
282 gradN(0) = 4.0*y;
283 gradN(1) = 4.0*x;
284 gradN(2) = 0.0;
285 break;
286 case 6:
287 gradN(0) = -4.0*y;
288 gradN(1) = 4.0-4.0*x-8.0*y-4.0*z;
289 gradN(2) = -4.0*y;
290 break;
291 case 7:
292 gradN(0) = -4.0*z;
293 gradN(1) = -4.0*z;
294 gradN(2) = 4.0-4.0*x-4.0*y-8.0*z;
295 break;
296 case 8:
297 gradN(0) = 4.0*z;
298 gradN(1) = 0;
299 gradN(2) = 4.0*x;
300 break;
301 case 9:
302 gradN(0) = 0;
303 gradN(1) = 4.0*z;
304 gradN(2) = 4.0*y;
305 break;
306 default:
307 NEVER_REACHED;
308 }
309 break;
310 }
311 return gradN;
312}
313
321template <unsigned ELEMENT_DIM>
322void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctions(const ChastePoint<ELEMENT_DIM>& rPoint,
323 c_vector<double, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
324{
325 assert(ELEMENT_DIM < 4 && ELEMENT_DIM >= 0);
326
327 for (unsigned i=0; i<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; i++)
328 {
329 rReturnValue(i) = ComputeBasisFunction(rPoint, i);
330 }
331}
332
342template <unsigned ELEMENT_DIM>
343void QuadraticBasisFunction<ELEMENT_DIM>::ComputeBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
344 c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
345{
346 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
347 for (unsigned j=0; j<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; j++)
348 {
349 matrix_column<c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2> > column(rReturnValue, j);
350 column = ComputeBasisFunctionDerivative(rPoint, j);
351 }
352}
353
367template <unsigned ELEMENT_DIM>
368void QuadraticBasisFunction<ELEMENT_DIM>::ComputeTransformedBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
369 const c_matrix<double, ELEMENT_DIM, ELEMENT_DIM>& rInverseJacobian,
370 c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
371{
372 assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
373
374 ComputeBasisFunctionDerivatives(rPoint, rReturnValue);
375 rReturnValue = prod(trans(rInverseJacobian), rReturnValue);
376}
377
378// Explicit instantiation
379template class QuadraticBasisFunction<1>;
380template class QuadraticBasisFunction<2>;
381template class QuadraticBasisFunction<3>;
NEVER_REACHED
#define NEVER_REACHED
Definition Exception.hpp:227
QuadraticBasisFunction::ComputeBasisFunction
static double ComputeBasisFunction(const ChastePoint< ELEMENT_DIM > &rPoint, unsigned basisIndex)
Definition QuadraticBasisFunction.cpp:78
QuadraticBasisFunction::ComputeBasisFunctionDerivative
static c_vector< double, ELEMENT_DIM > ComputeBasisFunctionDerivative(const ChastePoint< ELEMENT_DIM > &rPoint, unsigned basisIndex)
Definition QuadraticBasisFunction.cpp:190
QuadraticBasisFunction::ComputeBasisFunctions
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)
Definition QuadraticBasisFunction.cpp:322
QuadraticBasisFunction::ComputeTransformedBasisFunctionDerivatives
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)
Definition QuadraticBasisFunction.cpp:368
QuadraticBasisFunction::ComputeBasisFunctionDerivatives
static void ComputeBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, c_matrix< double, ELEMENT_DIM,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)
Definition QuadraticBasisFunction.cpp:343