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