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