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
43{
44namespace TetCollectionTests
45{
46SpatialDomains::SegGeomSharedPtr CreateSegGeom(
47 unsigned int id, SpatialDomains::PointGeomSharedPtr v0,
48 SpatialDomains::PointGeomSharedPtr v1)
49{
50 SpatialDomains::PointGeomSharedPtr vertices[] = {v0, v1};
51 SpatialDomains::SegGeomSharedPtr result(
52 new SpatialDomains::SegGeom(id, 3, vertices));
53 return result;
54}
55
56SpatialDomains::TetGeomSharedPtr CreateTet(
57 SpatialDomains::PointGeomSharedPtr v0,
58 SpatialDomains::PointGeomSharedPtr v1,
59 SpatialDomains::PointGeomSharedPtr v2,
60 SpatialDomains::PointGeomSharedPtr v3)
61{
62 Nektar::SpatialDomains::SegGeomSharedPtr e0 = CreateSegGeom(0, v0, v1);
63 Nektar::SpatialDomains::SegGeomSharedPtr e1 = CreateSegGeom(1, v1, v2);
64 Nektar::SpatialDomains::SegGeomSharedPtr e2 = CreateSegGeom(2, v2, v0);
65 Nektar::SpatialDomains::SegGeomSharedPtr e3 = CreateSegGeom(3, v0, v3);
66 Nektar::SpatialDomains::SegGeomSharedPtr e4 = CreateSegGeom(4, v1, v3);
67 Nektar::SpatialDomains::SegGeomSharedPtr e5 = CreateSegGeom(5, v2, v3);
68
69 Nektar::SpatialDomains::SegGeomSharedPtr
70 edgesF0[Nektar::SpatialDomains::TriGeom::kNedges] = {e0, e1, e2};
71 Nektar::SpatialDomains::SegGeomSharedPtr
72 edgesF1[Nektar::SpatialDomains::TriGeom::kNedges] = {e0, e3, e4};
73 Nektar::SpatialDomains::SegGeomSharedPtr
74 edgesF2[Nektar::SpatialDomains::TriGeom::kNedges] = {e1, e4, e5};
75 Nektar::SpatialDomains::SegGeomSharedPtr
76 edgesF3[Nektar::SpatialDomains::TriGeom::kNedges] = {e2, e3, e5};
77
78 Nektar::SpatialDomains::TriGeomSharedPtr face0(
79 new SpatialDomains::TriGeom(0, edgesF0));
80 Nektar::SpatialDomains::TriGeomSharedPtr face1(
81 new SpatialDomains::TriGeom(1, edgesF1));
82 Nektar::SpatialDomains::TriGeomSharedPtr face2(
83 new SpatialDomains::TriGeom(2, edgesF2));
84 Nektar::SpatialDomains::TriGeomSharedPtr face3(
85 new SpatialDomains::TriGeom(3, edgesF3));
86
87 Nektar::SpatialDomains::TriGeomSharedPtr tfaces[] = {face0, face1, face2,
88 face3};
89 SpatialDomains::TetGeomSharedPtr tetGeom(
90 new SpatialDomains::TetGeom(0, tfaces));
91 return tetGeom;
92}
93
94BOOST_AUTO_TEST_CASE(TestTetBwdTrans_IterPerExp_UniformP)
95{
96 SpatialDomains::PointGeomSharedPtr v0(
97 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
98 SpatialDomains::PointGeomSharedPtr v1(
99 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
100 SpatialDomains::PointGeomSharedPtr v2(
101 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
102 SpatialDomains::PointGeomSharedPtr v3(
103 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
104
105 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
106
107 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
108 Nektar::LibUtilities::eGaussLobattoLegendre;
109 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
110 triPointsTypeDir1);
111 Nektar::LibUtilities::BasisType basisTypeDir1 =
112 Nektar::LibUtilities::eModified_A;
113 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
114 triPointsKeyDir1);
115
116 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
117 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
118 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
119 triPointsTypeDir2);
120 Nektar::LibUtilities::BasisType basisTypeDir2 =
121 Nektar::LibUtilities::eModified_B;
122 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
123 triPointsKeyDir2);
124
125 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
126 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
127 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
128 triPointsTypeDir3);
129 Nektar::LibUtilities::BasisType basisTypeDir3 =
130 Nektar::LibUtilities::eModified_C;
131 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
132 triPointsKeyDir3);
133
134 Nektar::LocalRegions::TetExpSharedPtr Exp =
135 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
136 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
137
138 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
139 CollExp.push_back(Exp);
140
141 LibUtilities::SessionReaderSharedPtr dummySession;
142 Collections::CollectionOptimisation colOpt(dummySession, 3,
143 Collections::eIterPerExp);
144 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
145 Collections::Collection c(CollExp, impTypes);
146 c.Initialise(Collections::eBwdTrans);
147
148 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
149 Array<OneD, NekDouble> phys1(Exp->GetTotPoints());
150 Array<OneD, NekDouble> phys2(Exp->GetTotPoints());
151
152 Exp->BwdTrans(coeffs, phys1);
153 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
154
155 double epsilon = 1.0e-8;
156 for (int i = 0; i < phys1.size(); ++i)
157 {
158 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
159 }
160}
161
162BOOST_AUTO_TEST_CASE(TestTetBwdTrans_IterPerExp_VariableP_MultiElmt)
163{
164 SpatialDomains::PointGeomSharedPtr v0(
165 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
166 SpatialDomains::PointGeomSharedPtr v1(
167 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
168 SpatialDomains::PointGeomSharedPtr v2(
169 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
170 SpatialDomains::PointGeomSharedPtr v3(
171 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
172
173 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
174
175 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
176 Nektar::LibUtilities::eGaussLobattoLegendre;
177 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
178 triPointsTypeDir1);
179 Nektar::LibUtilities::BasisType basisTypeDir1 =
180 Nektar::LibUtilities::eModified_A;
181 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
182 triPointsKeyDir1);
183
184 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
185 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
186 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
187 triPointsTypeDir2);
188 Nektar::LibUtilities::BasisType basisTypeDir2 =
189 Nektar::LibUtilities::eModified_B;
190 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
191 triPointsKeyDir2);
192
193 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
194 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
195 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
196 triPointsTypeDir3);
197 Nektar::LibUtilities::BasisType basisTypeDir3 =
198 Nektar::LibUtilities::eModified_C;
199 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
200 triPointsKeyDir3);
201
202 Nektar::LocalRegions::TetExpSharedPtr Exp =
203 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
204 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
205
206 int nelmts = 10;
207
208 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
209 for (int i = 0; i < nelmts; ++i)
210 CollExp.push_back(Exp);
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 CollExp.push_back(Exp);
673
674 LibUtilities::SessionReaderSharedPtr dummySession;
675 Collections::CollectionOptimisation colOpt(dummySession, 3,
676 Collections::eSumFac);
677 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
678 Collections::Collection c(CollExp, impTypes);
679 c.Initialise(Collections::eBwdTrans);
680
681 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
682 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
683 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
684
685 for (int i = 0; i < nelmts; ++i)
686 {
687 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
688 tmp = phys1 + i * Exp->GetTotPoints());
689 }
690
691 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
692
693 double epsilon = 1.0e-8;
694 for (int i = 0; i < phys1.size(); ++i)
695 {
696 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
697 }
698}
699
700BOOST_AUTO_TEST_CASE(TestTetBwdTrans_SumFac_MultiElmt_VariableP)
701{
702 SpatialDomains::PointGeomSharedPtr v0(
703 new SpatialDomains::PointGeom(3u, 0u, -1.5, -1.5, -1.5));
704 SpatialDomains::PointGeomSharedPtr v1(
705 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
706 SpatialDomains::PointGeomSharedPtr v2(
707 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
708 SpatialDomains::PointGeomSharedPtr v3(
709 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
710
711 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
712
713 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
714 Nektar::LibUtilities::eGaussLobattoLegendre;
715 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
716 triPointsTypeDir1);
717 Nektar::LibUtilities::BasisType basisTypeDir1 =
718 Nektar::LibUtilities::eModified_A;
719 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
720 triPointsKeyDir1);
721
722 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
723 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
724 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(6,
725 triPointsTypeDir2);
726 Nektar::LibUtilities::BasisType basisTypeDir2 =
727 Nektar::LibUtilities::eModified_B;
728 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
729 triPointsKeyDir2);
730
731 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
732 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
733 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
734 triPointsTypeDir3);
735 Nektar::LibUtilities::BasisType basisTypeDir3 =
736 Nektar::LibUtilities::eModified_C;
737 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
738 triPointsKeyDir3);
739
740 Nektar::LocalRegions::TetExpSharedPtr Exp =
741 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
742 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
743
744 int nelmts = 1;
745
746 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
747 for (int i = 0; i < nelmts; ++i)
748 CollExp.push_back(Exp);
749
750 LibUtilities::SessionReaderSharedPtr dummySession;
751 Collections::CollectionOptimisation colOpt(dummySession, 3,
752 Collections::eSumFac);
753 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
754 Collections::Collection c(CollExp, impTypes);
755 c.Initialise(Collections::eBwdTrans);
756
757 Array<OneD, NekDouble> coeffs(nelmts * Exp->GetNcoeffs(), 1.0), tmp;
758 Array<OneD, NekDouble> phys1(nelmts * Exp->GetTotPoints());
759 Array<OneD, NekDouble> phys2(nelmts * Exp->GetTotPoints());
760
761 for (int i = 0; i < nelmts; ++i)
762 {
763 Exp->BwdTrans(coeffs + i * Exp->GetNcoeffs(),
764 tmp = phys1 + i * Exp->GetTotPoints());
765 }
766
767 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys2);
768
769 double epsilon = 1.0e-8;
770 for (int i = 0; i < phys1.size(); ++i)
771 {
772 BOOST_CHECK_CLOSE(phys1[i], phys2[i], epsilon);
773 }
774}
775
776BOOST_AUTO_TEST_CASE(TestTetBwdTrans_MatrixFree_UniformP)
777{
778 SpatialDomains::PointGeomSharedPtr v0(
779 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
780 SpatialDomains::PointGeomSharedPtr v1(
781 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
782 SpatialDomains::PointGeomSharedPtr v2(
783 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
784 SpatialDomains::PointGeomSharedPtr v3(
785 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
786
787 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
788
789 unsigned int numQuadPoints = 5;
790 unsigned int numModes = 4;
791
792 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
793 Nektar::LibUtilities::eGaussLobattoLegendre;
794 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
795 triPointsTypeDir1);
796 Nektar::LibUtilities::BasisType basisTypeDir1 =
797 Nektar::LibUtilities::eModified_A;
798 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
799 triPointsKeyDir1);
800
801 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
802 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
803 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
804 triPointsTypeDir2);
805 Nektar::LibUtilities::BasisType basisTypeDir2 =
806 Nektar::LibUtilities::eModified_B;
807 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
808 triPointsKeyDir2);
809
810 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
811 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
812 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
813 triPointsTypeDir3);
814 Nektar::LibUtilities::BasisType basisTypeDir3 =
815 Nektar::LibUtilities::eModified_C;
816 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
817 triPointsKeyDir3);
818
819 Nektar::LocalRegions::TetExpSharedPtr Exp =
820 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
821 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
822
823 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
824 CollExp.push_back(Exp);
825
826 LibUtilities::SessionReaderSharedPtr dummySession;
827 Collections::CollectionOptimisation colOpt(dummySession, 2,
828 Collections::eStdMat);
829 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
830 // ... only one op at the time ...
831 impTypes[Collections::eBwdTrans] = Collections::eMatrixFree;
832 Collections::Collection c(CollExp, impTypes);
833 c.Initialise(Collections::eBwdTrans);
834
835 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
836 Array<OneD, NekDouble> physRef(Exp->GetTotPoints());
837 Array<OneD, NekDouble> phys(Exp->GetTotPoints());
838
839 Exp->BwdTrans(coeffs, physRef);
840 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys);
841
842 double epsilon = 1.0e-8;
843 for (int i = 0; i < physRef.size(); ++i)
844 {
845 BOOST_CHECK_CLOSE(physRef[i], phys[i], epsilon);
846 }
847}
848
849BOOST_AUTO_TEST_CASE(TestTetBwdTrans_MatrixFree_UniformP_OverInt)
850{
851 SpatialDomains::PointGeomSharedPtr v0(
852 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
853 SpatialDomains::PointGeomSharedPtr v1(
854 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
855 SpatialDomains::PointGeomSharedPtr v2(
856 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
857 SpatialDomains::PointGeomSharedPtr v3(
858 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
859
860 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
861
862 unsigned int numQuadPoints = 8;
863 unsigned int numModes = 4;
864
865 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
866 Nektar::LibUtilities::eGaussLobattoLegendre;
867 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
868 triPointsTypeDir1);
869 Nektar::LibUtilities::BasisType basisTypeDir1 =
870 Nektar::LibUtilities::eModified_A;
871 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
872 triPointsKeyDir1);
873
874 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
875 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
876 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
877 triPointsTypeDir2);
878 Nektar::LibUtilities::BasisType basisTypeDir2 =
879 Nektar::LibUtilities::eModified_B;
880 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
881 triPointsKeyDir2);
882
883 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
884 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
885 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
886 triPointsTypeDir3);
887 Nektar::LibUtilities::BasisType basisTypeDir3 =
888 Nektar::LibUtilities::eModified_C;
889 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
890 triPointsKeyDir3);
891
892 Nektar::LocalRegions::TetExpSharedPtr Exp =
893 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
894 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
895
896 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
897 CollExp.push_back(Exp);
898
899 LibUtilities::SessionReaderSharedPtr dummySession;
900 Collections::CollectionOptimisation colOpt(dummySession, 2,
901 Collections::eStdMat);
902 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
903 // ... only one op at the time ...
904 impTypes[Collections::eBwdTrans] = Collections::eMatrixFree;
905 Collections::Collection c(CollExp, impTypes);
906 c.Initialise(Collections::eBwdTrans);
907
908 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs(), 1.0), tmp;
909 Array<OneD, NekDouble> physRef(Exp->GetTotPoints());
910 Array<OneD, NekDouble> phys(Exp->GetTotPoints());
911
912 Exp->BwdTrans(coeffs, physRef);
913 c.ApplyOperator(Collections::eBwdTrans, coeffs, phys);
914
915 double epsilon = 1.0e-8;
916 for (int i = 0; i < physRef.size(); ++i)
917 {
918 BOOST_CHECK_CLOSE(physRef[i], phys[i], epsilon);
919 }
920}
921
922BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_StdMat_UniformP)
923{
924 SpatialDomains::PointGeomSharedPtr v0(
925 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
926 SpatialDomains::PointGeomSharedPtr v1(
927 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
928 SpatialDomains::PointGeomSharedPtr v2(
929 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
930 SpatialDomains::PointGeomSharedPtr v3(
931 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
932
933 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
934
935 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
936 Nektar::LibUtilities::eGaussLobattoLegendre;
937 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
938 triPointsTypeDir1);
939 Nektar::LibUtilities::BasisType basisTypeDir1 =
940 Nektar::LibUtilities::eModified_A;
941 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
942 triPointsKeyDir1);
943
944 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
945 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
946 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
947 triPointsTypeDir2);
948 Nektar::LibUtilities::BasisType basisTypeDir2 =
949 Nektar::LibUtilities::eModified_B;
950 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
951 triPointsKeyDir2);
952
953 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
954 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
955 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
956 triPointsTypeDir3);
957 Nektar::LibUtilities::BasisType basisTypeDir3 =
958 Nektar::LibUtilities::eModified_C;
959 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
960 triPointsKeyDir3);
961
962 Nektar::LocalRegions::TetExpSharedPtr Exp =
963 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
964 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
965
966 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
967 CollExp.push_back(Exp);
968
969 LibUtilities::SessionReaderSharedPtr dummySession;
970 Collections::CollectionOptimisation colOpt(dummySession, 3,
971 Collections::eStdMat);
972 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
973 Collections::Collection c(CollExp, impTypes);
974 c.Initialise(Collections::eIProductWRTBase);
975
976 const int nq = Exp->GetTotPoints();
977 Array<OneD, NekDouble> phys(nq);
978 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
979 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
980
981 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
982
983 Exp->GetCoords(xc, yc, zc);
984
985 for (int i = 0; i < nq; ++i)
986 {
987 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
988 }
989
990 Exp->IProductWRTBase(phys, coeffs1);
991 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
992
993 double epsilon = 1.0e-8;
994 for (int i = 0; i < coeffs1.size(); ++i)
995 {
996 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
997 }
998}
999
1000BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_StdMat_VariableP_MultiElmt)
1001{
1002 SpatialDomains::PointGeomSharedPtr v0(
1003 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1004 SpatialDomains::PointGeomSharedPtr v1(
1005 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1006 SpatialDomains::PointGeomSharedPtr v2(
1007 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1008 SpatialDomains::PointGeomSharedPtr v3(
1009 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1010
1011 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1012
1013 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1014 Nektar::LibUtilities::eGaussLobattoLegendre;
1015 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1016 triPointsTypeDir1);
1017 Nektar::LibUtilities::BasisType basisTypeDir1 =
1018 Nektar::LibUtilities::eModified_A;
1019 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1020 triPointsKeyDir1);
1021
1022 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1023 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1024 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
1025 triPointsTypeDir2);
1026 Nektar::LibUtilities::BasisType basisTypeDir2 =
1027 Nektar::LibUtilities::eModified_B;
1028 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1029 triPointsKeyDir2);
1030
1031 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1032 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1033 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
1034 triPointsTypeDir3);
1035 Nektar::LibUtilities::BasisType basisTypeDir3 =
1036 Nektar::LibUtilities::eModified_C;
1037 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1038 triPointsKeyDir3);
1039
1040 Nektar::LocalRegions::TetExpSharedPtr Exp =
1041 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1042 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1043
1044 int nelmts = 10;
1045
1046 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1047 for (int i = 0; i < nelmts; ++i)
1048 {
1049 CollExp.push_back(Exp);
1050 }
1051
1052 LibUtilities::SessionReaderSharedPtr dummySession;
1053 Collections::CollectionOptimisation colOpt(dummySession, 3,
1054 Collections::eStdMat);
1055 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1056 Collections::Collection c(CollExp, impTypes);
1057 c.Initialise(Collections::eIProductWRTBase);
1058
1059 const int nq = Exp->GetTotPoints();
1060 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
1061 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
1062 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
1063
1064 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1065
1066 Exp->GetCoords(xc, yc, zc);
1067
1068 for (int i = 0; i < nq; ++i)
1069 {
1070 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1071 }
1072 Exp->IProductWRTBase(phys, coeffs1);
1073
1074 for (int i = 1; i < nelmts; ++i)
1075 {
1076 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1077 Exp->IProductWRTBase(phys + i * nq,
1078 tmp = coeffs1 + i * Exp->GetNcoeffs());
1079 }
1080 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1081
1082 double epsilon = 1.0e-4;
1083 for (int i = 0; i < coeffs1.size(); ++i)
1084 {
1085 // clamp values below 1e-14 to zero
1086 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1087 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1088 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1089 }
1090}
1091
1092BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_SumFac_UniformP)
1093{
1094 SpatialDomains::PointGeomSharedPtr v0(
1095 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1096 SpatialDomains::PointGeomSharedPtr v1(
1097 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1098 SpatialDomains::PointGeomSharedPtr v2(
1099 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1100 SpatialDomains::PointGeomSharedPtr v3(
1101 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1102
1103 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1104
1105 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1106 Nektar::LibUtilities::eGaussLobattoLegendre;
1107 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1108 triPointsTypeDir1);
1109 Nektar::LibUtilities::BasisType basisTypeDir1 =
1110 Nektar::LibUtilities::eModified_A;
1111 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1112 triPointsKeyDir1);
1113
1114 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1115 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1116 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
1117 triPointsTypeDir2);
1118 Nektar::LibUtilities::BasisType basisTypeDir2 =
1119 Nektar::LibUtilities::eModified_B;
1120 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1121 triPointsKeyDir2);
1122
1123 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1124 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1125 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
1126 triPointsTypeDir3);
1127 Nektar::LibUtilities::BasisType basisTypeDir3 =
1128 Nektar::LibUtilities::eModified_C;
1129 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1130 triPointsKeyDir3);
1131
1132 Nektar::LocalRegions::TetExpSharedPtr Exp =
1133 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1134 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1135
1136 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1137 CollExp.push_back(Exp);
1138
1139 LibUtilities::SessionReaderSharedPtr dummySession;
1140 Collections::CollectionOptimisation colOpt(dummySession, 3,
1141 Collections::eSumFac);
1142 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1143 Collections::Collection c(CollExp, impTypes);
1144 c.Initialise(Collections::eIProductWRTBase);
1145
1146 const int nq = Exp->GetTotPoints();
1147 Array<OneD, NekDouble> phys(nq);
1148 Array<OneD, NekDouble> coeffs1(Exp->GetNcoeffs());
1149 Array<OneD, NekDouble> coeffs2(Exp->GetNcoeffs());
1150
1151 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1152
1153 Exp->GetCoords(xc, yc, zc);
1154
1155 for (int i = 0; i < nq; ++i)
1156 {
1157 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1158 }
1159
1160 Exp->IProductWRTBase(phys, coeffs1);
1161 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1162
1163 double epsilon = 1.0e-8;
1164 for (int i = 0; i < coeffs1.size(); ++i)
1165 {
1166 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1167 }
1168}
1169
1170BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_SumFac_VariableP_MultiElmt)
1171{
1172 SpatialDomains::PointGeomSharedPtr v0(
1173 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1174 SpatialDomains::PointGeomSharedPtr v1(
1175 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1176 SpatialDomains::PointGeomSharedPtr v2(
1177 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1178 SpatialDomains::PointGeomSharedPtr v3(
1179 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1180
1181 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1182
1183 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1184 Nektar::LibUtilities::eGaussLobattoLegendre;
1185 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1186 triPointsTypeDir1);
1187 Nektar::LibUtilities::BasisType basisTypeDir1 =
1188 Nektar::LibUtilities::eModified_A;
1189 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1190 triPointsKeyDir1);
1191
1192 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1193 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1194 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
1195 triPointsTypeDir2);
1196 Nektar::LibUtilities::BasisType basisTypeDir2 =
1197 Nektar::LibUtilities::eModified_B;
1198 // const Nektar::LibUtilities::BasisKey
1199 // basisKeyDir2(basisTypeDir2,6,triPointsKeyDir2);
1200 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1201 triPointsKeyDir2);
1202
1203 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1204 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1205 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
1206 triPointsTypeDir3);
1207 Nektar::LibUtilities::BasisType basisTypeDir3 =
1208 Nektar::LibUtilities::eModified_C;
1209 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1210 triPointsKeyDir3);
1211
1212 Nektar::LocalRegions::TetExpSharedPtr Exp =
1213 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1214 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1215
1216 int nelmts = 10;
1217
1218 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1219 for (int i = 0; i < nelmts; ++i)
1220 {
1221 CollExp.push_back(Exp);
1222 }
1223
1224 LibUtilities::SessionReaderSharedPtr dummySession;
1225 Collections::CollectionOptimisation colOpt(dummySession, 3,
1226 Collections::eSumFac);
1227 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1228 Collections::Collection c(CollExp, impTypes);
1229 c.Initialise(Collections::eIProductWRTBase);
1230
1231 const int nq = Exp->GetTotPoints();
1232 Array<OneD, NekDouble> phys(nelmts * nq), tmp;
1233 Array<OneD, NekDouble> coeffs1(nelmts * Exp->GetNcoeffs());
1234 Array<OneD, NekDouble> coeffs2(nelmts * Exp->GetNcoeffs());
1235
1236 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1237
1238 Exp->GetCoords(xc, yc, zc);
1239
1240 for (int i = 0; i < nq; ++i)
1241 {
1242 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1243 }
1244 Exp->IProductWRTBase(phys, coeffs1);
1245
1246 for (int i = 1; i < nelmts; ++i)
1247 {
1248 Vmath::Vcopy(nq, &phys[0], 1, &phys[i * nq], 1);
1249 Exp->IProductWRTBase(phys + i * nq,
1250 tmp = coeffs1 + i * Exp->GetNcoeffs());
1251 }
1252 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs2);
1253
1254 double epsilon = 1.0e-4;
1255 for (int i = 0; i < coeffs1.size(); ++i)
1256 {
1257 // clamp values below 1e-14 to zero
1258 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
1259 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
1260 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
1261 }
1262}
1263
1264BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_MatrixFree_UniformP_Undeformed)
1265{
1266 SpatialDomains::PointGeomSharedPtr v0(
1267 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1268 SpatialDomains::PointGeomSharedPtr v1(
1269 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1270 SpatialDomains::PointGeomSharedPtr v2(
1271 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1272 SpatialDomains::PointGeomSharedPtr v3(
1273 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1274
1275 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1276
1277 unsigned int numQuadPoints = 5;
1278 unsigned int numModes = 4;
1279
1280 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1281 Nektar::LibUtilities::eGaussLobattoLegendre;
1282 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
1283 triPointsTypeDir1);
1284 Nektar::LibUtilities::BasisType basisTypeDir1 =
1285 Nektar::LibUtilities::eModified_A;
1286 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1287 triPointsKeyDir1);
1288
1289 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1290 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1291 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
1292 triPointsTypeDir2);
1293 Nektar::LibUtilities::BasisType basisTypeDir2 =
1294 Nektar::LibUtilities::eModified_B;
1295 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
1296 triPointsKeyDir2);
1297
1298 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1299 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1300 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
1301 triPointsTypeDir3);
1302 Nektar::LibUtilities::BasisType basisTypeDir3 =
1303 Nektar::LibUtilities::eModified_C;
1304 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
1305 triPointsKeyDir3);
1306
1307 Nektar::LocalRegions::TetExpSharedPtr Exp =
1308 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1309 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1310
1311 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1312 CollExp.push_back(Exp);
1313
1314 LibUtilities::SessionReaderSharedPtr dummySession;
1315 Collections::CollectionOptimisation colOpt(dummySession, 2,
1316 Collections::eStdMat);
1317 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1318 // ... only one op at the time ...
1319 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
1320 Collections::Collection c(CollExp, impTypes);
1321 c.Initialise(Collections::eIProductWRTBase);
1322
1323 const int nq = Exp->GetTotPoints();
1324 Array<OneD, NekDouble> phys(nq);
1325 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1326 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1327
1328 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1329
1330 Exp->GetCoords(xc, yc, zc);
1331
1332 for (int i = 0; i < nq; ++i)
1333 {
1334 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1335 }
1336
1337 Exp->IProductWRTBase(phys, coeffsRef);
1338 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1339
1340 double epsilon = 1.0e-8;
1341 for (int i = 0; i < coeffsRef.size(); ++i)
1342 {
1343 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1344 }
1345}
1346
1347BOOST_AUTO_TEST_CASE(TestTetIProductWRTBase_MatrixFree_UniformP_Deformed)
1348{
1349 SpatialDomains::PointGeomSharedPtr v0(
1350 new SpatialDomains::PointGeom(3u, 0u, -2.0, -3.0, -4.0));
1351 SpatialDomains::PointGeomSharedPtr v1(
1352 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1353 SpatialDomains::PointGeomSharedPtr v2(
1354 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1355 SpatialDomains::PointGeomSharedPtr v3(
1356 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1357
1358 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1359
1360 unsigned int numQuadPoints = 5;
1361 unsigned int numModes = 4;
1362
1363 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1364 Nektar::LibUtilities::eGaussLobattoLegendre;
1365 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
1366 triPointsTypeDir1);
1367 Nektar::LibUtilities::BasisType basisTypeDir1 =
1368 Nektar::LibUtilities::eModified_A;
1369 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1370 triPointsKeyDir1);
1371
1372 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1373 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1374 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
1375 triPointsTypeDir2);
1376 Nektar::LibUtilities::BasisType basisTypeDir2 =
1377 Nektar::LibUtilities::eModified_B;
1378 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
1379 triPointsKeyDir2);
1380
1381 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1382 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1383 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
1384 triPointsTypeDir3);
1385 Nektar::LibUtilities::BasisType basisTypeDir3 =
1386 Nektar::LibUtilities::eModified_C;
1387 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
1388 triPointsKeyDir3);
1389
1390 Nektar::LocalRegions::TetExpSharedPtr Exp =
1391 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1392 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1393
1394 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1395 CollExp.push_back(Exp);
1396
1397 LibUtilities::SessionReaderSharedPtr dummySession;
1398 Collections::CollectionOptimisation colOpt(dummySession, 2,
1399 Collections::eStdMat);
1400 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1401 // ... only one op at the time ...
1402 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
1403 Collections::Collection c(CollExp, impTypes);
1404 c.Initialise(Collections::eIProductWRTBase);
1405
1406 const int nq = Exp->GetTotPoints();
1407 Array<OneD, NekDouble> phys(nq);
1408 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1409 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1410
1411 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1412
1413 Exp->GetCoords(xc, yc, zc);
1414
1415 for (int i = 0; i < nq; ++i)
1416 {
1417 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1418 }
1419
1420 Exp->IProductWRTBase(phys, coeffsRef);
1421 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1422
1423 double epsilon = 1.0e-8;
1424 for (int i = 0; i < coeffsRef.size(); ++i)
1425 {
1426 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1427 }
1428}
1429
1430BOOST_AUTO_TEST_CASE(
1431 TestTetIProductWRTBase_MatrixFree_UniformP_Deformed_OverInt)
1432{
1433 SpatialDomains::PointGeomSharedPtr v0(
1434 new SpatialDomains::PointGeom(3u, 0u, -2.0, -3.0, -4.0));
1435 SpatialDomains::PointGeomSharedPtr v1(
1436 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1437 SpatialDomains::PointGeomSharedPtr v2(
1438 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1439 SpatialDomains::PointGeomSharedPtr v3(
1440 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1441
1442 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1443
1444 unsigned int numQuadPoints = 8;
1445 unsigned int numModes = 4;
1446
1447 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1448 Nektar::LibUtilities::eGaussLobattoLegendre;
1449 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
1450 triPointsTypeDir1);
1451 Nektar::LibUtilities::BasisType basisTypeDir1 =
1452 Nektar::LibUtilities::eModified_A;
1453 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
1454 triPointsKeyDir1);
1455
1456 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1457 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1458 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
1459 triPointsTypeDir2);
1460 Nektar::LibUtilities::BasisType basisTypeDir2 =
1461 Nektar::LibUtilities::eModified_B;
1462 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
1463 triPointsKeyDir2);
1464
1465 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1466 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1467 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
1468 triPointsTypeDir3);
1469 Nektar::LibUtilities::BasisType basisTypeDir3 =
1470 Nektar::LibUtilities::eModified_C;
1471 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
1472 triPointsKeyDir3);
1473
1474 Nektar::LocalRegions::TetExpSharedPtr Exp =
1475 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1476 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1477
1478 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1479 CollExp.push_back(Exp);
1480
1481 LibUtilities::SessionReaderSharedPtr dummySession;
1482 Collections::CollectionOptimisation colOpt(dummySession, 2,
1483 Collections::eStdMat);
1484 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1485 // ... only one op at the time ...
1486 impTypes[Collections::eIProductWRTBase] = Collections::eMatrixFree;
1487 Collections::Collection c(CollExp, impTypes);
1488 c.Initialise(Collections::eIProductWRTBase);
1489
1490 const int nq = Exp->GetTotPoints();
1491 Array<OneD, NekDouble> phys(nq);
1492 Array<OneD, NekDouble> coeffsRef(Exp->GetNcoeffs());
1493 Array<OneD, NekDouble> coeffs(Exp->GetNcoeffs());
1494
1495 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1496
1497 Exp->GetCoords(xc, yc, zc);
1498
1499 for (int i = 0; i < nq; ++i)
1500 {
1501 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1502 }
1503
1504 Exp->IProductWRTBase(phys, coeffsRef);
1505 c.ApplyOperator(Collections::eIProductWRTBase, phys, coeffs);
1506
1507 double epsilon = 1.0e-8;
1508 for (int i = 0; i < coeffsRef.size(); ++i)
1509 {
1510 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
1511 }
1512}
1513
1514BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_IterPerExp_UniformP)
1515{
1516 SpatialDomains::PointGeomSharedPtr v0(
1517 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1518 SpatialDomains::PointGeomSharedPtr v1(
1519 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1520 SpatialDomains::PointGeomSharedPtr v2(
1521 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1522 SpatialDomains::PointGeomSharedPtr v3(
1523 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1524
1525 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1526
1527 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1528 Nektar::LibUtilities::eGaussLobattoLegendre;
1529 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1530 triPointsTypeDir1);
1531 Nektar::LibUtilities::BasisType basisTypeDir1 =
1532 Nektar::LibUtilities::eModified_A;
1533 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1534 triPointsKeyDir1);
1535
1536 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1537 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1538 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
1539 triPointsTypeDir2);
1540 Nektar::LibUtilities::BasisType basisTypeDir2 =
1541 Nektar::LibUtilities::eModified_B;
1542 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1543 triPointsKeyDir2);
1544
1545 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1546 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1547 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
1548 triPointsTypeDir3);
1549 Nektar::LibUtilities::BasisType basisTypeDir3 =
1550 Nektar::LibUtilities::eModified_C;
1551 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1552 triPointsKeyDir3);
1553
1554 Nektar::LocalRegions::TetExpSharedPtr Exp =
1555 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1556 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1557
1558 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1559 CollExp.push_back(Exp);
1560
1561 LibUtilities::SessionReaderSharedPtr dummySession;
1562 Collections::CollectionOptimisation colOpt(dummySession, 3,
1563 Collections::eIterPerExp);
1564 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1565 Collections::Collection c(CollExp, impTypes);
1566 c.Initialise(Collections::ePhysDeriv);
1567
1568 const int nq = Exp->GetTotPoints();
1569 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1570 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
1571 Array<OneD, NekDouble> diff1(3 * nq);
1572 Array<OneD, NekDouble> diff2(3 * nq);
1573
1574 Exp->GetCoords(xc, yc, zc);
1575
1576 for (int i = 0; i < nq; ++i)
1577 {
1578 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1579 }
1580
1581 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
1582 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
1583 tmp2 = diff2 + 2 * nq);
1584
1585 double epsilon = 1.0e-8;
1586 for (int i = 0; i < diff1.size(); ++i)
1587 {
1588 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1589 }
1590}
1591
1592BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_IterPerExp_VariableP_MultiElmt)
1593{
1594 SpatialDomains::PointGeomSharedPtr v0(
1595 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1596 SpatialDomains::PointGeomSharedPtr v1(
1597 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1598 SpatialDomains::PointGeomSharedPtr v2(
1599 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1600 SpatialDomains::PointGeomSharedPtr v3(
1601 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1602
1603 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1604
1605 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1606 Nektar::LibUtilities::eGaussLobattoLegendre;
1607 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1608 triPointsTypeDir1);
1609 Nektar::LibUtilities::BasisType basisTypeDir1 =
1610 Nektar::LibUtilities::eModified_A;
1611 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1612 triPointsKeyDir1);
1613
1614 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1615 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1616 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
1617 triPointsTypeDir2);
1618 Nektar::LibUtilities::BasisType basisTypeDir2 =
1619 Nektar::LibUtilities::eModified_B;
1620 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1621 triPointsKeyDir2);
1622
1623 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1624 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1625 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
1626 triPointsTypeDir3);
1627 Nektar::LibUtilities::BasisType basisTypeDir3 =
1628 Nektar::LibUtilities::eModified_C;
1629 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1630 triPointsKeyDir3);
1631
1632 Nektar::LocalRegions::TetExpSharedPtr Exp =
1633 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1634 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1635
1636 int nelmts = 10;
1637
1638 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1639 for (int i = 0; i < nelmts; ++i)
1640 {
1641 CollExp.push_back(Exp);
1642 }
1643
1644 LibUtilities::SessionReaderSharedPtr dummySession;
1645 Collections::CollectionOptimisation colOpt(dummySession, 3,
1646 Collections::eIterPerExp);
1647 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1648 Collections::Collection c(CollExp, impTypes);
1649 c.Initialise(Collections::ePhysDeriv);
1650
1651 const int nq = Exp->GetTotPoints();
1652 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1653 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1654 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1655 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1656
1657 Exp->GetCoords(xc, yc, zc);
1658
1659 for (int i = 0; i < nq; ++i)
1660 {
1661 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1662 }
1663 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
1664 tmp2 = diff1 + (2 * nelmts) * nq);
1665
1666 for (int i = 1; i < nelmts; ++i)
1667 {
1668 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1669 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
1670 tmp1 = diff1 + (nelmts + i) * nq,
1671 tmp2 = diff1 + (2 * nelmts + i) * nq);
1672 }
1673
1674 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1675 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
1676
1677 double epsilon = 1.0e-8;
1678 for (int i = 0; i < diff1.size(); ++i)
1679 {
1680 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1681 }
1682}
1683
1684BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_StdMat_UniformP)
1685{
1686 SpatialDomains::PointGeomSharedPtr v0(
1687 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1688 SpatialDomains::PointGeomSharedPtr v1(
1689 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1690 SpatialDomains::PointGeomSharedPtr v2(
1691 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1692 SpatialDomains::PointGeomSharedPtr v3(
1693 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1694
1695 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1696
1697 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1698 Nektar::LibUtilities::eGaussLobattoLegendre;
1699 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1700 triPointsTypeDir1);
1701 Nektar::LibUtilities::BasisType basisTypeDir1 =
1702 Nektar::LibUtilities::eModified_A;
1703 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1704 triPointsKeyDir1);
1705
1706 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1707 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1708 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
1709 triPointsTypeDir2);
1710 Nektar::LibUtilities::BasisType basisTypeDir2 =
1711 Nektar::LibUtilities::eModified_B;
1712 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1713 triPointsKeyDir2);
1714
1715 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1716 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1717 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
1718 triPointsTypeDir3);
1719 Nektar::LibUtilities::BasisType basisTypeDir3 =
1720 Nektar::LibUtilities::eModified_C;
1721 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1722 triPointsKeyDir3);
1723
1724 Nektar::LocalRegions::TetExpSharedPtr Exp =
1725 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1726 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1727
1728 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1729 CollExp.push_back(Exp);
1730
1731 LibUtilities::SessionReaderSharedPtr dummySession;
1732 Collections::CollectionOptimisation colOpt(dummySession, 3,
1733 Collections::eStdMat);
1734 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1735 Collections::Collection c(CollExp, impTypes);
1736 c.Initialise(Collections::ePhysDeriv);
1737
1738 const int nq = Exp->GetTotPoints();
1739 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1740 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
1741 Array<OneD, NekDouble> diff1(3 * nq);
1742 Array<OneD, NekDouble> diff2(3 * nq);
1743
1744 Exp->GetCoords(xc, yc, zc);
1745
1746 for (int i = 0; i < nq; ++i)
1747 {
1748 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1749 }
1750
1751 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
1752 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
1753 tmp2 = diff2 + 2 * nq);
1754
1755 double epsilon = 1.0e-8;
1756 for (int i = 0; i < diff1.size(); ++i)
1757 {
1758 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1759 }
1760}
1761
1762BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_StdMat_VariableP_MultiElmt)
1763{
1764 SpatialDomains::PointGeomSharedPtr v0(
1765 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1766 SpatialDomains::PointGeomSharedPtr v1(
1767 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1768 SpatialDomains::PointGeomSharedPtr v2(
1769 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1770 SpatialDomains::PointGeomSharedPtr v3(
1771 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1772
1773 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1774
1775 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1776 Nektar::LibUtilities::eGaussLobattoLegendre;
1777 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1778 triPointsTypeDir1);
1779 Nektar::LibUtilities::BasisType basisTypeDir1 =
1780 Nektar::LibUtilities::eModified_A;
1781 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1782 triPointsKeyDir1);
1783
1784 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1785 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1786 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
1787 triPointsTypeDir2);
1788 Nektar::LibUtilities::BasisType basisTypeDir2 =
1789 Nektar::LibUtilities::eModified_B;
1790 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1791 triPointsKeyDir2);
1792
1793 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1794 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1795 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
1796 triPointsTypeDir3);
1797 Nektar::LibUtilities::BasisType basisTypeDir3 =
1798 Nektar::LibUtilities::eModified_C;
1799 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1800 triPointsKeyDir3);
1801
1802 Nektar::LocalRegions::TetExpSharedPtr Exp =
1803 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1804 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1805
1806 int nelmts = 10;
1807
1808 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1809 for (int i = 0; i < nelmts; ++i)
1810 {
1811 CollExp.push_back(Exp);
1812 }
1813
1814 LibUtilities::SessionReaderSharedPtr dummySession;
1815 Collections::CollectionOptimisation colOpt(dummySession, 3,
1816 Collections::eStdMat);
1817 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1818 Collections::Collection c(CollExp, impTypes);
1819 c.Initialise(Collections::ePhysDeriv);
1820
1821 const int nq = Exp->GetTotPoints();
1822 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1823 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1824 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1825 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1826
1827 Exp->GetCoords(xc, yc, zc);
1828
1829 for (int i = 0; i < nq; ++i)
1830 {
1831 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1832 }
1833 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
1834 tmp2 = diff1 + (2 * nelmts) * nq);
1835
1836 for (int i = 1; i < nelmts; ++i)
1837 {
1838 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
1839 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
1840 tmp1 = diff1 + (nelmts + i) * nq,
1841 tmp2 = diff1 + (2 * nelmts + i) * nq);
1842 }
1843
1844 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
1845 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
1846
1847 double epsilon = 1.0e-8;
1848 for (int i = 0; i < diff1.size(); ++i)
1849 {
1850 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1851 }
1852}
1853
1854BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_SumFac_UniformP)
1855{
1856 SpatialDomains::PointGeomSharedPtr v0(
1857 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1858 SpatialDomains::PointGeomSharedPtr v1(
1859 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1860 SpatialDomains::PointGeomSharedPtr v2(
1861 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1862 SpatialDomains::PointGeomSharedPtr v3(
1863 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1864
1865 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1866
1867 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1868 Nektar::LibUtilities::eGaussLobattoLegendre;
1869 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1870 triPointsTypeDir1);
1871 Nektar::LibUtilities::BasisType basisTypeDir1 =
1872 Nektar::LibUtilities::eModified_A;
1873 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1874 triPointsKeyDir1);
1875
1876 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1877 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1878 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
1879 triPointsTypeDir2);
1880 Nektar::LibUtilities::BasisType basisTypeDir2 =
1881 Nektar::LibUtilities::eModified_B;
1882 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
1883 triPointsKeyDir2);
1884
1885 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1886 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1887 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
1888 triPointsTypeDir3);
1889 Nektar::LibUtilities::BasisType basisTypeDir3 =
1890 Nektar::LibUtilities::eModified_C;
1891 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
1892 triPointsKeyDir3);
1893
1894 Nektar::LocalRegions::TetExpSharedPtr Exp =
1895 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1896 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1897
1898 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1899 CollExp.push_back(Exp);
1900
1901 LibUtilities::SessionReaderSharedPtr dummySession;
1902 Collections::CollectionOptimisation colOpt(dummySession, 3,
1903 Collections::eSumFac);
1904 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1905 Collections::Collection c(CollExp, impTypes);
1906 c.Initialise(Collections::ePhysDeriv);
1907
1908 const int nq = Exp->GetTotPoints();
1909 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1910 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
1911 Array<OneD, NekDouble> diff1(3 * nq);
1912 Array<OneD, NekDouble> diff2(3 * nq);
1913
1914 Exp->GetCoords(xc, yc, zc);
1915
1916 for (int i = 0; i < nq; ++i)
1917 {
1918 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
1919 }
1920
1921 Exp->PhysDeriv(phys, diff1, tmp = diff1 + nq, tmp1 = diff1 + 2 * nq);
1922 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2, tmp = diff2 + nq,
1923 tmp2 = diff2 + 2 * nq);
1924
1925 double epsilon = 1.0e-8;
1926 for (int i = 0; i < diff1.size(); ++i)
1927 {
1928 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
1929 }
1930}
1931
1932BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_SumFac_VariableP_MultiElmt)
1933{
1934 SpatialDomains::PointGeomSharedPtr v0(
1935 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
1936 SpatialDomains::PointGeomSharedPtr v1(
1937 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
1938 SpatialDomains::PointGeomSharedPtr v2(
1939 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
1940 SpatialDomains::PointGeomSharedPtr v3(
1941 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
1942
1943 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
1944
1945 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
1946 Nektar::LibUtilities::eGaussLobattoLegendre;
1947 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
1948 triPointsTypeDir1);
1949 Nektar::LibUtilities::BasisType basisTypeDir1 =
1950 Nektar::LibUtilities::eModified_A;
1951 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
1952 triPointsKeyDir1);
1953
1954 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
1955 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
1956 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(7,
1957 triPointsTypeDir2);
1958 Nektar::LibUtilities::BasisType basisTypeDir2 =
1959 Nektar::LibUtilities::eModified_B;
1960 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 6,
1961 triPointsKeyDir2);
1962
1963 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
1964 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
1965 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
1966 triPointsTypeDir3);
1967 Nektar::LibUtilities::BasisType basisTypeDir3 =
1968 Nektar::LibUtilities::eModified_C;
1969 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
1970 triPointsKeyDir3);
1971
1972 Nektar::LocalRegions::TetExpSharedPtr Exp =
1973 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
1974 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
1975
1976 int nelmts = 10;
1977
1978 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
1979 for (int i = 0; i < nelmts; ++i)
1980 {
1981 CollExp.push_back(Exp);
1982 }
1983
1984 LibUtilities::SessionReaderSharedPtr dummySession;
1985 Collections::CollectionOptimisation colOpt(dummySession, 3,
1986 Collections::eSumFac);
1987 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
1988 Collections::Collection c(CollExp, impTypes);
1989 c.Initialise(Collections::ePhysDeriv);
1990
1991 const int nq = Exp->GetTotPoints();
1992 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
1993 Array<OneD, NekDouble> phys(nelmts * nq), tmp, tmp1, tmp2;
1994 Array<OneD, NekDouble> diff1(3 * nelmts * nq);
1995 Array<OneD, NekDouble> diff2(3 * nelmts * nq);
1996
1997 Exp->GetCoords(xc, yc, zc);
1998
1999 for (int i = 0; i < nq; ++i)
2000 {
2001 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2002 }
2003 Exp->PhysDeriv(phys, diff1, tmp1 = diff1 + (nelmts)*nq,
2004 tmp2 = diff1 + (2 * nelmts) * nq);
2005
2006 for (int i = 1; i < nelmts; ++i)
2007 {
2008 Vmath::Vcopy(nq, phys, 1, tmp = phys + i * nq, 1);
2009 Exp->PhysDeriv(phys, tmp = diff1 + i * nq,
2010 tmp1 = diff1 + (nelmts + i) * nq,
2011 tmp2 = diff1 + (2 * nelmts + i) * nq);
2012 }
2013
2014 c.ApplyOperator(Collections::ePhysDeriv, phys, diff2,
2015 tmp = diff2 + nelmts * nq, tmp2 = diff2 + 2 * nelmts * nq);
2016
2017 double epsilon = 1.0e-8;
2018 for (int i = 0; i < diff1.size(); ++i)
2019 {
2020 BOOST_CHECK_CLOSE(diff1[i], diff2[i], epsilon);
2021 }
2022}
2023
2024BOOST_AUTO_TEST_CASE(TestTetPhysDeriv_MatrixFree_UniformP)
2025{
2026 SpatialDomains::PointGeomSharedPtr v0(
2027 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2028 SpatialDomains::PointGeomSharedPtr v1(
2029 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2030 SpatialDomains::PointGeomSharedPtr v2(
2031 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2032 SpatialDomains::PointGeomSharedPtr v3(
2033 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2034
2035 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2036
2037 unsigned int numQuadPoints = 5;
2038 unsigned int numModes = 4;
2039
2040 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2041 Nektar::LibUtilities::eGaussLobattoLegendre;
2042 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
2043 triPointsTypeDir1);
2044 Nektar::LibUtilities::BasisType basisTypeDir1 =
2045 Nektar::LibUtilities::eModified_A;
2046 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2047 triPointsKeyDir1);
2048
2049 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2050 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2051 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
2052 triPointsTypeDir2);
2053 Nektar::LibUtilities::BasisType basisTypeDir2 =
2054 Nektar::LibUtilities::eModified_B;
2055 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2056 triPointsKeyDir2);
2057
2058 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2059 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2060 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
2061 triPointsTypeDir3);
2062 Nektar::LibUtilities::BasisType basisTypeDir3 =
2063 Nektar::LibUtilities::eModified_C;
2064 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2065 triPointsKeyDir3);
2066
2067 Nektar::LocalRegions::TetExpSharedPtr Exp =
2068 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2069 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2070
2071 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2072 CollExp.push_back(Exp);
2073
2074 LibUtilities::SessionReaderSharedPtr dummySession;
2075 Collections::CollectionOptimisation colOpt(dummySession, 2,
2076 Collections::eStdMat);
2077 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2078 // ... only one op at the time ...
2079 impTypes[Collections::ePhysDeriv] = Collections::eMatrixFree;
2080 Collections::Collection c(CollExp, impTypes);
2081 c.Initialise(Collections::ePhysDeriv);
2082
2083 const int nq = Exp->GetTotPoints();
2084 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2085 Array<OneD, NekDouble> phys(nq), tmp, tmp1, tmp2;
2086 Array<OneD, NekDouble> diffRef(3 * nq);
2087 Array<OneD, NekDouble> diff(3 * nq);
2088
2089 Exp->GetCoords(xc, yc, zc);
2090
2091 for (int i = 0; i < nq; ++i)
2092 {
2093 phys[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2094 }
2095
2096 Exp->PhysDeriv(phys, diffRef, tmp = diffRef + nq, tmp1 = diffRef + 2 * nq);
2097 c.ApplyOperator(Collections::ePhysDeriv, phys, diff, tmp = diff + nq,
2098 tmp2 = diff + 2 * nq);
2099
2100 double epsilon = 1.0e-8;
2101 for (int i = 0; i < diffRef.size(); ++i)
2102 {
2103 BOOST_CHECK_CLOSE(diffRef[i], diff[i], epsilon);
2104 }
2105}
2106
2107BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_IterPerExp_UniformP)
2108{
2109 SpatialDomains::PointGeomSharedPtr v0(
2110 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2111 SpatialDomains::PointGeomSharedPtr v1(
2112 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2113 SpatialDomains::PointGeomSharedPtr v2(
2114 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2115 SpatialDomains::PointGeomSharedPtr v3(
2116 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2117
2118 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2119
2120 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2121 Nektar::LibUtilities::eGaussLobattoLegendre;
2122 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2123 triPointsTypeDir1);
2124 Nektar::LibUtilities::BasisType basisTypeDir1 =
2125 Nektar::LibUtilities::eModified_A;
2126 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2127 triPointsKeyDir1);
2128
2129 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2130 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2131 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
2132 triPointsTypeDir2);
2133 Nektar::LibUtilities::BasisType basisTypeDir2 =
2134 Nektar::LibUtilities::eModified_B;
2135 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2136 triPointsKeyDir2);
2137
2138 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2139 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2140 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
2141 triPointsTypeDir3);
2142 Nektar::LibUtilities::BasisType basisTypeDir3 =
2143 Nektar::LibUtilities::eModified_C;
2144 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2145 triPointsKeyDir3);
2146
2147 Nektar::LocalRegions::TetExpSharedPtr Exp =
2148 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2149 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2150
2151 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2152 CollExp.push_back(Exp);
2153
2154 LibUtilities::SessionReaderSharedPtr dummySession;
2155 Collections::CollectionOptimisation colOpt(dummySession, 3,
2156 Collections::eIterPerExp);
2157 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2158 Collections::Collection c(CollExp, impTypes);
2159 c.Initialise(Collections::eIProductWRTDerivBase);
2160
2161 const int nq = Exp->GetTotPoints();
2162 const int nm = Exp->GetNcoeffs();
2163 Array<OneD, NekDouble> phys1(nq);
2164 Array<OneD, NekDouble> phys2(nq);
2165 Array<OneD, NekDouble> phys3(nq);
2166 Array<OneD, NekDouble> coeffs1(nm);
2167 Array<OneD, NekDouble> coeffs2(nm);
2168
2169 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2170
2171 Exp->GetCoords(xc, yc, zc);
2172
2173 for (int i = 0; i < nq; ++i)
2174 {
2175 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2176 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2177 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2178 }
2179
2180 // Standard routines
2181 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
2182 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
2183 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2184 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
2185 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2186
2187 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2188 coeffs2);
2189
2190 double epsilon = 1.0e-8;
2191 for (int i = 0; i < coeffs1.size(); ++i)
2192 {
2193 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2194 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2195 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2196 }
2197}
2198
2199BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_IterPerExp_VariableP_MultiElmt)
2200{
2201 SpatialDomains::PointGeomSharedPtr v0(
2202 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2203 SpatialDomains::PointGeomSharedPtr v1(
2204 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2205 SpatialDomains::PointGeomSharedPtr v2(
2206 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2207 SpatialDomains::PointGeomSharedPtr v3(
2208 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2209
2210 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2211
2212 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2213 Nektar::LibUtilities::eGaussLobattoLegendre;
2214 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2215 triPointsTypeDir1);
2216 Nektar::LibUtilities::BasisType basisTypeDir1 =
2217 Nektar::LibUtilities::eModified_A;
2218 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2219 triPointsKeyDir1);
2220
2221 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2222 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2223 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(6,
2224 triPointsTypeDir2);
2225 Nektar::LibUtilities::BasisType basisTypeDir2 =
2226 Nektar::LibUtilities::eModified_B;
2227 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2228 triPointsKeyDir2);
2229 // const Nektar::LibUtilities::BasisKey
2230 // basisKeyDir2(basisTypeDir2,5,triPointsKeyDir2);
2231
2232 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2233 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2234 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
2235 triPointsTypeDir3);
2236 Nektar::LibUtilities::BasisType basisTypeDir3 =
2237 Nektar::LibUtilities::eModified_C;
2238 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
2239 triPointsKeyDir3);
2240
2241 Nektar::LocalRegions::TetExpSharedPtr Exp =
2242 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2243 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2244
2245 int nelmts = 1;
2246
2247 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2248 for (int i = 0; i < nelmts; ++i)
2249 {
2250 CollExp.push_back(Exp);
2251 }
2252
2253 LibUtilities::SessionReaderSharedPtr dummySession;
2254 Collections::CollectionOptimisation colOpt(dummySession, 3,
2255 Collections::eStdMat);
2256 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2257 Collections::Collection c(CollExp, impTypes);
2258 c.Initialise(Collections::eIProductWRTDerivBase);
2259
2260 const int nq = Exp->GetTotPoints();
2261 const int nm = Exp->GetNcoeffs();
2262 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2263 Array<OneD, NekDouble> phys2(nelmts * nq);
2264 Array<OneD, NekDouble> phys3(nelmts * nq);
2265 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2266 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2267
2268 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2269
2270 Exp->GetCoords(xc, yc, zc);
2271
2272 for (int i = 0; i < nq; ++i)
2273 {
2274 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2275 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2276 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2277 }
2278 for (int i = 1; i < nelmts; ++i)
2279 {
2280 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2281 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2282 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2283 }
2284
2285 // Standard routines
2286 for (int i = 0; i < nelmts; ++i)
2287 {
2288 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2289 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2290 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2291 tmp = coeffs1 + i * nm, 1);
2292 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2293 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2294 tmp = coeffs1 + i * nm, 1);
2295 }
2296
2297 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2298 coeffs2);
2299
2300 double epsilon = 1.0e-8;
2301 for (int i = 0; i < coeffs1.size(); ++i)
2302 {
2303 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2304 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2305 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2306 }
2307}
2308
2309BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_StdMat_UniformP)
2310{
2311 SpatialDomains::PointGeomSharedPtr v0(
2312 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2313 SpatialDomains::PointGeomSharedPtr v1(
2314 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2315 SpatialDomains::PointGeomSharedPtr v2(
2316 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2317 SpatialDomains::PointGeomSharedPtr v3(
2318 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2319
2320 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2321
2322 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2323 Nektar::LibUtilities::eGaussLobattoLegendre;
2324 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2325 triPointsTypeDir1);
2326 Nektar::LibUtilities::BasisType basisTypeDir1 =
2327 Nektar::LibUtilities::eModified_A;
2328 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2329 triPointsKeyDir1);
2330
2331 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2332 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2333 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
2334 triPointsTypeDir2);
2335 Nektar::LibUtilities::BasisType basisTypeDir2 =
2336 Nektar::LibUtilities::eModified_B;
2337 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2338 triPointsKeyDir2);
2339
2340 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2341 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2342 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
2343 triPointsTypeDir3);
2344 Nektar::LibUtilities::BasisType basisTypeDir3 =
2345 Nektar::LibUtilities::eModified_C;
2346 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2347 triPointsKeyDir3);
2348
2349 Nektar::LocalRegions::TetExpSharedPtr Exp =
2350 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2351 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2352
2353 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2354 CollExp.push_back(Exp);
2355
2356 LibUtilities::SessionReaderSharedPtr dummySession;
2357 Collections::CollectionOptimisation colOpt(dummySession, 3,
2358 Collections::eStdMat);
2359 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2360 Collections::Collection c(CollExp, impTypes);
2361 c.Initialise(Collections::eIProductWRTDerivBase);
2362
2363 const int nq = Exp->GetTotPoints();
2364 const int nm = Exp->GetNcoeffs();
2365 Array<OneD, NekDouble> phys1(nq);
2366 Array<OneD, NekDouble> phys2(nq);
2367 Array<OneD, NekDouble> phys3(nq);
2368 Array<OneD, NekDouble> coeffs1(nm);
2369 Array<OneD, NekDouble> coeffs2(nm);
2370
2371 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2372
2373 Exp->GetCoords(xc, yc, zc);
2374
2375 for (int i = 0; i < nq; ++i)
2376 {
2377 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2378 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2379 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2380 }
2381
2382 // Standard routines
2383 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
2384 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
2385 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2386 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
2387 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2388
2389 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2390 coeffs2);
2391
2392 double epsilon = 1.0e-8;
2393 for (int i = 0; i < coeffs1.size(); ++i)
2394 {
2395 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2396 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2397 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2398 }
2399}
2400
2401BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_StdMat_VariableP_MultiElmt)
2402{
2403 SpatialDomains::PointGeomSharedPtr v0(
2404 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2405 SpatialDomains::PointGeomSharedPtr v1(
2406 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2407 SpatialDomains::PointGeomSharedPtr v2(
2408 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2409 SpatialDomains::PointGeomSharedPtr v3(
2410 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2411
2412 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2413
2414 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2415 Nektar::LibUtilities::eGaussLobattoLegendre;
2416 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2417 triPointsTypeDir1);
2418 Nektar::LibUtilities::BasisType basisTypeDir1 =
2419 Nektar::LibUtilities::eModified_A;
2420 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2421 triPointsKeyDir1);
2422
2423 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2424 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2425 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(6,
2426 triPointsTypeDir2);
2427 Nektar::LibUtilities::BasisType basisTypeDir2 =
2428 Nektar::LibUtilities::eModified_B;
2429 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2430 triPointsKeyDir2);
2431 // const Nektar::LibUtilities::BasisKey
2432 // basisKeyDir2(basisTypeDir2,5,triPointsKeyDir2);
2433
2434 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2435 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2436 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
2437 triPointsTypeDir3);
2438 Nektar::LibUtilities::BasisType basisTypeDir3 =
2439 Nektar::LibUtilities::eModified_C;
2440 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
2441 triPointsKeyDir3);
2442
2443 Nektar::LocalRegions::TetExpSharedPtr Exp =
2444 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2445 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2446
2447 int nelmts = 1;
2448
2449 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2450 for (int i = 0; i < nelmts; ++i)
2451 {
2452 CollExp.push_back(Exp);
2453 }
2454
2455 LibUtilities::SessionReaderSharedPtr dummySession;
2456 Collections::CollectionOptimisation colOpt(dummySession, 3,
2457 Collections::eStdMat);
2458 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2459 Collections::Collection c(CollExp, impTypes);
2460 c.Initialise(Collections::eIProductWRTDerivBase);
2461
2462 const int nq = Exp->GetTotPoints();
2463 const int nm = Exp->GetNcoeffs();
2464 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2465 Array<OneD, NekDouble> phys2(nelmts * nq);
2466 Array<OneD, NekDouble> phys3(nelmts * nq);
2467 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2468 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2469
2470 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2471
2472 Exp->GetCoords(xc, yc, zc);
2473
2474 for (int i = 0; i < nq; ++i)
2475 {
2476 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2477 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2478 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2479 }
2480 for (int i = 1; i < nelmts; ++i)
2481 {
2482 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2483 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2484 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2485 }
2486
2487 // Standard routines
2488 for (int i = 0; i < nelmts; ++i)
2489 {
2490 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2491 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2492 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2493 tmp = coeffs1 + i * nm, 1);
2494 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2495 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2496 tmp = coeffs1 + i * nm, 1);
2497 }
2498
2499 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2500 coeffs2);
2501
2502 double epsilon = 1.0e-8;
2503 for (int i = 0; i < coeffs1.size(); ++i)
2504 {
2505 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2506 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2507 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2508 }
2509}
2510
2511BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_SumFac_UniformP)
2512{
2513 SpatialDomains::PointGeomSharedPtr v0(
2514 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2515 SpatialDomains::PointGeomSharedPtr v1(
2516 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2517 SpatialDomains::PointGeomSharedPtr v2(
2518 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2519 SpatialDomains::PointGeomSharedPtr v3(
2520 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2521
2522 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2523
2524 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2525 Nektar::LibUtilities::eGaussLobattoLegendre;
2526 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2527 triPointsTypeDir1);
2528 Nektar::LibUtilities::BasisType basisTypeDir1 =
2529 Nektar::LibUtilities::eModified_A;
2530 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2531 triPointsKeyDir1);
2532
2533 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2534 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2535 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(4,
2536 triPointsTypeDir2);
2537 Nektar::LibUtilities::BasisType basisTypeDir2 =
2538 Nektar::LibUtilities::eModified_B;
2539 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2540 triPointsKeyDir2);
2541
2542 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2543 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2544 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(4,
2545 triPointsTypeDir3);
2546 Nektar::LibUtilities::BasisType basisTypeDir3 =
2547 Nektar::LibUtilities::eModified_C;
2548 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 4,
2549 triPointsKeyDir3);
2550
2551 Nektar::LocalRegions::TetExpSharedPtr Exp =
2552 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2553 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2554
2555 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2556 CollExp.push_back(Exp);
2557
2558 LibUtilities::SessionReaderSharedPtr dummySession;
2559 Collections::CollectionOptimisation colOpt(dummySession, 3,
2560 Collections::eSumFac);
2561 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2562 Collections::Collection c(CollExp, impTypes);
2563 c.Initialise(Collections::eIProductWRTDerivBase);
2564
2565 const int nq = Exp->GetTotPoints();
2566 const int nm = Exp->GetNcoeffs();
2567 Array<OneD, NekDouble> phys1(nq);
2568 Array<OneD, NekDouble> phys2(nq);
2569 Array<OneD, NekDouble> phys3(nq);
2570 Array<OneD, NekDouble> coeffs1(nm);
2571 Array<OneD, NekDouble> coeffs2(nm);
2572
2573 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2574
2575 Exp->GetCoords(xc, yc, zc);
2576
2577 for (int i = 0; i < nq; ++i)
2578 {
2579 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2580 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2581 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2582 }
2583
2584 // Standard routines
2585 Exp->IProductWRTDerivBase(0, phys1, coeffs1);
2586 Exp->IProductWRTDerivBase(1, phys2, coeffs2);
2587 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2588 Exp->IProductWRTDerivBase(2, phys3, coeffs2);
2589 Vmath::Vadd(nm, coeffs1, 1, coeffs2, 1, coeffs1, 1);
2590
2591 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2592 coeffs2);
2593
2594 double epsilon = 1.0e-7;
2595 for (int i = 0; i < coeffs1.size(); ++i)
2596 {
2597 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2598 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2599 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2600 }
2601}
2602
2603BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_SumFac_VariableP_MultiElmt)
2604{
2605 SpatialDomains::PointGeomSharedPtr v0(
2606 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2607 SpatialDomains::PointGeomSharedPtr v1(
2608 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2609 SpatialDomains::PointGeomSharedPtr v2(
2610 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2611 SpatialDomains::PointGeomSharedPtr v3(
2612 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2613
2614 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2615
2616 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2617 Nektar::LibUtilities::eGaussLobattoLegendre;
2618 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(5,
2619 triPointsTypeDir1);
2620 Nektar::LibUtilities::BasisType basisTypeDir1 =
2621 Nektar::LibUtilities::eModified_A;
2622 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, 4,
2623 triPointsKeyDir1);
2624
2625 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2626 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2627 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(6,
2628 triPointsTypeDir2);
2629 Nektar::LibUtilities::BasisType basisTypeDir2 =
2630 Nektar::LibUtilities::eModified_B;
2631 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, 4,
2632 triPointsKeyDir2);
2633 // const Nektar::LibUtilities::BasisKey
2634 // basisKeyDir2(basisTypeDir2,5,triPointsKeyDir2);
2635
2636 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2637 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2638 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(9,
2639 triPointsTypeDir3);
2640 Nektar::LibUtilities::BasisType basisTypeDir3 =
2641 Nektar::LibUtilities::eModified_C;
2642 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, 8,
2643 triPointsKeyDir3);
2644
2645 Nektar::LocalRegions::TetExpSharedPtr Exp =
2646 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2647 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2648
2649 int nelmts = 1;
2650
2651 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2652 for (int i = 0; i < nelmts; ++i)
2653 {
2654 CollExp.push_back(Exp);
2655 }
2656
2657 LibUtilities::SessionReaderSharedPtr dummySession;
2658 Collections::CollectionOptimisation colOpt(dummySession, 3,
2659 Collections::eSumFac);
2660 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
2661 Collections::Collection c(CollExp, impTypes);
2662 c.Initialise(Collections::eIProductWRTDerivBase);
2663
2664 const int nq = Exp->GetTotPoints();
2665 const int nm = Exp->GetNcoeffs();
2666 Array<OneD, NekDouble> phys1(nelmts * nq), tmp;
2667 Array<OneD, NekDouble> phys2(nelmts * nq);
2668 Array<OneD, NekDouble> phys3(nelmts * nq);
2669 Array<OneD, NekDouble> coeffs1(nelmts * nm);
2670 Array<OneD, NekDouble> coeffs2(nelmts * nm);
2671
2672 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
2673
2674 Exp->GetCoords(xc, yc, zc);
2675
2676 for (int i = 0; i < nq; ++i)
2677 {
2678 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
2679 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
2680 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
2681 }
2682 for (int i = 1; i < nelmts; ++i)
2683 {
2684 Vmath::Vcopy(nq, phys1, 1, tmp = phys1 + i * nq, 1);
2685 Vmath::Vcopy(nq, phys2, 1, tmp = phys2 + i * nq, 1);
2686 Vmath::Vcopy(nq, phys3, 1, tmp = phys3 + i * nq, 1);
2687 }
2688
2689 // Standard routines
2690 for (int i = 0; i < nelmts; ++i)
2691 {
2692 Exp->IProductWRTDerivBase(0, phys1 + i * nq, tmp = coeffs1 + i * nm);
2693 Exp->IProductWRTDerivBase(1, phys2 + i * nq, tmp = coeffs2 + i * nm);
2694 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2695 tmp = coeffs1 + i * nm, 1);
2696 Exp->IProductWRTDerivBase(2, phys3 + i * nq, tmp = coeffs2 + i * nm);
2697 Vmath::Vadd(nm, coeffs1 + i * nm, 1, coeffs2 + i * nm, 1,
2698 tmp = coeffs1 + i * nm, 1);
2699 }
2700
2701 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
2702 coeffs2);
2703
2704 double epsilon = 1.0e-7;
2705 for (int i = 0; i < coeffs1.size(); ++i)
2706 {
2707 coeffs1[i] = (fabs(coeffs1[i]) < 1e-14) ? 0.0 : coeffs1[i];
2708 coeffs2[i] = (fabs(coeffs2[i]) < 1e-14) ? 0.0 : coeffs2[i];
2709 BOOST_CHECK_CLOSE(coeffs1[i], coeffs2[i], epsilon);
2710 }
2711}
2712
2713BOOST_AUTO_TEST_CASE(TestTetmHelmholtz_IterPerExp_UniformP_ConstVarDiff)
2714{
2715 SpatialDomains::PointGeomSharedPtr v0(
2716 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2717 SpatialDomains::PointGeomSharedPtr v1(
2718 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2719 SpatialDomains::PointGeomSharedPtr v2(
2720 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2721 SpatialDomains::PointGeomSharedPtr v3(
2722 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2723
2724 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2725
2726 unsigned int numQuadPoints = 5;
2727 unsigned int numModes = 4;
2728
2729 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2730 Nektar::LibUtilities::eGaussLobattoLegendre;
2731 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
2732 triPointsTypeDir1);
2733 Nektar::LibUtilities::BasisType basisTypeDir1 =
2734 Nektar::LibUtilities::eModified_A;
2735 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2736 triPointsKeyDir1);
2737
2738 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2739 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2740 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
2741 triPointsTypeDir2);
2742 Nektar::LibUtilities::BasisType basisTypeDir2 =
2743 Nektar::LibUtilities::eModified_B;
2744 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2745 triPointsKeyDir2);
2746
2747 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2748 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2749 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
2750 triPointsTypeDir3);
2751 Nektar::LibUtilities::BasisType basisTypeDir3 =
2752 Nektar::LibUtilities::eModified_C;
2753 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2754 triPointsKeyDir3);
2755
2756 Nektar::LocalRegions::TetExpSharedPtr Exp =
2757 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2758 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2759
2760 Nektar::StdRegions::StdTetExpSharedPtr stdExp =
2761 MemoryManager<Nektar::StdRegions::StdTetExp>::AllocateSharedPtr(
2762 basisKeyDir1, basisKeyDir1, basisKeyDir1);
2763
2764 int nelmts = 10;
2765
2766 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2767 for (int i = 0; i < nelmts; ++i)
2768 {
2769 CollExp.push_back(Exp);
2770 }
2771
2772 LibUtilities::SessionReaderSharedPtr dummySession;
2773 Collections::CollectionOptimisation colOpt(dummySession, 2,
2774 Collections::eIterPerExp);
2775 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(stdExp);
2776 Collections::Collection c(CollExp, impTypes);
2777 StdRegions::ConstFactorMap factors;
2778 factors[StdRegions::eFactorLambda] = 0.0;
2779 factors[StdRegions::eFactorCoeffD00] = 1.25;
2780 factors[StdRegions::eFactorCoeffD01] = 0.25;
2781 factors[StdRegions::eFactorCoeffD11] = 1.25;
2782 factors[StdRegions::eFactorCoeffD02] = 0.25;
2783 factors[StdRegions::eFactorCoeffD12] = 0.25;
2784 factors[StdRegions::eFactorCoeffD22] = 1.25;
2785
2786 c.Initialise(Collections::eHelmholtz, factors);
2787
2788 const int nm = Exp->GetNcoeffs();
2789 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
2790 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
2791 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
2792
2793 for (int i = 0; i < nm; ++i)
2794 {
2795 coeffsIn[i] = 1.0;
2796 }
2797
2798 for (int i = 1; i < nelmts; ++i)
2799 {
2800 Vmath::Vcopy(nm, coeffsIn, 1, tmp = coeffsIn + i * nm, 1);
2801 }
2802
2803 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
2804 *Exp, factors);
2805
2806 for (int i = 0; i < nelmts; ++i)
2807 {
2808 // Standard routines
2809 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
2810 }
2811
2812 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
2813
2814 double epsilon = 1.0e-8;
2815 for (int i = 0; i < coeffsRef.size(); ++i)
2816 {
2817 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
2818 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
2819 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
2820 }
2821}
2822
2823BOOST_AUTO_TEST_CASE(TestTetmHelmholtz_MatrixFree_UniformP)
2824{
2825 SpatialDomains::PointGeomSharedPtr v0(
2826 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
2827 SpatialDomains::PointGeomSharedPtr v1(
2828 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2829 SpatialDomains::PointGeomSharedPtr v2(
2830 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2831 SpatialDomains::PointGeomSharedPtr v3(
2832 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2833
2834 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2835
2836 unsigned int numQuadPoints = 5;
2837 unsigned int numModes = 4;
2838
2839 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2840 Nektar::LibUtilities::eGaussLobattoLegendre;
2841 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
2842 triPointsTypeDir1);
2843 Nektar::LibUtilities::BasisType basisTypeDir1 =
2844 Nektar::LibUtilities::eModified_A;
2845 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2846 triPointsKeyDir1);
2847
2848 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2849 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2850 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
2851 triPointsTypeDir2);
2852 Nektar::LibUtilities::BasisType basisTypeDir2 =
2853 Nektar::LibUtilities::eModified_B;
2854 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2855 triPointsKeyDir2);
2856
2857 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2858 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2859 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
2860 triPointsTypeDir3);
2861 Nektar::LibUtilities::BasisType basisTypeDir3 =
2862 Nektar::LibUtilities::eModified_C;
2863 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2864 triPointsKeyDir3);
2865
2866 Nektar::LocalRegions::TetExpSharedPtr Exp =
2867 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2868 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2869
2870 Nektar::StdRegions::StdTetExpSharedPtr stdExp =
2871 MemoryManager<Nektar::StdRegions::StdTetExp>::AllocateSharedPtr(
2872 basisKeyDir1, basisKeyDir1, basisKeyDir1);
2873
2874 int nelmts = 10;
2875
2876 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2877 for (int i = 0; i < nelmts; ++i)
2878 {
2879 CollExp.push_back(Exp);
2880 }
2881
2882 LibUtilities::SessionReaderSharedPtr dummySession;
2883 Collections::CollectionOptimisation colOpt(dummySession, 2,
2884 Collections::eMatrixFree);
2885 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(stdExp);
2886 Collections::Collection c(CollExp, impTypes);
2887 StdRegions::ConstFactorMap factors;
2888 factors[StdRegions::eFactorLambda] = 0.0;
2889
2890 c.Initialise(Collections::eHelmholtz, factors);
2891
2892 const int nm = Exp->GetNcoeffs();
2893 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
2894 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
2895 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
2896
2897 for (int i = 0; i < nm; ++i)
2898 {
2899 coeffsIn[i] = 1.0;
2900 }
2901
2902 for (int i = 1; i < nelmts; ++i)
2903 {
2904 Vmath::Vcopy(nm, coeffsIn, 1, tmp = coeffsIn + i * nm, 1);
2905 }
2906
2907 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
2908 *Exp, factors);
2909
2910 for (int i = 0; i < nelmts; ++i)
2911 {
2912 // Standard routines
2913 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
2914 }
2915
2916 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
2917
2918 double epsilon = 1.0e-8;
2919 for (int i = 0; i < coeffsRef.size(); ++i)
2920 {
2921 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
2922 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
2923 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
2924 }
2925}
2926
2927BOOST_AUTO_TEST_CASE(TestTetmHelmholtz_MatrixFree_Deformed_OverInt)
2928{
2929 SpatialDomains::PointGeomSharedPtr v0(
2930 new SpatialDomains::PointGeom(3u, 0u, -2.0, -3.0, -4.0));
2931 SpatialDomains::PointGeomSharedPtr v1(
2932 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
2933 SpatialDomains::PointGeomSharedPtr v2(
2934 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
2935 SpatialDomains::PointGeomSharedPtr v3(
2936 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
2937
2938 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
2939
2940 unsigned int numQuadPoints = 8;
2941 unsigned int numModes = 4;
2942
2943 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
2944 Nektar::LibUtilities::eGaussLobattoLegendre;
2945 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
2946 triPointsTypeDir1);
2947 Nektar::LibUtilities::BasisType basisTypeDir1 =
2948 Nektar::LibUtilities::eModified_A;
2949 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
2950 triPointsKeyDir1);
2951
2952 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
2953 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
2954 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
2955 triPointsTypeDir2);
2956 Nektar::LibUtilities::BasisType basisTypeDir2 =
2957 Nektar::LibUtilities::eModified_B;
2958 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
2959 triPointsKeyDir2);
2960
2961 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
2962 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
2963 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
2964 triPointsTypeDir3);
2965 Nektar::LibUtilities::BasisType basisTypeDir3 =
2966 Nektar::LibUtilities::eModified_C;
2967 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
2968 triPointsKeyDir3);
2969
2970 Nektar::LocalRegions::TetExpSharedPtr Exp =
2971 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
2972 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
2973
2974 Nektar::StdRegions::StdTetExpSharedPtr stdExp =
2975 MemoryManager<Nektar::StdRegions::StdTetExp>::AllocateSharedPtr(
2976 basisKeyDir1, basisKeyDir1, basisKeyDir1);
2977
2978 int nelmts = 10;
2979
2980 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
2981 for (int i = 0; i < nelmts; ++i)
2982 {
2983 CollExp.push_back(Exp);
2984 }
2985
2986 LibUtilities::SessionReaderSharedPtr dummySession;
2987 Collections::CollectionOptimisation colOpt(dummySession, 2,
2988 Collections::eMatrixFree);
2989 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(stdExp);
2990 Collections::Collection c(CollExp, impTypes);
2991 StdRegions::ConstFactorMap factors;
2992 factors[StdRegions::eFactorLambda] = 0.0;
2993
2994 c.Initialise(Collections::eHelmholtz, factors);
2995
2996 const int nm = Exp->GetNcoeffs();
2997 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
2998 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
2999 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3000
3001 for (int i = 0; i < nm; ++i)
3002 {
3003 coeffsIn[i] = 1.0;
3004 }
3005
3006 for (int i = 1; i < nelmts; ++i)
3007 {
3008 Vmath::Vcopy(nm, coeffsIn, 1, tmp = coeffsIn + i * nm, 1);
3009 }
3010
3011 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3012 *Exp, factors);
3013
3014 for (int i = 0; i < nelmts; ++i)
3015 {
3016 // Standard routines
3017 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3018 }
3019
3020 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3021
3022 double epsilon = 1.0e-8;
3023 for (int i = 0; i < coeffsRef.size(); ++i)
3024 {
3025 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3026 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3027 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3028 }
3029}
3030
3031BOOST_AUTO_TEST_CASE(TestTetIProductWRTDerivBase_MatrixFree_UniformP_Undeformed)
3032{
3033 SpatialDomains::PointGeomSharedPtr v0(
3034 new SpatialDomains::PointGeom(3u, 0u, -2.0, -3.0, -4.0));
3035 SpatialDomains::PointGeomSharedPtr v1(
3036 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3037 SpatialDomains::PointGeomSharedPtr v2(
3038 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
3039 SpatialDomains::PointGeomSharedPtr v3(
3040 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
3041
3042 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
3043
3044 unsigned int numQuadPoints = 5;
3045 unsigned int numModes = 4;
3046
3047 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
3048 Nektar::LibUtilities::eGaussLobattoLegendre;
3049 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
3050 triPointsTypeDir1);
3051 Nektar::LibUtilities::BasisType basisTypeDir1 =
3052 Nektar::LibUtilities::eModified_A;
3053 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3054 triPointsKeyDir1);
3055
3056 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
3057 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
3058 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
3059 triPointsTypeDir2);
3060 Nektar::LibUtilities::BasisType basisTypeDir2 =
3061 Nektar::LibUtilities::eModified_B;
3062 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3063 triPointsKeyDir2);
3064
3065 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
3066 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3067 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
3068 triPointsTypeDir3);
3069 Nektar::LibUtilities::BasisType basisTypeDir3 =
3070 Nektar::LibUtilities::eModified_C;
3071 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3072 triPointsKeyDir3);
3073
3074 Nektar::LocalRegions::TetExpSharedPtr Exp =
3075 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
3076 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
3077
3078 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
3079 CollExp.push_back(Exp);
3080
3081 LibUtilities::SessionReaderSharedPtr dummySession;
3082 Collections::CollectionOptimisation colOpt(dummySession, 2,
3083 Collections::eStdMat);
3084 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(Exp);
3085 // ... only one op at the time ...
3086 impTypes[Collections::eIProductWRTDerivBase] = Collections::eMatrixFree;
3087 Collections::Collection c(CollExp, impTypes);
3088 c.Initialise(Collections::eIProductWRTDerivBase);
3089
3090 const int nq = Exp->GetTotPoints();
3091 const int nm = Exp->GetNcoeffs();
3092 Array<OneD, NekDouble> phys1(nq);
3093 Array<OneD, NekDouble> phys2(nq);
3094 Array<OneD, NekDouble> phys3(nq);
3095 Array<OneD, NekDouble> coeffsRef(nm);
3096 Array<OneD, NekDouble> coeffs(nm);
3097
3098 Array<OneD, NekDouble> xc(nq), yc(nq), zc(nq);
3099
3100 Exp->GetCoords(xc, yc, zc);
3101
3102 for (int i = 0; i < nq; ++i)
3103 {
3104 phys1[i] = sin(xc[i]) * cos(yc[i]) * sin(zc[i]);
3105 phys2[i] = cos(xc[i]) * sin(yc[i]) * cos(zc[i]);
3106 phys3[i] = cos(xc[i]) * sin(yc[i]) * sin(zc[i]);
3107 }
3108
3109 // Standard routines
3110 Exp->IProductWRTDerivBase(0, phys1, coeffsRef);
3111 Exp->IProductWRTDerivBase(1, phys2, coeffs);
3112 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3113 Exp->IProductWRTDerivBase(2, phys3, coeffs);
3114 Vmath::Vadd(nm, coeffsRef, 1, coeffs, 1, coeffsRef, 1);
3115
3116 c.ApplyOperator(Collections::eIProductWRTDerivBase, phys1, phys2, phys3,
3117 coeffs);
3118
3119 double epsilon = 1.0e-7;
3120 for (int i = 0; i < coeffsRef.size(); ++i)
3121 {
3122 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3123 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3124 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3125 }
3126}
3127
3128BOOST_AUTO_TEST_CASE(TestTetmHelmholtz_MatrixFree_UniformP_ConstVarDiff)
3129{
3130 SpatialDomains::PointGeomSharedPtr v0(
3131 new SpatialDomains::PointGeom(3u, 0u, -1.0, -1.0, -1.0));
3132 SpatialDomains::PointGeomSharedPtr v1(
3133 new SpatialDomains::PointGeom(3u, 1u, 1.0, -1.0, -1.0));
3134 SpatialDomains::PointGeomSharedPtr v2(
3135 new SpatialDomains::PointGeom(3u, 2u, -1.0, 1.0, -1.0));
3136 SpatialDomains::PointGeomSharedPtr v3(
3137 new SpatialDomains::PointGeom(3u, 3u, -1.0, -1.0, 1.0));
3138
3139 SpatialDomains::TetGeomSharedPtr tetGeom = CreateTet(v0, v1, v2, v3);
3140
3141 unsigned int numQuadPoints = 5;
3142 unsigned int numModes = 4;
3143
3144 Nektar::LibUtilities::PointsType triPointsTypeDir1 =
3145 Nektar::LibUtilities::eGaussLobattoLegendre;
3146 const Nektar::LibUtilities::PointsKey triPointsKeyDir1(numQuadPoints,
3147 triPointsTypeDir1);
3148 Nektar::LibUtilities::BasisType basisTypeDir1 =
3149 Nektar::LibUtilities::eModified_A;
3150 const Nektar::LibUtilities::BasisKey basisKeyDir1(basisTypeDir1, numModes,
3151 triPointsKeyDir1);
3152
3153 Nektar::LibUtilities::PointsType triPointsTypeDir2 =
3154 Nektar::LibUtilities::eGaussRadauMAlpha1Beta0;
3155 const Nektar::LibUtilities::PointsKey triPointsKeyDir2(numQuadPoints - 1,
3156 triPointsTypeDir2);
3157 Nektar::LibUtilities::BasisType basisTypeDir2 =
3158 Nektar::LibUtilities::eModified_B;
3159 const Nektar::LibUtilities::BasisKey basisKeyDir2(basisTypeDir2, numModes,
3160 triPointsKeyDir2);
3161
3162 Nektar::LibUtilities::PointsType triPointsTypeDir3 =
3163 Nektar::LibUtilities::eGaussRadauMAlpha2Beta0;
3164 const Nektar::LibUtilities::PointsKey triPointsKeyDir3(numQuadPoints - 1,
3165 triPointsTypeDir3);
3166 Nektar::LibUtilities::BasisType basisTypeDir3 =
3167 Nektar::LibUtilities::eModified_C;
3168 const Nektar::LibUtilities::BasisKey basisKeyDir3(basisTypeDir3, numModes,
3169 triPointsKeyDir3);
3170
3171 Nektar::LocalRegions::TetExpSharedPtr Exp =
3172 MemoryManager<Nektar::LocalRegions::TetExp>::AllocateSharedPtr(
3173 basisKeyDir1, basisKeyDir2, basisKeyDir3, tetGeom);
3174
3175 Nektar::StdRegions::StdTetExpSharedPtr stdExp =
3176 MemoryManager<Nektar::StdRegions::StdTetExp>::AllocateSharedPtr(
3177 basisKeyDir1, basisKeyDir1, basisKeyDir1);
3178
3179 int nelmts = 10;
3180
3181 std::vector<StdRegions::StdExpansionSharedPtr> CollExp;
3182 for (int i = 0; i < nelmts; ++i)
3183 {
3184 CollExp.push_back(Exp);
3185 }
3186
3187 LibUtilities::SessionReaderSharedPtr dummySession;
3188 Collections::CollectionOptimisation colOpt(dummySession, 2,
3189 Collections::eMatrixFree);
3190 Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(stdExp);
3191 Collections::Collection c(CollExp, impTypes);
3192 StdRegions::ConstFactorMap factors;
3193 factors[StdRegions::eFactorLambda] = 0.0;
3194 factors[StdRegions::eFactorCoeffD00] = 1.25;
3195 factors[StdRegions::eFactorCoeffD01] = 0.25;
3196 factors[StdRegions::eFactorCoeffD11] = 1.25;
3197 factors[StdRegions::eFactorCoeffD02] = 0.25;
3198 factors[StdRegions::eFactorCoeffD12] = 0.25;
3199 factors[StdRegions::eFactorCoeffD22] = 1.25;
3200
3201 c.Initialise(Collections::eHelmholtz, factors);
3202
3203 const int nm = Exp->GetNcoeffs();
3204 Array<OneD, NekDouble> coeffsIn(nelmts * nm);
3205 Array<OneD, NekDouble> coeffsRef(nelmts * nm);
3206 Array<OneD, NekDouble> coeffs(nelmts * nm), tmp;
3207
3208 for (int i = 0; i < nm; ++i)
3209 {
3210 coeffsIn[i] = 1.0;
3211 }
3212
3213 for (int i = 1; i < nelmts; ++i)
3214 {
3215 Vmath::Vcopy(nm, coeffsIn, 1, tmp = coeffsIn + i * nm, 1);
3216 }
3217
3218 StdRegions::StdMatrixKey mkey(StdRegions::eHelmholtz, Exp->DetShapeType(),
3219 *Exp, factors);
3220
3221 for (int i = 0; i < nelmts; ++i)
3222 {
3223 // Standard routines
3224 Exp->GeneralMatrixOp(coeffsIn + i * nm, tmp = coeffsRef + i * nm, mkey);
3225 }
3226
3227 c.ApplyOperator(Collections::eHelmholtz, coeffsIn, coeffs);
3228
3229 double epsilon = 1.0e-8;
3230 for (int i = 0; i < coeffsRef.size(); ++i)
3231 {
3232 coeffsRef[i] = (std::abs(coeffsRef[i]) < 1e-14) ? 0.0 : coeffsRef[i];
3233 coeffs[i] = (std::abs(coeffs[i]) < 1e-14) ? 0.0 : coeffs[i];
3234 BOOST_CHECK_CLOSE(coeffsRef[i], coeffs[i], epsilon);
3235 }
3236}
3237} // namespace TetCollectionTests
3238} // namespace Nektar
Nektar::Collections::Collection::Initialise
COLLECTIONS_EXPORT void Initialise(const OperatorType opType, StdRegions::FactorMap factors=StdRegions::NullFactorMap)
Definition: Collection.cpp:63
Nektar::Collections::Collection::ApplyOperator
void ApplyOperator(const OperatorType &op, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &output)
Definition: Collection.h:116
Nektar::Collections::CollectionOptimisation::GetOperatorImpMap
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Definition: CollectionOptimisation.cpp:351
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:47
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:55
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
Nektar::SpatialDomains::TriGeom::kNedges
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:72
Nektar::Collections::OperatorImpMap
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:112
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LibUtilities::eGaussLobattoLegendre
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::LocalRegions::TetExpSharedPtr
std::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:214
Nektar::SpatialDomains::SegGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
Nektar::SpatialDomains::TetGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:85
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
Nektar::SpatialDomains::TriGeomSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::StdRegions::StdTetExpSharedPtr
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:255
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:408
Nektar::TetCollectionTests::BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(TestTetBwdTrans_IterPerExp_UniformP)
Definition: TestTetCollection.cpp:94
Nektar::TetCollectionTests::CreateTet
SpatialDomains::TetGeomSharedPtr CreateTet(SpatialDomains::PointGeomSharedPtr v0, SpatialDomains::PointGeomSharedPtr v1, SpatialDomains::PointGeomSharedPtr v2, SpatialDomains::PointGeomSharedPtr v3)
Definition: TestTetCollection.cpp:56
Nektar::TetCollectionTests::CreateSegGeom
SpatialDomains::SegGeomSharedPtr CreateSegGeom(unsigned int id, SpatialDomains::PointGeomSharedPtr v0, SpatialDomains::PointGeomSharedPtr v1)
Definition: TestTetCollection.cpp:46
Nektar::VarcoeffHashingTest::factors
StdRegions::ConstFactorMap factors
Definition: TestVarcoeffHashing.cpp:52
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:354
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:298
Generated on Wed Jan 31 2024 08:26:49 for Nektar++ by doxygen 1.9.4