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