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