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