TetGeom.cpp 26.3 KB
Newer Older
Mike Kirby's avatar
Mike Kirby committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
//  File: TetGeom.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: Tetrahedral geometry information.
Mike Kirby's avatar
Mike Kirby committed
33 34 35 36 37
//
////////////////////////////////////////////////////////////////////////////////

#include <SpatialDomains/TetGeom.h>

38 39 40 41
#include <SpatialDomains/Geometry1D.h>
#include <StdRegions/StdTetExp.h>
#include <SpatialDomains/SegGeom.h>

42 43
using namespace std;

Mike Kirby's avatar
Mike Kirby committed
44 45
namespace Nektar
{
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
namespace SpatialDomains
{
const unsigned int TetGeom::VertexEdgeConnectivity[4][3] = {
    {0, 2, 3}, {0, 1, 4}, {1, 2, 5}, {3, 4, 5}};
const unsigned int TetGeom::VertexFaceConnectivity[4][3] = {
    {0, 1, 3}, {0, 1, 2}, {0, 2, 3}, {1, 2, 3}};
const unsigned int TetGeom::EdgeFaceConnectivity[6][2] = {
    {0, 1}, {0, 2}, {0, 3}, {1, 3}, {1, 2}, {2, 3}};

TetGeom::TetGeom()
{
    m_shapeType = LibUtilities::eTetrahedron;
}

TetGeom::TetGeom(int id, const TriGeomSharedPtr faces[])
    : Geometry3D(faces[0]->GetEdge(0)->GetVertex(0)->GetCoordim())
{
    m_shapeType = LibUtilities::eTetrahedron;
    m_globalID = id;

    /// Copy the face shared pointers
    m_faces.insert(m_faces.begin(), faces, faces + TetGeom::kNfaces);

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

    SetUpLocalEdges();
    SetUpLocalVertices();
    SetUpEdgeOrientation();
    SetUpFaceOrientation();
}

TetGeom::~TetGeom()
{
}

/**
 * @brief Determines if a point specified in global coordinates is
 * located within this tetrahedral geometry and return local caretsian
 * coordinates
 */
bool TetGeom::v_ContainsPoint(const Array<OneD, const NekDouble> &gloCoord,
                              Array<OneD, NekDouble> &locCoord,
                              NekDouble tol,
                              NekDouble &resid)
{
    // Validation checks
    ASSERTL1(gloCoord.num_elements() == 3,
             "Three dimensional geometry expects three coordinates.");

    // find min, max point and check if within twice this
    // distance other false this is advisable since
    // GetLocCoord is expensive for non regular elements.
    if (GetMetricInfo()->GetGtype() != eRegular)
Mike Kirby's avatar
Mike Kirby committed
101
    {
102 103 104
        int i;
        Array<OneD, NekDouble> mincoord(3), maxcoord(3);
        NekDouble diff = 0.0;
105

106
        v_FillGeom();
107

108 109
        const int npts = m_xmap->GetTotPoints();
        Array<OneD, NekDouble> pts(npts);
110

111 112 113 114 115 116
        for (i = 0; i < 3; ++i)
        {
            m_xmap->BwdTrans(m_coeffs[i], pts);

            mincoord[i] = Vmath::Vmin(pts.num_elements(), pts, 1);
            maxcoord[i] = Vmath::Vmax(pts.num_elements(), pts, 1);
117

118
            diff = max(maxcoord[i] - mincoord[i], diff);
119
        }
120

121
        for (i = 0; i < 3; ++i)
122
        {
123 124 125 126 127
            if ((gloCoord[i] < mincoord[i] - 0.2 * diff) ||
                (gloCoord[i] > maxcoord[i] + 0.2 * diff))
            {
                return false;
            }
128
        }
129 130 131 132
    }

    // Convert to the local (eta) coordinates.
    resid = v_GetLocCoords(gloCoord, locCoord);
133

134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
    // Check local coordinate is within cartesian bounds.
    if (locCoord[0] >= -(1 + tol) && locCoord[1] >= -(1 + tol) &&
        locCoord[2] >= -(1 + tol) &&
        locCoord[0] + locCoord[1] + locCoord[2] <= -1 + tol)
    {
        return true;
    }

    // If out of range clamp locCoord to be within [-1,1]^3
    // since any larger value will be very oscillatory if
    // called by 'returnNearestElmt' option in
    // ExpList::GetExpIndex
    for (int i = 0; i < 3; ++i)
    {
        if (locCoord[i] < -(1 + tol))
149
        {
150
            locCoord[i] = -(1 + tol);
151 152
        }

153
        if (locCoord[i] > (1 + tol))
154
        {
155
            locCoord[i] = 1 + tol;
156
        }
157 158 159 160
    }

    return false;
}
161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
NekDouble TetGeom::v_GetLocCoords(const Array<OneD, const NekDouble> &coords,
                                  Array<OneD, NekDouble> &Lcoords)
{
    NekDouble ptdist = 1e6;

    // calculate local coordinates (eta) for coord
    if (GetMetricInfo()->GetGtype() == eRegular)
    {
        // Point inside tetrahedron
        PointGeom r(m_coordim, 0, coords[0], coords[1], coords[2]);

        // Edges
        PointGeom er0, e10, e20, e30;
        er0.Sub(r, *m_verts[0]);
        e10.Sub(*m_verts[1], *m_verts[0]);
        e20.Sub(*m_verts[2], *m_verts[0]);
        e30.Sub(*m_verts[3], *m_verts[0]);

        // Cross products (Normal times area)
        PointGeom cp1020, cp2030, cp3010;
        cp1020.Mult(e10, e20);
        cp2030.Mult(e20, e30);
        cp3010.Mult(e30, e10);

        // Barycentric coordinates (relative volume)
        NekDouble V = e30.dot(cp1020); // Tet Volume={(e30)dot(e10)x(e20)}/6
        NekDouble beta = er0.dot(cp2030) / V; // volume1={(er0)dot(e20)x(e30)}/6
        NekDouble gamma = er0.dot(cp3010) / V; // volume1={(er0)dot(e30)x(e10)}/6
        NekDouble delta = er0.dot(cp1020) / V; // volume1={(er0)dot(e10)x(e20)}/6

        // Make tet bigger
        Lcoords[0] = 2.0 * beta - 1.0;
        Lcoords[1] = 2.0 * gamma - 1.0;
        Lcoords[2] = 2.0 * delta - 1.0;

        // Set ptdist to distance to nearest vertex
        for (int i = 0; i < 4; ++i)
199
        {
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 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
            ptdist = min(ptdist, r.dist(*m_verts[i]));
        }
    }
    else
    {
        v_FillGeom();

        // Determine nearest point of coords  to values in m_xmap
        int npts = m_xmap->GetTotPoints();
        Array<OneD, NekDouble> ptsx(npts), ptsy(npts), ptsz(npts);
        Array<OneD, NekDouble> tmp1(npts), tmp2(npts);

        m_xmap->BwdTrans(m_coeffs[0], ptsx);
        m_xmap->BwdTrans(m_coeffs[1], ptsy);
        m_xmap->BwdTrans(m_coeffs[2], ptsz);

        const Array<OneD, const NekDouble> za = m_xmap->GetPoints(0);
        const Array<OneD, const NekDouble> zb = m_xmap->GetPoints(1);
        const Array<OneD, const NekDouble> zc = m_xmap->GetPoints(2);

        // guess the first local coords based on nearest point
        Vmath::Sadd(npts, -coords[0], ptsx, 1, tmp1, 1);
        Vmath::Vmul(npts, tmp1, 1, tmp1, 1, tmp1, 1);
        Vmath::Sadd(npts, -coords[1], ptsy, 1, tmp2, 1);
        Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);
        Vmath::Sadd(npts, -coords[2], ptsz, 1, tmp2, 1);
        Vmath::Vvtvp(npts, tmp2, 1, tmp2, 1, tmp1, 1, tmp1, 1);

        int min_i = Vmath::Imin(npts, tmp1, 1);

        // distance from coordinate to nearest point for return value.
        ptdist = sqrt(tmp1[min_i]);

        // Get collapsed coordinate
        int qa = za.num_elements(), qb = zb.num_elements();
        Lcoords[2] = zc[min_i / (qa * qb)];
        min_i = min_i % (qa * qb);
        Lcoords[1] = zb[min_i / qa];
        Lcoords[0] = za[min_i % qa];

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

        // Perform newton iteration to find local coordinates
        NekDouble resid = 0.0;
        NewtonIterationForLocCoord(coords, ptsx, ptsy, ptsz, Lcoords, resid);
    }
    return ptdist;
}

int TetGeom::v_GetDir(const int faceidx, const int facedir) const
{
    if (faceidx == 0)
    {
        return facedir;
    }
    else if (faceidx == 1)
    {
        return 2 * facedir;
    }
    else
    {
        return 1 + facedir;
    }
}
266

267 268 269 270
int TetGeom::v_GetVertexEdgeMap(const int i, const int j) const
{
    return VertexEdgeConnectivity[i][j];
}
271

272 273 274 275
int TetGeom::v_GetVertexFaceMap(const int i, const int j) const
{
    return VertexFaceConnectivity[i][j];
}
276

277 278 279 280
int TetGeom::v_GetEdgeFaceMap(const int i, const int j) const
{
    return EdgeFaceConnectivity[i][j];
}
281

282 283
void TetGeom::SetUpLocalEdges()
{
284

285 286 287
    // find edge 0
    int i, j;
    unsigned int check;
288

289
    SegGeomSharedPtr edge;
290

291
    // First set up the 3 bottom edges
292

293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
    if (m_faces[0]->GetEid(0) != m_faces[1]->GetEid(0))
    {
        std::ostringstream errstrm;
        errstrm << "Local edge 0 (eid=" << m_faces[0]->GetEid(0);
        errstrm << ") on face " << m_faces[0]->GetGlobalID();
        errstrm << " must be the same as local edge 0 (eid="
                << m_faces[1]->GetEid(0);
        errstrm << ") on face " << m_faces[1]->GetGlobalID();
        ASSERTL0(false, errstrm.str());
    }

    int faceConnected;
    for (faceConnected = 1; faceConnected < 4; faceConnected++)
    {
        check = 0;
        for (i = 0; i < 3; i++)
        {
            if ((m_faces[0])->GetEid(i) == (m_faces[faceConnected])->GetEid(0))
311
            {
312 313 314 315
                edge = dynamic_pointer_cast<SegGeom>(
                    (m_faces[0])->GetEdge(i));
                m_edges.push_back(edge);
                check++;
316
            }
317
        }
318

319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336
        if (check < 1)
        {
            std::ostringstream errstrm;
            errstrm << "Face 0 does not share an edge with first edge of "
                       "adjacent face. Faces ";
            errstrm << (m_faces[0])->GetGlobalID() << ", "
                    << (m_faces[faceConnected])->GetGlobalID();
            ASSERTL0(false, errstrm.str());
        }
        else if (check > 1)
        {
            std::ostringstream errstrm;
            errstrm << "Connected faces share more than one edge. Faces ";
            errstrm << (m_faces[0])->GetGlobalID() << ", "
                    << (m_faces[faceConnected])->GetGlobalID();
            ASSERTL0(false, errstrm.str());
        }
    }
337

338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 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
    // Then, set up the 3 vertical edges
    check = 0;
    for (i = 0; i < 3; i++) // Set up the vertical edge :face(1) and face(3)
    {
        for (j = 0; j < 3; j++)
        {
            if ((m_faces[1])->GetEid(i) == (m_faces[3])->GetEid(j))
            {
                edge = 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])->GetGlobalID() << ", "
                << (m_faces[3])->GetGlobalID();
        ASSERTL0(false, errstrm.str());
    }
    else if (check > 1)
    {
        std::ostringstream errstrm;
        errstrm << "Connected faces share more than one edge. Faces ";
        errstrm << (m_faces[1])->GetGlobalID() << ", "
                << (m_faces[3])->GetGlobalID();
        ASSERTL0(false, errstrm.str());
    }
    // Set up vertical edges: face(1) through face(3)
    for (faceConnected = 1; faceConnected < 3; faceConnected++)
    {
        check = 0;
        for (i = 0; i < 3; i++)
        {
            for (j = 0; j < 3; j++)
            {
                if ((m_faces[faceConnected])->GetEid(i) ==
                    (m_faces[faceConnected + 1])->GetEid(j))
379
                {
380 381 382 383
                    edge = dynamic_pointer_cast<SegGeom>(
                        (m_faces[faceConnected])->GetEdge(i));
                    m_edges.push_back(edge);
                    check++;
384 385
                }
            }
386
        }
387

388
        if (check < 1)
389
        {
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
            std::ostringstream errstrm;
            errstrm << "Connected faces do not share an edge. Faces ";
            errstrm << (m_faces[faceConnected])->GetGlobalID() << ", "
                    << (m_faces[faceConnected + 1])->GetGlobalID();
            ASSERTL0(false, errstrm.str());
        }
        else if (check > 1)
        {
            std::ostringstream errstrm;
            errstrm << "Connected faces share more than one edge. Faces ";
            errstrm << (m_faces[faceConnected])->GetGlobalID() << ", "
                    << (m_faces[faceConnected + 1])->GetGlobalID();
            ASSERTL0(false, errstrm.str());
        }
    }
};
406

407 408
void TetGeom::SetUpLocalVertices()
{
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 438 439 440 441 442 443 444 445 446 447 448 449 450 451
    // Set up the first 2 vertices (i.e. vertex 0,1)
    if ((m_edges[0]->GetVid(0) == m_edges[1]->GetVid(0)) ||
        (m_edges[0]->GetVid(0) == m_edges[1]->GetVid(1)))
    {
        m_verts.push_back(m_edges[0]->GetVertex(1));
        m_verts.push_back(m_edges[0]->GetVertex(0));
    }
    else if ((m_edges[0]->GetVid(1) == m_edges[1]->GetVid(0)) ||
             (m_edges[0]->GetVid(1) == m_edges[1]->GetVid(1)))
    {
        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]->GetGlobalID() << ", "
                << m_edges[1]->GetGlobalID();
        ASSERTL0(false, errstrm.str());
    }

    // set up the other bottom vertices (i.e. vertex 2)
    for (int i = 1; i < 2; i++)
    {
        if (m_edges[i]->GetVid(0) == m_verts[i]->GetGlobalID())
        {
            m_verts.push_back(m_edges[i]->GetVertex(1));
        }
        else if (m_edges[i]->GetVid(1) == m_verts[i]->GetGlobalID())
        {
            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]->GetGlobalID() << ", "
                    << m_edges[i - 1]->GetGlobalID();
            ASSERTL0(false, errstrm.str());
        }
    }
452

453 454 455 456 457 458 459 460 461
    // set up top vertex
    if (m_edges[3]->GetVid(0) == m_verts[0]->GetGlobalID())
    {
        m_verts.push_back(m_edges[3]->GetVertex(1));
    }
    else
    {
        m_verts.push_back(m_edges[3]->GetVertex(0));
    }
462

463 464 465 466 467 468 469 470 471 472
    // Check the other edges match up.
    int check = 0;
    for (int i = 4; i < 6; ++i)
    {
        if ((m_edges[i]->GetVid(0) == m_verts[i - 3]->GetGlobalID() &&
             m_edges[i]->GetVid(1) == m_verts[3]->GetGlobalID()) ||
            (m_edges[i]->GetVid(1) == m_verts[i - 3]->GetGlobalID() &&
             m_edges[i]->GetVid(0) == m_verts[3]->GetGlobalID()))
        {
            check++;
473
        }
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492
    }
    if (check != 2)
    {
        std::ostringstream errstrm;
        errstrm << "Connected edges do not share a vertex. Edges ";
        errstrm << m_edges[3]->GetGlobalID() << ", "
                << m_edges[2]->GetGlobalID();
        ASSERTL0(false, errstrm.str());
    }
};

void TetGeom::SetUpEdgeOrientation()
{

    // 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
    const unsigned int edgeVerts[kNedges][2] = {
        {0, 1}, {1, 2}, {0, 2}, {0, 3}, {1, 3}, {2, 3}};
493

494 495 496 497
    int i;
    for (i = 0; i < kNedges; i++)
    {
        if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][0]]->GetGlobalID())
498
        {
499
            m_eorient[i] = StdRegions::eForwards;
500
        }
501
        else if (m_edges[i]->GetVid(0) == m_verts[edgeVerts[i][1]]->GetGlobalID())
502
        {
503
            m_eorient[i] = StdRegions::eBackwards;
504
        }
505
        else
506
        {
507
            ASSERTL0(false, "Could not find matching vertex for the edge");
508
        }
509 510 511 512 513
    }
};

void TetGeom::SetUpFaceOrientation()
{
514

515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 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 576 577 578 579 580 581 582
    int f, i;

    // 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)
    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);

    // 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;

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

    // 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][TriGeom::kNverts] = {
        {0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};

    NekDouble dotproduct1 = 0.0;
    NekDouble dotproduct2 = 0.0;

    unsigned int orientation;

    // Loop over all the faces to set up the orientation
    for (f = 0; f < kNqfaces + kNtfaces; f++)
    {
        // initialisation
        elementAaxis_length = 0.0;
        elementBaxis_length = 0.0;
        faceAaxis_length = 0.0;
        faceBaxis_length = 0.0;

        dotproduct1 = 0.0;
        dotproduct2 = 0.0;

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

        // 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...

        // Compute the length of edges on a base-face
        if (baseVertex == m_verts[faceVerts[f][0]]->GetGlobalID())
583
        {
584
            for (i = 0; i < m_coordim; i++)
585
            {
586 587 588 589
                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][0]])[i];
590
            }
591 592 593 594
        }
        else if (baseVertex == m_verts[faceVerts[f][1]]->GetGlobalID())
        {
            for (i = 0; i < m_coordim; i++)
595
            {
596 597 598 599
                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];
600
            }
601 602 603 604
        }
        else if (baseVertex == m_verts[faceVerts[f][2]]->GetGlobalID())
        {
            for (i = 0; i < m_coordim; i++)
605
            {
606 607 608 609
                elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
                                  (*m_verts[faceVerts[f][2]])[i];
                elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
                                  (*m_verts[faceVerts[f][0]])[i];
610 611
            }
        }
612 613 614 615
        else
        {
            ASSERTL0(false, "Could not find matching vertex for the face");
        }
616

617 618 619
        // Now, construct the edge-vectors of the local coordinates of
        // the Geometry2D-representation of the face
        for (i = 0; i < m_coordim; i++)
620
        {
621 622 623 624 625 626 627 628 629 630
            faceAaxis[i] =
                (*m_faces[f]->GetVertex(1))[i] - (*m_faces[f]->GetVertex(0))[i];
            faceBaxis[i] =
                (*m_faces[f]->GetVertex(2))[i] - (*m_faces[f]->GetVertex(0))[i];

            elementAaxis_length += pow(elementAaxis[i], 2);
            elementBaxis_length += pow(elementBaxis[i], 2);
            faceAaxis_length += pow(faceAaxis[i], 2);
            faceBaxis_length += pow(faceBaxis[i], 2);
        }
631

632 633 634 635
        elementAaxis_length = sqrt(elementAaxis_length);
        elementBaxis_length = sqrt(elementBaxis_length);
        faceAaxis_length = sqrt(faceAaxis_length);
        faceBaxis_length = sqrt(faceBaxis_length);
636

637 638 639 640 641 642
        // 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];
        }
643

644 645
        NekDouble norm = fabs(dotproduct1) / elementAaxis_length / faceAaxis_length;
        orientation = 0;
646

647 648 649 650 651 652 653 654
        // if the innerproduct is equal to the (absolute value of the ) products
        // of the lengths of both vectors, then, the coordinate systems will NOT
        // be transposed
        if (fabs(norm - 1.0) < NekConstants::kNekZeroTol)
        {
            // if the inner product is negative, both A-axis point
            // in reverse direction
            if (dotproduct1 < 0.0)
655
            {
656
                orientation += 2;
657 658
            }

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

665
            norm = fabs(dotproduct2) / elementBaxis_length / faceBaxis_length;
666

667 668 669 670 671 672 673
            // check that both these axis are indeed parallel
            ASSERTL1(fabs(norm - 1.0) < NekConstants::kNekZeroTol,
                     "These vectors should be parallel");

            // if the inner product is negative, both B-axis point
            // in reverse direction
            if (dotproduct2 < 0.0)
674
            {
675
                orientation++;
676
            }
677 678 679 680 681 682 683 684 685 686
        }
        // 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++)
687
            {
688
                dotproduct1 += elementAaxis[i] * faceBaxis[i];
689
            }
690

691
            norm = fabs(dotproduct1) / elementAaxis_length / faceBaxis_length;
692

693 694 695
            // check that both these axis are indeed parallel
            ASSERTL1(fabs(norm - 1.0) < NekConstants::kNekZeroTol,
                     "These vectors should be parallel");
696

697 698 699
            // if the result is negative, both axis point in reverse
            // directions
            if (dotproduct1 < 0.0)
700
            {
701
                orientation += 2;
702 703
            }

704 705 706
            // Do the same for the other two corresponding axis
            dotproduct2 = 0.0;
            for (i = 0; i < m_coordim; i++)
707
            {
708
                dotproduct2 += elementBaxis[i] * faceAaxis[i];
709 710
            }

711
            norm = fabs(dotproduct2) / elementBaxis_length / faceAaxis_length;
712

713 714 715
            // check that both these axis are indeed parallel
            ASSERTL1(fabs(norm - 1.0) < NekConstants::kNekZeroTol,
                     "These vectors should be parallel");
716

717
            if (dotproduct2 < 0.0)
718
            {
719 720 721
                orientation++;
            }
        }
722

723
        orientation = orientation + 5;
724

725 726 727 728
        // Fill the m_forient array
        m_forient[f] = (StdRegions::Orientation)orientation;
    }
}
729

730 731 732
void TetGeom::v_Reset(CurveMap &curvedEdges, CurveMap &curvedFaces)
{
    Geometry::v_Reset(curvedEdges, curvedFaces);
733

734 735 736 737 738
    for (int i = 0; i < 4; ++i)
    {
        m_faces[i]->Reset(curvedEdges, curvedFaces);
    }
}
Michael Turner's avatar
Michael Turner committed
739

740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
void TetGeom::v_Setup()
{
    if(!m_setupState)
    {
        for (int i = 0; i < 4; ++i)
        {
            m_faces[i]->Setup();
        }
        SetUpXmap();
        SetUpCoeffs(m_xmap->GetNcoeffs());
        m_setupState = true;
    }
}

/**
 * Generate the geometry factors for this element.
 */
void TetGeom::v_GenGeomFactors()
{
    if(!m_setupState)
    {
        TetGeom::v_Setup();
    }
763

764 765 766
    if (m_geomFactorsState != ePtsFilled)
    {
        GeomType Gtype = eRegular;
767

768
        v_FillGeom();
769

770 771
        // check to see if expansions are linear
        for (int i = 0; i < m_coordim; ++i)
772
        {
773 774 775
            if (m_xmap->GetBasisNumModes(0) != 2 ||
                m_xmap->GetBasisNumModes(1) != 2 ||
                m_xmap->GetBasisNumModes(2) != 2)
776
            {
777
                Gtype = eDeformed;
778 779
            }
        }
Dave Moxey's avatar
Dave Moxey committed
780

781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830
        m_geomFactors = MemoryManager<GeomFactors>::AllocateSharedPtr(
            Gtype, m_coordim, m_xmap, m_coeffs);
        m_geomFactorsState = ePtsFilled;
    }
}

/**
 * @brief Set up the #m_xmap object by determining the order of each
 * direction from derived faces.
 */
void TetGeom::SetUpXmap()
{
    vector<int> tmp;
    tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(0));
    int order0 = *max_element(tmp.begin(), tmp.end());

    tmp.clear();
    tmp.push_back(order0);
    tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(1));
    tmp.push_back(m_faces[0]->GetXmap()->GetEdgeNcoeffs(2));
    int order1 = *max_element(tmp.begin(), tmp.end());

    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));
    int order2 = *max_element(tmp.begin(), tmp.end());

    const LibUtilities::BasisKey A(
        LibUtilities::eModified_A,
        order0,
        LibUtilities::PointsKey(order0+1, LibUtilities::eGaussLobattoLegendre));
    const LibUtilities::BasisKey B(
        LibUtilities::eModified_B,
        order1,
        LibUtilities::PointsKey(order1,
                                LibUtilities::eGaussRadauMAlpha1Beta0));
    const LibUtilities::BasisKey C(
        LibUtilities::eModified_C,
        order2,
        LibUtilities::PointsKey(order2,
                                LibUtilities::eGaussRadauMAlpha2Beta0));

    m_xmap = MemoryManager<StdRegions::StdTetExp>::AllocateSharedPtr(A, B, C);
}

}
}