Chaste  Release::3.4
GaussianQuadratureRule.cpp
1 /*
2 
3 Copyright (c) 2005-2016, 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 <cmath>
37 
38 #include "GaussianQuadratureRule.hpp"
39 #include "Exception.hpp"
40 #include "UblasCustomFunctions.hpp"
41 
42 template<unsigned ELEMENT_DIM>
44 {
45  assert(index < mNumQuadPoints);
46  return mPoints[index];
47 }
48 
49 template<unsigned ELEMENT_DIM>
51 {
52  assert(index < mNumQuadPoints);
53  return mWeights[index];
54 }
55 
56 template<unsigned ELEMENT_DIM>
58 {
59  return mNumQuadPoints;
60 }
61 
67 template<>
69 {
70  mNumQuadPoints = 1;
71  mWeights.push_back(1);
72  mPoints.push_back(ChastePoint<0>());
73 }
74 
80 template<>
82 {
83  switch (quadratureOrder)
84  {
85  case 0:
86  case 1: // 1d, 1st order
87  // 1 point rule
88  mWeights.push_back(1);
89  mPoints.push_back(ChastePoint<1>(0.5));
90  break;
91  case 2:
92  case 3: // 1d, 3rd order
93  // 2 point rule
94  mWeights.push_back(0.5);
95  mWeights.push_back(0.5);
96  {
97  double sqrt_one_third = sqrt(1.0/3.0);
98  mPoints.push_back(ChastePoint<1>((-sqrt_one_third+1.0)/2.0));
99  mPoints.push_back(ChastePoint<1>((sqrt_one_third+1.0)/2.0));
100  }
101  break;
102  case 4:
103  case 5: // 1d, 5th order
104  // 3 point rule
105  mWeights.push_back(5.0/18.0);
106  mWeights.push_back(4.0/9.0);
107  mWeights.push_back(5.0/18.0);
108 
109  {
110  double sqrt_three_fifths = sqrt(3.0/5.0);
111  mPoints.push_back(ChastePoint<1>((-sqrt_three_fifths+1.0)/2.0));
112  mPoints.push_back(ChastePoint<1>(0.5));
113  mPoints.push_back(ChastePoint<1>((sqrt_three_fifths+1.0)/2.0));
114  }
115  break;
116  default:
117  EXCEPTION("Gauss quadrature order not supported.");
118  }
119  assert(mPoints.size() == mWeights.size());
120  mNumQuadPoints = mPoints.size();
121 }
122 
129 template<>
131 {
132  double one_third = 1.0/3.0;
133  double one_sixth = 1.0/6.0;
134  double two_thirds = 2.0/3.0;
135  switch (quadratureOrder)
136  {
137  case 0: // 2d, 0th order
138  case 1: // 2d, 1st order
139  // 1 point rule
140  mWeights.push_back(0.5);
141  mPoints.push_back(ChastePoint<2>(one_third, one_third));
142  break;
143 
144  case 2: // 2d, 2nd order
145  // 3 point rule
146  mWeights.push_back(one_sixth);
147  mWeights.push_back(one_sixth);
148  mWeights.push_back(one_sixth);
149 
150  mPoints.push_back(ChastePoint<2>(two_thirds, one_sixth));
151  mPoints.push_back(ChastePoint<2>(one_sixth, one_sixth));
152  mPoints.push_back(ChastePoint<2>(one_sixth, two_thirds));
153  break;
154 
155  case 3: // 2d, 3rd order - derived by hand and using a Macsyma script to solve the cubic
156  // 60*x^3 - 60*x^2 + 15*x - 1;
157  // 6 point rule
158  {
159  double w = 1.0/12.0;
160  mWeights.push_back(w);
161  mWeights.push_back(w);
162  mWeights.push_back(w);
163  mWeights.push_back(w);
164  mWeights.push_back(w);
165  mWeights.push_back(w);
166 
167  double inverse_tan = atan(0.75);
168  double cos_third = cos(inverse_tan/3.0);
169  double sin_third = sin(inverse_tan/3.0);
170  // a = 0.23193336461755
171  double a = sin_third/(2*sqrt(3.0)) - cos_third/6.0 + 1.0/3.0;
172  // b = 0.10903901046622
173  double b = - sin_third/(2*sqrt(3.0)) - cos_third/6.0 + 1.0/3.0;
174  // c = 0.659028
175  double c = cos_third/3.0 + 1.0/3.0;
176 
177  mPoints.push_back(ChastePoint<2>(a, b));
178  mPoints.push_back(ChastePoint<2>(a, c));
179  mPoints.push_back(ChastePoint<2>(b, a));
180  mPoints.push_back(ChastePoint<2>(b, c));
181  mPoints.push_back(ChastePoint<2>(c, a));
182  mPoints.push_back(ChastePoint<2>(c, b));
183  }
184  break;
185  default:
186  EXCEPTION("Gauss quadrature order not supported.");
187  }
188  assert(mPoints.size() == mWeights.size());
189  mNumQuadPoints = mPoints.size();
190 }
191 
197 template<>
199 {
200  switch (quadratureOrder)
201  {
202  case 0: // 3d, 0th order
203  case 1: // 3d, 1st order
204  // 1 point rule
205  mWeights.push_back(1.0/6.0);
206  mPoints.push_back(ChastePoint<3>(0.25, 0.25, 0.25));
207  break;
208 
209  case 2: //2nd order
210  // 4 point rule
211  {
212  double sqrt_fifth = 1.0/sqrt(5.0);
213  double a = (1.0 + 3.0*sqrt_fifth)/4.0; //0.585410196624969;
214  double b = (1.0 - sqrt_fifth)/4.0; //0.138196601125011;
215  double w = 1.0/24.0;
216 
217  mWeights.push_back(w);
218  mWeights.push_back(w);
219  mWeights.push_back(w);
220  mWeights.push_back(w);
221 
222  mPoints.push_back(ChastePoint<3>(a,b,b));
223  mPoints.push_back(ChastePoint<3>(b,a,b));
224  mPoints.push_back(ChastePoint<3>(b,b,a));
225  mPoints.push_back(ChastePoint<3>(b,b,b));
226  }
227  break;
228 
229  case 3: // 3d, 3rd order
230  // 8 point rule
231  /* The main options were
232  * 5-point rule. Commonly published rule has four symmetric points and
233  * a negative weight in the centre. We would like to avoid
234  * negative weight (certainly for interpolation.
235  * 8-point rule. Uses two sets of symmetric points (as 4 point rule with a,b and then with c,d).
236  * This one is hard to derive a closed form solution to.
237  */
238  {
239  double root_seventeen = sqrt(17.0);
240  double root_term = sqrt(1022.0-134.0*root_seventeen);
241  double b = (55.0 - 3.0*root_seventeen + root_term)/196; //b = 0.328055
242  double d = (55.0 - 3.0*root_seventeen - root_term)/196; //d = 0.106952
243 
244  double a = 1.0 - 3.0*b; // a = 0.0158359
245  double c = 1.0 - 3.0*d; // c = 0.679143
246 
247  // w1 = 0.023088 (= 0.138528/6)
248  double w1 = (20.0*d*d - 10.0*d + 1.0)/(240.0*(2.0*d*d - d - 2.0*b*b + b)); // w1 = 0.0362942
249  double w2 = 1.0/24.0 - w1; // w2 = 0.0185787 (=0.111472/6)
250 
251  mWeights.push_back(w1);
252  mWeights.push_back(w1);
253  mWeights.push_back(w1);
254  mWeights.push_back(w1);
255 
256  mWeights.push_back(w2);
257  mWeights.push_back(w2);
258  mWeights.push_back(w2);
259  mWeights.push_back(w2);
260 
261  mPoints.push_back(ChastePoint<3>(a, b, b));
262  mPoints.push_back(ChastePoint<3>(b, a, b));
263  mPoints.push_back(ChastePoint<3>(b, b, a));
264  mPoints.push_back(ChastePoint<3>(b, b, b));
265 
266  mPoints.push_back(ChastePoint<3>(c, d, d));
267  mPoints.push_back(ChastePoint<3>(d, c, d));
268  mPoints.push_back(ChastePoint<3>(d, d, c));
269  mPoints.push_back(ChastePoint<3>(d, d, d));
270 
271  }
272 break;
273 
274  default:
275  EXCEPTION("Gauss quadrature order not supported.");
276  }
277  assert(mPoints.size() == mWeights.size());
278  mNumQuadPoints = mPoints.size();
279 }
280 
281 template<unsigned ELEMENT_DIM>
283 {
284  EXCEPTION("Gauss quadrature rule not available for this dimension.");
285 }
286 
288 // Explicit instantiation
290 
291 template class GaussianQuadratureRule<0>;
292 template class GaussianQuadratureRule<1>;
293 template class GaussianQuadratureRule<2>;
294 template class GaussianQuadratureRule<3>;
295 template class GaussianQuadratureRule<4>;
GaussianQuadratureRule(unsigned quadratureOrder)
std::vector< double > mWeights
#define EXCEPTION(message)
Definition: Exception.hpp:143
unsigned GetNumQuadPoints() const
double GetWeight(unsigned index) const
std::vector< ChastePoint< ELEMENT_DIM > > mPoints
const ChastePoint< ELEMENT_DIM > & rGetQuadPoint(unsigned index) const