Nektar++
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
Nektar::SpatialDomains::LinearMeshGraph Class Reference

#include <LinearMeshGraph.hpp>

Public Member Functions

 LinearMeshGraph (SpatialDomains::MeshGraphSharedPtr graph)
 
SpatialDomains::MeshGraphSharedPtr GetLinearGraph ()
 
SpatialDomains::MeshGraphSharedPtr CreateLinearGraph (int nsplit, std::map< int, std::pair< int, std::vector< int > > > &LinCoeffMap, bool UseGLL, bool useSimplex)
 

Private Member Functions

void LinMeshSetUp1DGeom (int nsplit, int voffset, bool UseGLL)
 
void AddSplitEdge (int nsplit, int vid0, int vid1, int maxvertid, int edgeid, Array< OneD, Array< OneD, NekDouble > > &Coords)
 
void LinMeshSetUp2DGeom (int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, int voffset, bool UseGLL, bool useSimplex)
 
void AddSplitTri (int nsplit, Array< TwoD, int > &vids, Array< TwoD, int > &eids_h, Array< TwoD, int > &eids_v, Array< TwoD, int > &eids_d, int maxedgeid, int edgeoffset, int triid, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap)
 
void LinMeshSetUpTetGeom (int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, bool UseGLL)
 
void LinMeshSetUpPrismGeom (int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, bool UseGLL)
 
void LinMeshSetUpCompositesDomain (int nsplit, std::map< int, std::pair< int, std::vector< int > > > &LinCoeffMap, std::map< int, std::map< int, int > > &CoeffMap2D, std::map< int, std::map< int, int > > &CoeffMap3D, bool useSimplex)
 
TiXmlElement * SetupLinearExpansionType (void)
 
int TriShuffleEdges (std::array< SegGeom *, 3 > &edges, std::array< SegGeom *, 3 > *edgesSort=nullptr)
 
void TetShuffleFaces (std::array< TriGeom *, 4 > &faces)
 
void PrismShuffleFaces (std::array< Geometry2D *, 5 > &faces)
 

Private Attributes

const unsigned int TetFaceOrder [4][4][2]
 
SpatialDomains::MeshGraphSharedPtr m_graph
 
SpatialDomains::MeshGraphSharedPtr m_linMesh
 
LibUtilities::SessionReaderSharedPtr m_session
 
int m_meshDimension
 
int m_spaceDimension
 

Detailed Description

Definition at line 93 of file LinearMeshGraph.hpp.

Constructor & Destructor Documentation

◆ LinearMeshGraph()

Nektar::SpatialDomains::LinearMeshGraph::LinearMeshGraph ( SpatialDomains::MeshGraphSharedPtr  graph)
inline

Definition at line 96 of file LinearMeshGraph.hpp.

97 : m_graph(graph), m_session(m_graph->GetSession()),
98 m_meshDimension(m_graph->GetMeshDimension()),
99 m_spaceDimension(m_graph->GetSpaceDimension())
100 {
101 }
Nektar::SpatialDomains::LinearMeshGraph::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Definition LinearMeshGraph.hpp:336
Nektar::SpatialDomains::LinearMeshGraph::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition LinearMeshGraph.hpp:338

Member Function Documentation

◆ AddSplitEdge()

void Nektar::SpatialDomains::LinearMeshGraph::AddSplitEdge ( int  nsplit,
int  vid0,
int  vid1,
int  maxvertid,
int  edgeid,
Array< OneD, Array< OneD, NekDouble > > &  Coords 
)
private

Definition at line 515 of file LinearMeshGraph.hpp.

518{
519 for (int i = 1; i < nsplit; ++i)
520 {
521 int vid = maxvertid + edgeid * (nsplit - 1) + i - 1;
522 PointGeomUniquePtr vert = ObjPoolManager<PointGeom>::AllocateUniquePtr(
523 m_spaceDimension, vid, Coords[0][i], Coords[1][i], Coords[2][i]);
524 m_linMesh->AddGeom(vid, std::move(vert));
525 }
526
527 // add SegGeoms
528 for (int i = 0; i < nsplit; ++i)
529 {
530 int edgid = edgeid * nsplit + i;
531 std::array<PointGeom *, 2> vert;
532
533 if (i == 0)
534 {
535 vert[0] = m_linMesh->GetPointGeom(vid0);
536
537 vert[1] = (nsplit == 1)
538 ? m_linMesh->GetPointGeom(vid1)
539 : m_linMesh->GetPointGeom(maxvertid +
540 edgeid * (nsplit - 1) + i);
541 }
542 else if (i == nsplit - 1)
543 {
544 vert[0] = m_linMesh->GetPointGeom(maxvertid +
545 edgeid * (nsplit - 1) + i - 1);
546 vert[1] = m_linMesh->GetPointGeom(vid1);
547 }
548 else
549 {
550 vert[0] = m_linMesh->GetPointGeom(maxvertid +
551 edgeid * (nsplit - 1) + i - 1);
552 vert[1] =
553 m_linMesh->GetPointGeom(maxvertid + edgeid * (nsplit - 1) + i);
554 }
555
556 SegGeomUniquePtr edge = ObjPoolManager<SegGeom>::AllocateUniquePtr(
557 edgid, m_spaceDimension, vert);
558
559 m_linMesh->AddGeom(edgid, std::move(edge));
560 }
561}
Nektar::ObjPoolManager::AllocateUniquePtr
static std::unique_ptr< DataType, UniquePtrDeleter > AllocateUniquePtr(const Args &...args)
Definition ObjectPool.hpp:86
Nektar::SpatialDomains::LinearMeshGraph::m_linMesh
SpatialDomains::MeshGraphSharedPtr m_linMesh
Definition LinearMeshGraph.hpp:337
Nektar::SpatialDomains::SegGeomUniquePtr
unique_ptr_objpool< SegGeom > SegGeomUniquePtr
Definition MeshGraph.h:102
Nektar::SpatialDomains::PointGeomUniquePtr
unique_ptr_objpool< PointGeom > PointGeomUniquePtr
Definition MeshGraph.h:99

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_linMesh, and m_spaceDimension.

Referenced by LinMeshSetUp1DGeom().

◆ AddSplitTri()

void Nektar::SpatialDomains::LinearMeshGraph::AddSplitTri ( int  nsplit,
Array< TwoD, int > &  vids,
Array< TwoD, int > &  eids_h,
Array< TwoD, int > &  eids_v,
Array< TwoD, int > &  eids_d,
int  maxedgeid,
int  edgeoffset,
int  triid,
std::map< int, int > &  FceEdgOffset,
std::map< int, std::map< int, int > > &  CoeffMap 
)
private

Definition at line 952 of file LinearMeshGraph.hpp.

960{
961 // Add in horizontal edges
962 int cnt = 0;
963 for (int j = 1; j < nsplit; ++j)
964 {
965 for (int i = 0; i < nsplit - j; ++i)
966 {
967 int edgid = maxedgeid + triid * edgeoffset + cnt++;
968 eids_h[j][i] = edgid;
969
970 std::array<PointGeom *, 2> vert;
971 vert[0] = m_linMesh->GetPointGeom(vids[j][i]);
972 vert[1] = m_linMesh->GetPointGeom(vids[j][i + 1]);
973
974 SegGeomUniquePtr edge = ObjPoolManager<SegGeom>::AllocateUniquePtr(
975 edgid, m_spaceDimension, vert);
976
977 m_linMesh->AddGeom(edgid, std::move(edge));
978 }
979 }
980
981 // Add in vertical and diagonal edges
982 for (int j = 0; j < nsplit; ++j)
983 {
984 for (int i = 1; i < nsplit - j; ++i)
985 {
986 int edgid = maxedgeid + triid * edgeoffset + cnt++;
987
988 eids_v[j][i] = edgid;
989
990 std::array<PointGeom *, 2> vert;
991 vert[0] = m_linMesh->GetPointGeom(vids[j][i]);
992 vert[1] = m_linMesh->GetPointGeom(vids[j + 1][i]);
993 SegGeomUniquePtr edge = ObjPoolManager<SegGeom>::AllocateUniquePtr(
994 edgid, m_spaceDimension, vert);
995
996 m_linMesh->AddGeom(edgid, std::move(edge));
997
998 edgid = maxedgeid + triid * edgeoffset + cnt++;
999
1000 eids_d[j][i] = edgid;
1001
1002 vert[0] = m_linMesh->GetPointGeom(vids[j][i]);
1003 vert[1] = m_linMesh->GetPointGeom(vids[j + 1][i - 1]);
1004 SegGeomUniquePtr edge_d =
1005 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1006 edgid, m_spaceDimension, vert);
1007
1008 m_linMesh->AddGeom(edgid, std::move(edge_d));
1009 }
1010 }
1011
1012 // add in tri elements
1013 cnt = 0;
1014 for (int j = 0; j < nsplit; ++j)
1015 {
1016 // lower triangular elements.
1017 for (int i = 0; i < nsplit - j; ++i)
1018 {
1019 int newtriid = 2 * triid * nsplit * nsplit + cnt++;
1020
1021 std::array<SegGeom *, 3> edges;
1022 edges[0] = m_linMesh->GetSegGeom(eids_h[j][i]);
1023 edges[1] = m_linMesh->GetSegGeom(eids_d[j][i + 1]);
1024 edges[2] = m_linMesh->GetSegGeom(eids_v[j][i]);
1025 // shuffle so lowest vid is between edges 1,2
1026 FceEdgOffset[newtriid] = TriShuffleEdges(edges);
1027
1028 TriGeomUniquePtr tri =
1029 ObjPoolManager<SpatialDomains::TriGeom>::AllocateUniquePtr(
1030 newtriid, edges);
1031
1032 m_linMesh->AddGeom(newtriid, std::move(tri));
1033 }
1034
1035 // upper triangular elements.
1036 for (int i = 1; i < nsplit - j; ++i)
1037 {
1038 int newtriid = 2 * triid * nsplit * nsplit + cnt++;
1039
1040 std::array<SegGeom *, 3> edges;
1041 edges[0] = m_linMesh->GetSegGeom(eids_d[j][i]);
1042 edges[1] = m_linMesh->GetSegGeom(eids_v[j][i]);
1043 edges[2] = m_linMesh->GetSegGeom(eids_h[j + 1][i - 1]);
1044 // shuffle so lowest vid is on edges 1,2
1045 FceEdgOffset[newtriid] = TriShuffleEdges(edges);
1046 TriGeomUniquePtr tri =
1047 ObjPoolManager<SpatialDomains::TriGeom>::AllocateUniquePtr(
1048 newtriid, edges);
1049 m_linMesh->AddGeom(newtriid, std::move(tri));
1050 }
1051 }
1052
1053 // collect a mapping of Vid to point vertex location
1054 cnt = 0;
1055 for (int j = 0; j < nsplit + 1; ++j)
1056 {
1057 for (int i = 0; i < nsplit + 1 - j; ++i)
1058 {
1059 CoeffMap[triid][vids[j][i]] = cnt++;
1060 }
1061 }
1062}
Nektar::SpatialDomains::LinearMeshGraph::TriShuffleEdges
int TriShuffleEdges(std::array< SegGeom *, 3 > &edges, std::array< SegGeom *, 3 > *edgesSort=nullptr)
Definition LinearMeshGraph.hpp:199
Nektar::SpatialDomains::TriGeomUniquePtr
unique_ptr_objpool< TriGeom > TriGeomUniquePtr
Definition MeshGraph.h:103

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_linMesh, m_spaceDimension, and TriShuffleEdges().

Referenced by LinMeshSetUp2DGeom().

◆ CreateLinearGraph()

SpatialDomains::MeshGraphSharedPtr Nektar::SpatialDomains::LinearMeshGraph::CreateLinearGraph ( int  nsplit,
std::map< int, std::pair< int, std::vector< int > > > &  LinCoeffMap,
bool  UseGLL,
bool  useSimplex 
)

Definition at line 343 of file LinearMeshGraph.hpp.

346{
347 //--------------------------------------------------------
348 // initial setup
349 //--------------------------------------------------------
350 LibUtilities::CommSharedPtr comm = m_session->GetComm();
351 ASSERTL0(comm.get(), "Communication not initialised.");
352
353 // Populate SessionReader. This should be done only on the root process
354 // so that we can partition appropriately without all processes having
355 // to read in the input file.
356 const bool isRoot = comm->TreatAsRankZero();
357 std::string geomType;
358
359 if (isRoot)
360 {
361 // Get geometry type, i.e. XML (compressed/uncompressed) or HDF5.
362 geomType = m_session->GetGeometryType();
363
364 // Convert to a vector of chars so that we can broadcast.
365 std::vector<char> v(geomType.c_str(),
366 geomType.c_str() + geomType.length());
367
368 size_t length = v.size();
369 comm->Bcast(length, 0);
370 comm->Bcast(v, 0);
371 }
372 else
373 {
374 size_t length;
375 comm->Bcast(length, 0);
376
377 std::vector<char> v(length);
378 comm->Bcast(v, 0);
379
380 geomType = std::string(v.begin(), v.end());
381 }
382
383 m_linMesh = MemoryManager<SpatialDomains::MeshGraph>::AllocateSharedPtr();
384
385 // copy starting info
386 m_linMesh->SetPartitionNumber(m_graph->GetPartitionNumber());
387 m_linMesh->SetMeshPartitioned(m_graph->GetMeshPartitioned());
388 m_linMesh->SetMeshDimension(m_meshDimension);
389 m_linMesh->SetSpaceDimension(m_spaceDimension);
390 m_linMesh->SetSession(m_graph->GetSession());
391
392 //--------------------------------------------------------
393 // Need to offset vertices so that new face vertids are lower than
394 // original vertices, edges so that new orientation are aligned
395 // (particularly for prisms).
396 //--------------------------------------------------------
397 int voffset = 0;
398 auto &prismGeoms = m_graph->GetGeomMap<PrismGeom>();
399 if (m_meshDimension == 3 && prismGeoms.size() > 0)
400 {
401 // check for alignment of triangular faces or edge id of bottom edges
402 for (auto [id, prism] : prismGeoms)
403 {
404 if (prism->GetEorient(1) != prism->GetEorient(3))
405 {
406 NEKERROR(ErrorUtil::ewarning,
407 "triangular faces are not aligned ");
408 }
409 }
410
411 ASSERTL0(useSimplex == false, "Prisms are not currently "
412 "setup for simplex LOR scheme");
413 }
414
415 //--------------------------------------------------------
416 // Set up 1D elements
417 //--------------------------------------------------------
418 LinMeshSetUp1DGeom(nsplit, voffset, UseGLL);
419
420 //--------------------------------------------------------
421 // Set up 2D elements
422 //--------------------------------------------------------
423 std::map<int, std::map<int, int>> CoeffMap2D;
424 std::map<int, int> FceEdgOffset; // map of face id and offset when shuffled
425
426 LinMeshSetUp2DGeom(nsplit, FceEdgOffset, CoeffMap2D, voffset, UseGLL,
427 useSimplex);
428
429 //--------------------------------------------------------
430 // Set up 3D elements
431 //--------------------------------------------------------
432 std::map<int, std::map<int, int>> CoeffMap3D;
433 if (m_meshDimension == 3)
434 {
435 LinMeshSetUpTetGeom(nsplit, FceEdgOffset, CoeffMap3D, UseGLL);
436 LinMeshSetUpPrismGeom(nsplit, FceEdgOffset, CoeffMap3D, UseGLL);
437 }
438
439 //--------------------------------------------------------
440 // Reset composites for BCs and domain definitions
441 //--------------------------------------------------------
442 LinMeshSetUpCompositesDomain(nsplit, LinCoeffMap, CoeffMap2D, CoeffMap3D,
443 useSimplex);
444
445 return m_linMesh;
446}
ASSERTL0
#define ASSERTL0(condition, msg)
Definition ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition ErrorUtil.hpp:210
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition NekMemoryManager.hpp:167
Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUp1DGeom
void LinMeshSetUp1DGeom(int nsplit, int voffset, bool UseGLL)
Definition LinearMeshGraph.hpp:448
Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUp2DGeom
void LinMeshSetUp2DGeom(int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, int voffset, bool UseGLL, bool useSimplex)
Definition LinearMeshGraph.hpp:563
Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpCompositesDomain
void LinMeshSetUpCompositesDomain(int nsplit, std::map< int, std::pair< int, std::vector< int > > > &LinCoeffMap, std::map< int, std::map< int, int > > &CoeffMap2D, std::map< int, std::map< int, int > > &CoeffMap3D, bool useSimplex)
Definition LinearMeshGraph.hpp:2405
Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpPrismGeom
void LinMeshSetUpPrismGeom(int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, bool UseGLL)
Definition LinearMeshGraph.hpp:1819
Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpTetGeom
void LinMeshSetUpTetGeom(int nsplit, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap, bool UseGLL)
Definition LinearMeshGraph.hpp:1064
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition Comm.h:55

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::ErrorUtil::ewarning, LinMeshSetUp1DGeom(), LinMeshSetUp2DGeom(), LinMeshSetUpCompositesDomain(), LinMeshSetUpPrismGeom(), LinMeshSetUpTetGeom(), m_graph, m_linMesh, m_meshDimension, m_session, m_spaceDimension, and NEKERROR.

Referenced by Nektar::MultiRegions::PreconditionerLOR::v_BuildPreconditioner().

◆ GetLinearGraph()

SpatialDomains::MeshGraphSharedPtr Nektar::SpatialDomains::LinearMeshGraph::GetLinearGraph ( )
inline

Definition at line 103 of file LinearMeshGraph.hpp.

104 {
105 return m_linMesh;
106 }

References m_linMesh.

Referenced by Nektar::MultiRegions::PreconditionerLOR::v_BuildPreconditioner().

◆ LinMeshSetUp1DGeom()

void Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUp1DGeom ( int  nsplit,
int  voffset,
bool  UseGLL 
)
private

set up a deep copy of existing vertices with voffset

Definition at line 448 of file LinearMeshGraph.hpp.

449{
450 LibUtilities::CommSharedPtr comm = m_session->GetComm();
451
452 /** set up a deep copy of existing vertices with voffset **/
453 for (auto [id, vert] : m_graph->GetGeomMap<PointGeom>())
454 {
455 int vid = voffset + vert->GetGlobalID();
456 PointGeomUniquePtr v = ObjPoolManager<PointGeom>::AllocateUniquePtr(
457 m_spaceDimension, vid, (*vert)[0], (*vert)[1], (*vert)[2]);
458
459 m_linMesh->AddGeom(vid, std::move(v));
460 }
461
462 int maxvertid = m_linMesh->GetGeomMap<PointGeom>().rbegin()->first + 1;
463 comm->AllReduce(maxvertid, LibUtilities::ReduceMax);
464
465 Array<OneD, Array<OneD, NekDouble>> Coords(3);
466 Coords[0] = Array<OneD, NekDouble>(nsplit + 1);
467 Coords[1] = Array<OneD, NekDouble>(nsplit + 1);
468 Coords[2] = Array<OneD, NekDouble>(nsplit + 1);
469
470 //--------------------------------------------------------
471 // Split edges of macro mesh into new edges and vertices
472 //--------------------------------------------------------
473
474 for (auto [id, edgegeom] : m_graph->GetGeomMap<SegGeom>())
475 {
476 StdRegions::StdExpansionSharedPtr Xmap = edgegeom->GetXmap();
477
478 Array<OneD, NekDouble> tmp(Xmap->GetTotPoints());
479 if (UseGLL) // get GLL points
480 {
481 Xmap->BwdTrans(edgegeom->GetCoeffs(0), tmp);
482 Xmap->PhysInterpToGLL(tmp, Coords[0], nsplit + 1);
483 Xmap->BwdTrans(edgegeom->GetCoeffs(1), tmp);
484 Xmap->PhysInterpToGLL(tmp, Coords[1], nsplit + 1);
485
486 if (m_meshDimension == 3)
487 {
488 Xmap->BwdTrans(edgegeom->GetCoeffs(2), tmp);
489 Xmap->PhysInterpToGLL(tmp, Coords[2], nsplit + 1);
490 }
491 }
492 else // get equispaced points
493 {
494 Xmap->BwdTrans(edgegeom->GetCoeffs(0), tmp);
495 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Coords[0], nsplit + 1);
496 Xmap->BwdTrans(edgegeom->GetCoeffs(1), tmp);
497 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Coords[1], nsplit + 1);
498
499 if (m_meshDimension == 3)
500 {
501 Xmap->BwdTrans(edgegeom->GetCoeffs(2), tmp);
502 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Coords[2], nsplit + 1);
503 }
504 }
505
506 // replace with method with api EdgeAdd(nsplit,vid0,vid1,
507 // maxvertid, edgeid)
508 int edgeid = edgegeom->GetGlobalID();
509 AddSplitEdge(nsplit, edgegeom->GetVertex(0)->GetGlobalID(),
510 edgegeom->GetVertex(1)->GetGlobalID(), maxvertid, edgeid,
511 Coords);
512 }
513}
Nektar::SpatialDomains::LinearMeshGraph::AddSplitEdge
void AddSplitEdge(int nsplit, int vid0, int vid1, int maxvertid, int edgeid, Array< OneD, Array< OneD, NekDouble > > &Coords)
Definition LinearMeshGraph.hpp:515
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition StdExpansion.h:1817

References AddSplitEdge(), Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), m_graph, m_linMesh, m_meshDimension, m_session, m_spaceDimension, and Nektar::LibUtilities::ReduceMax.

Referenced by CreateLinearGraph().

◆ LinMeshSetUp2DGeom()

void Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUp2DGeom ( int  nsplit,
std::map< int, int > &  FceEdgOffset,
std::map< int, std::map< int, int > > &  CoeffMap,
int  voffset,
bool  UseGLL,
bool  useSimplex 
)
private

Definition at line 563 of file LinearMeshGraph.hpp.

567{
568 LibUtilities::CommSharedPtr comm = m_session->GetComm();
569
570 Array<OneD, NekDouble> Xpts, Ypts, Zpts;
571 Xpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1));
572 Ypts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1));
573 Zpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1));
574
575 // use max vertid
576 int maxvertid =
577 m_graph->GetGeomMap<PointGeom>().rbegin()->first + 1 + voffset;
578 comm->AllReduce(maxvertid, LibUtilities::ReduceMax);
579
580 int maxedgeid = m_linMesh->GetGeomMap<SegGeom>().rbegin()->first + 1;
581 comm->AllReduce(maxedgeid, LibUtilities::ReduceMax);
582
583 //--------------------------------------------------------
584 // Split tri faces of macro mesh into new edges and vertices
585 //--------------------------------------------------------
586 Array<TwoD, int> vids(nsplit + 1, nsplit + 1, 0);
587 // horirzotnal edges
588 Array<TwoD, int> eids_h(nsplit + 1, nsplit, 0);
589 // vertical edges
590 Array<TwoD, int> eids_v(nsplit, nsplit + 1, 0);
591 // diagonal edges
592 Array<TwoD, int> eids_d(nsplit + 1, nsplit + 1, 0);
593
594 // number of interior edges or quad if split into simplices ~ 3
595 // npslit*nsplit so need to use that here too
596 int edgeoffset = nsplit * (3 * nsplit - 2);
597
598 // add vertices, edges for tri geoms
599 for (auto [id, geom] : m_graph->GetGeomMap<TriGeom>())
600 {
601 StdRegions::StdExpansionSharedPtr Xmap = geom->GetXmap();
602
603 Array<OneD, NekDouble> tmp(Xmap->GetTotPoints());
604
605 if (UseGLL) // get GLL points
606 {
607 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
608 Xmap->PhysInterpToGLL(tmp, Xpts, nsplit + 1);
609 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
610 Xmap->PhysInterpToGLL(tmp, Ypts, nsplit + 1);
611
612 if (m_meshDimension == 3)
613 {
614 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
615 Xmap->PhysInterpToGLL(tmp, Zpts, nsplit + 1);
616 }
617 }
618 else // get equispaced points
619 {
620 // get equispaced points
621 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
622 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Xpts, nsplit + 1);
623 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
624 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Ypts, nsplit + 1);
625 if (m_meshDimension == 3)
626 {
627 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
628 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Zpts, nsplit + 1);
629 }
630 }
631
632 int triid = geom->GetGlobalID();
633
634 int edgeid0 = geom->GetEdge(0)->GetGlobalID();
635 int edgeid1 = geom->GetEdge(1)->GetGlobalID();
636 int edgeid2 = geom->GetEdge(2)->GetGlobalID();
637
638 // Determine if edges are cartesian aligned:
639 bool fwd0 = geom->GetEorient(0) == StdRegions::eForwards;
640 bool fwd1 = geom->GetEorient(1) == StdRegions::eForwards;
641 bool fwd2 = geom->GetEorient(2) == StdRegions::eForwards;
642
643 // fill in vertex ids along boundary of element //
644 for (int i = 1; i < nsplit; ++i)
645 {
646 vids[0][i] =
647 fwd0 ? maxvertid + edgeid0 * (nsplit - 1) + i - 1
648 : maxvertid + edgeid0 * (nsplit - 1) + nsplit - i - 1;
649 vids[i][nsplit - i] =
650 fwd1 ? maxvertid + edgeid1 * (nsplit - 1) + i - 1
651 : maxvertid + edgeid1 * (nsplit - 1) + nsplit - i - 1;
652 vids[i][0] =
653 fwd2 ? maxvertid + edgeid2 * (nsplit - 1) + i - 1
654 : maxvertid + edgeid2 * (nsplit - 1) + nsplit - i - 1;
655 }
656
657 // set up corner ids
658 vids[0][0] =
659 fwd0 ? m_graph->GetSegGeom(edgeid0)->GetVertex(0)->GetGlobalID()
660 : m_graph->GetSegGeom(edgeid0)->GetVertex(1)->GetGlobalID();
661 vids[0][nsplit] =
662 fwd0 ? m_graph->GetSegGeom(edgeid0)->GetVertex(1)->GetGlobalID()
663 : m_graph->GetSegGeom(edgeid0)->GetVertex(0)->GetGlobalID();
664 vids[nsplit][0] =
665 fwd2 ? m_graph->GetSegGeom(edgeid2)->GetVertex(1)->GetGlobalID()
666 : m_graph->GetSegGeom(edgeid2)->GetVertex(0)->GetGlobalID();
667
668 // fill in edge ids along boundary of element
669 for (int i = 0; i < nsplit; ++i)
670 {
671 eids_h[0][i] =
672 fwd0 ? edgeid0 * nsplit + i : edgeid0 * nsplit + nsplit - i - 1;
673 eids_v[i][0] =
674 fwd2 ? edgeid2 * nsplit + i : edgeid2 * nsplit + nsplit - i - 1;
675 eids_d[i][nsplit - i] =
676 fwd1 ? edgeid1 * nsplit + i : edgeid1 * nsplit + nsplit - i - 1;
677 }
678
679 // add in vertices on interior of element
680 int cnt = nsplit + 1;
681 int vid;
682 for (int j = 1; j < nsplit; ++j)
683 {
684 for (int i = 1; i < nsplit - j; ++i)
685 {
686 vid = triid * (nsplit - 1) * (nsplit - 1) +
687 (j - 1) * (nsplit - 1) + (i - 1);
688 if (voffset == 0)
689 {
690 vid += maxvertid + maxedgeid * (nsplit - 1) / nsplit;
691 }
692
693 vids[j][i] = vid;
694 PointGeomUniquePtr vert =
695 ObjPoolManager<PointGeom>::AllocateUniquePtr(
696 m_spaceDimension, vid, Xpts[cnt + i], Ypts[cnt + i],
697 Zpts[cnt + i]);
698 m_linMesh->AddGeom(vid, std::move(vert));
699 }
700 cnt += nsplit + 1 - j;
701 }
702
703 AddSplitTri(nsplit, vids, eids_h, eids_v, eids_d, maxedgeid, edgeoffset,
704 triid, FceEdgOffset, CoeffMap);
705 }
706
707 //--------------------------------------------------------
708 // Split quad faces of macro mesh into new edges and vertices
709 //--------------------------------------------------------
710
711 // add vertices, edges for quad geoms
712 for (auto [id, quadgeom] : m_graph->GetGeomMap<QuadGeom>())
713 {
714 StdRegions::StdExpansionSharedPtr Xmap = quadgeom->GetXmap();
715
716 Array<OneD, NekDouble> tmp(Xmap->GetTotPoints());
717
718 if (UseGLL) // get GLL points
719 {
720 Xmap->BwdTrans(quadgeom->GetCoeffs(0), tmp);
721 Xmap->PhysInterpToGLL(tmp, Xpts, nsplit + 1);
722 Xmap->BwdTrans(quadgeom->GetCoeffs(1), tmp);
723 Xmap->PhysInterpToGLL(tmp, Ypts, nsplit + 1);
724
725 if (m_meshDimension == 3)
726 {
727 Xmap->BwdTrans(quadgeom->GetCoeffs(2), tmp);
728 Xmap->PhysInterpToGLL(tmp, Zpts, nsplit + 1);
729 }
730 }
731 else // get equispaced points
732 {
733 // get equispaced points
734 Xmap->BwdTrans(quadgeom->GetCoeffs(0), tmp);
735 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Xpts, nsplit + 1);
736 Xmap->BwdTrans(quadgeom->GetCoeffs(1), tmp);
737 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Ypts, nsplit + 1);
738
739 if (m_meshDimension == 3)
740 {
741 Xmap->BwdTrans(quadgeom->GetCoeffs(2), tmp);
742 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Zpts, nsplit + 1);
743 }
744 }
745
746 int quadid = quadgeom->GetGlobalID();
747
748 int edgeid0 = quadgeom->GetEdge(0)->GetGlobalID();
749 int edgeid1 = quadgeom->GetEdge(1)->GetGlobalID();
750 int edgeid2 = quadgeom->GetEdge(2)->GetGlobalID();
751 int edgeid3 = quadgeom->GetEdge(3)->GetGlobalID();
752
753 // Determine if edges are cartesian aligned:
754 bool fwd0 = quadgeom->GetEorient(0) == StdRegions::eForwards;
755 bool fwd1 = quadgeom->GetEorient(1) == StdRegions::eForwards;
756 bool fwd2 = quadgeom->GetEorient(2) == StdRegions::eForwards;
757 bool fwd3 = quadgeom->GetEorient(3) == StdRegions::eForwards;
758
759 // fill in vertex ids along boundary of element //
760 for (int i = 1; i < nsplit; ++i)
761 {
762 vids[0][i] =
763 fwd0 ? maxvertid + edgeid0 * (nsplit - 1) + i - 1
764 : maxvertid + edgeid0 * (nsplit - 1) + nsplit - i - 1;
765 vids[i][nsplit] =
766 fwd1 ? maxvertid + edgeid1 * (nsplit - 1) + i - 1
767 : maxvertid + edgeid1 * (nsplit - 1) + nsplit - i - 1;
768
769 // Need to check if fwd is anticlockwise or cartesian -
770 // set up for cartesian
771 vids[nsplit][i] =
772 fwd2 ? maxvertid + edgeid2 * (nsplit - 1) + i - 1
773 : maxvertid + edgeid2 * (nsplit - 1) + nsplit - i - 1;
774 vids[i][0] =
775 fwd3 ? maxvertid + edgeid3 * (nsplit - 1) + i - 1
776 : maxvertid + edgeid3 * (nsplit - 1) + nsplit - i - 1;
777 }
778
779 // set up corner ids
780 vids[0][0] =
781 fwd0 ? m_graph->GetSegGeom(edgeid0)->GetVertex(0)->GetGlobalID()
782 : m_graph->GetSegGeom(edgeid0)->GetVertex(1)->GetGlobalID();
783 vids[0][nsplit] =
784 fwd0 ? m_graph->GetSegGeom(edgeid0)->GetVertex(1)->GetGlobalID()
785 : m_graph->GetSegGeom(edgeid0)->GetVertex(0)->GetGlobalID();
786 vids[nsplit][0] =
787 fwd2 ? m_graph->GetSegGeom(edgeid2)->GetVertex(0)->GetGlobalID()
788 : m_graph->GetSegGeom(edgeid2)->GetVertex(1)->GetGlobalID();
789 vids[nsplit][nsplit] =
790 fwd2 ? m_graph->GetSegGeom(edgeid2)->GetVertex(1)->GetGlobalID()
791 : m_graph->GetSegGeom(edgeid2)->GetVertex(0)->GetGlobalID();
792
793 // fill in edge ids along boundary of element
794 for (int i = 0; i < nsplit; ++i)
795 {
796 eids_h[0][i] =
797 fwd0 ? edgeid0 * nsplit + i : edgeid0 * nsplit + nsplit - i - 1;
798 eids_h[nsplit][i] =
799 fwd2 ? edgeid2 * nsplit + i : edgeid2 * nsplit + nsplit - i - 1;
800 eids_v[i][nsplit] =
801 fwd1 ? edgeid1 * nsplit + i : edgeid1 * nsplit + nsplit - i - 1;
802 eids_v[i][0] =
803 fwd3 ? edgeid3 * nsplit + i : edgeid3 * nsplit + nsplit - i - 1;
804 }
805
806 // add in vertices on interior of element
807 int vid;
808 for (int j = 1; j < nsplit; ++j)
809 {
810 for (int i = 1; i < nsplit; ++i)
811 {
812 vid = quadid * (nsplit - 1) * (nsplit - 1) +
813 (j - 1) * (nsplit - 1) + (i - 1);
814
815 if (voffset == 0)
816 {
817 vid += maxvertid + maxedgeid * (nsplit - 1) / nsplit;
818 }
819
820 vids[j][i] = vid;
821 PointGeomUniquePtr vert =
822 ObjPoolManager<PointGeom>::AllocateUniquePtr(
823 m_spaceDimension, vid, Xpts[j * (nsplit + 1) + i],
824 Ypts[j * (nsplit + 1) + i], Zpts[j * (nsplit + 1) + i]);
825 m_linMesh->AddGeom(vid, std::move(vert));
826 }
827 }
828
829 // Add in horizontal edges
830 // KK - we might not need to create explicitly the edges now
831 int cnt = 0;
832 for (int j = 1; j < nsplit; ++j)
833 {
834 for (int i = 0; i < nsplit; ++i)
835 {
836 int edgid = maxedgeid + quadid * edgeoffset + cnt++;
837
838 eids_h[j][i] = edgid;
839
840 std::array<PointGeom *, 2> vert;
841 vert[0] = m_linMesh->GetPointGeom(vids[j][i]);
842 vert[1] = m_linMesh->GetPointGeom(vids[j][i + 1]);
843
844 SegGeomUniquePtr edge =
845 ObjPoolManager<SegGeom>::AllocateUniquePtr(
846 edgid, m_spaceDimension, vert);
847
848 m_linMesh->AddGeom(edgid, std::move(edge));
849 }
850 }
851
852 // Add in vertical edges
853 // KK - we might not need to create explicitly the edges now
854 for (int j = 0; j < nsplit; ++j)
855 {
856 for (int i = 1; i < nsplit; ++i)
857 {
858 int edgid = maxedgeid + quadid * edgeoffset + cnt++;
859
860 eids_v[j][i] = edgid;
861
862 std::array<PointGeom *, 2> vert;
863 vert[0] = m_linMesh->GetPointGeom(vids[j][i]);
864 vert[1] = m_linMesh->GetPointGeom(vids[j + 1][i]);
865 SegGeomUniquePtr edge =
866 ObjPoolManager<SegGeom>::AllocateUniquePtr(
867 edgid, m_spaceDimension, vert);
868
869 m_linMesh->AddGeom(edgid, std::move(edge));
870 }
871 }
872
873 if (useSimplex)
874 {
875 // Add in vertical and diagonal edges
876 for (int j = 0; j < nsplit; ++j)
877 {
878 for (int i = 0; i < nsplit; ++i)
879 {
880 int edgid = maxedgeid + quadid * edgeoffset + cnt++;
881
882 eids_d[j][i] = edgid;
883
884 std::array<PointGeom *, 2> vert;
885 vert[0] = m_linMesh->GetPointGeom(vids[j][i + 1]);
886 vert[1] = m_linMesh->GetPointGeom(vids[j + 1][i]);
887 SegGeomUniquePtr edge_d =
888 ObjPoolManager<SegGeom>::AllocateUniquePtr(
889 edgid, m_spaceDimension, vert);
890 m_linMesh->AddGeom(edgid, std::move(edge_d));
891 }
892 }
893
894 // add in quad elements
895 for (int j = 0; j < nsplit; ++j)
896 {
897 for (int i = 0; i < nsplit; ++i)
898 {
899 int newtriid =
900 2 * (quadid * nsplit * nsplit + j * nsplit + i);
901
902 std::array<SegGeom *, 3> edges;
903 edges[0] = m_linMesh->GetSegGeom(eids_h[j][i]);
904 edges[1] = m_linMesh->GetSegGeom(eids_d[j][i]);
905 edges[2] = m_linMesh->GetSegGeom(eids_v[j][i]);
906 TriGeomUniquePtr tri = ObjPoolManager<
907 SpatialDomains::TriGeom>::AllocateUniquePtr(newtriid,
908 edges);
909 m_linMesh->AddGeom(newtriid, std::move(tri));
910
911 newtriid =
912 2 * (quadid * nsplit * nsplit + j * nsplit + i) + 1;
913
914 edges[0] = m_linMesh->GetSegGeom(eids_h[j + 1][i]);
915 edges[1] = m_linMesh->GetSegGeom(eids_d[j][i]);
916 edges[2] = m_linMesh->GetSegGeom(eids_v[j][i + 1]);
917
918 TriGeomUniquePtr tri1 = ObjPoolManager<
919 SpatialDomains::TriGeom>::AllocateUniquePtr(newtriid,
920 edges);
921
922 m_linMesh->AddGeom(newtriid, std::move(tri1));
923 }
924 }
925 }
926 else
927 {
928 // add in elements
929 for (int j = 0; j < nsplit; ++j)
930 {
931 for (int i = 0; i < nsplit; ++i)
932 {
933 int newquadid =
934 2 * quadid * nsplit * nsplit + j * nsplit + i;
935
936 std::array<SegGeom *, 4> edges;
937 edges[0] = m_linMesh->GetSegGeom(eids_h[j][i]);
938 edges[2] = m_linMesh->GetSegGeom(eids_h[j + 1][i]);
939 edges[1] = m_linMesh->GetSegGeom(eids_v[j][i + 1]);
940 edges[3] = m_linMesh->GetSegGeom(eids_v[j][i]);
941
942 QuadGeomUniquePtr quad = ObjPoolManager<
943 SpatialDomains::QuadGeom>::AllocateUniquePtr(newquadid,
944 edges);
945 m_linMesh->AddGeom(newquadid, std::move(quad));
946 }
947 }
948 }
949 }
950}
Nektar::SpatialDomains::LinearMeshGraph::AddSplitTri
void AddSplitTri(int nsplit, Array< TwoD, int > &vids, Array< TwoD, int > &eids_h, Array< TwoD, int > &eids_v, Array< TwoD, int > &eids_d, int maxedgeid, int edgeoffset, int triid, std::map< int, int > &FceEdgOffset, std::map< int, std::map< int, int > > &CoeffMap)
Definition LinearMeshGraph.hpp:952
Nektar::SpatialDomains::QuadGeomUniquePtr
unique_ptr_objpool< QuadGeom > QuadGeomUniquePtr
Definition MeshGraph.h:104

References AddSplitTri(), Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), Nektar::StdRegions::eForwards, m_graph, m_linMesh, m_meshDimension, m_session, m_spaceDimension, and Nektar::LibUtilities::ReduceMax.

Referenced by CreateLinearGraph().

◆ LinMeshSetUpCompositesDomain()

void Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpCompositesDomain ( int  nsplit,
std::map< int, std::pair< int, std::vector< int > > > &  LinCoeffMap,
std::map< int, std::map< int, int > > &  CoeffMap2D,
std::map< int, std::map< int, int > > &  CoeffMap3D,
bool  useSimplex 
)
private

Definition at line 2405 of file LinearMeshGraph.hpp.

2409{
2410 auto &linMeshComp = m_linMesh->GetComposites();
2411
2412 for (auto &c : m_graph->GetComposites())
2413 {
2414 CompositeSharedPtr curVector =
2415 MemoryManager<Composite>::AllocateSharedPtr();
2416 linMeshComp[c.first] = curVector;
2417
2418 for (auto &g : c.second->m_geomVec)
2419 {
2420 switch (g->GetShapeType())
2421 {
2422 case LibUtilities::eSegment:
2423 {
2424 int id = g->GetGlobalID();
2425 for (int i = 0; i < nsplit; ++i)
2426 {
2427 int newid = id * nsplit + i;
2428
2429 auto it = m_linMesh->GetGeomMap<SegGeom>().find(newid);
2430 ASSERTL1(it != m_linMesh->GetGeomMap<SegGeom>().end(),
2431 "Failed to find LinMesh edge id");
2432 linMeshComp[c.first]->m_geomVec.push_back(it->second);
2433 }
2434 }
2435 break;
2436 case LibUtilities::eTriangle:
2437 {
2438 int id = g->GetGlobalID();
2439 for (int i = 0; i < nsplit * nsplit; ++i)
2440 {
2441 int newid = 2 * id * nsplit * nsplit + i;
2442
2443 auto it = m_linMesh->GetGeomMap<TriGeom>().find(newid);
2444 ASSERTL1(it != m_linMesh->GetGeomMap<TriGeom>().end(),
2445 "Failed to find LinMesh tri id");
2446 linMeshComp[c.first]->m_geomVec.push_back(it->second);
2447 }
2448 }
2449 break;
2450 case LibUtilities::eQuadrilateral:
2451 {
2452 int id = g->GetGlobalID();
2453
2454 if (useSimplex)
2455 {
2456 for (int i = 0; i < nsplit * nsplit; ++i)
2457 {
2458 int newid = 2 * (id * nsplit * nsplit + i);
2459
2460 auto it1 =
2461 m_linMesh->GetGeomMap<TriGeom>().find(newid);
2462 ASSERTL1(
2463 it1 != m_linMesh->GetGeomMap<TriGeom>().end(),
2464 "Failed to find LinMesh tri id (from quad)");
2465 linMeshComp[c.first]->m_geomVec.push_back(
2466 it1->second);
2467
2468 auto it2 = m_linMesh->GetGeomMap<TriGeom>().find(
2469 newid + 1);
2470 ASSERTL1(
2471 it2 != m_linMesh->GetGeomMap<TriGeom>().end(),
2472 "Failed to find LinMesh tri id (from quad)");
2473 linMeshComp[c.first]->m_geomVec.push_back(
2474 it2->second);
2475 }
2476 }
2477 else
2478 {
2479 for (int i = 0; i < nsplit * nsplit; ++i)
2480 {
2481 int newid = 2 * id * nsplit * nsplit + i;
2482
2483 auto it =
2484 m_linMesh->GetGeomMap<QuadGeom>().find(newid);
2485 ASSERTL1(
2486 it != m_linMesh->GetGeomMap<QuadGeom>().end(),
2487 "Failed to find LinMesh quad id");
2488 linMeshComp[c.first]->m_geomVec.push_back(
2489 it->second);
2490 }
2491 }
2492 }
2493 break;
2494 case LibUtilities::eTetrahedron:
2495 {
2496 int id = g->GetGlobalID();
2497 for (int i = 0; i < nsplit * nsplit * nsplit; ++i)
2498 {
2499 int newid = id * nsplit * nsplit * nsplit + i;
2500
2501 auto it = m_linMesh->GetGeomMap<TetGeom>().find(newid);
2502 ASSERTL1(it != m_linMesh->GetGeomMap<TetGeom>().end(),
2503 "Failed to find LinMesh tet id");
2504 linMeshComp[c.first]->m_geomVec.push_back(it->second);
2505 }
2506 }
2507 break;
2508 case LibUtilities::ePrism:
2509 {
2510 int id = g->GetGlobalID();
2511 for (int i = 0; i < nsplit * nsplit * nsplit; ++i)
2512 {
2513 int newid = id * nsplit * nsplit * nsplit + i;
2514
2515 auto it =
2516 m_linMesh->GetGeomMap<PrismGeom>().find(newid);
2517 ASSERTL1(it != m_linMesh->GetGeomMap<PrismGeom>().end(),
2518 "Failed to find LinMesh prism id");
2519 linMeshComp[c.first]->m_geomVec.push_back(it->second);
2520 }
2521 }
2522 break;
2523 default:
2524 NEKERROR(ErrorUtil::ewarning, "Shape Not set up ");
2525 break;
2526 }
2527 }
2528 }
2529
2530 // reset domain information based on new composites;
2531 for (auto &d : m_graph->GetDomain())
2532 {
2533 // make a list of composite ids
2534 std::string compstr;
2535 for (auto sd = d.second.begin(); sd != d.second.end();)
2536 {
2537 compstr += std::to_string(sd->first);
2538 compstr += (++sd != d.second.end()) ? "," : " ";
2539 }
2540 if (compstr.size()) // put in this check since in hdf5 format have empty
2541 // string cases
2542 {
2543 std::map<int, CompositeSharedPtr> unrollDomain;
2544 m_linMesh->GetCompositeList(compstr, unrollDomain);
2545 m_linMesh->GetDomain()[d.first] = unrollDomain;
2546 }
2547 }
2548
2549 // Setup Epxansion info as linear expansion
2550 m_linMesh->ReadExpansionInfo(SetupLinearExpansionType());
2551
2552 // loop over first expansion entry
2553 for (auto &expIt : *m_graph->GetExpansionInfoMap().begin()->second)
2554 {
2555 Geometry *geom = expIt.second->m_geomPtr;
2556 // expIt.second->m_geomPtr
2557 switch (expIt.second->m_geomPtr->GetShapeType())
2558 {
2559 case LibUtilities::Tri:
2560 {
2561 int elmtid = geom->GetGlobalID();
2562 ASSERTL1(CoeffMap2D[elmtid].size() ==
2563 (nsplit + 1) * (nsplit + 2) / 2,
2564 "CoeffMap2D is not the correct length");
2565
2566 for (int i = 0; i < nsplit * nsplit; ++i)
2567 {
2568 TriGeom *tri =
2569 m_linMesh->GetTriGeom(2 * elmtid * nsplit * nsplit + i);
2570
2571 std::vector<int> offset(3);
2572
2573 offset[0] = CoeffMap2D[elmtid][tri->GetVid(0)];
2574 offset[1] = CoeffMap2D[elmtid][tri->GetVid(2)];
2575 offset[2] = CoeffMap2D[elmtid][tri->GetVid(1)];
2576
2577 LinCoeffMap[2 * elmtid * nsplit * nsplit + i] =
2578 std::make_pair(elmtid, offset);
2579 }
2580 }
2581 break;
2582 case LibUtilities::Quad:
2583 if (useSimplex)
2584 {
2585 int elmtid = geom->GetGlobalID();
2586 for (int i = 0, e = 0; i < nsplit; ++i)
2587 {
2588 for (int j = 0; j < nsplit; ++j, ++e)
2589 {
2590 std::vector<int> offset(3);
2591 offset[0] = i * (nsplit + 1) + j;
2592 offset[1] = (i + 1) * (nsplit + 1) + j;
2593 offset[2] = i * (nsplit + 1) + j + 1;
2594
2595 LinCoeffMap[2 * (elmtid * nsplit * nsplit + e)] =
2596 std::make_pair(elmtid, offset);
2597
2598 std::vector<int> offset1(3);
2599 offset1[0] = (i + 1) * (nsplit + 1) + j + 1;
2600 offset1[1] = i * (nsplit + 1) + j + 1;
2601 offset1[2] = (i + 1) * (nsplit + 1) + j;
2602
2603 LinCoeffMap[2 * (elmtid * nsplit * nsplit + e) +
2604 1] = std::make_pair(elmtid, offset1);
2605 }
2606 }
2607 }
2608 else
2609 {
2610 int elmtid = geom->GetGlobalID();
2611 for (int i = 0, e = 0; i < nsplit; ++i)
2612 {
2613 for (int j = 0; j < nsplit; ++j, ++e)
2614 {
2615 std::vector<int> offset(4);
2616 offset[0] = i * (nsplit + 1) + j;
2617 offset[1] = i * (nsplit + 1) + j + 1;
2618 offset[2] = (i + 1) * (nsplit + 1) + j;
2619 offset[3] = (i + 1) * (nsplit + 1) + j + 1;
2620
2621 LinCoeffMap[2 * elmtid * nsplit * nsplit + e] =
2622 std::make_pair(elmtid, offset);
2623 }
2624 }
2625 }
2626 break;
2627 case LibUtilities::Tet:
2628 {
2629 int elmtid = geom->GetGlobalID();
2630 ASSERTL1(CoeffMap3D[elmtid].size() ==
2631 (nsplit + 1) * (nsplit + 2) * (nsplit + 3) / 6,
2632 "CoeffMap3D (Tet) is not the correct length");
2633
2634 for (int i = 0; i < nsplit * nsplit * nsplit; ++i)
2635 {
2636 TetGeom *tet = m_linMesh->GetTetGeom(
2637 elmtid * nsplit * nsplit * nsplit + i);
2638
2639 std::vector<int> offset(4);
2640
2641 offset[0] = CoeffMap3D[elmtid][tet->GetVid(0)];
2642 offset[1] = CoeffMap3D[elmtid][tet->GetVid(3)];
2643 offset[2] = CoeffMap3D[elmtid][tet->GetVid(2)];
2644 offset[3] = CoeffMap3D[elmtid][tet->GetVid(1)];
2645
2646 LinCoeffMap[elmtid * nsplit * nsplit * nsplit + i] =
2647 std::make_pair(elmtid, offset);
2648 }
2649 }
2650 break;
2651 case LibUtilities::Prism:
2652 {
2653 int elmtid = geom->GetGlobalID();
2654 ASSERTL1(CoeffMap3D[elmtid].size() ==
2655 (nsplit + 1) * (nsplit + 1) * (nsplit + 2) / 2,
2656 "CoeffMap3D (Prism) is not the correct length");
2657
2658 for (int i = 0; i < nsplit * nsplit * nsplit; ++i)
2659 {
2660 PrismGeom *prism = m_linMesh->GetPrismGeom(
2661 elmtid * nsplit * nsplit * nsplit + i);
2662
2663 std::vector<int> offset(6);
2664
2665 offset[0] = CoeffMap3D[elmtid][prism->GetVid(0)];
2666 offset[1] = CoeffMap3D[elmtid][prism->GetVid(4)];
2667 offset[2] = CoeffMap3D[elmtid][prism->GetVid(3)];
2668 offset[3] = CoeffMap3D[elmtid][prism->GetVid(5)];
2669 offset[4] = CoeffMap3D[elmtid][prism->GetVid(1)];
2670 offset[5] = CoeffMap3D[elmtid][prism->GetVid(2)];
2671
2672 LinCoeffMap[elmtid * nsplit * nsplit * nsplit + i] =
2673 std::make_pair(elmtid, offset);
2674 }
2675 }
2676 break;
2677 default:
2678 NEKERROR(ErrorUtil::efatal, "Need to setup local mapping");
2679 break;
2680 }
2681 }
2682}
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition ErrorUtil.hpp:250
Nektar::SpatialDomains::CompositeSharedPtr
std::shared_ptr< Composite > CompositeSharedPtr
Definition MeshGraph.h:185
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition StdRegions.hpp:539
Nektar::UnitTests::d
std::vector< double > d(NPUPPER *NPUPPER)

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Nektar::ErrorUtil::efatal, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, Nektar::ErrorUtil::ewarning, Nektar::SpatialDomains::Geometry::GetGlobalID(), Nektar::SpatialDomains::Geometry::GetVid(), m_graph, m_linMesh, NEKERROR, Nektar::LibUtilities::Prism, Nektar::LibUtilities::Quad, SetupLinearExpansionType(), Nektar::LibUtilities::Tet, and Nektar::LibUtilities::Tri.

Referenced by CreateLinearGraph().

◆ LinMeshSetUpPrismGeom()

void Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpPrismGeom ( int  nsplit,
std::map< int, int > &  FceEdgOffset,
std::map< int, std::map< int, int > > &  CoeffMap,
bool  UseGLL 
)
private

Definition at line 1819 of file LinearMeshGraph.hpp.

1822{
1823 // get vertex ids on faces of 3D shape
1824 Array3D<int> E_vid(nsplit + 1, nsplit + 1, nsplit + 1);
1825
1826 // X-aligned edge
1827 Array3D<int> E_eidx(nsplit + 1, nsplit + 1, nsplit + 1);
1828 // Y-aligned edge
1829 Array3D<int> E_eidy(nsplit + 1, nsplit + 1, nsplit + 1);
1830 // Z-aligned edge
1831 Array3D<int> E_eidz(nsplit + 1, nsplit + 1, nsplit + 1);
1832
1833 // XZ-diagonal edge
1834 Array3D<int> E_eidxz(nsplit + 1, nsplit + 1, nsplit + 1);
1835
1836 // xy face ids
1837 Array3D<int> E_fidxy(nsplit + 1, nsplit + 1, nsplit + 1);
1838 // xz face ids (lower triangular if simplex)
1839 Array3D<int> E_fidxz1(nsplit + 1, nsplit + 1, nsplit + 1);
1840 // xz face ids (upper triangular if simplex)
1841 Array3D<int> E_fidxz2(nsplit + 1, nsplit + 1, nsplit + 1);
1842 // yz face ids
1843 Array3D<int> E_fidyz(nsplit + 1, nsplit + 1, nsplit + 1);
1844
1845 // Fill with yz diagonal face ids
1846 Array3D<int> E_fidyzd(nsplit + 1, nsplit + 1, nsplit + 1);
1847
1848 int maxvertid, maxedgeid;
1849 LibUtilities::CommSharedPtr comm = m_session->GetComm();
1850
1851 Array<OneD, NekDouble> Xpts, Ypts, Zpts;
1852 Xpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1853 Ypts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1854 Zpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1855
1856 /* reset maxedgeid and maxvertid */
1857 maxvertid = m_linMesh->GetGeomMap<PointGeom>().rbegin()->first + 1;
1858 comm->AllReduce(maxvertid, LibUtilities::ReduceMax);
1859 maxedgeid = m_linMesh->GetGeomMap<SegGeom>().rbegin()->first + 1;
1860 comm->AllReduce(maxedgeid, LibUtilities::ReduceMax);
1861
1862 int maxfaceid = 0;
1863 if (m_linMesh->GetGeomMap<TriGeom>().size())
1864 {
1865 maxfaceid = std::max(
1866 maxfaceid, m_linMesh->GetGeomMap<TriGeom>().rbegin()->first + 1);
1867 }
1868 if (m_linMesh->GetGeomMap<QuadGeom>().size())
1869 {
1870 maxfaceid = std::max(
1871 maxfaceid, m_linMesh->GetGeomMap<QuadGeom>().rbegin()->first + 1);
1872 }
1873
1874 comm->AllReduce(maxfaceid, LibUtilities::ReduceMax);
1875
1876 // define verte reorientaiton - might need for all square faces
1877 std::map<StdRegions::Orientation, std::array<int, 4>> vertMap{
1878 {StdRegions::eDir1FwdDir1_Dir2FwdDir2, {{0, 1, 2, 3}}},
1879 {StdRegions::eDir1BwdDir1_Dir2FwdDir2, {{1, 0, 3, 2}}},
1880 {StdRegions::eDir1FwdDir1_Dir2BwdDir2, {{3, 2, 1, 0}}},
1881 {StdRegions::eDir1BwdDir1_Dir2BwdDir2, {{2, 3, 0, 1}}},
1882 {StdRegions::eDir1FwdDir2_Dir2FwdDir1, {{0, 3, 2, 1}}},
1883 {StdRegions::eDir1BwdDir2_Dir2FwdDir1, {{3, 0, 1, 2}}},
1884 {StdRegions::eDir1FwdDir2_Dir2BwdDir1, {{1, 2, 3, 0}}},
1885 {StdRegions::eDir1BwdDir2_Dir2BwdDir1, {{2, 1, 0, 3}}}};
1886
1887 // define edge reorientaiton - might need for all square faces
1888 std::map<StdRegions::Orientation, std::array<int, 4>> edgMap{
1889 {StdRegions::eDir1FwdDir1_Dir2FwdDir2, {{0, 1, 2, 3}}},
1890 {StdRegions::eDir1BwdDir1_Dir2FwdDir2, {{0, 3, 2, 1}}},
1891 {StdRegions::eDir1FwdDir1_Dir2BwdDir2, {{2, 1, 0, 3}}},
1892 {StdRegions::eDir1BwdDir1_Dir2BwdDir2, {{2, 3, 0, 1}}},
1893 {StdRegions::eDir1FwdDir2_Dir2FwdDir1, {{3, 2, 1, 0}}},
1894 {StdRegions::eDir1BwdDir2_Dir2FwdDir1, {{3, 0, 1, 2}}},
1895 {StdRegions::eDir1FwdDir2_Dir2BwdDir1, {{1, 2, 3, 0}}},
1896 {StdRegions::eDir1BwdDir2_Dir2BwdDir1, {{1, 0, 3, 2}}}};
1897
1898 //--------------------------------------------------------
1899 // Split prism elements of macro mesh into new edges and vertices
1900 //--------------------------------------------------------
1901
1902 int edgeoffset = 3 * nsplit * (nsplit - 1) * (nsplit - 1);
1903 int faceoffset = 3 * nsplit * nsplit * (nsplit - 1);
1904
1905 // For every original elmt keep a map between the vid of the split
1906 // tet and the location of the points in the macro nodes to create
1907 // the LinCoeffMap when required. Needed since we reshuffle tri faces
1908 // orientaiton so cannot infer pattern directly
1909 for (auto [id, geom] : m_graph->GetGeomMap<PrismGeom>())
1910 {
1911 int vcnt = 0;
1912 int ecnt = 0;
1913 int fcnt = 0;
1914 int elmtid = geom->GetGlobalID();
1915
1916 StdRegions::StdExpansionSharedPtr Xmap = geom->GetXmap();
1917
1918 /* face 0 - Quad */
1919 int fid = geom->GetFace(0)->GetGlobalID();
1920 int fce_offset = 2 * fid * nsplit * nsplit;
1921
1922 StdRegions::Orientation orient = geom->GetForient(0);
1923 Array<OneD, int> idmap;
1924 // macro split reorientation
1925 Xmap->ReOrientTracePhysMap(orient, idmap, nsplit, nsplit, false);
1926
1927 int cnt = 0;
1928 for (int j = 0; j < nsplit; ++j)
1929 {
1930 for (int i = 0; i < nsplit; ++i)
1931 {
1932 int Fce_id = fce_offset + idmap[i + j * nsplit];
1933
1934 E_fidxy(0, j, i) = Fce_id;
1935
1936 QuadGeom *face = m_linMesh->GetQuadGeom(Fce_id);
1937
1938 // fill in vertices & edges from faces1
1939 E_vid(0, j, i) = face->GetVid(vertMap[orient][0]);
1940 E_vid(0, j, i + 1) = face->GetVid(vertMap[orient][1]);
1941 E_vid(0, j + 1, i) = face->GetVid(vertMap[orient][3]);
1942 E_vid(0, j + 1, i + 1) = face->GetVid(vertMap[orient][2]);
1943
1944 E_eidx(0, j, i) = face->GetEid(edgMap[orient][0]);
1945 E_eidx(0, j + 1, i) = face->GetEid(edgMap[orient][2]);
1946 E_eidy(0, j, i) = face->GetEid(edgMap[orient][3]);
1947 E_eidy(0, j, i + 1) = face->GetEid(edgMap[orient][1]);
1948 }
1949 }
1950
1951 /* face 1 - Tri */
1952 fid = geom->GetFace(1)->GetGlobalID();
1953 fce_offset = 2 * fid * nsplit * nsplit;
1954 bool FceFwd =
1955 geom->GetForient(1) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
1956 cnt = 0;
1957 for (int j = 0; j < nsplit; ++j)
1958 {
1959 for (int i = 0; i < nsplit - j; ++i)
1960 {
1961 int Fce_id = FceFwd ? fce_offset + cnt + i
1962 : fce_offset + cnt + nsplit - j - 1 - i;
1963
1964 E_fidxz1(j, 0, i) = Fce_id;
1965
1966 // due to edge shuffling earlier need to identify
1967 // know
1968 // orientation (singlar vertex at top) orientation to
1969 // fill in stored edges and vertices
1970 int oset = FceEdgOffset[Fce_id];
1971
1972 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
1973 // fill in vertices & edges from faces
1974 if (FceFwd)
1975 {
1976 E_vid(j, 0, i) = face->GetVid(oset);
1977 E_vid(j, 0, i + 1) = face->GetVid((oset + 1) % 3);
1978
1979 E_eidz(j, 0, i) = face->GetEid((oset + 2) % 3);
1980 E_eidxz(j, 0, i) = face->GetEid((oset + 1) % 3);
1981 }
1982 else
1983 {
1984 E_vid(j, 0, i) = face->GetVid((oset + 1) % 3);
1985 E_vid(j, 0, i + 1) = face->GetVid(oset);
1986
1987 E_eidz(j, 0, i) = face->GetEid((oset + 1) % 3);
1988 E_eidxz(j, 0, i) = face->GetEid((oset + 2) % 3);
1989 }
1990 E_vid(j + 1, 0, i) = face->GetVid((oset + 2) % 3);
1991 E_eidx(j, 0, i) = face->GetEid(oset);
1992 }
1993 cnt += nsplit - j;
1994
1995 for (int i = 0; i < nsplit - j - 1; ++i)
1996 {
1997 int Fce_id = FceFwd ? fce_offset + cnt + i
1998 : fce_offset + cnt + nsplit - j - 2 - i;
1999 E_fidxz2(j, 0, i) = Fce_id;
2000 }
2001 cnt += nsplit - j - 1;
2002 }
2003
2004 /* face 2 - Quad on diagonal */
2005 fid = geom->GetFace(2)->GetGlobalID();
2006 fce_offset = 2 * fid * nsplit * nsplit;
2007 orient = geom->GetForient(2);
2008
2009 // macro split reorientation
2010 Xmap->ReOrientTracePhysMap(orient, idmap, nsplit, nsplit, false);
2011
2012 cnt = 0;
2013 for (int j = 0; j < nsplit; ++j)
2014 {
2015 for (int i = 0; i < nsplit; ++i)
2016 {
2017 int Fce_id = fce_offset + idmap[i + j * nsplit];
2018
2019 E_fidyzd(j, i, nsplit - 1 - j) = Fce_id;
2020
2021 QuadGeom *face = m_linMesh->GetQuadGeom(Fce_id);
2022
2023 // fill in vertices & edges from faces1
2024 E_vid(j, i, nsplit - j) = face->GetVid(vertMap[orient][0]);
2025 E_vid(j, i + 1, nsplit - j) = face->GetVid(vertMap[orient][1]);
2026 E_vid(j + 1, i, nsplit - 1 - j) =
2027 face->GetVid(vertMap[orient][3]);
2028 E_vid(j + 1, i + 1, nsplit - 1 - j) =
2029 face->GetVid(vertMap[orient][2]);
2030
2031 E_eidy(j, i, nsplit - j) = face->GetEid(edgMap[orient][0]);
2032 E_eidy(j + 1, i, nsplit - 1 - j) =
2033 face->GetEid(edgMap[orient][2]);
2034 E_eidxz(j, i, nsplit - 1 - j) = face->GetEid(edgMap[orient][3]);
2035 E_eidxz(j, i + 1, nsplit - 1 - j) =
2036 face->GetEid(edgMap[orient][1]);
2037 }
2038 }
2039 /* face 3 - Tri */
2040 fid = geom->GetFace(3)->GetGlobalID();
2041 fce_offset = 2 * fid * nsplit * nsplit;
2042 FceFwd = geom->GetForient(3) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
2043 cnt = 0;
2044 for (int j = 0; j < nsplit; ++j)
2045 {
2046 for (int i = 0; i < nsplit - j; ++i)
2047 {
2048 int Fce_id = FceFwd ? fce_offset + cnt + i
2049 : fce_offset + cnt + nsplit - j - 1 - i;
2050
2051 E_fidxz1(j, nsplit, i) = Fce_id;
2052
2053 // due to edge shuffling earlier need to identify
2054 // know
2055 // orientation (singlar vertex at top) orientation to
2056 // fill in stored edges and vertices
2057 int oset = FceEdgOffset[Fce_id];
2058
2059 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
2060 // fill in vertices & edges from faces
2061 if (FceFwd)
2062 {
2063 E_vid(j, nsplit, i) = face->GetVid(oset);
2064 E_vid(j, nsplit, i + 1) = face->GetVid((oset + 1) % 3);
2065
2066 E_eidz(j, nsplit, i) = face->GetEid((oset + 2) % 3);
2067 E_eidxz(j, nsplit, i) = face->GetEid((oset + 1) % 3);
2068 }
2069 else
2070 {
2071 E_vid(j, nsplit, i) = face->GetVid((oset + 1) % 3);
2072 E_vid(j, nsplit, i + 1) = face->GetVid(oset);
2073
2074 E_eidz(j, nsplit, i) = face->GetEid((oset + 1) % 3);
2075 E_eidxz(j, nsplit, i) = face->GetEid((oset + 2) % 3);
2076 }
2077 E_vid(j + 1, nsplit, i) = face->GetVid((oset + 2) % 3);
2078 E_eidx(j, nsplit, i) = face->GetEid(oset);
2079 }
2080 cnt += nsplit - j;
2081
2082 for (int i = 0; i < nsplit - j - 1; ++i)
2083 {
2084 int Fce_id = FceFwd ? fce_offset + cnt + i
2085 : fce_offset + cnt + nsplit - j - 2 - i;
2086 E_fidxz2(j, nsplit, i) = Fce_id;
2087 }
2088 cnt += nsplit - j - 1;
2089 }
2090
2091 /* face 4 - Quad */
2092 fid = geom->GetFace(4)->GetGlobalID();
2093 fce_offset = 2 * fid * nsplit * nsplit;
2094
2095 orient = geom->GetForient(4);
2096 // macro split reorientation
2097 Xmap->ReOrientTracePhysMap(orient, idmap, nsplit, nsplit, false);
2098
2099 cnt = 0;
2100 for (int j = 0; j < nsplit; ++j)
2101 {
2102 for (int i = 0; i < nsplit; ++i)
2103 {
2104 int Fce_id = fce_offset + idmap[i + j * nsplit];
2105
2106 E_fidyz(j, i, 0) = Fce_id;
2107
2108 QuadGeom *face = m_linMesh->GetQuadGeom(Fce_id);
2109
2110 // fill in vertices & edges from faces1
2111 E_vid(j, i, 0) = face->GetVid(vertMap[orient][0]);
2112 E_vid(j, i + 1, 0) = face->GetVid(vertMap[orient][1]);
2113 E_vid(j + 1, i, 0) = face->GetVid(vertMap[orient][3]);
2114 E_vid(j + 1, i + 1, 0) = face->GetVid(vertMap[orient][2]);
2115
2116 E_eidy(j, i, 0) = face->GetEid(edgMap[orient][0]);
2117 E_eidy(j + 1, i, 0) = face->GetEid(edgMap[orient][2]);
2118 E_eidz(j, i, 0) = face->GetEid(edgMap[orient][3]);
2119 E_eidz(j, i + 1, 0) = face->GetEid(edgMap[orient][1]);
2120 }
2121 }
2122
2123 Array<OneD, NekDouble> tmp(Xmap->GetTotPoints());
2124
2125 // fill in interior vids if exist
2126 if (UseGLL) // get GLL points
2127 {
2128 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
2129 Xmap->PhysInterpToGLL(tmp, Xpts, nsplit + 1);
2130 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
2131 Xmap->PhysInterpToGLL(tmp, Ypts, nsplit + 1);
2132 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
2133 Xmap->PhysInterpToGLL(tmp, Zpts, nsplit + 1);
2134 }
2135 else /* get equispaced points */
2136 {
2137 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
2138 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Xpts, nsplit + 1);
2139 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
2140 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Ypts, nsplit + 1);
2141 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
2142 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Zpts, nsplit + 1);
2143 }
2144
2145 cnt = (nsplit + 1) * (nsplit + 1);
2146 for (int k = 1; k < nsplit; ++k)
2147 {
2148 cnt += nsplit + 1 - k; // first row
2149 for (int j = 1; j < nsplit; ++j)
2150 {
2151 cnt++;
2152 for (int i = 1; i < nsplit - k; ++i, cnt++)
2153 {
2154 int vid =
2155 maxvertid +
2156 elmtid * (nsplit - 1) * (nsplit - 1) * (nsplit - 1) +
2157 vcnt++;
2158
2159 E_vid(k, j, i) = vid;
2160
2161 PointGeomUniquePtr vert =
2162 ObjPoolManager<PointGeom>::AllocateUniquePtr(
2163 m_spaceDimension, vid, Xpts[cnt], Ypts[cnt],
2164 Zpts[cnt]);
2165
2166 m_linMesh->AddGeom(vid, std::move(vert));
2167 }
2168 cnt++;
2169 }
2170 cnt += nsplit + 1 - k; // last row
2171 }
2172
2173 // declare Prisms
2174 cnt = 0;
2175 for (int k = 0; k < nsplit; ++k)
2176 {
2177 for (int j = 0; j < nsplit; ++j)
2178 {
2179 for (int i = 0; i < nsplit - k; ++i)
2180 {
2181
2182 if ((i < nsplit - 1 - k) && (j < nsplit - 1))
2183 {
2184 /* add three new edges */
2185
2186 /* new diagonals */
2187 int edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
2188 std::array<PointGeom *, 2> vert;
2189
2190 vert[0] =
2191 m_linMesh->GetPointGeom(E_vid(k, j + 1, i + 1));
2192 vert[1] =
2193 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i));
2194 SegGeomUniquePtr edge0 =
2195 ObjPoolManager<SpatialDomains::SegGeom>::
2196 AllocateUniquePtr(edgid, m_spaceDimension,
2197 vert);
2198 m_linMesh->AddGeom(edgid, std::move(edge0));
2199 E_eidxz(k, j + 1, i) = edgid;
2200
2201 // add in vertical edge
2202 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
2203 vert[0] =
2204 m_linMesh->GetPointGeom(E_vid(k, j + 1, i + 1));
2205 vert[1] =
2206 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
2207 SegGeomUniquePtr edge1 =
2208 ObjPoolManager<SpatialDomains::SegGeom>::
2209 AllocateUniquePtr(edgid, m_spaceDimension,
2210 vert);
2211 m_linMesh->AddGeom(edgid, std::move(edge1));
2212 E_eidz(k, j + 1, i + 1) = edgid;
2213
2214 // horizontal edge
2215 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
2216 vert[0] =
2217 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
2218 vert[1] =
2219 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i));
2220 SegGeomUniquePtr edge2 =
2221 ObjPoolManager<SpatialDomains::SegGeom>::
2222 AllocateUniquePtr(edgid, m_spaceDimension,
2223 vert);
2224 m_linMesh->AddGeom(edgid, std::move(edge2));
2225 E_eidx(k + 1, j + 1, i) = edgid;
2226
2227 /* Add in new faces */
2228 /* Far end Tri */
2229 int newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
2230 std::array<SegGeom *, 3>
2231 edges; // BUG KK - it was 4 originally ???
2232 std::array<SegGeom *, 3> edgessort;
2233
2234 edgessort[0] =
2235 m_linMesh->GetSegGeom(E_eidx(k + 1, 0, i));
2236 edgessort[1] = m_linMesh->GetSegGeom(E_eidxz(k, 0, i));
2237 edgessort[2] =
2238 m_linMesh->GetSegGeom(E_eidz(k, 0, i + 1));
2239 edges[0] =
2240 m_linMesh->GetSegGeom(E_eidx(k + 1, j + 1, i));
2241 edges[1] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
2242 edges[2] =
2243 m_linMesh->GetSegGeom(E_eidz(k, j + 1, i + 1));
2244 // shuffle so lowest vid is
2245 // on edges 1, 2
2246 TriShuffleEdges(edges, &edgessort); // BUG KK ?
2247 TriGeomUniquePtr tri1 =
2248 ObjPoolManager<SpatialDomains::TriGeom>::
2249 AllocateUniquePtr(newtriid, edges);
2250 m_linMesh->AddGeom(newtriid, std::move(tri1));
2251
2252 E_fidxz2(k, j + 1, i) = newtriid;
2253 }
2254
2255 if (i < (nsplit - 2 - k))
2256 {
2257 /* new y-parallel edge */
2258 int edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
2259 std::array<PointGeom *, 2> vert;
2260
2261 vert[0] =
2262 m_linMesh->GetPointGeom(E_vid(k + 1, j, i + 1));
2263 vert[1] =
2264 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
2265 SegGeomUniquePtr edge0 =
2266 ObjPoolManager<SpatialDomains::SegGeom>::
2267 AllocateUniquePtr(edgid, m_spaceDimension,
2268 vert);
2269 m_linMesh->AddGeom(edgid, std::move(edge0));
2270 E_eidy(k + 1, j, i + 1) = edgid;
2271 }
2272
2273 if (i < (nsplit - 1 - k))
2274 {
2275 std::array<SegGeom *, 4> edges;
2276 /* diagonal face */
2277 int newquadid =
2278 elmtid * faceoffset + maxfaceid + fcnt++;
2279 edges[0] = m_linMesh->GetSegGeom(E_eidy(k, j, i + 1));
2280 edges[1] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
2281 edges[2] = m_linMesh->GetSegGeom(E_eidy(k + 1, j, i));
2282 edges[3] = m_linMesh->GetSegGeom(E_eidxz(k, j, i));
2283
2284 QuadGeomUniquePtr quad0 =
2285 ObjPoolManager<SpatialDomains::QuadGeom>::
2286 AllocateUniquePtr(newquadid, edges);
2287 m_linMesh->AddGeom(newquadid, std::move(quad0));
2288 E_fidyzd(k, j, i) = newquadid;
2289
2290 /* vertical face */
2291 newquadid = elmtid * faceoffset + maxfaceid + fcnt++;
2292 edges[0] = m_linMesh->GetSegGeom(E_eidy(k, j, i + 1));
2293 edges[1] =
2294 m_linMesh->GetSegGeom(E_eidz(k, j + 1, i + 1));
2295 edges[2] =
2296 m_linMesh->GetSegGeom(E_eidy(k + 1, j, i + 1));
2297 edges[3] = m_linMesh->GetSegGeom(E_eidz(k, j, i + 1));
2298
2299 QuadGeomUniquePtr quad1 =
2300 ObjPoolManager<SpatialDomains::QuadGeom>::
2301 AllocateUniquePtr(newquadid, edges);
2302 m_linMesh->AddGeom(newquadid, std::move(quad1));
2303 E_fidyz(k, j, i + 1) = newquadid;
2304
2305 /* horizontal face */
2306 newquadid = elmtid * faceoffset + maxfaceid + fcnt++;
2307 edges[0] =
2308 m_linMesh->GetSegGeom(E_eidy(k + 1, j, i + 1));
2309 edges[1] =
2310 m_linMesh->GetSegGeom(E_eidx(k + 1, j + 1, i));
2311 edges[2] = m_linMesh->GetSegGeom(E_eidy(k + 1, j, i));
2312 edges[3] = m_linMesh->GetSegGeom(E_eidx(k + 1, j, i));
2313
2314 QuadGeomUniquePtr quad2 =
2315 ObjPoolManager<SpatialDomains::QuadGeom>::
2316 AllocateUniquePtr(newquadid, edges);
2317 m_linMesh->AddGeom(newquadid, std::move(quad2));
2318 E_fidxy(k + 1, j, i) = newquadid;
2319
2320 // /* top Corner Prism */
2321 int newprismid =
2322 elmtid * nsplit * nsplit * nsplit + cnt++;
2323
2324 std::array<Geometry2D *, 5> faces;
2325 faces[0] = m_linMesh->GetQuadGeom(E_fidxy(k + 1, j, i));
2326 faces[1] = m_linMesh->GetTriGeom(E_fidxz2(k, j, i));
2327 faces[2] = m_linMesh->GetQuadGeom(E_fidyzd(k, j, i));
2328 faces[3] = m_linMesh->GetTriGeom(E_fidxz2(k, j + 1, i));
2329 faces[4] = m_linMesh->GetQuadGeom(E_fidyz(k, j, i + 1));
2330 PrismShuffleFaces(faces); // ensure faces are
2331 // ordered correctly
2332 PrismGeomUniquePtr prism =
2333 ObjPoolManager<SpatialDomains::PrismGeom>::
2334 AllocateUniquePtr(newprismid, faces);
2335 m_linMesh->AddGeom(newprismid, std::move(prism));
2336 m_linMesh->PopulateFaceToElMap(
2337 m_linMesh->GetPrismGeom(newprismid), 5);
2338 }
2339
2340 if (j < nsplit - 1)
2341 {
2342 int newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
2343 std::array<SegGeom *, 3> edges;
2344 std::array<SegGeom *, 3> edgessort;
2345
2346 edgessort[0] = m_linMesh->GetSegGeom(E_eidx(k, 0, i));
2347 edgessort[1] = m_linMesh->GetSegGeom(E_eidxz(k, 0, i));
2348 edgessort[2] = m_linMesh->GetSegGeom(E_eidz(k, 0, i));
2349 edges[0] = m_linMesh->GetSegGeom(E_eidx(k, j + 1, i));
2350 edges[1] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
2351 edges[2] = m_linMesh->GetSegGeom(E_eidz(k, j + 1, i));
2352 // shuffle so lowest vid is on
2353 // edges 1, 2
2354 TriShuffleEdges(edges, &edgessort);
2355 TriGeomUniquePtr tri0 =
2356 ObjPoolManager<SpatialDomains::TriGeom>::
2357 AllocateUniquePtr(newtriid, edges);
2358 m_linMesh->AddGeom(newtriid, std::move(tri0));
2359 E_fidxz1(k, j + 1, i) = newtriid;
2360 }
2361
2362 /* bottom Corner Prism */
2363 int newprismid = elmtid * nsplit * nsplit * nsplit + cnt++;
2364
2365 std::array<Geometry2D *, 5> faces;
2366 faces[0] = m_linMesh->GetQuadGeom(E_fidxy(k, j, i));
2367 faces[1] = m_linMesh->GetTriGeom(E_fidxz1(k, j, i));
2368 faces[2] = m_linMesh->GetQuadGeom(E_fidyzd(k, j, i));
2369 faces[3] = m_linMesh->GetTriGeom(E_fidxz1(k, j + 1, i));
2370 faces[4] = m_linMesh->GetQuadGeom(E_fidyz(k, j, i));
2371 PrismShuffleFaces(faces); // ensure faces are
2372 // ordered correctly
2373 PrismGeomUniquePtr prism =
2374 ObjPoolManager<SpatialDomains::PrismGeom>::
2375 AllocateUniquePtr(newprismid, faces);
2376 m_linMesh->AddGeom(newprismid, std::move(prism));
2377
2378 m_linMesh->PopulateFaceToElMap(
2379 m_linMesh->GetPrismGeom(newprismid), 5);
2380 }
2381 }
2382 }
2383 ASSERTL1(
2384 ecnt <= edgeoffset,
2385 "This should not happen, edgeoffset might need to be increased");
2386 ASSERTL1(
2387 fcnt <= faceoffset,
2388 "This should not happen, faceoffset might need to be increased");
2389
2390 // collect a mapping of Vid to point vertex location
2391 cnt = 0;
2392 for (int k = 0; k < nsplit + 1; ++k)
2393 {
2394 for (int j = 0; j < nsplit + 1; ++j)
2395 {
2396 for (int i = 0; i < nsplit + 1 - k; ++i)
2397 {
2398 CoeffMap[elmtid][E_vid(k, j, i)] = cnt++;
2399 }
2400 }
2401 }
2402 }
2403}
Nektar::SpatialDomains::LinearMeshGraph::PrismShuffleFaces
void PrismShuffleFaces(std::array< Geometry2D *, 5 > &faces)
Definition LinearMeshGraph.hpp:294
Nektar::SpatialDomains::PrismGeomUniquePtr
unique_ptr_objpool< PrismGeom > PrismGeomUniquePtr
Definition MeshGraph.h:107

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), ASSERTL1, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eDir1FwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1FwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1FwdDir2_Dir2FwdDir1, Nektar::SpatialDomains::Geometry::GetEid(), Nektar::SpatialDomains::Geometry::GetVid(), m_graph, m_linMesh, m_session, m_spaceDimension, PrismShuffleFaces(), Nektar::LibUtilities::ReduceMax, and TriShuffleEdges().

Referenced by CreateLinearGraph().

◆ LinMeshSetUpTetGeom()

void Nektar::SpatialDomains::LinearMeshGraph::LinMeshSetUpTetGeom ( int  nsplit,
std::map< int, int > &  FceEdgOffset,
std::map< int, std::map< int, int > > &  CoeffMap,
bool  UseGLL 
)
private

Definition at line 1064 of file LinearMeshGraph.hpp.

1067{
1068 // get vertex ids on faces of 3D shape
1069 Array3D<int> E_vid(nsplit + 1, nsplit + 1, nsplit + 1);
1070 // X-aligned edge
1071 Array3D<int> E_eidx(nsplit + 1, nsplit + 1, nsplit + 1);
1072 // Y-aligned edge
1073 Array3D<int> E_eidy(nsplit + 1, nsplit + 1, nsplit + 1);
1074 // Z-aligned edge
1075 Array3D<int> E_eidz(nsplit + 1, nsplit + 1, nsplit + 1);
1076
1077 // XY-diagonal edge
1078 Array3D<int> E_eidxy(nsplit + 1, nsplit + 1, nsplit + 1);
1079 // ZY-diagonal edge
1080 Array3D<int> E_eidyz(nsplit + 1, nsplit + 1, nsplit + 1);
1081 // XZ-diagonal edge
1082 Array3D<int> E_eidxz(nsplit + 1, nsplit + 1, nsplit + 1);
1083
1084 // Fill with xy face ids (lower triangular if simplex)
1085 Array3D<int> E_fidxy1(nsplit + 1, nsplit + 1, nsplit + 1);
1086 // Fill with xy face ids (upperr triangular if simplex)
1087 Array3D<int> E_fidxy2(nsplit + 1, nsplit + 1, nsplit + 1);
1088
1089 // Fill with xz face ids (lower triangular if simplex)
1090 Array3D<int> E_fidxz1(nsplit + 1, nsplit + 1, nsplit + 1);
1091 // Fill with xz face ids (upperr triangular if simplex)
1092 Array3D<int> E_fidxz2(nsplit + 1, nsplit + 1, nsplit + 1);
1093
1094 // Fill with yz face ids (lower triangular if simplex)
1095 Array3D<int> E_fidyz1(nsplit + 1, nsplit + 1, nsplit + 1);
1096 // Fill with yz face ids (upperr triangular if simplex)
1097 Array3D<int> E_fidyz2(nsplit + 1, nsplit + 1, nsplit + 1);
1098
1099 // Fill with yz diagonal face ids (lower triangular if simplex)
1100 Array3D<int> E_fidyzd1(nsplit + 1, nsplit + 1, nsplit + 1);
1101 // Fill with yz diagonal face ids (upperr triangular if simplex)
1102 Array3D<int> E_fidyzd2(nsplit + 1, nsplit + 1, nsplit + 1);
1103
1104 int maxvertid, maxedgeid;
1105 LibUtilities::CommSharedPtr comm = m_session->GetComm();
1106
1107 Array<OneD, NekDouble> Xpts, Ypts, Zpts;
1108 Xpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1109 Ypts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1110 Zpts = Array<OneD, NekDouble>((nsplit + 1) * (nsplit + 1) * (nsplit + 1));
1111
1112 /* reset maxedgeid and maxvertid */
1113 maxvertid = m_linMesh->GetGeomMap<PointGeom>().rbegin()->first + 1;
1114 comm->AllReduce(maxvertid, LibUtilities::ReduceMax);
1115 maxedgeid = m_linMesh->GetGeomMap<SegGeom>().rbegin()->first + 1;
1116 comm->AllReduce(maxedgeid, LibUtilities::ReduceMax);
1117
1118 int maxfaceid = 0;
1119
1120 if (m_linMesh->GetGeomMap<TriGeom>().size())
1121 {
1122 maxfaceid = std::max(
1123 maxfaceid, m_linMesh->GetGeomMap<TriGeom>().rbegin()->first + 1);
1124 }
1125 if (m_linMesh->GetGeomMap<QuadGeom>().size())
1126 {
1127 maxfaceid = std::max(
1128 maxfaceid, m_linMesh->GetGeomMap<QuadGeom>().rbegin()->first + 1);
1129 }
1130
1131 comm->AllReduce(maxfaceid, LibUtilities::ReduceMax);
1132 //--------------------------------------------------------
1133 // Split tet elements of macro mesh into new edges and vertices
1134 //--------------------------------------------------------
1135
1136 int edgeoffset = 3 * nsplit * (nsplit - 1) * (nsplit - 1);
1137 int faceoffset = 3 * nsplit * nsplit * (nsplit - 1);
1138
1139 // For every original tet keep a map between the vid of the split
1140 // tet and the location of the points in the macro nodes to create
1141 // the LinCoeffMap when required. Needed since we reshuffle tet
1142 // orientaiton so cannot infer pattern directly
1143 for (auto [id, geom] : m_graph->GetGeomMap<TetGeom>())
1144 {
1145 int vcnt = 0;
1146 int ecnt = 0;
1147 int fcnt = 0;
1148 int elmtid = geom->GetGlobalID();
1149
1150 /* face 0 */
1151 int fid = geom->GetFace(0)->GetGlobalID();
1152 int fce_offset = 2 * fid * nsplit * nsplit;
1153 bool FceFwd =
1154 geom->GetForient(0) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
1155
1156 int cnt = 0;
1157 for (int j = 0; j < nsplit; ++j)
1158 {
1159 for (int i = 0; i < nsplit - j; ++i)
1160 {
1161 int Fce_id = FceFwd ? fce_offset + cnt + i
1162 : fce_offset + cnt + nsplit - j - 1 - i;
1163 E_fidxy1(0, j, i) = Fce_id;
1164
1165 // due to edge shuffling earlier need to find new
1166 // orientation (singlar vertex at top) to fill in
1167 // stored edges and vertices
1168 int oset = FceEdgOffset[Fce_id];
1169
1170 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
1171 // fill in vertices & edges from faces1
1172 if (FceFwd)
1173 {
1174 E_vid(0, j, i) = face->GetVid(oset);
1175 E_vid(0, j, i + 1) = face->GetVid((oset + 1) % 3);
1176
1177 E_eidy(0, j, i) = face->GetEid((oset + 2) % 3);
1178 E_eidxy(0, j, i) = face->GetEid((oset + 1) % 3);
1179 }
1180 else
1181 {
1182 E_vid(0, j, i) = face->GetVid((oset + 1) % 3);
1183 E_vid(0, j, i + 1) = face->GetVid(oset);
1184
1185 E_eidy(0, j, i) = face->GetEid((oset + 1) % 3);
1186 E_eidxy(0, j, i) = face->GetEid((oset + 2) % 3);
1187 }
1188 E_vid(0, j + 1, i) = face->GetVid((oset + 2) % 3);
1189 E_eidx(0, j, i) = face->GetEid(oset);
1190 }
1191 cnt += nsplit - j;
1192
1193 for (int i = 0; i < nsplit - j - 1; ++i)
1194 {
1195 int Fce_id = FceFwd ? fce_offset + cnt + i
1196 : fce_offset + cnt + nsplit - j - 2 - i;
1197 E_fidxy2(0, j, i) = Fce_id;
1198 }
1199 cnt += nsplit - j - 1;
1200 }
1201
1202 /* face 1 */
1203 fid = geom->GetFace(1)->GetGlobalID();
1204 fce_offset = 2 * fid * nsplit * nsplit;
1205 FceFwd = geom->GetForient(1) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
1206 cnt = 0;
1207 for (int j = 0; j < nsplit; ++j)
1208 {
1209 for (int i = 0; i < nsplit - j; ++i)
1210 {
1211 int Fce_id = FceFwd ? fce_offset + cnt + i
1212 : fce_offset + cnt + nsplit - j - 1 - i;
1213 E_fidxz1(j, 0, i) = Fce_id;
1214
1215 // due to edge shuffling earlier need to identify know
1216 // orientation (singlar vertex at top) orientation to
1217 // fill in stored edges and vertices
1218 int oset = FceEdgOffset[Fce_id];
1219
1220 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
1221 // fill in vertices & edges from faces
1222 if (FceFwd)
1223 {
1224 E_vid(j, 0, i) = face->GetVid(oset);
1225 E_vid(j, 0, i + 1) = face->GetVid((oset + 1) % 3);
1226
1227 E_eidz(j, 0, i) = face->GetEid((oset + 2) % 3);
1228 E_eidxz(j, 0, i) = face->GetEid((oset + 1) % 3);
1229 }
1230 else
1231 {
1232 E_vid(j, 0, i) = face->GetVid((oset + 1) % 3);
1233 E_vid(j, 0, i + 1) = face->GetVid(oset);
1234
1235 E_eidz(j, 0, i) = face->GetEid((oset + 1) % 3);
1236 E_eidxz(j, 0, i) = face->GetEid((oset + 2) % 3);
1237 }
1238 E_vid(j + 1, 0, i) = face->GetVid((oset + 2) % 3);
1239 E_eidx(j, 0, i) = face->GetEid(oset);
1240 }
1241 cnt += nsplit - j;
1242
1243 for (int i = 0; i < nsplit - j - 1; ++i)
1244 {
1245 int Fce_id = FceFwd ? fce_offset + cnt + i
1246 : fce_offset + cnt + nsplit - j - 2 - i;
1247 E_fidxz2(j, 0, i) = Fce_id;
1248 }
1249 cnt += nsplit - j - 1;
1250 }
1251
1252 /* face 2 */
1253 fid = geom->GetFace(2)->GetGlobalID();
1254 fce_offset = 2 * fid * nsplit * nsplit;
1255 FceFwd = geom->GetForient(2) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
1256
1257 cnt = 0;
1258 for (int j = 0; j < nsplit; ++j)
1259 {
1260 for (int i = 0; i < nsplit - j; ++i)
1261 {
1262 int Fce_id = FceFwd ? fce_offset + cnt + i
1263 : fce_offset + cnt + nsplit - j - 1 - i;
1264
1265 E_fidyzd1(j, i, nsplit - j - i - 1) = Fce_id;
1266
1267 // due to edge shuffling earlier need to identify know
1268 // orientation (singlar vertex at top) orientation to
1269 // fill in stored edges and vertices
1270 int oset = FceEdgOffset[Fce_id];
1271
1272 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
1273 // fill in vertices & edges from faces
1274 if (FceFwd)
1275 {
1276 E_vid(j, i, nsplit - j - i) = face->GetVid(oset);
1277 E_vid(j, i + 1, nsplit - j - i - 1) =
1278 face->GetVid((oset + 1) % 3);
1279
1280 E_eidxz(j, i, nsplit - j - i - 1) =
1281 face->GetEid((oset + 2) % 3);
1282 E_eidyz(j, i, nsplit - j - i - 1) =
1283 face->GetEid((oset + 1) % 3);
1284 }
1285 else
1286 {
1287 E_vid(j, i, nsplit - j - i) = face->GetVid((oset + 1) % 3);
1288 E_vid(j, i + 1, nsplit - j - i - 1) = face->GetVid(oset);
1289
1290 E_eidxz(j, i, nsplit - j - i - 1) =
1291 face->GetEid((oset + 1) % 3);
1292 E_eidyz(j, i, nsplit - j - i - 1) =
1293 face->GetEid((oset + 2) % 3);
1294 }
1295
1296 E_vid(j + 1, i, nsplit - j - i - 1) =
1297 face->GetVid((oset + 2) % 3);
1298 E_eidxy(j, i, nsplit - j - i - 1) = face->GetEid(oset);
1299
1300 // reorient face if required. Note cannto do earlier
1301 // since need to setup the arrays above
1302 // LinMesh->m_triGeoms[Fce_id] = ReOrientFace(face);
1303 }
1304 cnt += nsplit - j;
1305
1306 for (int i = 0; i < nsplit - j - 1; ++i)
1307 {
1308 int Fce_id = FceFwd ? fce_offset + cnt + i
1309 : fce_offset + cnt + nsplit - j - 2 - i;
1310 E_fidyzd2(j, i, nsplit - j - i - 2) = Fce_id;
1311 }
1312 cnt += nsplit - j - 1;
1313 }
1314
1315 /* face 3 */
1316 fid = geom->GetFace(3)->GetGlobalID();
1317 fce_offset = 2 * fid * nsplit * nsplit;
1318 FceFwd = geom->GetForient(3) == StdRegions::eDir1FwdDir1_Dir2FwdDir2;
1319 cnt = 0;
1320 for (int j = 0; j < nsplit; ++j)
1321 {
1322 for (int i = 0; i < nsplit - j; ++i)
1323 {
1324 int Fce_id = FceFwd ? fce_offset + cnt + i
1325 : fce_offset + cnt + nsplit - j - 1 - i;
1326 E_fidyz1(j, i, 0) = Fce_id;
1327
1328 // due to edge shuffling earlier need to identify know
1329 // orientation (singlar vertex at top) orientation to
1330 // fill in stored edges and vertices
1331 int oset = FceEdgOffset[Fce_id];
1332
1333 TriGeom *face = m_linMesh->GetTriGeom(Fce_id);
1334 // fill in vertices & edges from faces
1335 if (FceFwd)
1336 {
1337 E_vid(j, i, 0) = face->GetVid(oset);
1338 E_vid(j, i + 1, 0) = face->GetVid((oset + 1) % 3);
1339
1340 E_eidz(j, i, 0) = face->GetEid((oset + 2) % 3);
1341 E_eidyz(j, i, 0) = face->GetEid((oset + 1) % 3);
1342 }
1343 else
1344 {
1345 E_vid(j, i, 0) = face->GetVid((oset + 1) % 3);
1346 E_vid(j, i + 1, 0) = face->GetVid(oset);
1347
1348 E_eidz(j, i, 0) = face->GetEid((oset + 1) % 3);
1349 E_eidyz(j, i, 0) = face->GetEid((oset + 2) % 3);
1350 }
1351 E_vid(j + 1, i, 0) = face->GetVid((oset + 2) % 3);
1352 E_eidy(j, i, 0) = face->GetEid(oset);
1353 }
1354 cnt += nsplit - j;
1355
1356 for (int i = 0; i < nsplit - j - 1; ++i)
1357 {
1358 int Fce_id = FceFwd ? fce_offset + cnt + i
1359 : fce_offset + cnt + nsplit - j - 2 - i;
1360 E_fidyz2(j, i, 0) = Fce_id;
1361 }
1362 cnt += nsplit - j - 1;
1363 }
1364
1365 StdRegions::StdExpansionSharedPtr Xmap = geom->GetXmap();
1366
1367 Array<OneD, NekDouble> tmp(Xmap->GetTotPoints());
1368
1369 // fill in interior vids if exist
1370 if (UseGLL) // get GLL points
1371 {
1372 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
1373 Xmap->PhysInterpToGLL(tmp, Xpts, nsplit + 1);
1374 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
1375 Xmap->PhysInterpToGLL(tmp, Ypts, nsplit + 1);
1376 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
1377 Xmap->PhysInterpToGLL(tmp, Zpts, nsplit + 1);
1378 }
1379 else // get equispaced points
1380 {
1381 Xmap->BwdTrans(geom->GetCoeffs(0), tmp);
1382 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Xpts, nsplit + 1);
1383 Xmap->BwdTrans(geom->GetCoeffs(1), tmp);
1384 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Ypts, nsplit + 1);
1385 Xmap->BwdTrans(geom->GetCoeffs(2), tmp);
1386 Xmap->PhysInterpToSimplexEquiSpaced(tmp, Zpts, nsplit + 1);
1387 }
1388
1389 cnt = (nsplit + 2) * (nsplit + 1) / 2;
1390 for (int k = 1; k < nsplit; ++k)
1391 {
1392 cnt += nsplit + 1 - k; // first row
1393 for (int j = 1; j < nsplit - k; ++j)
1394 {
1395 cnt++;
1396 for (int i = 1; i < nsplit - k - j; ++i, cnt++)
1397 {
1398 int vid =
1399 maxvertid +
1400 elmtid * (nsplit - 1) * (nsplit - 1) * (nsplit - 1) +
1401 vcnt++;
1402
1403 E_vid(k, j, i) = vid;
1404
1405 PointGeomUniquePtr vert =
1406 ObjPoolManager<PointGeom>::AllocateUniquePtr(
1407 m_spaceDimension, vid, Xpts[cnt], Ypts[cnt],
1408 Zpts[cnt]);
1409
1410 m_linMesh->AddGeom(vid, std::move(vert));
1411 }
1412 cnt++;
1413 }
1414 cnt++;
1415 }
1416
1417 // declare Tets
1418 cnt = 0;
1419 for (int k = 0; k < nsplit; ++k)
1420 {
1421 for (int j = 0; j < nsplit - k; ++j)
1422 {
1423 for (int i = 0; i < nsplit - k - j; ++i)
1424 {
1425 if (i < nsplit - k - j - 2)
1426 {
1427 /* addd in 6 new edges */
1428 // edge 0
1429 int edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1430 std::array<PointGeom *, 2> vert;
1431 vert[0] =
1432 m_linMesh->GetPointGeom(E_vid(k, j + 1, i + 1));
1433 vert[1] =
1434 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
1435 SegGeomUniquePtr edge0 =
1436 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1437 edgid, m_spaceDimension, vert);
1438 m_linMesh->AddGeom(edgid, std::move(edge0));
1439 E_eidz(k, j + 1, i + 1) = edgid;
1440
1441 // edge 1
1442 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1443 vert[0] =
1444 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i));
1445 vert[1] =
1446 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
1447 SegGeomUniquePtr edge1 =
1448 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1449 edgid, m_spaceDimension, vert);
1450 m_linMesh->AddGeom(edgid, std::move(edge1));
1451 E_eidx(k + 1, j + 1, i) = edgid;
1452
1453 // edge 2
1454 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1455 vert[0] =
1456 m_linMesh->GetPointGeom(E_vid(k + 1, j, i + 1));
1457 vert[1] =
1458 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i + 1));
1459 SegGeomUniquePtr edge2 =
1460 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1461 edgid, m_spaceDimension, vert);
1462 m_linMesh->AddGeom(edgid, std::move(edge2));
1463 E_eidy(k + 1, j, i + 1) = edgid;
1464
1465 // edge 3
1466 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1467 vert[0] =
1468 m_linMesh->GetPointGeom(E_vid(k, j + 1, i + 1));
1469 vert[1] =
1470 m_linMesh->GetPointGeom(E_vid(k + 1, j, i + 1));
1471 SegGeomUniquePtr edge3 =
1472 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1473 edgid, m_spaceDimension, vert);
1474 m_linMesh->AddGeom(edgid, std::move(edge3));
1475 E_eidyz(k, j, i + 1) = edgid;
1476
1477 // edge 4
1478 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1479 vert[0] =
1480 m_linMesh->GetPointGeom(E_vid(k, j + 1, i + 1));
1481 vert[1] =
1482 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i));
1483 SegGeomUniquePtr edge4 =
1484 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1485 edgid, m_spaceDimension, vert);
1486 m_linMesh->AddGeom(edgid, std::move(edge4));
1487 E_eidxz(k, j + 1, i) = edgid;
1488
1489 // edge 5
1490 edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1491 vert[0] =
1492 m_linMesh->GetPointGeom(E_vid(k + 1, j, i + 1));
1493 vert[1] =
1494 m_linMesh->GetPointGeom(E_vid(k + 1, j + 1, i));
1495 SegGeomUniquePtr edge5 =
1496 ObjPoolManager<SegGeom>::AllocateUniquePtr(
1497 edgid, m_spaceDimension, vert);
1498 m_linMesh->AddGeom(edgid, std::move(edge5));
1499 E_eidxy(k + 1, j, i) = edgid;
1500
1501 // /* Add in 4 new faces */
1502
1503 /* Tri0 = bottom left Tri */
1504 int newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
1505 std::array<SegGeom *, 3> edges;
1506 edges[0] = m_linMesh->GetSegGeom(E_eidyz(k, j, i + 1));
1507 edges[1] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
1508 edges[2] = m_linMesh->GetSegGeom(E_eidxy(k + 1, j, i));
1509 // shuffle so lowest vid is on edges 1,2
1510 TriShuffleEdges(edges);
1511 TriGeomUniquePtr tri0 =
1512 ObjPoolManager<SpatialDomains::TriGeom>::
1513 AllocateUniquePtr(newtriid, edges);
1514
1515 m_linMesh->AddGeom(newtriid, std::move(tri0));
1516 E_fidyzd2(k, j, i) = newtriid;
1517
1518 /* Tri1 = */
1519 newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
1520 edges[0] = m_linMesh->GetSegGeom(E_eidyz(k, j, i + 1));
1521 edges[1] =
1522 m_linMesh->GetSegGeom(E_eidz(k, j + 1, i + 1));
1523 edges[2] =
1524 m_linMesh->GetSegGeom(E_eidy(k + 1, j, i + 1));
1525 // shuffle so lowest vid is on edges 1, 2
1526 TriShuffleEdges(edges);
1527 TriGeomUniquePtr tri1 =
1528 ObjPoolManager<SpatialDomains::TriGeom>::
1529 AllocateUniquePtr(newtriid, edges);
1530 m_linMesh->AddGeom(newtriid, std::move(tri1));
1531 E_fidyz2(k, j, i + 1) = newtriid;
1532
1533 /* Tri2 = */
1534 newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
1535 edges[0] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
1536 edges[1] =
1537 m_linMesh->GetSegGeom(E_eidz(k, j + 1, i + 1));
1538 edges[2] =
1539 m_linMesh->GetSegGeom(E_eidx(k + 1, j + 1, i));
1540 // shuffle so lowest vid is on edges 1,2
1541 TriShuffleEdges(edges);
1542 TriGeomUniquePtr tri2 =
1543 ObjPoolManager<SpatialDomains::TriGeom>::
1544 AllocateUniquePtr(newtriid, edges);
1545 m_linMesh->AddGeom(newtriid, std::move(tri2));
1546 E_fidxz2(k, j + 1, i) = newtriid;
1547
1548 /* Tri3 = */
1549 newtriid = elmtid * faceoffset + maxfaceid + fcnt++;
1550 edges[0] = m_linMesh->GetSegGeom(E_eidxy(k + 1, j, i));
1551 edges[1] =
1552 m_linMesh->GetSegGeom(E_eidy(k + 1, j, i + 1));
1553 edges[2] =
1554 m_linMesh->GetSegGeom(E_eidx(k + 1, j + 1, i));
1555 // shuffle so lowest vid is on
1556 TriShuffleEdges(edges);
1557 // edges 1, 2
1558 TriGeomUniquePtr tri3 =
1559 ObjPoolManager<SpatialDomains::TriGeom>::
1560 AllocateUniquePtr(newtriid, edges);
1561 m_linMesh->AddGeom(newtriid, std::move(tri3));
1562 E_fidxy2(k + 1, j, i) = newtriid;
1563
1564 /* add new top Corner Tet */
1565 int newtetid =
1566 elmtid * nsplit * nsplit * nsplit + cnt++;
1567 std::array<TriGeom *, 4> faces;
1568 faces[0] =
1569 m_linMesh->GetTriGeom(E_fidyzd2(k, j, i)); // tri0
1570 faces[1] = m_linMesh->GetTriGeom(
1571 E_fidyz2(k, j, i + 1)); // tri1
1572 faces[2] = m_linMesh->GetTriGeom(
1573 E_fidxz2(k, j + 1, i)); // tri2
1574 faces[3] = m_linMesh->GetTriGeom(
1575 E_fidxy2(k + 1, j, i)); // tri3
1576 TetShuffleFaces(faces); // shuffle so max
1577 // vert id at top
1578 TetGeomUniquePtr tet =
1579 ObjPoolManager<SpatialDomains::TetGeom>::
1580 AllocateUniquePtr(newtetid, faces);
1581 m_linMesh->AddGeom(newtetid, std::move(tet));
1582 m_linMesh->PopulateFaceToElMap(
1583 m_linMesh->GetTetGeom(newtetid), TetGeom::kNfaces);
1584 }
1585
1586 if (i < nsplit - k - j - 1)
1587 {
1588 /* addd in new edge (off
1589 diagonal) */
1590 int edgid = elmtid * edgeoffset + maxedgeid + ecnt++;
1591
1592 std::array<PointGeom *, 2> vert;
1593 vert[0] = m_linMesh->GetPointGeom(E_vid(k, j + 1, i));
1594 vert[1] =
1595 m_linMesh->GetPointGeom(E_vid(k + 1, j, i + 1));
1596
1597 SegGeomUniquePtr edge =
1598 ObjPoolManager<SpatialDomains::SegGeom>::
1599 AllocateUniquePtr(edgid, m_spaceDimension,
1600 vert);
1601 m_linMesh->AddGeom(edgid, std::move(edge));
1602
1603 /* Add in 8 new faces */
1604
1605 /* Tri0 = bottom left Tri */
1606 int newtriid0 =
1607 elmtid * faceoffset + maxfaceid + fcnt++;
1608
1609 std::array<SegGeom *, 3> edges;
1610 edges[0] = m_linMesh->GetSegGeom(E_eidxy(k, j, i));
1611 edges[1] = m_linMesh->GetSegGeom(E_eidyz(k, j, i));
1612 edges[2] = m_linMesh->GetSegGeom(E_eidxz(k, j, i));
1613 // shuffle so lowest vid is on edges 1,2
1614 TriShuffleEdges(edges);
1615 TriGeomUniquePtr tri0 =
1616 ObjPoolManager<SpatialDomains::TriGeom>::
1617 AllocateUniquePtr(newtriid0, edges);
1618 m_linMesh->AddGeom(newtriid0, std::move(tri0));
1619 E_fidyzd1(k, j, i) = newtriid0;
1620
1621 /* Tri1 temp Tri - yz diagonal
1622 */
1623 int newtriid1 =
1624 elmtid * faceoffset + maxfaceid + fcnt++;
1625 edges[0] = m_linMesh->GetSegGeom(E_eidyz(k, j, i));
1626 edges[1] = m_linMesh->GetSegGeom(edgid);
1627 edges[2] = m_linMesh->GetSegGeom(E_eidx(k + 1, j, i));
1628 // shuffle so lowest vid is on edges 1,2
1629 TriShuffleEdges(edges);
1630 TriGeomUniquePtr tri1 =
1631 ObjPoolManager<SpatialDomains::TriGeom>::
1632 AllocateUniquePtr(newtriid1, edges);
1633 m_linMesh->AddGeom(newtriid1, std::move(tri1));
1634
1635 /* Tri2 temp Tri yz diagonal */
1636 int newtriid2 =
1637 elmtid * faceoffset + maxfaceid + fcnt++;
1638 edges[0] = m_linMesh->GetSegGeom(E_eidx(k, j + 1, i));
1639 edges[1] = m_linMesh->GetSegGeom(edgid);
1640 edges[2] = m_linMesh->GetSegGeom(E_eidyz(k, j, i + 1));
1641 // shuffle so lowest vid is on edges 1, 2
1642 TriShuffleEdges(edges);
1643 TriGeomUniquePtr tri2 =
1644 ObjPoolManager<SpatialDomains::TriGeom>::
1645 AllocateUniquePtr(newtriid2, edges);
1646 m_linMesh->AddGeom(newtriid2, std::move(tri2));
1647
1648 /* Tri3 temp Tri xz diagonal */
1649 int newtriid3 =
1650 elmtid * faceoffset + maxfaceid + fcnt++;
1651 edges[0] = m_linMesh->GetSegGeom(E_eidxy(k, j, i));
1652 edges[1] = m_linMesh->GetSegGeom(edgid);
1653 edges[2] = m_linMesh->GetSegGeom(E_eidz(k, j, i + 1));
1654 // shuffle so lowest vid is on edges 1, 2
1655 TriShuffleEdges(edges);
1656 TriGeomUniquePtr tri3 =
1657 ObjPoolManager<SpatialDomains::TriGeom>::
1658 AllocateUniquePtr(newtriid3, edges);
1659 m_linMesh->AddGeom(newtriid3, std::move(tri3));
1660
1661 /* Tri4 temp Tri xz diagonal */
1662 int newtriid4 =
1663 elmtid * faceoffset + maxfaceid + fcnt++;
1664 edges[0] = m_linMesh->GetSegGeom(edgid);
1665 edges[1] = m_linMesh->GetSegGeom(E_eidxy(k + 1, j, i));
1666 edges[2] = m_linMesh->GetSegGeom(E_eidz(k, j + 1, i));
1667 // shuffle so lowest vid is on edges 1,2
1668 TriShuffleEdges(edges);
1669 TriGeomUniquePtr tri4 =
1670 ObjPoolManager<SpatialDomains::TriGeom>::
1671 AllocateUniquePtr(newtriid4, edges);
1672 m_linMesh->AddGeom(newtriid4, std::move(tri4));
1673
1674 // /* Tri5 Tri xz vertical */
1675 int newtriid5 =
1676 elmtid * faceoffset + maxfaceid + fcnt++;
1677 edges[0] = m_linMesh->GetSegGeom(E_eidx(k, j + 1, i));
1678 edges[1] = m_linMesh->GetSegGeom(E_eidxz(k, j + 1, i));
1679 edges[2] = m_linMesh->GetSegGeom(E_eidz(k, j + 1, i));
1680 // shuffle so lowest vid is on edges 1,2
1681 TriShuffleEdges(edges);
1682 TriGeomUniquePtr tri5 =
1683 ObjPoolManager<SpatialDomains::TriGeom>::
1684 AllocateUniquePtr(newtriid5, edges);
1685 m_linMesh->AddGeom(newtriid5, std::move(tri5));
1686 E_fidxz1(k, j + 1, i) = newtriid5;
1687
1688 /* Tri6 Tri yz vertical */
1689 int newtriid6 =
1690 elmtid * faceoffset + maxfaceid + fcnt++;
1691 edges[0] = m_linMesh->GetSegGeom(E_eidy(k, j, i + 1));
1692 edges[1] = m_linMesh->GetSegGeom(E_eidyz(k, j, i + 1));
1693 edges[2] = m_linMesh->GetSegGeom(E_eidz(k, j, i + 1));
1694 // shuffle so lowest vid is on edges 1,2
1695 TriShuffleEdges(edges);
1696 TriGeomUniquePtr tri6 =
1697 ObjPoolManager<SpatialDomains::TriGeom>::
1698 AllocateUniquePtr(newtriid6, edges);
1699 m_linMesh->AddGeom(newtriid6, std::move(tri6));
1700 E_fidyz1(k, j, i + 1) = newtriid6;
1701
1702 // /* Tri7 Tri xy horizontal */
1703 int newtriid7 =
1704 elmtid * faceoffset + maxfaceid + fcnt++;
1705 edges[0] = m_linMesh->GetSegGeom(E_eidx(k + 1, j, i));
1706 edges[1] = m_linMesh->GetSegGeom(E_eidy(k + 1, j, i));
1707 edges[2] = m_linMesh->GetSegGeom(E_eidxy(k + 1, j, i));
1708 // shuffle so lowest vid is on edges 1,2
1709 TriShuffleEdges(edges);
1710 TriGeomUniquePtr tri7 =
1711 ObjPoolManager<SpatialDomains::TriGeom>::
1712 AllocateUniquePtr(newtriid7, edges);
1713 m_linMesh->AddGeom(newtriid7, std::move(tri7));
1714 E_fidxy1(k + 1, j, i) = newtriid7;
1715
1716 /* add in four new Tets */
1717 // Tet 0
1718 int newtetid =
1719 elmtid * nsplit * nsplit * nsplit + cnt++;
1720 std::array<TriGeom *, 4> faces;
1721 faces[0] = m_linMesh->GetTriGeom(E_fidyzd1(k, j, i));
1722 faces[1] = m_linMesh->GetTriGeom(newtriid3);
1723 faces[2] = m_linMesh->GetTriGeom(newtriid1);
1724 faces[3] = m_linMesh->GetTriGeom(E_fidxz2(k, j, i));
1725 // store location of vertid
1726 TetShuffleFaces(faces); // shuffle so maxvert id at top
1727 TetGeomUniquePtr tet0 =
1728 ObjPoolManager<SpatialDomains::TetGeom>::
1729 AllocateUniquePtr(newtetid, faces);
1730 m_linMesh->AddGeom(newtetid, std::move(tet0));
1731 m_linMesh->PopulateFaceToElMap(
1732 m_linMesh->GetTetGeom(newtetid), 4);
1733
1734 // Tet 1
1735 newtetid = elmtid * nsplit * nsplit * nsplit + cnt++;
1736 faces[0] = m_linMesh->GetTriGeom(newtriid1);
1737 faces[1] = m_linMesh->GetTriGeom(E_fidyz2(k, j, i));
1738 faces[2] = m_linMesh->GetTriGeom(E_fidxy1(k + 1, j, i));
1739 faces[3] = m_linMesh->GetTriGeom(newtriid4);
1740 TetShuffleFaces(faces);
1741 // shuffle so max vert id at top
1742 TetGeomUniquePtr tet1 =
1743 ObjPoolManager<SpatialDomains::TetGeom>::
1744 AllocateUniquePtr(newtetid, faces);
1745 m_linMesh->AddGeom(newtetid, std::move(tet1));
1746 m_linMesh->PopulateFaceToElMap(
1747 m_linMesh->GetTetGeom(newtetid), 4);
1748
1749 // Tet 2
1750 newtetid = elmtid * nsplit * nsplit * nsplit + cnt++;
1751 faces[0] = m_linMesh->GetTriGeom(E_fidxy2(k, j, i));
1752 faces[1] = m_linMesh->GetTriGeom(newtriid3);
1753 faces[2] = m_linMesh->GetTriGeom(E_fidyz1(k, j, i + 1));
1754 faces[3] = m_linMesh->GetTriGeom(newtriid2);
1755 TetShuffleFaces(faces);
1756 // shuffle so max vert id at top
1757 TetGeomUniquePtr tet2 =
1758 ObjPoolManager<SpatialDomains::TetGeom>::
1759 AllocateUniquePtr(newtetid, faces);
1760 m_linMesh->AddGeom(newtetid, std::move(tet2));
1761 m_linMesh->PopulateFaceToElMap(
1762 m_linMesh->GetTetGeom(newtetid), 4);
1763
1764 // Tet 3
1765 newtetid = elmtid * nsplit * nsplit * nsplit + cnt++;
1766 faces[0] = m_linMesh->GetTriGeom(newtriid2);
1767 faces[1] = m_linMesh->GetTriGeom(E_fidxz1(k, j + 1, i));
1768 faces[2] = m_linMesh->GetTriGeom(newtriid4);
1769 faces[3] = m_linMesh->GetTriGeom(E_fidyzd2(k, j, i));
1770 TetShuffleFaces(faces); // shuffle so max
1771 // vert id at top
1772 TetGeomUniquePtr tet3 =
1773 ObjPoolManager<SpatialDomains::TetGeom>::
1774 AllocateUniquePtr(newtetid, faces);
1775 m_linMesh->AddGeom(newtetid, std::move(tet3));
1776 m_linMesh->PopulateFaceToElMap(
1777 m_linMesh->GetTetGeom(newtetid), 4);
1778 }
1779
1780 /* bottom Corner Tet */
1781 int newtetid = elmtid * nsplit * nsplit * nsplit + cnt++;
1782
1783 std::array<TriGeom *, 4> faces;
1784 faces[0] = m_linMesh->GetTriGeom(E_fidxy1(k, j, i));
1785 faces[1] = m_linMesh->GetTriGeom(E_fidxz1(k, j, i));
1786 faces[2] = m_linMesh->GetTriGeom(E_fidyzd1(k, j, i));
1787 faces[3] = m_linMesh->GetTriGeom(E_fidyz1(k, j, i));
1788 TetShuffleFaces(faces); // shuffle
1789 // so max vert id at top
1790 TetGeomUniquePtr tet = ObjPoolManager<
1791 SpatialDomains::TetGeom>::AllocateUniquePtr(newtetid,
1792 faces);
1793 m_linMesh->AddGeom(newtetid, std::move(tet));
1794 m_linMesh->PopulateFaceToElMap(
1795 m_linMesh->GetTetGeom(newtetid), 4);
1796 }
1797 }
1798 }
1799 ASSERTL1(ecnt <= edgeoffset, "This should not happen, "
1800 "edgeoffset might need to be increased");
1801 ASSERTL1(fcnt <= faceoffset, "This should not happen, "
1802 "faceoffset might need to be increased");
1803
1804 // collect a mapping of Vid to point vertex location
1805 cnt = 0;
1806 for (int k = 0; k < nsplit + 1; ++k)
1807 {
1808 for (int j = 0; j < nsplit + 1 - k; ++j)
1809 {
1810 for (int i = 0; i < nsplit + 1 - j - k; ++i)
1811 {
1812 CoeffMap[elmtid][E_vid(k, j, i)] = cnt++;
1813 }
1814 }
1815 }
1816 }
1817}
Nektar::SpatialDomains::LinearMeshGraph::TetShuffleFaces
void TetShuffleFaces(std::array< TriGeom *, 4 > &faces)
Definition LinearMeshGraph.hpp:242
Nektar::SpatialDomains::TetGeom::kNfaces
static const int kNfaces
Definition TetGeom.h:55
Nektar::SpatialDomains::TetGeomUniquePtr
unique_ptr_objpool< TetGeom > TetGeomUniquePtr
Definition MeshGraph.h:105

References Nektar::ObjPoolManager< DataType >::AllocateUniquePtr(), ASSERTL1, Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, Nektar::SpatialDomains::Geometry::GetEid(), Nektar::SpatialDomains::Geometry::GetVid(), Nektar::SpatialDomains::TetGeom::kNfaces, m_graph, m_linMesh, m_session, m_spaceDimension, Nektar::LibUtilities::ReduceMax, TetShuffleFaces(), and TriShuffleEdges().

Referenced by CreateLinearGraph().

◆ PrismShuffleFaces()

void Nektar::SpatialDomains::LinearMeshGraph::PrismShuffleFaces ( std::array< Geometry2D *, 5 > &  faces)
inlineprivate

Definition at line 294 of file LinearMeshGraph.hpp.

295 {
296 // assume second face defined singular vertex
297 int edg0 = faces[1]->GetEdge(0)->GetGlobalID();
298
299 // find face with this same edge
300 int i, j;
301 for (i = 0; i < 5; i += 2)
302 {
303 for (j = 0; j < 4; ++j)
304 {
305 if (faces[i]->GetEdge(j)->GetGlobalID() == edg0)
306 {
307 break;
308 }
309 }
310 if (j != 4)
311 {
312 break;
313 }
314 }
315 ASSERTL1(i != 6, "Failed to find common edge");
316
317 /* rotate faces */
318
319 if (i == 2)
320 {
321 Geometry2D *save0 = std::move(faces[0]);
322 faces[0] = std::move(faces[2]);
323 faces[2] = std::move(faces[4]);
324 faces[4] = std::move(save0);
325 }
326 else if (i == 4)
327 {
328 Geometry2D *save0 = std::move(faces[0]);
329 faces[0] = std::move(faces[4]);
330 faces[4] = std::move(faces[2]);
331 faces[2] = std::move(save0);
332 }
333 }

References ASSERTL1.

Referenced by LinMeshSetUpPrismGeom().

◆ SetupLinearExpansionType()

TiXmlElement * Nektar::SpatialDomains::LinearMeshGraph::SetupLinearExpansionType ( void  )
inlineprivate

Definition at line 143 of file LinearMeshGraph.hpp.

144 {
145 TiXmlElement *expansionTypes =
146 m_session->GetElement("NEKTAR/EXPANSIONS");
147
148 TiXmlElement *newExpType = new TiXmlElement("NEKTAR/EXPANSIONS");
149
150 // Find the Expansions tag
151 ASSERTL0(expansionTypes, "Unable to find EXPANSIONS tag in file.");
152
153 if (expansionTypes)
154 {
155 TiXmlElement *expansion = expansionTypes->FirstChildElement();
156 ASSERTL0(expansion, "Unable to find entries in EXPANSIONS tag in"
157 "file.");
158 std::string expType = expansion->Value();
159 std::vector<std::string> vars = m_session->GetVariables();
160
161 if (expType == "E")
162 {
163 while (expansion)
164 {
165 TiXmlElement *exp = new TiXmlElement("E");
166
167 std::string expstr = expansion->Attribute("COMPOSITE");
168 exp->SetAttribute("COMPOSITE", expstr);
169 exp->SetAttribute("NUMMODES", 2);
170 const char *rType = expansion->Attribute("TYPE");
171 if (rType) // use same type if exists
172 {
173 std::string typestr = expansion->Attribute("TYPE");
174 exp->SetAttribute("TYPE", typestr);
175 }
176 else // otherwise specify as modified. Not currently
177 // using BasisType
178 {
179 exp->SetAttribute("TYPE", "MODIFIED");
180 }
181 std::string fieldstr = expansion->Attribute("FIELDS");
182 exp->SetAttribute("FIELDS", fieldstr);
183
184 newExpType->LinkEndChild(exp);
185 expansion = expansion->NextSiblingElement("E");
186 }
187 }
188 else
189 {
190 ASSERTL0(false, "Expansion type not defined");
191 }
192 }
193
194 return newExpType;
195 }

References ASSERTL0, and m_session.

Referenced by LinMeshSetUpCompositesDomain().

◆ TetShuffleFaces()

void Nektar::SpatialDomains::LinearMeshGraph::TetShuffleFaces ( std::array< TriGeom *, 4 > &  faces)
inlineprivate

reshuffles tet to point to singular vertices assuming faces are ordered in conventiaional manner

We know all faces are orientated so singular vertex on vertex 2

Definition at line 242 of file LinearMeshGraph.hpp.

243 {
244 std::array<int, 4> faceid = {-1, -1, -1, -1};
245 std::array<TriGeom *, 4> save;
246
247 for (int i = 0; i < 4; ++i)
248 {
249 ASSERTL0(faces[i] != nullptr,
250 "Null face pointer at index " + std::to_string(i));
251 }
252
253 /** We know all faces are orientated so singular vertex on vertex 2
254 **/
255 faceid[0] = 0;
256 /* face 0 */
257 for (int i = 1; i < 4; ++i)
258 {
259 if (faces[i]->GetVid(2) < faces[faceid[0]]->GetVid(2))
260 {
261 faceid[0] = i;
262 }
263 }
264
265 // /* face 1 */
266 int eid = faces[faceid[0]]->GetEid(0);
267 int i;
268 for (i = 0; i < 4; ++i)
269 {
270 if ((i != faceid[0]) && (faces[i]->GetEid(0) == eid))
271 {
272 faceid[1] = i;
273 break;
274 }
275 }
276 ASSERTL0(i != 4, "Failed to find face 1 id in shuffle_faces: faces not "
277 "set up correctly");
278
279 // // assign last two faces according to predefined map above.
280 faceid[2] = TetFaceOrder[faceid[0]][faceid[1]][0];
281 faceid[3] = TetFaceOrder[faceid[0]][faceid[1]][1];
282
283 for (int i = 0; i < 4; ++i)
284 {
285 save[i] = std::move(faces[i]);
286 }
287
288 for (int i = 0; i < 4; ++i)
289 {
290 faces[i] = std::move(save[faceid[i]]);
291 }
292 }
Nektar::SpatialDomains::LinearMeshGraph::TetFaceOrder
const unsigned int TetFaceOrder[4][4][2]
Definition LinearMeshGraph.hpp:234

References ASSERTL0, and TetFaceOrder.

Referenced by LinMeshSetUpTetGeom().

◆ TriShuffleEdges()

int Nektar::SpatialDomains::LinearMeshGraph::TriShuffleEdges ( std::array< SegGeom *, 3 > &  edges,
std::array< SegGeom *, 3 > *  edgesSort = nullptr 
)
inlineprivate

Definition at line 199 of file LinearMeshGraph.hpp.

201 {
202 std::array<SegGeom *, 3> save;
203 int maxVid0, maxVid1, startEd;
204
205 if (edgesSort == nullptr)
206 {
207 maxVid0 = std::max(edges[0]->GetVid(0), edges[0]->GetVid(1));
208 maxVid1 = std::max(edges[1]->GetVid(0), edges[1]->GetVid(1));
209
210 startEd = (maxVid0 == maxVid1) ? 2 : (maxVid0 > maxVid1) ? 1 : 0;
211 }
212 else
213 {
214 maxVid0 = std::max((*edgesSort)[0]->GetVid(0),
215 (*edgesSort)[0]->GetVid(1));
216 maxVid1 = std::max((*edgesSort)[1]->GetVid(0),
217 (*edgesSort)[1]->GetVid(1));
218 startEd = (maxVid0 == maxVid1) ? 2 : (maxVid0 > maxVid1) ? 1 : 0;
219 }
220
221 // // swap start and 0
222 save[0] = std::move(edges[0]);
223 save[1] = std::move(edges[1]);
224 save[2] = std::move(edges[2]);
225
226 edges[0] = std::move(save[startEd]);
227 edges[1] = std::move(save[(startEd + 1) % 3]);
228 edges[2] = std::move(save[(startEd + 2) % 3]);
229
230 return (3 - startEd) % 3;
231 }

Referenced by AddSplitTri(), LinMeshSetUpPrismGeom(), and LinMeshSetUpTetGeom().

Member Data Documentation

◆ m_graph

SpatialDomains::MeshGraphSharedPtr Nektar::SpatialDomains::LinearMeshGraph::m_graph
private

Definition at line 336 of file LinearMeshGraph.hpp.

Referenced by CreateLinearGraph(), LinMeshSetUp1DGeom(), LinMeshSetUp2DGeom(), LinMeshSetUpCompositesDomain(), LinMeshSetUpPrismGeom(), and LinMeshSetUpTetGeom().

◆ m_linMesh

SpatialDomains::MeshGraphSharedPtr Nektar::SpatialDomains::LinearMeshGraph::m_linMesh
private

Definition at line 337 of file LinearMeshGraph.hpp.

Referenced by AddSplitEdge(), AddSplitTri(), CreateLinearGraph(), GetLinearGraph(), LinMeshSetUp1DGeom(), LinMeshSetUp2DGeom(), LinMeshSetUpCompositesDomain(), LinMeshSetUpPrismGeom(), and LinMeshSetUpTetGeom().

◆ m_meshDimension

int Nektar::SpatialDomains::LinearMeshGraph::m_meshDimension
private

Definition at line 339 of file LinearMeshGraph.hpp.

Referenced by CreateLinearGraph(), LinMeshSetUp1DGeom(), and LinMeshSetUp2DGeom().

◆ m_session

LibUtilities::SessionReaderSharedPtr Nektar::SpatialDomains::LinearMeshGraph::m_session
private

Definition at line 338 of file LinearMeshGraph.hpp.

Referenced by CreateLinearGraph(), LinMeshSetUp1DGeom(), LinMeshSetUp2DGeom(), LinMeshSetUpPrismGeom(), LinMeshSetUpTetGeom(), and SetupLinearExpansionType().

◆ m_spaceDimension

int Nektar::SpatialDomains::LinearMeshGraph::m_spaceDimension
private

Definition at line 340 of file LinearMeshGraph.hpp.

Referenced by AddSplitEdge(), AddSplitTri(), CreateLinearGraph(), LinMeshSetUp1DGeom(), LinMeshSetUp2DGeom(), LinMeshSetUpPrismGeom(), and LinMeshSetUpTetGeom().

◆ TetFaceOrder

const unsigned int Nektar::SpatialDomains::LinearMeshGraph::TetFaceOrder[4][4][2]
private
Initial value:
= {
{{0, 0}, {2, 3}, {3, 1}, {1, 2}},
{{3, 2}, {0, 0}, {0, 3}, {2, 0}},
{{1, 3}, {3, 0}, {0, 0}, {0, 1}},
{{2, 1}, {0, 2}, {1, 0}, {0, 0}}}

Definition at line 234 of file LinearMeshGraph.hpp.

234 {
235 {{0, 0}, {2, 3}, {3, 1}, {1, 2}},
236 {{3, 2}, {0, 0}, {0, 3}, {2, 0}},
237 {{1, 3}, {3, 0}, {0, 0}, {0, 1}},
238 {{2, 1}, {0, 2}, {1, 0}, {0, 0}}};

Referenced by TetShuffleFaces().

Generated on Sat Jul 4 2026 01:02:31 for Nektar++ by doxygen 1.9.8