Chaste  Release::2017.1
QuadraticBasisFunction.cpp
1 /*
2 
3 Copyright (c) 2005-2017, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
49 double QuadraticBasisFunction<0>::ComputeBasisFunction(const ChastePoint<0>& rPoint, unsigned basisIndex)
50 {
51  assert(basisIndex == 0);
52  return 1.0;
53 }
54 
62 void QuadraticBasisFunction<0>::ComputeBasisFunctions(const ChastePoint<0>& rPoint,
63  c_vector<double, 1>& rReturnValue)
64 {
65  rReturnValue(0) = ComputeBasisFunction(rPoint, 0);
66 }
67 
77 template <unsigned ELEMENT_DIM>
78 double 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 
189 template <unsigned ELEMENT_DIM>
190 c_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 
321 template <unsigned ELEMENT_DIM>
322 void 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 
342 template <unsigned ELEMENT_DIM>
343 void 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 
348  for (unsigned j=0; j<(ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2; j++)
349  {
350  matrix_column<c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2> > column(rReturnValue, j);
351  column = ComputeBasisFunctionDerivative(rPoint, j);
352  }
353 }
354 
368 template <unsigned ELEMENT_DIM>
369 void QuadraticBasisFunction<ELEMENT_DIM>::ComputeTransformedBasisFunctionDerivatives(const ChastePoint<ELEMENT_DIM>& rPoint,
370  const c_matrix<double, ELEMENT_DIM, ELEMENT_DIM>& rInverseJacobian,
371  c_matrix<double, ELEMENT_DIM, (ELEMENT_DIM+1)*(ELEMENT_DIM+2)/2>& rReturnValue)
372 {
373  assert(ELEMENT_DIM < 4 && ELEMENT_DIM > 0);
374 
375  ComputeBasisFunctionDerivatives(rPoint, rReturnValue);
376  rReturnValue = prod(trans(rInverseJacobian), rReturnValue);
377 }
378 
379 // Explicit instantiation
380 template class QuadraticBasisFunction<1>;
381 template class QuadraticBasisFunction<2>;
382 template class QuadraticBasisFunction<3>;
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
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:369
NEVER_REACHED
#define NEVER_REACHED
Definition: Exception.hpp:206
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::ComputeBasisFunctionDerivative
static c_vector< double, ELEMENT_DIM > ComputeBasisFunctionDerivative(const ChastePoint< ELEMENT_DIM > &rPoint, unsigned basisIndex)
Definition: QuadraticBasisFunction.cpp:190
QuadraticBasisFunction::ComputeBasisFunction
static double ComputeBasisFunction(const ChastePoint< ELEMENT_DIM > &rPoint, unsigned basisIndex)
Definition: QuadraticBasisFunction.cpp:78