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