Nektar++
Loading...
Searching...
No Matches
TestEquiSpaced.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestEquiSpaced.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description:
32//
33///////////////////////////////////////////////////////////////////////////////
34
41#include <boost/test/tools/floating_point_comparison.hpp>
42#include <boost/test/unit_test.hpp>
43
45{
46
47BOOST_AUTO_TEST_CASE(TestQuadExpInterpPhysToEquiSpaced)
48{
49
50 LibUtilities::PointsType PointsTypeDir1 =
53 unsigned int numPoints = 6;
54 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numPoints,
55 PointsTypeDir1);
56 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
57 PointsKeyDir1);
58
61 basisKeyDir1, basisKeyDir1);
62
63 unsigned int numEQ = 4;
64 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
65 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
66 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
68
69 double epsilon = 1.0e-8;
70
71 Exp->GetCoords(c0, c1, c2);
72
73 Exp->PhysInterpToSimplexEquiSpaced(c0, out, numEQ);
74 for (int j = 0, cnt = 0; j < numEQ; ++j)
75 {
76 for (int i = 0; i < numEQ; ++i)
77 {
78 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * i / (numEQ - 1.0),
79 epsilon);
80 }
81 }
82
83 Exp->PhysInterpToSimplexEquiSpaced(c1, out, numEQ);
84 for (int j = 0, cnt = 0; j < numEQ; ++j)
85 {
86 for (int i = 0; i < numEQ; ++i)
87 {
88 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * j / (numEQ - 1.0),
89 epsilon);
90 }
91 }
92}
93
94BOOST_AUTO_TEST_CASE(TestTriExpInterpPhysToEquiSpaced)
95{
96 using namespace LibUtilities;
97
98 PointsType PointsTypeDir1 = eGaussLobattoLegendre;
99 PointsType PointsTypeDir2 = eGaussRadauMAlpha1Beta0;
100 BasisType basisTypeDir1 = eModified_A;
101 BasisType basisTypeDir2 = eModified_B;
102
103 unsigned int numPoints = 6;
104 unsigned int numEQ = 4;
105
106 const PointsKey PointsKeyDir1(numPoints, PointsTypeDir1);
107 const BasisKey basisKeyDir1(basisTypeDir1, numEQ, PointsKeyDir1);
108 const PointsKey PointsKeyDir2(numPoints - 1, PointsTypeDir2);
109 const BasisKey basisKeyDir2(basisTypeDir2, numEQ, PointsKeyDir2);
110
113 basisKeyDir1, basisKeyDir2);
114
115 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
116 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
117 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
119
120 double epsilon = 1.0e-8;
121
122 Exp->GetCoords(c0, c1, c2);
123
124 Exp->PhysInterpToSimplexEquiSpaced(c0, out, numEQ);
125 for (int j = 0, cnt = 0; j < numEQ; ++j)
126 {
127 for (int i = 0; i < numEQ - j; ++i)
128 {
129 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * i / (numEQ - 1.0),
130 epsilon);
131 }
132 }
133
134 Exp->PhysInterpToSimplexEquiSpaced(c1, out, numEQ);
135 for (int j = 0, cnt = 0; j < numEQ; ++j)
136 {
137 for (int i = 0; i < numEQ - j; ++i)
138 {
139 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * j / (numEQ - 1.0),
140 epsilon);
141 }
142 }
143}
144
145BOOST_AUTO_TEST_CASE(TestTetExpInterpPhysToEquiSpaced)
146{
147 using namespace LibUtilities;
148
149 PointsType PointsTypeDir1 = eGaussLobattoLegendre;
150 PointsType PointsTypeDir2 = eGaussRadauMAlpha1Beta0;
151 PointsType PointsTypeDir3 = eGaussRadauMAlpha2Beta0;
152 BasisType basisTypeDir1 = eModified_A;
153 BasisType basisTypeDir2 = eModified_B;
154 BasisType basisTypeDir3 = eModified_C;
155
156 unsigned int numPoints = 10;
157 unsigned int numEQ = 8;
158
159 // Set up standard element.
160 const PointsKey PointsKeyDir1(numPoints, PointsTypeDir1);
161 const PointsKey PointsKeyDir2(numPoints - 1, PointsTypeDir2);
162 const PointsKey PointsKeyDir3(numPoints - 1, PointsTypeDir3);
163 const BasisKey basisKeyDir2(basisTypeDir2, numEQ, PointsKeyDir2);
164 const BasisKey basisKeyDir1(basisTypeDir1, numEQ, PointsKeyDir1);
165 const BasisKey basisKeyDir3(basisTypeDir3, numEQ, PointsKeyDir3);
166
169 basisKeyDir1, basisKeyDir2, basisKeyDir3);
170
171 // define an equispaced Tet points
172 PointsType PointsTypeEq = eNodalTetEvenlySpaced;
173
176 basisKeyDir1, basisKeyDir2, basisKeyDir3, PointsTypeEq);
177
178 // Get coordinates at quadrature points
179 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
180 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
181 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
182 Exp->GetCoords(c0, c1, c2);
183
184 // Get coordinates at equispaced points
185 Array<OneD, NekDouble> ceq0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
186 Array<OneD, NekDouble> ceq1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
187 Array<OneD, NekDouble> ceq2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
188 ExpEq->GetCoords(ceq0, ceq1, ceq2);
189
190 Array<OneD, NekDouble> out = Array<OneD, NekDouble>(Exp->GetNcoeffs());
191 Array<OneD, NekDouble> outeq = Array<OneD, NekDouble>(Exp->GetNcoeffs());
192 double epsilon = 1.0e-8;
193
194 Array<OneD, int> sorted;
196
197 // compare x-coordinates
198 Exp->PhysInterpToSimplexEquiSpaced(c0, out, numEQ);
199 int cnt = 0;
200 for (int k = 0; k < numEQ; ++k)
201 {
202 for (int j = 0; j < numEQ - k; ++j)
203 {
204 for (int i = 0; i < numEQ - k - j; ++i)
205 {
206 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * i / (numEQ - 1.0),
207 epsilon);
208 }
209 }
210 }
211
212 // compare y-coordinates
213 Exp->PhysInterpToSimplexEquiSpaced(c1, out, numEQ);
214 cnt = 0;
215 for (int k = 0; k < numEQ; ++k)
216 {
217 for (int j = 0; j < numEQ - k; ++j)
218 {
219 for (int i = 0; i < numEQ - k - j; ++i)
220 {
221 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * j / (numEQ - 1.0),
222 epsilon);
223 }
224 }
225 }
226
227 // compare z-coordinates
228 Exp->PhysInterpToSimplexEquiSpaced(c2, out, numEQ);
229 cnt = 0;
230 for (int k = 0; k < numEQ; ++k)
231 {
232 for (int j = 0; j < numEQ - k; ++j)
233 {
234 for (int i = 0; i < numEQ - k - j; ++i)
235 {
236 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * k / (numEQ - 1.0),
237 epsilon);
238 }
239 }
240 }
241}
242
243BOOST_AUTO_TEST_CASE(TestPrismExpInterpPhysToEquiSpaced)
244{
245 using namespace LibUtilities;
246
247 PointsType PointsTypeDir1 = eGaussLobattoLegendre;
248 PointsType PointsTypeDir2 = eGaussLobattoLegendre;
249 PointsType PointsTypeDir3 = eGaussRadauMAlpha1Beta0;
250 BasisType basisTypeDir1 = eModified_A;
251 BasisType basisTypeDir2 = eModified_A;
252 BasisType basisTypeDir3 = eModified_B;
253
254 unsigned int numPoints = 10;
255 unsigned int numEQ = 8;
256
257 // Set up standard element.
258 const PointsKey PointsKeyDir1(numPoints, PointsTypeDir1);
259 const PointsKey PointsKeyDir2(numPoints, PointsTypeDir2);
260 const PointsKey PointsKeyDir3(numPoints - 1, PointsTypeDir3);
261 const BasisKey basisKeyDir1(basisTypeDir1, numEQ, PointsKeyDir1);
262 const BasisKey basisKeyDir2(basisTypeDir2, numEQ, PointsKeyDir2);
263 const BasisKey basisKeyDir3(basisTypeDir3, numEQ, PointsKeyDir3);
264
267 basisKeyDir1, basisKeyDir2, basisKeyDir3);
268
269 // define an equispaced Tet points
270 PointsType PointsTypeEq = eNodalPrismEvenlySpaced;
271
274 basisKeyDir1, basisKeyDir2, basisKeyDir3, PointsTypeEq);
275
276 // Get coordinates at quadrature points
277 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
278 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
279 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
280 Exp->GetCoords(c0, c1, c2);
281
282 // Get coordinates at equispaced points
283 Array<OneD, NekDouble> ceq0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
284 Array<OneD, NekDouble> ceq1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
285 Array<OneD, NekDouble> ceq2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
286 ExpEq->GetCoords(ceq0, ceq1, ceq2);
287
288 Array<OneD, NekDouble> out = Array<OneD, NekDouble>(Exp->GetNcoeffs());
289 Array<OneD, NekDouble> outeq = Array<OneD, NekDouble>(Exp->GetNcoeffs());
290 double epsilon = 1.0e-8;
291
292 Array<OneD, int> sorted;
294
295 // compare x-coordinates
296 Exp->PhysInterpToSimplexEquiSpaced(c0, out, numEQ);
297 int cnt = 0;
298 for (int k = 0; k < numEQ; ++k)
299 {
300 for (int j = 0; j < numEQ; ++j)
301 {
302 for (int i = 0; i < numEQ - k; ++i)
303 {
304 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * i / (numEQ - 1.0),
305 epsilon);
306 }
307 }
308 }
309
310 // compare y-coordinates
311 Exp->PhysInterpToSimplexEquiSpaced(c1, out, numEQ);
312 cnt = 0;
313 for (int k = 0; k < numEQ; ++k)
314 {
315 for (int j = 0; j < numEQ; ++j)
316 {
317 for (int i = 0; i < numEQ - k; ++i)
318 {
319 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * j / (numEQ - 1.0),
320 epsilon);
321 }
322 }
323 }
324
325 // compare z-coordinates
326 Exp->PhysInterpToSimplexEquiSpaced(c2, out, numEQ);
327 cnt = 0;
328 for (int k = 0; k < numEQ; ++k)
329 {
330 for (int j = 0; j < numEQ; ++j)
331 {
332 for (int i = 0; i < numEQ - k; ++i)
333 {
334 BOOST_CHECK_CLOSE(out[cnt++], -1.0 + 2.0 * k / (numEQ - 1.0),
335 epsilon);
336 }
337 }
338 }
339}
340
341BOOST_AUTO_TEST_CASE(TestQuadExpEquiSpacedToPhys)
342{
343 LibUtilities::PointsType PointsTypeDir1 =
346 unsigned int numPoints = 10;
347 unsigned int numModes = 8;
348 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numPoints,
349 PointsTypeDir1);
350 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
351 PointsKeyDir1);
352
355 basisKeyDir1, basisKeyDir1);
356
357 unsigned int numEQ = 4;
358 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
359 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
360 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
361 Array<OneD, NekDouble> phys = Array<OneD, NekDouble>(Exp->GetTotPoints());
363
364 double epsilon = 1.0e-8;
365
366 Exp->GetCoords(c0, c1, c2);
367
368 for (int j = 0, cnt = 0; j < numEQ; ++j)
369 {
370 for (int i = 0; i < numEQ; ++i)
371 {
372 in[cnt++] = -1.0 + 2.0 * i / (numEQ - 1.0);
373 }
374 }
375
376 Exp->EquiSpacedToPhys(numEQ, in, phys);
377 for (int j = 0, cnt = 0; j < numPoints; ++j)
378 {
379 for (int i = 0; i < numPoints; ++i)
380 {
381 BOOST_CHECK_CLOSE(phys[cnt], c0[cnt], epsilon);
382 cnt++;
383 }
384 }
385
386 for (int j = 0, cnt = 0; j < numEQ; ++j)
387 {
388 for (int i = 0; i < numEQ; ++i)
389 {
390 in[cnt++] = -1.0 + 2.0 * j / (numEQ - 1.0);
391 }
392 }
393 Exp->EquiSpacedToPhys(numEQ, in, phys);
394 for (int j = 0, cnt = 0; j < numPoints; ++j)
395 {
396 for (int i = 0; i < numPoints; ++i)
397 {
398 BOOST_CHECK_CLOSE(phys[cnt], c1[cnt], epsilon);
399 cnt++;
400 }
401 }
402}
403
404BOOST_AUTO_TEST_CASE(TestTriExpEquiSpacedToPhys)
405{
406
407 LibUtilities::PointsType PointsTypeDir1 =
409 LibUtilities::PointsType PointsTypeDir2 =
410 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
411 Nektar::LibUtilities::BasisType basisTypeDir1 =
413 Nektar::LibUtilities::BasisType basisTypeDir2 =
415
416 unsigned int numPoints = 6;
417
418 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numPoints,
419 PointsTypeDir1);
420 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
421 PointsKeyDir1);
422 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numPoints - 1,
423 PointsTypeDir2);
424 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
425 PointsKeyDir2);
426
429 basisKeyDir1, basisKeyDir2);
430
431 unsigned int numEQ = 4;
432 Array<OneD, NekDouble> c0 = Array<OneD, NekDouble>(Exp->GetTotPoints());
433 Array<OneD, NekDouble> c1 = Array<OneD, NekDouble>(Exp->GetTotPoints());
434 Array<OneD, NekDouble> c2 = Array<OneD, NekDouble>(Exp->GetTotPoints());
435 Array<OneD, NekDouble> phys = Array<OneD, NekDouble>(Exp->GetTotPoints());
437
438 double epsilon = 1.0e-8;
439
440 Exp->GetCoords(c0, c1, c2);
441
442 for (int j = 0, cnt = 0; j < numEQ; ++j)
443 {
444 for (int i = 0; i < numEQ - j; ++i)
445 {
446 in[cnt++] = -1.0 + 2.0 * i / (numEQ - 1.0);
447 }
448 }
449 Exp->EquiSpacedToPhys(numEQ, in, phys);
450 for (int j = 0, cnt = 0; j < numPoints - 1; ++j)
451 {
452 for (int i = 0; i < numPoints; ++i)
453 {
454 BOOST_CHECK_CLOSE(phys[cnt], c0[cnt], epsilon);
455 cnt++;
456 }
457 }
458
459 for (int j = 0, cnt = 0; j < numEQ; ++j)
460 {
461 for (int i = 0; i < numEQ - j; ++i)
462 {
463 in[cnt++] = -1.0 + 2.0 * j / (numEQ - 1.0);
464 }
465 }
466 Exp->EquiSpacedToPhys(numEQ, in, phys);
467 for (int j = 0, cnt = 0; j < numPoints - 1; ++j)
468 {
469 for (int i = 0; i < numPoints; ++i)
470 {
471 BOOST_CHECK_CLOSE(phys[cnt], c1[cnt], epsilon);
472 cnt++;
473 }
474 }
475}
476
477} // namespace Nektar::QuadEquiSpaced
Describes the specification for a Basis.
Definition Basis.h:45
static void CartesianOrdering(const int nq, Array< OneD, int > &sorted)
Definition NodalUtil.h:363
static void CartesianOrdering(const int nq, Array< OneD, int > &sorted)
Definition NodalUtil.h:256
Defines a specification for a set of points.
Definition Points.h:50
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
@ eNodalPrismEvenlySpaced
3D Evenly-spaced points on a Prism
Definition PointsType.h:86
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition PointsType.h:51
@ eNodalTetEvenlySpaced
3D Evenly-spaced points on a Tetrahedron
Definition PointsType.h:84
@ eModified_B
Principle Modified Functions .
Definition BasisType.h:49
@ eModified_C
Principle Modified Functions .
Definition BasisType.h:50
@ eModified_A
Principle Modified Functions .
Definition BasisType.h:48
BOOST_AUTO_TEST_CASE(TestQuadExpInterpPhysToEquiSpaced)
std::shared_ptr< StdPrismExp > StdPrismExpSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition StdTetExp.h:187
std::shared_ptr< StdQuadExp > StdQuadExpSharedPtr
Definition StdQuadExp.h:181
std::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition StdTriExp.h:177