Nektar++
Loading...
Searching...
No Matches
TestPyrCollection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestPyrCollection.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
35#include <Collections/Collection.h>
36#include <Collections/CollectionOptimisation.h>
37#include <LibUtilities/BasicUtils/Timer.h>
38#include <LibUtilities/Communication/CommSerial.h>
39#include <LocalRegions/PyrExp.h>
40#include <SpatialDomains/MeshGraph.h>
41#include <boost/test/tools/floating_point_comparison.hpp>
42#include <boost/test/unit_test.hpp>
43
44namespace Nektar::PyrCollectionTests
45{
46#define NELMTS 10
47
48SpatialDomains::SegGeomUniquePtr CreateSegGeom(unsigned int id,
49 SpatialDomains::PointGeom *v0,
50 SpatialDomains::PointGeom *v1)
51{
52 std::array<SpatialDomains::PointGeom *, 2> vertices = {v0, v1};
53 SpatialDomains::SegGeomUniquePtr result(
54 new SpatialDomains::SegGeom(id, 3, vertices));
55 return result;
56}
57
58SpatialDomains::PyrGeomUniquePtr CreatePyr(
59 std::array<SpatialDomains::PointGeom *, 5> v,
60 std::array<SpatialDomains::SegGeomUniquePtr, 8> &segVec,
61 std::array<SpatialDomains::TriGeomUniquePtr, 4> &triVec,
62 std::array<SpatialDomains::QuadGeomUniquePtr, 1> &quadVec)
63{
64 std::array<std::array<int, 2>, 8> edgeVerts = {{{{0, 1}},
65 {{1, 2}},
66 {{3, 2}},
67 {{0, 3}},
68 {{0, 4}},
69 {{1, 4}},
70 {{2, 4}},
71 {{3, 4}}}};
72 std::array<std::array<int, 4>, 5> faceEdges = {{{{0, 1, 2, 3}},
73 {{0, 5, 4, -1}},
74 {{1, 6, 5, -1}},
75 {{2, 7, 6, -1}},
76 {{3, 7, 4, -1}}}};
77
78 // Create segments from vertices
79 for (int i = 0; i < 8; ++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::Geometry2D *, 5> faces;
86 for (int i = 0; i < 5; ++i)
87 {
88 if (i == 0)
89 {
90 // Quad faces
91 std::array<SpatialDomains::SegGeom *, 4> face;
92 for (int j = 0; j < 4; ++j)
93 {
94 face[j] = segVec[faceEdges[i][j]].get();
95 }
96 quadVec[0] = SpatialDomains::QuadGeomUniquePtr(
97 new SpatialDomains::QuadGeom(i, face));
98 faces[i] = quadVec[0].get();
99 }
100 else
101 {
102 // Tri faces
103 std::array<SpatialDomains::SegGeom *, 3> face;
104 for (int j = 0; j < 3; ++j)
105 {
106 face[j] = segVec[faceEdges[i][j]].get();
107 }
108 triVec[i - 1] = SpatialDomains::TriGeomUniquePtr(
109 new SpatialDomains::TriGeom(i, face));
110 faces[i] = triVec[i - 1].get();
111 }
112 }
113
114 SpatialDomains::PyrGeomUniquePtr pyrGeom(
115 new SpatialDomains::PyrGeom(0, faces));
116 return pyrGeom;
117}
118
119BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_IterPerExp_UniformP_MultiElmt)
120{
121 SpatialDomains::PointGeomUniquePtr v0(
122 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
123 SpatialDomains::PointGeomUniquePtr v1(
124 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
125 SpatialDomains::PointGeomUniquePtr v2(
126 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
127 SpatialDomains::PointGeomUniquePtr v3(
128 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
129 SpatialDomains::PointGeomUniquePtr v4(
130 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
131
132 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
133 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
134 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
135 std::array<SpatialDomains::PointGeom *, 5> v = {
136 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
137 SpatialDomains::PyrGeomUniquePtr pyrGeom =
138 CreatePyr(v, segVec, triVec, quadVec);
139
140 Nektar::LibUtilities::PointsType PointsTypeDir1 =
141 Nektar::LibUtilities::eGaussLobattoLegendre;
142 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
143 Nektar::LibUtilities::BasisType basisTypeDir1 =
144 Nektar::LibUtilities::eModified_A;
145 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
146 PointsKeyDir1);
147
148 Nektar::LibUtilities::PointsType PointsTypeDir2 =
149 Nektar::LibUtilities::eGaussLobattoLegendre;
150 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
151 Nektar::LibUtilities::BasisType basisTypeDir2 =
152 Nektar::LibUtilities::eModified_A;
153 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
154 PointsKeyDir2);
155
156 Nektar::LibUtilities::PointsType PointsTypeDir3 =
157 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
158 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
159 Nektar::LibUtilities::BasisType basisTypeDir3 =
160 Nektar::LibUtilities::eModifiedPyr_C;
161 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
162 PointsKeyDir3);
163
164 Nektar::LocalRegions::PyrExpSharedPtr Exp =
165 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
166 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
167
168 int nelmts = NELMTS;
169
170 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
171 for (int i = 0; i < nelmts; ++i)
172 {
173 CollExp.push_back(Exp);
174 }
175
176 LibUtilities::SessionReaderSharedPtr dummySession;
177 Collections::CollectionOptimisation colOpt(dummySession, 3,
178 Collections::eIterPerExp);
179 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
180 Collections::Collection c(CollExp, impTypes);
181 c.Initialise(Collections::eBwdTrans);
182
183 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
184 for (int i = 0; i < coeffs.size(); ++i)
185 {
186 coeffs[i] = i + 1;
187 }
188 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
189 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
190
191 for (int i = 0; i < nelmts; ++i)
192 {
193 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
194 tmp = phys1 + i * Exp->GetTotPoints());
195 }
196
197 LibUtilities::Timer timer;
198 timer.Start();
199 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
200 timer.Stop();
201 timer.AccumulateRegion("Pyr BwdTrans IterPerExp");
202
203 double epsilon = 1.0e-8;
204 for (int i = 0; i < phys1.size(); ++i)
205 {
206 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
207 }
208}
209
210BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_StdMat_UniformP_MultiElmt)
211{
212 SpatialDomains::PointGeomUniquePtr v0(
213 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
214 SpatialDomains::PointGeomUniquePtr v1(
215 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
216 SpatialDomains::PointGeomUniquePtr v2(
217 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
218 SpatialDomains::PointGeomUniquePtr v3(
219 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
220 SpatialDomains::PointGeomUniquePtr v4(
221 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
222
223 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
224 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
225 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
226 std::array<SpatialDomains::PointGeom *, 5> v = {
227 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
228 SpatialDomains::PyrGeomUniquePtr pyrGeom =
229 CreatePyr(v, segVec, triVec, quadVec);
230
231 Nektar::LibUtilities::PointsType PointsTypeDir1 =
232 Nektar::LibUtilities::eGaussLobattoLegendre;
233 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
234 Nektar::LibUtilities::BasisType basisTypeDir1 =
235 Nektar::LibUtilities::eModified_A;
236 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
237 PointsKeyDir1);
238
239 Nektar::LibUtilities::PointsType PointsTypeDir2 =
240 Nektar::LibUtilities::eGaussLobattoLegendre;
241 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
242 Nektar::LibUtilities::BasisType basisTypeDir2 =
243 Nektar::LibUtilities::eModified_A;
244 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
245 PointsKeyDir2);
246
247 Nektar::LibUtilities::PointsType PointsTypeDir3 =
248 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
249 const Nektar::LibUtilities::PointsKey PointsKeyDir3(5, PointsTypeDir3);
250 Nektar::LibUtilities::BasisType basisTypeDir3 =
251 Nektar::LibUtilities::eModifiedPyr_C;
252 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
253 PointsKeyDir3);
254
255 Nektar::LocalRegions::PyrExpSharedPtr Exp =
256 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
257 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
258
259 int nelmts = 10;
260
261 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
262 for (int i = 0; i < nelmts; ++i)
263 {
264 CollExp.push_back(Exp);
265 }
266
267 LibUtilities::SessionReaderSharedPtr dummySession;
268 Collections::CollectionOptimisation colOpt(dummySession, 3,
269 Collections::eStdMat);
270 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
271 Collections::Collection c(CollExp, impTypes);
272 c.Initialise(Collections::eBwdTrans);
273
274 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
275 for (int i = 0; i < coeffs.size(); ++i)
276 {
277 coeffs[i] = i + 1;
278 }
279 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints(), 0.0);
280 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints(), 0.0);
281 Array<OneD, NekDouble> tmp;
282
283 for (int i = 0; i < nelmts; ++i)
284 {
285 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
286 tmp = phys1 + i * Exp->GetTotPoints());
287 }
288 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
289
290 double epsilon = 1.0e-8;
291 for (int i = 0; i < phys1.size(); ++i)
292 {
293 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
294 }
295}
296
297BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_SumFac_UniformP_MultiElmt)
298{
299 SpatialDomains::PointGeomUniquePtr v0(
300 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
301 SpatialDomains::PointGeomUniquePtr v1(
302 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
303 SpatialDomains::PointGeomUniquePtr v2(
304 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
305 SpatialDomains::PointGeomUniquePtr v3(
306 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
307 SpatialDomains::PointGeomUniquePtr v4(
308 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
309
310 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
311 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
312 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
313 std::array<SpatialDomains::PointGeom *, 5> v = {
314 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
315 SpatialDomains::PyrGeomUniquePtr pyrGeom =
316 CreatePyr(v, segVec, triVec, quadVec);
317
318 Nektar::LibUtilities::PointsType PointsTypeDir1 =
319 Nektar::LibUtilities::eGaussLobattoLegendre;
320 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
321 Nektar::LibUtilities::BasisType basisTypeDir1 =
322 Nektar::LibUtilities::eModified_A;
323 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
324 PointsKeyDir1);
325
326 Nektar::LibUtilities::PointsType PointsTypeDir2 =
327 Nektar::LibUtilities::eGaussLobattoLegendre;
328 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
329 Nektar::LibUtilities::BasisType basisTypeDir2 =
330 Nektar::LibUtilities::eModified_A;
331 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
332 PointsKeyDir2);
333
334 Nektar::LibUtilities::PointsType PointsTypeDir3 =
335 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
336 const Nektar::LibUtilities::PointsKey PointsKeyDir3(5, PointsTypeDir3);
337 Nektar::LibUtilities::BasisType basisTypeDir3 =
338 Nektar::LibUtilities::eModifiedPyr_C;
339 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
340 PointsKeyDir3);
341
342 Nektar::LocalRegions::PyrExpSharedPtr Exp =
343 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
344 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
345
346 int nelmts = NELMTS;
347
348 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
349 for (int i = 0; i < nelmts; ++i)
350 {
351 CollExp.push_back(Exp);
352 }
353
354 LibUtilities::SessionReaderSharedPtr dummySession;
355 Collections::CollectionOptimisation colOpt(dummySession, 3,
356 Collections::eSumFac);
357 Collections::OperatorImpMap impTypes;
358 impTypes[Collections::eBwdTrans] = Collections::eSumFac;
359 Collections::Collection c(CollExp, impTypes);
360 c.Initialise(Collections::eBwdTrans);
361
362 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
363 for (int i = 0; i < coeffs.size(); ++i)
364 {
365 coeffs[i] = i + 1;
366 }
367 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints(), 0.0);
368 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints(), 0.0);
369 Array<OneD, NekDouble> tmp;
370
371 for (int i = 0; i < nelmts; ++i)
372 {
373 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
374 tmp = phys1 + i * Exp->GetTotPoints());
375 }
376
377 LibUtilities::Timer timer;
378 timer.Start();
379 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
380 timer.Stop();
381 timer.AccumulateRegion("Pyr BwdTrans SumFactor");
382
383 double epsilon = 1.0e-8;
384 for (int i = 0; i < phys1.size(); ++i)
385 {
386 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
387 }
388}
389
390BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_IterPerExp_VariableP_MultiElmt)
391{
392 SpatialDomains::PointGeomUniquePtr v0(
393 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
394 SpatialDomains::PointGeomUniquePtr v1(
395 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
396 SpatialDomains::PointGeomUniquePtr v2(
397 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
398 SpatialDomains::PointGeomUniquePtr v3(
399 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
400 SpatialDomains::PointGeomUniquePtr v4(
401 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
402
403 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
404 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
405 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
406 std::array<SpatialDomains::PointGeom *, 5> v = {
407 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
408 SpatialDomains::PyrGeomUniquePtr pyrGeom =
409 CreatePyr(v, segVec, triVec, quadVec);
410
411 Nektar::LibUtilities::PointsType PointsTypeDir1 =
412 Nektar::LibUtilities::eGaussLobattoLegendre;
413 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
414 Nektar::LibUtilities::BasisType basisTypeDir1 =
415 Nektar::LibUtilities::eModified_A;
416 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
417 PointsKeyDir1);
418
419 Nektar::LibUtilities::PointsType PointsTypeDir2 =
420 Nektar::LibUtilities::eGaussLobattoLegendre;
421 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
422 Nektar::LibUtilities::BasisType basisTypeDir2 =
423 Nektar::LibUtilities::eModified_A;
424 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
425 PointsKeyDir2);
426
427 Nektar::LibUtilities::PointsType PointsTypeDir3 =
428 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
429 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
430 Nektar::LibUtilities::BasisType basisTypeDir3 =
431 Nektar::LibUtilities::eModifiedPyr_C;
432 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
433 PointsKeyDir3);
434
435 Nektar::LocalRegions::PyrExpSharedPtr Exp =
436 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
437 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
438 int nelmts = 10;
439
440 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
441 for (int i = 0; i < nelmts; ++i)
442 {
443 CollExp.push_back(Exp);
444 }
445
446 LibUtilities::SessionReaderSharedPtr dummySession;
447 Collections::CollectionOptimisation colOpt(dummySession, 3,
448 Collections::eIterPerExp);
449 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
450 Collections::Collection c(CollExp, impTypes);
451 c.Initialise(Collections::eBwdTrans);
452
453 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
454 for (int i = 0; i < coeffs.size(); ++i)
455 {
456 coeffs[i] = i + 1;
457 }
458 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints(), 0.0);
459 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints(), 0.0);
460 Array<OneD, NekDouble> tmp;
461
462 for (int i = 0; i < nelmts; ++i)
463 {
464 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
465 tmp = phys1 + i * Exp->GetTotPoints());
466 }
467 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
468
469 double epsilon = 1.0e-8;
470 for (int i = 0; i < phys1.size(); ++i)
471 {
472 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
473 }
474}
475
476BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_StdMat_VariableP_MultiElmt)
477{
478 SpatialDomains::PointGeomUniquePtr v0(
479 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
480 SpatialDomains::PointGeomUniquePtr v1(
481 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
482 SpatialDomains::PointGeomUniquePtr v2(
483 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
484 SpatialDomains::PointGeomUniquePtr v3(
485 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
486 SpatialDomains::PointGeomUniquePtr v4(
487 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
488
489 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
490 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
491 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
492 std::array<SpatialDomains::PointGeom *, 5> v = {
493 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
494 SpatialDomains::PyrGeomUniquePtr pyrGeom =
495 CreatePyr(v, segVec, triVec, quadVec);
496
497 Nektar::LibUtilities::PointsType PointsTypeDir1 =
498 Nektar::LibUtilities::eGaussLobattoLegendre;
499 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
500 Nektar::LibUtilities::BasisType basisTypeDir1 =
501 Nektar::LibUtilities::eModified_A;
502 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
503 PointsKeyDir1);
504
505 Nektar::LibUtilities::PointsType PointsTypeDir2 =
506 Nektar::LibUtilities::eGaussLobattoLegendre;
507 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
508 Nektar::LibUtilities::BasisType basisTypeDir2 =
509 Nektar::LibUtilities::eModified_A;
510 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
511 PointsKeyDir2);
512
513 Nektar::LibUtilities::PointsType PointsTypeDir3 =
514 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
515 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
516 Nektar::LibUtilities::BasisType basisTypeDir3 =
517 Nektar::LibUtilities::eModifiedPyr_C;
518 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
519 PointsKeyDir3);
520
521 Nektar::LocalRegions::PyrExpSharedPtr Exp =
522 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
523 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
524
525 int nelmts = 10;
526
527 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
528 for (int i = 0; i < nelmts; ++i)
529 {
530 CollExp.push_back(Exp);
531 }
532
533 LibUtilities::SessionReaderSharedPtr dummySession;
534 Collections::CollectionOptimisation colOpt(dummySession, 3,
535 Collections::eStdMat);
536 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
537 Collections::Collection c(CollExp, impTypes);
538 c.Initialise(Collections::eBwdTrans);
539
540 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
541 for (int i = 0; i < coeffs.size(); ++i)
542 {
543 coeffs[i] = i + 1;
544 }
545 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints(), 0.0);
546 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints(), 0.0);
547 Array<OneD, NekDouble> tmp;
548
549 for (int i = 0; i < nelmts; ++i)
550 {
551 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
552 tmp = phys1 + i * Exp->GetTotPoints());
553 }
554 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
555
556 double epsilon = 1.0e-8;
557 for (int i = 0; i < phys1.size(); ++i)
558 {
559 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
560 }
561}
562
563BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_SumFac_VariableP_MultiElmt)
564{
565 SpatialDomains::PointGeomUniquePtr v0(
566 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
567 SpatialDomains::PointGeomUniquePtr v1(
568 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
569 SpatialDomains::PointGeomUniquePtr v2(
570 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
571 SpatialDomains::PointGeomUniquePtr v3(
572 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
573 SpatialDomains::PointGeomUniquePtr v4(
574 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
575
576 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
577 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
578 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
579 std::array<SpatialDomains::PointGeom *, 5> v = {
580 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
581 SpatialDomains::PyrGeomUniquePtr pyrGeom =
582 CreatePyr(v, segVec, triVec, quadVec);
583
584 Nektar::LibUtilities::PointsType PointsTypeDir1 =
585 Nektar::LibUtilities::eGaussLobattoLegendre;
586 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
587 Nektar::LibUtilities::BasisType basisTypeDir1 =
588 Nektar::LibUtilities::eModified_A;
589 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
590 PointsKeyDir1);
591
592 Nektar::LibUtilities::PointsType PointsTypeDir2 =
593 Nektar::LibUtilities::eGaussLobattoLegendre;
594 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
595 Nektar::LibUtilities::BasisType basisTypeDir2 =
596 Nektar::LibUtilities::eModified_A;
597 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
598 PointsKeyDir2);
599
600 Nektar::LibUtilities::PointsType PointsTypeDir3 =
601 Nektar::LibUtilities::eGaussLobattoLegendre;
602 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
603 Nektar::LibUtilities::BasisType basisTypeDir3 =
604 Nektar::LibUtilities::eModifiedPyr_C;
605 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
606 PointsKeyDir3);
607
608 Nektar::LocalRegions::PyrExpSharedPtr Exp =
609 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
610 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
611
612 int nelmts = 10;
613
614 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
615 for (int i = 0; i < nelmts; ++i)
616 {
617 CollExp.push_back(Exp);
618 }
619
620 LibUtilities::SessionReaderSharedPtr dummySession;
621 Collections::CollectionOptimisation colOpt(dummySession, 3,
622 Collections::eSumFac);
623 Collections::OperatorImpMap impTypes;
624 impTypes[Collections::eBwdTrans] = Collections::eSumFac;
625 Collections::Collection c(CollExp, impTypes);
626 c.Initialise(Collections::eBwdTrans);
627
628 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
629 for (int i = 0; i < coeffs.size(); ++i)
630 {
631 coeffs[i] = i + 1;
632 }
633 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints(), 0.0);
634 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints(), 0.0);
635 Array<OneD, NekDouble> tmp;
636
637 for (int i = 0; i < nelmts; ++i)
638 {
639 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
640 tmp = phys1 + i * Exp->GetTotPoints());
641 }
642 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
643
644 double epsilon = 1.0e-8;
645 for (int i = 0; i < phys1.size(); ++i)
646 {
647 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
648 }
649}
650
651BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_IterPerExp_UniformP_MultiElmt)
652{
653 SpatialDomains::PointGeomUniquePtr v0(
654 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
655 SpatialDomains::PointGeomUniquePtr v1(
656 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
657 SpatialDomains::PointGeomUniquePtr v2(
658 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
659 SpatialDomains::PointGeomUniquePtr v3(
660 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
661 SpatialDomains::PointGeomUniquePtr v4(
662 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
663
664 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
665 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
666 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
667 std::array<SpatialDomains::PointGeom *, 5> v = {
668 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
669 SpatialDomains::PyrGeomUniquePtr pyrGeom =
670 CreatePyr(v, segVec, triVec, quadVec);
671
672 Nektar::LibUtilities::PointsType PointsTypeDir1 =
673 Nektar::LibUtilities::eGaussLobattoLegendre;
674 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
675 Nektar::LibUtilities::BasisType basisTypeDir1 =
676 Nektar::LibUtilities::eModified_A;
677 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
678 PointsKeyDir1);
679
680 Nektar::LibUtilities::PointsType PointsTypeDir2 =
681 Nektar::LibUtilities::eGaussLobattoLegendre;
682 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
683 Nektar::LibUtilities::BasisType basisTypeDir2 =
684 Nektar::LibUtilities::eModified_A;
685 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
686 PointsKeyDir2);
687
688 Nektar::LibUtilities::PointsType PointsTypeDir3 =
689 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
690 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
691 Nektar::LibUtilities::BasisType basisTypeDir3 =
692 Nektar::LibUtilities::eModifiedPyr_C;
693 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
694 PointsKeyDir3);
695
696 Nektar::LocalRegions::PyrExpSharedPtr Exp =
697 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
698 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
699 int nelmts = NELMTS;
700
701 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
702 for (int i = 0; i < nelmts; ++i)
703 {
704 CollExp.push_back(Exp);
705 }
706
707 LibUtilities::SessionReaderSharedPtr dummySession;
708 Collections::CollectionOptimisation colOpt(dummySession, 3,
709 Collections::eIterPerExp);
710 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
711 Collections::Collection c(CollExp, impTypes);
712 c.Initialise(Collections::ePhysDeriv);
713
714 const int nq = Exp->GetTotPoints();
715 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
716 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
717 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
718 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
719
720 Exp->GetCoords(xc, yc, zc);
721
722 for (int i = 0; i < nq; ++i)
723 {
724 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
725 }
726 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
727 tmp2 = diff1 + (2 * nelmts) * nq);
728
729 for (int i = 1; i < nelmts; ++i)
730 {
731 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
732 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
733 tmp1 = diff1 + (nelmts + i) * nq,
734 tmp2 = diff1 + (2 * nelmts + i) * nq);
735 }
736
737 LibUtilities::Timer timer;
738 timer.Start();
739 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
740 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
741 timer.Stop();
742 timer.AccumulateRegion("Pyr PhysDeriv IterPerExp");
743
744 double epsilon = 1.0e-8;
745 for (int i = 0; i < diff1.size(); ++i)
746 {
747 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
748 }
749}
750
751BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_StdMat_UniformP_MultiElmt)
752{
753 SpatialDomains::PointGeomUniquePtr v0(
754 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
755 SpatialDomains::PointGeomUniquePtr v1(
756 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
757 SpatialDomains::PointGeomUniquePtr v2(
758 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
759 SpatialDomains::PointGeomUniquePtr v3(
760 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
761 SpatialDomains::PointGeomUniquePtr v4(
762 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
763
764 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
765 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
766 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
767 std::array<SpatialDomains::PointGeom *, 5> v = {
768 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
769 SpatialDomains::PyrGeomUniquePtr pyrGeom =
770 CreatePyr(v, segVec, triVec, quadVec);
771
772 Nektar::LibUtilities::PointsType PointsTypeDir1 =
773 Nektar::LibUtilities::eGaussLobattoLegendre;
774 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
775 Nektar::LibUtilities::BasisType basisTypeDir1 =
776 Nektar::LibUtilities::eModified_A;
777 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
778 PointsKeyDir1);
779
780 Nektar::LibUtilities::PointsType PointsTypeDir2 =
781 Nektar::LibUtilities::eGaussLobattoLegendre;
782 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
783 Nektar::LibUtilities::BasisType basisTypeDir2 =
784 Nektar::LibUtilities::eModified_A;
785 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
786 PointsKeyDir2);
787
788 Nektar::LibUtilities::PointsType PointsTypeDir3 =
789 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
790 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
791 Nektar::LibUtilities::BasisType basisTypeDir3 =
792 Nektar::LibUtilities::eModifiedPyr_C;
793 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
794 PointsKeyDir3);
795
796 Nektar::LocalRegions::PyrExpSharedPtr Exp =
797 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
798 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
799
800 int nelmts = 10;
801
802 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
803 for (int i = 0; i < nelmts; ++i)
804 {
805 CollExp.push_back(Exp);
806 }
807
808 LibUtilities::SessionReaderSharedPtr dummySession;
809 Collections::CollectionOptimisation colOpt(dummySession, 3,
810 Collections::eStdMat);
811 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
812 Collections::Collection c(CollExp, impTypes);
813 c.Initialise(Collections::ePhysDeriv);
814
815 const int nq = Exp->GetTotPoints();
816 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
817 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
818 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
819 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
820
821 Exp->GetCoords(xc, yc, zc);
822
823 for (int i = 0; i < nq; ++i)
824 {
825 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
826 }
827 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
828 tmp2 = diff1 + (2 * nelmts) * nq);
829
830 for (int i = 1; i < nelmts; ++i)
831 {
832 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
833 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
834 tmp1 = diff1 + (nelmts + i) * nq,
835 tmp2 = diff1 + (2 * nelmts + i) * nq);
836 }
837
838 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
839 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
840
841 double epsilon = 1.0e-8;
842 for (int i = 0; i < diff1.size(); ++i)
843 {
844 // clamp values below 1e-14 to zero
845 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
846 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
847 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
848 }
849}
850
851BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_SumFac_UniformP_MultiElmt)
852{
853 SpatialDomains::PointGeomUniquePtr v0(
854 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
855 SpatialDomains::PointGeomUniquePtr v1(
856 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
857 SpatialDomains::PointGeomUniquePtr v2(
858 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
859 SpatialDomains::PointGeomUniquePtr v3(
860 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
861 SpatialDomains::PointGeomUniquePtr v4(
862 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
863
864 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
865 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
866 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
867 std::array<SpatialDomains::PointGeom *, 5> v = {
868 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
869 SpatialDomains::PyrGeomUniquePtr pyrGeom =
870 CreatePyr(v, segVec, triVec, quadVec);
871
872 Nektar::LibUtilities::PointsType PointsTypeDir1 =
873 Nektar::LibUtilities::eGaussLobattoLegendre;
874 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
875 Nektar::LibUtilities::BasisType basisTypeDir1 =
876 Nektar::LibUtilities::eModified_A;
877 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
878 PointsKeyDir1);
879
880 Nektar::LibUtilities::PointsType PointsTypeDir2 =
881 Nektar::LibUtilities::eGaussLobattoLegendre;
882 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
883 Nektar::LibUtilities::BasisType basisTypeDir2 =
884 Nektar::LibUtilities::eModified_A;
885 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
886 PointsKeyDir2);
887
888 Nektar::LibUtilities::PointsType PointsTypeDir3 =
889 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
890 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
891 Nektar::LibUtilities::BasisType basisTypeDir3 =
892 Nektar::LibUtilities::eModifiedPyr_C;
893 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
894 PointsKeyDir3);
895
896 Nektar::LocalRegions::PyrExpSharedPtr Exp =
897 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
898 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
899 int nelmts = NELMTS;
900
901 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
902 for (int i = 0; i < nelmts; ++i)
903 {
904 CollExp.push_back(Exp);
905 }
906
907 LibUtilities::SessionReaderSharedPtr dummySession;
908 Collections::CollectionOptimisation colOpt(dummySession, 3,
909 Collections::eSumFac);
910 Collections::OperatorImpMap impTypes;
911 impTypes[Collections::ePhysDeriv] = Collections::eSumFac;
912 Collections::Collection c(CollExp, impTypes);
913 c.Initialise(Collections::ePhysDeriv);
914
915 const int nq = Exp->GetTotPoints();
916 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
917 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
918 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
919 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
920
921 Exp->GetCoords(xc, yc, zc);
922
923 for (int i = 0; i < nq; ++i)
924 {
925 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
926 }
927 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
928 tmp2 = diff1 + (2 * nelmts) * nq);
929
930 for (int i = 1; i < nelmts; ++i)
931 {
932 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
933 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
934 tmp1 = diff1 + (nelmts + i) * nq,
935 tmp2 = diff1 + (2 * nelmts + i) * nq);
936 }
937
938 LibUtilities::Timer timer;
939 timer.Start();
940 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
941 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
942 timer.Stop();
943 timer.AccumulateRegion("Pyr PhysDeriv SumFactor");
944
945 double epsilon = 1.0e-8;
946 for (int i = 0; i < diff1.size(); ++i)
947 {
948 // clamp values below 1e-14 to zero
949 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
950 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
951 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
952 }
953}
954
955BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_IterPerExp_VariableP_MultiElmt)
956{
957 SpatialDomains::PointGeomUniquePtr v0(
958 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
959 SpatialDomains::PointGeomUniquePtr v1(
960 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
961 SpatialDomains::PointGeomUniquePtr v2(
962 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
963 SpatialDomains::PointGeomUniquePtr v3(
964 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
965 SpatialDomains::PointGeomUniquePtr v4(
966 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
967
968 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
969 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
970 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
971 std::array<SpatialDomains::PointGeom *, 5> v = {
972 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
973 SpatialDomains::PyrGeomUniquePtr pyrGeom =
974 CreatePyr(v, segVec, triVec, quadVec);
975
976 Nektar::LibUtilities::PointsType PointsTypeDir1 =
977 Nektar::LibUtilities::eGaussLobattoLegendre;
978 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
979 Nektar::LibUtilities::BasisType basisTypeDir1 =
980 Nektar::LibUtilities::eModified_A;
981 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
982 PointsKeyDir1);
983
984 Nektar::LibUtilities::PointsType PointsTypeDir2 =
985 Nektar::LibUtilities::eGaussLobattoLegendre;
986 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
987 Nektar::LibUtilities::BasisType basisTypeDir2 =
988 Nektar::LibUtilities::eModified_A;
989 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
990 PointsKeyDir2);
991
992 Nektar::LibUtilities::PointsType PointsTypeDir3 =
993 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
994 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
995 Nektar::LibUtilities::BasisType basisTypeDir3 =
996 Nektar::LibUtilities::eModifiedPyr_C;
997 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
998 PointsKeyDir3);
999
1000 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1001 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1002 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1003
1004 int nelmts = 10;
1005
1006 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1007 for (int i = 0; i < nelmts; ++i)
1008 {
1009 CollExp.push_back(Exp);
1010 }
1011
1012 LibUtilities::SessionReaderSharedPtr dummySession;
1013 Collections::CollectionOptimisation colOpt(dummySession, 3,
1014 Collections::eIterPerExp);
1015 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1016 Collections::Collection c(CollExp, impTypes);
1017 c.Initialise(Collections::ePhysDeriv);
1018
1019 const int nq = Exp->GetTotPoints();
1020 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1021 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1022 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1023 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1024
1025 Exp->GetCoords(xc, yc, zc);
1026
1027 for (int i = 0; i < nq; ++i)
1028 {
1029 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1030 }
1031 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
1032 tmp2 = diff1 + (2 * nelmts) * nq);
1033
1034 for (int i = 1; i < nelmts; ++i)
1035 {
1036 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1037 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
1038 tmp1 = diff1 + (nelmts + i) * nq,
1039 tmp2 = diff1 + (2 * nelmts + i) * nq);
1040 }
1041
1042 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1043 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
1044
1045 double epsilon = 1.0e-8;
1046 for (int i = 0; i < diff1.size(); ++i)
1047 {
1048 // clamp values below 1e-14 to zero
1049 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1050 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1051 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1052 }
1053}
1054
1055BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_SumFac_VariableP_MultiElmt)
1056{
1057 SpatialDomains::PointGeomUniquePtr v0(
1058 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1059 SpatialDomains::PointGeomUniquePtr v1(
1060 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1061 SpatialDomains::PointGeomUniquePtr v2(
1062 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1063 SpatialDomains::PointGeomUniquePtr v3(
1064 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1065 SpatialDomains::PointGeomUniquePtr v4(
1066 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1067
1068 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1069 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1070 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1071 std::array<SpatialDomains::PointGeom *, 5> v = {
1072 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1073 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1074 CreatePyr(v, segVec, triVec, quadVec);
1075
1076 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1077 Nektar::LibUtilities::eGaussLobattoLegendre;
1078 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1079 Nektar::LibUtilities::BasisType basisTypeDir1 =
1080 Nektar::LibUtilities::eModified_A;
1081 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1082 PointsKeyDir1);
1083
1084 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1085 Nektar::LibUtilities::eGaussLobattoLegendre;
1086 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
1087 Nektar::LibUtilities::BasisType basisTypeDir2 =
1088 Nektar::LibUtilities::eModified_A;
1089 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1090 PointsKeyDir2);
1091
1092 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1093 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1094 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
1095 Nektar::LibUtilities::BasisType basisTypeDir3 =
1096 Nektar::LibUtilities::eModifiedPyr_C;
1097 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1098 PointsKeyDir3);
1099
1100 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1101 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1102 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1103
1104 int nelmts = 10;
1105
1106 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1107 for (int i = 0; i < nelmts; ++i)
1108 {
1109 CollExp.push_back(Exp);
1110 }
1111
1112 LibUtilities::SessionReaderSharedPtr dummySession;
1113 Collections::CollectionOptimisation colOpt(dummySession, 3,
1114 Collections::eSumFac);
1115 Collections::OperatorImpMap impTypes;
1116 impTypes[Collections::ePhysDeriv] = Collections::eSumFac;
1117
1118 Collections::Collection c(CollExp, impTypes);
1119 c.Initialise(Collections::ePhysDeriv);
1120
1121 const int nq = Exp->GetTotPoints();
1122 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1123 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1124 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1125 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1126
1127 Exp->GetCoords(xc, yc, zc);
1128
1129 for (int i = 0; i < nq; ++i)
1130 {
1131 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1132 }
1133 Exp->PhysDeriv(phys, tmp = diff1, tmp1 = diff1 + (nelmts)*nq,
1134 tmp2 = diff1 + (2 * nelmts) * nq);
1135 for (int i = 1; i < nelmts; ++i)
1136 {
1137 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1138 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
1139 tmp1 = diff1 + (nelmts + i) * nq,
1140 tmp2 = diff1 + (2 * nelmts + i) * nq);
1141 }
1142
1143 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1144 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
1145
1146 double epsilon = 1.0e-8;
1147 for (int i = 0; i < diff1.size(); ++i)
1148 {
1149 // clamp values below 1e-14 to zero
1150 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1151 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1152 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1153 }
1154}
1155
1156BOOST_AUTO_TEST_CASE(TestPyrPhysDeriv_MatrixFree_UniformP_MultiElmt)
1157{
1158 SpatialDomains::PointGeomUniquePtr v0(
1159 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1160 SpatialDomains::PointGeomUniquePtr v1(
1161 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1162 SpatialDomains::PointGeomUniquePtr v2(
1163 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1164 SpatialDomains::PointGeomUniquePtr v3(
1165 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1166 SpatialDomains::PointGeomUniquePtr v4(
1167 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1168
1169 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1170 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1171 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1172 std::array<SpatialDomains::PointGeom *, 5> v = {
1173 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1174 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1175 CreatePyr(v, segVec, triVec, quadVec);
1176
1177 unsigned int numQuadPoints = 5;
1178 unsigned int numModes = 2;
1179
1180 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1181 Nektar::LibUtilities::eGaussLobattoLegendre;
1182 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
1183 PointsTypeDir1);
1184 Nektar::LibUtilities::BasisType basisTypeDir1 =
1185 Nektar::LibUtilities::eModified_A;
1186 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1187 PointsKeyDir1);
1188
1189 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1190 Nektar::LibUtilities::eGaussLobattoLegendre;
1191 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
1192 PointsTypeDir2);
1193 Nektar::LibUtilities::BasisType basisTypeDir2 =
1194 Nektar::LibUtilities::eModified_A;
1195 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
1196 PointsKeyDir2);
1197
1198 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1199 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1200 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
1201 PointsTypeDir3);
1202 Nektar::LibUtilities::BasisType basisTypeDir3 =
1203 Nektar::LibUtilities::eModifiedPyr_C;
1204 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
1205 PointsKeyDir3);
1206
1207 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1208 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1209 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1210
1211 int nelmts = NELMTS;
1212
1213 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1214 for (int i = 0; i < nelmts; ++i)
1215 {
1216 CollExp.push_back(Exp);
1217 }
1218
1219 LibUtilities::SessionReaderSharedPtr dummySession;
1220 Collections::CollectionOptimisation colOpt(dummySession, 2,
1221 Collections::eIterPerExp);
1222 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1223
1224 impTypes[Collections::ePhysDeriv] = Collections::eMatrixFree;
1225 Collections::Collection c(CollExp, impTypes);
1226 c.Initialise(Collections::ePhysDeriv);
1227
1228 const int nq = Exp->GetTotPoints();
1229 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1230 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1231 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1232 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1233
1234 Exp->GetCoords(xc, yc, zc);
1235
1236 for (int i = 0; i < nq; ++i)
1237 {
1238 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1239 }
1240 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
1241 tmp2 = diff1 + (2 * nelmts) * nq);
1242
1243 for (int i = 1; i < nelmts; ++i)
1244 {
1245 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1246 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
1247 tmp1 = diff1 + (nelmts + i) * nq,
1248 tmp2 = diff1 + (2 * nelmts + i) * nq);
1249 }
1250
1251 LibUtilities::Timer timer;
1252 timer.Start();
1253 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1254 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
1255 timer.Stop();
1256 timer.AccumulateRegion("Pyr PhysDeriv MatrixFree");
1257
1258 double epsilon = 1.0e-8;
1259 for (int i = 0; i < diff1.size(); ++i)
1260 {
1261 // clamp values below 1e-14 to zero
1262 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1263 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1264 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1265 }
1266}
1267
1268BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_IterPerExp_UniformP_MultiElmt)
1269{
1270 SpatialDomains::PointGeomUniquePtr v0(
1271 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1272 SpatialDomains::PointGeomUniquePtr v1(
1273 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1274 SpatialDomains::PointGeomUniquePtr v2(
1275 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1276 SpatialDomains::PointGeomUniquePtr v3(
1277 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1278 SpatialDomains::PointGeomUniquePtr v4(
1279 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1280
1281 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1282 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1283 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1284 std::array<SpatialDomains::PointGeom *, 5> v = {
1285 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1286 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1287 CreatePyr(v, segVec, triVec, quadVec);
1288
1289 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1290 Nektar::LibUtilities::eGaussLobattoLegendre;
1291 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1292 Nektar::LibUtilities::BasisType basisTypeDir1 =
1293 Nektar::LibUtilities::eModified_A;
1294 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1295 PointsKeyDir1);
1296
1297 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1298 Nektar::LibUtilities::eGaussLobattoLegendre;
1299 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
1300 Nektar::LibUtilities::BasisType basisTypeDir2 =
1301 Nektar::LibUtilities::eModified_A;
1302 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1303 PointsKeyDir2);
1304
1305 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1306 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1307 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
1308 Nektar::LibUtilities::BasisType basisTypeDir3 =
1309 Nektar::LibUtilities::eModifiedPyr_C;
1310 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1311 PointsKeyDir3);
1312
1313 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1314 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1315 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1316
1317 int nelmts = NELMTS;
1318
1319 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1320 for (int i = 0; i < nelmts; ++i)
1321 {
1322 CollExp.push_back(Exp);
1323 }
1324
1325 LibUtilities::SessionReaderSharedPtr dummySession;
1326 Collections::CollectionOptimisation colOpt(dummySession, 3,
1327 Collections::eIterPerExp);
1328 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1329 Collections::Collection c(CollExp, impTypes);
1330 c.Initialise(Collections::eIProductWRTBase);
1331
1332 const int nq = Exp->GetTotPoints();
1333 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1334 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1335 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1336 Array<OneD, NekDouble> tmp;
1337 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1338
1339 Exp->GetCoords(xc, yc, zc);
1340
1341 for (int i = 0; i < nq; ++i)
1342 {
1343 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1344 }
1345 Exp->IProductWRTBase(phys, coeffs1);
1346
1347 for (int i = 1; i < nelmts; ++i)
1348 {
1349 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1350 Exp->IProductWRTBase(phys + i * nq,
1351 tmp = coeffs1 + i * Exp->GetNcoeffs());
1352 }
1353
1354 LibUtilities::Timer timer;
1355 timer.Start();
1356 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1357 timer.Stop();
1358 timer.AccumulateRegion("Pyr IProdWRTB IterPerExp");
1359
1360 double epsilon = 1.0e-8;
1361 for (int i = 0; i < coeffs1.size(); ++i)
1362 {
1363 // clamp values below 1e-14 to zero
1364 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1365 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1366 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1367 }
1368}
1369
1370BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_StdMat_UniformP_MultiElmt)
1371{
1372 SpatialDomains::PointGeomUniquePtr v0(
1373 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1374 SpatialDomains::PointGeomUniquePtr v1(
1375 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1376 SpatialDomains::PointGeomUniquePtr v2(
1377 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1378 SpatialDomains::PointGeomUniquePtr v3(
1379 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1380 SpatialDomains::PointGeomUniquePtr v4(
1381 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1382
1383 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1384 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1385 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1386 std::array<SpatialDomains::PointGeom *, 5> v = {
1387 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1388 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1389 CreatePyr(v, segVec, triVec, quadVec);
1390
1391 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1392 Nektar::LibUtilities::eGaussLobattoLegendre;
1393 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1394 Nektar::LibUtilities::BasisType basisTypeDir1 =
1395 Nektar::LibUtilities::eModified_A;
1396 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1397 PointsKeyDir1);
1398
1399 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1400 Nektar::LibUtilities::eGaussLobattoLegendre;
1401 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
1402 Nektar::LibUtilities::BasisType basisTypeDir2 =
1403 Nektar::LibUtilities::eModified_A;
1404 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1405 PointsKeyDir2);
1406
1407 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1408 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1409 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
1410 Nektar::LibUtilities::BasisType basisTypeDir3 =
1411 Nektar::LibUtilities::eModifiedPyr_C;
1412 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1413 PointsKeyDir3);
1414
1415 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1416 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1417 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1418
1419 int nelmts = 10;
1420
1421 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1422 for (int i = 0; i < nelmts; ++i)
1423 {
1424 CollExp.push_back(Exp);
1425 }
1426
1427 LibUtilities::SessionReaderSharedPtr dummySession;
1428 Collections::CollectionOptimisation colOpt(dummySession, 3,
1429 Collections::eStdMat);
1430 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1431 Collections::Collection c(CollExp, impTypes);
1432 c.Initialise(Collections::eIProductWRTBase);
1433
1434 const int nq = Exp->GetTotPoints();
1435 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1436 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1437 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1438 Array<OneD, NekDouble> tmp;
1439 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1440
1441 Exp->GetCoords(xc, yc, zc);
1442
1443 for (int i = 0; i < nq; ++i)
1444 {
1445 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1446 }
1447 Exp->IProductWRTBase(phys, coeffs1);
1448
1449 for (int i = 1; i < nelmts; ++i)
1450 {
1451 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1452 Exp->IProductWRTBase(phys + i * nq,
1453 tmp = coeffs1 + i * Exp->GetNcoeffs());
1454 }
1455
1456 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1457
1458 double epsilon = 1.0e-8;
1459 for (int i = 0; i < coeffs1.size(); ++i)
1460 {
1461 // clamp values below 1e-14 to zero
1462 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1463 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1464 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1465 }
1466}
1467
1468BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_SumFac_UniformP_MultiElmt)
1469{
1470 SpatialDomains::PointGeomUniquePtr v0(
1471 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1472 SpatialDomains::PointGeomUniquePtr v1(
1473 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1474 SpatialDomains::PointGeomUniquePtr v2(
1475 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1476 SpatialDomains::PointGeomUniquePtr v3(
1477 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1478 SpatialDomains::PointGeomUniquePtr v4(
1479 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1480
1481 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1482 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1483 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1484 std::array<SpatialDomains::PointGeom *, 5> v = {
1485 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1486 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1487 CreatePyr(v, segVec, triVec, quadVec);
1488
1489 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1490 Nektar::LibUtilities::eGaussLobattoLegendre;
1491 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1492 Nektar::LibUtilities::BasisType basisTypeDir1 =
1493 Nektar::LibUtilities::eModified_A;
1494 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1495 PointsKeyDir1);
1496
1497 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1498 Nektar::LibUtilities::eGaussLobattoLegendre;
1499 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
1500 Nektar::LibUtilities::BasisType basisTypeDir2 =
1501 Nektar::LibUtilities::eModified_A;
1502 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1503 PointsKeyDir2);
1504
1505 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1506 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1507 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
1508 Nektar::LibUtilities::BasisType basisTypeDir3 =
1509 Nektar::LibUtilities::eModifiedPyr_C;
1510 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1511 PointsKeyDir3);
1512
1513 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1514 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1515 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1516
1517 int nelmts = NELMTS;
1518
1519 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1520 for (int i = 0; i < nelmts; ++i)
1521 {
1522 CollExp.push_back(Exp);
1523 }
1524
1525 LibUtilities::SessionReaderSharedPtr dummySession;
1526 Collections::CollectionOptimisation colOpt(dummySession, 3,
1527 Collections::eSumFac);
1528 Collections::OperatorImpMap impTypes;
1529 impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
1530 Collections::Collection c(CollExp, impTypes);
1531 c.Initialise(Collections::eIProductWRTBase);
1532
1533 const int nq = Exp->GetTotPoints();
1534 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1535 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1536 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1537 Array<OneD, NekDouble> tmp;
1538 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1539
1540 Exp->GetCoords(xc, yc, zc);
1541
1542 for (int i = 0; i < nq; ++i)
1543 {
1544 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1545 }
1546 Exp->IProductWRTBase(phys, coeffs1);
1547
1548 for (int i = 1; i < nelmts; ++i)
1549 {
1550 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1551 Exp->IProductWRTBase(phys + i * nq,
1552 tmp = coeffs1 + i * Exp->GetNcoeffs());
1553 }
1554
1555 LibUtilities::Timer timer;
1556 timer.Start();
1557 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1558 timer.Stop();
1559 timer.AccumulateRegion("Pyr IProdWRTB SumFactor");
1560
1561 double epsilon = 1.0e-8;
1562 for (int i = 0; i < coeffs1.size(); ++i)
1563 {
1564 // clamp values below 1e-14 to zero
1565 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1566 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1567 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1568 }
1569}
1570
1571BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_IterPerExp_VariableP_MultiElmt)
1572{
1573 SpatialDomains::PointGeomUniquePtr v0(
1574 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1575 SpatialDomains::PointGeomUniquePtr v1(
1576 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1577 SpatialDomains::PointGeomUniquePtr v2(
1578 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1579 SpatialDomains::PointGeomUniquePtr v3(
1580 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1581 SpatialDomains::PointGeomUniquePtr v4(
1582 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1583
1584 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1585 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1586 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1587 std::array<SpatialDomains::PointGeom *, 5> v = {
1588 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1589 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1590 CreatePyr(v, segVec, triVec, quadVec);
1591
1592 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1593 Nektar::LibUtilities::eGaussLobattoLegendre;
1594 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1595 Nektar::LibUtilities::BasisType basisTypeDir1 =
1596 Nektar::LibUtilities::eModified_A;
1597 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1598 PointsKeyDir1);
1599
1600 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1601 Nektar::LibUtilities::eGaussLobattoLegendre;
1602 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
1603 Nektar::LibUtilities::BasisType basisTypeDir2 =
1604 Nektar::LibUtilities::eModified_A;
1605 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1606 PointsKeyDir2);
1607
1608 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1609 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1610 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
1611 Nektar::LibUtilities::BasisType basisTypeDir3 =
1612 Nektar::LibUtilities::eModifiedPyr_C;
1613 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1614 PointsKeyDir3);
1615
1616 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1617 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1618 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1619
1620 int nelmts = 10;
1621
1622 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1623 for (int i = 0; i < nelmts; ++i)
1624 {
1625 CollExp.push_back(Exp);
1626 }
1627
1628 LibUtilities::SessionReaderSharedPtr dummySession;
1629 Collections::CollectionOptimisation colOpt(dummySession, 3,
1630 Collections::eIterPerExp);
1631 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1632 Collections::Collection c(CollExp, impTypes);
1633 c.Initialise(Collections::eIProductWRTBase);
1634
1635 const int nq = Exp->GetTotPoints();
1636 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1637 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1638 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1639 Array<OneD, NekDouble> tmp;
1640 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1641
1642 Exp->GetCoords(xc, yc, zc);
1643
1644 for (int i = 0; i < nq; ++i)
1645 {
1646 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1647 }
1648 Exp->IProductWRTBase(phys, coeffs1);
1649
1650 for (int i = 1; i < nelmts; ++i)
1651 {
1652 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1653 Exp->IProductWRTBase(phys + i * nq,
1654 tmp = coeffs1 + i * Exp->GetNcoeffs());
1655 }
1656
1657 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1658
1659 double epsilon = 1.0e-8;
1660 for (int i = 0; i < coeffs1.size(); ++i)
1661 {
1662 // clamp values below 1e-14 to zero
1663 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1664 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1665 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1666 }
1667}
1668
1669BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_StdMat_VariableP_MultiElmt)
1670{
1671 SpatialDomains::PointGeomUniquePtr v0(
1672 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1673 SpatialDomains::PointGeomUniquePtr v1(
1674 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1675 SpatialDomains::PointGeomUniquePtr v2(
1676 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1677 SpatialDomains::PointGeomUniquePtr v3(
1678 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1679 SpatialDomains::PointGeomUniquePtr v4(
1680 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1681
1682 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1683 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1684 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1685 std::array<SpatialDomains::PointGeom *, 5> v = {
1686 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1687 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1688 CreatePyr(v, segVec, triVec, quadVec);
1689
1690 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1691 Nektar::LibUtilities::eGaussLobattoLegendre;
1692 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1693 Nektar::LibUtilities::BasisType basisTypeDir1 =
1694 Nektar::LibUtilities::eModified_A;
1695 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1696 PointsKeyDir1);
1697
1698 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1699 Nektar::LibUtilities::eGaussLobattoLegendre;
1700 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
1701 Nektar::LibUtilities::BasisType basisTypeDir2 =
1702 Nektar::LibUtilities::eModified_A;
1703 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1704 PointsKeyDir2);
1705
1706 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1707 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1708 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
1709 Nektar::LibUtilities::BasisType basisTypeDir3 =
1710 Nektar::LibUtilities::eModifiedPyr_C;
1711 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1712 PointsKeyDir3);
1713
1714 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1715 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1716 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1717
1718 int nelmts = 10;
1719
1720 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1721 for (int i = 0; i < nelmts; ++i)
1722 {
1723 CollExp.push_back(Exp);
1724 }
1725
1726 LibUtilities::SessionReaderSharedPtr dummySession;
1727 Collections::CollectionOptimisation colOpt(dummySession, 3,
1728 Collections::eStdMat);
1729 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1730 Collections::Collection c(CollExp, impTypes);
1731 c.Initialise(Collections::eIProductWRTBase);
1732
1733 const int nq = Exp->GetTotPoints();
1734 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1735 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1736 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1737 Array<OneD, NekDouble> tmp;
1738 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1739
1740 Exp->GetCoords(xc, yc, zc);
1741
1742 for (int i = 0; i < nq; ++i)
1743 {
1744 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1745 }
1746 Exp->IProductWRTBase(phys, coeffs1);
1747
1748 for (int i = 1; i < nelmts; ++i)
1749 {
1750 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1751 Exp->IProductWRTBase(phys + i * nq,
1752 tmp = coeffs1 + i * Exp->GetNcoeffs());
1753 }
1754
1755 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1756
1757 double epsilon = 1.0e-8;
1758 for (int i = 0; i < coeffs1.size(); ++i)
1759 {
1760 // clamp values below 1e-14 to zero
1761 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1762 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1763 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1764 }
1765}
1766
1767BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_SumFac_VariableP_MultiElmt)
1768{
1769 SpatialDomains::PointGeomUniquePtr v0(
1770 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1771 SpatialDomains::PointGeomUniquePtr v1(
1772 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1773 SpatialDomains::PointGeomUniquePtr v2(
1774 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1775 SpatialDomains::PointGeomUniquePtr v3(
1776 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1777 SpatialDomains::PointGeomUniquePtr v4(
1778 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1779
1780 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1781 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1782 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1783 std::array<SpatialDomains::PointGeom *, 5> v = {
1784 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1785 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1786 CreatePyr(v, segVec, triVec, quadVec);
1787
1788 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1789 Nektar::LibUtilities::eGaussLobattoLegendre;
1790 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
1791 Nektar::LibUtilities::BasisType basisTypeDir1 =
1792 Nektar::LibUtilities::eModified_A;
1793 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1794 PointsKeyDir1);
1795
1796 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1797 Nektar::LibUtilities::eGaussLobattoLegendre;
1798 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
1799 Nektar::LibUtilities::BasisType basisTypeDir2 =
1800 Nektar::LibUtilities::eModified_A;
1801 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1802 PointsKeyDir2);
1803
1804 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1805 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1806 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
1807 Nektar::LibUtilities::BasisType basisTypeDir3 =
1808 Nektar::LibUtilities::eModifiedPyr_C;
1809 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1810 PointsKeyDir3);
1811
1812 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1813 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1814 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1815
1816 int nelmts = 10;
1817
1818 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1819 for (int i = 0; i < nelmts; ++i)
1820 {
1821 CollExp.push_back(Exp);
1822 }
1823
1824 LibUtilities::SessionReaderSharedPtr dummySession;
1825 Collections::CollectionOptimisation colOpt(dummySession, 3,
1826 Collections::eSumFac);
1827 Collections::OperatorImpMap impTypes;
1828 impTypes[Collections::eIProductWRTBase] = Collections::eSumFac;
1829 Collections::Collection c(CollExp, impTypes);
1830 c.Initialise(Collections::eIProductWRTBase);
1831
1832 const int nq = Exp->GetTotPoints();
1833 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1834 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1835 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1836 Array<OneD, NekDouble> tmp;
1837 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1838
1839 Exp->GetCoords(xc, yc, zc);
1840
1841 for (int i = 0; i < nq; ++i)
1842 {
1843 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1844 }
1845 Exp->IProductWRTBase(phys, coeffs1);
1846
1847 for (int i = 1; i < nelmts; ++i)
1848 {
1849 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1850 Exp->IProductWRTBase(phys + i * nq,
1851 tmp = coeffs1 + i * Exp->GetNcoeffs());
1852 }
1853
1854 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1855
1856 double epsilon = 1.0e-8;
1857 for (int i = 0; i < coeffs1.size(); ++i)
1858 {
1859 // clamp values below 1e-14 to zero
1860 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1861 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1862 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1863 }
1864}
1865
1866BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_MatrixFree_UniformP_MultiElmt)
1867{
1868 SpatialDomains::PointGeomUniquePtr v0(
1869 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1870 SpatialDomains::PointGeomUniquePtr v1(
1871 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1872 SpatialDomains::PointGeomUniquePtr v2(
1873 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1874 SpatialDomains::PointGeomUniquePtr v3(
1875 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1876 SpatialDomains::PointGeomUniquePtr v4(
1877 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1878
1879 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1880 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1881 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1882 std::array<SpatialDomains::PointGeom *, 5> v = {
1883 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1884 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1885 CreatePyr(v, segVec, triVec, quadVec);
1886
1887 unsigned int numQuadPoints = 5;
1888 unsigned int numModes = 4;
1889
1890 Nektar::LibUtilities::PointsType PointsTypeDir1 =
1891 Nektar::LibUtilities::eGaussLobattoLegendre;
1892 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
1893 PointsTypeDir1);
1894 Nektar::LibUtilities::BasisType basisTypeDir1 =
1895 Nektar::LibUtilities::eModified_A;
1896 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1897 PointsKeyDir1);
1898
1899 Nektar::LibUtilities::PointsType PointsTypeDir2 =
1900 Nektar::LibUtilities::eGaussLobattoLegendre;
1901 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
1902 PointsTypeDir2);
1903 Nektar::LibUtilities::BasisType basisTypeDir2 =
1904 Nektar::LibUtilities::eModified_A;
1905 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
1906 PointsKeyDir2);
1907
1908 Nektar::LibUtilities::PointsType PointsTypeDir3 =
1909 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1910 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
1911 PointsTypeDir3);
1912 Nektar::LibUtilities::BasisType basisTypeDir3 =
1913 Nektar::LibUtilities::eModifiedPyr_C;
1914 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
1915 PointsKeyDir3);
1916
1917 Nektar::LocalRegions::PyrExpSharedPtr Exp =
1918 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
1919 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
1920
1921 int nelmts = NELMTS;
1922
1923 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
1924 for (int i = 0; i < nelmts; ++i)
1925 {
1926 CollExp.push_back(Exp);
1927 }
1928
1929 LibUtilities::SessionReaderSharedPtr dummySession;
1930 Collections::CollectionOptimisation colOpt(dummySession, 2,
1931 Collections::eIterPerExp);
1932 Collections::OperatorImpMap impTypes;
1933 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
1934 Collections::Collection c(CollExp, impTypes);
1935 c.Initialise(Collections::eIProductWRTBase);
1936
1937 const int nq = Exp->GetTotPoints();
1938 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
1939 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
1940 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
1941 Array<OneD, NekDouble> tmp;
1942 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1943
1944 Exp->GetCoords(xc, yc, zc);
1945
1946 for (int i = 0; i < nq; ++i)
1947 {
1948 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1949 }
1950 Exp->IProductWRTBase(phys, coeffs1);
1951
1952 for (int i = 1; i < nelmts; ++i)
1953 {
1954 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1955 Exp->IProductWRTBase(phys + i * nq,
1956 tmp = coeffs1 + i * Exp->GetNcoeffs());
1957 }
1958
1959 LibUtilities::Timer timer;
1960 timer.Start();
1961 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1962 timer.Stop();
1963 timer.AccumulateRegion("Pyr IProdWRTB MatrixFree");
1964
1965 double epsilon = 1.0e-8;
1966 for (int i = 0; i < coeffs1.size(); ++i)
1967 {
1968 // clamp values below 1e-14 to zero
1969 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1970 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1971 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1972 }
1973}
1974
1975BOOST_AUTO_TEST_CASE(TestPyrIProductWRTBase_MatrixFree_Deformed_MultiElmt)
1976{
1977 SpatialDomains::PointGeomUniquePtr v0(
1978 new SpatialDomains::PointGeom(3u, 0u, -2.0, -3.0, -4.0));
1979 SpatialDomains::PointGeomUniquePtr v1(
1980 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1981 SpatialDomains::PointGeomUniquePtr v2(
1982 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
1983 SpatialDomains::PointGeomUniquePtr v3(
1984 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
1985 SpatialDomains::PointGeomUniquePtr v4(
1986 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
1987
1988 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
1989 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
1990 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
1991 std::array<SpatialDomains::PointGeom *, 5> v = {
1992 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
1993 SpatialDomains::PyrGeomUniquePtr pyrGeom =
1994 CreatePyr(v, segVec, triVec, quadVec);
1995
1996 unsigned int numQuadPoints = 5;
1997 unsigned int numModes = 4;
1998
1999 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2000 Nektar::LibUtilities::eGaussLobattoLegendre;
2001 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
2002 PointsTypeDir1);
2003 Nektar::LibUtilities::BasisType basisTypeDir1 =
2004 Nektar::LibUtilities::eModified_A;
2005 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2006 PointsKeyDir1);
2007
2008 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2009 Nektar::LibUtilities::eGaussLobattoLegendre;
2010 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
2011 PointsTypeDir2);
2012 Nektar::LibUtilities::BasisType basisTypeDir2 =
2013 Nektar::LibUtilities::eModified_A;
2014 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2015 PointsKeyDir2);
2016
2017 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2018 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2019 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
2020 PointsTypeDir3);
2021 Nektar::LibUtilities::BasisType basisTypeDir3 =
2022 Nektar::LibUtilities::eModifiedPyr_C;
2023 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2024 PointsKeyDir3);
2025
2026 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2027 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2028 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2029
2030 int nelmts = 5;
2031
2032 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2033 for (int i = 0; i < nelmts; ++i)
2034 {
2035 CollExp.push_back(Exp);
2036 }
2037
2038 LibUtilities::SessionReaderSharedPtr dummySession;
2039 Collections::CollectionOptimisation colOpt(dummySession, 2,
2040 Collections::eIterPerExp);
2041 Collections::OperatorImpMap impTypes;
2042
2043 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
2044 Collections::Collection c(CollExp, impTypes);
2045 c.Initialise(Collections::eIProductWRTBase);
2046
2047 const int nq = Exp->GetTotPoints();
2048 Array<OneD, NekDouble> phys(nelmts * nq, 0.0);
2049 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs(), 0.0);
2050 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs(), 0.0);
2051 Array<OneD, NekDouble> tmp;
2052 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2053
2054 Exp->GetCoords(xc, yc, zc);
2055
2056 for (int i = 0; i < nq; ++i)
2057 {
2058 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2059 }
2060 Exp->IProductWRTBase(phys, coeffs1);
2061
2062 for (int i = 1; i < nelmts; ++i)
2063 {
2064 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
2065 Exp->IProductWRTBase(phys + i * nq,
2066 tmp = coeffs1 + i * Exp->GetNcoeffs());
2067 }
2068
2069 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
2070
2071 double epsilon = 1.0e-8;
2072 for (int i = 0; i < coeffs1.size(); ++i)
2073 {
2074 // clamp values below 1e-14 to zero
2075 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2076 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2077 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2078 }
2079}
2080
2081BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_IterPerExp_UniformP_MultiElmt)
2082{
2083 SpatialDomains::PointGeomUniquePtr v0(
2084 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2085 SpatialDomains::PointGeomUniquePtr v1(
2086 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2087 SpatialDomains::PointGeomUniquePtr v2(
2088 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2089 SpatialDomains::PointGeomUniquePtr v3(
2090 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2091 SpatialDomains::PointGeomUniquePtr v4(
2092 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2093
2094 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2095 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2096 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2097 std::array<SpatialDomains::PointGeom *, 5> v = {
2098 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2099 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2100 CreatePyr(v, segVec, triVec, quadVec);
2101
2102 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2103 Nektar::LibUtilities::eGaussLobattoLegendre;
2104 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
2105 Nektar::LibUtilities::BasisType basisTypeDir1 =
2106 Nektar::LibUtilities::eModified_A;
2107 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2108 PointsKeyDir1);
2109
2110 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2111 Nektar::LibUtilities::eGaussLobattoLegendre;
2112 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
2113 Nektar::LibUtilities::BasisType basisTypeDir2 =
2114 Nektar::LibUtilities::eModified_A;
2115 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2116 PointsKeyDir2);
2117
2118 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2119 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2120 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
2121 Nektar::LibUtilities::BasisType basisTypeDir3 =
2122 Nektar::LibUtilities::eModifiedPyr_C;
2123 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2124 PointsKeyDir3);
2125
2126 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2127 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2128 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2129
2130 int nelmts = NELMTS;
2131
2132 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2133 for (int i = 0; i < nelmts; ++i)
2134 {
2135 CollExp.push_back(Exp);
2136 }
2137
2138 LibUtilities::SessionReaderSharedPtr dummySession;
2139 Collections::CollectionOptimisation colOpt(dummySession, 3,
2140 Collections::eIterPerExp);
2141 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2142 Collections::Collection c(CollExp, impTypes);
2143 c.Initialise(Collections::eIProductWRTDerivBase);
2144
2145 const int nq = Exp->GetTotPoints();
2146 const int nm = Exp->GetNcoeffs();
2147 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2148 Array<OneD, NekDouble> phys2(nelmts * nq);
2149 Array<OneD, NekDouble> phys3(nelmts * nq);
2150 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2151 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2152
2153 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2154
2155 Exp->GetCoords(xc, yc, zc);
2156
2157 for (int i = 0; i < nq; ++i)
2158 {
2159 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2160 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2161 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2162 }
2163 for (int i = 1; i < nelmts; ++i)
2164 {
2165 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2166 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2167 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2168 }
2169
2170 // Standard routines
2171 for (int i = 0; i < nelmts; ++i)
2172 {
2173 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2174 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2175 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2176 tmp = coeffs1 + i * nm, 1);
2177 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2178 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2179 tmp = coeffs1 + i * nm, 1);
2180 }
2181
2182 LibUtilities::Timer timer;
2183 timer.Start();
2184 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2185 coeffs2);
2186 timer.Stop();
2187 timer.AccumulateRegion("Pyr IProdWRTDB IterPerExp");
2188
2189 double epsilon = 1.0e-8;
2190 for (int i = 0; i < coeffs1.size(); ++i)
2191 {
2192 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2193 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2194 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2195 }
2196}
2197
2198BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_StdMat_UniformP_MultiElmt)
2199{
2200 SpatialDomains::PointGeomUniquePtr v0(
2201 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2202 SpatialDomains::PointGeomUniquePtr v1(
2203 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2204 SpatialDomains::PointGeomUniquePtr v2(
2205 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2206 SpatialDomains::PointGeomUniquePtr v3(
2207 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2208 SpatialDomains::PointGeomUniquePtr v4(
2209 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2210
2211 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2212 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2213 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2214 std::array<SpatialDomains::PointGeom *, 5> v = {
2215 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2216 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2217 CreatePyr(v, segVec, triVec, quadVec);
2218
2219 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2220 Nektar::LibUtilities::eGaussLobattoLegendre;
2221 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
2222 Nektar::LibUtilities::BasisType basisTypeDir1 =
2223 Nektar::LibUtilities::eModified_A;
2224 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2225 PointsKeyDir1);
2226
2227 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2228 Nektar::LibUtilities::eGaussLobattoLegendre;
2229 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
2230 Nektar::LibUtilities::BasisType basisTypeDir2 =
2231 Nektar::LibUtilities::eModified_A;
2232 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2233 PointsKeyDir2);
2234
2235 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2236 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2237 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
2238 Nektar::LibUtilities::BasisType basisTypeDir3 =
2239 Nektar::LibUtilities::eModifiedPyr_C;
2240 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2241 PointsKeyDir3);
2242
2243 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2244 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2245 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2246
2247 int nelmts = 10;
2248
2249 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2250 for (int i = 0; i < nelmts; ++i)
2251 {
2252 CollExp.push_back(Exp);
2253 }
2254
2255 LibUtilities::SessionReaderSharedPtr dummySession;
2256 Collections::CollectionOptimisation colOpt(dummySession, 3,
2257 Collections::eStdMat);
2258 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2259 Collections::Collection c(CollExp, impTypes);
2260 c.Initialise(Collections::eIProductWRTDerivBase);
2261
2262 const int nq = Exp->GetTotPoints();
2263 const int nm = Exp->GetNcoeffs();
2264 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2265 Array<OneD, NekDouble> phys2(nelmts * nq);
2266 Array<OneD, NekDouble> phys3(nelmts * nq);
2267 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2268 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2269
2270 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2271
2272 Exp->GetCoords(xc, yc, zc);
2273
2274 for (int i = 0; i < nq; ++i)
2275 {
2276 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2277 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2278 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2279 }
2280 for (int i = 1; i < nelmts; ++i)
2281 {
2282 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2283 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2284 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2285 }
2286
2287 // Standard routines
2288 for (int i = 0; i < nelmts; ++i)
2289 {
2290 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2291 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2292 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2293 tmp = coeffs1 + i * nm, 1);
2294 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2295 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2296 tmp = coeffs1 + i * nm, 1);
2297 }
2298
2299 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2300 coeffs2);
2301
2302 double epsilon = 1.0e-8;
2303 for (int i = 0; i < coeffs1.size(); ++i)
2304 {
2305 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2306 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2307 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2308 }
2309}
2310
2311BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_SumFac_UniformP_MultiElmt)
2312{
2313 SpatialDomains::PointGeomUniquePtr v0(
2314 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2315 SpatialDomains::PointGeomUniquePtr v1(
2316 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2317 SpatialDomains::PointGeomUniquePtr v2(
2318 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2319 SpatialDomains::PointGeomUniquePtr v3(
2320 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2321 SpatialDomains::PointGeomUniquePtr v4(
2322 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2323
2324 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2325 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2326 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2327 std::array<SpatialDomains::PointGeom *, 5> v = {
2328 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2329 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2330 CreatePyr(v, segVec, triVec, quadVec);
2331
2332 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2333 Nektar::LibUtilities::eGaussLobattoLegendre;
2334 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
2335 Nektar::LibUtilities::BasisType basisTypeDir1 =
2336 Nektar::LibUtilities::eModified_A;
2337 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2338 PointsKeyDir1);
2339
2340 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2341 Nektar::LibUtilities::eGaussLobattoLegendre;
2342 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
2343 Nektar::LibUtilities::BasisType basisTypeDir2 =
2344 Nektar::LibUtilities::eModified_A;
2345 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2346 PointsKeyDir2);
2347
2348 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2349 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2350 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
2351 Nektar::LibUtilities::BasisType basisTypeDir3 =
2352 Nektar::LibUtilities::eModifiedPyr_C;
2353 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2354 PointsKeyDir3);
2355
2356 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2357 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2358 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2359
2360 int nelmts = NELMTS;
2361
2362 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2363 for (int i = 0; i < nelmts; ++i)
2364 {
2365 CollExp.push_back(Exp);
2366 }
2367
2368 LibUtilities::SessionReaderSharedPtr dummySession;
2369 Collections::CollectionOptimisation colOpt(dummySession, 3,
2370 Collections::eSumFac);
2371 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2372 Collections::Collection c(CollExp, impTypes);
2373 c.Initialise(Collections::eIProductWRTDerivBase);
2374
2375 const int nq = Exp->GetTotPoints();
2376 const int nm = Exp->GetNcoeffs();
2377 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2378 Array<OneD, NekDouble> phys2(nelmts * nq);
2379 Array<OneD, NekDouble> phys3(nelmts * nq);
2380 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2381 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2382
2383 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2384
2385 Exp->GetCoords(xc, yc, zc);
2386
2387 for (int i = 0; i < nq; ++i)
2388 {
2389 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2390 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2391 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2392 }
2393 for (int i = 1; i < nelmts; ++i)
2394 {
2395 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2396 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2397 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2398 }
2399
2400 // Standard routines
2401 for (int i = 0; i < nelmts; ++i)
2402 {
2403 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2404 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2405 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2406 tmp = coeffs1 + i * nm, 1);
2407 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2408 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2409 tmp = coeffs1 + i * nm, 1);
2410 }
2411
2412 LibUtilities::Timer timer;
2413 timer.Start();
2414 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2415 coeffs2);
2416 timer.Stop();
2417 timer.AccumulateRegion("Pyr IProdWRTDB SumFactor");
2418
2419 double epsilon = 1.0e-8;
2420 for (int i = 0; i < coeffs1.size(); ++i)
2421 {
2422 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2423 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2424 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2425 }
2426}
2427
2428BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_IterPerExp_VariableP_MultiElmt)
2429{
2430 SpatialDomains::PointGeomUniquePtr v0(
2431 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2432 SpatialDomains::PointGeomUniquePtr v1(
2433 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2434 SpatialDomains::PointGeomUniquePtr v2(
2435 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2436 SpatialDomains::PointGeomUniquePtr v3(
2437 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2438 SpatialDomains::PointGeomUniquePtr v4(
2439 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2440
2441 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2442 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2443 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2444 std::array<SpatialDomains::PointGeom *, 5> v = {
2445 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2446 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2447 CreatePyr(v, segVec, triVec, quadVec);
2448
2449 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2450 Nektar::LibUtilities::eGaussLobattoLegendre;
2451 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
2452 Nektar::LibUtilities::BasisType basisTypeDir1 =
2453 Nektar::LibUtilities::eModified_A;
2454 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2455 PointsKeyDir1);
2456
2457 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2458 Nektar::LibUtilities::eGaussLobattoLegendre;
2459 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
2460 Nektar::LibUtilities::BasisType basisTypeDir2 =
2461 Nektar::LibUtilities::eModified_A;
2462 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
2463 PointsKeyDir2);
2464
2465 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2466 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2467 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
2468 Nektar::LibUtilities::BasisType basisTypeDir3 =
2469 Nektar::LibUtilities::eModifiedPyr_C;
2470 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
2471 PointsKeyDir3);
2472
2473 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2474 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2475 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2476
2477 int nelmts = 10;
2478
2479 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2480 for (int i = 0; i < nelmts; ++i)
2481 {
2482 CollExp.push_back(Exp);
2483 }
2484
2485 LibUtilities::SessionReaderSharedPtr dummySession;
2486 Collections::CollectionOptimisation colOpt(dummySession, 3,
2487 Collections::eIterPerExp);
2488 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2489 Collections::Collection c(CollExp, impTypes);
2490 c.Initialise(Collections::eIProductWRTDerivBase);
2491
2492 const int nq = Exp->GetTotPoints();
2493 const int nm = Exp->GetNcoeffs();
2494 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2495 Array<OneD, NekDouble> phys2(nelmts * nq);
2496 Array<OneD, NekDouble> phys3(nelmts * nq);
2497 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2498 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2499
2500 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2501
2502 Exp->GetCoords(xc, yc, zc);
2503
2504 for (int i = 0; i < nq; ++i)
2505 {
2506 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2507 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2508 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2509 }
2510 for (int i = 1; i < nelmts; ++i)
2511 {
2512 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2513 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2514 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2515 }
2516
2517 // Standard routines
2518 for (int i = 0; i < nelmts; ++i)
2519 {
2520 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2521 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2522 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2523 tmp = coeffs1 + i * nm, 1);
2524 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2525 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2526 tmp = coeffs1 + i * nm, 1);
2527 }
2528
2529 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2530 coeffs2);
2531
2532 double epsilon = 1.0e-8;
2533 for (int i = 0; i < coeffs1.size(); ++i)
2534 {
2535 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2536 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2537 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2538 }
2539}
2540
2541BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_MatrixFree_Deformed_MultiElmt)
2542{
2543 SpatialDomains::PointGeomUniquePtr v0(
2544 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2545 SpatialDomains::PointGeomUniquePtr v1(
2546 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2547 SpatialDomains::PointGeomUniquePtr v2(
2548 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2549 SpatialDomains::PointGeomUniquePtr v3(
2550 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2551 SpatialDomains::PointGeomUniquePtr v4(
2552 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2553
2554 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2555 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2556 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2557 std::array<SpatialDomains::PointGeom *, 5> v = {
2558 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2559 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2560 CreatePyr(v, segVec, triVec, quadVec);
2561
2562 unsigned int numQuadPoints = 5;
2563 unsigned int numModes = 4;
2564
2565 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2566 Nektar::LibUtilities::eGaussLobattoLegendre;
2567 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
2568 PointsTypeDir1);
2569 Nektar::LibUtilities::BasisType basisTypeDir1 =
2570 Nektar::LibUtilities::eModified_A;
2571 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2572 PointsKeyDir1);
2573
2574 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2575 Nektar::LibUtilities::eGaussLobattoLegendre;
2576 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
2577 PointsTypeDir2);
2578 Nektar::LibUtilities::BasisType basisTypeDir2 =
2579 Nektar::LibUtilities::eModified_A;
2580 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2581 PointsKeyDir2);
2582
2583 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2584 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2585 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
2586 PointsTypeDir3);
2587 Nektar::LibUtilities::BasisType basisTypeDir3 =
2588 Nektar::LibUtilities::eModifiedPyr_C;
2589 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2590 PointsKeyDir3);
2591
2592 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2593 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2594 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2595
2596 int nelmts = NELMTS;
2597
2598 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2599 for (int i = 0; i < nelmts; ++i)
2600 {
2601 CollExp.push_back(Exp);
2602 }
2603
2604 LibUtilities::SessionReaderSharedPtr dummySession;
2605 Collections::CollectionOptimisation colOpt(dummySession, 2,
2606 Collections::eMatrixFree);
2607 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2608 Collections::Collection c(CollExp, impTypes);
2609 c.Initialise(Collections::eIProductWRTDerivBase);
2610
2611 const int nq = Exp->GetTotPoints();
2612 const int nm = Exp->GetNcoeffs();
2613 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2614 Array<OneD, NekDouble> phys2(nelmts * nq);
2615 Array<OneD, NekDouble> phys3(nelmts * nq);
2616 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2617 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2618
2619 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2620 Exp->GetCoords(xc, yc, zc);
2621
2622 for (int i = 0; i < nq; ++i)
2623 {
2624 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2625 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2626 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2627 }
2628 for (int i = 1; i < nelmts; ++i)
2629 {
2630 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2631 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2632 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2633 }
2634
2635 // Standard routines
2636 for (int i = 0; i < nelmts; ++i)
2637 {
2638 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2639 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2640 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2641 tmp = coeffs1 + i * nm, 1);
2642 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2643 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2644 tmp = coeffs1 + i * nm, 1);
2645 }
2646
2647 LibUtilities::Timer timer;
2648 timer.Start();
2649 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2650 coeffs2);
2651 timer.Stop();
2652 timer.AccumulateRegion("Pyr IProdWRTDB MatrixFree");
2653
2654 double epsilon = 1.0e-8;
2655 for (int i = 0; i < coeffs1.size(); ++i)
2656 {
2657 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2658 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2659 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2660 }
2661}
2662BOOST_AUTO_TEST_CASE(
2663 TestPyrIProductWRTDerivBase_MatrixFree_Deformed_MultiElmt_OverInt)
2664{
2665 SpatialDomains::PointGeomUniquePtr v0(
2666 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2667 SpatialDomains::PointGeomUniquePtr v1(
2668 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2669 SpatialDomains::PointGeomUniquePtr v2(
2670 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2671 SpatialDomains::PointGeomUniquePtr v3(
2672 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2673 SpatialDomains::PointGeomUniquePtr v4(
2674 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2675
2676 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2677 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2678 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2679 std::array<SpatialDomains::PointGeom *, 5> v = {
2680 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2681 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2682 CreatePyr(v, segVec, triVec, quadVec);
2683
2684 unsigned int numQuadPoints = 8;
2685 unsigned int numModes = 4;
2686
2687 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2688 Nektar::LibUtilities::eGaussLobattoLegendre;
2689 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
2690 PointsTypeDir1);
2691 Nektar::LibUtilities::BasisType basisTypeDir1 =
2692 Nektar::LibUtilities::eModified_A;
2693 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2694 PointsKeyDir1);
2695
2696 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2697 Nektar::LibUtilities::eGaussLobattoLegendre;
2698 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
2699 PointsTypeDir2);
2700 Nektar::LibUtilities::BasisType basisTypeDir2 =
2701 Nektar::LibUtilities::eModified_A;
2702 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2703 PointsKeyDir2);
2704
2705 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2706 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2707 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
2708 PointsTypeDir3);
2709 Nektar::LibUtilities::BasisType basisTypeDir3 =
2710 Nektar::LibUtilities::eModifiedPyr_C;
2711 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2712 PointsKeyDir3);
2713
2714 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2715 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2716 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2717
2718 int nelmts = 5;
2719
2720 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2721 for (int i = 0; i < nelmts; ++i)
2722 {
2723 CollExp.push_back(Exp);
2724 }
2725
2726 LibUtilities::SessionReaderSharedPtr dummySession;
2727 Collections::CollectionOptimisation colOpt(dummySession, 2,
2728 Collections::eMatrixFree);
2729 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2730 Collections::Collection c(CollExp, impTypes);
2731 c.Initialise(Collections::eIProductWRTDerivBase);
2732
2733 const int nq = Exp->GetTotPoints();
2734 const int nm = Exp->GetNcoeffs();
2735 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2736 Array<OneD, NekDouble> phys2(nelmts * nq);
2737 Array<OneD, NekDouble> phys3(nelmts * nq);
2738 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2739 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2740
2741 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2742 Exp->GetCoords(xc, yc, zc);
2743
2744 for (int i = 0; i < nq; ++i)
2745 {
2746 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2747 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2748 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2749 }
2750 for (int i = 1; i < nelmts; ++i)
2751 {
2752 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2753 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2754 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2755 }
2756
2757 // Standard routines
2758 for (int i = 0; i < nelmts; ++i)
2759 {
2760 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2761 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2762 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2763 tmp = coeffs1 + i * nm, 1);
2764 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2765 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2766 tmp = coeffs1 + i * nm, 1);
2767 }
2768
2769 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2770 coeffs2);
2771
2772 double epsilon = 1.0e-8;
2773 for (int i = 0; i < coeffs1.size(); ++i)
2774 {
2775 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2776 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2777 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2778 }
2779}
2780
2781BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_MatrixFree_UniformP_MultiElmt)
2782{
2783 SpatialDomains::PointGeomUniquePtr v0(
2784 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2785 SpatialDomains::PointGeomUniquePtr v1(
2786 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2787 SpatialDomains::PointGeomUniquePtr v2(
2788 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2789 SpatialDomains::PointGeomUniquePtr v3(
2790 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2791 SpatialDomains::PointGeomUniquePtr v4(
2792 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2793
2794 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2795 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2796 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2797 std::array<SpatialDomains::PointGeom *, 5> v = {
2798 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2799 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2800 CreatePyr(v, segVec, triVec, quadVec);
2801
2802 unsigned int numQuadPoints = 5;
2803 unsigned int numModes = 4;
2804
2805 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2806 Nektar::LibUtilities::eGaussLobattoLegendre;
2807 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
2808 PointsTypeDir1);
2809 Nektar::LibUtilities::BasisType basisTypeDir1 =
2810 Nektar::LibUtilities::eModified_A;
2811 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2812 PointsKeyDir1);
2813
2814 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2815 Nektar::LibUtilities::eGaussLobattoLegendre;
2816 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
2817 PointsTypeDir2);
2818 Nektar::LibUtilities::BasisType basisTypeDir2 =
2819 Nektar::LibUtilities::eModified_A;
2820 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2821 PointsKeyDir2);
2822
2823 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2824 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2825 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
2826 PointsTypeDir3);
2827 Nektar::LibUtilities::BasisType basisTypeDir3 =
2828 Nektar::LibUtilities::eModifiedPyr_C;
2829 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2830 PointsKeyDir3);
2831
2832 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2833 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2834 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2835
2836 unsigned int nelmts = NELMTS;
2837
2838 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2839 for (unsigned int i = 0; i < nelmts; ++i)
2840 {
2841 CollExp.push_back(Exp);
2842 }
2843
2844 LibUtilities::SessionReaderSharedPtr dummySession;
2845 Collections::CollectionOptimisation colOpt(dummySession, 2,
2846 Collections::eIterPerExp);
2847 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2848
2849 // ... only one op at the time ...
2850 impTypes[Collections::eBwdTrans] = Collections::eMatrixFree;
2851 Collections::Collection c(CollExp, impTypes);
2852 c.Initialise(Collections::eBwdTrans);
2853
2854 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
2855 for (int i = 0; i < coeffs.size(); ++i)
2856 {
2857 coeffs[i] = i + 1;
2858 }
2859 Array<OneD, NekDouble> physRef(nelmts * Exp->GetTotPoints(), 0.0);
2860 Array<OneD, NekDouble> phys(nelmts * Exp->GetTotPoints(), 0.0);
2861 Array<OneD, NekDouble> tmp;
2862
2863 for (unsigned int i = 0; i < nelmts; ++i)
2864 {
2865 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
2866 tmp = physRef + i * Exp->GetTotPoints());
2867 }
2868 LibUtilities::Timer timer;
2869 timer.Start();
2870 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys);
2871 timer.Stop();
2872 timer.AccumulateRegion("Pyr BwdTrans MatrixFree");
2873
2874 double epsilon = 1.0e-8;
2875 for (unsigned int i = 0; i < physRef.size(); ++i)
2876 {
2877 BOOST_CHECK_CLOSE(physRef[i], phys[i], epsilon);
2878 }
2879}
2880
2881BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_MatrixFree_UniformP_OverInt_MultiElmt)
2882{
2883 SpatialDomains::PointGeomUniquePtr v0(
2884 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2885 SpatialDomains::PointGeomUniquePtr v1(
2886 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2887 SpatialDomains::PointGeomUniquePtr v2(
2888 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2889 SpatialDomains::PointGeomUniquePtr v3(
2890 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2891 SpatialDomains::PointGeomUniquePtr v4(
2892 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2893
2894 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2895 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2896 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2897 std::array<SpatialDomains::PointGeom *, 5> v = {
2898 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2899 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2900 CreatePyr(v, segVec, triVec, quadVec);
2901
2902 unsigned int numQuadPoints = 8;
2903 unsigned int numModes = 4;
2904
2905 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2906 Nektar::LibUtilities::eGaussLobattoLegendre;
2907 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
2908 PointsTypeDir1);
2909 Nektar::LibUtilities::BasisType basisTypeDir1 =
2910 Nektar::LibUtilities::eModified_A;
2911 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2912 PointsKeyDir1);
2913
2914 Nektar::LibUtilities::PointsType PointsTypeDir2 =
2915 Nektar::LibUtilities::eGaussLobattoLegendre;
2916 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
2917 PointsTypeDir2);
2918 Nektar::LibUtilities::BasisType basisTypeDir2 =
2919 Nektar::LibUtilities::eModified_A;
2920 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2921 PointsKeyDir2);
2922
2923 Nektar::LibUtilities::PointsType PointsTypeDir3 =
2924 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2925 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
2926 PointsTypeDir3);
2927 Nektar::LibUtilities::BasisType basisTypeDir3 =
2928 Nektar::LibUtilities::eModifiedPyr_C;
2929 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2930 PointsKeyDir3);
2931
2932 Nektar::LocalRegions::PyrExpSharedPtr Exp =
2933 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
2934 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
2935
2936 unsigned int nelmts = 10;
2937
2938 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
2939 for (unsigned int i = 0; i < nelmts; ++i)
2940 {
2941 CollExp.push_back(Exp);
2942 }
2943
2944 LibUtilities::SessionReaderSharedPtr dummySession;
2945 Collections::CollectionOptimisation colOpt(dummySession, 2,
2946 Collections::eIterPerExp);
2947 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2948
2949 // ... only one op at the time ...
2950 impTypes[Collections::eBwdTrans] = Collections::eMatrixFree;
2951 Collections::Collection c(CollExp, impTypes);
2952 c.Initialise(Collections::eBwdTrans);
2953
2954 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0);
2955 for (int i = 0; i < coeffs.size(); ++i)
2956 {
2957 coeffs[i] = i + 1;
2958 }
2959 Array<OneD, NekDouble> physRef(nelmts * Exp->GetTotPoints(), 0.0);
2960 Array<OneD, NekDouble> phys(nelmts * Exp->GetTotPoints(), 0.0);
2961 Array<OneD, NekDouble> tmp;
2962
2963 for (unsigned int i = 0; i < nelmts; ++i)
2964 {
2965 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
2966 tmp = physRef + i * Exp->GetTotPoints());
2967 }
2968 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys);
2969
2970 double epsilon = 1.0e-8;
2971 for (unsigned int i = 0; i < physRef.size(); ++i)
2972 {
2973 BOOST_CHECK_CLOSE(physRef[i], phys[i], epsilon);
2974 }
2975}
2976
2977BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_StdMat_VariableP_MultiElmt)
2978{
2979 SpatialDomains::PointGeomUniquePtr v0(
2980 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
2981 SpatialDomains::PointGeomUniquePtr v1(
2982 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2983 SpatialDomains::PointGeomUniquePtr v2(
2984 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
2985 SpatialDomains::PointGeomUniquePtr v3(
2986 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
2987 SpatialDomains::PointGeomUniquePtr v4(
2988 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
2989
2990 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
2991 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
2992 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
2993 std::array<SpatialDomains::PointGeom *, 5> v = {
2994 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
2995 SpatialDomains::PyrGeomUniquePtr pyrGeom =
2996 CreatePyr(v, segVec, triVec, quadVec);
2997
2998 Nektar::LibUtilities::PointsType PointsTypeDir1 =
2999 Nektar::LibUtilities::eGaussLobattoLegendre;
3000 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
3001 Nektar::LibUtilities::BasisType basisTypeDir1 =
3002 Nektar::LibUtilities::eModified_A;
3003 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3004 PointsKeyDir1);
3005
3006 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3007 Nektar::LibUtilities::eGaussLobattoLegendre;
3008 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
3009 Nektar::LibUtilities::BasisType basisTypeDir2 =
3010 Nektar::LibUtilities::eModified_A;
3011 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
3012 PointsKeyDir2);
3013
3014 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3015 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3016 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
3017 Nektar::LibUtilities::BasisType basisTypeDir3 =
3018 Nektar::LibUtilities::eModifiedPyr_C;
3019 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
3020 PointsKeyDir3);
3021
3022 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3023 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3024 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3025
3026 int nelmts = 10;
3027
3028 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3029 for (int i = 0; i < nelmts; ++i)
3030 {
3031 CollExp.push_back(Exp);
3032 }
3033
3034 LibUtilities::SessionReaderSharedPtr dummySession;
3035 Collections::CollectionOptimisation colOpt(dummySession, 3,
3036 Collections::eStdMat);
3037 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3038 Collections::Collection c(CollExp, impTypes);
3039 c.Initialise(Collections::eIProductWRTDerivBase);
3040
3041 const int nq = Exp->GetTotPoints();
3042 const int nm = Exp->GetNcoeffs();
3043 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
3044 Array<OneD, NekDouble> phys2(nelmts * nq);
3045 Array<OneD, NekDouble> phys3(nelmts * nq);
3046 Array<OneD, NekDouble> coeffs1(nelmts * nm);
3047 Array<OneD, NekDouble> coeffs2(nelmts * nm);
3048
3049 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3050
3051 Exp->GetCoords(xc, yc, zc);
3052
3053 for (int i = 0; i < nq; ++i)
3054 {
3055 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3056 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3057 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3058 }
3059 for (int i = 1; i < nelmts; ++i)
3060 {
3061 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
3062 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
3063 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
3064 }
3065
3066 // Standard routines
3067 for (int i = 0; i < nelmts; ++i)
3068 {
3069 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
3070 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
3071 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3072 tmp = coeffs1 + i * nm, 1);
3073 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
3074 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3075 tmp = coeffs1 + i * nm, 1);
3076 }
3077
3078 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3079 coeffs2);
3080
3081 double epsilon = 1.0e-8;
3082 for (int i = 0; i < coeffs1.size(); ++i)
3083 {
3084 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3085 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3086 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3087 }
3088}
3089
3090BOOST_AUTO_TEST_CASE(TestPyrIProductWRTDerivBase_SumFac_VariableP_MultiElmt)
3091{
3092 SpatialDomains::PointGeomUniquePtr v0(
3093 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3094 SpatialDomains::PointGeomUniquePtr v1(
3095 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3096 SpatialDomains::PointGeomUniquePtr v2(
3097 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3098 SpatialDomains::PointGeomUniquePtr v3(
3099 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3100 SpatialDomains::PointGeomUniquePtr v4(
3101 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3102
3103 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3104 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3105 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3106 std::array<SpatialDomains::PointGeom *, 5> v = {
3107 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3108 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3109 CreatePyr(v, segVec, triVec, quadVec);
3110
3111 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3112 Nektar::LibUtilities::eGaussLobattoLegendre;
3113 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
3114 Nektar::LibUtilities::BasisType basisTypeDir1 =
3115 Nektar::LibUtilities::eModified_A;
3116 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3117 PointsKeyDir1);
3118
3119 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3120 Nektar::LibUtilities::eGaussLobattoLegendre;
3121 const Nektar::LibUtilities::PointsKey PointsKeyDir2(7, PointsTypeDir2);
3122 Nektar::LibUtilities::BasisType basisTypeDir2 =
3123 Nektar::LibUtilities::eModified_A;
3124 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
3125 PointsKeyDir2);
3126
3127 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3128 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3129 const Nektar::LibUtilities::PointsKey PointsKeyDir3(8, PointsTypeDir3);
3130 Nektar::LibUtilities::BasisType basisTypeDir3 =
3131 Nektar::LibUtilities::eModifiedPyr_C;
3132 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
3133 PointsKeyDir3);
3134
3135 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3136 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3137 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3138
3139 int nelmts = 10;
3140
3141 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3142 for (int i = 0; i < nelmts; ++i)
3143 {
3144 CollExp.push_back(Exp);
3145 }
3146
3147 LibUtilities::SessionReaderSharedPtr dummySession;
3148 Collections::CollectionOptimisation colOpt(dummySession, 3,
3149 Collections::eSumFac);
3150 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3151 Collections::Collection c(CollExp, impTypes);
3152 c.Initialise(Collections::eIProductWRTDerivBase);
3153
3154 const int nq = Exp->GetTotPoints();
3155 const int nm = Exp->GetNcoeffs();
3156 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
3157 Array<OneD, NekDouble> phys2(nelmts * nq);
3158 Array<OneD, NekDouble> phys3(nelmts * nq);
3159 Array<OneD, NekDouble> coeffs1(nelmts * nm);
3160 Array<OneD, NekDouble> coeffs2(nelmts * nm);
3161
3162 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3163
3164 Exp->GetCoords(xc, yc, zc);
3165
3166 for (int i = 0; i < nq; ++i)
3167 {
3168 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3169 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3170 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3171 }
3172 for (int i = 1; i < nelmts; ++i)
3173 {
3174 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
3175 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
3176 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
3177 }
3178
3179 // Standard routines
3180 for (int i = 0; i < nelmts; ++i)
3181 {
3182 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
3183 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
3184 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3185 tmp = coeffs1 + i * nm, 1);
3186 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
3187 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
3188 tmp = coeffs1 + i * nm, 1);
3189 }
3190
3191 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3192 coeffs2);
3193
3194 double epsilon = 1.0e-8;
3195 for (int i = 0; i < coeffs1.size(); ++i)
3196 {
3197 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
3198 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
3199 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
3200 }
3201}
3202
3203BOOST_AUTO_TEST_CASE(TestPyrHelmholtz_IterPerExp_UniformP_ConstVarDiff)
3204{
3205 SpatialDomains::PointGeomUniquePtr v0(
3206 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3207 SpatialDomains::PointGeomUniquePtr v1(
3208 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3209 SpatialDomains::PointGeomUniquePtr v2(
3210 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3211 SpatialDomains::PointGeomUniquePtr v3(
3212 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3213 SpatialDomains::PointGeomUniquePtr v4(
3214 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3215
3216 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3217 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3218 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3219 std::array<SpatialDomains::PointGeom *, 5> v = {
3220 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3221 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3222 CreatePyr(v, segVec, triVec, quadVec);
3223
3224 unsigned int numQuadPoints = 5;
3225 unsigned int numModes = 4;
3226
3227 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3228 Nektar::LibUtilities::eGaussLobattoLegendre;
3229 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3230 PointsTypeDir1);
3231 Nektar::LibUtilities::BasisType basisTypeDir1 =
3232 Nektar::LibUtilities::eModified_A;
3233 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3234 PointsKeyDir1);
3235
3236 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3237 Nektar::LibUtilities::eGaussLobattoLegendre;
3238 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3239 PointsTypeDir2);
3240 Nektar::LibUtilities::BasisType basisTypeDir2 =
3241 Nektar::LibUtilities::eModified_A;
3242 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3243 PointsKeyDir2);
3244
3245 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3246 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3247 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3248 PointsTypeDir3);
3249 Nektar::LibUtilities::BasisType basisTypeDir3 =
3250 Nektar::LibUtilities::eModifiedPyr_C;
3251 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3252 PointsKeyDir3);
3253
3254 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3255 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3256 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3257
3258 int nelmts = 10;
3259
3260 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3261 for (int i = 0; i < nelmts; ++i)
3262 {
3263 CollExp.push_back(Exp);
3264 }
3265
3266 LibUtilities::SessionReaderSharedPtr dummySession;
3267 Collections::CollectionOptimisation colOpt(dummySession, 2,
3268 Collections::eIterPerExp);
3269 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3270 Collections::Collection c(CollExp, impTypes);
3271 StdRegions::ConstFactorMap factors;
3272 factors[StdRegions::eFactorLambda] = 1.5;
3273 factors[StdRegions::eFactorCoeffD00] = 1.25;
3274 factors[StdRegions::eFactorCoeffD01] = 0.25;
3275 factors[StdRegions::eFactorCoeffD11] = 1.25;
3276 factors[StdRegions::eFactorCoeffD02] = 0.25;
3277 factors[StdRegions::eFactorCoeffD12] = 0.25;
3278 factors[StdRegions::eFactorCoeffD22] = 1.25;
3279
3280 c.Initialise(Collections::eHelmholtz, factors);
3281
3282 const int nm = Exp->GetNcoeffs();
3283 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3284 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3285 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3286
3287 for (int i = 0; i < coeffsIn.size(); ++i)
3288 {
3289 coeffsIn[i] = i + 1.0;
3290 }
3291
3292 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3293 *Exp, factors);
3294
3295 for (int i = 0; i < nelmts; ++i)
3296 {
3297 // Standard routines
3298 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3299 }
3300
3301 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3302
3303 double epsilon = 1.0e-8;
3304 for (int i = 0; i < coeffsRef.size(); ++i)
3305 {
3306 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3307 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3308 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3309 }
3310}
3311
3312BOOST_AUTO_TEST_CASE(TestPyrHelmholtz_MatrixFree_UniformP)
3313{
3314 SpatialDomains::PointGeomUniquePtr v0(
3315 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3316 SpatialDomains::PointGeomUniquePtr v1(
3317 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3318 SpatialDomains::PointGeomUniquePtr v2(
3319 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3320 SpatialDomains::PointGeomUniquePtr v3(
3321 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3322 SpatialDomains::PointGeomUniquePtr v4(
3323 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3324
3325 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3326 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3327 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3328 std::array<SpatialDomains::PointGeom *, 5> v = {
3329 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3330 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3331 CreatePyr(v, segVec, triVec, quadVec);
3332
3333 unsigned int numQuadPoints = 5;
3334 unsigned int numModes = 4;
3335
3336 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3337 Nektar::LibUtilities::eGaussLobattoLegendre;
3338 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3339 PointsTypeDir1);
3340 Nektar::LibUtilities::BasisType basisTypeDir1 =
3341 Nektar::LibUtilities::eModified_A;
3342 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3343 PointsKeyDir1);
3344
3345 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3346 Nektar::LibUtilities::eGaussLobattoLegendre;
3347 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3348 PointsTypeDir2);
3349 Nektar::LibUtilities::BasisType basisTypeDir2 =
3350 Nektar::LibUtilities::eModified_A;
3351 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3352 PointsKeyDir2);
3353
3354 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3355 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3356 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3357 PointsTypeDir3);
3358 Nektar::LibUtilities::BasisType basisTypeDir3 =
3359 Nektar::LibUtilities::eModifiedPyr_C;
3360 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3361 PointsKeyDir3);
3362
3363 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3364 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3365 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3366
3367 int nelmts = 10;
3368
3369 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3370 for (int i = 0; i < nelmts; ++i)
3371 {
3372 CollExp.push_back(Exp);
3373 }
3374
3375 LibUtilities::SessionReaderSharedPtr dummySession;
3376 Collections::CollectionOptimisation colOpt(dummySession, 2,
3377 Collections::eMatrixFree);
3378 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3379 Collections::Collection c(CollExp, impTypes);
3380 StdRegions::ConstFactorMap factors;
3381 factors[StdRegions::eFactorLambda] = 1.5;
3382
3383 c.Initialise(Collections::eHelmholtz, factors);
3384
3385 const int nm = Exp->GetNcoeffs();
3386 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3387 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3388 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3389
3390 for (int i = 0; i < coeffsIn.size(); ++i)
3391 {
3392 coeffsIn[i] = i + 1.0;
3393 }
3394
3395 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3396 *Exp, factors);
3397
3398 for (int i = 0; i < nelmts; ++i)
3399 {
3400 // Standard routines
3401 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3402 }
3403
3404 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3405
3406 double epsilon = 1.0e-8;
3407 for (int i = 0; i < coeffsRef.size(); ++i)
3408 {
3409 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3410 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3411 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3412 }
3413}
3414
3415BOOST_AUTO_TEST_CASE(TestPyrHelmholtz_MatrixFree_Deformed_overInt)
3416{
3417 SpatialDomains::PointGeomUniquePtr v0(
3418 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
3419 SpatialDomains::PointGeomUniquePtr v1(
3420 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3421 SpatialDomains::PointGeomUniquePtr v2(
3422 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3423 SpatialDomains::PointGeomUniquePtr v3(
3424 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3425 SpatialDomains::PointGeomUniquePtr v4(
3426 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3427
3428 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3429 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3430 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3431 std::array<SpatialDomains::PointGeom *, 5> v = {
3432 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3433 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3434 CreatePyr(v, segVec, triVec, quadVec);
3435
3436 unsigned int numQuadPoints = 8;
3437 unsigned int numModes = 4;
3438
3439 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3440 Nektar::LibUtilities::eGaussLobattoLegendre;
3441 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3442 PointsTypeDir1);
3443 Nektar::LibUtilities::BasisType basisTypeDir1 =
3444 Nektar::LibUtilities::eModified_A;
3445 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3446 PointsKeyDir1);
3447
3448 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3449 Nektar::LibUtilities::eGaussLobattoLegendre;
3450 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3451 PointsTypeDir2);
3452 Nektar::LibUtilities::BasisType basisTypeDir2 =
3453 Nektar::LibUtilities::eModified_A;
3454 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3455 PointsKeyDir2);
3456
3457 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3458 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3459 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3460 PointsTypeDir3);
3461 Nektar::LibUtilities::BasisType basisTypeDir3 =
3462 Nektar::LibUtilities::eModifiedPyr_C;
3463 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3464 PointsKeyDir3);
3465
3466 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3467 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3468 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3469
3470 int nelmts = 10;
3471
3472 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3473 for (int i = 0; i < nelmts; ++i)
3474 {
3475 CollExp.push_back(Exp);
3476 }
3477
3478 LibUtilities::SessionReaderSharedPtr dummySession;
3479 Collections::CollectionOptimisation colOpt(dummySession, 2,
3480 Collections::eMatrixFree);
3481 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3482 Collections::Collection c(CollExp, impTypes);
3483 StdRegions::ConstFactorMap factors;
3484 factors[StdRegions::eFactorLambda] = 1.5;
3485
3486 c.Initialise(Collections::eHelmholtz, factors);
3487
3488 const int nm = Exp->GetNcoeffs();
3489 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3490 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3491 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3492
3493 for (int i = 0; i < coeffsIn.size(); ++i)
3494 {
3495 coeffsIn[i] = i + 1.0;
3496 }
3497
3498 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3499 *Exp, factors);
3500
3501 for (int i = 0; i < nelmts; ++i)
3502 {
3503 // Standard routines
3504 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3505 }
3506
3507 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3508
3509 double epsilon = 1.0e-8;
3510 for (int i = 0; i < coeffsRef.size(); ++i)
3511 {
3512 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3513 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3514 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3515 }
3516}
3517
3518BOOST_AUTO_TEST_CASE(TestPyrHelmholtz_MatrixFree_UniformP_ConstVarDiff)
3519{
3520 SpatialDomains::PointGeomUniquePtr v0(
3521 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3522 SpatialDomains::PointGeomUniquePtr v1(
3523 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3524 SpatialDomains::PointGeomUniquePtr v2(
3525 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3526 SpatialDomains::PointGeomUniquePtr v3(
3527 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3528 SpatialDomains::PointGeomUniquePtr v4(
3529 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3530
3531 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3532 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3533 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3534 std::array<SpatialDomains::PointGeom *, 5> v = {
3535 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3536 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3537 CreatePyr(v, segVec, triVec, quadVec);
3538
3539 unsigned int numQuadPoints = 5;
3540 unsigned int numModes = 4;
3541
3542 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3543 Nektar::LibUtilities::eGaussLobattoLegendre;
3544 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3545 PointsTypeDir1);
3546 Nektar::LibUtilities::BasisType basisTypeDir1 =
3547 Nektar::LibUtilities::eModified_A;
3548 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3549 PointsKeyDir1);
3550
3551 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3552 Nektar::LibUtilities::eGaussLobattoLegendre;
3553 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3554 PointsTypeDir2);
3555 Nektar::LibUtilities::BasisType basisTypeDir2 =
3556 Nektar::LibUtilities::eModified_A;
3557 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3558 PointsKeyDir2);
3559
3560 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3561 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3562 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3563 PointsTypeDir3);
3564 Nektar::LibUtilities::BasisType basisTypeDir3 =
3565 Nektar::LibUtilities::eModifiedPyr_C;
3566 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3567 PointsKeyDir3);
3568
3569 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3570 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3571 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3572
3573 int nelmts = 10;
3574
3575 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3576 for (int i = 0; i < nelmts; ++i)
3577 {
3578 CollExp.push_back(Exp);
3579 }
3580
3581 LibUtilities::SessionReaderSharedPtr dummySession;
3582 Collections::CollectionOptimisation colOpt(dummySession, 3,
3583 Collections::eMatrixFree);
3584 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3585 Collections::Collection c(CollExp, impTypes);
3586 StdRegions::ConstFactorMap factors;
3587 factors[StdRegions::eFactorLambda] = 1.5;
3588 factors[StdRegions::eFactorCoeffD00] = 1.25;
3589 factors[StdRegions::eFactorCoeffD01] = 0.25;
3590 factors[StdRegions::eFactorCoeffD11] = 1.25;
3591 factors[StdRegions::eFactorCoeffD02] = 0.25;
3592 factors[StdRegions::eFactorCoeffD12] = 0.25;
3593 factors[StdRegions::eFactorCoeffD22] = 1.25;
3594
3595 c.Initialise(Collections::eHelmholtz, factors);
3596
3597 const int nm = Exp->GetNcoeffs();
3598 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3599 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3600 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3601
3602 for (int i = 0; i < coeffsIn.size(); ++i)
3603 {
3604 coeffsIn[i] = i + 1.0;
3605 }
3606
3607 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3608 *Exp, factors);
3609
3610 for (int i = 0; i < nelmts; ++i)
3611 {
3612 // Standard routines
3613 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3614 }
3615
3616 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3617
3618 double epsilon = 1.0e-8;
3619 for (int i = 0; i < coeffsRef.size(); ++i)
3620 {
3621 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3622 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3623 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3624 }
3625}
3626
3627BOOST_AUTO_TEST_CASE(TestPyrPhysInterp1DScaled_NoCollection_UniformP_MultiElmt)
3628{
3629 SpatialDomains::PointGeomUniquePtr v0(
3630 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3631 SpatialDomains::PointGeomUniquePtr v1(
3632 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3633 SpatialDomains::PointGeomUniquePtr v2(
3634 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3635 SpatialDomains::PointGeomUniquePtr v3(
3636 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3637 SpatialDomains::PointGeomUniquePtr v4(
3638 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3639
3640 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3641 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3642 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3643 std::array<SpatialDomains::PointGeom *, 5> v = {
3644 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3645 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3646 CreatePyr(v, segVec, triVec, quadVec);
3647
3648 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3649 Nektar::LibUtilities::eGaussLobattoLegendre;
3650 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
3651 Nektar::LibUtilities::BasisType basisTypeDir1 =
3652 Nektar::LibUtilities::eModified_A;
3653 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3654 PointsKeyDir1);
3655
3656 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3657 Nektar::LibUtilities::eGaussLobattoLegendre;
3658 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
3659 Nektar::LibUtilities::BasisType basisTypeDir2 =
3660 Nektar::LibUtilities::eModified_A;
3661 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
3662 PointsKeyDir2);
3663
3664 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3665 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3666 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
3667 Nektar::LibUtilities::BasisType basisTypeDir3 =
3668 Nektar::LibUtilities::eModifiedPyr_C;
3669 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
3670 PointsKeyDir3);
3671
3672 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3673 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3674 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3675
3676 int nelmts = 1;
3677
3678 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3679 for (int i = 0; i < nelmts; ++i)
3680 {
3681 CollExp.push_back(Exp);
3682 }
3683
3684 LibUtilities::SessionReaderSharedPtr dummySession;
3685 Collections::CollectionOptimisation colOpt(dummySession, 3,
3686 Collections::eNoCollection);
3687 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3688 Collections::Collection c(CollExp, impTypes);
3689 StdRegions::ConstFactorMap factors;
3690 factors[StdRegions::eFactorConst] = 1.5;
3691 c.Initialise(Collections::ePhysInterp1DScaled, factors);
3692
3693 const int nq = Exp->GetTotPoints();
3694 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3695 Array<OneD, NekDouble> phys(nq);
3696
3697 Exp->GetCoords(xc, yc, zc);
3698
3699 for (int i = 0; i < nq; ++i)
3700 {
3701 phys[i] = pow(xc[i], 3) + pow(yc[i], 3) + pow(zc[i], 3);
3702 }
3703
3704 const int nq1 = c.GetOutputSize(Collections::ePhysInterp1DScaled);
3705 Array<OneD, NekDouble> xc1(nq1);
3706 Array<OneD, NekDouble> yc1(nq1);
3707 Array<OneD, NekDouble> zc1(nq1);
3708 Array<OneD, NekDouble> phys1(nq1);
3709
3710 c.ApplyOperator(Collections::ePhysInterp1DScaled, xc, xc1);
3711 c.ApplyOperator(Collections::ePhysInterp1DScaled, yc, yc1);
3712 c.ApplyOperator(Collections::ePhysInterp1DScaled, zc, zc1);
3713 c.ApplyOperator(Collections::ePhysInterp1DScaled, phys, phys1);
3714
3715 double epsilon = 2.0e-8;
3716 for (int i = 0; i < nq1; ++i)
3717 {
3718 NekDouble exact = pow(xc1[i], 3) + pow(yc1[i], 3) + pow(zc1[i], 3);
3719 phys1[i] = (fabs(phys1[i]) < 1e-14) ? 0.0 : phys1[i];
3720 exact = (fabs(exact) < 1e-14) ? 0.0 : exact;
3721 BOOST_CHECK_CLOSE(phys1[i], exact, epsilon);
3722 }
3723}
3724
3725BOOST_AUTO_TEST_CASE(TestPyrPhysInterp1DScaled_MatrixFree_UniformP_MultiElmt)
3726{
3727 SpatialDomains::PointGeomUniquePtr v0(
3728 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3729 SpatialDomains::PointGeomUniquePtr v1(
3730 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3731 SpatialDomains::PointGeomUniquePtr v2(
3732 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3733 SpatialDomains::PointGeomUniquePtr v3(
3734 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3735 SpatialDomains::PointGeomUniquePtr v4(
3736 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3737
3738 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3739 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3740 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3741 std::array<SpatialDomains::PointGeom *, 5> v = {
3742 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3743 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3744 CreatePyr(v, segVec, triVec, quadVec);
3745
3746 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3747 Nektar::LibUtilities::eGaussLobattoLegendre;
3748 const Nektar::LibUtilities::PointsKey PointsKeyDir1(5, PointsTypeDir1);
3749 Nektar::LibUtilities::BasisType basisTypeDir1 =
3750 Nektar::LibUtilities::eModified_A;
3751 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
3752 PointsKeyDir1);
3753
3754 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3755 Nektar::LibUtilities::eGaussLobattoLegendre;
3756 const Nektar::LibUtilities::PointsKey PointsKeyDir2(5, PointsTypeDir2);
3757 Nektar::LibUtilities::BasisType basisTypeDir2 =
3758 Nektar::LibUtilities::eModified_A;
3759 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
3760 PointsKeyDir2);
3761
3762 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3763 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3764 const Nektar::LibUtilities::PointsKey PointsKeyDir3(4, PointsTypeDir3);
3765 Nektar::LibUtilities::BasisType basisTypeDir3 =
3766 Nektar::LibUtilities::eModifiedPyr_C;
3767 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
3768 PointsKeyDir3);
3769
3770 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3771 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3772 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3773
3774 int nelmts = 1;
3775
3776 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3777 for (int i = 0; i < nelmts; ++i)
3778 {
3779 CollExp.push_back(Exp);
3780 }
3781
3782 LibUtilities::SessionReaderSharedPtr dummySession;
3783 Collections::CollectionOptimisation colOpt(dummySession, 3,
3784 Collections::eMatrixFree);
3785 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3786 Collections::Collection c(CollExp, impTypes);
3787 StdRegions::ConstFactorMap factors;
3788 factors[StdRegions::eFactorConst] = 1.5;
3789 c.Initialise(Collections::ePhysInterp1DScaled, factors);
3790
3791 const int nq = Exp->GetTotPoints();
3792 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3793 Array<OneD, NekDouble> phys(nq);
3794
3795 Exp->GetCoords(xc, yc, zc);
3796
3797 for (int i = 0; i < nq; ++i)
3798 {
3799 phys[i] = pow(xc[i], 3) + pow(yc[i], 3) + pow(zc[i], 3);
3800 }
3801
3802 const int nq1 = c.GetOutputSize(Collections::ePhysInterp1DScaled);
3803 Array<OneD, NekDouble> xc1(nq1);
3804 Array<OneD, NekDouble> yc1(nq1);
3805 Array<OneD, NekDouble> zc1(nq1);
3806 Array<OneD, NekDouble> phys1(nq1);
3807
3808 c.ApplyOperator(Collections::ePhysInterp1DScaled, xc, xc1);
3809 c.ApplyOperator(Collections::ePhysInterp1DScaled, yc, yc1);
3810 c.ApplyOperator(Collections::ePhysInterp1DScaled, zc, zc1);
3811 c.ApplyOperator(Collections::ePhysInterp1DScaled, phys, phys1);
3812
3813 double epsilon = 2.0e-8;
3814 for (int i = 0; i < nq1; ++i)
3815 {
3816 NekDouble exact = pow(xc1[i], 3) + pow(yc1[i], 3) + pow(zc1[i], 3);
3817 phys1[i] = (fabs(phys1[i]) < 1e-14) ? 0.0 : phys1[i];
3818 exact = (fabs(exact) < 1e-14) ? 0.0 : exact;
3819 BOOST_CHECK_CLOSE(phys1[i], exact, epsilon);
3820 }
3821}
3822
3823BOOST_AUTO_TEST_CASE(
3824 TestPyrLinearAdvectionDiffusionReaction_IterPerExp_UniformP_ConstVarDiff)
3825{
3826 SpatialDomains::PointGeomUniquePtr v0(
3827 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3828 SpatialDomains::PointGeomUniquePtr v1(
3829 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3830 SpatialDomains::PointGeomUniquePtr v2(
3831 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3832 SpatialDomains::PointGeomUniquePtr v3(
3833 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3834 SpatialDomains::PointGeomUniquePtr v4(
3835 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3836
3837 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3838 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3839 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3840 std::array<SpatialDomains::PointGeom *, 5> v = {
3841 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3842 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3843 CreatePyr(v, segVec, triVec, quadVec);
3844
3845 unsigned int numQuadPoints = 5;
3846 unsigned int numModes = 4;
3847
3848 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3849 Nektar::LibUtilities::eGaussLobattoLegendre;
3850 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3851 PointsTypeDir1);
3852 Nektar::LibUtilities::BasisType basisTypeDir1 =
3853 Nektar::LibUtilities::eModified_A;
3854 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3855 PointsKeyDir1);
3856
3857 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3858 Nektar::LibUtilities::eGaussLobattoLegendre;
3859 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3860 PointsTypeDir2);
3861 Nektar::LibUtilities::BasisType basisTypeDir2 =
3862 Nektar::LibUtilities::eModified_A;
3863 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3864 PointsKeyDir2);
3865
3866 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3867 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3868 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3869 PointsTypeDir3);
3870 Nektar::LibUtilities::BasisType basisTypeDir3 =
3871 Nektar::LibUtilities::eModifiedPyr_C;
3872 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3873 PointsKeyDir3);
3874
3875 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3876 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3877 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3878
3879 int nelmts = 10;
3880
3881 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
3882 for (int i = 0; i < nelmts; ++i)
3883 {
3884 CollExp.push_back(Exp);
3885 }
3886
3887 LibUtilities::SessionReaderSharedPtr dummySession;
3888 Collections::CollectionOptimisation colOpt(dummySession, 2,
3889 Collections::eIterPerExp);
3890 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3891 Collections::Collection c(CollExp, impTypes);
3892 StdRegions::ConstFactorMap factors;
3893 factors[StdRegions::eFactorLambda] = 1.5;
3894
3895 c.Initialise(Collections::eLinearAdvectionDiffusionReaction, factors);
3896
3897 // Add advection velocities via varcoeffs
3898 int npoints = Exp->GetTotPoints() * nelmts;
3899 StdRegions::VarCoeffMap varcoeffs;
3900 StdRegions::VarCoeffType varcoefftypes[] = {StdRegions::eVarCoeffVelX,
3901 StdRegions::eVarCoeffVelY,
3902 StdRegions::eVarCoeffVelZ};
3903 for (int i = 0; i < Exp->GetShapeDimension(); i++)
3904 {
3905 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
3906 }
3907 c.UpdateVarcoeffs(Collections::eLinearAdvectionDiffusionReaction,
3908 varcoeffs);
3909
3910 const int nm = Exp->GetNcoeffs();
3911 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3912 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3913 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3914
3915 for (int i = 0; i < coeffsIn.size(); ++i)
3916 {
3917 coeffsIn[i] = i + 1.0;
3918 }
3919
3920 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
3921 Exp->DetShapeType(), *Exp, factors,
3922 varcoeffs);
3923
3924 for (int i = 0; i < nelmts; ++i)
3925 {
3926 // Standard routines
3927 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3928 }
3929
3930 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
3931 coeffs);
3932
3933 double epsilon = 1.0e-8;
3934 for (int i = 0; i < coeffsRef.size(); ++i)
3935 {
3936 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3937 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3938 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3939 }
3940}
3941
3942BOOST_AUTO_TEST_CASE(
3943 TestPyrLinearAdvectionDiffusionReaction_MatrixFree_UniformP)
3944{
3945 SpatialDomains::PointGeomUniquePtr v0(
3946 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3947 SpatialDomains::PointGeomUniquePtr v1(
3948 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3949 SpatialDomains::PointGeomUniquePtr v2(
3950 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
3951 SpatialDomains::PointGeomUniquePtr v3(
3952 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
3953 SpatialDomains::PointGeomUniquePtr v4(
3954 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
3955
3956 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
3957 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
3958 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
3959 std::array<SpatialDomains::PointGeom *, 5> v = {
3960 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
3961 SpatialDomains::PyrGeomUniquePtr pyrGeom =
3962 CreatePyr(v, segVec, triVec, quadVec);
3963
3964 unsigned int numQuadPoints = 5;
3965 unsigned int numModes = 4;
3966
3967 Nektar::LibUtilities::PointsType PointsTypeDir1 =
3968 Nektar::LibUtilities::eGaussLobattoLegendre;
3969 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
3970 PointsTypeDir1);
3971 Nektar::LibUtilities::BasisType basisTypeDir1 =
3972 Nektar::LibUtilities::eModified_A;
3973 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3974 PointsKeyDir1);
3975
3976 Nektar::LibUtilities::PointsType PointsTypeDir2 =
3977 Nektar::LibUtilities::eGaussLobattoLegendre;
3978 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
3979 PointsTypeDir2);
3980 Nektar::LibUtilities::BasisType basisTypeDir2 =
3981 Nektar::LibUtilities::eModified_A;
3982 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3983 PointsKeyDir2);
3984
3985 Nektar::LibUtilities::PointsType PointsTypeDir3 =
3986 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3987 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
3988 PointsTypeDir3);
3989 Nektar::LibUtilities::BasisType basisTypeDir3 =
3990 Nektar::LibUtilities::eModifiedPyr_C;
3991 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3992 PointsKeyDir3);
3993
3994 Nektar::LocalRegions::PyrExpSharedPtr Exp =
3995 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
3996 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
3997
3998 int nelmts = 10;
3999
4000 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4001 for (int i = 0; i < nelmts; ++i)
4002 {
4003 CollExp.push_back(Exp);
4004 }
4005
4006 LibUtilities::SessionReaderSharedPtr dummySession;
4007 Collections::CollectionOptimisation colOpt(dummySession, 2,
4008 Collections::eMatrixFree);
4009 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4010 Collections::Collection c(CollExp, impTypes);
4011 StdRegions::ConstFactorMap factors;
4012 factors[StdRegions::eFactorLambda] = 1.5;
4013
4014 c.Initialise(Collections::eLinearAdvectionDiffusionReaction, factors);
4015
4016 // Add advection velocities via varcoeffs
4017 int npoints = Exp->GetTotPoints() * nelmts;
4018 StdRegions::VarCoeffMap varcoeffs;
4019 StdRegions::VarCoeffType varcoefftypes[] = {StdRegions::eVarCoeffVelX,
4020 StdRegions::eVarCoeffVelY,
4021 StdRegions::eVarCoeffVelZ};
4022 for (int i = 0; i < Exp->GetShapeDimension(); i++)
4023 {
4024 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
4025 }
4026 c.UpdateVarcoeffs(Collections::eLinearAdvectionDiffusionReaction,
4027 varcoeffs);
4028
4029 const int nm = Exp->GetNcoeffs();
4030 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4031 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4032 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4033
4034 for (int i = 0; i < coeffsIn.size(); ++i)
4035 {
4036 coeffsIn[i] = i + 1.0;
4037 }
4038
4039 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
4040 Exp->DetShapeType(), *Exp, factors,
4041 varcoeffs);
4042
4043 for (int i = 0; i < nelmts; ++i)
4044 {
4045 // Standard routines
4046 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4047 }
4048
4049 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
4050 coeffs);
4051
4052 double epsilon = 1.0e-8;
4053 for (int i = 0; i < coeffsRef.size(); ++i)
4054 {
4055 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4056 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4057 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4058 }
4059}
4060
4061BOOST_AUTO_TEST_CASE(
4062 TestPyrLinearAdvectionDiffusionReaction_MatrixFree_Deformed_overInt)
4063{
4064 SpatialDomains::PointGeomUniquePtr v0(
4065 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
4066 SpatialDomains::PointGeomUniquePtr v1(
4067 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
4068 SpatialDomains::PointGeomUniquePtr v2(
4069 new SpatialDomains::PointGeom(3u, 2u, 1.0, 1.0, -1.0));
4070 SpatialDomains::PointGeomUniquePtr v3(
4071 new SpatialDomains::PointGeom(3u, 3u, -1.0, 1.0, -1.0));
4072 SpatialDomains::PointGeomUniquePtr v4(
4073 new SpatialDomains::PointGeom(3u, 4u, -1.0, -1.0, 1.0));
4074
4075 std::array<SpatialDomains::SegGeomUniquePtr, 8> segVec;
4076 std::array<SpatialDomains::QuadGeomUniquePtr, 1> quadVec;
4077 std::array<SpatialDomains::TriGeomUniquePtr, 4> triVec;
4078 std::array<SpatialDomains::PointGeom *, 5> v = {
4079 v0.get(), v1.get(), v2.get(), v3.get(), v4.get()};
4080 SpatialDomains::PyrGeomUniquePtr pyrGeom =
4081 CreatePyr(v, segVec, triVec, quadVec);
4082
4083 unsigned int numQuadPoints = 8;
4084 unsigned int numModes = 4;
4085
4086 Nektar::LibUtilities::PointsType PointsTypeDir1 =
4087 Nektar::LibUtilities::eGaussLobattoLegendre;
4088 const Nektar::LibUtilities::PointsKey PointsKeyDir1(numQuadPoints,
4089 PointsTypeDir1);
4090 Nektar::LibUtilities::BasisType basisTypeDir1 =
4091 Nektar::LibUtilities::eModified_A;
4092 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
4093 PointsKeyDir1);
4094
4095 Nektar::LibUtilities::PointsType PointsTypeDir2 =
4096 Nektar::LibUtilities::eGaussLobattoLegendre;
4097 const Nektar::LibUtilities::PointsKey PointsKeyDir2(numQuadPoints,
4098 PointsTypeDir2);
4099 Nektar::LibUtilities::BasisType basisTypeDir2 =
4100 Nektar::LibUtilities::eModified_A;
4101 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
4102 PointsKeyDir2);
4103
4104 Nektar::LibUtilities::PointsType PointsTypeDir3 =
4105 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
4106 const Nektar::LibUtilities::PointsKey PointsKeyDir3(numQuadPoints - 1,
4107 PointsTypeDir3);
4108 Nektar::LibUtilities::BasisType basisTypeDir3 =
4109 Nektar::LibUtilities::eModifiedPyr_C;
4110 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
4111 PointsKeyDir3);
4112
4113 Nektar::LocalRegions::PyrExpSharedPtr Exp =
4114 MemoryManager<Nektar::LocalRegions::PyrExp>::AllocateSharedPtr(
4115 basisKeyDir1, basisKeyDir2, basisKeyDir3, pyrGeom.get());
4116
4117 int nelmts = 10;
4118
4119 std::vector<LocalRegions::ExpansionSharedPtr> CollExp;
4120 for (int i = 0; i < nelmts; ++i)
4121 {
4122 CollExp.push_back(Exp);
4123 }
4124
4125 LibUtilities::SessionReaderSharedPtr dummySession;
4126 Collections::CollectionOptimisation colOpt(dummySession, 2,
4127 Collections::eMatrixFree);
4128 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
4129 Collections::Collection c(CollExp, impTypes);
4130 StdRegions::ConstFactorMap factors;
4131 factors[StdRegions::eFactorLambda] = 1.5;
4132
4133 c.Initialise(Collections::eLinearAdvectionDiffusionReaction, factors);
4134
4135 // Add advection velocities via varcoeffs
4136 int npoints = Exp->GetTotPoints() * nelmts;
4137 StdRegions::VarCoeffMap varcoeffs;
4138 StdRegions::VarCoeffType varcoefftypes[] = {StdRegions::eVarCoeffVelX,
4139 StdRegions::eVarCoeffVelY,
4140 StdRegions::eVarCoeffVelZ};
4141 for (int i = 0; i < Exp->GetShapeDimension(); i++)
4142 {
4143 varcoeffs[varcoefftypes[i]] = Array<OneD, NekDouble>(npoints, 1.0);
4144 }
4145 c.UpdateVarcoeffs(Collections::eLinearAdvectionDiffusionReaction,
4146 varcoeffs);
4147
4148 const int nm = Exp->GetNcoeffs();
4149 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
4150 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
4151 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
4152
4153 for (int i = 0; i < coeffsIn.size(); ++i)
4154 {
4155 coeffsIn[i] = i + 1.0;
4156 }
4157
4158 StdRegions::StdMatrixKey mkey(StdRegions::eLinearAdvectionDiffusionReaction,
4159 Exp->DetShapeType(), *Exp, factors,
4160 varcoeffs);
4161
4162 for (int i = 0; i < nelmts; ++i)
4163 {
4164 // Standard routines
4165 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
4166 }
4167
4168 c.ApplyOperator(Collections::eLinearAdvectionDiffusionReaction, coeffsIn,
4169 coeffs);
4170
4171 double epsilon = 1.0e-8;
4172 for (int i = 0; i < coeffsRef.size(); ++i)
4173 {
4174 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
4175 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
4176 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
4177 }
4178}
4179
4180} // namespace Nektar::PyrCollectionTests
NELMTS
#define NELMTS
Definition TestPyrCollection.cpp:46
Nektar::Collections::Collection::UpdateVarcoeffs
COLLECTIONS_EXPORT void UpdateVarcoeffs(const OperatorType opType, StdRegions::VarCoeffMap &varcoeffs=StdRegions::NullVarCoeffMap)
Definition Collection.cpp:60
Nektar::Collections::Collection::GetOutputSize
int GetOutputSize(const OperatorType &op)
Definition Collection.h:112
Nektar::Collections::Collection::Initialise
COLLECTIONS_EXPORT void Initialise(const OperatorType opType, StdRegions::FactorMap factors=StdRegions::NullFactorMap)
Definition Collection.cpp:66
Nektar::Collections::Collection::ApplyOperator
void ApplyOperator(const OperatorType &op, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &output)
Definition Collection.h:148
Nektar::Collections::CollectionOptimisation::GetOperatorImpMap
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(LocalRegions::ExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Definition CollectionOptimisation.cpp:356
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition Basis.h:45
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition Points.h:50
Nektar::LibUtilities::Timer::AccumulateRegion
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition Timer.cpp:70
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition NekMemoryManager.hpp:167
Nektar::Collections::OperatorImpMap
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition Operator.h:131
Nektar::Collections::eLinearAdvectionDiffusionReaction
@ eLinearAdvectionDiffusionReaction
Definition Operator.h:66
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition SessionReader.h:119
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition PointsType.h:51
Nektar::LibUtilities::eModifiedPyr_C
@ eModifiedPyr_C
Principle Modified Functions.
Definition BasisType.h:53
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition BasisType.h:48
Nektar::LocalRegions::PyrExpSharedPtr
std::shared_ptr< PyrExp > PyrExpSharedPtr
Definition PyrExp.h:139
Nektar::PyrCollectionTests::CreateSegGeom
SpatialDomains::SegGeomUniquePtr CreateSegGeom(unsigned int id, SpatialDomains::PointGeom *v0, SpatialDomains::PointGeom *v1)
Definition TestPyrCollection.cpp:48
Nektar::PyrCollectionTests::CreatePyr
SpatialDomains::PyrGeomUniquePtr CreatePyr(std::array< SpatialDomains::PointGeom *, 5 > v, std::array< SpatialDomains::SegGeomUniquePtr, 8 > &segVec, std::array< SpatialDomains::TriGeomUniquePtr, 4 > &triVec, std::array< SpatialDomains::QuadGeomUniquePtr, 1 > &quadVec)
Definition TestPyrCollection.cpp:58
Nektar::PyrCollectionTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestPyrBwdTrans_IterPerExp_UniformP_MultiElmt)
Definition TestPyrCollection.cpp:119
Nektar::SpatialDomains::QuadGeomUniquePtr
unique_ptr_objpool< QuadGeom > QuadGeomUniquePtr
Definition MeshGraph.h:104
Nektar::SpatialDomains::PyrGeomUniquePtr
unique_ptr_objpool< PyrGeom > PyrGeomUniquePtr
Definition MeshGraph.h:106
Nektar::SpatialDomains::SegGeomUniquePtr
unique_ptr_objpool< SegGeom > SegGeomUniquePtr
Definition MeshGraph.h:102
Nektar::SpatialDomains::PointGeomUniquePtr
unique_ptr_objpool< PointGeom > PointGeomUniquePtr
Definition MeshGraph.h:99
Nektar::SpatialDomains::TriGeomUniquePtr
unique_ptr_objpool< TriGeom > TriGeomUniquePtr
Definition MeshGraph.h:103
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition StdRegions.hpp:489
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition StdRegions.hpp:432
Vmath::Vadd
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
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition Vmath.hpp:825
Generated on Sat Jul 4 2026 01:02:24 for Nektar++ by doxygen 1.9.8