PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Quadrature.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2023 PISM Authors
2  *
3  * This file is part of PISM.
4  *
5  * PISM is free software; you can redistribute it and/or modify it under the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13  * details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with PISM; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #include "pism/util/fem/Quadrature.hh"
21 
22 namespace pism {
23 namespace fem {
24 
25 const std::vector<QuadPoint>& Quadrature::points() const {
26  return m_points;
27 }
28 
29 const std::vector<double>& Quadrature::weights() const {
30  return m_weights;
31 }
32 
33 //! Build quadrature points and weights for a tensor product quadrature based on a 1D quadrature
34 //! rule. Uses the same 1D quadrature in both directions.
35 /**
36  @param[in] n 1D quadrature size (the resulting quadrature has size n*n)
37  @param[in] points1 1D quadrature points
38  @param[in] weights1 1D quadrature weights
39  @param[out] points resulting 2D quadrature points
40  @param[out] weights resulting 2D quadrature weights
41  */
42 static void tensor_product_quadrature(unsigned int n,
43  const double *points1,
44  const double *weights1,
45  std::vector<QuadPoint>& points,
46  std::vector<double>& weights) {
47  points.resize(n * n);
48  weights.resize(n * n);
49 
50  unsigned int q = 0;
51  for (unsigned int j = 0; j < n; ++j) {
52  for (unsigned int i = 0; i < n; ++i) {
53  points[q] = {points1[i], points1[j], 0.0};
54  weights[q] = weights1[i] * weights1[j];
55 
56  ++q;
57  }
58  }
59 }
61 
62  // coordinates and weights of the 2-point 1D Gaussian quadrature
63  double A = 1.0 / std::sqrt(3.0);
64 
65  m_points = {{-A, 0.0, 0.0}, {A, 0.0, 0.0}};
66  m_weights = {0.5 * D, 0.5 * D};
67 }
68 
69 //! One-point quadrature on a rectangle.
71  m_points = {{0.0, 0.0, 0.0}};
72  m_weights = {4.0};
73 }
74 
75 //! Two-by-two Gaussian quadrature on a rectangle.
77 
78  // coordinates and weights of the 2-point 1D Gaussian quadrature
79  const double
80  A = 1.0 / sqrt(3.0),
81  points2[2] = {-A, A},
82  weights2[2] = {1.0, 1.0};
83 
84  tensor_product_quadrature(2, points2, weights2, m_points, m_weights);
85 }
86 
88  const double
89  A = 0.0,
90  B = sqrt(0.6),
91  points3[3] = {-B, A, B};
92 
93  const double
94  w1 = 5.0 / 9.0,
95  w2 = 8.0 / 9.0,
96  weights3[3] = {w1, w2, w1};
97 
98  tensor_product_quadrature(3, points3, weights3, m_points, m_weights);
99 }
100 
102  const double
103  A = sqrt(3.0 / 7.0 - (2.0 / 7.0) * sqrt(6.0 / 5.0)), // smaller magnitude
104  B = sqrt(3.0 / 7.0 + (2.0 / 7.0) * sqrt(6.0 / 5.0)), // larger magnitude
105  points4[4] = {-B, -A, A, B};
106 
107  // The weights w_i for Gaussian quadrature on the reference element with these
108  // quadrature points
109  const double
110  w1 = (18.0 + sqrt(30.0)) / 36.0, // larger
111  w2 = (18.0 - sqrt(30.0)) / 36.0, // smaller
112  weights4[4] = {w2, w1, w1, w2};
113 
114  tensor_product_quadrature(4, points4, weights4, m_points, m_weights);
115 }
116 
117 
118 //! @brief N*N-point uniform (*not* Gaussian) quadrature for integrating discontinuous
119 //! functions.
121 
122  std::vector<double> xi(N), w(N);
123  const double dxi = 2.0 / N;
124  for (unsigned int k = 0; k < N; ++k) {
125  xi[k] = -1.0 + dxi*(k + 0.5);
126  w[k] = 2.0 / N;
127  }
128 
129  tensor_product_quadrature(N, xi.data(), w.data(), m_points, m_weights);
130 }
131 
132 /*!
133  * 3-point Gaussian quadrature on the triangle (0,0)-(1,0)-(0,1)
134  */
136 
137  const double
138  one_over_six = 1.0 / 6.0,
139  two_over_three = 2.0 / 3.0;
140 
141  m_points = {{two_over_three, one_over_six, 0.0},
142  {one_over_six, two_over_three, 0.0},
143  {one_over_six, one_over_six, 0.0}};
144 
145  m_weights = {one_over_six, one_over_six, one_over_six};
146 }
147 
149  double xis[8] = {-1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0};
150  double etas[8] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0};
151  double zetas[8] = {-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0};
152 
153  m_points.resize(8);
154  m_weights.resize(8);
155 
156  double C = 1.0 / sqrt(3.0);
157 
158  for (int k = 0; k < 8; ++k) {
159  m_points[k] = {C * xis[k], C * etas[k], C * zetas[k]};
160  m_weights[k] = 1.0;
161  }
162 }
163 
165  m_points = {{0.0, 0.0, 0.0}};
166  m_weights = {8.0};
167 }
168 
170  const double
171  A = sqrt(3.0 / 7.0 - (2.0 / 7.0) * sqrt(6.0 / 5.0)), // smaller magnitude
172  B = sqrt(3.0 / 7.0 + (2.0 / 7.0) * sqrt(6.0 / 5.0)), // larger magnitude
173  pt[4] = {-B, -A, A, B};
174 
175  // The weights w_i for Gaussian quadrature on the reference element with these
176  // quadrature points
177  const double
178  w1 = (18.0 + sqrt(30.0)) / 36.0, // larger
179  w2 = (18.0 - sqrt(30.0)) / 36.0, // smaller
180  w[4] = {w2, w1, w1, w2};
181 
182  m_points.resize(64);
183  m_weights.resize(64);
184  int q = 0;
185  for (int i = 0; i < 4; ++i) {
186  for (int j = 0; j < 4; ++j) {
187  for (int k = 0; k < 4; ++k) {
188  m_points[q] = {pt[i], pt[j], pt[k]};
189  m_weights[q] = w[i] * w[j] * w[k];
190  ++q;
191  }
192  }
193  }
194 }
195 
196 } // end of namespace fem
197 } // end of namespace pism
Q1Quadrature1()
One-point quadrature on a rectangle.
Definition: Quadrature.cc:70
Q1Quadrature4()
Two-by-two Gaussian quadrature on a rectangle.
Definition: Quadrature.cc:76
Q1QuadratureN(unsigned int n)
N*N-point uniform (not Gaussian) quadrature for integrating discontinuous functions.
Definition: Quadrature.cc:120
std::vector< QuadPoint > m_points
Definition: Quadrature.hh:75
const std::vector< double > & weights() const
Definition: Quadrature.cc:29
std::vector< double > m_weights
Definition: Quadrature.hh:76
const std::vector< QuadPoint > & points() const
Definition: Quadrature.cc:25
#define n
Definition: exactTestM.c:37
static void tensor_product_quadrature(unsigned int n, const double *points1, const double *weights1, std::vector< QuadPoint > &points, std::vector< double > &weights)
Definition: Quadrature.cc:42
static const double k
Definition: exactTestP.cc:42
const int C[]
Definition: ssafd_code.cc:44