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