29 #ifndef CORE_BASIS_TETSAMPLINGSCHEMES_H
30 #define CORE_BASIS_TETSAMPLINGSCHEMES_H 1
46 template <
class ARRAY1,
class ARRAY2>
49 typedef typename ARRAY1::value_type coords_type;
50 if (order == 1 || order == 0)
52 const double gaussian_weights[1] = { 1.0};
53 const double gaussian_coords[1][3] = { {1./4., 1./4., 1./4.}};
54 const unsigned int num_coords = 3;
55 const unsigned int num_points = 1;
57 coords.resize(num_points);
58 weights.resize(num_points);
59 for (
unsigned int i=0; i<num_points; i++)
61 coords[i].resize(num_coords);
62 for (
unsigned int j=0; j<num_coords; j++)
63 coords[i][j] = static_cast<typename coords_type::value_type>(gaussian_coords[i][j]);
64 weights[i] =
static_cast<typename ARRAY2::value_type
>(gaussian_weights[i]);
69 const double gaussian_weights[4] = {.25, .25, .25, .25};
70 const double gaussian_coords[4][3] = {
71 {0.138196601125011, 0.138196601125011, 0.138196601125011},
72 {0.585410196624969, 0.138196601125011, 0.138196601125011},
73 {0.138196601125011, 0.585410196624969, 0.138196601125011},
74 {0.138196601125011, 0.138196601125011, 0.585410196624969}};
75 const unsigned int num_coords = 3;
79 const unsigned int num_points = 4;
81 coords.resize(num_points);
82 weights.resize(num_points);
83 for (
unsigned int i=0; i<num_points; i++)
85 coords[i].resize(num_coords);
86 for (
unsigned int j=0; j<num_coords; j++)
87 coords[i][j] = static_cast<typename coords_type::value_type>(gaussian_coords[i][j]);
88 weights[i] =
static_cast<typename ARRAY2::value_type
>(gaussian_weights[i]);
93 const double gaussian_weights[11] = {
108 const double gaussian_coords[11][3] = {
109 {0.2500000, 0.2500000, 0.2500000},
110 {0.7857143, 0.07142857, 0.07142857},
111 {0.07142857, 0.7857143, 0.07142857},
112 {0.07142857, 0.07142857, 0.7857143},
113 {0.07142857, 0.07142857, 0.07142857},
114 {0.1005964, 0.1005964, 0.3994034},
115 {0.1005964, 0.3994034, 0.1005964},
116 {0.1005964, 0.3994034, 0.3994034},
117 {0.3994034, 0.1005964, 0.1005964},
118 {0.3994034, 0.1005964, 0.3994034},
119 {0.3994034, 0.3994034, 0.1005964}};
120 const unsigned int num_coords = 3;
126 const unsigned int num_points = 11;
128 coords.resize(num_points);
129 weights.resize(num_points);
130 for (
unsigned int i=0; i<num_points; i++)
132 coords[i].resize(num_coords);
133 for (
unsigned int j=0; j<num_coords; j++)
134 coords[i][j] = static_cast<typename coords_type::value_type>(gaussian_coords[i][j]);
135 weights[i] =
static_cast<typename ARRAY2::value_type
>(gaussian_weights[i]);
144 template <
class ARRAY1,
class ARRAY2>
147 typedef typename ARRAY1::value_type coords_type;
152 for (
int p=0; p<order;p++)
153 for (
int q=0; q<order;q++)
154 for (
int r=0; r<order;r++)
156 if (p+q+r == (order-1))
159 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/4.0)/
static_cast<double>(order));
160 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/4.0)/
static_cast<double>(order));
161 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/4.0)/
static_cast<double>(order));
162 coords.push_back(c); m++;
164 if (p+q+r == (order-2))
167 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+2.0/4.0)/
static_cast<double>(order));
168 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/4.0)/
static_cast<double>(order));
169 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/4.0)/
static_cast<double>(order));
170 coords.push_back(c); m++;
171 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/4.0)/
static_cast<double>(order));
172 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+2.0/4.0)/
static_cast<double>(order));
173 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/4.0)/
static_cast<double>(order));
174 coords.push_back(c); m++;
175 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/4.0)/
static_cast<double>(order));
176 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/4.0)/
static_cast<double>(order));
177 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+2.0/4.0)/
static_cast<double>(order));
178 coords.push_back(c); m++;
179 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/3.0)/
static_cast<double>(order));
180 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/3.0)/
static_cast<double>(order));
181 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/3.0)/
static_cast<double>(order));
182 coords.push_back(c); m++;
184 if (p+q+r < (order-2))
187 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+2.0/4.0)/
static_cast<double>(order));
188 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/4.0)/
static_cast<double>(order));
189 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/4.0)/
static_cast<double>(order));
190 coords.push_back(c); m++;
191 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/4.0)/
static_cast<double>(order));
192 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+2.0/4.0)/
static_cast<double>(order));
193 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/4.0)/
static_cast<double>(order));
194 coords.push_back(c); m++;
195 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/4.0)/
static_cast<double>(order));
196 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/4.0)/
static_cast<double>(order));
197 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+2.0/4.0)/
static_cast<double>(order));
198 coords.push_back(c); m++;
199 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+1.0/3.0)/
static_cast<double>(order));
200 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+1.0/3.0)/
static_cast<double>(order));
201 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+1.0/3.0)/
static_cast<double>(order));
202 coords.push_back(c); m++;
203 c[0] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(p)+2.0/4.0)/
static_cast<double>(order));
204 c[1] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(q)+2.0/4.0)/
static_cast<double>(order));
205 c[2] =
static_cast<typename coords_type::value_type
>((
static_cast<double>(r)+2.0/4.0)/
static_cast<double>(order));
206 coords.push_back(c); m++;
211 for (
int p=0; p<m;p++) weights[p] = static_cast<typename ARRAY2::value_type>(1.0/
static_cast<double>(m));
void get_regular_scheme(ARRAY1 &coords, ARRAY2 &weights, int order)
Definition: TetSamplingSchemes.h:145
#define SCISHARE
Definition: share.h:39
Utility for specifying data invariants (Assertions)
void get_gaussian_scheme(ARRAY1 &coords, ARRAY2 &weights, int order)
Definition: TetSamplingSchemes.h:47
#define REPORT_NOT_IMPLEMENTED(message)
Definition: Exception.h:106
Definition: TetSamplingSchemes.h:42