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