TriGeom.cpp 32.5 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1
2
////////////////////////////////////////////////////////////////////////////////
//
Peter Vos's avatar
Peter Vos committed
3
//  File: TriGeom.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:
Mike Kirby's avatar
Mike Kirby committed
33
34
35
36
37
//
//
////////////////////////////////////////////////////////////////////////////////

#include <SpatialDomains/TriGeom.h>
38
39
#include <StdRegions/StdNodalTriExp.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
40
#include <LibUtilities/Foundations/Interp.h>
Mike Kirby's avatar
Mike Kirby committed
41

42
43
44
45
#include <SpatialDomains/SegGeom.h>
#include <SpatialDomains/GeomFactors2D.h>


Mike Kirby's avatar
Mike Kirby committed
46
47
48
49
namespace Nektar
{
    namespace SpatialDomains
    {
50
51
52
        /**
         *
         */
Sehun Chun's avatar
Sehun Chun committed
53
        TriGeom::TriGeom()
Mike Kirby's avatar
Mike Kirby committed
54
        {
55
            m_shapeType = LibUtilities::eTriangle;
Mike Kirby's avatar
Mike Kirby committed
56
57
        }

58
59
60
61

        /**
         *
         */
62
        TriGeom::TriGeom(int id, const int coordim):
63
                                   Geometry2D(coordim), m_fid(id)
64
65
        {
            const LibUtilities::BasisKey B0(LibUtilities::eModified_A, 2,
66
                    LibUtilities::PointsKey(3,LibUtilities::eGaussLobattoLegendre));
67
            const LibUtilities::BasisKey B1(LibUtilities::eModified_B, 2,
68
                    LibUtilities::PointsKey(3,LibUtilities::eGaussRadauMAlpha1Beta0));
69
70

            m_xmap = Array<OneD, StdRegions::StdExpansion2DSharedPtr>(coordim);
71

72
73
74
75
            m_fid = id;

            for(int i = 0; i < m_coordim; ++i)
            {
76
                m_xmap[i] = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(B0,B1);
77
78
79
80
            }

        }

81
82
83
84

        /**
         *
         */
85
        TriGeom::TriGeom(const int id,
86
87
                const VertexComponentSharedPtr verts[],
                const SegGeomSharedPtr edges[],
88
                const StdRegions::Orientation eorient[]):
89
90
                Geometry2D(verts[0]->GetCoordim()),
                m_fid(id)
91
        {
92
            m_shapeType = LibUtilities::eTriangle;
93

94
            /// Copy the vert shared pointers.
95
            m_verts.insert(m_verts.begin(), verts, verts+TriGeom::kNverts);
96
97
98
99
100
101
102
103
104
105
106

            /// Copy the edge shared pointers.
            m_edges.insert(m_edges.begin(), edges, edges+TriGeom::kNedges);

            for (int j=0; j<kNedges; ++j)
            {
                m_eorient[j] = eorient[j];
            }

            m_coordim = verts[0]->GetCoordim();
            ASSERTL0(m_coordim > 1,
107
                    "Cannot call function with dim == 1");
108

Spencer Sherwin's avatar
Spencer Sherwin committed
109
            int order0  = edges[0]->GetBasis(0,0)->GetNumModes();
110
            int points0 = edges[0]->GetBasis(0,0)->GetNumPoints();
111
            int order1  = max(order0,
112
113
                    max(edges[1]->GetBasis(0,0)->GetNumModes(),
                            edges[2]->GetBasis(0,0)->GetNumModes()));
Spencer Sherwin's avatar
Spencer Sherwin committed
114
            int points1 = max(points0,
115
116
117
                    max(edges[1]->GetBasis(0,0)->GetNumPoints(),
                            edges[2]->GetBasis(0,0)->GetNumPoints()));

Spencer Sherwin's avatar
Spencer Sherwin committed
118
119

            const LibUtilities::BasisKey B0(LibUtilities::eModified_A, order0,
120
                    LibUtilities::PointsKey(points0,LibUtilities::eGaussLobattoLegendre));
Spencer Sherwin's avatar
Spencer Sherwin committed
121
            const LibUtilities::BasisKey B1(LibUtilities::eModified_B, order1,
122
                    LibUtilities::PointsKey(points1,LibUtilities::eGaussRadauMAlpha1Beta0));
123
124
125
126
127

            m_xmap = Array<OneD, StdRegions::StdExpansion2DSharedPtr>(m_coordim);

            for(int i = 0; i < m_coordim; ++i)
            {
Spencer Sherwin's avatar
Spencer Sherwin committed
128
                m_xmap[i] = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(B0,B1);
129
            }
130
131
        }

132
133
134
135

        /**
         *
         */
136
        TriGeom::TriGeom(const int id, const SegGeomSharedPtr edges[],
137
                const StdRegions::Orientation eorient[]):
138
139
                Geometry2D(edges[0]->GetVertex(0)->GetCoordim()),
                m_fid(id)
Mike Kirby's avatar
Mike Kirby committed
140
        {
141
            m_shapeType = LibUtilities::eTriangle;
142

Mike Kirby's avatar
Mike Kirby committed
143
144
145
            /// Copy the edge shared pointers.
            m_edges.insert(m_edges.begin(), edges, edges+TriGeom::kNedges);

146
147
148
149
150
151
152
153
154
155
156
            for(int j=0; j <kNedges; ++j)
            {
                if(eorient[j] == StdRegions::eForwards)
                {
                    m_verts.push_back(edges[j]->GetVertex(0));
                }
                else
                {
                    m_verts.push_back(edges[j]->GetVertex(1));
                }
            }
157

158
            for (int j=0; j<kNedges; ++j)
Spencer Sherwin's avatar
   
Spencer Sherwin committed
159
160
161
162
            {
                m_eorient[j] = eorient[j];
            }

Mike Kirby's avatar
Mike Kirby committed
163
            m_coordim = edges[0]->GetVertex(0)->GetCoordim();
164
            ASSERTL0(m_coordim > 1,"Cannot call function with dim == 1");
165

Spencer Sherwin's avatar
Spencer Sherwin committed
166
            int order0  = edges[0]->GetBasis(0,0)->GetNumModes();
167
168
            int points0 = edges[0]->GetBasis(0,0)->GetNumPoints();
            int order1  = max(order0,
169
170
                    max(edges[1]->GetBasis(0,0)->GetNumModes(),
                            edges[2]->GetBasis(0,0)->GetNumModes()));
Spencer Sherwin's avatar
Spencer Sherwin committed
171
            int points1 = max(points0,
172
173
174
                    max(edges[1]->GetBasis(0,0)->GetNumPoints(),
                            edges[2]->GetBasis(0,0)->GetNumPoints()));

175
            const LibUtilities::BasisKey B0(LibUtilities::eModified_A, order0,
176
                    LibUtilities::PointsKey(points0,LibUtilities::eGaussLobattoLegendre));
177
            const LibUtilities::BasisKey B1(LibUtilities::eModified_B, order1,
178
                    LibUtilities::PointsKey(points1,LibUtilities::eGaussRadauMAlpha1Beta0));
179
180
181
182
183
184
185
186
187

            m_xmap = Array<OneD, StdRegions::StdExpansion2DSharedPtr>(m_coordim);

            for(int i = 0; i < m_coordim; ++i)
            {
                m_xmap[i] = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(B0,B1);
            }
        }

188
189
190
191

        /**
         *
         */
192
        TriGeom::TriGeom(const int id,
193
                const SegGeomSharedPtr edges[],
194
                const StdRegions::Orientation eorient[],
195
196
197
                const CurveSharedPtr &curve) :
                Geometry2D(edges[0]->GetVertex(0)->GetCoordim()),
                m_fid(id)
198
        {
199
            m_shapeType =  LibUtilities::eTriangle;
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226

            /// Copy the edge shared pointers.
            m_edges.insert(m_edges.begin(), edges, edges+TriGeom::kNedges);

            for(int j=0; j <kNedges; ++j)
            {
                if(eorient[j] == StdRegions::eForwards)
                {
                    m_verts.push_back(edges[j]->GetVertex(0));
                }
                else
                {
                    m_verts.push_back(edges[j]->GetVertex(1));
                }
            }

            for (int j=0; j<kNedges; ++j)
            {
                m_eorient[j] = eorient[j];
            }

            m_coordim = edges[0]->GetVertex(0)->GetCoordim();
            ASSERTL0(m_coordim > 1,"Cannot call function with dim == 1");

            int order0  = edges[0]->GetBasis(0,0)->GetNumModes();
            int points0 = edges[0]->GetBasis(0,0)->GetNumPoints();
            int order1  = max(order0,
227
228
                    max(edges[1]->GetBasis(0,0)->GetNumModes(),
                            edges[2]->GetBasis(0,0)->GetNumModes()));
229
            int points1 = max(points0,
230
231
                    max(edges[1]->GetBasis(0,0)->GetNumPoints(),
                            edges[2]->GetBasis(0,0)->GetNumPoints()));
Spencer Sherwin's avatar
Spencer Sherwin committed
232
233

            const LibUtilities::BasisKey B0(LibUtilities::eModified_A, order0,
234
                    LibUtilities::PointsKey(points0,LibUtilities::eGaussLobattoLegendre));
Spencer Sherwin's avatar
Spencer Sherwin committed
235
            const LibUtilities::BasisKey B1(LibUtilities::eModified_B, order1,
236
                    LibUtilities::PointsKey(points1,LibUtilities::eGaussRadauMAlpha1Beta0));
237

238
            m_xmap = Array<OneD, StdRegions::StdExpansion2DSharedPtr>(m_coordim);
239

240
241
            for(int i = 0; i < m_coordim; ++i)
            {
Spencer Sherwin's avatar
Spencer Sherwin committed
242
                m_xmap[i] = MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(B0,B1);
243

244
                int pdim = LibUtilities::PointsManager()[LibUtilities::PointsKey(2, curve->m_ptype)]->GetPointsDim();
245

246
                // Deal with 2D points type separately (e.g. electrostatic or
247
248
249
250
                // Fekete points).
                if (pdim == 2)
                {
                    int N = curve->m_points.size();
251
                    int nEdgePts = (-1+(int)sqrt(static_cast<NekDouble>(8*N+1)))/2;
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
                    
                    ASSERTL0(nEdgePts*(nEdgePts+1)/2 == N,
                             "NUMPOINTS must be a triangle number for 2D basis.");
                    
                    for (int j = 0; j < kNedges; ++j)
                    {
                        ASSERTL0(edges[j]->GetXmap(i)->GetNcoeffs() == nEdgePts,
                                 "Number of edge points does not correspond "
                                 "to number of face points.");
                    }
                    
                    // Create a StdNodalTriExp.
                    const LibUtilities::PointsKey P0(
                        nEdgePts, LibUtilities::eGaussLobattoLegendre);
                    const LibUtilities::PointsKey P1(
                        nEdgePts, LibUtilities::eGaussRadauMAlpha1Beta0);
                    const LibUtilities::BasisKey  T0(
                        LibUtilities::eOrtho_A, nEdgePts, P0);
                    const LibUtilities::BasisKey  T1(
                        LibUtilities::eOrtho_B, nEdgePts, P1);
                    
                    StdRegions::StdNodalTriExpSharedPtr t =
                        MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
                            T0,T1,curve->m_ptype);
                    
                    for (int j = 0; j < N; ++j)
                    {
                        t->UpdatePhys()[j] = (curve->m_points[j]->GetPtr())[i];
                    }
                    
                    Array<OneD, NekDouble> tmp(nEdgePts*nEdgePts);
                    t->BwdTrans(t->GetPhys(), tmp);
                    // Interpolate points to standard region.
                    LibUtilities::Interp2D(P0, P1, tmp,
                                           B0.GetPointsKey(),B1.GetPointsKey(),
                                           m_xmap[i]->UpdatePhys());
                    
                    // Forwards transform to get coefficient space.
                    m_xmap[i]->FwdTrans(m_xmap[i]->GetPhys(), m_xmap[i]->UpdateCoeffs());
                }
                else if (pdim == 1)
                {
                    int npts = curve->m_points.size();
295
                    int nEdgePts = (int)sqrt(static_cast<NekDouble>(npts));
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
                    Array<OneD,NekDouble> tmp(npts);
                    LibUtilities::PointsKey curveKey(nEdgePts, curve->m_ptype);
                    
                    // Sanity checks:
                    // - Curved faces should have square number of points;
                    // - Each edge should have sqrt(npts) points.
                    ASSERTL0(nEdgePts*nEdgePts == npts,
                             "NUMPOINTS should be a square number");
                
                    for (int j = 0; j < kNedges; ++j)
                    {
                        ASSERTL0(edges[j]->GetXmap(i)->GetNcoeffs() == nEdgePts,
                                 "Number of edge points does not correspond "
                                 "to number of face points.");
                    }
                    
                    for (int j = 0; j < npts; ++j)
                    {
                        tmp[j] = (curve->m_points[j]->GetPtr())[i];
                    }
                    
                    // Interpolate curve points to standard triangle points.
                    LibUtilities::Interp2D(curveKey,curveKey,tmp,
                                           B0.GetPointsKey(),B1.GetPointsKey(),
                                           m_xmap[i]->UpdatePhys());
                    
                    // Forwards transform to get coefficient space.
                    m_xmap[i]->FwdTrans(m_xmap[i]->GetPhys(),m_xmap[i]->UpdateCoeffs());
                }
                else
                {
                    ASSERTL0(false, "Only 1D/2D points distributions supported.");
                }
329
            }
Mike Kirby's avatar
Mike Kirby committed
330
331
        }

332

333
334
335
        /**
         *
         */
Spencer Sherwin's avatar
Spencer Sherwin committed
336
337
338
        TriGeom::TriGeom(const TriGeom &in)
        {
            // From Geomtry
339
            m_shapeType = in.m_shapeType;
340

Spencer Sherwin's avatar
Spencer Sherwin committed
341
342
            // From TriFaceComponent
            m_fid = in.m_fid;
343
            m_ownVerts = in.m_ownVerts;
344
            std::list<CompToElmt>::const_iterator def;
345
            for(def = in.m_elmtMap.begin(); def != in.m_elmtMap.end(); def++)
346
            {
347
                m_elmtMap.push_back(*def);
348
            }
349

Spencer Sherwin's avatar
Spencer Sherwin committed
350
351
352
            // From TriGeom
            m_verts = in.m_verts;
            m_edges = in.m_edges;
353
            for (int i = 0; i < kNedges; i++)
Spencer Sherwin's avatar
Spencer Sherwin committed
354
355
356
            {
                m_eorient[i] = in.m_eorient[i];
            }
357
            m_ownData = in.m_ownData;
Spencer Sherwin's avatar
Spencer Sherwin committed
358
        }
359

360
361
362
363

        /**
         *
         */
Mike Kirby's avatar
Mike Kirby committed
364
365
366
        TriGeom::~TriGeom()
        {
        }
367

368
369
370
371
372
373
374

        /**
         * Given local collapsed coordinate Lcoord return the value of physical
         * coordinate in direction i
         */
        NekDouble TriGeom::GetCoord(const int i,
                const Array<OneD, const NekDouble> &Lcoord)
Mike Kirby's avatar
Mike Kirby committed
375
        {
376
377
378
379
            ASSERTL1(m_state == ePtsFilled,
                    "Goemetry is not in physical space");
            return m_xmap[i]->PhysEvaluate(Lcoord);
        }
380

Mike Kirby's avatar
Mike Kirby committed
381

382
383
384
        /**
         * TODO: implement eight different case of face orientation
         */
385
        StdRegions::Orientation TriGeom::GetFaceOrientation(
386
387
388
                const TriGeom &face1,
                const TriGeom &face2)
        {
389
            StdRegions::Orientation returnval = StdRegions::eNoOrientation;
390
391
            
            int i, j, map[3] = {-1,-1,-1};
392
            NekDouble x, y, z, x1, y1, z1, cx = 0.0, cy = 0.0, cz = 0.0;
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
           
            // For periodic faces, we calculate the vector between the centre
            // points of the two faces. (For connected faces this will be
            // zero). We can then use this to determine alignment later in the
            // algorithm.
            for (i = 0; i < 3; ++i)
            {
                cx += (*face2.m_verts[i])(0) - (*face1.m_verts[i])(0);
                cy += (*face2.m_verts[i])(1) - (*face1.m_verts[i])(1);
                cz += (*face2.m_verts[i])(2) - (*face1.m_verts[i])(2);
            }
            cx /= 3;
            cy /= 3;
            cz /= 3;
           
            /*
            cout << "From face " << face1.GetGlobalID() << " -> " 
                 << face2.GetGlobalID() << endl;
            cout << "Vector c = " << cx << " " << cy << " " << cz << endl;
            */
            
            // Now construct a mapping which takes us from the vertices of one
            // face to the other. That is, vertex j of face2 corresponds to
            // vertex map[j] of face1.
            for (i = 0; i < 3; ++i)
            {
                x = (*face1.m_verts[i])(0);
                y = (*face1.m_verts[i])(1);
                z = (*face1.m_verts[i])(2);
                for (j = 0; j < 3; ++j)
                {
                    x1 = (*face2.m_verts[j])(0)-cx;
                    y1 = (*face2.m_verts[j])(1)-cy;
                    z1 = (*face2.m_verts[j])(2)-cz;
                    if (sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z1-z)*(z1-z)) < 1e-5)
                    {
                        map[j] = i;
                        break;
                    }
                }
            }
Mike Kirby's avatar
Mike Kirby committed
434

435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
            for (i = 0; i < 3; ++i)
            {
                ASSERTL0(map[i] != -1, "Unable to determine face mapping.");
            }
            
            // Use the mapping to determine the eight alignment options between
            // faces.
            if (map[1] == (map[0]+1) % 3)
            {
                switch (map[0])
                {
                    case 0:
                        returnval = StdRegions::eDir1FwdDir1_Dir2FwdDir2;
                        break;
                    case 1:
                        returnval = StdRegions::eDir1FwdDir2_Dir2BwdDir1;
                        break;
                    case 2:
                        returnval = StdRegions::eDir1BwdDir1_Dir2BwdDir2;
                        break;
                }
            }
            else 
            {
                switch (map[0])
                {
                    case 0:
                        returnval = StdRegions::eDir1FwdDir2_Dir2FwdDir1;
                        break;
                    case 1:
                        returnval = StdRegions::eDir1BwdDir1_Dir2FwdDir2;
                        break;
                    case 2:
                        returnval = StdRegions::eDir1BwdDir2_Dir2BwdDir1;
                        break;
                }
            }
            
473
            return returnval;
Mike Kirby's avatar
Mike Kirby committed
474
        }
475
        
476
477
478
479
480

        /**
         *
         */
        void TriGeom::v_AddElmtConnected(int gvo_id, int locid)
481
482
        {
            CompToElmt ee(gvo_id,locid);
483
            m_elmtMap.push_back(ee);
484
485
486
        }


487
488
489
490
        /**
         *
         */
        int  TriGeom::v_NumElmtConnected() const
491
        {
492
            return int(m_elmtMap.size());
493
494
495
        }


496
497
498
499
        /**
         *
         */
        bool TriGeom::v_IsElmtConnected(int gvo_id, int locid) const
500
501
502
503
        {
            std::list<CompToElmt>::const_iterator def;
            CompToElmt ee(gvo_id,locid);

504
            def = find(m_elmtMap.begin(),m_elmtMap.end(),ee);
505
506

            // Found the element connectivity object in the list
507
            return (def != m_elmtMap.end());
508
509
510
        }


511
512
513
514
515
516
517
        /**
         *
         */
        int  TriGeom::v_GetFid() const
        {
            return m_fid;
        }
518

519
520
521
522

        /**
         *
         */
523
        int  TriGeom::v_GetCoordim() const
524
        {
525
526
            return m_coordim;
        }
527

528
529
530
531
532
533
534

        /**
         *
         */
        const LibUtilities::BasisSharedPtr TriGeom::v_GetBasis(const int i, const int j)
        {
            return m_xmap[i]->GetBasis(j);
535
536
        }

Mike Kirby's avatar
Mike Kirby committed
537

538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
        /**
         *
         */
        const LibUtilities::BasisSharedPtr TriGeom::v_GetEdgeBasis(const int i, const int j)
        {
            ASSERTL1(j <=2,"edge is out of range");
            if(j == 0)
            {
                return m_xmap[i]->GetBasis(0);
            }
            else
            {
                return m_xmap[i]->GetBasis(1);
            }
        }


        /**
         *
         */
        Array<OneD,NekDouble> &TriGeom::v_UpdatePhys(const int i)
        {
            return m_xmap[i]->UpdatePhys();
        }


        /**
         *
         */
        NekDouble TriGeom::v_GetCoord(const int i, const Array<OneD,const NekDouble> &Lcoord)
        {
            return GetCoord(i,Lcoord);
        }


        /**
         * Set up GeoFac for this geometry using Coord quadrature distribution
         */
576
577
        void TriGeom::v_GenGeomFactors(
                const Array<OneD, const LibUtilities::BasisSharedPtr> &tbasis)
578
        {
579
580
581
            if (m_geomFactorsState != ePtsFilled)
            {
                GeomType Gtype = eRegular;
582

583
                TriGeom::v_FillGeom();
584

585
586
                // check to see if expansions are linear
                for(int i = 0; i < m_coordim; ++i)
587
                {
588
589
590
591
592
                    if((m_xmap[i]->GetBasisNumModes(0) != 2)||
                            (m_xmap[i]->GetBasisNumModes(1) != 2))
                    {
                        Gtype = eDeformed;
                    }
593
594
                }

595
                m_geomFactors = MemoryManager<GeomFactors2D>::AllocateSharedPtr(
596
                    Gtype, m_coordim, m_xmap, tbasis, true);
597
598
599

                m_geomFactorsState = ePtsFilled;
            }
600
601
602
603
604
605
606
607
608
609
        }


        /**
         *
         */
        void TriGeom::v_SetOwnData()
        {
            m_ownData = true;
        }
Mike Kirby's avatar
Mike Kirby committed
610
611


612
613
614
615
616
617
        /**
         * Note verts and edges are listed according to anticlockwise
         * convention but points in _coeffs have to be in array format from
         * left to right.
         */
        void TriGeom::v_FillGeom()
Mike Kirby's avatar
Mike Kirby committed
618
619
        {
            // check to see if geometry structure is already filled
Peter Vos's avatar
Peter Vos committed
620
            if(m_state != ePtsFilled)
Mike Kirby's avatar
Mike Kirby committed
621
            {
Peter Vos's avatar
Peter Vos committed
622
623
                int i,j,k;
                int nEdgeCoeffs = m_xmap[0]->GetEdgeNcoeffs(0);
Mike Kirby's avatar
Mike Kirby committed
624

Peter Vos's avatar
Peter Vos committed
625
626
                Array<OneD, unsigned int> mapArray (nEdgeCoeffs);
                Array<OneD, int>          signArray(nEdgeCoeffs);
Mike Kirby's avatar
Mike Kirby committed
627

Peter Vos's avatar
Peter Vos committed
628
                for(i = 0; i < kNedges; i++)
Mike Kirby's avatar
Mike Kirby committed
629
                {
Peter Vos's avatar
Peter Vos committed
630
631
632
                    m_edges[i]->FillGeom();
                    m_xmap[0]->GetEdgeToElementMap(i,m_eorient[i],mapArray,signArray);

Spencer Sherwin's avatar
Spencer Sherwin committed
633
634
                    //nEdgeCoeffs = m_xmap[0]->GetEdgeNcoeffs(i);
                    nEdgeCoeffs = (*m_edges[i])[0]->GetNcoeffs();
Peter Vos's avatar
Peter Vos committed
635
636
637
638
639

                    for(j = 0 ; j < m_coordim; j++)
                    {
                        for(k = 0; k < nEdgeCoeffs; k++)
                        {
640
                            (m_xmap[j]->UpdateCoeffs())[ mapArray[k] ]
641
642
                                                         = signArray[k]*
                                                         ((*m_edges[i])[j]->GetCoeffs())[k];
Peter Vos's avatar
Peter Vos committed
643
644
                        }
                    }
Mike Kirby's avatar
Mike Kirby committed
645
                }
646

Peter Vos's avatar
Peter Vos committed
647
                for(i = 0; i < m_coordim; ++i)
Mike Kirby's avatar
Mike Kirby committed
648
                {
Peter Vos's avatar
Peter Vos committed
649
                    m_xmap[i]->BwdTrans(m_xmap[i]->GetCoeffs(),
650
                            m_xmap[i]->UpdatePhys());
Mike Kirby's avatar
Mike Kirby committed
651
                }
652

Peter Vos's avatar
Peter Vos committed
653
                m_state = ePtsFilled;
Mike Kirby's avatar
Mike Kirby committed
654
655
656
            }
        }

657
658
659
660
661

        /**
         *
         */
        void TriGeom::v_GetLocCoords(const Array<OneD,const NekDouble> &coords, Array<OneD,NekDouble> &Lcoords)
Mike Kirby's avatar
Mike Kirby committed
662
        {
663
            TriGeom::v_FillGeom();
Mike Kirby's avatar
Mike Kirby committed
664

665
            // calculate local coordinate for coord
Spencer Sherwin's avatar
Spencer Sherwin committed
666
            if(GetGtype() == eRegular)
667
668
            { 
                NekDouble coords2 = (m_coordim == 3)? coords[2]: 0.0; 
Mike Kirby's avatar
Mike Kirby committed
669
                VertexComponent dv1, dv2, norm, orth1, orth2;
670
                VertexComponent xin(m_coordim,0,coords[0],coords[1],coords2);
Mike Kirby's avatar
Mike Kirby committed
671

672
                // Calculate edge vectors from 0-1 and 0-2 edges. 
Mike Kirby's avatar
Mike Kirby committed
673
674
675
                dv1.Sub(*m_verts[1],*m_verts[0]);
                dv2.Sub(*m_verts[2],*m_verts[0]);

676
                // Obtain normal to plane in which dv1 and dv2 lie
Mike Kirby's avatar
Mike Kirby committed
677
                norm.Mult(dv1,dv2);
678
679
                
                // Obtain vector which are proportional to normal of dv1 and dv2. 
Mike Kirby's avatar
Mike Kirby committed
680
681
                orth1.Mult(norm,dv1);
                orth2.Mult(norm,dv2);
682
683
684
685
686
687
688
689
690
691
                
                // Start with vector of desired points minus vertex_0
                xin -= *m_verts[0];

                // Calculate length using L/|dv1| = (x-v0).n1/(dv1.n1) for coordiante 1
                // Then rescale to [-1,1]. 
                Lcoords[0] = xin.dot(orth2)/dv1.dot(orth2);
                Lcoords[0] = 2*Lcoords[0]-1;
                Lcoords[1] = xin.dot(orth1)/dv2.dot(orth1);
                Lcoords[1] = 2*Lcoords[1]-1;
Mike Kirby's avatar
Mike Kirby committed
692
693
694
            }
            else
            {
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
                // Determine nearest point of coords  to values in m_xmap
                Array<OneD, NekDouble> ptsx = m_xmap[0]->GetPhys();
                Array<OneD, NekDouble> ptsy = m_xmap[1]->GetPhys();
                int npts = ptsx.num_elements();
                Array<OneD, NekDouble> tmpx(npts), tmpy(npts);
                const Array<OneD, const NekDouble> za = m_xmap[0]->GetPoints(0);
                const Array<OneD, const NekDouble> zb = m_xmap[0]->GetPoints(1);
                
                
                //guess the first local coords based on nearest point
                Vmath::Sadd(npts, -coords[0], ptsx,1,tmpx,1);
                Vmath::Sadd(npts, -coords[1], ptsy,1,tmpy,1);
                Vmath::Vmul (npts, tmpx,1,tmpx,1,tmpx,1);
                Vmath::Vvtvp(npts, tmpy,1,tmpy,1,tmpx,1,tmpx,1);
                          
                int min_i = Vmath::Imin(npts,tmpx,1);
                
                Lcoords[0] = za[min_i%za.num_elements()];
                Lcoords[1] = zb[min_i/za.num_elements()];

                // recover cartesian coordinate from collapsed coordinate. 
                Lcoords[0] = (1.0+Lcoords[0])*(1.0-Lcoords[1])/2 -1.0;

                // Perform newton iteration to find local coordinates 
                NewtonIterationForLocCoord(coords,Lcoords);

                
#if 0 // Original verion

Andrea Isoni's avatar
Andrea Isoni committed
724
725
726
727
728
729
730
731
                //NEKERROR(ErrorUtil::efatal,
                //        "inverse mapping must be set up to use this call");

                int i;
                Array<OneD, const NekDouble> pts;

                Array<OneD, NekDouble> ptsx;
                Array<OneD, NekDouble> ptsy;  
732
733
734
735
                Array<OneD, NekDouble> za;
                Array<OneD, NekDouble> zb;  
                Array<OneD, NekDouble> wa;
                Array<OneD, NekDouble> wb;  
Andrea Isoni's avatar
Andrea Isoni committed
736
737
738
                NekDouble xmap,ymap, F1,F2;
                NekDouble jac, derx_1k, derx_2k, dery_1k, dery_2k ;
                NekDouble invderx_1k, invderx_2k, invdery_1k, invdery_2k;
739
740
                F1= F2 = 2000;

Andrea Isoni's avatar
Andrea Isoni committed
741
742
743
744
745
746
747
748

                ptsx = m_xmap[0]->GetPhys();
                ptsy = m_xmap[1]->GetPhys();
                Array<OneD, NekDouble> derx_1 (ptsx.num_elements());
                Array<OneD, NekDouble> derx_2 (ptsx.num_elements());                 
                Array<OneD, NekDouble> dery_1 (ptsy.num_elements());
                Array<OneD, NekDouble> dery_2 (ptsy.num_elements());
                m_xmap[0]->StdPhysDeriv(ptsx, derx_1, derx_2);                  
749
750
751
752
753
                m_xmap[1]->StdPhysDeriv(ptsy, dery_1, dery_2);    

                //guess the first local coords
                //Lcoords[0]=0.0;
                //Lcoords[1]=0.0; 
754
                
755
756
757
758
                boost::shared_ptr<StdRegions::StdTriExp> trimap0 = 
                        boost::dynamic_pointer_cast<StdRegions::StdTriExp>(m_xmap[0]);
                boost::shared_ptr<StdRegions::StdTriExp> trimap1 = 
                        boost::dynamic_pointer_cast<StdRegions::StdTriExp>(m_xmap[1]);
759
                int ic;
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780

                int trinp = m_xmap[0]->GetTotPoints();
                Array<OneD, NekDouble> ltrix(trinp);
                Array<OneD, NekDouble> ltriy(trinp);
                m_xmap[0]->GetCoords(ltrix,ltriy);
        

                NekDouble tmp = 1e10;

                for(i = 0; i < ptsx.num_elements()*ptsy.num_elements(); ++i)
                {
                    Lcoords[0] = (coords[0]-ptsx[i])*(coords[0]-ptsx[i]) 
                               + (coords[1]-ptsy[i])*(coords[1]-ptsy[i]);
                    if(Lcoords[0] < tmp){
                    tmp = Lcoords[0];
                    ic = i;
                    }
                }

                Lcoords[0] = ltrix[ic];
                Lcoords[1] = ltriy[ic];
781
               	 
Andrea Isoni's avatar
Andrea Isoni committed
782
783
784
                //int offset=0;              
                //determine y
                int cnt=0;
785
                while( abs(F2) > 0.00001 || abs(F1)> 0.00001        )
Andrea Isoni's avatar
Andrea Isoni committed
786
                {
787

Andrea Isoni's avatar
Andrea Isoni committed
788
789
790
791
792
                    //calculate the gradient tensor at Lcoords
                    derx_1k = m_xmap[0]->PhysEvaluate(Lcoords, derx_1);
                    derx_2k = m_xmap[0]->PhysEvaluate(Lcoords, derx_2);
                    dery_1k = m_xmap[1]->PhysEvaluate(Lcoords, dery_1);
                    dery_2k = m_xmap[1]->PhysEvaluate(Lcoords, dery_2);                  
793

Andrea Isoni's avatar
Andrea Isoni committed
794
795
                    jac = (derx_1k*dery_2k - derx_2k*dery_1k);//determinant IS the jac

796

Andrea Isoni's avatar
Andrea Isoni committed
797
                    //invert matrix:
798
                    invderx_1k =  dery_2k/jac;
Andrea Isoni's avatar
Andrea Isoni committed
799
800
                    invderx_2k = -derx_2k/jac;
                    invdery_1k = -dery_1k/jac;
801
                    invdery_2k =  derx_1k/jac;
802

Andrea Isoni's avatar
Andrea Isoni committed
803
                    //calculate the global point corresponding to Lcoords
804
805
                    xmap = trimap0->PhysEvaluate(Lcoords, ptsx);
                    ymap = trimap1->PhysEvaluate(Lcoords, ptsy);
Andrea Isoni's avatar
Andrea Isoni committed
806
807
808
809
810
                    Lcoords[0] = Lcoords[0] + invderx_1k*(coords[0]-xmap) + invderx_2k*(coords[1]-ymap);
                    Lcoords[1] = Lcoords[1] + invdery_1k*(coords[0]-xmap) + invdery_2k*(coords[1]-ymap);

                    F1 = coords[0] - xmap;
                    F2 = coords[1] - ymap;
811
     
Andrea Isoni's avatar
Andrea Isoni committed
812
813
814
815
816
817
818
                    cnt++;
         
                    if( cnt >= 21)
                    {   
                             Lcoords[0] = Lcoords[1] = 2.0;    
                             break;
                    }
819
820
821
822
823
824
825
826
827
                }
                //cout<<"it finished"<<endl;
                if(Lcoords[1]>1.01 && Lcoords[0]>1.01)
                {
                    Lcoords[0] = Lcoords[1] = 2.0;    
                }
                
                //cout<<elmtid<<"Locx="<<Lcoords[0]<<"  Locy="<<Lcoords[1]<<endl;
#endif
Mike Kirby's avatar
Mike Kirby committed
828
            }
829
            
Mike Kirby's avatar
Mike Kirby committed
830
        }
831

832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875

        /**
         *
         */
        int TriGeom::v_GetEid(int i) const
        {
            ASSERTL2((i >=0) && (i <= 2),"Edge id must be between 0 and 2");
            return m_edges[i]->GetEid();
        }


        /**
         *
         */
        int TriGeom::v_GetVid(int i) const
        {
            ASSERTL2((i >=0) && (i <= 2),"Vertex id must be between 0 and 2");
            return m_verts[i]->GetVid();
        }


        /**
         *
         */
        const VertexComponentSharedPtr TriGeom::v_GetVertex(int i) const
        {
            ASSERTL2((i >=0) && (i <= 2),"Vertex id must be between 0 and 2");
            return m_verts[i];
        }


        /**
         *
         */
        const Geometry1DSharedPtr TriGeom::v_GetEdge(int i) const
        {
            ASSERTL2((i >=0) && (i <= 2),"Edge id must be between 0 and 3");
            return m_edges[i];
        }


        /**
         *
         */
876
        StdRegions::Orientation TriGeom::v_GetEorient(const int i) const
877
878
879
880
881
882
883
884
885
        {
            ASSERTL2((i >=0) && (i <= 2),"Edge id must be between 0 and 2");
            return m_eorient[i];
        }


        /**
         *
         */
886
        StdRegions::Orientation TriGeom::v_GetCartesianEorient(const int i) const
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
        {
            ASSERTL2((i >=0) && (i <= 3),"Edge id must be between 0 and 3");
            if(i < 2)
            {
                return m_eorient[i];
            }
            else
            {
                if(m_eorient[i] == StdRegions::eForwards)
                {
                    return StdRegions::eBackwards;
                }
                else
                {
                    return StdRegions::eForwards;
                }
            }
        }


        /**
         *
         */
        int TriGeom::v_WhichEdge(SegGeomSharedPtr edge)
        {
            int returnval = -1;

            SegGeomVector::iterator edgeIter;
            int i;

            for (i=0,edgeIter = m_edges.begin(); edgeIter != m_edges.end(); ++edgeIter,++i)
            {
                if (*edgeIter == edge)
                {
                    returnval = i;
                    break;
                }
            }
925
926
927
928

            return returnval;
        }

929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950

        /**
         *
         */
        int TriGeom::v_GetNumVerts() const
        {
            return kNverts;
        }


        /**
         *
         */
        int TriGeom::v_GetNumEdges() const
        {
            return kNedges;
        }


        /**
         *
         */
951
        bool TriGeom::v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord, NekDouble tol)
952
953
        {
            ASSERTL1(gloCoord.num_elements() >= 2,
954
                    "Two dimensional geometry expects at least two coordinates.");
955

956
957
            Array<OneD,NekDouble> stdCoord(GetCoordim(),0.0);
            GetLocCoords(gloCoord, stdCoord);
958
            if (stdCoord[0] >= -(1+tol) && stdCoord[1] >= -(1+tol)
959
                    && stdCoord[0] + stdCoord[1] <= tol)
960
961
962
963
964
            {
                return true;
            }
            return false;
        }
Mike Kirby's avatar
Mike Kirby committed
965
966
    }; //end of namespace
}; //end of namespace