AssemblyMapDG.cpp 39.8 KB
Newer Older
1 2
///////////////////////////////////////////////////////////////////////////////
//
3
// File AssemblyMapDG.cpp
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// 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.
//
// Description: Local to Global Base Class mapping routines
//
///////////////////////////////////////////////////////////////////////////////

36
#include <LibUtilities/BasicUtils/HashUtils.hpp>
37
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
38 39 40 41
#include <MultiRegions/ExpList.h>
#include <LocalRegions/SegExp.h>
#include <LocalRegions/TriExp.h>
#include <LocalRegions/QuadExp.h>
42 43 44 45 46
#include <LocalRegions/HexExp.h>
#include <LocalRegions/TetExp.h>
#include <LocalRegions/PrismExp.h>
#include <LocalRegions/PyrExp.h>
#include <LocalRegions/PointExp.h>
47 48
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
49 50 51 52 53 54 55

#include <boost/config.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/cuthill_mckee_ordering.hpp>
#include <boost/graph/properties.hpp>
#include <boost/graph/bandwidth.hpp>

56
using namespace std;
57 58 59 60 61 62 63 64 65 66 67 68 69 70

namespace Nektar
{
    namespace MultiRegions
    {
        AssemblyMapDG::AssemblyMapDG():
            m_numDirichletBndPhys(0)
        {
        }

        AssemblyMapDG::~AssemblyMapDG()
        {
        }

71
        AssemblyMapDG::AssemblyMapDG(
72
            const LibUtilities::SessionReaderSharedPtr                &pSession,
Dave Moxey's avatar
Dave Moxey committed
73 74
            const SpatialDomains::MeshGraphSharedPtr                  &graph,
            const ExpListSharedPtr                                    &trace,
75
            const ExpList                                             &locExp,
Dave Moxey's avatar
Dave Moxey committed
76 77 78
            const Array<OneD, const MultiRegions::ExpListSharedPtr>         &bndCondExp,
            const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &bndCond,
            const PeriodicMap                                         &periodicTrace,
79 80
            const std::string variable):
            AssemblyMap(pSession,variable)
81
        {
82 83
            int i, j, k, cnt, eid, id, id1, gid;
            int order_e   = 0;
Dave Moxey's avatar
Dave Moxey committed
84
            int nTraceExp = trace->GetExpSize();
85
            int nbnd      = bndCondExp.num_elements();
86

Dave Moxey's avatar
Dave Moxey committed
87 88 89 90 91 92
            LocalRegions::ExpansionSharedPtr  exp;
            LocalRegions::ExpansionSharedPtr  bndExp;
            SpatialDomains::GeometrySharedPtr traceGeom;

            const LocalRegions::ExpansionVector expList = *(locExp.GetExp());
            int nel = expList.size();
93

Dave Moxey's avatar
Dave Moxey committed
94
            map<int, int> meshTraceId;
95 96 97 98

            m_signChange = true;

            // determine mapping from geometry edges to trace
Dave Moxey's avatar
Dave Moxey committed
99
            for(i = 0; i < nTraceExp; ++i)
100
            {
101
                meshTraceId[trace->GetExp(i)->GetGeom()->GetGlobalID()] = i;
102 103
            }

Dave Moxey's avatar
Dave Moxey committed
104
            // Count total number of trace elements
105 106 107
            cnt = 0;
            for(i = 0; i < nel; ++i)
            {
Dave Moxey's avatar
Dave Moxey committed
108
                cnt += expList[i]->GetNtrace();
109 110
            }

Dave Moxey's avatar
Dave Moxey committed
111 112 113
            Array<OneD, LocalRegions::ExpansionSharedPtr> tracemap(cnt);
            m_elmtToTrace = Array<
                OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >(nel);
114

Dave Moxey's avatar
Dave Moxey committed
115
            // set up trace expansions links;
116
            for(cnt = i = 0; i < nel; ++i)
117
            {
Dave Moxey's avatar
Dave Moxey committed
118 119 120
                m_elmtToTrace[i] = tracemap + cnt;

                for(j = 0; j < expList[i]->GetNtrace(); ++j)
121
                {
Dave Moxey's avatar
Dave Moxey committed
122 123 124
                    id = expList[i]->GetGeom()->GetTid(j);

                    if(meshTraceId.count(id) > 0)
125
                    {
Dave Moxey's avatar
Dave Moxey committed
126 127
                        m_elmtToTrace[i][j] =
                            trace->GetExp(meshTraceId.find(id)->second);
128
                    }
129
                    else
130
                    {
Dave Moxey's avatar
Dave Moxey committed
131
                        ASSERTL0(false, "Failed to find trace map");
132 133
                    }
                }
134

Dave Moxey's avatar
Dave Moxey committed
135
                cnt += expList[i]->GetNtrace();
136 137 138 139 140 141 142 143 144
            }

            // Set up boundary mapping
            cnt = 0;
            for(i = 0; i < nbnd; ++i)
            {
                cnt += bndCondExp[i]->GetExpSize();
            }

Dave Moxey's avatar
Dave Moxey committed
145
            set<int> dirTrace;
146

147 148 149 150 151 152 153 154
            m_numLocalDirBndCoeffs = 0;
            m_numDirichletBndPhys  = 0;

            cnt = 0;
            for(i = 0; i < bndCondExp.num_elements(); ++i)
            {
                for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
                {
Dave Moxey's avatar
Dave Moxey committed
155 156 157 158 159
                    bndExp    = bndCondExp[i]->GetExp(j);
                    traceGeom = bndExp->GetGeom();
                    id        = traceGeom->GetGlobalID();

                    if(bndCond[i]->GetBoundaryConditionType() ==
160
                           SpatialDomains::eDirichlet)
161
                    {
Dave Moxey's avatar
Dave Moxey committed
162 163 164
                        m_numLocalDirBndCoeffs += bndExp->GetNcoeffs();
                        m_numDirichletBndPhys  += bndExp->GetTotPoints();
                        dirTrace.insert(id);
165 166
                    }
                }
167

168 169 170
                cnt += j;
            }

Dave Moxey's avatar
Dave Moxey committed
171 172
            // Set up integer mapping array and sign change for each degree of
            // freedom + initialise some more data members.
173
            m_staticCondLevel           = 0;
Dave Moxey's avatar
Dave Moxey committed
174
            m_lowestStaticCondLevel     = 0;
175 176 177 178
            m_numPatches                = nel;
            m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
            m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);

179 180 181
            int nbndry = 0;
            for(i = 0; i < nel; ++i) // count number of elements in array
            {
182
                eid     = i;
Dave Moxey's avatar
Dave Moxey committed
183 184 185 186
                nbndry += expList[eid]->NumDGBndryCoeffs();
                m_numLocalIntCoeffsPerPatch[i] = 0;
                m_numLocalBndCoeffsPerPatch[i] =
                    (unsigned int) expList[eid]->NumDGBndryCoeffs();
187 188 189
            }

            m_numGlobalDirBndCoeffs = m_numLocalDirBndCoeffs;
190 191 192 193
            m_numLocalBndCoeffs     = nbndry;
            m_numLocalCoeffs        = nbndry;
            m_localToGlobalBndMap   = Array<OneD, int>       (nbndry);
            m_localToGlobalBndSign  = Array<OneD, NekDouble> (nbndry,1);
194 195

            // Set up array for potential mesh optimsation
Dave Moxey's avatar
Dave Moxey committed
196
            Array<OneD,int> traceElmtGid(nTraceExp, -1);
197 198 199
            int nDir = 0;
            cnt = 0;

Dave Moxey's avatar
Dave Moxey committed
200 201 202 203
            // We are now going to construct a graph of the mesh which can be
            // reordered depending on the type of solver we would like to use.
            typedef boost::adjacency_list<
                boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
204 205

            BoostGraph boostGraphObj;
Dave Moxey's avatar
Dave Moxey committed
206
            int trace_id, trace_id1;
207
            int dirOffset = 0;
208

Dave Moxey's avatar
Dave Moxey committed
209 210 211
            // make trace trace renumbering map where first solved trace starts
            // at 0 so we can set up graph.
            for(i = 0; i < nTraceExp; ++i)
212
            {
Dave Moxey's avatar
Dave Moxey committed
213 214 215
                id = trace->GetExp(i)->GetGeom()->GetGlobalID();

                if (dirTrace.count(id) == 0)
216
                {
Dave Moxey's avatar
Dave Moxey committed
217 218
                    // Initial put in element ordering (starting from zero) into
                    // traceElmtGid
219
                    boost::add_vertex(boostGraphObj);
Dave Moxey's avatar
Dave Moxey committed
220
                    traceElmtGid[i] = cnt++;
221 222 223 224
                }
                else
                {
                    // Use existing offset for Dirichlet edges
Dave Moxey's avatar
Dave Moxey committed
225 226
                    traceElmtGid[i] = dirOffset;
                    dirOffset      += trace->GetExp(i)->GetNcoeffs();
227 228 229 230 231 232 233
                    nDir++;
                }
            }

            // Set up boost Graph
            for(i = 0; i < nel; ++i)
            {
234
                eid = i;
235

Dave Moxey's avatar
Dave Moxey committed
236
                for(j = 0; j < expList[eid]->GetNtrace(); ++j)
237
                {
Dave Moxey's avatar
Dave Moxey committed
238 239 240 241
                    // Add trace to boost graph for non-Dirichlet Boundary
                    traceGeom = m_elmtToTrace[eid][j]->GetGeom();
                    id        = traceGeom->GetGlobalID();
                    trace_id  = meshTraceId.find(id)->second;
242

Dave Moxey's avatar
Dave Moxey committed
243
                    if(dirTrace.count(id) == 0)
244
                    {
Dave Moxey's avatar
Dave Moxey committed
245
                        for(k = j+1; k < expList[eid]->GetNtrace(); ++k)
246
                        {
Dave Moxey's avatar
Dave Moxey committed
247 248 249 250 251
                            traceGeom = m_elmtToTrace[eid][k]->GetGeom();
                            id1       = traceGeom->GetGlobalID();
                            trace_id1 = meshTraceId.find(id1)->second;

                            if(dirTrace.count(id1) == 0)
252
                            {
Dave Moxey's avatar
Dave Moxey committed
253 254
                                boost::add_edge((size_t)traceElmtGid[trace_id],
                                                (size_t)traceElmtGid[trace_id1],
255
                                                boostGraphObj);
256 257 258 259 260 261
                            }
                        }
                    }
                }
            }

Dave Moxey's avatar
Dave Moxey committed
262
            int                                 nGraphVerts = nTraceExp - nDir;
263 264 265
            Array<OneD, int>                    perm (nGraphVerts);
            Array<OneD, int>                    iperm(nGraphVerts);
            Array<OneD, int>                    vwgts(nGraphVerts);
266
            BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
Dave Moxey's avatar
Dave Moxey committed
267

268 269 270 271 272 273 274 275 276
            for(i = 0; i < nGraphVerts; ++i)
            {
                vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
            }

            if(nGraphVerts)
            {
                switch(m_solnType)
                {
277 278 279
                    case eDirectFullMatrix:
                    case eIterativeFull:
                    case eIterativeStaticCond:
Dave Moxey's avatar
Dave Moxey committed
280 281
                    case eXxtFullMatrix:
                    case eXxtStaticCond:
Dave Moxey's avatar
Dave Moxey committed
282 283
                    case ePETScFullMatrix:
                    case ePETScStaticCond:
284 285
                    {
                        NoReordering(boostGraphObj,perm,iperm);
286
                        break;
287
                    }
288
                    case eDirectStaticCond:
289 290
                    {
                        CuthillMckeeReordering(boostGraphObj,perm,iperm);
291
                        break;
292
                    }
293
                    case eDirectMultiLevelStaticCond:
294
                    case eIterativeMultiLevelStaticCond:
Dave Moxey's avatar
Dave Moxey committed
295 296
                    case eXxtMultiLevelStaticCond:
                    case ePETScMultiLevelStaticCond:
297
                    {
298 299
                        MultiLevelBisectionReordering(boostGraphObj,perm,iperm,
                                                      bottomUpGraph);
300
                        break;
301
                    }
302
                    default:
303 304 305 306 307 308
                    {
                        ASSERTL0(false,"Unrecognised solution type");
                    }
                }
            }

Dave Moxey's avatar
Dave Moxey committed
309 310
            // Recast the permutation so that it can be used as a map from old
            // trace ID to new trace ID
311
            cnt = m_numLocalDirBndCoeffs;
Dave Moxey's avatar
Dave Moxey committed
312
            for(i = 0; i < nTraceExp - nDir; ++i)
313
            {
Dave Moxey's avatar
Dave Moxey committed
314
                traceElmtGid[perm[i]+nDir] = cnt;
315 316 317 318
                cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
            }

            // Now have trace edges Gid position
Dave Moxey's avatar
Dave Moxey committed
319

320 321 322
            cnt = 0;
            for(i = 0; i < nel; ++i)
            {
323
                eid = i;
Dave Moxey's avatar
Dave Moxey committed
324
                exp = expList[eid];
325

Dave Moxey's avatar
Dave Moxey committed
326
                for(j = 0; j < exp->GetNtrace(); ++j)
327
                {
Dave Moxey's avatar
Dave Moxey committed
328 329 330 331 332 333 334 335
                    traceGeom = m_elmtToTrace[eid][j]->GetGeom();
                    id        = traceGeom->GetGlobalID();
                    gid       = traceElmtGid[meshTraceId.find(id)->second];

                    const int nDim = expList[eid]->GetNumBases();

                    if (nDim == 1)
                    {
Dave Moxey's avatar
Dave Moxey committed
336 337
                        order_e = 1;
                        m_localToGlobalBndMap[cnt] = gid;
Dave Moxey's avatar
Dave Moxey committed
338 339
                    }
                    else if (nDim == 2)
340
                    {
Dave Moxey's avatar
Dave Moxey committed
341 342 343 344 345 346 347 348 349 350
                        order_e = expList[eid]->GetEdgeNcoeffs(j);
                    
                        if(expList[eid]->GetEorient(j) == StdRegions::eForwards)
                        {
                            for(k = 0; k < order_e; ++k)
                            {
                                m_localToGlobalBndMap[k+cnt] = gid + k;
                            }
                        }
                        else
351
                        {
Dave Moxey's avatar
Dave Moxey committed
352
                            switch(m_elmtToTrace[eid][j]->GetBasisType(0))
353
                            {
Dave Moxey's avatar
Dave Moxey committed
354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
                                case LibUtilities::eModified_A:
                                {
                                    // reverse vertex order
                                    m_localToGlobalBndMap[cnt]   = gid + 1;
                                    m_localToGlobalBndMap[cnt+1] = gid;
                                    for (k = 2; k < order_e; ++k)
                                    {
                                        m_localToGlobalBndMap[k+cnt] = gid + k;
                                    }

                                    // negate odd modes
                                    for(k = 3; k < order_e; k+=2)
                                    {
                                        m_localToGlobalBndSign[cnt+k] = -1.0;
                                    }
                                    break;
                                }
                                case LibUtilities::eGLL_Lagrange:
                                {
                                    // reverse  order
                                    for(k = 0; k < order_e; ++k)
                                    {
                                        m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
                                    }
                                    break;
                                }
                                case LibUtilities::eGauss_Lagrange:
                                {
                                    // reverse  order
                                    for(k = 0; k < order_e; ++k)
                                    {
                                        m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
                                    }
                                    break;
                                }
                                default:
                                {
                                    ASSERTL0(false,"Boundary type not permitted");
                                }
393 394 395
                            }
                        }
                    }
Dave Moxey's avatar
Dave Moxey committed
396
                    else if (nDim == 3)
397
                    {
Dave Moxey's avatar
Dave Moxey committed
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 434 435 436 437
                        order_e = expList[eid]->GetFaceNcoeffs(j);

                        std::map<int, int> orientMap;

                        Array<OneD, unsigned int> elmMap1 (order_e);
                        Array<OneD,          int> elmSign1(order_e);
                        Array<OneD, unsigned int> elmMap2 (order_e);
                        Array<OneD,          int> elmSign2(order_e);

                        StdRegions::Orientation fo = expList[eid]->GetForient(j);

                        // Construct mapping which will permute global IDs
                        // according to face orientations.
                        expList[eid]->GetFaceToElementMap(j,fo,elmMap1,elmSign1);
                        expList[eid]->GetFaceToElementMap(
                            j,StdRegions::eDir1FwdDir1_Dir2FwdDir2,elmMap2,elmSign2);

                        for (k = 0; k < elmMap1.num_elements(); ++k)
                        {
                            // Find the elemental co-efficient in the original
                            // mapping.
                            int idx = -1;
                            for (int l = 0; l < elmMap2.num_elements(); ++l)
                            {
                                if (elmMap1[k] == elmMap2[l])
                                {
                                    idx = l;
                                    break;
                                }
                            }

                            ASSERTL2(idx != -1, "Problem with face to element map!");
                            orientMap[k] = idx;
                        }

                        for(k = 0; k < order_e; ++k)
                        {
                            m_localToGlobalBndMap [k+cnt] = gid + orientMap[k];
                            m_localToGlobalBndSign[k+cnt] = elmSign2[orientMap[k]];
                        }
438 439 440 441 442 443 444 445 446 447 448 449 450
                    }

                    cnt += order_e;
                }
            }

            // set up m_bndCondCoeffsToGlobalCoeffsMap to align with map
            cnt = 0;
            for(i = 0; i < nbnd; ++i)
            {
                cnt += bndCondExp[i]->GetNcoeffs();
            }

Dave Moxey's avatar
Dave Moxey committed
451
            m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int>(cnt);
452 453

            // Number of boundary expansions
454
            int nbndexp = 0, bndOffset, bndTotal = 0;
455 456 457 458
            for(cnt = i = 0; i < nbnd; ++i)
            {
                for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
                {
Dave Moxey's avatar
Dave Moxey committed
459 460 461
                    bndExp    = bndCondExp[i]->GetExp(j);
                    id        = bndExp->GetGeom()->GetGlobalID();
                    gid       = traceElmtGid[meshTraceId.find(id)->second];
462
                    bndOffset = bndCondExp[i]->GetCoeff_Offset(j) + bndTotal;
Dave Moxey's avatar
Dave Moxey committed
463

464 465 466
                    // Since boundary information is defined to be aligned with
                    // the geometry just use forward/forward (both coordinate
                    // directions) defintiion for gid's
Dave Moxey's avatar
Dave Moxey committed
467
                    for(k = 0; k < bndExp->GetNcoeffs(); ++k)
468
                    {
469
                        m_bndCondCoeffsToGlobalCoeffsMap[bndOffset+k] = gid + k;
470 471
                    }
                }
472

473 474
                nbndexp  += bndCondExp[i]->GetExpSize();
                bndTotal += bndCondExp[i]->GetNcoeffs();
475
            }
Dave Moxey's avatar
Dave Moxey committed
476

477
            m_numGlobalBndCoeffs = trace->GetNcoeffs();
478
            m_numGlobalCoeffs    = m_numGlobalBndCoeffs;
479 480 481 482

            CalculateBndSystemBandWidth();

            cnt = 0;
483
            m_bndCondTraceToGlobalTraceMap = Array<OneD, int>(nbndexp);
484 485 486 487
            for(i = 0; i < bndCondExp.num_elements(); ++i)
            {
                for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
                {
Dave Moxey's avatar
Dave Moxey committed
488 489 490 491 492
                    bndExp    = bndCondExp[i]->GetExp(j);
                    traceGeom = bndExp->GetGeom();
                    id        = traceGeom->GetGlobalID();
                    m_bndCondTraceToGlobalTraceMap[cnt++] =
                        meshTraceId.find(id)->second;
493 494 495 496
                }
            }

            // Now set up mapping from global coefficients to universal.
497
            ExpListSharedPtr tr = std::dynamic_pointer_cast<ExpList>(trace);
498
            SetUpUniversalDGMap   (locExp);
Dave Moxey's avatar
Dave Moxey committed
499
            SetUpUniversalTraceMap(locExp, tr, periodicTrace);
500

Dave Moxey's avatar
Dave Moxey committed
501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520
            if ((m_solnType == eDirectMultiLevelStaticCond ||
                 m_solnType == eIterativeMultiLevelStaticCond ||
                 m_solnType == eXxtMultiLevelStaticCond ||
                 m_solnType == ePETScMultiLevelStaticCond)
                && nGraphVerts)
            {
                if (m_staticCondLevel < (bottomUpGraph->GetNlevels() - 1))
                {
                    Array<OneD, int> vwgts_perm(nGraphVerts);

                    for (int i = 0; i < nGraphVerts; i++)
                    {
                        vwgts_perm[i] = vwgts[perm[i]];
                    }

                    bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
                    m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
                        AllocateSharedPtr(this, bottomUpGraph);
                }
            }
Dave Moxey's avatar
Dave Moxey committed
521

522 523
            m_hash = hash_range(m_localToGlobalBndMap.begin(),
                                m_localToGlobalBndMap.end());
524 525 526 527 528 529 530
        }

        /**
         * Constructs a mapping between the process-local global numbering and
         * a universal numbering of the trace space expansion. The universal
         * numbering is defined by the mesh edge IDs to enforce consistency
         * across processes.
Dave Moxey's avatar
Dave Moxey committed
531
         *
532 533 534 535
         * @param       locExp  List of local elemental expansions.
         */
        void AssemblyMapDG::SetUpUniversalDGMap(const ExpList &locExp)
        {
Dave Moxey's avatar
Dave Moxey committed
536
            LocalRegions::ExpansionSharedPtr locExpansion;
537 538 539 540
            int eid       = 0;
            int cnt       = 0;
            int id        = 0;
            int order_e   = 0;
541
            int vGlobalId = 0;
542 543 544 545
            int maxDof    = 0;
            int dof       = 0;
            int nDim      = 0;
            int i,j,k;
546

547
            const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
548 549 550 551 552

            // Initialise the global to universal maps.
            m_globalToUniversalBndMap = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
            m_globalToUniversalBndMapUnique = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);

Dave Moxey's avatar
Dave Moxey committed
553
            // Loop over all the elements in the domain and compute max
554 555 556
            // DOF. Reduce across all processes to get universal maximum.
            for(i = 0; i < locExpVector.size(); ++i)
            {
Dave Moxey's avatar
Dave Moxey committed
557
                locExpansion = locExpVector[i];
558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
                nDim = locExpansion->GetShapeDimension();

                // Loop over all edges of element i
                if (nDim == 1)
                {
                    maxDof = (1 > maxDof ? 1 : maxDof);
                }
                else if (nDim == 2)
                {
                    for (j = 0; j < locExpansion->GetNedges(); ++j)
                    {
                        dof    = locExpansion->GetEdgeNcoeffs(j);
                        maxDof = (dof > maxDof ? dof : maxDof);
                    }
                }
                else if (nDim == 3)
                {
                    for (j = 0; j < locExpansion->GetNfaces(); ++j)
                    {
                        dof    = locExpansion->GetFaceNcoeffs(j);
                        maxDof = (dof > maxDof ? dof : maxDof);
                    }
                }
            }
            m_comm->AllReduce(maxDof, LibUtilities::ReduceMax);

            // Now have trace edges Gid position
            cnt = 0;
            for(i = 0; i < locExpVector.size(); ++i)
            {
588
                eid = i;
Douglas Serson's avatar
Douglas Serson committed
589 590
                locExpansion = locExpVector[eid];
                nDim = locExpansion->GetShapeDimension();
591 592 593 594

                // Populate mapping for each edge of the element.
                if (nDim == 1)
                {
595 596
                    int nverts = locExpansion->GetNverts();
                    for(j = 0; j < nverts; ++j)
597
                    {
Dave Moxey's avatar
Dave Moxey committed
598
                        LocalRegions::PointExpSharedPtr locPointExp =
599
                            m_elmtToTrace[eid][j]->as<LocalRegions::PointExp>();
600
                        id = locPointExp->GetGeom()->GetGlobalID();
601 602 603 604
                        vGlobalId = m_localToGlobalBndMap[cnt+j];
                        m_globalToUniversalBndMap[vGlobalId]
                            = id * maxDof + j + 1;
                    }
605
                    cnt += nverts;
Dave Moxey's avatar
Dave Moxey committed
606
                }
607 608 609 610
                else if (nDim == 2)
                {
                    for(j = 0; j < locExpansion->GetNedges(); ++j)
                    {
Dave Moxey's avatar
Dave Moxey committed
611
                        LocalRegions::SegExpSharedPtr locSegExp =
612
                            m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
613 614

                        id  = locSegExp->GetGeom1D()->GetEid();
615
                        order_e = locExpansion->GetEdgeNcoeffs(j);
Dave Moxey's avatar
Dave Moxey committed
616

617 618 619
                        map<int,int> orientMap;
                        Array<OneD, unsigned int> map1(order_e), map2(order_e);
                        Array<OneD, int> sign1(order_e), sign2(order_e);
Dave Moxey's avatar
Dave Moxey committed
620

621 622
                        locExpansion->GetEdgeToElementMap(j, StdRegions::eForwards, map1, sign1);
                        locExpansion->GetEdgeToElementMap(j, locExpansion->GetEorient(j), map2, sign2);
Dave Moxey's avatar
Dave Moxey committed
623

624 625 626 627 628 629 630 631 632 633 634 635 636
                        for (k = 0; k < map1.num_elements(); ++k)
                        {
                            // Find the elemental co-efficient in the original
                            // mapping.
                            int idx = -1;
                            for (int l = 0; l < map2.num_elements(); ++l)
                            {
                                if (map1[k] == map2[l])
                                {
                                    idx = l;
                                    break;
                                }
                            }
Dave Moxey's avatar
Dave Moxey committed
637

638 639 640
                            ASSERTL2(idx != -1, "Problem with face to element map!");
                            orientMap[k] = idx;
                        }
641 642 643 644 645

                        for(k = 0; k < order_e; ++k)
                        {
                            vGlobalId = m_localToGlobalBndMap[k+cnt];
                            m_globalToUniversalBndMap[vGlobalId]
646
                                = id * maxDof + orientMap[k] + 1;
647 648 649 650 651 652 653 654
                        }
                        cnt += order_e;
                    }
                }
                else if (nDim == 3)
                {
                    for(j = 0; j < locExpansion->GetNfaces(); ++j)
                    {
Dave Moxey's avatar
Dave Moxey committed
655
                        LocalRegions::Expansion2DSharedPtr locFaceExp =
656 657
                                m_elmtToTrace[eid][j]
                                           ->as<LocalRegions::Expansion2D>();
658 659

                        id  = locFaceExp->GetGeom2D()->GetFid();
660
                        order_e = locExpansion->GetFaceNcoeffs(j);
661

662 663 664
                        map<int,int> orientMap;
                        Array<OneD, unsigned int> map1(order_e), map2(order_e);
                        Array<OneD, int> sign1(order_e), sign2(order_e);
Dave Moxey's avatar
Dave Moxey committed
665

666 667
                        locExpansion->GetFaceToElementMap(j, StdRegions::eDir1FwdDir1_Dir2FwdDir2, map1, sign1);
                        locExpansion->GetFaceToElementMap(j, locExpansion->GetForient(j), map2, sign2);
Dave Moxey's avatar
Dave Moxey committed
668

669 670 671 672 673 674 675 676 677 678 679 680 681
                        for (k = 0; k < map1.num_elements(); ++k)
                        {
                            // Find the elemental co-efficient in the original
                            // mapping.
                            int idx = -1;
                            for (int l = 0; l < map2.num_elements(); ++l)
                            {
                                if (map1[k] == map2[l])
                                {
                                    idx = l;
                                    break;
                                }
                            }
Dave Moxey's avatar
Dave Moxey committed
682

683 684 685 686
                            ASSERTL2(idx != -1, "Problem with face to element map!");
                            orientMap[k] = idx;
                        }

687 688 689 690
                        for(k = 0; k < order_e; ++k)
                        {
                            vGlobalId = m_localToGlobalBndMap[k+cnt];
                            m_globalToUniversalBndMap[vGlobalId]
691
                                = id * maxDof + orientMap[k] + 1;
692 693 694 695 696
                        }
                        cnt += order_e;
                    }
                }
            }
Dave Moxey's avatar
Dave Moxey committed
697

698 699 700 701 702 703 704 705 706 707 708 709
            // Initialise GSlib and populate the unique map.
            Array<OneD, long> tmp(m_globalToUniversalBndMap.num_elements());
            for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
            {
                tmp[i] = m_globalToUniversalBndMap[i];
            }
            m_bndGsh = m_gsh = Gs::Init(tmp, m_comm);
            Gs::Unique(tmp, m_comm);
            for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
            {
                m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
            }
710 711
        }

712 713 714 715
        void AssemblyMapDG::SetUpUniversalTraceMap(
            const ExpList         &locExp,
            const ExpListSharedPtr trace,
            const PeriodicMap     &perMap)
716
        {
717
            Array<OneD, int> tmp;
Dave Moxey's avatar
Dave Moxey committed
718
            LocalRegions::ExpansionSharedPtr locExpansion;
719
            int i;
720 721
            int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;

722
            const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
723 724 725 726

            int nTracePhys = trace->GetTotPoints();

            // Initialise the trace to universal maps.
727
            m_traceToUniversalMap       =
728
                Nektar::Array<OneD, int>(nTracePhys, -1);
729
            m_traceToUniversalMapUnique =
730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756
                Nektar::Array<OneD, int>(nTracePhys, -1);

            // Assume that each element of the expansion is of the same
            // dimension.
            nDim = locExpVector[0]->GetShapeDimension();

            if (nDim == 1)
            {
                maxQuad = (1 > maxQuad ? 1 : maxQuad);
            }
            else
            {
                for (i = 0; i < trace->GetExpSize(); ++i)
                {
                    quad = trace->GetExp(i)->GetTotPoints();
                    if (quad > maxQuad)
                    {
                        maxQuad = quad;
                    }
                }
            }
            m_comm->AllReduce(maxQuad, LibUtilities::ReduceMax);

            if (nDim == 1)
            {
                for (int i = 0; i < trace->GetExpSize(); ++i)
                {
757
                    eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
758
                    offset = trace->GetPhys_Offset(i);
759 760

                    // Check to see if this vert is periodic. If it is, then we
Dave Moxey's avatar
Dave Moxey committed
761
                    // need use the unique eid of the two points
762
                    auto it = perMap.find(eid);
763 764 765 766 767 768 769 770 771
                    if (perMap.count(eid) > 0)
                    {
                        PeriodicEntity ent = it->second[0];
                        if (ent.isLocal == false) // Not sure if true in 1D
                        {
                            eid = min(eid, ent.id);
                        }
                    }

772 773 774
                    m_traceToUniversalMap[offset] = eid*maxQuad+1;
                }
            }
775
            else
776 777 778
            {
                for (int i = 0; i < trace->GetExpSize(); ++i)
                {
779
                    eid    = trace->GetExp(i)->GetGeom()->GetGlobalID();
780 781 782
                    offset = trace->GetPhys_Offset(i);
                    quad   = trace->GetExp(i)->GetTotPoints();

783 784 785 786
                    // Check to see if this edge is periodic. If it is, then we
                    // need to reverse the trace order of one edge only in the
                    // universal map so that the data are reversed w.r.t each
                    // other. We do this by using the minimum of the two IDs.
787
                    auto it = perMap.find(eid);
788
                    bool realign = false;
789 790
                    if (perMap.count(eid) > 0)
                    {
791
                        PeriodicEntity ent = it->second[0];
792
                        if (ent.isLocal == false)
793
                        {
794
                            realign = eid == min(eid, ent.id);
795
                            eid = min(eid, ent.id);
796 797 798
                        }
                    }

799
                    for (int j = 0; j < quad; ++j)
800 801 802
                    {
                        m_traceToUniversalMap[j+offset] = eid*maxQuad+j+1;
                    }
803 804

                    if (realign)
805
                    {
806 807
                        if (nDim == 2)
                        {
808 809 810
                            RealignTraceElement(
                                tmp = m_traceToUniversalMap+offset,
                                it->second[0].orient, quad);
811 812
                        }
                        else
813
                        {
814 815 816 817 818
                            RealignTraceElement(
                                tmp = m_traceToUniversalMap+offset,
                                it->second[0].orient,
                                trace->GetExp(i)->GetNumPoints(0),
                                trace->GetExp(i)->GetNumPoints(1));
819
                        }
820 821 822 823
                    }
                }
            }

824
            Array<OneD, long> tmp2(nTracePhys);
825 826
            for (int i = 0; i < nTracePhys; ++i)
            {
827
                tmp2[i] = m_traceToUniversalMap[i];
828
            }
829 830
            m_traceGsh = Gs::Init(tmp2, m_comm);
            Gs::Unique(tmp2, m_comm);
831 832
            for (int i = 0; i < nTracePhys; ++i)
            {
833
                m_traceToUniversalMapUnique[i] = tmp2[i];
834 835 836
            }
        }

837
        void AssemblyMapDG::RealignTraceElement(
838 839 840 841
            Array<OneD, int>        &toAlign,
            StdRegions::Orientation  orient,
            int                      nquad1,
            int                      nquad2)
842 843 844 845 846 847 848 849 850
        {
            if (orient == StdRegions::eBackwards)
            {
                ASSERTL1(nquad2 == 0, "nquad2 != 0 for reorienation");
                for (int i = 0; i < nquad1/2; ++i)
                {
                    swap(toAlign[i], toAlign[nquad1-1-i]);
                }
            }
851
            else if (orient != StdRegions::eForwards)
852 853 854 855 856
            {
                ASSERTL1(nquad2 != 0, "nquad2 == 0 for reorienation");

                Array<OneD, int> tmp(nquad1*nquad2);

857
                // Copy transpose.
858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 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
                if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
                    orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                    orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                    orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                {
                    for (int i = 0; i < nquad2; ++i)
                    {
                        for (int j = 0; j < nquad1; ++j)
                        {
                            tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
                        }
                    }
                }
                else
                {
                    for (int i = 0; i < nquad2; ++i)
                    {
                        for (int j = 0; j < nquad1; ++j)
                        {
                            tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
                        }
                    }
                }

                if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
                    orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                    orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
                    orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                {
                    // Reverse x direction
                    for (int i = 0; i < nquad2; ++i)
                    {
                        for (int j = 0; j < nquad1/2; ++j)
                        {
                            swap(tmp[i*nquad1 + j],
                                 tmp[i*nquad1 + nquad1-j-1]);
                        }
                    }
                }

                if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
                    orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
                    orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
                    orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
                {
                    // Reverse y direction
                    for (int j = 0; j < nquad1; ++j)
                    {
                        for (int i = 0; i < nquad2/2; ++i)
                        {
                            swap(tmp[i*nquad1 + j],
                                 tmp[(nquad2-i-1)*nquad1 + j]);
                        }
                    }
                }
                Vmath::Vcopy(nquad1*nquad2, tmp, 1, toAlign, 1);
914 915 916 917 918 919
            }
        }

        void AssemblyMapDG::UniversalTraceAssemble(
            Array<OneD, NekDouble> &pGlobal) const
        {
920
            Gs::Gather(pGlobal, Gs::gs_add, m_traceGsh);
921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958
        }

        int AssemblyMapDG::v_GetLocalToGlobalMap(const int i) const
        {
            return m_localToGlobalBndMap[i];
        }

        int AssemblyMapDG::v_GetGlobalToUniversalMap(const int i) const
        {
            return m_globalToUniversalBndMap[i];
        }

        int AssemblyMapDG::v_GetGlobalToUniversalMapUnique(const int i) const
        {
            return m_globalToUniversalBndMapUnique[i];
        }

        const Array<OneD,const int>& AssemblyMapDG::v_GetLocalToGlobalMap()
        {
            return m_localToGlobalBndMap;
        }

        const Array<OneD,const int>& AssemblyMapDG::v_GetGlobalToUniversalMap()
        {
            return m_globalToUniversalBndMap;
        }

        const Array<OneD,const int>& AssemblyMapDG::v_GetGlobalToUniversalMapUnique()
        {
            return m_globalToUniversalBndMapUnique;
        }

        NekDouble AssemblyMapDG::v_GetLocalToGlobalSign(
                    const int i) const
        {
            return GetLocalToGlobalBndSign(i);
        }

959
        void AssemblyMapDG::v_LocalToGlobal(
960
                    const Array<OneD, const NekDouble>& loc,
961 962
                    Array<OneD,       NekDouble>& global,
                    bool useComm ) const
963 964 965 966
        {
            AssembleBnd(loc,global);
        }

967
        void AssemblyMapDG::v_LocalToGlobal(
968
                    const NekVector<NekDouble>& loc,
969 970
                    NekVector<      NekDouble>& global,
                    bool useComm) const
971 972 973 974
        {
            AssembleBnd(loc,global);
        }

975
        void AssemblyMapDG::v_GlobalToLocal(
976 977 978 979 980 981
                    const Array<OneD, const NekDouble>& global,
                          Array<OneD,       NekDouble>& loc) const
        {
            GlobalToLocalBnd(global,loc);
        }

982
        void AssemblyMapDG::v_GlobalToLocal(
983 984 985 986 987 988
                    const NekVector<NekDouble>& global,
                          NekVector<      NekDouble>& loc) const
        {
            GlobalToLocalBnd(global,loc);
        }

989
        void AssemblyMapDG::v_Assemble(
990 991 992 993 994 995
                    const Array<OneD, const NekDouble> &loc,
                          Array<OneD,       NekDouble> &global) const
        {
            AssembleBnd(loc,global);
        }

996
        void AssemblyMapDG::v_Assemble(
997 998 999 1000 1001 1002
                    const NekVector<NekDouble>& loc,
                          NekVector<      NekDouble>& global) const
        {
            AssembleBnd(loc,global);
        }

1003
        void AssemblyMapDG::v_UniversalAssemble(
1004 1005 1006 1007 1008
                      Array<OneD,     NekDouble>& pGlobal) const
        {
            Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
        }

1009
        void AssemblyMapDG::v_UniversalAssemble(
1010 1011 1012 1013 1014
                      NekVector<      NekDouble>& pGlobal) const
        {
            UniversalAssemble(pGlobal.GetPtr());
        }

1015
        int AssemblyMapDG::v_GetFullSystemBandWidth() const
1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034
        {
            return GetBndSystemBandWidth();
        }

        int AssemblyMapDG::GetTraceToUniversalMap(int i)
        {
            return m_traceToUniversalMap[i];
        }

        int AssemblyMapDG::GetTraceToUniversalMapUnique(int i)
        {
            return m_traceToUniversalMapUnique[i];
        }

        int AssemblyMapDG::GetNumDirichletBndPhys()
        {
            return m_numDirichletBndPhys;
        }

Dave Moxey's avatar
Dave Moxey committed
1035
        Array<OneD, LocalRegions::ExpansionSharedPtr>&
1036 1037 1038 1039 1040 1041 1042
                    AssemblyMapDG::GetElmtToTrace(const int i)
        {
            ASSERTL1(i >= 0 && i < m_elmtToTrace.num_elements(),
                     "i is out of range");
            return m_elmtToTrace[i];
        }

Dave Moxey's avatar
Dave Moxey committed
1043
        Array<OneD, Array< OneD, LocalRegions::ExpansionSharedPtr> >&
1044 1045 1046 1047 1048 1049
                    AssemblyMapDG::GetElmtToTrace()
        {
            return m_elmtToTrace;
        }
    } //namespace
} // namespace