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