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