Nektar++
TestSegCollection.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TestSegCollection.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/SegExp.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::SegCollectionTests
43{
44SpatialDomains::SegGeomSharedPtr CreateSegGeom(
45 unsigned int id, SpatialDomains::PointGeomSharedPtr v0,
46 SpatialDomains::PointGeomSharedPtr v1, int coordim = 1)
47{
48 SpatialDomains::PointGeomSharedPtr vertices[] = {v0, v1};
49 SpatialDomains::SegGeomSharedPtr result(
50 new SpatialDomains::SegGeom(id, coordim, vertices));
51 return result;
52}
53
54BOOST_AUTO_TEST_CASE(TestSegBwdTrans_StdMat_UniformP)
55{
56 SpatialDomains::PointGeomSharedPtr v0(
57 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
58 SpatialDomains::PointGeomSharedPtr v1(
59 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
60
61 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
62
63 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
64 Nektar::LibUtilities::eGaussLobattoLegendre;
65 Nektar::LibUtilities::BasisType basisTypeDir1 =
66 Nektar::LibUtilities::eModified_A;
67 unsigned int numSegPoints = 6;
68 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
69 segPointsTypeDir1);
70 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
71 segPointsKeyDir1);
72
73 Nektar::LocalRegions::SegExpSharedPtr Exp =
74 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
75 basisKeyDir1, segGeom);
76
77 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
78 CollExp.push_back(Exp);
79
80 LibUtilities::SessionReaderSharedPtr dummySession;
81 Collections::CollectionOptimisation colOpt(dummySession, 1,
82 Collections::eStdMat);
83 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
84 Collections::Collection c(CollExp, impTypes);
85 c.Initialise(Collections::eBwdTrans);
86
87 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
88 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
89 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
90
91 Exp->BwdTrans(coeffs, phys1);
92 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
93
94 double epsilon = 1.0e-8;
95 for (int i = 0; i < phys1.size(); ++i)
96 {
97 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
98 }
99}
100
101BOOST_AUTO_TEST_CASE(TestSegBwdTrans_StdMat_UniformP_MultiElmt)
102{
103 SpatialDomains::PointGeomSharedPtr v0(
104 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
105 SpatialDomains::PointGeomSharedPtr v1(
106 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
107
108 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
109
110 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
111 Nektar::LibUtilities::eGaussLobattoLegendre;
112 Nektar::LibUtilities::BasisType basisTypeDir1 =
113 Nektar::LibUtilities::eModified_A;
114 unsigned int numSegPoints = 6;
115 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
116 segPointsTypeDir1);
117 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
118 segPointsKeyDir1);
119
120 Nektar::LocalRegions::SegExpSharedPtr Exp =
121 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
122 basisKeyDir1, segGeom);
123
124 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
125
126 int nelmts = 10;
127 for (int i = 0; i < nelmts; ++i)
128 {
129 CollExp.push_back(Exp);
130 }
131
132 LibUtilities::SessionReaderSharedPtr dummySession;
133 Collections::CollectionOptimisation colOpt(dummySession, 1,
134 Collections::eStdMat);
135 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
136 Collections::Collection c(CollExp, impTypes);
137 c.Initialise(Collections::eBwdTrans);
138
139 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
140 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
141 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
142
143 for (int i = 0; i < nelmts; ++i)
144 {
145 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
146 tmp = phys1 + i * Exp->GetTotPoints());
147 }
148 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
149
150 double epsilon = 1.0e-8;
151 for (int i = 0; i < phys1.size(); ++i)
152 {
153 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
154 }
155}
156
157BOOST_AUTO_TEST_CASE(TestSegBwdTrans_IterPerExp_UniformP)
158{
159 SpatialDomains::PointGeomSharedPtr v0(
160 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
161 SpatialDomains::PointGeomSharedPtr v1(
162 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
163
164 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
165
166 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
167 Nektar::LibUtilities::eGaussLobattoLegendre;
168 Nektar::LibUtilities::BasisType basisTypeDir1 =
169 Nektar::LibUtilities::eModified_A;
170 unsigned int numSegPoints = 6;
171 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
172 segPointsTypeDir1);
173 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
174 segPointsKeyDir1);
175
176 Nektar::LocalRegions::SegExpSharedPtr Exp =
177 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
178 basisKeyDir1, segGeom);
179
180 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
181 CollExp.push_back(Exp);
182
183 LibUtilities::SessionReaderSharedPtr dummySession;
184 Collections::CollectionOptimisation colOpt(dummySession, 1,
185 Collections::eIterPerExp);
186 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
187 Collections::Collection c(CollExp, impTypes);
188 c.Initialise(Collections::eBwdTrans);
189
190 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
191 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
192 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
193
194 Exp->BwdTrans(coeffs, phys1);
195 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
196
197 double epsilon = 1.0e-8;
198 for (int i = 0; i < phys1.size(); ++i)
199 {
200 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
201 }
202}
203
204BOOST_AUTO_TEST_CASE(TestSegBwdTrans_SumFac_UniformP)
205{
206 SpatialDomains::PointGeomSharedPtr v0(
207 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
208 SpatialDomains::PointGeomSharedPtr v1(
209 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
210
211 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
212
213 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
214 Nektar::LibUtilities::eGaussLobattoLegendre;
215 Nektar::LibUtilities::BasisType basisTypeDir1 =
216 Nektar::LibUtilities::eModified_A;
217 unsigned int numSegPoints = 6;
218 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
219 segPointsTypeDir1);
220 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
221 segPointsKeyDir1);
222
223 Nektar::LocalRegions::SegExpSharedPtr Exp =
224 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
225 basisKeyDir1, segGeom);
226
227 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
228
229 int nelmts = 1;
230 for (int i = 0; i < nelmts; ++i)
231 {
232 CollExp.push_back(Exp);
233 }
234
235 LibUtilities::SessionReaderSharedPtr dummySession;
236 Collections::CollectionOptimisation colOpt(dummySession, 1,
237 Collections::eSumFac);
238 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
239 Collections::Collection c(CollExp, impTypes);
240 c.Initialise(Collections::eBwdTrans);
241
242 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
243 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
244 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
245
246 for (int i = 0; i < nelmts; ++i)
247 {
248 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
249 tmp = phys1 + i * Exp->GetTotPoints());
250 }
251 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
252
253 double epsilon = 1.0e-8;
254 for (int i = 0; i < phys1.size(); ++i)
255 {
256 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
257 }
258}
259
260BOOST_AUTO_TEST_CASE(TestSegBwdTrans_SumFac_UniformP_MultiElmt)
261{
262 SpatialDomains::PointGeomSharedPtr v0(
263 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
264 SpatialDomains::PointGeomSharedPtr v1(
265 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
266
267 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
268
269 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
270 Nektar::LibUtilities::eGaussLobattoLegendre;
271 Nektar::LibUtilities::BasisType basisTypeDir1 =
272 Nektar::LibUtilities::eModified_A;
273 unsigned int numSegPoints = 6;
274 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
275 segPointsTypeDir1);
276 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
277 segPointsKeyDir1);
278
279 Nektar::LocalRegions::SegExpSharedPtr Exp =
280 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
281 basisKeyDir1, segGeom);
282
283 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
284
285 int nelmts = 10;
286 for (int i = 0; i < nelmts; ++i)
287 {
288 CollExp.push_back(Exp);
289 }
290
291 LibUtilities::SessionReaderSharedPtr dummySession;
292 Collections::CollectionOptimisation colOpt(dummySession, 1,
293 Collections::eSumFac);
294 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
295 Collections::Collection c(CollExp, impTypes);
296 c.Initialise(Collections::eBwdTrans);
297
298 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
299 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
300 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
301
302 for (int i = 0; i < nelmts; ++i)
303 {
304 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
305 tmp = phys1 + i * Exp->GetTotPoints());
306 }
307 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
308
309 double epsilon = 1.0e-8;
310 for (int i = 0; i < phys1.size(); ++i)
311 {
312 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
313 }
314}
315
316BOOST_AUTO_TEST_CASE(TestSegBwdTrans_MatrixFree_UniformP_MultiElmt)
317{
318 SpatialDomains::PointGeomSharedPtr v0(
319 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
320 SpatialDomains::PointGeomSharedPtr v1(
321 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
322
323 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
324
325 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
326 Nektar::LibUtilities::eGaussLobattoLegendre;
327 Nektar::LibUtilities::BasisType basisTypeDir1 =
328 Nektar::LibUtilities::eModified_A;
329 unsigned int numSegPoints = 6;
330 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
331 segPointsTypeDir1);
332 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
333 segPointsKeyDir1);
334
335 Nektar::LocalRegions::SegExpSharedPtr Exp =
336 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
337 basisKeyDir1, segGeom);
338
339 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
340
341 int nelmts = 10;
342 for (int i = 0; i < nelmts; ++i)
343 {
344 CollExp.push_back(Exp);
345 }
346
347 LibUtilities::SessionReaderSharedPtr dummySession;
348 Collections::CollectionOptimisation colOpt(dummySession, 1,
349 Collections::eIterPerExp);
350 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
351 impTypes[Collections::eBwdTrans] = Collections::eMatrixFree;
352 Collections::Collection c(CollExp, impTypes);
353 c.Initialise(Collections::eBwdTrans);
354
355 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
356 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
357 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
358
359 for (int i = 0; i < nelmts; ++i)
360 {
361 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
362 tmp = phys1 + i * Exp->GetTotPoints());
363 }
364 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
365
366 double epsilon = 1.0e-8;
367 for (int i = 0; i < phys1.size(); ++i)
368 {
369 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
370 }
371}
372
373BOOST_AUTO_TEST_CASE(TestSegIProductWRTBase_IterPerExp_UniformP_MultiElmt)
374{
375 SpatialDomains::PointGeomSharedPtr v0(
376 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
377 SpatialDomains::PointGeomSharedPtr v1(
378 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
379
380 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
381
382 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
383 Nektar::LibUtilities::eGaussLobattoLegendre;
384 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
385 segPointsTypeDir1);
386 Nektar::LibUtilities::BasisType basisTypeDir1 =
387 Nektar::LibUtilities::eModified_A;
388 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
389 segPointsKeyDir1);
390
391 Nektar::LocalRegions::SegExpSharedPtr Exp =
392 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
393 basisKeyDir1, segGeom);
394
395 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
396
397 int nelmts = 10;
398 for (int i = 0; i < nelmts; ++i)
399 {
400 CollExp.push_back(Exp);
401 }
402
403 LibUtilities::SessionReaderSharedPtr dummySession;
404 Collections::CollectionOptimisation colOpt(dummySession, 1,
405 Collections::eIterPerExp);
406 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
407 Collections::Collection c(CollExp, impTypes);
408 c.Initialise(Collections::eIProductWRTBase);
409
410 const int nq = Exp->GetTotPoints();
411 const int nc = Exp->GetNcoeffs();
412 Array<OneD, NekDouble> xc(nq), yc(nq);
413 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
414 Array<OneD, NekDouble> coeffs1(nelmts * nc);
415 Array<OneD, NekDouble> coeffs2(nelmts * nc);
416
417 Exp->GetCoords(xc, yc);
418
419 for (int i = 0; i < nq; ++i)
420 {
421 phys[i] = sin(xc[i]);
422 }
423 Exp->IProductWRTBase(phys, coeffs1);
424
425 for (int i = 1; i < nelmts; ++i)
426 {
427 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
428 Exp->IProductWRTBase(phys + i * nq, tmp = coeffs1 + i * nc);
429 }
430 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
431
432 double epsilon = 1.0e-8;
433 for (int i = 0; i < coeffs1.size(); ++i)
434 {
435 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
436 }
437}
438
439BOOST_AUTO_TEST_CASE(TestSegIProductWRTBase_StdMat_UniformP_MultiElmt)
440{
441 SpatialDomains::PointGeomSharedPtr v0(
442 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
443 SpatialDomains::PointGeomSharedPtr v1(
444 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
445
446 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
447
448 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
449 Nektar::LibUtilities::eGaussLobattoLegendre;
450 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
451 segPointsTypeDir1);
452 Nektar::LibUtilities::BasisType basisTypeDir1 =
453 Nektar::LibUtilities::eModified_A;
454 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
455 segPointsKeyDir1);
456
457 Nektar::LocalRegions::SegExpSharedPtr Exp =
458 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
459 basisKeyDir1, segGeom);
460
461 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
462
463 int nelmts = 10;
464 for (int i = 0; i < nelmts; ++i)
465 {
466 CollExp.push_back(Exp);
467 }
468
469 LibUtilities::SessionReaderSharedPtr dummySession;
470 Collections::CollectionOptimisation colOpt(dummySession, 1,
471 Collections::eStdMat);
472 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
473 Collections::Collection c(CollExp, impTypes);
474 c.Initialise(Collections::eIProductWRTBase);
475
476 const int nq = Exp->GetTotPoints();
477 const int nc = Exp->GetNcoeffs();
478 Array<OneD, NekDouble> xc(nq), yc(nq);
479 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
480 Array<OneD, NekDouble> coeffs1(nelmts * nc);
481 Array<OneD, NekDouble> coeffs2(nelmts * nc);
482
483 Exp->GetCoords(xc, yc);
484
485 for (int i = 0; i < nq; ++i)
486 {
487 phys[i] = sin(xc[i]);
488 }
489 Exp->IProductWRTBase(phys, coeffs1);
490
491 for (int i = 1; i < nelmts; ++i)
492 {
493 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
494 Exp->IProductWRTBase(phys + i * nq, tmp = coeffs1 + i * nc);
495 }
496 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
497
498 double epsilon = 1.0e-8;
499 for (int i = 0; i < coeffs1.size(); ++i)
500 {
501 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
502 }
503}
504
505BOOST_AUTO_TEST_CASE(TestSegIProductWRTBase_SumFac_UniformP_MultiElmt)
506{
507 SpatialDomains::PointGeomSharedPtr v0(
508 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
509 SpatialDomains::PointGeomSharedPtr v1(
510 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
511
512 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
513
514 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
515 Nektar::LibUtilities::eGaussLobattoLegendre;
516 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
517 segPointsTypeDir1);
518 Nektar::LibUtilities::BasisType basisTypeDir1 =
519 Nektar::LibUtilities::eModified_A;
520 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
521 segPointsKeyDir1);
522
523 Nektar::LocalRegions::SegExpSharedPtr Exp =
524 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
525 basisKeyDir1, segGeom);
526
527 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
528
529 int nelmts = 10;
530 for (int i = 0; i < nelmts; ++i)
531 {
532 CollExp.push_back(Exp);
533 }
534
535 LibUtilities::SessionReaderSharedPtr dummySession;
536 Collections::CollectionOptimisation colOpt(dummySession, 1,
537 Collections::eSumFac);
538 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
539 Collections::Collection c(CollExp, impTypes);
540 c.Initialise(Collections::eIProductWRTBase);
541
542 const int nq = Exp->GetTotPoints();
543 const int nc = Exp->GetNcoeffs();
544 Array<OneD, NekDouble> xc(nq), yc(nq);
545 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
546 Array<OneD, NekDouble> coeffs1(nelmts * nc);
547 Array<OneD, NekDouble> coeffs2(nelmts * nc);
548
549 Exp->GetCoords(xc, yc);
550
551 for (int i = 0; i < nq; ++i)
552 {
553 phys[i] = sin(xc[i]);
554 }
555 Exp->IProductWRTBase(phys, coeffs1);
556
557 for (int i = 1; i < nelmts; ++i)
558 {
559 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
560 Exp->IProductWRTBase(phys + i * nq, tmp = coeffs1 + i * nc);
561 }
562 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
563
564 double epsilon = 1.0e-8;
565 for (int i = 0; i < coeffs1.size(); ++i)
566 {
567 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
568 }
569}
570
571BOOST_AUTO_TEST_CASE(TestSegIProductWRTBase_MatrixFree_UniformP_MultiElmt)
572{
573 SpatialDomains::PointGeomSharedPtr v0(
574 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
575 SpatialDomains::PointGeomSharedPtr v1(
576 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
577
578 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
579
580 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
581 Nektar::LibUtilities::eGaussLobattoLegendre;
582 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
583 segPointsTypeDir1);
584 Nektar::LibUtilities::BasisType basisTypeDir1 =
585 Nektar::LibUtilities::eModified_A;
586 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
587 segPointsKeyDir1);
588
589 Nektar::LocalRegions::SegExpSharedPtr Exp =
590 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
591 basisKeyDir1, segGeom);
592
593 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
594
595 int nelmts = 10;
596 for (int i = 0; i < nelmts; ++i)
597 {
598 CollExp.push_back(Exp);
599 }
600
601 LibUtilities::SessionReaderSharedPtr dummySession;
602 Collections::CollectionOptimisation colOpt(dummySession, 1,
603 Collections::eIterPerExp);
604 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
605 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
606 Collections::Collection c(CollExp, impTypes);
607 c.Initialise(Collections::eIProductWRTBase);
608
609 const int nq = Exp->GetTotPoints();
610 const int nc = Exp->GetNcoeffs();
611 Array<OneD, NekDouble> xc(nq), yc(nq);
612 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
613 Array<OneD, NekDouble> coeffs1(nelmts * nc);
614 Array<OneD, NekDouble> coeffs2(nelmts * nc);
615
616 Exp->GetCoords(xc, yc);
617
618 for (int i = 0; i < nq; ++i)
619 {
620 phys[i] = sin(xc[i]);
621 }
622 Exp->IProductWRTBase(phys, coeffs1);
623
624 for (int i = 1; i < nelmts; ++i)
625 {
626 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
627 Exp->IProductWRTBase(phys + i * nq, tmp = coeffs1 + i * nc);
628 }
629 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
630
631 double epsilon = 1.0e-8;
632 for (int i = 0; i < coeffs1.size(); ++i)
633 {
634 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
635 }
636}
637
638BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_IterPerExp_UniformP)
639{
640 SpatialDomains::PointGeomSharedPtr v0(
641 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
642 SpatialDomains::PointGeomSharedPtr v1(
643 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
644
645 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
646
647 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
648 Nektar::LibUtilities::eGaussLobattoLegendre;
649 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
650 segPointsTypeDir1);
651 Nektar::LibUtilities::BasisType basisTypeDir1 =
652 Nektar::LibUtilities::eModified_A;
653 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
654 segPointsKeyDir1);
655
656 Nektar::LocalRegions::SegExpSharedPtr Exp =
657 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
658 basisKeyDir1, segGeom);
659
660 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
661 CollExp.push_back(Exp);
662
663 LibUtilities::SessionReaderSharedPtr dummySession;
664 Collections::CollectionOptimisation colOpt(dummySession, 1,
665 Collections::eIterPerExp);
666 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
667 Collections::Collection c(CollExp, impTypes);
668 c.Initialise(Collections::ePhysDeriv);
669
670 const int nq = Exp->GetTotPoints();
671 Array<OneD, NekDouble> xc(nq), yc(nq);
672 Array<OneD, NekDouble> phys(nq), tmp;
673 Array<OneD, NekDouble> diff1(nq);
674 Array<OneD, NekDouble> diff2(nq);
675
676 Exp->GetCoords(xc, yc);
677
678 for (int i = 0; i < nq; ++i)
679 {
680 phys[i] = sin(xc[i]);
681 }
682
683 Exp->PhysDeriv(phys, diff1);
684 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
685
686 double epsilon = 1.0e-8;
687 for (int i = 0; i < diff1.size(); ++i)
688 {
689 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
690 }
691}
692
693BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_IterPerExp_UniformP_MultiElmt)
694{
695 SpatialDomains::PointGeomSharedPtr v0(
696 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
697 SpatialDomains::PointGeomSharedPtr v1(
698 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
699
700 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
701
702 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
703 Nektar::LibUtilities::eGaussLobattoLegendre;
704 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
705 segPointsTypeDir1);
706 Nektar::LibUtilities::BasisType basisTypeDir1 =
707 Nektar::LibUtilities::eModified_A;
708 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
709 segPointsKeyDir1);
710
711 Nektar::LocalRegions::SegExpSharedPtr Exp =
712 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
713 basisKeyDir1, segGeom);
714
715 int nelmts = 10;
716
717 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
718 for (int i = 0; i < nelmts; ++i)
719 {
720 CollExp.push_back(Exp);
721 }
722
723 LibUtilities::SessionReaderSharedPtr dummySession;
724 Collections::CollectionOptimisation colOpt(dummySession, 1,
725 Collections::eIterPerExp);
726 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
727 Collections::Collection c(CollExp, impTypes);
728 c.Initialise(Collections::ePhysDeriv);
729
730 const int nq = Exp->GetTotPoints();
731 Array<OneD, NekDouble> xc(nq);
732 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
733 Array<OneD, NekDouble> diff1(nelmts * nq);
734 Array<OneD, NekDouble> diff2(nelmts * nq);
735
736 Exp->GetCoords(xc);
737
738 for (int i = 0; i < nq; ++i)
739 {
740 phys[i] = sin(xc[i]);
741 }
742 Exp->PhysDeriv(phys, diff1);
743
744 for (int i = 1; i < nelmts; ++i)
745 {
746 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
747 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq);
748 }
749
750 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
751
752 double epsilon = 1.0e-8;
753 for (int i = 0; i < diff1.size(); ++i)
754 {
755 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
756 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
757 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
758 }
759}
760
761BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_StdMat_UniformP)
762{
763 SpatialDomains::PointGeomSharedPtr v0(
764 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
765 SpatialDomains::PointGeomSharedPtr v1(
766 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
767
768 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
769
770 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
771 Nektar::LibUtilities::eGaussLobattoLegendre;
772 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
773 segPointsTypeDir1);
774 Nektar::LibUtilities::BasisType basisTypeDir1 =
775 Nektar::LibUtilities::eModified_A;
776 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
777 segPointsKeyDir1);
778
779 Nektar::LocalRegions::SegExpSharedPtr Exp =
780 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
781 basisKeyDir1, segGeom);
782
783 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
784 CollExp.push_back(Exp);
785
786 LibUtilities::SessionReaderSharedPtr dummySession;
787 Collections::CollectionOptimisation colOpt(dummySession, 1,
788 Collections::eStdMat);
789 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
790 Collections::Collection c(CollExp, impTypes);
791 c.Initialise(Collections::ePhysDeriv);
792
793 const int nq = Exp->GetTotPoints();
794 Array<OneD, NekDouble> xc(nq);
795 Array<OneD, NekDouble> phys(nq), tmp;
796 Array<OneD, NekDouble> diff1(nq);
797 Array<OneD, NekDouble> diff2(nq);
798
799 Exp->GetCoords(xc);
800
801 for (int i = 0; i < nq; ++i)
802 {
803 phys[i] = sin(xc[i]);
804 }
805
806 Exp->PhysDeriv(phys, diff1);
807 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
808
809 double epsilon = 1.0e-8;
810 for (int i = 0; i < diff1.size(); ++i)
811 {
812 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
813 }
814}
815
816BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_StdMat_UniformP_MultiElmt)
817{
818 SpatialDomains::PointGeomSharedPtr v0(
819 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
820 SpatialDomains::PointGeomSharedPtr v1(
821 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
822
823 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
824
825 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
826 Nektar::LibUtilities::eGaussLobattoLegendre;
827 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
828 segPointsTypeDir1);
829 Nektar::LibUtilities::BasisType basisTypeDir1 =
830 Nektar::LibUtilities::eModified_A;
831 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
832 segPointsKeyDir1);
833
834 Nektar::LocalRegions::SegExpSharedPtr Exp =
835 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
836 basisKeyDir1, segGeom);
837
838 int nelmts = 10;
839
840 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
841 for (int i = 0; i < nelmts; ++i)
842 {
843 CollExp.push_back(Exp);
844 }
845
846 LibUtilities::SessionReaderSharedPtr dummySession;
847 Collections::CollectionOptimisation colOpt(dummySession, 1,
848 Collections::eStdMat);
849 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
850 Collections::Collection c(CollExp, impTypes);
851 c.Initialise(Collections::ePhysDeriv);
852
853 const int nq = Exp->GetTotPoints();
854 Array<OneD, NekDouble> xc(nq);
855 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
856 Array<OneD, NekDouble> diff1(nelmts * nq);
857 Array<OneD, NekDouble> diff2(nelmts * nq);
858
859 Exp->GetCoords(xc);
860
861 for (int i = 0; i < nq; ++i)
862 {
863 phys[i] = sin(xc[i]);
864 }
865 Exp->PhysDeriv(phys, diff1);
866
867 for (int i = 1; i < nelmts; ++i)
868 {
869 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
870 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq);
871 }
872
873 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
874
875 double epsilon = 1.0e-8;
876 for (int i = 0; i < diff1.size(); ++i)
877 {
878 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
879 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
880 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
881 }
882}
883
884BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_SumFac_UniformP_MultiElmt)
885{
886 SpatialDomains::PointGeomSharedPtr v0(
887 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
888 SpatialDomains::PointGeomSharedPtr v1(
889 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
890
891 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
892
893 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
894 Nektar::LibUtilities::eGaussLobattoLegendre;
895 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
896 segPointsTypeDir1);
897 Nektar::LibUtilities::BasisType basisTypeDir1 =
898 Nektar::LibUtilities::eModified_A;
899 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
900 segPointsKeyDir1);
901
902 Nektar::LocalRegions::SegExpSharedPtr Exp =
903 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
904 basisKeyDir1, segGeom);
905
906 int nelmts = 10;
907
908 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
909 for (int i = 0; i < nelmts; ++i)
910 {
911 CollExp.push_back(Exp);
912 }
913
914 LibUtilities::SessionReaderSharedPtr dummySession;
915 Collections::CollectionOptimisation colOpt(dummySession, 1,
916 Collections::eSumFac);
917 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
918 Collections::Collection c(CollExp, impTypes);
919 c.Initialise(Collections::ePhysDeriv);
920
921 const int nq = Exp->GetTotPoints();
922 Array<OneD, NekDouble> xc(nq);
923 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
924 Array<OneD, NekDouble> diff1(nelmts * nq);
925 Array<OneD, NekDouble> diff2(nelmts * nq);
926
927 Exp->GetCoords(xc);
928
929 for (int i = 0; i < nq; ++i)
930 {
931 phys[i] = sin(xc[i]);
932 }
933 Exp->PhysDeriv(phys, diff1);
934
935 for (int i = 1; i < nelmts; ++i)
936 {
937 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
938 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq);
939 }
940
941 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
942
943 double epsilon = 1.0e-8;
944 for (int i = 0; i < diff1.size(); ++i)
945 {
946 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
947 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
948 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
949 }
950}
951
952BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_MatrixFree_UniformP_MultiElmt_1D)
953{
954 SpatialDomains::PointGeomSharedPtr v0(
955 new SpatialDomains::PointGeom(1u, 0u, -1.5, 0.0, 0.0));
956 SpatialDomains::PointGeomSharedPtr v1(
957 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
958
959 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
960
961 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
962 Nektar::LibUtilities::eGaussLobattoLegendre;
963 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
964 segPointsTypeDir1);
965 Nektar::LibUtilities::BasisType basisTypeDir1 =
966 Nektar::LibUtilities::eModified_A;
967 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
968 segPointsKeyDir1);
969
970 Nektar::LocalRegions::SegExpSharedPtr Exp =
971 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
972 basisKeyDir1, segGeom);
973
974 int nelmts = 10;
975
976 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
977 for (int i = 0; i < nelmts; ++i)
978 {
979 CollExp.push_back(Exp);
980 }
981
982 LibUtilities::SessionReaderSharedPtr dummySession;
983 Collections::CollectionOptimisation colOpt(dummySession, 1,
984 Collections::eIterPerExp);
985 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
986 impTypes[Collections::ePhysDeriv] = Collections::eMatrixFree;
987 Collections::Collection c(CollExp, impTypes);
988 c.Initialise(Collections::ePhysDeriv);
989
990 const int nq = Exp->GetTotPoints();
991 Array<OneD, NekDouble> xc(nq);
992 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
993 Array<OneD, NekDouble> diff1(nelmts * nq);
994 Array<OneD, NekDouble> diff2(nelmts * nq);
995
996 Exp->GetCoords(xc);
997
998 for (int i = 0; i < nq; ++i)
999 {
1000 phys[i] = sin(xc[i]);
1001 }
1002 Exp->PhysDeriv(phys, diff1);
1003
1004 for (int i = 1; i < nelmts; ++i)
1005 {
1006 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1007 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq);
1008 }
1009
1010 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2);
1011
1012 double epsilon = 1.0e-8;
1013 for (int i = 0; i < diff1.size(); ++i)
1014 {
1015 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1016 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1017 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1018 }
1019}
1020
1021BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_MatrixFree_UniformP_MultiElmt_2D)
1022{
1023 SpatialDomains::PointGeomSharedPtr v0(
1024 new SpatialDomains::PointGeom(2u, 0u, -1.5, 0.0, 0.0));
1025 SpatialDomains::PointGeomSharedPtr v1(
1026 new SpatialDomains::PointGeom(2u, 1u, 1.0, 1.0, 0.0));
1027
1028 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1, 2u);
1029
1030 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1031 Nektar::LibUtilities::eGaussLobattoLegendre;
1032 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
1033 segPointsTypeDir1);
1034 Nektar::LibUtilities::BasisType basisTypeDir1 =
1035 Nektar::LibUtilities::eModified_A;
1036 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1037 segPointsKeyDir1);
1038
1039 Nektar::LocalRegions::SegExpSharedPtr Exp =
1040 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1041 basisKeyDir1, segGeom);
1042
1043 int nelmts = 10;
1044
1045 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1046 for (int i = 0; i < nelmts; ++i)
1047 {
1048 CollExp.push_back(Exp);
1049 }
1050
1051 LibUtilities::SessionReaderSharedPtr dummySession;
1052 Collections::CollectionOptimisation colOpt(dummySession, 1,
1053 Collections::eIterPerExp);
1054 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1055 impTypes[Collections::ePhysDeriv] = Collections::eMatrixFree;
1056 Collections::Collection c(CollExp, impTypes);
1057 c.Initialise(Collections::ePhysDeriv);
1058
1059 const int nq = Exp->GetTotPoints();
1060 Array<OneD, NekDouble> xc(nq), yc(nq);
1061 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1;
1062 Array<OneD, NekDouble> diff1(2 * nelmts * nq);
1063 Array<OneD, NekDouble> diff2(2 * nelmts * nq);
1064
1065 Exp->GetCoords(xc, yc);
1066
1067 for (int i = 0; i < nq; ++i)
1068 {
1069 phys[i] = sin(xc[i]) * cos(yc[i]);
1070 }
1071 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nelmts * nq);
1072
1073 for (int i = 1; i < nelmts; ++i)
1074 {
1075 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1076 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq,
1077 tmp1 = diff1 + (nelmts + i) * nq);
1078 }
1079
1080 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1081 tmp = diff2 + nelmts * nq);
1082
1083 double epsilon = 1.0e-8;
1084 for (int i = 0; i < diff1.size(); ++i)
1085 {
1086 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1087 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1088 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1089 }
1090}
1091
1092BOOST_AUTO_TEST_CASE(TestSegPhysDeriv_MatrixFree_UniformP_MultiElmt_3D)
1093{
1094 SpatialDomains::PointGeomSharedPtr v0(
1095 new SpatialDomains::PointGeom(3u, 0u, -1.5, 0.0, 0.0));
1096 SpatialDomains::PointGeomSharedPtr v1(
1097 new SpatialDomains::PointGeom(3u, 1u, 1.0, 1.0, 1.0));
1098
1099 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1, 3u);
1100
1101 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1102 Nektar::LibUtilities::eGaussLobattoLegendre;
1103 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(5,
1104 segPointsTypeDir1);
1105 Nektar::LibUtilities::BasisType basisTypeDir1 =
1106 Nektar::LibUtilities::eModified_A;
1107 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1108 segPointsKeyDir1);
1109
1110 Nektar::LocalRegions::SegExpSharedPtr Exp =
1111 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1112 basisKeyDir1, segGeom);
1113
1114 int nelmts = 10;
1115
1116 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1117 for (int i = 0; i < nelmts; ++i)
1118 {
1119 CollExp.push_back(Exp);
1120 }
1121
1122 LibUtilities::SessionReaderSharedPtr dummySession;
1123 Collections::CollectionOptimisation colOpt(dummySession, 1,
1124 Collections::eIterPerExp);
1125 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1126 impTypes[Collections::ePhysDeriv] = Collections::eMatrixFree;
1127 Collections::Collection c(CollExp, impTypes);
1128 c.Initialise(Collections::ePhysDeriv);
1129
1130 const int nq = Exp->GetTotPoints();
1131 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1132 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1133 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1134 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1135
1136 Exp->GetCoords(xc, yc, zc);
1137
1138 for (int i = 0; i < nq; ++i)
1139 {
1140 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1141 }
1142 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nelmts * nq,
1143 tmp1 = diff1 + 2 * nelmts * nq);
1144
1145 for (int i = 1; i < nelmts; ++i)
1146 {
1147 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1148 Exp->PhysDeriv(phys + i * nq, tmp = diff1 + i * nq,
1149 tmp1 = diff1 + (nelmts + i) * nq,
1150 tmp2 = diff1 + (2 * nelmts + i) * nq);
1151 }
1152
1153 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1154 tmp = diff2 + nelmts * nq, tmp1 = diff2 + 2 * nelmts * nq);
1155
1156 double epsilon = 1.0e-8;
1157 for (int i = 0; i < diff1.size(); ++i)
1158 {
1159 diff1[i] = (fabs(diff1[i]) < 1e-14) ? 0.0 : diff1[i];
1160 diff2[i] = (fabs(diff2[i]) < 1e-14) ? 0.0 : diff2[i];
1161 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1162 }
1163}
1164
1165BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_IterPerExp_UniformP)
1166{
1167 SpatialDomains::PointGeomSharedPtr v0(
1168 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1169 SpatialDomains::PointGeomSharedPtr v1(
1170 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1171
1172 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1173
1174 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1175 Nektar::LibUtilities::eGaussLobattoLegendre;
1176 Nektar::LibUtilities::BasisType basisTypeDir1 =
1177 Nektar::LibUtilities::eModified_A;
1178 unsigned int numSegPoints = 6;
1179 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1180 segPointsTypeDir1);
1181 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1182 segPointsKeyDir1);
1183
1184 Nektar::LocalRegions::SegExpSharedPtr Exp =
1185 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1186 basisKeyDir1, segGeom);
1187
1188 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1189 CollExp.push_back(Exp);
1190
1191 LibUtilities::SessionReaderSharedPtr dummySession;
1192 Collections::CollectionOptimisation colOpt(dummySession, 1,
1193 Collections::eIterPerExp);
1194 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1195 Collections::Collection c(CollExp, impTypes);
1196 c.Initialise(Collections::eIProductWRTDerivBase);
1197
1198 const int nq = Exp->GetTotPoints();
1199 const int nm = Exp->GetNcoeffs();
1200 Array<OneD, NekDouble> phys1(nq);
1201 Array<OneD, NekDouble> coeffs1(nm);
1202 Array<OneD, NekDouble> coeffs2(nm);
1203
1204 Array<OneD, NekDouble> xc(nq);
1205
1206 Exp->GetCoords(xc);
1207
1208 for (int i = 0; i < nq; ++i)
1209 {
1210 phys1[i] = sin(xc[i]);
1211 }
1212
1213 // Standard routines
1214 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1215
1216 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1217
1218 double epsilon = 1.0e-8;
1219 for (int i = 0; i < coeffs1.size(); ++i)
1220 {
1221 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1222 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1223 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1224 }
1225}
1226
1227BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_IterPerExp_UniformP_MultiElmt)
1228{
1229 SpatialDomains::PointGeomSharedPtr v0(
1230 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1231 SpatialDomains::PointGeomSharedPtr v1(
1232 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1233
1234 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1235
1236 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1237 Nektar::LibUtilities::eGaussLobattoLegendre;
1238 Nektar::LibUtilities::BasisType basisTypeDir1 =
1239 Nektar::LibUtilities::eModified_A;
1240 unsigned int numSegPoints = 6;
1241 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1242 segPointsTypeDir1);
1243 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1244 segPointsKeyDir1);
1245
1246 Nektar::LocalRegions::SegExpSharedPtr Exp =
1247 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1248 basisKeyDir1, segGeom);
1249
1250 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1251
1252 int nelmts = 10;
1253 for (int i = 0; i < nelmts; ++i)
1254 {
1255 CollExp.push_back(Exp);
1256 }
1257
1258 LibUtilities::SessionReaderSharedPtr dummySession;
1259 Collections::CollectionOptimisation colOpt(dummySession, 1,
1260 Collections::eIterPerExp);
1261 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1262 Collections::Collection c(CollExp, impTypes);
1263 c.Initialise(Collections::eIProductWRTDerivBase);
1264
1265 const int nq = Exp->GetTotPoints();
1266 const int nm = Exp->GetNcoeffs();
1267 Array<OneD, NekDouble> xc(nq), tmp, tmp1;
1268 Array<OneD, NekDouble> phys1(nelmts * nq);
1269 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1270 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1271
1272 Exp->GetCoords(xc);
1273
1274 for (int i = 0; i < nq; ++i)
1275 {
1276 phys1[i] = sin(xc[i]);
1277 }
1278 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1279
1280 for (int i = 1; i < nelmts; ++i)
1281 {
1282 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1283
1284 // Standard routines
1285 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1286 }
1287
1288 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1289
1290 double epsilon = 1.0e-8;
1291 for (int i = 0; i < coeffs1.size(); ++i)
1292 {
1293 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1294 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1295 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1296 }
1297}
1298
1299BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_StdMat_UniformP)
1300{
1301 SpatialDomains::PointGeomSharedPtr v0(
1302 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1303 SpatialDomains::PointGeomSharedPtr v1(
1304 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1305
1306 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1307
1308 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1309 Nektar::LibUtilities::eGaussLobattoLegendre;
1310 Nektar::LibUtilities::BasisType basisTypeDir1 =
1311 Nektar::LibUtilities::eModified_A;
1312 unsigned int numSegPoints = 6;
1313 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1314 segPointsTypeDir1);
1315 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1316 segPointsKeyDir1);
1317
1318 Nektar::LocalRegions::SegExpSharedPtr Exp =
1319 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1320 basisKeyDir1, segGeom);
1321
1322 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1323 CollExp.push_back(Exp);
1324
1325 LibUtilities::SessionReaderSharedPtr dummySession;
1326 Collections::CollectionOptimisation colOpt(dummySession, 1,
1327 Collections::eStdMat);
1328 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1329 Collections::Collection c(CollExp, impTypes);
1330 c.Initialise(Collections::eIProductWRTDerivBase);
1331
1332 const int nq = Exp->GetTotPoints();
1333 const int nm = Exp->GetNcoeffs();
1334 Array<OneD, NekDouble> phys1(nq);
1335 Array<OneD, NekDouble> coeffs1(nm);
1336 Array<OneD, NekDouble> coeffs2(nm);
1337
1338 Array<OneD, NekDouble> xc(nq);
1339
1340 Exp->GetCoords(xc);
1341
1342 for (int i = 0; i < nq; ++i)
1343 {
1344 phys1[i] = sin(xc[i]);
1345 }
1346
1347 // Standard routines
1348 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1349
1350 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1351
1352 double epsilon = 1.0e-8;
1353 for (int i = 0; i < coeffs1.size(); ++i)
1354 {
1355 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1356 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1357 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1358 }
1359}
1360
1361BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_StdMat_UniformP_MultiElmt)
1362{
1363 SpatialDomains::PointGeomSharedPtr v0(
1364 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1365 SpatialDomains::PointGeomSharedPtr v1(
1366 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1367
1368 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1369
1370 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1371 Nektar::LibUtilities::eGaussLobattoLegendre;
1372 Nektar::LibUtilities::BasisType basisTypeDir1 =
1373 Nektar::LibUtilities::eModified_A;
1374 unsigned int numSegPoints = 6;
1375 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1376 segPointsTypeDir1);
1377 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1378 segPointsKeyDir1);
1379
1380 Nektar::LocalRegions::SegExpSharedPtr Exp =
1381 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1382 basisKeyDir1, segGeom);
1383
1384 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1385
1386 int nelmts = 10;
1387 for (int i = 0; i < nelmts; ++i)
1388 {
1389 CollExp.push_back(Exp);
1390 }
1391
1392 LibUtilities::SessionReaderSharedPtr dummySession;
1393 Collections::CollectionOptimisation colOpt(dummySession, 1,
1394 Collections::eStdMat);
1395 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1396 Collections::Collection c(CollExp, impTypes);
1397 c.Initialise(Collections::eIProductWRTDerivBase);
1398
1399 const int nq = Exp->GetTotPoints();
1400 const int nm = Exp->GetNcoeffs();
1401 Array<OneD, NekDouble> xc(nq), tmp, tmp1;
1402 Array<OneD, NekDouble> phys1(nelmts * nq);
1403 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1404 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1405
1406 Exp->GetCoords(xc);
1407
1408 for (int i = 0; i < nq; ++i)
1409 {
1410 phys1[i] = sin(xc[i]);
1411 }
1412 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1413
1414 for (int i = 1; i < nelmts; ++i)
1415 {
1416 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1417
1418 // Standard routines
1419 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1420 }
1421
1422 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1423
1424 double epsilon = 1.0e-8;
1425 for (int i = 0; i < coeffs1.size(); ++i)
1426 {
1427 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1428 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1429 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1430 }
1431}
1432
1433BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_SumFac_UniformP)
1434{
1435 SpatialDomains::PointGeomSharedPtr v0(
1436 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1437 SpatialDomains::PointGeomSharedPtr v1(
1438 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1439
1440 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1441
1442 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1443 Nektar::LibUtilities::eGaussLobattoLegendre;
1444 Nektar::LibUtilities::BasisType basisTypeDir1 =
1445 Nektar::LibUtilities::eModified_A;
1446 unsigned int numSegPoints = 6;
1447 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1448 segPointsTypeDir1);
1449 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1450 segPointsKeyDir1);
1451
1452 Nektar::LocalRegions::SegExpSharedPtr Exp =
1453 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1454 basisKeyDir1, segGeom);
1455
1456 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1457 CollExp.push_back(Exp);
1458
1459 LibUtilities::SessionReaderSharedPtr dummySession;
1460 Collections::CollectionOptimisation colOpt(dummySession, 1,
1461 Collections::eSumFac);
1462 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1463 Collections::Collection c(CollExp, impTypes);
1464 c.Initialise(Collections::eIProductWRTDerivBase);
1465
1466 const int nq = Exp->GetTotPoints();
1467 const int nm = Exp->GetNcoeffs();
1468 Array<OneD, NekDouble> phys1(nq);
1469 Array<OneD, NekDouble> coeffs1(nm);
1470 Array<OneD, NekDouble> coeffs2(nm);
1471
1472 Array<OneD, NekDouble> xc(nq);
1473
1474 Exp->GetCoords(xc);
1475
1476 for (int i = 0; i < nq; ++i)
1477 {
1478 phys1[i] = sin(xc[i]);
1479 }
1480
1481 // Standard routines
1482 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1483
1484 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1485
1486 double epsilon = 1.0e-8;
1487 for (int i = 0; i < coeffs1.size(); ++i)
1488 {
1489 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1490 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1491 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1492 }
1493}
1494
1495BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_SumFac_UniformP_MultiElmt)
1496{
1497 SpatialDomains::PointGeomSharedPtr v0(
1498 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1499 SpatialDomains::PointGeomSharedPtr v1(
1500 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1501
1502 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1503
1504 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1505 Nektar::LibUtilities::eGaussLobattoLegendre;
1506 Nektar::LibUtilities::BasisType basisTypeDir1 =
1507 Nektar::LibUtilities::eModified_A;
1508 unsigned int numSegPoints = 6;
1509 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1510 segPointsTypeDir1);
1511 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1512 segPointsKeyDir1);
1513
1514 Nektar::LocalRegions::SegExpSharedPtr Exp =
1515 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1516 basisKeyDir1, segGeom);
1517
1518 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1519
1520 int nelmts = 10;
1521 for (int i = 0; i < nelmts; ++i)
1522 {
1523 CollExp.push_back(Exp);
1524 }
1525
1526 LibUtilities::SessionReaderSharedPtr dummySession;
1527 Collections::CollectionOptimisation colOpt(dummySession, 1,
1528 Collections::eSumFac);
1529 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1530 Collections::Collection c(CollExp, impTypes);
1531 c.Initialise(Collections::eIProductWRTDerivBase);
1532
1533 const int nq = Exp->GetTotPoints();
1534 const int nm = Exp->GetNcoeffs();
1535 Array<OneD, NekDouble> xc(nq), tmp, tmp1;
1536 Array<OneD, NekDouble> phys1(nelmts * nq);
1537 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1538 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1539
1540 Exp->GetCoords(xc);
1541
1542 for (int i = 0; i < nq; ++i)
1543 {
1544 phys1[i] = sin(xc[i]);
1545 }
1546 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1547
1548 for (int i = 1; i < nelmts; ++i)
1549 {
1550 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1551
1552 // Standard routines
1553 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1554 }
1555
1556 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1557
1558 double epsilon = 1.0e-8;
1559 for (int i = 0; i < coeffs1.size(); ++i)
1560 {
1561 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1562 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1563 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1564 }
1565}
1566
1567BOOST_AUTO_TEST_CASE(
1568 TestSegIProductWRTDerivBase_SumFac_UniformP_MultiElmt_CoordimTwo)
1569{
1570 SpatialDomains::PointGeomSharedPtr v0(
1571 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1572 SpatialDomains::PointGeomSharedPtr v1(
1573 new SpatialDomains::PointGeom(1u, 1u, 1.0, 1.0, 0.0));
1574
1575 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1, 2);
1576
1577 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1578 Nektar::LibUtilities::eGaussLobattoLegendre;
1579 Nektar::LibUtilities::BasisType basisTypeDir1 =
1580 Nektar::LibUtilities::eModified_A;
1581 unsigned int numSegPoints = 6;
1582 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1583 segPointsTypeDir1);
1584 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1585 segPointsKeyDir1);
1586
1587 Nektar::LocalRegions::SegExpSharedPtr Exp =
1588 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1589 basisKeyDir1, segGeom);
1590
1591 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1592
1593 int nelmts = 10;
1594 for (int i = 0; i < nelmts; ++i)
1595 {
1596 CollExp.push_back(Exp);
1597 }
1598
1599 LibUtilities::SessionReaderSharedPtr dummySession;
1600 Collections::CollectionOptimisation colOpt(dummySession, 1,
1601 Collections::eSumFac);
1602 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1603 Collections::Collection c(CollExp, impTypes);
1604 c.Initialise(Collections::eIProductWRTDerivBase);
1605
1606 const int nq = Exp->GetTotPoints();
1607 const int nm = Exp->GetNcoeffs();
1608 Array<OneD, NekDouble> xc(nq), yc(nq), tmp, tmp1;
1609 Array<OneD, NekDouble> phys1(nelmts * nq);
1610 Array<OneD, NekDouble> phys2(nelmts * nq);
1611 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1612 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1613
1614 Exp->GetCoords(xc, yc);
1615
1616 for (int i = 0; i < nq; ++i)
1617 {
1618 phys1[i] = sin(xc[i]);
1619 phys2[i] = cos(yc[i]);
1620 }
1621 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1622 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
1623
1624 for (int i = 1; i < nelmts; ++i)
1625 {
1626 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1627 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
1628
1629 // Standard routines
1630 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1631 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
1632 }
1633 Vmath::Vadd(nelmts * nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
1634
1635 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, coeffs2);
1636
1637 double epsilon = 1.0e-8;
1638 for (int i = 0; i < coeffs1.size(); ++i)
1639 {
1640 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1641 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1642 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1643 }
1644}
1645
1646BOOST_AUTO_TEST_CASE(TestSegIProductWRTDerivBase_MatrixFree_UniformP_MultiElmt)
1647{
1648 SpatialDomains::PointGeomSharedPtr v0(
1649 new SpatialDomains::PointGeom(1u, 0u, -1.0, 0.0, 0.0));
1650 SpatialDomains::PointGeomSharedPtr v1(
1651 new SpatialDomains::PointGeom(1u, 1u, 1.0, 0.0, 0.0));
1652
1653 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1);
1654
1655 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1656 Nektar::LibUtilities::eGaussLobattoLegendre;
1657 Nektar::LibUtilities::BasisType basisTypeDir1 =
1658 Nektar::LibUtilities::eModified_A;
1659 unsigned int numSegPoints = 6;
1660 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1661 segPointsTypeDir1);
1662 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1663 segPointsKeyDir1);
1664
1665 Nektar::LocalRegions::SegExpSharedPtr Exp =
1666 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1667 basisKeyDir1, segGeom);
1668
1669 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1670
1671 int nelmts = 10;
1672 for (int i = 0; i < nelmts; ++i)
1673 {
1674 CollExp.push_back(Exp);
1675 }
1676
1677 LibUtilities::SessionReaderSharedPtr dummySession;
1678 Collections::CollectionOptimisation colOpt(dummySession, 1,
1679 Collections::eIterPerExp);
1680 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1681 impTypes[Collections::eIProductWRTDerivBase] = Collections::eMatrixFree;
1682 Collections::Collection c(CollExp, impTypes);
1683 c.Initialise(Collections::eIProductWRTDerivBase);
1684
1685 const int nq = Exp->GetTotPoints();
1686 const int nm = Exp->GetNcoeffs();
1687 Array<OneD, NekDouble> xc(nq), tmp, tmp1;
1688 Array<OneD, NekDouble> phys1(nelmts * nq);
1689 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1690 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1691
1692 Exp->GetCoords(xc);
1693
1694 for (int i = 0; i < nq; ++i)
1695 {
1696 phys1[i] = sin(xc[i]);
1697 }
1698 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1699
1700 for (int i = 1; i < nelmts; ++i)
1701 {
1702 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1703
1704 // Standard routines
1705 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1706 }
1707
1708 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, coeffs2);
1709
1710 double epsilon = 1.0e-8;
1711 for (int i = 0; i < coeffs1.size(); ++i)
1712 {
1713 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1714 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1715 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1716 }
1717}
1718
1719BOOST_AUTO_TEST_CASE(
1720 TestSegIProductWRTDerivBase_MatrixFree_UniformP_MultiElmt_CoordimTwo)
1721{
1722 SpatialDomains::PointGeomSharedPtr v0(
1723 new SpatialDomains::PointGeom(2u, 0u, -1.0, 0.0, 0.0));
1724 SpatialDomains::PointGeomSharedPtr v1(
1725 new SpatialDomains::PointGeom(2u, 1u, 1.0, 1.0, 0.0));
1726
1727 SpatialDomains::SegGeomSharedPtr segGeom = CreateSegGeom(0, v0, v1, 2);
1728
1729 Nektar::LibUtilities::PointsType segPointsTypeDir1 =
1730 Nektar::LibUtilities::eGaussLobattoLegendre;
1731 Nektar::LibUtilities::BasisType basisTypeDir1 =
1732 Nektar::LibUtilities::eModified_A;
1733 unsigned int numSegPoints = 6;
1734 const Nektar::LibUtilities::PointsKey segPointsKeyDir1(numSegPoints,
1735 segPointsTypeDir1);
1736 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1737 segPointsKeyDir1);
1738
1739 Nektar::LocalRegions::SegExpSharedPtr Exp =
1740 MemoryManager<Nektar::LocalRegions::SegExp>::AllocateSharedPtr(
1741 basisKeyDir1, segGeom);
1742
1743 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1744
1745 int nelmts = 10;
1746 for (int i = 0; i < nelmts; ++i)
1747 {
1748 CollExp.push_back(Exp);
1749 }
1750
1751 LibUtilities::SessionReaderSharedPtr dummySession;
1752 Collections::CollectionOptimisation colOpt(dummySession, 1,
1753 Collections::eIterPerExp);
1754 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1755 impTypes[Collections::eIProductWRTDerivBase] = Collections::eMatrixFree;
1756 Collections::Collection c(CollExp, impTypes);
1757 c.Initialise(Collections::eIProductWRTDerivBase);
1758
1759 const int nq = Exp->GetTotPoints();
1760 const int nm = Exp->GetNcoeffs();
1761 Array<OneD, NekDouble> xc(nq), yc(nq), tmp, tmp1;
1762 Array<OneD, NekDouble> phys1(nelmts * nq);
1763 Array<OneD, NekDouble> phys2(nelmts * nq);
1764 Array<OneD, NekDouble> coeffs1(nelmts * nm);
1765 Array<OneD, NekDouble> coeffs2(nelmts * nm);
1766
1767 Exp->GetCoords(xc, yc);
1768
1769 for (int i = 0; i < nq; ++i)
1770 {
1771 phys1[i] = sin(xc[i]);
1772 phys2[i] = cos(yc[i]);
1773 }
1774 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
1775 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
1776
1777 for (int i = 1; i < nelmts; ++i)
1778 {
1779 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
1780 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
1781
1782 // Standard routines
1783 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
1784 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
1785 }
1786 Vmath::Vadd(nelmts * nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
1787
1788 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, coeffs2);
1789
1790 double epsilon = 1.0e-8;
1791 for (int i = 0; i < coeffs1.size(); ++i)
1792 {
1793 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1794 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1795 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1796 }
1797}
1798} // namespace Nektar::SegCollectionTests
Nektar::Collections::Collection::Initialise
COLLECTIONS_EXPORT void Initialise(const OperatorType opType, StdRegions::FactorMap factors=StdRegions::NullFactorMap)
Definition: Collection.cpp:61
Nektar::Collections::Collection::ApplyOperator
void ApplyOperator(const OperatorType &op, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &output)
Definition: Collection.h:124
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::Collections::OperatorImpMap
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:126
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::SegExpSharedPtr
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:248
Nektar::SegCollectionTests::CreateSegGeom
SpatialDomains::SegGeomSharedPtr CreateSegGeom(unsigned int id, SpatialDomains::PointGeomSharedPtr v0, SpatialDomains::PointGeomSharedPtr v1, int coordim=1)
Definition: TestSegCollection.cpp:44
Nektar::SegCollectionTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestSegBwdTrans_StdMat_UniformP)
Definition: TestSegCollection.cpp:54
Nektar::SpatialDomains::SegGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:59
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
Generated on Wed Jan 31 2024 09:35:57 for Nektar++ by doxygen 1.9.4