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