PyrGeom.cpp 27.5 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1
2
////////////////////////////////////////////////////////////////////////////////
//
3
//  File: PyrGeom.cpp
Mike Kirby's avatar
Mike Kirby committed
4
5
6
7
8
9
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10
//  Department of Aeronautics, Imperial College London (UK), and Scientific
Mike Kirby's avatar
Mike Kirby committed
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
32
//  Description: Pyramidic geometry information.
Mike Kirby's avatar
Mike Kirby committed
33
34
35
36
37
//
////////////////////////////////////////////////////////////////////////////////

#include <SpatialDomains/PyrGeom.h>

38
39
40
41
42
43
44
#include <SpatialDomains/Geometry1D.h>
#include <SpatialDomains/Geometry2D.h>
#include <StdRegions/StdPyrExp.h>
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/TriGeom.h>
#include <SpatialDomains/MeshComponents.h>

45
using namespace std;
46

Mike Kirby's avatar
Mike Kirby committed
47
48
49
50
namespace Nektar
{
    namespace SpatialDomains
    {
51
52
        PyrGeom::PyrGeom()
        {
53
            m_shapeType = LibUtilities::ePyramid;
54
55
        }

56
57
        PyrGeom::PyrGeom(const Geometry2DSharedPtr faces[]) :
            Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
58
        {
59
            m_shapeType = LibUtilities::ePyramid;
60

61
62
63
            /// Copy the face shared pointers
            m_faces.insert(m_faces.begin(), faces, faces+PyrGeom::kNfaces);

64
65
66
67
            /// Set up orientation vectors with correct amount of elements.
            m_eorient.resize(kNedges);
            m_forient.resize(kNfaces);

68
69
70
71
            SetUpLocalEdges();
            SetUpLocalVertices();
            SetUpEdgeOrientation();
            SetUpFaceOrientation();
72
            SetUpXmap();
73
            SetUpCoeffs(m_xmap->GetNcoeffs());
74
        }
75

76
        PyrGeom::~PyrGeom()
Mike Kirby's avatar
Mike Kirby committed
77
        {
78

79
        }
80

81
82
83
84
        void PyrGeom::v_GenGeomFactors()
        {
            if (m_geomFactorsState != ePtsFilled)
            {
Dave Moxey's avatar
Dave Moxey committed
85
                int i;
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
                GeomType Gtype = eRegular;

                v_FillGeom();

                // check to see if expansions are linear
                for(i = 0; i < m_coordim; ++i)
                {
                    if (m_xmap->GetBasisNumModes(0) != 2 ||
                        m_xmap->GetBasisNumModes(1) != 2 ||
                        m_xmap->GetBasisNumModes(2) != 2 )
                    {
                        Gtype = eDeformed;
                    }
                }

                // check to see if all quadrilateral faces are parallelograms
                if(Gtype == eRegular)
                {
                    // Ensure each face is a parallelogram? Check this.
                    for (i = 0; i < m_coordim; i++)
                    {
                        if( fabs( (*m_verts[0])(i) -
                                  (*m_verts[1])(i) +
                                  (*m_verts[2])(i) -
                                  (*m_verts[3])(i) )
                            > NekConstants::kNekZeroTol )
                        {
                            Gtype = eDeformed;
                            break;
                        }
                    }
                }

                m_geomFactors = MemoryManager<GeomFactors>::AllocateSharedPtr(
                    Gtype, m_coordim, m_xmap, m_coeffs);
                m_geomFactorsState = ePtsFilled;
            }
        }

125
        NekDouble PyrGeom::v_GetLocCoords(
126
127
128
            const Array<OneD, const NekDouble> &coords,
                  Array<OneD,       NekDouble> &Lcoords)
        {
129
            NekDouble ptdist = 1e6;
Chris Cantwell's avatar
Chris Cantwell committed
130

131
            v_FillGeom();
132

133
            // calculate local coordinate for coord
134
            if(GetMetricInfo()->GetGtype() == eRegular)
135
            {   // Based on Spen's book, page 99
136

137
                // Point inside tetrahedron
138
                PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);
139

140
                // Edges
141
                PointGeom er0, e10, e30, e40;
142
143
144
145
                er0.Sub(r,*m_verts[0]);
                e10.Sub(*m_verts[1],*m_verts[0]);
                e30.Sub(*m_verts[3],*m_verts[0]);
                e40.Sub(*m_verts[4],*m_verts[0]);
146

147

148
                // Cross products (Normal times area)
149
                PointGeom cp1030, cp3040, cp4010;
150
151
152
                cp1030.Mult(e10,e30);
                cp3040.Mult(e30,e40);
                cp4010.Mult(e40,e10);
153
154


155
156
157
158
159
160
161
162
                // Barycentric coordinates (relative volume)
                NekDouble V = e40.dot(cp1030); // Pyramid Volume = {(e40)dot(e10)x(e30)}/4
                NekDouble scaleFactor = 2.0/3.0;
                NekDouble v1 = er0.dot(cp3040) / V; // volume1 = {(er0)dot(e30)x(e40)}/6
                NekDouble v2 = er0.dot(cp4010) / V; // volume2 = {(er0)dot(e40)x(e10)}/6
                NekDouble beta  = v1 * scaleFactor;
                NekDouble gamma = v2 * scaleFactor;
                NekDouble delta = er0.dot(cp1030) / V; // volume3 = {(er0)dot(e10)x(e30)}/4
163

164
165
166
167
                // Make Pyramid bigger
                Lcoords[0] = 2.0*beta  - 1.0;
                Lcoords[1] = 2.0*gamma - 1.0;
                Lcoords[2] = 2.0*delta - 1.0;
168

169
                // Set ptdist to distance to nearest vertex
170
171
172
173
174
                for(int i = 0; i < 5; ++i)
                {
                    ptdist = min(ptdist,r.dist(*m_verts[i]));
                }

175
            }
176
            else
177
            {
178
179
                NEKERROR(ErrorUtil::efatal,
                         "inverse mapping must be set up to use this call");
180
            }
181
            return ptdist;
Mike Kirby's avatar
Mike Kirby committed
182
183
        }

184
        int PyrGeom::v_GetNumVerts() const
185
        {
186
            return 5;
187
        }
188

189
        int PyrGeom::v_GetNumEdges() const
Mike Kirby's avatar
Mike Kirby committed
190
        {
191
            return 8;
Mike Kirby's avatar
Mike Kirby committed
192
        }
193

194
195
196
197
198
        int PyrGeom::v_GetNumFaces() const
        {
            return 5;
        }

199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
        int PyrGeom::v_GetDir(const int faceidx, const int facedir) const
        {
            if (faceidx == 0)
            {
                return facedir;
            }
            else if (faceidx == 1 || faceidx == 3)
            {
                return 2 * facedir;
            }
            else
            {
                return 1 + facedir;
            }
        }

215
216
        void PyrGeom::SetUpLocalEdges()
        {
217
218
219
220
221
222
223
            // find edge 0
            int i,j;
            unsigned int check;

            SegGeomSharedPtr edge;

            // First set up the 4 bottom edges
224
225
            int f;
            for (f = 1; f < 5; f++)
226
            {
227
                int nEdges = m_faces[f]->GetNumEdges();
228
                check = 0;
229
                for (i = 0; i < 4; i++)
230
                {
231
                    for (j = 0; j < nEdges; j++)
232
                    {
233
                        if (m_faces[0]->GetEid(i) == m_faces[f]->GetEid(j))
234
235
236
237
238
239
240
                        {
                            edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[0])->GetEdge(i));
                            m_edges.push_back(edge);
                            check++;
                        }
                    }
                }
241

242
                if (check < 1)
243
244
245
                {
                    std::ostringstream errstrm;
                    errstrm << "Connected faces do not share an edge. Faces ";
246
                    errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
247
248
                    ASSERTL0(false, errstrm.str());
                }
249
                else if (check > 1)
250
251
                {
                    std::ostringstream errstrm;
Blake Nelson's avatar
Blake Nelson committed
252
                    errstrm << "Connected faces share more than one edge. Faces ";
253
                    errstrm << (m_faces[0])->GetFid() << ", " << (m_faces[f])->GetFid();
254
                    ASSERTL0(false, errstrm.str());
255
256
257
                }
            }

258
259
            // Then, set up the 4 vertical edges
            check = 0;
260
            for(i = 0; i < 3; i++) //Set up the vertical edge :face(1) and face(4)
261
            {
262
                for(j = 0; j < 3; j++)
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
                {
                    if( (m_faces[1])->GetEid(i) == (m_faces[4])->GetEid(j) )
                    {
                        edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[1])->GetEdge(i));
                        m_edges.push_back(edge);
                        check++;
                    }
                }
            }
            if( check < 1 )
            {
                std::ostringstream errstrm;
                errstrm << "Connected faces do not share an edge. Faces ";
                errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
                ASSERTL0(false, errstrm.str());
            }
            else if( check > 1)
            {
                std::ostringstream errstrm;
Blake Nelson's avatar
Blake Nelson committed
282
                errstrm << "Connected faces share more than one edge. Faces ";
283
284
285
                errstrm << (m_faces[1])->GetFid() << ", " << (m_faces[4])->GetFid();
                ASSERTL0(false, errstrm.str());
            }
286

287
            // Set up vertical edges: face(1) through face(4)
288
            for (f = 1; f < 4; f++)
289
290
            {
                check = 0;
291
                for(i = 0; i < m_faces[f]->GetNumEdges(); i++)
292
                {
293
                    for(j = 0; j < m_faces[f+1]->GetNumEdges(); j++)
294
                    {
295
                        if( (m_faces[f])->GetEid(i) == (m_faces[f+1])->GetEid(j))
296
                        {
297
                            edge = boost::dynamic_pointer_cast<SegGeom>((m_faces[f])->GetEdge(i));
298
299
300
301
302
                            m_edges.push_back(edge);
                            check++;
                        }
                    }
                }
303

304
305
306
307
                if( check < 1 )
                {
                    std::ostringstream errstrm;
                    errstrm << "Connected faces do not share an edge. Faces ";
308
                    errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
309
310
311
312
313
                    ASSERTL0(false, errstrm.str());
                }
                else if( check > 1)
                {
                    std::ostringstream errstrm;
Blake Nelson's avatar
Blake Nelson committed
314
                    errstrm << "Connected faces share more than one edge. Faces ";
315
                    errstrm << (m_faces[f])->GetFid() << ", " << (m_faces[f+1])->GetFid();
316
                    ASSERTL0(false, errstrm.str());
317
                }
318
            }
319
        }
320

321
322
        void PyrGeom::SetUpLocalVertices()
        {
323
            // Set up the first 2 vertices (i.e. vertex 0,1)
324
325
            if (m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0) ||
                m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1))
326
327
328
329
            {
                m_verts.push_back(m_edges[0]->GetVertex(1));
                m_verts.push_back(m_edges[0]->GetVertex(0));
            }
330
331
            else if (m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0) ||
                     m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1))
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
            {
                m_verts.push_back(m_edges[0]->GetVertex(0));
                m_verts.push_back(m_edges[0]->GetVertex(1));
            }
            else
            {
                std::ostringstream errstrm;
                errstrm << "Connected edges do not share a vertex. Edges ";
                errstrm << m_edges[0]->GetEid() << ", " << m_edges[1]->GetEid();
                ASSERTL0(false, errstrm.str());
            }

            // set up the other bottom vertices (i.e. vertex 2,3)
            for(int i = 1; i < 3; i++)
            {
347
                if (m_edges[i]->GetVid(0) == m_verts[i]->GetVid())
348
349
350
                {
                    m_verts.push_back(m_edges[i]->GetVertex(1));
                }
351
                else if (m_edges[i]->GetVid(1) == m_verts[i]->GetVid())
352
353
354
355
356
357
358
359
360
361
362
363
                {
                    m_verts.push_back(m_edges[i]->GetVertex(0));
                }
                else
                {
                    std::ostringstream errstrm;
                    errstrm << "Connected edges do not share a vertex. Edges ";
                    errstrm << m_edges[i]->GetEid() << ", " << m_edges[i-1]->GetEid();
                    ASSERTL0(false, errstrm.str());
                }
            }

364
            // set up top vertex
365
            if (m_edges[4]->GetVid(0) == m_verts[0]->GetVid())
366
367
368
369
370
            {
                m_verts.push_back(m_edges[4]->GetVertex(1));
            }
            else
            {
371
372
                m_verts.push_back(m_edges[4]->GetVertex(0));
            }
373

374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
            int check = 0;
            for (int i = 5; i < 8; ++i)
            {
                if( (m_edges[i]->GetVid(0) == m_verts[i-4]->GetVid()
                    && m_edges[i]->GetVid(1) == m_verts[4]->GetVid())
                    ||(m_edges[i]->GetVid(1) == m_verts[i-4]->GetVid()
                    && m_edges[i]->GetVid(0) == m_verts[4]->GetVid()))
                {
                    check++;
                }
            }
            if (check != 3) {
                std::ostringstream errstrm;
                errstrm << "Connected edges do not share a vertex. Edges ";
                errstrm << m_edges[3]->GetEid() << ", " << m_edges[2]->GetEid();
389
390
                ASSERTL0(false, errstrm.str());
            }
391
        }
392

393
394
        void PyrGeom::SetUpEdgeOrientation()
        {
395
396
397
            // This 2D array holds the local id's of all the vertices for every
            // edge. For every edge, they are ordered to what we define as being
            // Forwards.
398
            const unsigned int edgeVerts[kNedges][2] =
399
                { {0,1}, {1,2}, {3,2}, {0,3}, {0,4}, {1,4}, {2,4}, {3,4} };
400

401
            int i;
402
            for (i = 0; i < kNedges; i++)
403
            {
404
                if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][0]]->GetVid())
405
406
407
                {
                    m_eorient[i] = StdRegions::eForwards;
                }
408
                else if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][1]]->GetVid())
409
410
411
412
413
414
415
416
                {
                    m_eorient[i] = StdRegions::eBackwards;
                }
                else
                {
                    ASSERTL0(false,"Could not find matching vertex for the edge");
                }
            }
417
        }
418

419
420
        void PyrGeom::SetUpFaceOrientation()
        {
421
422
            int f,i;

423
424
425
426
427
            // These arrays represent the vector of the A and B coordinate of
            // the local elemental coordinate system where A corresponds with
            // the coordinate direction xi_i with the lowest index i (for that
            // particular face) Coordinate 'B' then corresponds to the other
            // local coordinate (i.e. with the highest index)
428
429
430
431
432
433
434
435
436
437
            Array<OneD,NekDouble> elementAaxis(m_coordim);
            Array<OneD,NekDouble> elementBaxis(m_coordim);

            // These arrays correspond to the local coordinate
            // system of the face itself (i.e. the Geometry2D)
            // faceAaxis correspond to the xi_0 axis
            // faceBaxis correspond to the xi_1 axis
            Array<OneD,NekDouble> faceAaxis(m_coordim);
            Array<OneD,NekDouble> faceBaxis(m_coordim);

438
439
440
            // This is the base vertex of the face (i.e. the Geometry2D) This
            // corresponds to thevertex with local ID 0 of the Geometry2D
            unsigned int baseVertex;
441
442
443
444
445
446
447

            // The lenght of the vectors above
            NekDouble elementAaxis_length;
            NekDouble elementBaxis_length;
            NekDouble faceAaxis_length;
            NekDouble faceBaxis_length;

448
449
450
451
452
453
454
455
456
457
            // This 2D array holds the local id's of all the vertices for every
            // face. For every face, they are ordered in such a way that the
            // implementation below allows a unified approach for all faces.
            const unsigned int faceVerts[kNfaces][4] = {
                {0,1,2,3},
                {0,1,4,0}, // Last four elements are triangles which only
                {1,2,4,0}, // require three vertices.
                {3,2,4,0},
                {0,3,4,0}
            };
458
459
460
461
462
463

            NekDouble dotproduct1 = 0.0;
            NekDouble dotproduct2 = 0.0;

            unsigned int orientation;

464
            // Loop over all the faces to set up the orientation
465
466
467
468
469
470
471
            for(f = 0; f < kNqfaces + kNtfaces; f++)
            {
                // initialisation
                elementAaxis_length = 0.0;
                elementBaxis_length = 0.0;
                faceAaxis_length = 0.0;
                faceBaxis_length = 0.0;
472

473
474
475
                dotproduct1 = 0.0;
                dotproduct2 = 0.0;

476
                baseVertex = m_faces[f]->GetVid(0);
477

478
479
480
481
482
483
484
                // We are going to construct the vectors representing the A and
                // B axis of every face. These vectors will be constructed as a
                // vector-representation of the edges of the face. However, for
                // both coordinate directions, we can represent the vectors by
                // two different edges. That's why we need to make sure that we
                // pick the edge to which the baseVertex of the
                // Geometry2D-representation of the face belongs...
485
486

                // Compute the length of edges on a base-face
487
                if (f > 0)
488
489
                {
                    for(i = 0; i < m_coordim; i++)
490
                    {
491
                        elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
492
                        elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
493
                    }
494
                }
495
                else
496
                {
497
                    if( baseVertex == m_verts[ faceVerts[f][0] ]->GetVid() )
498
                    {
499
500
501
502
503
                        for(i = 0; i < m_coordim; i++)
                        {
                            elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                            elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                        }
504
                    }
505
                    else if( baseVertex == m_verts[ faceVerts[f][1] ]->GetVid() )
506
                    {
507
508
509
510
511
                        for(i = 0; i < m_coordim; i++)
                        {
                            elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                            elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
                        }
512
                    }
513
                    else if( baseVertex == m_verts[ faceVerts[f][2] ]->GetVid() )
514
                    {
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
                        for(i = 0; i < m_coordim; i++)
                        {
                            elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
                            elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][1] ])[i];
                        }
                    }
                    else if( baseVertex == m_verts[ faceVerts[f][3] ]->GetVid() )
                    {
                        for(i = 0; i < m_coordim; i++)
                        {
                            elementAaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][3] ])[i];
                            elementBaxis[i] = (*m_verts[ faceVerts[f][3] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
                        }
                    }
                    else
                    {
                        ASSERTL0(false, "Could not find matching vertex for the face");
532
                    }
533
                }
534

535
                // Now, construct the edge-vectors of the local coordinates of
536
537
                // the Geometry2D-representation of the face
                for(i = 0; i < m_coordim; i++)
538
                {
539
                    int v = m_faces[f]->GetNumVerts()-1;
540
                    faceAaxis[i] = (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
541
                    faceBaxis[i] = (*m_faces[f]->GetVertex(v))[i] - (*m_faces[f]->GetVertex(0))[i];
542

543
544
545
546
547
                    elementAaxis_length += pow(elementAaxis[i],2);
                    elementBaxis_length += pow(elementBaxis[i],2);
                    faceAaxis_length += pow(faceAaxis[i],2);
                    faceBaxis_length += pow(faceBaxis[i],2);
                }
548

549
550
551
552
553
554
555
556
557
558
559
560
561
                elementAaxis_length = sqrt(elementAaxis_length);
                elementBaxis_length = sqrt(elementBaxis_length);
                faceAaxis_length = sqrt(faceAaxis_length);
                faceBaxis_length = sqrt(faceBaxis_length);

                // Calculate the inner product of both the A-axis
                // (i.e. Elemental A axis and face A axis)
                for(i = 0 ; i < m_coordim; i++)
                {
                    dotproduct1 += elementAaxis[i]*faceAaxis[i];
                }

                orientation = 0;
562
                // if the innerproduct is equal to the (absolute value of the ) products of the lengths
563
                // of both vectors, then, the coordinate systems will NOT be transposed
Peter Vos's avatar
Peter Vos committed
564
                if( fabs(elementAaxis_length*faceAaxis_length - fabs(dotproduct1)) < NekConstants::kNekZeroTol )
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
                {
                    // if the inner product is negative, both A-axis point
                    // in reverse direction
                    if(dotproduct1 < 0.0)
                    {
                        orientation += 2;
                    }

                    // calculate the inner product of both B-axis
                    for(i = 0 ; i < m_coordim; i++)
                    {
                        dotproduct2 += elementBaxis[i]*faceBaxis[i];
                    }

                    // if the inner product is negative, both B-axis point
                    // in reverse direction
                    if( dotproduct2 < 0.0 )
                    {
                        orientation++;
                    }
                }
                // The coordinate systems are transposed
                else
                {
                    orientation = 4;

                    // Calculate the inner product between the elemental A-axis
                    // and the B-axis of the face (which are now the corresponding axis)
                    dotproduct1 = 0.0;
                    for(i = 0 ; i < m_coordim; i++)
                    {
                        dotproduct1 += elementAaxis[i]*faceBaxis[i];
                    }
598

599
                    // check that both these axis are indeed parallel
600
601
602
                    if (fabs(elementAaxis_length*faceBaxis_length
                            - fabs(dotproduct1)) > NekConstants::kNekZeroTol)
                    {
603
                        cout << "Warning: Pyramid axes not parallel" << endl;
604
                    }
605

606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
                    // if the result is negative, both axis point in reverse
                    // directions
                    if(dotproduct1 < 0.0)
                    {
                        orientation += 2;
                    }

                    // Do the same for the other two corresponding axis
                    dotproduct2 = 0.0;
                    for(i = 0 ; i < m_coordim; i++)
                    {
                        dotproduct2 += elementBaxis[i]*faceAaxis[i];
                    }

                    // check that both these axis are indeed parallel
621
622
623
                    if (fabs(elementBaxis_length*faceAaxis_length
                            - fabs(dotproduct2)) > NekConstants::kNekZeroTol)
                    {
624
                        cout << "Warning: Pyramid axes not parallel" << endl;
625
                    }
626
627
628
629
630
631

                    if( dotproduct2 < 0.0 )
                    {
                        orientation++;
                    }
                }
632

633
                orientation = orientation + 5;
634
635

                // Fill the m_forient array
636
                m_forient[f] = (StdRegions::Orientation) orientation;
637
            }
638
        }
639
640

        void PyrGeom::v_Reset(
641
642
            CurveMap &curvedEdges,
            CurveMap &curvedFaces)
643
644
645
646
647
648
649
650
651
652
653
        {
            Geometry::v_Reset(curvedEdges, curvedFaces);

            for (int i = 0; i < 5; ++i)
            {
                m_faces[i]->Reset(curvedEdges, curvedFaces);
            }

            SetUpXmap();
            SetUpCoeffs(m_xmap->GetNcoeffs());
        }
654

655
656
657
658
659
660
661
        /**
         * @brief Set up the #m_xmap object by determining the order of each
         * direction from derived faces.
         */
        void PyrGeom::SetUpXmap()
        {
            vector<int> tmp;
662
            int order0, order1;
663
664
665
666
667
668
669
670
671
672
673
674
675

            if (m_forient[0] < 9)
            {
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
                order0 = *max_element(tmp.begin(), tmp.end());
            }
            else
            {
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
                order0 = *max_element(tmp.begin(), tmp.end());
            }
676

677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
            if (m_forient[0] < 9)
            {
                tmp.clear();
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(3));
                tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
                order1 = *max_element(tmp.begin(), tmp.end());
            }
            else
            {
                tmp.clear();
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
                tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
                tmp.push_back(m_faces[2]->GetXmap()->GetEdgeNcoeffs(2));
                order1 = *max_element(tmp.begin(), tmp.end());
            }
693

694
695
696
697
698
699
700
701
            tmp.clear();
            tmp.push_back(order0);
            tmp.push_back(order1);
            tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(1));
            tmp.push_back(m_faces[1]->GetXmap()->GetEdgeNcoeffs(2));
            tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(1));
            tmp.push_back(m_faces[3]->GetXmap()->GetEdgeNcoeffs(2));
            int order2 = *max_element(tmp.begin(), tmp.end());
702
703


704
705
706
            const LibUtilities::BasisKey A1(
                LibUtilities::eModified_A, order0,
                LibUtilities::PointsKey(
707
                    order0+1, LibUtilities::eGaussLobattoLegendre));
708
709
710
            const LibUtilities::BasisKey A2(
                LibUtilities::eModified_A, order1,
                LibUtilities::PointsKey(
711
                    order1+1, LibUtilities::eGaussLobattoLegendre));
712
            const LibUtilities::BasisKey C(
713
                LibUtilities::eModifiedPyr_C, order2,
714
                LibUtilities::PointsKey(
715
                    order2, LibUtilities::eGaussRadauMAlpha2Beta0));
716
717
718
719

            m_xmap = MemoryManager<StdRegions::StdPyrExp>::AllocateSharedPtr(
                A1, A2, C);
        }
Mike Kirby's avatar
Mike Kirby committed
720
721
    }; //end of namespace
}; //end of namespace