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