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