Nektar++
Loading...
Searching...
No Matches
TestHexCollection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestHexCollection.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
37#include <LocalRegions/HexExp.h>
39#include <boost/test/tools/floating_point_comparison.hpp>
40#include <boost/test/unit_test.hpp>
42{
43
47{
48 std::array<SpatialDomains::PointGeom *, 2> vertices = {v0, v1};
50 new SpatialDomains::SegGeom(id, 3, vertices));
51 return result;
52}
53
55 std::array<SpatialDomains::PointGeom *, 8> v,
56 std::array<SpatialDomains::SegGeomUniquePtr, 12> &segVec,
57 std::array<SpatialDomains::QuadGeomUniquePtr, 6> &faceVec)
58{
59 std::array<std::array<int, 2>, 12> edgeVerts = {{{{0, 1}},
60 {{1, 2}},
61 {{2, 3}},
62 {{3, 0}},
63 {{0, 4}},
64 {{1, 5}},
65 {{2, 6}},
66 {{3, 7}},
67 {{4, 5}},
68 {{5, 6}},
69 {{6, 7}},
70 {{7, 4}}}};
71 std::array<std::array<int, 4>, 6> faceEdges = {{{{0, 1, 2, 3}},
72 {{0, 5, 8, 4}},
73 {{1, 6, 9, 5}},
74 {{2, 6, 10, 7}},
75 {{3, 7, 11, 4}},
76 {{8, 9, 10, 11}}}};
77
78 // Create segments from vertices
79 for (int i = 0; i < 12; ++i)
80 {
81 segVec[i] = CreateSegGeom(i, v[edgeVerts[i][0]], v[edgeVerts[i][1]]);
82 }
83
84 // Create faces from edges
85 std::array<SpatialDomains::QuadGeom *, 6> faces;
86 for (int i = 0; i < 6; ++i)
87 {
88 std::array<SpatialDomains::SegGeom *, 4> face;
89 for (int j = 0; j < 4; ++j)
90 {
91 face[j] = segVec[faceEdges[i][j]].get();
92 }
94 new SpatialDomains::QuadGeom(i, face));
95 faces[i] = faceVec[i].get();
96 }
97
99 new SpatialDomains::HexGeom(0, faces));
100 return hexGeom;
101}
102
103BOOST_AUTO_TEST_CASE(TestHexBwdTrans_StdMat_UniformP)
104{
106 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
108 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
110 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
112 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
114 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
116 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
118 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
120 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
121
122 std::array<SpatialDomains::PointGeom *, 8> v = {
123 v0.get(), v1.get(), v2.get(), v3.get(),
124 v4.get(), v5.get(), v6.get(), v7.get()};
125 std::array<SpatialDomains::SegGeomUniquePtr, 12> segVec;
126 std::array<SpatialDomains::QuadGeomUniquePtr, 6> faceVec;
127 SpatialDomains::HexGeomUniquePtr hexGeom = CreateHex(v, segVec, faceVec);
128
129 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
131 Nektar::LibUtilities::BasisType basisTypeDir1 =
133 unsigned int numQuadPoints = 6;
134 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
135 quadPointsTypeDir1);
136 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
137 quadPointsKeyDir1);
138
141 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
142
143 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
144 CollExp.push_back(Exp);
145
147 Collections::CollectionOptimisation colOpt(dummySession, 3,
149 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
150 Collections::Collection c(CollExp, impTypes);
152
153 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
154 for (int i = 0; i < coeffs.size(); ++i)
155 {
156 coeffs[i] = i + 1;
157 }
158 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
159 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
160
161 Exp->BwdTrans(coeffs, phys1);
162 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
163
164 double epsilon = 1.0e-8;
165 for (int i = 0; i < phys1.size(); ++i)
166 {
167 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
168 }
169}
170#if 0
171BOOST_AUTO_TEST_CASE(TestHexBwdTrans_StdMat_VariableP)
172{
174 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
176 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
178 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
180 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
182 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
184 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
186 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
188 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
189
190 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
191 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
193 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
194 v6.get(), v7.get(), segVec, faceVec);
195
196 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
198 Nektar::LibUtilities::BasisType basisTypeDir1 =
200 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
201 quadPointsTypeDir1);
202 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
203 quadPointsTypeDir1);
204 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
205 quadPointsTypeDir1);
206 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
207 quadPointsKeyDir1);
208 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
209 quadPointsKeyDir2);
210 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
211 quadPointsKeyDir3);
212
215 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
216
217 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
218 CollExp.push_back(Exp);
219
221 Collections::CollectionOptimisation colOpt(dummySession, 3,
223 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
224 Collections::Collection c(CollExp, impTypes);
225 c.Initialise(Collections::eBwdTrans);
226
227 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
228 for (int i = 0; i < coeffs.size(); ++i)
229 {
230 coeffs[i] = i + 1;
231 }
232 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
233 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
234
235 Exp->BwdTrans(coeffs, phys1);
236 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
237
238 double epsilon = 1.0e-8;
239 for (int i = 0; i < phys1.size(); ++i)
240 {
241 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
242 }
243}
244
245BOOST_AUTO_TEST_CASE(TestHexBwdTrans_IterPerExp_UniformP)
246{
248 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
250 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
252 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
254 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
256 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
258 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
260 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
262 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
263
264 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
265 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
267 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
268 v6.get(), v7.get(), segVec, faceVec);
269
270 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
272 Nektar::LibUtilities::BasisType basisTypeDir1 =
274 unsigned int numQuadPoints = 6;
275 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
276 quadPointsTypeDir1);
277 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
278 quadPointsKeyDir1);
279
282 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
283
284 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
285 CollExp.push_back(Exp);
286
288 Collections::CollectionOptimisation colOpt(dummySession, 3,
290 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
291 Collections::Collection c(CollExp, impTypes);
292 c.Initialise(Collections::eBwdTrans);
293
294 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
295 for (int i = 0; i < coeffs.size(); ++i)
296 {
297 coeffs[i] = i + 1;
298 }
299 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
300 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
301
302 Exp->BwdTrans(coeffs, phys1);
303 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
304
305 double epsilon = 1.0e-8;
306 for (int i = 0; i < phys1.size(); ++i)
307 {
308 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
309 }
310}
311
312BOOST_AUTO_TEST_CASE(TestHexBwdTrans_IterPerExp_VariableP)
313{
315 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
317 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
319 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
321 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
323 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
325 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
327 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
329 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
330
331 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
332 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
334 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
335 v6.get(), v7.get(), segVec, faceVec);
336
337 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
339 Nektar::LibUtilities::BasisType basisTypeDir1 =
341 unsigned int numQuadPoints = 6;
342 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
343 quadPointsTypeDir1);
344 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
345 quadPointsKeyDir1);
346 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
347 quadPointsKeyDir1);
348 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
349 quadPointsKeyDir1);
350
353 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
354
355 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
356 CollExp.push_back(Exp);
357
359 Collections::CollectionOptimisation colOpt(dummySession, 3,
361 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
362 Collections::Collection c(CollExp, impTypes);
363 c.Initialise(Collections::eBwdTrans);
364
365 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
366 for (int i = 0; i < coeffs.size(); ++i)
367 {
368 coeffs[i] = i + 1;
369 }
370 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
371 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
372
373 Exp->BwdTrans(coeffs, phys1);
374 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
375 c.Initialise(Collections::eBwdTrans);
376
377 double epsilon = 1.0e-8;
378 for (int i = 0; i < phys1.size(); ++i)
379 {
380 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
381 }
382}
383
384BOOST_AUTO_TEST_CASE(TestHexBwdTrans_IterPerExp_VariableP_MultiElmt)
385{
387 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
389 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
391 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
393 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
395 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
397 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
399 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
401 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
402
403 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
404 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
406 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
407 v6.get(), v7.get(), segVec, faceVec);
408
409 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
411 Nektar::LibUtilities::BasisType basisTypeDir1 =
413 unsigned int numQuadPoints = 6;
414 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
415 quadPointsTypeDir1);
416 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
417 quadPointsKeyDir1);
418 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
419 quadPointsKeyDir1);
420 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
421 quadPointsKeyDir1);
422
425 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
426
427 int nelmts = 10;
428
429 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
430 for (int i = 0; i < nelmts; ++i)
431 {
432 CollExp.push_back(Exp);
433 }
434
436 Collections::CollectionOptimisation colOpt(dummySession, 3,
438 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
439 Collections::Collection c(CollExp, impTypes);
440 c.Initialise(Collections::eBwdTrans);
441
442 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
443 for (int i = 0; i < coeffs.size(); ++i)
444 {
445 coeffs[i] = i + 1;
446 }
447 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
448 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
449
450 for (int i = 0; i < nelmts; ++i)
451 {
452 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
453 tmp = phys1 + i * Exp->GetTotPoints());
454 }
455
456 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
457
458 double epsilon = 1.0e-8;
459 for (int i = 0; i < phys1.size(); ++i)
460 {
461 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
462 }
463}
464
465BOOST_AUTO_TEST_CASE(TestHexBwdTrans_NoCollection_VariableP)
466{
468 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
470 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
472 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
474 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
476 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
478 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
480 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
482 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
483
484 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
485 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
487 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
488 v6.get(), v7.get(), segVec, faceVec);
489
490 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
492 Nektar::LibUtilities::BasisType basisTypeDir1 =
494 unsigned int numQuadPoints = 6;
495 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
496 quadPointsTypeDir1);
497 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
498 quadPointsKeyDir1);
499 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
500 quadPointsKeyDir1);
501 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
502 quadPointsKeyDir1);
503
506 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
507
508 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
509 CollExp.push_back(Exp);
510
512 Collections::CollectionOptimisation colOpt(dummySession, 3,
514 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
515 Collections::Collection c(CollExp, impTypes);
516 c.Initialise(Collections::eBwdTrans);
517
518 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
519 for (int i = 0; i < coeffs.size(); ++i)
520 {
521 coeffs[i] = i + 1;
522 }
523 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
524 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
525
526 Exp->BwdTrans(coeffs, phys1);
527
528 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
529
530 double epsilon = 1.0e-8;
531 for (int i = 0; i < phys1.size(); ++i)
532 {
533 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
534 }
535}
536
537BOOST_AUTO_TEST_CASE(TestHexBwdTrans_SumFac_UniformP)
538{
540 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
542 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
544 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
546 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
548 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
550 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
552 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
554 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
555
556 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
557 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
559 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
560 v6.get(), v7.get(), segVec, faceVec);
561
562 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
564 Nektar::LibUtilities::BasisType basisTypeDir1 =
566 unsigned int numQuadPoints = 6;
567 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
568 quadPointsTypeDir1);
569 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
570 quadPointsKeyDir1);
571
574 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
575
576 int nelmts = 1;
577
578 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
579 for (int i = 0; i < nelmts; ++i)
580 {
581 CollExp.push_back(Exp);
582 }
583
585 Collections::CollectionOptimisation colOpt(dummySession, 3,
587 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
588 Collections::Collection c(CollExp, impTypes);
589 c.Initialise(Collections::eBwdTrans);
590
591 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
592 for (int i = 0; i < coeffs.size(); ++i)
593 {
594 coeffs[i] = i + 1;
595 }
596 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
597 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
598
599 for (int i = 0; i < nelmts; ++i)
600 {
601 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
602 tmp = phys1 + i * Exp->GetTotPoints());
603 }
604 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
605
606 double epsilon = 1.0e-8;
607 for (int i = 0; i < phys1.size(); ++i)
608 {
609 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
610 }
611}
612
613BOOST_AUTO_TEST_CASE(TestHexBwdTrans_SumFac_UniformP_MultiElmt)
614{
616 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
618 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
620 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
622 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
624 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
626 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
628 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
630 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
631
632 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
633 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
635 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
636 v6.get(), v7.get(), segVec, faceVec);
637
638 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
640 Nektar::LibUtilities::BasisType basisTypeDir1 =
642 unsigned int numQuadPoints = 6;
643 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
644 quadPointsTypeDir1);
645 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
646 quadPointsKeyDir1);
647
650 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
651
652 int nelmts = 10;
653
654 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
655 for (int i = 0; i < nelmts; ++i)
656 {
657 CollExp.push_back(Exp);
658 }
659
661 Collections::CollectionOptimisation colOpt(dummySession, 3,
663 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
664 Collections::Collection c(CollExp, impTypes);
665 c.Initialise(Collections::eBwdTrans);
666
667 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
668 for (int i = 0; i < coeffs.size(); ++i)
669 {
670 coeffs[i] = i + 1;
671 }
672 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
673 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
674
675 for (int i = 0; i < nelmts; ++i)
676 {
677 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
678 tmp = phys1 + i * Exp->GetTotPoints());
679 }
680 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
681
682 double epsilon = 1.0e-8;
683 for (int i = 0; i < phys1.size(); ++i)
684 {
685 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
686 }
687}
688
689BOOST_AUTO_TEST_CASE(TestHexBwdTrans_SumFac_VariableP)
690{
692 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
694 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
696 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
698 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
700 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
702 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
704 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
706 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
707
708 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
709 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
711 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
712 v6.get(), v7.get(), segVec, faceVec);
713
714 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
716 Nektar::LibUtilities::BasisType basisTypeDir1 =
718 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
719 quadPointsTypeDir1);
720 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
721 quadPointsTypeDir1);
722 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
723 quadPointsTypeDir1);
724 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
725 quadPointsKeyDir1);
726 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
727 quadPointsKeyDir2);
728 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
729 quadPointsKeyDir3);
730
733 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
734
735 int nelmts = 1;
736
737 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
738 for (int i = 0; i < nelmts; ++i)
739 {
740 CollExp.push_back(Exp);
741 }
742
744 Collections::CollectionOptimisation colOpt(dummySession, 3,
746 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
747 Collections::Collection c(CollExp, impTypes);
748 c.Initialise(Collections::eBwdTrans);
749
750 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
751 for (int i = 0; i < coeffs.size(); ++i)
752 {
753 coeffs[i] = i + 1;
754 }
755 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
756 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
757
758 for (int i = 0; i < nelmts; ++i)
759 {
760 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
761 tmp = phys1 + i * Exp->GetTotPoints());
762 }
763 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
764
765 double epsilon = 1.0e-8;
766 for (int i = 0; i < phys1.size(); ++i)
767 {
768 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
769 }
770}
771
772BOOST_AUTO_TEST_CASE(TestHexBwdTrans_SumFac_VariableP_MultiElmt)
773{
775 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
777 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
779 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
781 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
783 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
785 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
787 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
789 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
790
791 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
792 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
794 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
795 v6.get(), v7.get(), segVec, faceVec);
796
797 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
799 Nektar::LibUtilities::BasisType basisTypeDir1 =
801 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
802 quadPointsTypeDir1);
803 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
804 quadPointsTypeDir1);
805 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
806 quadPointsTypeDir1);
807 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
808 quadPointsKeyDir1);
809 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
810 quadPointsKeyDir2);
811 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
812 quadPointsKeyDir3);
813
816 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
817
818 int nelmts = 10;
819
820 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
821 for (int i = 0; i < nelmts; ++i)
822 {
823 CollExp.push_back(Exp);
824 }
825
827 Collections::CollectionOptimisation colOpt(dummySession, 3,
829 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
830 Collections::Collection c(CollExp, impTypes);
831 c.Initialise(Collections::eBwdTrans);
832
833 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
834 for (int i = 0; i < coeffs.size(); ++i)
835 {
836 coeffs[i] = i + 1;
837 }
838 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
839 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
840
841 for (int i = 0; i < nelmts; ++i)
842 {
843 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
844 tmp = phys1 + i * Exp->GetTotPoints());
845 }
846 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
847
848 double epsilon = 1.0e-8;
849 for (int i = 0; i < phys1.size(); ++i)
850 {
851 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
852 }
853}
854
855BOOST_AUTO_TEST_CASE(TestHexBwdTrans_MatrixFree_UniformP)
856{
858 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
860 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
862 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
864 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
866 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
868 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
870 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
872 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
873
874 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
875 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
877 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
878 v6.get(), v7.get(), segVec, faceVec);
879
880 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
882 Nektar::LibUtilities::BasisType basisTypeDir1 =
884 unsigned int numQuadPoints = 6;
885 unsigned int numModes = 4;
886 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
887 quadPointsTypeDir1);
888 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
889 quadPointsKeyDir1);
890
893 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
894
895 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
896 CollExp.push_back(Exp);
897
899 Collections::CollectionOptimisation colOpt(dummySession, 3,
901 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
902 Collections::Collection c(CollExp, impTypes);
903 c.Initialise(Collections::eBwdTrans);
904
905 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
906 for (int i = 0; i < coeffs.size(); ++i)
907 {
908 coeffs[i] = i + 1;
909 }
910 Array<OneD, NekDouble> physRef(Exp->GetTotPoints());
911 Array<OneD, NekDouble> phys(Exp->GetTotPoints());
912
913 Exp->BwdTrans(coeffs, physRef);
914 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys);
915
916 double epsilon = 1.0e-8;
917 for (int i = 0; i < physRef.size(); ++i)
918 {
919 BOOST_CHECK_CLOSE(physRef[i], phys[i], epsilon);
920 }
921}
922
923BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_StdMat_UniformP)
924{
926 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
928 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
930 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
932 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
934 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
936 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
938 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
940 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
941
942 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
943 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
945 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
946 v6.get(), v7.get(), segVec, faceVec);
947
948 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
950 Nektar::LibUtilities::BasisType basisTypeDir1 =
952 unsigned int numQuadPoints = 6;
953 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
954 quadPointsTypeDir1);
955 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
956 quadPointsKeyDir1);
957
960 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
961
962 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
963 CollExp.push_back(Exp);
964
966 Collections::CollectionOptimisation colOpt(dummySession, 3,
968 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
969 Collections::Collection c(CollExp, impTypes);
970 c.Initialise(Collections::eIProductWRTBase);
971
972 const int nq = Exp->GetTotPoints();
973 Array<OneD, NekDouble> phys(nq);
974 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
975 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
976
977 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
978
979 Exp->GetCoords(xc, yc, zc);
980
981 for (int i = 0; i < nq; ++i)
982 {
983 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
984 }
985
986 Exp->IProductWRTBase(phys, coeffs1);
987 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
988
989 double epsilon = 1.0e-8;
990 for (int i = 0; i < coeffs1.size(); ++i)
991 {
992 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
993 }
994}
995
996BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_MatrixFree_UniformP_Undeformed)
997{
999 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1001 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1003 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1005 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1007 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1009 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1011 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1013 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1014
1015 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1016 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1018 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1019 v6.get(), v7.get(), segVec, faceVec);
1020
1021 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1023 Nektar::LibUtilities::BasisType basisTypeDir1 =
1025 unsigned int numQuadPoints = 5;
1026 unsigned int numModes = 4;
1027 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
1028 quadPointsTypeDir1);
1029 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1030 quadPointsKeyDir1);
1031
1034 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
1035
1036 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1037 CollExp.push_back(Exp);
1038
1040 Collections::CollectionOptimisation colOpt(dummySession, 3,
1042 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1043 Collections::Collection c(CollExp, impTypes);
1044 c.Initialise(Collections::eIProductWRTBase);
1045
1046 const int nq = Exp->GetTotPoints();
1047 Array<OneD, NekDouble> phys(nq);
1048 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1049 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1050
1051 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1052
1053 Exp->GetCoords(xc, yc, zc);
1054
1055 for (int i = 0; i < nq; ++i)
1056 {
1057 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1058 }
1059
1060 Exp->IProductWRTBase(phys, coeffsRef);
1061 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1062
1063 double epsilon = 1.0e-8;
1064 for (int i = 0; i < coeffsRef.size(); ++i)
1065 {
1066 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
1067 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
1068 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1069 }
1070}
1071
1072BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_MatrixFree_UniformP_Deformed)
1073{
1075 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1077 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1079 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1081 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1083 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1085 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1087 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
1089 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1090
1091 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1092 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1094 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1095 v6.get(), v7.get(), segVec, faceVec);
1096
1097 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1099 Nektar::LibUtilities::BasisType basisTypeDir1 =
1101 unsigned int numQuadPoints = 5;
1102 unsigned int numModes = 4;
1103 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
1104 quadPointsTypeDir1);
1105 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1106 quadPointsKeyDir1);
1107
1110 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
1111
1112 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1113 CollExp.push_back(Exp);
1114
1116 Collections::CollectionOptimisation colOpt(dummySession, 3,
1118 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1119 Collections::Collection c(CollExp, impTypes);
1120 c.Initialise(Collections::eIProductWRTBase);
1121
1122 const int nq = Exp->GetTotPoints();
1123 Array<OneD, NekDouble> phys(nq);
1124 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1125 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1126
1127 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1128
1129 Exp->GetCoords(xc, yc, zc);
1130
1131 for (int i = 0; i < nq; ++i)
1132 {
1133 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1134 }
1135
1136 Exp->IProductWRTBase(phys, coeffsRef);
1137 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1138
1139 double epsilon = 1.0e-8;
1140 for (int i = 0; i < coeffsRef.size(); ++i)
1141 {
1142 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
1143 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
1144 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1145 }
1146}
1147
1149 TestHexIProductWRTBase_MatrixFree_UniformP_Deformed_OverInt)
1150{
1152 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1154 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1156 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1158 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1160 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1162 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1164 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
1166 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1167
1168 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1169 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1171 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1172 v6.get(), v7.get(), segVec, faceVec);
1173
1174 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1176 Nektar::LibUtilities::BasisType basisTypeDir1 =
1178 unsigned int numQuadPoints = 8;
1179 unsigned int numModes = 4;
1180 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
1181 quadPointsTypeDir1);
1182 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1183 quadPointsKeyDir1);
1184
1187 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
1188
1189 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1190 CollExp.push_back(Exp);
1191
1193 Collections::CollectionOptimisation colOpt(dummySession, 3,
1195 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1196 Collections::Collection c(CollExp, impTypes);
1197 c.Initialise(Collections::eIProductWRTBase);
1198
1199 const int nq = Exp->GetTotPoints();
1200 Array<OneD, NekDouble> phys(nq);
1201 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1202 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1203
1204 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1205
1206 Exp->GetCoords(xc, yc, zc);
1207
1208 for (int i = 0; i < nq; ++i)
1209 {
1210 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1211 }
1212
1213 Exp->IProductWRTBase(phys, coeffsRef);
1214 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1215
1216 double epsilon = 1.0e-8;
1217 for (int i = 0; i < coeffsRef.size(); ++i)
1218 {
1219 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
1220 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
1221 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1222 }
1223}
1224
1225BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_StdMat_VariableP)
1226{
1228 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1230 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1232 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1234 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1236 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1238 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1240 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1242 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1243
1244 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1245 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1247 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1248 v6.get(), v7.get(), segVec, faceVec);
1249
1250 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1252 Nektar::LibUtilities::BasisType basisTypeDir1 =
1254 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1255 quadPointsTypeDir1);
1256 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1257 quadPointsTypeDir1);
1258 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
1259 quadPointsTypeDir1);
1260 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1261 quadPointsKeyDir1);
1262 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1263 quadPointsKeyDir2);
1264 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1265 quadPointsKeyDir3);
1266
1269 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1270
1271 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1272 CollExp.push_back(Exp);
1273
1275 Collections::CollectionOptimisation colOpt(dummySession, 3,
1277 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1278 Collections::Collection c(CollExp, impTypes);
1279 c.Initialise(Collections::eIProductWRTBase);
1280
1281 const int nq = Exp->GetTotPoints();
1282 Array<OneD, NekDouble> phys(nq);
1283 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1284 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1285
1286 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1287
1288 Exp->GetCoords(xc, yc, zc);
1289
1290 for (int i = 0; i < nq; ++i)
1291 {
1292 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1293 }
1294
1295 Exp->IProductWRTBase(phys, coeffs1);
1296 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1297 c.Initialise(Collections::eIProductWRTBase);
1298
1299 double epsilon = 1.0e-8;
1300 for (int i = 0; i < coeffs1.size(); ++i)
1301 {
1302 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1303 }
1304}
1305
1306BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_NoCollection_VariableP)
1307{
1309 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1311 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1313 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1315 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1317 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1319 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1321 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1323 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1324
1325 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1326 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1328 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1329 v6.get(), v7.get(), segVec, faceVec);
1330
1331 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1333 Nektar::LibUtilities::BasisType basisTypeDir1 =
1335 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1336 quadPointsTypeDir1);
1337 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1338 quadPointsTypeDir1);
1339 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
1340 quadPointsTypeDir1);
1341 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1342 quadPointsKeyDir1);
1343 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1344 quadPointsKeyDir2);
1345 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1346 quadPointsKeyDir3);
1347
1350 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1351
1352 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1353 CollExp.push_back(Exp);
1354
1356 Collections::CollectionOptimisation colOpt(dummySession, 3,
1358 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1359 Collections::Collection c(CollExp, impTypes);
1360 c.Initialise(Collections::eIProductWRTBase);
1361
1362 const int nq = Exp->GetTotPoints();
1363 Array<OneD, NekDouble> phys(nq);
1364 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1365 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1366
1367 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1368
1369 Exp->GetCoords(xc, yc, zc);
1370
1371 for (int i = 0; i < nq; ++i)
1372 {
1373 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1374 }
1375
1376 Exp->IProductWRTBase(phys, coeffs1);
1377 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1378 c.Initialise(Collections::eIProductWRTBase);
1379
1380 double epsilon = 1.0e-8;
1381 for (int i = 0; i < coeffs1.size(); ++i)
1382 {
1383 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1384 }
1385}
1386
1387BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_VariableP_CollAll)
1388{
1390 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1392 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1394 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1396 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1398 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1400 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1402 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1404 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1405
1406 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1407 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1409 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1410 v6.get(), v7.get(), segVec, faceVec);
1411
1412 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1414 Nektar::LibUtilities::BasisType basisTypeDir1 =
1416 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(4,
1417 quadPointsTypeDir1);
1418 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
1419 quadPointsTypeDir1);
1420 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
1421 quadPointsTypeDir1);
1422 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1423 quadPointsKeyDir1);
1424 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1425 quadPointsKeyDir2);
1426 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1427 quadPointsKeyDir3);
1428
1431 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1432
1433 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1434 CollExp.push_back(Exp);
1435
1437 Collections::CollectionOptimisation colOpt(dummySession, 3,
1439 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1440 Collections::Collection c(CollExp, impTypes);
1441 c.Initialise(Collections::eIProductWRTBase);
1442
1443 const int nq = Exp->GetTotPoints();
1444 Array<OneD, NekDouble> phys(nq);
1445 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1446 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1447
1448 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1449
1450 Exp->GetCoords(xc, yc, zc);
1451
1452 for (int i = 0; i < nq; ++i)
1453 {
1454 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1455 }
1456
1457 Exp->IProductWRTBase(phys, coeffs1);
1458 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1459
1460 double epsilon = 1.0e-8;
1461 for (int i = 0; i < coeffs1.size(); ++i)
1462 {
1463 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1464 }
1465}
1466
1467BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_VariableP_CollDir02)
1468{
1470 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1472 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1474 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1476 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1478 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1480 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1482 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1484 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1485
1486 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1487 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1489 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1490 v6.get(), v7.get(), segVec, faceVec);
1491
1492 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1494 Nektar::LibUtilities::BasisType basisTypeDir1 =
1496 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(4,
1497 quadPointsTypeDir1);
1498 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1499 quadPointsTypeDir1);
1500 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
1501 quadPointsTypeDir1);
1502 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1503 quadPointsKeyDir1);
1504 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1505 quadPointsKeyDir2);
1506 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1507 quadPointsKeyDir3);
1508
1511 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1512
1513 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1514 CollExp.push_back(Exp);
1515
1517 Collections::CollectionOptimisation colOpt(dummySession, 3,
1519 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1520 Collections::Collection c(CollExp, impTypes);
1521 c.Initialise(Collections::eIProductWRTBase);
1522
1523 const int nq = Exp->GetTotPoints();
1524 Array<OneD, NekDouble> phys(nq);
1525 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1526 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1527
1528 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1529
1530 Exp->GetCoords(xc, yc, zc);
1531
1532 for (int i = 0; i < nq; ++i)
1533 {
1534 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1535 }
1536
1537 Exp->IProductWRTBase(phys, coeffs1);
1538 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1539
1540 double epsilon = 1.0e-8;
1541 for (int i = 0; i < coeffs1.size(); ++i)
1542 {
1543 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1544 }
1545}
1546
1547BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_VariableP_CollDir12)
1548{
1550 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1552 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1554 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1556 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1558 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1560 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1562 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1564 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1565
1566 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1567 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1569 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1570 v6.get(), v7.get(), segVec, faceVec);
1571
1572 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1574 Nektar::LibUtilities::BasisType basisTypeDir1 =
1576 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1577 quadPointsTypeDir1);
1578 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
1579 quadPointsTypeDir1);
1580 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
1581 quadPointsTypeDir1);
1582 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1583 quadPointsKeyDir1);
1584 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1585 quadPointsKeyDir2);
1586 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1587 quadPointsKeyDir3);
1588
1591 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1592
1593 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1594 CollExp.push_back(Exp);
1595
1597 Collections::CollectionOptimisation colOpt(dummySession, 3,
1599 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1600 Collections::Collection c(CollExp, impTypes);
1601 c.Initialise(Collections::eIProductWRTBase);
1602
1603 const int nq = Exp->GetTotPoints();
1604 Array<OneD, NekDouble> phys(nq);
1605 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1606 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1607
1608 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1609
1610 Exp->GetCoords(xc, yc, zc);
1611
1612 for (int i = 0; i < nq; ++i)
1613 {
1614 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1615 }
1616
1617 Exp->IProductWRTBase(phys, coeffs1);
1618 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1619
1620 double epsilon = 1.0e-8;
1621 for (int i = 0; i < coeffs1.size(); ++i)
1622 {
1623 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1624 }
1625}
1626
1627BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_StdMat_VariableP_MultiElmt)
1628{
1630 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1632 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1634 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1636 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1638 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1640 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1642 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1644 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1645
1646 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1647 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1649 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1650 v6.get(), v7.get(), segVec, faceVec);
1651
1652 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1654 Nektar::LibUtilities::BasisType basisTypeDir1 =
1656 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1657 quadPointsTypeDir1);
1658 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1659 quadPointsTypeDir1);
1660 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
1661 quadPointsTypeDir1);
1662 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1663 quadPointsKeyDir1);
1664 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1665 quadPointsKeyDir2);
1666 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1667 quadPointsKeyDir3);
1668
1671 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1672
1673 int nelmts = 10;
1674
1675 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1676 for (int i = 0; i < nelmts; ++i)
1677 {
1678 CollExp.push_back(Exp);
1679 }
1680
1682 Collections::CollectionOptimisation colOpt(dummySession, 3,
1684 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1685 Collections::Collection c(CollExp, impTypes);
1686 c.Initialise(Collections::eIProductWRTBase);
1687
1688 const int nq = Exp->GetTotPoints();
1689 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
1690 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
1691 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
1692
1693 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1694
1695 Exp->GetCoords(xc, yc, zc);
1696
1697 for (int i = 0; i < nq; ++i)
1698 {
1699 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1700 }
1701 Exp->IProductWRTBase(phys, coeffs1);
1702
1703 for (int i = 1; i < nelmts; ++i)
1704 {
1705 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1706 Exp->IProductWRTBase(phys + i * nq,
1707 tmp = coeffs1 + i * Exp->GetNcoeffs());
1708 }
1709
1710 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1711
1712 double epsilon = 1.0e-8;
1713 for (int i = 0; i < coeffs1.size(); ++i)
1714 {
1715 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1716 }
1717}
1718
1719BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_IterPerExp_VariableP_MultiElmt)
1720{
1722 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1724 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1726 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1728 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1730 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1732 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1734 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1736 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1737
1738 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1739 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1741 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1742 v6.get(), v7.get(), segVec, faceVec);
1743
1744 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1746 Nektar::LibUtilities::BasisType basisTypeDir1 =
1748 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1749 quadPointsTypeDir1);
1750 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1751 quadPointsTypeDir1);
1752 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
1753 quadPointsTypeDir1);
1754 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1755 quadPointsKeyDir1);
1756 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1757 quadPointsKeyDir2);
1758 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1759 quadPointsKeyDir3);
1760
1763 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1764
1765 int nelmts = 10;
1766
1767 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1768 for (int i = 0; i < nelmts; ++i)
1769 {
1770 CollExp.push_back(Exp);
1771 }
1772
1774 Collections::CollectionOptimisation colOpt(dummySession, 3,
1776 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1777 Collections::Collection c(CollExp, impTypes);
1778 c.Initialise(Collections::eIProductWRTBase);
1779
1780 const int nq = Exp->GetTotPoints();
1781 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
1782 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
1783 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
1784
1785 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1786
1787 Exp->GetCoords(xc, yc, zc);
1788
1789 for (int i = 0; i < nq; ++i)
1790 {
1791 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1792 }
1793 Exp->IProductWRTBase(phys, coeffs1);
1794
1795 for (int i = 1; i < nelmts; ++i)
1796 {
1797 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1798 Exp->IProductWRTBase(phys + i * nq,
1799 tmp = coeffs1 + i * Exp->GetNcoeffs());
1800 }
1801
1802 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1803
1804 double epsilon = 1.0e-8;
1805 for (int i = 0; i < coeffs1.size(); ++i)
1806 {
1807 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1808 }
1809}
1810
1811BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_UniformP)
1812{
1814 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1816 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1818 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1820 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1822 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1824 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1826 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1828 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1829
1830 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1831 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1833 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1834 v6.get(), v7.get(), segVec, faceVec);
1835
1836 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1838 Nektar::LibUtilities::BasisType basisTypeDir1 =
1840 unsigned int numQuadPoints = 6;
1841 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
1842 quadPointsTypeDir1);
1843 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1844 quadPointsKeyDir1);
1845
1848 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
1849
1850 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1851 CollExp.push_back(Exp);
1852
1854 Collections::CollectionOptimisation colOpt(dummySession, 3,
1856 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1857 Collections::Collection c(CollExp, impTypes);
1858 c.Initialise(Collections::eIProductWRTBase);
1859
1860 const int nq = Exp->GetTotPoints();
1861 Array<OneD, NekDouble> phys(nq);
1862 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1863 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1864
1865 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1866
1867 Exp->GetCoords(xc, yc, zc);
1868
1869 for (int i = 0; i < nq; ++i)
1870 {
1871 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1872 }
1873
1874 Exp->IProductWRTBase(phys, coeffs1);
1875 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1876
1877 double epsilon = 1.0e-6;
1878 for (int i = 0; i < coeffs1.size(); ++i)
1879 {
1880 // clamp values below 1e-16 to zero
1881 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-16) ? 0.0 : coeffs1[i];
1882 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-16) ? 0.0 : coeffs2[i];
1883 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1884 }
1885}
1886
1887BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_VariableP)
1888{
1890 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1892 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1894 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1896 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1898 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1900 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1902 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1904 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1905
1906 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1907 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1909 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1910 v6.get(), v7.get(), segVec, faceVec);
1911
1912 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1914 Nektar::LibUtilities::BasisType basisTypeDir1 =
1916 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
1917 quadPointsTypeDir1);
1918 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
1919 quadPointsTypeDir1);
1920 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
1921 quadPointsTypeDir1);
1922 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1923 quadPointsKeyDir1);
1924 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
1925 quadPointsKeyDir2);
1926 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
1927 quadPointsKeyDir3);
1928
1931 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
1932
1933 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1934 CollExp.push_back(Exp);
1935
1937 Collections::CollectionOptimisation colOpt(dummySession, 3,
1939 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1940 Collections::Collection c(CollExp, impTypes);
1941 c.Initialise(Collections::eIProductWRTBase);
1942
1943 const int nq = Exp->GetTotPoints();
1944 Array<OneD, NekDouble> phys(nq);
1945 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1946 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1947
1948 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1949
1950 Exp->GetCoords(xc, yc, zc);
1951
1952 for (int i = 0; i < nq; ++i)
1953 {
1954 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1955 }
1956
1957 Exp->IProductWRTBase(phys, coeffs1);
1958 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1959
1960 double epsilon = 1.0e-6;
1961 for (int i = 0; i < coeffs1.size(); ++i)
1962 {
1963 // clamp values below 1e-16 to zero
1964 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-16) ? 0.0 : coeffs1[i];
1965 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-16) ? 0.0 : coeffs2[i];
1966 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1967 }
1968}
1969
1970BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_UniformP_MultiElmt)
1971{
1973 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
1975 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1977 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1979 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1981 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1983 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
1985 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
1987 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
1988
1989 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
1990 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
1992 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
1993 v6.get(), v7.get(), segVec, faceVec);
1994
1995 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
1997 Nektar::LibUtilities::BasisType basisTypeDir1 =
1999 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2000 quadPointsTypeDir1);
2001 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2002 quadPointsKeyDir1);
2003
2006 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2007
2008 int nelmts = 10;
2009
2010 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2011 for (int i = 0; i < nelmts; ++i)
2012 {
2013 CollExp.push_back(Exp);
2014 }
2015
2017 Collections::CollectionOptimisation colOpt(dummySession, 3,
2019 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2020 Collections::Collection c(CollExp, impTypes);
2021 c.Initialise(Collections::eIProductWRTBase);
2022
2023 const int nq = Exp->GetTotPoints();
2024 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
2025 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
2026 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
2027
2028 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2029
2030 Exp->GetCoords(xc, yc, zc);
2031
2032 for (int i = 0; i < nq; ++i)
2033 {
2034 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2035 }
2036 Exp->IProductWRTBase(phys, coeffs1);
2037
2038 for (int i = 1; i < nelmts; ++i)
2039 {
2040 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
2041 Exp->IProductWRTBase(phys + i * nq,
2042 tmp = coeffs1 + i * Exp->GetNcoeffs());
2043 }
2044 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
2045
2046 double epsilon = 1.0e-6;
2047 for (int i = 0; i < coeffs1.size(); ++i)
2048 {
2049 // clamp values below 1e-16 to zero
2050 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-16) ? 0.0 : coeffs1[i];
2051 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-16) ? 0.0 : coeffs2[i];
2052 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2053 }
2054}
2055
2056BOOST_AUTO_TEST_CASE(TestHexIProductWRTBase_SumFac_VariableP_MultiElmt)
2057{
2059 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2061 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2063 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2065 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2067 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2069 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2071 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2073 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2074
2075 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2076 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2078 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2079 v6.get(), v7.get(), segVec, faceVec);
2080
2081 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2083 Nektar::LibUtilities::BasisType basisTypeDir1 =
2085 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2086 quadPointsTypeDir1);
2087 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(7,
2088 quadPointsTypeDir1);
2089 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(9,
2090 quadPointsTypeDir1);
2091 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2092 quadPointsKeyDir1);
2093 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2094 quadPointsKeyDir2);
2095 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2096 quadPointsKeyDir3);
2097
2100 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2101
2102 int nelmts = 10;
2103
2104 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2105 for (int i = 0; i < nelmts; ++i)
2106 {
2107 CollExp.push_back(Exp);
2108 }
2109
2111 Collections::CollectionOptimisation colOpt(dummySession, 3,
2113 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2114 Collections::Collection c(CollExp, impTypes);
2115 c.Initialise(Collections::eIProductWRTBase);
2116
2117 const int nq = Exp->GetTotPoints();
2118 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
2119 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
2120 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
2121
2122 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2123
2124 Exp->GetCoords(xc, yc, zc);
2125
2126 for (int i = 0; i < nq; ++i)
2127 {
2128 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2129 }
2130 Exp->IProductWRTBase(phys, coeffs1);
2131
2132 for (int i = 1; i < nelmts; ++i)
2133 {
2134 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
2135 Exp->IProductWRTBase(phys + i * nq,
2136 tmp = coeffs1 + i * Exp->GetNcoeffs());
2137 }
2138 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
2139
2140 double epsilon = 1.0e-4;
2141 for (int i = 0; i < coeffs1.size(); ++i)
2142 {
2143 // clamp values below 1e-14 to zero
2144 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2145 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2146 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2147 }
2148}
2149
2151 TestHexIProductWRTBase_SumFac_VariableP_MultiElmt_CollDir02)
2152{
2154 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2156 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2158 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2160 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2162 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2164 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2166 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2168 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2169
2170 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2171 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2173 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2174 v6.get(), v7.get(), segVec, faceVec);
2175
2176 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2178 Nektar::LibUtilities::BasisType basisTypeDir1 =
2180 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(4,
2181 quadPointsTypeDir1);
2182 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
2183 quadPointsTypeDir1);
2184 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
2185 quadPointsTypeDir1);
2186 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2187 quadPointsKeyDir1);
2188 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2189 quadPointsKeyDir2);
2190 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2191 quadPointsKeyDir3);
2192
2195 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2196
2197 int nelmts = 10;
2198
2199 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2200 for (int i = 0; i < nelmts; ++i)
2201 {
2202 CollExp.push_back(Exp);
2203 }
2204
2206 Collections::CollectionOptimisation colOpt(dummySession, 3,
2208 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2209 Collections::Collection c(CollExp, impTypes);
2210 c.Initialise(Collections::eIProductWRTBase);
2211
2212 const int nq = Exp->GetTotPoints();
2213 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
2214 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
2215 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
2216
2217 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2218
2219 Exp->GetCoords(xc, yc, zc);
2220
2221 for (int i = 0; i < nq; ++i)
2222 {
2223 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2224 }
2225 Exp->IProductWRTBase(phys, coeffs1);
2226
2227 for (int i = 1; i < nelmts; ++i)
2228 {
2229 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
2230 Exp->IProductWRTBase(phys + i * nq,
2231 tmp = coeffs1 + i * Exp->GetNcoeffs());
2232 }
2233 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
2234
2235 double epsilon = 1.0e-4;
2236 for (int i = 0; i < coeffs1.size(); ++i)
2237 {
2238 // clamp values below 1e-14 to zero
2239 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2240 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2241 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2242 }
2243}
2244
2246 TestHexIProductWRTBase_SumFac_VariableP_MultiElmt_CollDir12)
2247{
2249 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2251 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2253 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2255 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2257 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2259 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2261 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2263 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2264
2265 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2266 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2268 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2269 v6.get(), v7.get(), segVec, faceVec);
2270
2271 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2273 Nektar::LibUtilities::BasisType basisTypeDir1 =
2275 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2276 quadPointsTypeDir1);
2277 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
2278 quadPointsTypeDir1);
2279 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
2280 quadPointsTypeDir1);
2281 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2282 quadPointsKeyDir1);
2283 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2284 quadPointsKeyDir2);
2285 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2286 quadPointsKeyDir3);
2287
2290 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2291
2292 int nelmts = 10;
2293
2294 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2295 for (int i = 0; i < nelmts; ++i)
2296 {
2297 CollExp.push_back(Exp);
2298 }
2299
2301 Collections::CollectionOptimisation colOpt(dummySession, 3,
2303 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2304 Collections::Collection c(CollExp, impTypes);
2305 c.Initialise(Collections::eIProductWRTBase);
2306
2307 const int nq = Exp->GetTotPoints();
2308 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
2309 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
2310 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
2311
2312 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2313
2314 Exp->GetCoords(xc, yc, zc);
2315
2316 for (int i = 0; i < nq; ++i)
2317 {
2318 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2319 }
2320 Exp->IProductWRTBase(phys, coeffs1);
2321
2322 for (int i = 1; i < nelmts; ++i)
2323 {
2324 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
2325 Exp->IProductWRTBase(phys + i * nq,
2326 tmp = coeffs1 + i * Exp->GetNcoeffs());
2327 }
2328 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
2329
2330 double epsilon = 1.0e-4;
2331 for (int i = 0; i < coeffs1.size(); ++i)
2332 {
2333 // clamp values below 1e-14 to zero
2334 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2335 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2336 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2337 }
2338}
2339
2340BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_IterPerExp_UniformP)
2341{
2343 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2345 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2347 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2349 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2351 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2353 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2355 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2357 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2358
2359 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2360 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2362 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2363 v6.get(), v7.get(), segVec, faceVec);
2364
2365 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2367 Nektar::LibUtilities::BasisType basisTypeDir1 =
2369 unsigned int numQuadPoints = 6;
2370 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
2371 quadPointsTypeDir1);
2372 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2373 quadPointsKeyDir1);
2374
2377 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2378
2379 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2380 CollExp.push_back(Exp);
2381
2383 Collections::CollectionOptimisation colOpt(dummySession, 3,
2385 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2386 Collections::Collection c(CollExp, impTypes);
2387 c.Initialise(Collections::ePhysDeriv);
2388
2389 const int nq = Exp->GetTotPoints();
2390 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2391 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2392 Array<OneD, NekDouble> diff1(3 * nq);
2393 Array<OneD, NekDouble> diff2(3 * nq);
2394
2395 Exp->GetCoords(xc, yc, zc);
2396
2397 for (int i = 0; i < nq; ++i)
2398 {
2399 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2400 }
2401
2402 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
2403 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
2404 tmp2 = diff2 + 2 * nq);
2405
2406 double epsilon = 1.0e-8;
2407 for (int i = 0; i < diff1.size(); ++i)
2408 {
2409 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2410 }
2411}
2412
2413BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_MatrixFree_UniformP_Undeformed)
2414{
2416 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2418 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2420 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2422 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2424 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2426 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2428 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2430 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2431
2432 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2433 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2435 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2436 v6.get(), v7.get(), segVec, faceVec);
2437
2438 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2440 Nektar::LibUtilities::BasisType basisTypeDir1 =
2442 unsigned int numQuadPoints = 6;
2443 unsigned int numModes = 2;
2444 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
2445 quadPointsTypeDir1);
2446 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2447 quadPointsKeyDir1);
2448
2451 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2452
2453 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2454 CollExp.push_back(Exp);
2455
2457 Collections::CollectionOptimisation colOpt(dummySession, 2,
2459 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2460 Collections::Collection c(CollExp, impTypes);
2461 c.Initialise(Collections::ePhysDeriv);
2462
2463 const int nq = Exp->GetTotPoints();
2464 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2465 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2466 Array<OneD, NekDouble> diffRef(3 * nq);
2467 Array<OneD, NekDouble> diff(3 * nq);
2468
2469 Exp->GetCoords(xc, yc, zc);
2470
2471 for (int i = 0; i < nq; ++i)
2472 {
2473 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2474 }
2475
2476 Exp->PhysDeriv(phys, diffRef, tmp = diffRef + nq, tmp1 = diffRef + 2 * nq);
2477 c.ApplyOperator(Collections::ePhysDeriv, phys, diff, tmp = diff + nq,
2478 tmp2 = diff + 2 * nq);
2479
2480 double epsilon = 1.0e-8;
2481 for (int i = 0; i < diffRef.size(); ++i)
2482 {
2483 diffRef[i] = (std::abs(diffRef[i]) < 1e-14) ? 0.0 : diffRef[i];
2484 diff[i] = (std::abs(diff[i]) < 1e-14) ? 0.0 : diff[i];
2485 BOOST_CHECK_CLOSE(diffRef[i], diff[i], epsilon);
2486 }
2487}
2488
2489BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_MatrixFree_UniformP_Deformed)
2490{
2492 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2494 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2496 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2498 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2500 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2502 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2504 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
2506 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2507
2508 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2509 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2511 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2512 v6.get(), v7.get(), segVec, faceVec);
2513
2514 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2516 Nektar::LibUtilities::BasisType basisTypeDir1 =
2518 unsigned int numQuadPoints = 5;
2519 unsigned int numModes = 2;
2520 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
2521 quadPointsTypeDir1);
2522 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2523 quadPointsKeyDir1);
2524
2527 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2528
2529 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2530 CollExp.push_back(Exp);
2531
2533 Collections::CollectionOptimisation colOpt(dummySession, 2,
2535 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2536 Collections::Collection c(CollExp, impTypes);
2537 c.Initialise(Collections::ePhysDeriv);
2538
2539 const int nq = Exp->GetTotPoints();
2540 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2541 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2542 Array<OneD, NekDouble> diffRef(3 * nq);
2543 Array<OneD, NekDouble> diff(3 * nq);
2544
2545 Exp->GetCoords(xc, yc, zc);
2546
2547 for (int i = 0; i < nq; ++i)
2548 {
2549 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2550 }
2551
2552 Exp->PhysDeriv(phys, diffRef, tmp = diffRef + nq, tmp1 = diffRef + 2 * nq);
2553 c.ApplyOperator(Collections::ePhysDeriv, phys, diff, tmp = diff + nq,
2554 tmp2 = diff + 2 * nq);
2555
2556 double epsilon = 1.0e-8;
2557 for (int i = 0; i < diffRef.size(); ++i)
2558 {
2559 diffRef[i] = (std::abs(diffRef[i]) < 1e-14) ? 0.0 : diffRef[i];
2560 diff[i] = (std::abs(diff[i]) < 1e-14) ? 0.0 : diff[i];
2561 BOOST_CHECK_CLOSE(diffRef[i], diff[i], epsilon);
2562 }
2563}
2564
2565BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_IterPerExp_VariableP_MultiElmt)
2566{
2568 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2570 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2572 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2574 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2576 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2578 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2580 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2582 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2583
2584 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2585 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2587 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2588 v6.get(), v7.get(), segVec, faceVec);
2589
2590 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2592 Nektar::LibUtilities::BasisType basisTypeDir1 =
2594 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2595 quadPointsTypeDir1);
2596 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
2597 quadPointsTypeDir1);
2598 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
2599 quadPointsTypeDir1);
2600 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2601 quadPointsKeyDir1);
2602 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2603 quadPointsKeyDir2);
2604 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2605 quadPointsKeyDir3);
2606
2609 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2610
2611 int nelmts = 10;
2612
2613 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2614 for (int i = 0; i < nelmts; ++i)
2615 {
2616 CollExp.push_back(Exp);
2617 }
2618
2620 Collections::CollectionOptimisation colOpt(dummySession, 3,
2622 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2623 Collections::Collection c(CollExp, impTypes);
2624 c.Initialise(Collections::ePhysDeriv);
2625
2626 const int nq = Exp->GetTotPoints();
2627 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2628 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
2629 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
2630 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
2631
2632 Exp->GetCoords(xc, yc, zc);
2633
2634 for (int i = 0; i < nq; ++i)
2635 {
2636 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2637 }
2638 Exp->PhysDeriv(phys, tmp = diff1, tmp1 = diff1 + (nelmts)*nq,
2639 tmp2 = diff1 + (2 * nelmts) * nq);
2640 for (int i = 1; i < nelmts; ++i)
2641 {
2642 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
2643 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
2644 tmp1 = diff1 + (nelmts + i) * nq,
2645 tmp2 = diff1 + (2 * nelmts + i) * nq);
2646 }
2647
2648 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
2649 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
2650
2651 double epsilon = 1.0e-8;
2652 for (int i = 0; i < diff1.size(); ++i)
2653 {
2654 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2655 }
2656}
2657
2658BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_NoCollection_VariableP_MultiElmt)
2659{
2661 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2663 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2665 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2667 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2669 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2671 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2673 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2675 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2676
2677 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2678 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2680 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2681 v6.get(), v7.get(), segVec, faceVec);
2682
2683 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2685 Nektar::LibUtilities::BasisType basisTypeDir1 =
2687 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2688 quadPointsTypeDir1);
2689 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
2690 quadPointsTypeDir1);
2691 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
2692 quadPointsTypeDir1);
2693 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2694 quadPointsKeyDir1);
2695 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2696 quadPointsKeyDir2);
2697 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2698 quadPointsKeyDir3);
2699
2702 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2703
2704 int nelmts = 10;
2705
2706 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2707 for (int i = 0; i < nelmts; ++i)
2708 {
2709 CollExp.push_back(Exp);
2710 }
2711
2713 Collections::CollectionOptimisation colOpt(dummySession, 3,
2715 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2716 Collections::Collection c(CollExp, impTypes);
2717 c.Initialise(Collections::ePhysDeriv);
2718
2719 const int nq = Exp->GetTotPoints();
2720 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2721 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
2722 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
2723 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
2724
2725 Exp->GetCoords(xc, yc, zc);
2726
2727 for (int i = 0; i < nq; ++i)
2728 {
2729 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2730 }
2731 Exp->PhysDeriv(phys, tmp = diff1, tmp1 = diff1 + (nelmts)*nq,
2732 tmp2 = diff1 + (2 * nelmts) * nq);
2733
2734 for (int i = 1; i < nelmts; ++i)
2735 {
2736 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
2737 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
2738 tmp1 = diff1 + (nelmts + i) * nq,
2739 tmp2 = diff1 + (2 * nelmts + i) * nq);
2740 }
2741
2742 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
2743 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
2744
2745 double epsilon = 1.0e-8;
2746 for (int i = 0; i < diff1.size(); ++i)
2747 {
2748 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2749 }
2750}
2751
2752BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_StdMat_UniformP)
2753{
2755 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2757 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2759 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2761 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2763 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2765 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2767 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2769 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2770
2771 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2772 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2774 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2775 v6.get(), v7.get(), segVec, faceVec);
2776
2777 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2779 Nektar::LibUtilities::BasisType basisTypeDir1 =
2781 unsigned int numQuadPoints = 4;
2782 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
2783 quadPointsTypeDir1);
2784 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 3,
2785 quadPointsKeyDir1);
2786
2789 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2790
2791 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2792 CollExp.push_back(Exp);
2793
2795 Collections::CollectionOptimisation colOpt(dummySession, 3,
2797 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2798 Collections::Collection c(CollExp, impTypes);
2799 c.Initialise(Collections::ePhysDeriv);
2800
2801 const int nq = Exp->GetTotPoints();
2802 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2803 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2804 Array<OneD, NekDouble> diff1(3 * nq);
2805 Array<OneD, NekDouble> diff2(3 * nq);
2806
2807 Exp->GetCoords(xc, yc, zc);
2808
2809 for (int i = 0; i < nq; ++i)
2810 {
2811 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2812 }
2813
2814 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
2815 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
2816 tmp2 = diff2 + 2 * nq);
2817
2818 double epsilon = 1.0e-8;
2819 for (int i = 0; i < diff1.size(); ++i)
2820 {
2821 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2822 }
2823}
2824
2825BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_StdMat_VariableP_MultiElmt)
2826{
2828 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2830 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2832 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2834 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2836 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2838 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2840 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2842 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2843
2844 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2845 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2847 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2848 v6.get(), v7.get(), segVec, faceVec);
2849
2850 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2852 Nektar::LibUtilities::BasisType basisTypeDir1 =
2854 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
2855 quadPointsTypeDir1);
2856 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
2857 quadPointsTypeDir1);
2858 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
2859 quadPointsTypeDir1);
2860 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2861 quadPointsKeyDir1);
2862 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
2863 quadPointsKeyDir2);
2864 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
2865 quadPointsKeyDir3);
2866
2869 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
2870
2871 int nelmts = 10;
2872
2873 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2874 for (int i = 0; i < nelmts; ++i)
2875 {
2876 CollExp.push_back(Exp);
2877 }
2878
2880 Collections::CollectionOptimisation colOpt(dummySession, 3,
2882 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2883 Collections::Collection c(CollExp, impTypes);
2884 c.Initialise(Collections::ePhysDeriv);
2885
2886 const int nq = Exp->GetTotPoints();
2887 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2888 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
2889 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
2890 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
2891
2892 Exp->GetCoords(xc, yc, zc);
2893
2894 for (int i = 0; i < nq; ++i)
2895 {
2896 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2897 }
2898 Exp->PhysDeriv(phys, tmp = diff1, tmp1 = diff1 + (nelmts)*nq,
2899 tmp2 = diff1 + (2 * nelmts) * nq);
2900 for (int i = 1; i < nelmts; ++i)
2901 {
2902 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
2903 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
2904 tmp1 = diff1 + (nelmts + i) * nq,
2905 tmp2 = diff1 + (2 * nelmts + i) * nq);
2906 }
2907
2908 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
2909 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
2910
2911 double epsilon = 1.0e-8;
2912 for (int i = 0; i < diff1.size(); ++i)
2913 {
2914 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2915 }
2916}
2917
2918BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_SumFac_UniformP)
2919{
2921 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2923 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2925 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2927 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2929 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2931 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
2933 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
2935 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
2936
2937 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
2938 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
2940 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
2941 v6.get(), v7.get(), segVec, faceVec);
2942
2943 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
2945 Nektar::LibUtilities::BasisType basisTypeDir1 =
2947 unsigned int numQuadPoints = 6;
2948 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
2949 quadPointsTypeDir1);
2950 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2951 quadPointsKeyDir1);
2952
2955 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
2956
2957 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2958 CollExp.push_back(Exp);
2959
2961 Collections::CollectionOptimisation colOpt(dummySession, 3,
2963 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2964 Collections::Collection c(CollExp, impTypes);
2965 c.Initialise(Collections::ePhysDeriv);
2966
2967 const int nq = Exp->GetTotPoints();
2968 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2969 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2970 Array<OneD, NekDouble> diff1(3 * nq);
2971 Array<OneD, NekDouble> diff2(3 * nq);
2972
2973 Exp->GetCoords(xc, yc, zc);
2974
2975 for (int i = 0; i < nq; ++i)
2976 {
2977 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2978 }
2979
2980 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
2981 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
2982 tmp2 = diff2 + 2 * nq);
2983
2984 double epsilon = 1.0e-8;
2985 for (int i = 0; i < diff1.size(); ++i)
2986 {
2987 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2988 }
2989}
2990
2991BOOST_AUTO_TEST_CASE(TestHexPhysDeriv_SumFac_VariableP_MultiElmt)
2992{
2994 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2996 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2998 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3000 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3002 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3004 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3006 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3008 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3009
3010 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3011 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3013 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3014 v6.get(), v7.get(), segVec, faceVec);
3015
3016 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3018 Nektar::LibUtilities::BasisType basisTypeDir1 =
3020 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
3021 quadPointsTypeDir1);
3022 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
3023 quadPointsTypeDir1);
3024 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
3025 quadPointsTypeDir1);
3026 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3027 quadPointsKeyDir1);
3028 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
3029 quadPointsKeyDir2);
3030 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
3031 quadPointsKeyDir3);
3032
3035 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
3036
3037 int nelmts = 10;
3038
3039 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3040 for (int i = 0; i < nelmts; ++i)
3041 {
3042 CollExp.push_back(Exp);
3043 }
3044
3046 Collections::CollectionOptimisation colOpt(dummySession, 3,
3048 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3049 Collections::Collection c(CollExp, impTypes);
3050 c.Initialise(Collections::ePhysDeriv);
3051
3052 const int nq = Exp->GetTotPoints();
3053 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3054 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
3055 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
3056 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
3057
3058 Exp->GetCoords(xc, yc, zc);
3059
3060 for (int i = 0; i < nq; ++i)
3061 {
3062 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3063 }
3064 Exp->PhysDeriv(phys, tmp = diff1, tmp1 = diff1 + (nelmts)*nq,
3065 tmp2 = diff1 + (2 * nelmts) * nq);
3066 for (int i = 1; i < nelmts; ++i)
3067 {
3068 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
3069 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
3070 tmp1 = diff1 + (nelmts + i) * nq,
3071 tmp2 = diff1 + (2 * nelmts + i) * nq);
3072 }
3073
3074 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
3075 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
3076
3077 double epsilon = 1.0e-8;
3078 for (int i = 0; i < diff1.size(); ++i)
3079 {
3080 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
3081 }
3082}
3083
3084BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_Iterperexp_UniformP)
3085{
3087 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3089 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3091 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3093 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3095 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3097 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3099 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3101 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3102
3103 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3104 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3106 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3107 v6.get(), v7.get(), segVec, faceVec);
3108
3109 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3111 Nektar::LibUtilities::BasisType basisTypeDir1 =
3113 unsigned int numQuadPoints = 6;
3114 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3115 quadPointsTypeDir1);
3116 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3117 quadPointsKeyDir1);
3118
3121 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3122
3123 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3124 CollExp.push_back(Exp);
3125
3127 Collections::CollectionOptimisation colOpt(dummySession, 3,
3129 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3130 Collections::Collection c(CollExp, impTypes);
3132
3133 const int nq = Exp->GetTotPoints();
3134 const int nm = Exp->GetNcoeffs();
3135 Array<OneD, NekDouble> phys1(nq, 0.0);
3136 Array<OneD, NekDouble> phys2(nq, 0.0);
3137 Array<OneD, NekDouble> phys3(nq, 0.0);
3138 Array<OneD, NekDouble> coeffs1(nm, 0.0);
3139 Array<OneD, NekDouble> coeffs2(nm, 0.0);
3140
3141 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3142
3143 Exp->GetCoords(xc, yc, zc);
3144
3145 for (int i = 0; i < nq; ++i)
3146 {
3147 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3148 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3149 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3150 }
3151
3152 // Standard routines
3153 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
3154 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
3155 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3156 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
3157 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3158
3159 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3160 coeffs2);
3161
3162 double epsilon = 1.0e-8;
3163 for (int i = 0; i < nm; ++i)
3164 {
3165 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3166 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3167 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3168 }
3169}
3170
3171BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_MatrixFree_UniformP_Undeformed)
3172{
3174 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3176 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3178 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3180 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3182 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3184 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3186 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3188 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3189
3190 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3191 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3193 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3194 v6.get(), v7.get(), segVec, faceVec);
3195
3196 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3198 Nektar::LibUtilities::BasisType basisTypeDir1 =
3200 unsigned int numQuadPoints = 5;
3201 unsigned int numModes = 4;
3202 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3203 quadPointsTypeDir1);
3204 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3205 quadPointsKeyDir1);
3206
3209 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3210
3211 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3212 CollExp.push_back(Exp);
3213
3215 Collections::CollectionOptimisation colOpt(dummySession, 2,
3217 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3218 Collections::Collection c(CollExp, impTypes);
3220
3221 const int nq = Exp->GetTotPoints();
3222 const int nm = Exp->GetNcoeffs();
3223 Array<OneD, NekDouble> phys1(nq, 0.0);
3224 Array<OneD, NekDouble> phys2(nq, 0.0);
3225 Array<OneD, NekDouble> phys3(nq, 0.0);
3226 Array<OneD, NekDouble> coeffsRef(nm, 0.0);
3227 Array<OneD, NekDouble> coeffs(nm, 0.0);
3228
3229 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3230
3231 Exp->GetCoords(xc, yc, zc);
3232
3233 for (int i = 0; i < nq; ++i)
3234 {
3235 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3236 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3237 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3238 }
3239
3240 // Standard routines
3241 Exp->IProductWRTDerivBase(0, phys1, coeffsRef);
3242 Exp->IProductWRTDerivBase(1, phys2, coeffs);
3243 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3244 Exp->IProductWRTDerivBase(2, phys3, coeffs);
3245 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3246
3247 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3248 coeffs);
3249
3250 double epsilon = 1.0e-8;
3251 for (int i = 0; i < nm; ++i)
3252 {
3253 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3254 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3255 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3256 }
3257}
3258
3259BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_MatrixFree_UniformP_Deformed)
3260{
3262 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3264 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3266 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3268 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3270 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3272 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3274 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
3276 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3277
3278 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3279 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3281 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3282 v6.get(), v7.get(), segVec, faceVec);
3283
3284 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3286 Nektar::LibUtilities::BasisType basisTypeDir1 =
3288 unsigned int numQuadPoints = 5;
3289 unsigned int numModes = 4;
3290 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3291 quadPointsTypeDir1);
3292 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3293 quadPointsKeyDir1);
3294
3297 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3298
3299 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3300 CollExp.push_back(Exp);
3301
3303 Collections::CollectionOptimisation colOpt(dummySession, 2,
3305 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3306 Collections::Collection c(CollExp, impTypes);
3308
3309 const int nq = Exp->GetTotPoints();
3310 const int nm = Exp->GetNcoeffs();
3311 Array<OneD, NekDouble> phys1(nq, 0.0);
3312 Array<OneD, NekDouble> phys2(nq, 0.0);
3313 Array<OneD, NekDouble> phys3(nq, 0.0);
3314 Array<OneD, NekDouble> coeffsRef(nm, 0.0);
3315 Array<OneD, NekDouble> coeffs(nm, 0.0);
3316
3317 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3318
3319 Exp->GetCoords(xc, yc, zc);
3320
3321 for (int i = 0; i < nq; ++i)
3322 {
3323 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3324 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3325 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3326 }
3327
3328 // Standard routines
3329 Exp->IProductWRTDerivBase(0, phys1, coeffsRef);
3330 Exp->IProductWRTDerivBase(1, phys2, coeffs);
3331 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3332 Exp->IProductWRTDerivBase(2, phys3, coeffs);
3333 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3334
3335 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3336 coeffs);
3337
3338 double epsilon = 1.0e-8;
3339 for (int i = 0; i < nm; ++i)
3340 {
3341 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3342 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3343 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3344 }
3345}
3346
3348 TestHexIProductWRTDerivBase_MatrixFree_UniformP_Deformed_OverInt)
3349{
3351 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3353 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3355 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3357 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3359 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3361 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3363 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
3365 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3366
3367 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3368 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3370 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3371 v6.get(), v7.get(), segVec, faceVec);
3372
3373 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3375 Nektar::LibUtilities::BasisType basisTypeDir1 =
3377 unsigned int numQuadPoints = 8;
3378 unsigned int numModes = 4;
3379 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3380 quadPointsTypeDir1);
3381 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3382 quadPointsKeyDir1);
3383
3386 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3387
3388 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3389 CollExp.push_back(Exp);
3390
3392 Collections::CollectionOptimisation colOpt(dummySession, 2,
3394 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3395 Collections::Collection c(CollExp, impTypes);
3397
3398 const int nq = Exp->GetTotPoints();
3399 const int nm = Exp->GetNcoeffs();
3400 Array<OneD, NekDouble> phys1(nq, 0.0);
3401 Array<OneD, NekDouble> phys2(nq, 0.0);
3402 Array<OneD, NekDouble> phys3(nq, 0.0);
3403 Array<OneD, NekDouble> coeffsRef(nm, 0.0);
3404 Array<OneD, NekDouble> coeffs(nm, 0.0);
3405
3406 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3407
3408 Exp->GetCoords(xc, yc, zc);
3409
3410 for (int i = 0; i < nq; ++i)
3411 {
3412 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3413 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3414 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3415 }
3416
3417 // Standard routines
3418 Exp->IProductWRTDerivBase(0, phys1, coeffsRef);
3419 Exp->IProductWRTDerivBase(1, phys2, coeffs);
3420 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3421 Exp->IProductWRTDerivBase(2, phys3, coeffs);
3422 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3423
3424 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3425 coeffs);
3426
3427 double epsilon = 1.0e-8;
3428 for (int i = 0; i < nm; ++i)
3429 {
3430 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3431 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3432 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3433 }
3434}
3435
3436BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_IterPerExp_VariableP_MultiElmt)
3437{
3439 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3441 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3443 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3445 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3447 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3449 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3451 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3453 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3454
3455 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3456 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3458 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3459 v6.get(), v7.get(), segVec, faceVec);
3460
3461 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3463 Nektar::LibUtilities::BasisType basisTypeDir1 =
3465 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
3466 quadPointsTypeDir1);
3467 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
3468 quadPointsTypeDir1);
3469 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
3470 quadPointsTypeDir1);
3471 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3472 quadPointsKeyDir1);
3473 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
3474 quadPointsKeyDir2);
3475 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
3476 quadPointsKeyDir3);
3477
3480 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
3481
3482 int nelmts = 10;
3483
3484 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3485 for (int i = 0; i < nelmts; ++i)
3486 {
3487 CollExp.push_back(Exp);
3488 }
3489
3491 Collections::CollectionOptimisation colOpt(dummySession, 3,
3493 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3494 Collections::Collection c(CollExp, impTypes);
3496
3497 const int nq = Exp->GetTotPoints();
3498 const int nm = Exp->GetNcoeffs();
3499 Array<OneD, NekDouble> phys1(nelmts * nq, 0.0);
3500 Array<OneD, NekDouble> phys2(nelmts * nq, 0.0);
3501 Array<OneD, NekDouble> phys3(nelmts * nq, 0.0);
3502 Array<OneD, NekDouble> coeffs1(nelmts * nm, 0.0);
3503 Array<OneD, NekDouble> coeffs2(nelmts * nm, 0.0);
3504 Array<OneD, NekDouble> tmp;
3505
3506 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3507
3508 Exp->GetCoords(xc, yc, zc);
3509
3510 for (int i = 0; i < nq; ++i)
3511 {
3512 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3513 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3514 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3515 }
3516
3517 for (int i = 1; i < nelmts; ++i)
3518 {
3519 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
3520 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
3521 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
3522 }
3523
3524 // Standard routines
3525 for (int i = 0; i < nelmts; ++i)
3526 {
3527
3528 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
3529 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
3530 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3531 tmp = coeffs1 + i * nm, 1);
3532 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
3533 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3534 tmp = coeffs1 + i * nm, 1);
3535 }
3536
3537 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3538 coeffs2);
3539
3540 double epsilon = 1.0e-6;
3541 for (int i = 0; i < coeffs1.size(); ++i)
3542 {
3543 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3544 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3545 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3546 }
3547}
3548
3549BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_StdMat_UniformP)
3550{
3552 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3554 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3556 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3558 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3560 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3562 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3564 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3566 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3567
3568 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3569 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3571 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3572 v6.get(), v7.get(), segVec, faceVec);
3573
3574 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3576 Nektar::LibUtilities::BasisType basisTypeDir1 =
3578 unsigned int numQuadPoints = 6;
3579 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3580 quadPointsTypeDir1);
3581 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3582 quadPointsKeyDir1);
3583
3586 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3587
3588 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3589 CollExp.push_back(Exp);
3590
3592 Collections::CollectionOptimisation colOpt(dummySession, 3,
3594 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3595 Collections::Collection c(CollExp, impTypes);
3597
3598 const int nq = Exp->GetTotPoints();
3599 const int nm = Exp->GetNcoeffs();
3600 Array<OneD, NekDouble> phys1(nq, 0.0);
3601 Array<OneD, NekDouble> phys2(nq, 0.0);
3602 Array<OneD, NekDouble> phys3(nq, 0.0);
3603 Array<OneD, NekDouble> coeffs1(nm, 0.0);
3604 Array<OneD, NekDouble> coeffs2(nm, 0.0);
3605
3606 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3607
3608 Exp->GetCoords(xc, yc, zc);
3609
3610 for (int i = 0; i < nq; ++i)
3611 {
3612 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3613 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3614 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3615 }
3616
3617 // Standard routines
3618 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
3619 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
3620 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3621 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
3622 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3623
3624 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3625 coeffs2);
3626
3627 double epsilon = 1.0e-8;
3628 for (int i = 0; i < coeffs1.size(); ++i)
3629 {
3630 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3631 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3632 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3633 }
3634}
3635
3636BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_StdMat_VariableP_MultiElmt)
3637{
3639 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3641 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3643 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3645 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3647 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3649 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3651 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3653 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3654
3655 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3656 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3658 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3659 v6.get(), v7.get(), segVec, faceVec);
3660
3661 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3663 Nektar::LibUtilities::BasisType basisTypeDir1 =
3665 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
3666 quadPointsTypeDir1);
3667 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
3668 quadPointsTypeDir1);
3669 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
3670 quadPointsTypeDir1);
3671 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3672 quadPointsKeyDir1);
3673 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
3674 quadPointsKeyDir2);
3675 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
3676 quadPointsKeyDir3);
3677
3680 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
3681
3682 int nelmts = 10;
3683
3684 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3685 for (int i = 0; i < nelmts; ++i)
3686 {
3687 CollExp.push_back(Exp);
3688 }
3689
3691 Collections::CollectionOptimisation colOpt(dummySession, 3,
3693 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3694 Collections::Collection c(CollExp, impTypes);
3696
3697 const int nq = Exp->GetTotPoints();
3698 const int nm = Exp->GetNcoeffs();
3699 Array<OneD, NekDouble> phys1(nelmts * nq, 0.0);
3700 Array<OneD, NekDouble> phys2(nelmts * nq, 0.0);
3701 Array<OneD, NekDouble> phys3(nelmts * nq, 0.0);
3702 Array<OneD, NekDouble> coeffs1(nelmts * nm, 0.0);
3703 Array<OneD, NekDouble> coeffs2(nelmts * nm, 0.0);
3704 Array<OneD, NekDouble> tmp;
3705
3706 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3707
3708 Exp->GetCoords(xc, yc, zc);
3709
3710 for (int i = 0; i < nq; ++i)
3711 {
3712 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3713 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3714 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3715 }
3716
3717 for (int i = 1; i < nelmts; ++i)
3718 {
3719 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
3720 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
3721 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
3722 }
3723
3724 // Standard routines
3725 for (int i = 0; i < nelmts; ++i)
3726 {
3727
3728 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
3729 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
3730 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3731 tmp = coeffs1 + i * nm, 1);
3732 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
3733 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3734 tmp = coeffs1 + i * nm, 1);
3735 }
3736
3737 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3738 coeffs2);
3739
3740 double epsilon = 1.0e-8;
3741 for (int i = 0; i < coeffs1.size(); ++i)
3742 {
3743 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3744 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3745 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3746 }
3747}
3748
3750 TestHexIProductWRTDerivBase_NoCollection_VariableP_MultiElmt)
3751{
3753 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3755 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3757 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3759 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3761 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3763 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3765 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3767 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3768
3769 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3770 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3772 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3773 v6.get(), v7.get(), segVec, faceVec);
3774
3775 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3777 Nektar::LibUtilities::BasisType basisTypeDir1 =
3779 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
3780 quadPointsTypeDir1);
3781 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
3782 quadPointsTypeDir1);
3783 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
3784 quadPointsTypeDir1);
3785 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3786 quadPointsKeyDir1);
3787 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
3788 quadPointsKeyDir2);
3789 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
3790 quadPointsKeyDir3);
3791
3794 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
3795
3796 int nelmts = 10;
3797
3798 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3799 for (int i = 0; i < nelmts; ++i)
3800 {
3801 CollExp.push_back(Exp);
3802 }
3803
3805 Collections::CollectionOptimisation colOpt(dummySession, 3,
3807 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3808 Collections::Collection c(CollExp, impTypes);
3810
3811 const int nq = Exp->GetTotPoints();
3812 const int nm = Exp->GetNcoeffs();
3813 Array<OneD, NekDouble> phys1(nelmts * nq, 0.0);
3814 Array<OneD, NekDouble> phys2(nelmts * nq, 0.0);
3815 Array<OneD, NekDouble> phys3(nelmts * nq, 0.0);
3816 Array<OneD, NekDouble> coeffs1(nelmts * nm, 0.0);
3817 Array<OneD, NekDouble> coeffs2(nelmts * nm, 0.0);
3818 Array<OneD, NekDouble> tmp;
3819 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3820
3821 Exp->GetCoords(xc, yc, zc);
3822
3823 for (int i = 0; i < nq; ++i)
3824 {
3825 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3826 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3827 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3828 }
3829
3830 for (int i = 1; i < nelmts; ++i)
3831 {
3832 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
3833 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
3834 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
3835 }
3836
3837 // Standard routines
3838 for (int i = 0; i < nelmts; ++i)
3839 {
3840
3841 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
3842 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
3843 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3844 tmp = coeffs1 + i * nm, 1);
3845 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
3846 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3847 tmp = coeffs1 + i * nm, 1);
3848 }
3849
3850 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3851 coeffs2);
3852
3853 double epsilon = 1.0e-6;
3854 for (int i = 0; i < coeffs1.size(); ++i)
3855 {
3856 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3857 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3858 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3859 }
3860}
3861
3862BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_SumFac_UniformP)
3863{
3865 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3867 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3869 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3871 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3873 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3875 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3877 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3879 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3880
3881 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3882 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3884 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3885 v6.get(), v7.get(), segVec, faceVec);
3886
3887 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3889 Nektar::LibUtilities::BasisType basisTypeDir1 =
3891 unsigned int numQuadPoints = 6;
3892 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
3893 quadPointsTypeDir1);
3894 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3895 quadPointsKeyDir1);
3896
3899 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
3900
3901 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3902 CollExp.push_back(Exp);
3903
3905 Collections::CollectionOptimisation colOpt(dummySession, 3,
3907 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3908 Collections::Collection c(CollExp, impTypes);
3910
3911 const int nq = Exp->GetTotPoints();
3912 const int nm = Exp->GetNcoeffs();
3913 Array<OneD, NekDouble> phys1(nq, 0.0);
3914 Array<OneD, NekDouble> phys2(nq, 0.0);
3915 Array<OneD, NekDouble> phys3(nq, 0.0);
3916 Array<OneD, NekDouble> coeffs1(nm, 0.0);
3917 Array<OneD, NekDouble> coeffs2(nm, 0.0);
3918
3919 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3920
3921 Exp->GetCoords(xc, yc, zc);
3922
3923 for (int i = 0; i < nq; ++i)
3924 {
3925 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3926 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3927 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3928 }
3929
3930 // Standard routines
3931 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
3932 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
3933 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3934 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
3935 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
3936
3937 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3938 coeffs2);
3939
3940 double epsilon = 1.0e-8;
3941 for (int i = 0; i < coeffs1.size(); ++i)
3942 {
3943 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3944 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3945 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3946 }
3947}
3948
3949BOOST_AUTO_TEST_CASE(TestHexIProductWRTDerivBase_SumFac_VariableP_MultiElmt)
3950{
3952 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3954 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3956 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3958 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3960 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3962 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
3964 new SpatialDomains::PointGeom(3u, 6u, 1.0, 1.0, 1.0));
3966 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
3967
3968 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
3969 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
3971 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
3972 v6.get(), v7.get(), segVec, faceVec);
3973
3974 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
3976 Nektar::LibUtilities::BasisType basisTypeDir1 =
3978 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(5,
3979 quadPointsTypeDir1);
3980 const Nektar::LibUtilities::PointsKey quadPointsKeyDir2(6,
3981 quadPointsTypeDir1);
3982 const Nektar::LibUtilities::PointsKey quadPointsKeyDir3(8,
3983 quadPointsTypeDir1);
3984 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3985 quadPointsKeyDir1);
3986 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir1, 6,
3987 quadPointsKeyDir2);
3988 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir1, 8,
3989 quadPointsKeyDir3);
3990
3993 basisKeyDir1, basisKeyDir2, basisKeyDir3, hexGeom.get());
3994
3995 int nelmts = 10;
3996
3997 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3998 for (int i = 0; i < nelmts; ++i)
3999 {
4000 CollExp.push_back(Exp);
4001 }
4002
4004 Collections::CollectionOptimisation colOpt(dummySession, 3,
4006 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4007 Collections::Collection c(CollExp, impTypes);
4009
4010 const int nq = Exp->GetTotPoints();
4011 const int nm = Exp->GetNcoeffs();
4012 Array<OneD, NekDouble> phys1(nelmts * nq, 0.0);
4013 Array<OneD, NekDouble> phys2(nelmts * nq, 0.0);
4014 Array<OneD, NekDouble> phys3(nelmts * nq, 0.0);
4015 Array<OneD, NekDouble> coeffs1(nelmts * nm, 0.0);
4016 Array<OneD, NekDouble> coeffs2(nelmts * nm, 0.0);
4017 Array<OneD, NekDouble> tmp;
4018
4019 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
4020
4021 Exp->GetCoords(xc, yc, zc);
4022
4023 for (int i = 0; i < nq; ++i)
4024 {
4025 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
4026 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
4027 phys2[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
4028 }
4029
4030 for (int i = 1; i < nelmts; ++i)
4031 {
4032 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
4033 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
4034 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
4035 }
4036
4037 // Standard routines
4038 for (int i = 0; i < nelmts; ++i)
4039 {
4040 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
4041 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
4042 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
4043 tmp = coeffs1 + i * nm, 1);
4044 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
4045 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
4046 tmp = coeffs1 + i * nm, 1);
4047 }
4048
4049 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
4050 coeffs2);
4051
4052 double epsilon = 1.0e-8;
4053 for (int i = 0; i < coeffs1.size(); ++i)
4054 {
4055 coeffs1[i] = (std::abs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
4056 coeffs2[i] = (std::abs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
4057 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
4058 }
4059}
4060
4061BOOST_AUTO_TEST_CASE(TestHexHelmholtz_NoCollection_UniformP)
4062{
4064 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4066 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4068 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4070 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4072 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4074 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4076 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4078 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4079
4080 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4081 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4083 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4084 v6.get(), v7.get(), segVec, faceVec);
4085
4086 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4088 Nektar::LibUtilities::BasisType basisTypeDir1 =
4090 unsigned int numQuadPoints = 5;
4091 unsigned int numModes = 4;
4092
4093 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4094 quadPointsTypeDir1);
4095 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4096 quadPointsKeyDir1);
4097
4100 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4101
4102 int nelmts = 10;
4103
4104 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4105 for (int i = 0; i < nelmts; ++i)
4106 {
4107 CollExp.push_back(Exp);
4108 }
4109
4111 Collections::CollectionOptimisation colOpt(dummySession, 2,
4113 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4114 Collections::Collection c(CollExp, impTypes);
4117
4118 c.Initialise(Collections::eHelmholtz, factors);
4119
4120 const int nm = Exp->GetNcoeffs();
4121 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4122 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4123 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4124
4125 for (int i = 0; i < coeffsIn.size(); ++i)
4126 {
4127 coeffsIn[i] = i + 1.0;
4128 }
4129
4130 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4131 *Exp, factors);
4132
4133 for (int i = 0; i < nelmts; ++i)
4134 {
4135 // Standard routines
4136 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4137 }
4138
4139 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4140
4141 double epsilon = 1.0e-8;
4142 for (int i = 0; i < coeffsRef.size(); ++i)
4143 {
4144 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4145 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4146 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4147 }
4148}
4149
4150BOOST_AUTO_TEST_CASE(TestHexHelmholtz_IterPerExp_UniformP)
4151{
4153 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4155 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4157 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4159 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4161 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4163 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4165 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4167 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4168
4169 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4170 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4172 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4173 v6.get(), v7.get(), segVec, faceVec);
4174
4175 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4177 Nektar::LibUtilities::BasisType basisTypeDir1 =
4179 unsigned int numQuadPoints = 5;
4180 unsigned int numModes = 4;
4181
4182 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4183 quadPointsTypeDir1);
4184 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4185 quadPointsKeyDir1);
4186
4189 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4190
4191 int nelmts = 10;
4192
4193 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4194 for (int i = 0; i < nelmts; ++i)
4195 {
4196 CollExp.push_back(Exp);
4197 }
4198
4200 Collections::CollectionOptimisation colOpt(dummySession, 2,
4202 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4203 Collections::Collection c(CollExp, impTypes);
4206
4207 c.Initialise(Collections::eHelmholtz, factors);
4208
4209 const int nm = Exp->GetNcoeffs();
4210 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4211 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4212 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4213
4214 for (int i = 0; i < coeffsIn.size(); ++i)
4215 {
4216 coeffsIn[i] = i + 1.0;
4217 }
4218
4219 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4220 *Exp, factors);
4221
4222 for (int i = 0; i < nelmts; ++i)
4223 {
4224 // Standard routines
4225 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4226 }
4227
4228 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4229
4230 double epsilon = 1.0e-8;
4231 for (int i = 0; i < coeffsRef.size(); ++i)
4232 {
4233 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4234 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4235 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4236 }
4237}
4238
4239BOOST_AUTO_TEST_CASE(TestHexHelmholtz_IterPerExp_UniformP_ConstVarDiff)
4240{
4242 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4244 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4246 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4248 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4250 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4252 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4254 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4256 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4257
4258 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4259 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4261 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4262 v6.get(), v7.get(), segVec, faceVec);
4263
4264 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4266 Nektar::LibUtilities::BasisType basisTypeDir1 =
4268 unsigned int numQuadPoints = 5;
4269 unsigned int numModes = 4;
4270
4271 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4272 quadPointsTypeDir1);
4273 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4274 quadPointsKeyDir1);
4275
4278 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4279
4280 int nelmts = 10;
4281
4282 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4283 for (int i = 0; i < nelmts; ++i)
4284 {
4285 CollExp.push_back(Exp);
4286 }
4287
4289 Collections::CollectionOptimisation colOpt(dummySession, 2,
4291 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4292 Collections::Collection c(CollExp, impTypes);
4301
4302 c.Initialise(Collections::eHelmholtz, factors);
4303
4304 const int nm = Exp->GetNcoeffs();
4305 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4306 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4307 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4308
4309 for (int i = 0; i < coeffsIn.size(); ++i)
4310 {
4311 coeffsIn[i] = i + 1.0;
4312 }
4313
4314 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4315 *Exp, factors);
4316
4317 for (int i = 0; i < nelmts; ++i)
4318 {
4319 // Standard routines
4320 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4321 }
4322
4323 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4324
4325 double epsilon = 1.0e-8;
4326 for (int i = 0; i < coeffsRef.size(); ++i)
4327 {
4328 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4329 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4330 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4331 }
4332}
4333
4334BOOST_AUTO_TEST_CASE(TestHexHelmholtz_MatrixFree_UniformP)
4335{
4337 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4339 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4341 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4343 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4345 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4347 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4349 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4351 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4352
4353 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4354 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4356 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4357 v6.get(), v7.get(), segVec, faceVec);
4358
4359 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4361 Nektar::LibUtilities::BasisType basisTypeDir1 =
4363 unsigned int numQuadPoints = 5;
4364 unsigned int numModes = 4;
4365
4366 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4367 quadPointsTypeDir1);
4368 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4369 quadPointsKeyDir1);
4370
4373 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4374
4375 int nelmts = 10;
4376
4377 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4378 for (int i = 0; i < nelmts; ++i)
4379 {
4380 CollExp.push_back(Exp);
4381 }
4382
4384 Collections::CollectionOptimisation colOpt(dummySession, 2,
4386 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4387 Collections::Collection c(CollExp, impTypes);
4390
4391 c.Initialise(Collections::eHelmholtz, factors);
4392
4393 const int nm = Exp->GetNcoeffs();
4394 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4395 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4396 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4397
4398 for (int i = 0; i < coeffsIn.size(); ++i)
4399 {
4400 coeffsIn[i] = i + 1.0;
4401 }
4402
4403 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4404 *Exp, factors);
4405
4406 for (int i = 0; i < nelmts; ++i)
4407 {
4408 // Standard routines
4409 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4410 }
4411
4412 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4413
4414 double epsilon = 1.0e-8;
4415 for (int i = 0; i < coeffsRef.size(); ++i)
4416 {
4417 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4418 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4419 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4420 }
4421}
4422
4423BOOST_AUTO_TEST_CASE(TestHexHelmholtz_MatrixFree_UniformP_Deformed_OverInt)
4424{
4426 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4428 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4430 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4432 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4434 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4436 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4438 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4440 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4441
4442 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4443 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4445 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4446 v6.get(), v7.get(), segVec, faceVec);
4447
4448 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4450 Nektar::LibUtilities::BasisType basisTypeDir1 =
4452 unsigned int numQuadPoints = 8;
4453 unsigned int numModes = 4;
4454
4455 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4456 quadPointsTypeDir1);
4457 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4458 quadPointsKeyDir1);
4459
4462 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4463
4464 int nelmts = 10;
4465
4466 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4467 for (int i = 0; i < nelmts; ++i)
4468 {
4469 CollExp.push_back(Exp);
4470 }
4471
4473 Collections::CollectionOptimisation colOpt(dummySession, 2,
4475 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4476 Collections::Collection c(CollExp, impTypes);
4479
4480 c.Initialise(Collections::eHelmholtz, factors);
4481
4482 const int nm = Exp->GetNcoeffs();
4483 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4484 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4485 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4486
4487 for (int i = 0; i < coeffsIn.size(); ++i)
4488 {
4489 coeffsIn[i] = i + 1.0;
4490 }
4491
4492 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4493 *Exp, factors);
4494
4495 for (int i = 0; i < nelmts; ++i)
4496 {
4497 // Standard routines
4498 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4499 }
4500
4501 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4502
4503 double epsilon = 1.0e-8;
4504 for (int i = 0; i < coeffsRef.size(); ++i)
4505 {
4506 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4507 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4508 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4509 }
4510}
4511
4512BOOST_AUTO_TEST_CASE(TestHexHelmholtz_MatrixFree_UniformP_ConstVarDiff)
4513{
4515 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4517 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4519 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4521 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4523 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4525 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4527 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4529 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4530
4531 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4532 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4534 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(),
4535 v6.get(), v7.get(), segVec, faceVec);
4536
4537 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4539 Nektar::LibUtilities::BasisType basisTypeDir1 =
4541 unsigned int numQuadPoints = 5;
4542 unsigned int numModes = 4;
4543
4544 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4545 quadPointsTypeDir1);
4546 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4547 quadPointsKeyDir1);
4548
4551 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4552
4553 int nelmts = 10;
4554
4555 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4556 for (int i = 0; i < nelmts; ++i)
4557 {
4558 CollExp.push_back(Exp);
4559 }
4560
4562 Collections::CollectionOptimisation colOpt(dummySession, 2,
4564 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4565 Collections::Collection c(CollExp, impTypes);
4574
4575 c.Initialise(Collections::eHelmholtz, factors);
4576
4577 const int nm = Exp->GetNcoeffs();
4578 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4579 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4580 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4581
4582 for (int i = 0; i < coeffsIn.size(); ++i)
4583 {
4584 coeffsIn[i] = i + 1.0;
4585 }
4586
4587 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
4588 *Exp, factors);
4589
4590 for (int i = 0; i < nelmts; ++i)
4591 {
4592 // Standard routines
4593 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4594 }
4595
4596 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
4597
4598 double epsilon = 1.0e-8;
4599 for (int i = 0; i < coeffsRef.size(); ++i)
4600 {
4601 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4602 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4603 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4604 }
4605}
4606
4607BOOST_AUTO_TEST_CASE(TestHexPhysInterp1D_NoCollection_UniformP)
4608{
4610 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4612 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4614 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4616 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4618 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4620 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4622 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4624 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4625
4626 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4627 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4629 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(), v6.get(), v7.get(), segVec, faceVec);
4630
4631 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4633 Nektar::LibUtilities::BasisType basisTypeDir1 =
4635 unsigned int numQuadPoints = 5;
4636 unsigned int numModes = 4;
4637
4638 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4639 quadPointsTypeDir1);
4640 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4641 quadPointsKeyDir1);
4642
4645 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4646
4647 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4648 CollExp.push_back(Exp);
4649
4651 Collections::CollectionOptimisation colOpt(dummySession, 3,
4653 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4654 Collections::Collection c(CollExp, impTypes);
4655
4658 c.Initialise(Collections::ePhysInterp1DScaled, factors);
4659
4660 const int nq = Exp->GetTotPoints();
4661
4662 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
4663 Array<OneD, NekDouble> phys(nq), tmp;
4664
4665 Exp->GetCoords(xc, yc, zc);
4666
4667 for (int i = 0; i < nq; ++i)
4668 {
4669 phys[i] = pow(xc[i], 3) + pow(yc[i], 3) + pow(zc[i], 3);
4670 }
4671
4672 const int nq1 = c.GetOutputSize(Collections::ePhysInterp1DScaled);
4673 Array<OneD, NekDouble> xc1(nq1);
4674 Array<OneD, NekDouble> yc1(nq1);
4675 Array<OneD, NekDouble> zc1(nq1);
4676 Array<OneD, NekDouble> phys1(nq1);
4677
4678 c.ApplyOperator(Collections::ePhysInterp1DScaled, xc, xc1);
4679 c.ApplyOperator(Collections::ePhysInterp1DScaled, yc, yc1);
4680 c.ApplyOperator(Collections::ePhysInterp1DScaled, zc, zc1);
4681 c.ApplyOperator(Collections::ePhysInterp1DScaled, phys, phys1);
4682
4683 double epsilon = 1.0e-8;
4684 // since solution is a polynomial should be able to compare soln directly
4685 for (int i = 0; i < nq1; ++i)
4686 {
4687 NekDouble exact = pow(xc1[i], 3) + pow(yc1[i], 3) + pow(zc1[i], 3);
4688 phys1[i] = (fabs(phys1[i]) < 1e-14) ? 0.0 : phys1[i];
4689 exact = (fabs(exact) < 1e-14) ? 0.0 : exact;
4690 BOOST_CHECK_CLOSE(phys1[i], exact, epsilon);
4691 }
4692}
4693
4694BOOST_AUTO_TEST_CASE(TestHexPhysInterp1D_MatrixFree_UniformP)
4695{
4697 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4699 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4701 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4703 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4705 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4707 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4709 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4711 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4712
4713 std::vector<SpatialDomains::SegGeomUniquePtr> segVec;
4714 std::vector<SpatialDomains::QuadGeomUniquePtr> faceVec;
4716 CreateHex(v0.get(), v1.get(), v2.get(), v3.get(), v4.get(), v5.get(), v6.get(), v7.get(), segVec, faceVec);
4717
4718 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4720 Nektar::LibUtilities::BasisType basisTypeDir1 =
4722 unsigned int numQuadPoints = 5;
4723 unsigned int numModes = 4;
4724
4725 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4726 quadPointsTypeDir1);
4727 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4728 quadPointsKeyDir1);
4729
4732 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4733
4734 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4735 CollExp.push_back(Exp);
4736
4738 Collections::CollectionOptimisation colOpt(dummySession, 3,
4740 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4741 Collections::Collection c(CollExp, impTypes);
4742
4745 c.Initialise(Collections::ePhysInterp1DScaled, factors);
4746
4747 const int nq = Exp->GetTotPoints();
4748
4749 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
4750 Array<OneD, NekDouble> phys(nq), tmp;
4751
4752 Exp->GetCoords(xc, yc, zc);
4753
4754 for (int i = 0; i < nq; ++i)
4755 {
4756 phys[i] = pow(xc[i], 3) + pow(yc[i], 3) + pow(zc[i], 3);
4757 }
4758
4759 const int nq1 = c.GetOutputSize(Collections::ePhysInterp1DScaled);
4760 Array<OneD, NekDouble> xc1(nq1);
4761 Array<OneD, NekDouble> yc1(nq1);
4762 Array<OneD, NekDouble> zc1(nq1);
4763 Array<OneD, NekDouble> phys1(nq1);
4764
4765 c.ApplyOperator(Collections::ePhysInterp1DScaled, xc, xc1);
4766 c.ApplyOperator(Collections::ePhysInterp1DScaled, yc, yc1);
4767 c.ApplyOperator(Collections::ePhysInterp1DScaled, zc, zc1);
4768 c.ApplyOperator(Collections::ePhysInterp1DScaled, phys, phys1);
4769
4770 double epsilon = 1.0e-8;
4771 // since solution is a polynomial should be able to compare soln directly
4772 for (int i = 0; i < nq1; ++i)
4773 {
4774 NekDouble exact = pow(xc1[i], 3) + pow(yc1[i], 3) + pow(zc1[i], 3);
4775 phys1[i] = (fabs(phys1[i]) < 1e-14) ? 0.0 : phys1[i];
4776 exact = (fabs(exact) < 1e-14) ? 0.0 : exact;
4777 BOOST_CHECK_CLOSE(phys1[i], exact, epsilon);
4778 }
4779}
4780
4782 TestHexLinearAdvectionDiffusionReaction_NoCollection_UniformP)
4783{
4785 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4787 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4789 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4791 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4793 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4795 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4797 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4799 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4800
4801 std::array<SpatialDomains::PointGeom *, 8> v = {
4802 v0.get(), v1.get(), v2.get(), v3.get(),
4803 v4.get(), v5.get(), v6.get(), v7.get()};
4804 std::array<SpatialDomains::SegGeomUniquePtr, 12> segVec;
4805 std::array<SpatialDomains::QuadGeomUniquePtr, 6> faceVec;
4806 SpatialDomains::HexGeomUniquePtr hexGeom = CreateHex(v, segVec, faceVec);
4807
4808 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4810 Nektar::LibUtilities::BasisType basisTypeDir1 =
4812 unsigned int numQuadPoints = 5;
4813 unsigned int numModes = 4;
4814
4815 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4816 quadPointsTypeDir1);
4817 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4818 quadPointsKeyDir1);
4819
4822 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4823
4824 int nelmts = 10;
4825
4826 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4827 for (int i = 0; i < nelmts; ++i)
4828 {
4829 CollExp.push_back(Exp);
4830 }
4831
4833 Collections::CollectionOptimisation colOpt(dummySession, 2,
4835 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4836 Collections::Collection c(CollExp, impTypes);
4839
4841
4842 // Add advection velocities via varcoeffs
4843 int npoints = Exp->GetTotPoints() * nelmts;
4844 StdRegions::VarCoeffMap varcoeffs;
4848 for (int i = 0; i < Exp->GetShapeDimension(); i++)
4849 {
4850 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
4851 }
4853 varcoeffs);
4854
4855 const int nm = Exp->GetNcoeffs();
4856 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4857 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4858 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4859
4860 for (int i = 0; i < coeffsIn.size(); ++i)
4861 {
4862 coeffsIn[i] = i + 1.0;
4863 }
4864
4865 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
4866 Exp->DetShapeType(), *Exp, factors,
4867 varcoeffs);
4868
4869 for (int i = 0; i < nelmts; ++i)
4870 {
4871 // Standard routines
4872 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4873 }
4874
4875 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
4876 coeffs);
4877
4878 double epsilon = 1.0e-8;
4879 for (int i = 0; i < coeffsRef.size(); ++i)
4880 {
4881 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4882 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4883 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4884 }
4885}
4886
4888 TestHexLinearAdvectionDiffusionReaction_IterPerExp_UniformP)
4889{
4891 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4893 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4895 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4897 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4899 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4901 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
4903 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
4905 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
4906
4907 std::array<SpatialDomains::PointGeom *, 8> v = {
4908 v0.get(), v1.get(), v2.get(), v3.get(),
4909 v4.get(), v5.get(), v6.get(), v7.get()};
4910 std::array<SpatialDomains::SegGeomUniquePtr, 12> segVec;
4911 std::array<SpatialDomains::QuadGeomUniquePtr, 6> faceVec;
4912 SpatialDomains::HexGeomUniquePtr hexGeom = CreateHex(v, segVec, faceVec);
4913
4914 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
4916 Nektar::LibUtilities::BasisType basisTypeDir1 =
4918 unsigned int numQuadPoints = 5;
4919 unsigned int numModes = 4;
4920
4921 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
4922 quadPointsTypeDir1);
4923 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4924 quadPointsKeyDir1);
4925
4928 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
4929
4930 int nelmts = 10;
4931
4932 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4933 for (int i = 0; i < nelmts; ++i)
4934 {
4935 CollExp.push_back(Exp);
4936 }
4937
4939 Collections::CollectionOptimisation colOpt(dummySession, 2,
4941 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4942 Collections::Collection c(CollExp, impTypes);
4945
4947
4948 // Add advection velocities via varcoeffs
4949 int npoints = Exp->GetTotPoints() * nelmts;
4950 StdRegions::VarCoeffMap varcoeffs;
4954 for (int i = 0; i < Exp->GetShapeDimension(); i++)
4955 {
4956 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
4957 }
4959 varcoeffs);
4960
4961 const int nm = Exp->GetNcoeffs();
4962 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4963 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4964 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4965
4966 for (int i = 0; i < coeffsIn.size(); ++i)
4967 {
4968 coeffsIn[i] = i + 1.0;
4969 }
4970
4971 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
4972 Exp->DetShapeType(), *Exp, factors,
4973 varcoeffs);
4974
4975 for (int i = 0; i < nelmts; ++i)
4976 {
4977 // Standard routines
4978 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4979 }
4980
4981 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
4982 coeffs);
4983
4984 double epsilon = 1.0e-8;
4985 for (int i = 0; i < coeffsRef.size(); ++i)
4986 {
4987 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4988 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4989 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4990 }
4991}
4992
4994 TestHexLinearAdvectionDiffusionReaction_MatrixFree_UniformP)
4995{
4997 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
4999 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
5001 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
5003 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
5005 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
5007 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
5009 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
5011 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
5012
5013 std::array<SpatialDomains::PointGeom *, 8> v = {
5014 v0.get(), v1.get(), v2.get(), v3.get(),
5015 v4.get(), v5.get(), v6.get(), v7.get()};
5016 std::array<SpatialDomains::SegGeomUniquePtr, 12> segVec;
5017 std::array<SpatialDomains::QuadGeomUniquePtr, 6> faceVec;
5018 SpatialDomains::HexGeomUniquePtr hexGeom = CreateHex(v, segVec, faceVec);
5019
5020 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
5022 Nektar::LibUtilities::BasisType basisTypeDir1 =
5024 unsigned int numQuadPoints = 5;
5025 unsigned int numModes = 4;
5026
5027 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
5028 quadPointsTypeDir1);
5029 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
5030 quadPointsKeyDir1);
5031
5034 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
5035
5036 int nelmts = 10;
5037
5038 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
5039 for (int i = 0; i < nelmts; ++i)
5040 {
5041 CollExp.push_back(Exp);
5042 }
5043
5045 Collections::CollectionOptimisation colOpt(dummySession, 2,
5047 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
5048 Collections::Collection c(CollExp, impTypes);
5051
5053
5054 // Add advection velocities via varcoeffs
5055 int npoints = Exp->GetTotPoints() * nelmts;
5056 StdRegions::VarCoeffMap varcoeffs;
5060 for (int i = 0; i < Exp->GetShapeDimension(); i++)
5061 {
5062 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
5063 }
5065 varcoeffs);
5066
5067 const int nm = Exp->GetNcoeffs();
5068 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
5069 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
5070 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
5071
5072 for (int i = 0; i < coeffsIn.size(); ++i)
5073 {
5074 coeffsIn[i] = i + 1.0;
5075 }
5076
5077 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
5078 Exp->DetShapeType(), *Exp, factors,
5079 varcoeffs);
5080
5081 for (int i = 0; i < nelmts; ++i)
5082 {
5083 // Standard routines
5084 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
5085 }
5086
5087 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
5088 coeffs);
5089
5090 double epsilon = 1.0e-8;
5091 for (int i = 0; i < coeffsRef.size(); ++i)
5092 {
5093 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
5094 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
5095 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
5096 }
5097}
5098
5100 TestHexLinearAdvectionDiffusionReaction_MatrixFree_UniformP_Deformed_OverInt)
5101{
5103 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
5105 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
5107 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
5109 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
5111 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
5113 new SpatialDomains::PointGeom(3u, 5u, 1.0, -1.0, 1.0));
5115 new SpatialDomains::PointGeom(3u, 6u, 2.0, 3.0, 4.0));
5117 new SpatialDomains::PointGeom(3u, 7u, -1.0, 1.0, 1.0));
5118
5119 std::array<SpatialDomains::PointGeom *, 8> v = {
5120 v0.get(), v1.get(), v2.get(), v3.get(),
5121 v4.get(), v5.get(), v6.get(), v7.get()};
5122 std::array<SpatialDomains::SegGeomUniquePtr, 12> segVec;
5123 std::array<SpatialDomains::QuadGeomUniquePtr, 6> faceVec;
5124 SpatialDomains::HexGeomUniquePtr hexGeom = CreateHex(v, segVec, faceVec);
5125
5126 Nektar::LibUtilities::PointsType quadPointsTypeDir1 =
5128 Nektar::LibUtilities::BasisType basisTypeDir1 =
5130 unsigned int numQuadPoints = 8;
5131 unsigned int numModes = 4;
5132
5133 const Nektar::LibUtilities::PointsKey quadPointsKeyDir1(numQuadPoints,
5134 quadPointsTypeDir1);
5135 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
5136 quadPointsKeyDir1);
5137
5140 basisKeyDir1, basisKeyDir1, basisKeyDir1, hexGeom.get());
5141
5142 int nelmts = 10;
5143
5144 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
5145 for (int i = 0; i < nelmts; ++i)
5146 {
5147 CollExp.push_back(Exp);
5148 }
5149
5151 Collections::CollectionOptimisation colOpt(dummySession, 2,
5153 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
5154 Collections::Collection c(CollExp, impTypes);
5157
5159
5160 // Add advection velocities via varcoeffs
5161 int npoints = Exp->GetTotPoints() * nelmts;
5162 StdRegions::VarCoeffMap varcoeffs;
5166 for (int i = 0; i < Exp->GetShapeDimension(); i++)
5167 {
5168 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
5169 }
5171 varcoeffs);
5172
5173 const int nm = Exp->GetNcoeffs();
5174 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
5175 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
5176 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
5177
5178 for (int i = 0; i < coeffsIn.size(); ++i)
5179 {
5180 coeffsIn[i] = i + 1.0;
5181 }
5182
5183 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
5184 Exp->DetShapeType(), *Exp, factors,
5185 varcoeffs);
5186
5187 for (int i = 0; i < nelmts; ++i)
5188 {
5189 // Standard routines
5190 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
5191 }
5192
5193 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
5194 coeffs);
5195
5196 double epsilon = 1.0e-8;
5197 for (int i = 0; i < coeffsRef.size(); ++i)
5198 {
5199 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
5200 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
5201 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
5202 }
5203}
5204
5205#endif
5206} // namespace Nektar::HexCollectionTests
COLLECTIONS_EXPORT void Initialise(const OperatorType opType, StdRegions::FactorMap factors=StdRegions::NullFactorMap)
void ApplyOperator(const OperatorType &op, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &output)
Definition Collection.h:148
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(LocalRegions::ExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Describes the specification for a Basis.
Definition Basis.h:45
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.
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition Operator.h:131
@ eLinearAdvectionDiffusionReaction
Definition Operator.h:66
The above copyright notice and this permission notice shall be included.
BOOST_AUTO_TEST_CASE(TestHexBwdTrans_StdMat_UniformP)
SpatialDomains::SegGeomUniquePtr CreateSegGeom(unsigned int id, SpatialDomains::PointGeom *v0, SpatialDomains::PointGeom *v1)
SpatialDomains::HexGeomUniquePtr CreateHex(std::array< SpatialDomains::PointGeom *, 8 > v, std::array< SpatialDomains::SegGeomUniquePtr, 12 > &segVec, std::array< SpatialDomains::QuadGeomUniquePtr, 6 > &faceVec)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition PointsType.h:51
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition BasisType.h:56
@ eModified_A
Principle Modified Functions .
Definition BasisType.h:48
std::shared_ptr< HexExp > HexExpSharedPtr
Definition HexExp.h:189
unique_ptr_objpool< HexGeom > HexGeomUniquePtr
Definition MeshGraph.h:108
unique_ptr_objpool< QuadGeom > QuadGeomUniquePtr
Definition MeshGraph.h:104
unique_ptr_objpool< SegGeom > SegGeomUniquePtr
Definition MeshGraph.h:102
unique_ptr_objpool< PointGeom > PointGeomUniquePtr
Definition MeshGraph.h:99
std::map< ConstFactorType, NekDouble > ConstFactorMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
StdRegions::ConstFactorMap factors
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition Vmath.hpp:180
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825