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