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