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