Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | List of all members
Nektar::FieldUtils::ProcessEquiSpacedOutput Class Reference

This processing module interpolates one field to another. More...

#include <ProcessEquiSpacedOutput.h>

Inheritance diagram for Nektar::FieldUtils::ProcessEquiSpacedOutput:
[legend]

Public Member Functions

 ProcessEquiSpacedOutput (FieldSharedPtr f)
 
virtual ~ProcessEquiSpacedOutput ()
 
- Public Member Functions inherited from Nektar::FieldUtils::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
- Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT Module (FieldSharedPtr p_f)
 
virtual ~Module ()=default
 
void Process (po::variables_map &vm)
 
std::string GetModuleName ()
 
std::string GetModuleDescription ()
 
const ConfigOption & GetConfigOption (const std::string &key) const
 
ModulePriority GetModulePriority ()
 
FIELD_UTILS_EXPORT void RegisterConfig (std::string key, std::string value="")
 Register a configuration option with a module. More...
 
FIELD_UTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
FIELD_UTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
FIELD_UTILS_EXPORT void AddFile (std::string fileType, std::string fileName)
 
FIELD_UTILS_EXPORT void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 

Static Public Member Functions

static std::shared_ptr< Module > create (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Protected Member Functions

virtual void v_Process (po::variables_map &vm) override
 Write mesh to output file. More...
 
virtual std::string v_GetModuleName () override
 
virtual std::string v_GetModuleDescription () override
 
virtual ModulePriority v_GetModulePriority () override
 
void SetHomogeneousConnectivity (void)
 
void GenOrthoModes (int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
 
- Protected Member Functions inherited from Nektar::FieldUtils::Module
 Module ()
 

Additional Inherited Members

- Public Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr m_f
 Field object. More...
 
- Protected Attributes inherited from Nektar::FieldUtils::Module
std::map< std::string, ConfigOption > m_config
 List of configuration values. More...
 
std::set< std::string > m_allowedFiles
 List of allowed file formats. More...
 

Detailed Description

This processing module interpolates one field to another.

Definition at line 48 of file ProcessEquiSpacedOutput.h.

Constructor & Destructor Documentation

◆ ProcessEquiSpacedOutput()

Nektar::FieldUtils::ProcessEquiSpacedOutput::ProcessEquiSpacedOutput ( FieldSharedPtr  f)

Definition at line 60 of file ProcessEquiSpacedOutput.cpp.

61  : ProcessModule(f)
62 {
63  m_config["tetonly"] =
64  ConfigOption(true, "0", "Only process tetrahedral elements");
65 
66  m_config["modalenergy"] =
67  ConfigOption(true, "0", "Write output as modal energy");
68 }
Nektar::FieldUtils::Module::m_config
std::map< std::string, ConfigOption > m_config
List of configuration values.
Definition: Module.h:263

References Nektar::FieldUtils::Module::m_config.

◆ ~ProcessEquiSpacedOutput()

Nektar::FieldUtils::ProcessEquiSpacedOutput::~ProcessEquiSpacedOutput ( )
virtual

Definition at line 70 of file ProcessEquiSpacedOutput.cpp.

71 {
72 }

Member Function Documentation

◆ create()

static std::shared_ptr<Module> Nektar::FieldUtils::ProcessEquiSpacedOutput::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 52 of file ProcessEquiSpacedOutput.h.

53  {
54  return MemoryManager<ProcessEquiSpacedOutput>::AllocateSharedPtr(f);
55  }
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

◆ GenOrthoModes()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::GenOrthoModes ( int  n,
const Array< OneD, const NekDouble > &  phys,
Array< OneD, NekDouble > &  coeffs 
)
protected

Definition at line 814 of file ProcessEquiSpacedOutput.cpp.

817 {
818  LocalRegions::ExpansionSharedPtr e;
819  e = m_f->m_exp[0]->GetExp(n);
820 
821  switch (e->DetShapeType())
822  {
823  case LibUtilities::eTriangle:
824  {
825  int np0 = e->GetBasis(0)->GetNumPoints();
826  int np1 = e->GetBasis(1)->GetNumPoints();
827  int np = max(np0, np1);
828 
829  // to ensure points are correctly projected to np need
830  // to increase the order slightly of coordinates
831  LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
832  LibUtilities::PointsKey pb(np, e->GetPointsType(1));
833  Array<OneD, NekDouble> tophys(np * (np + 1));
834 
835  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
836  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, np, pb);
837  StdRegions::StdTriExp OrthoExp(Ba, Bb);
838 
839  // interpolate points to new phys points!
840  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
841  e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
842  tophys);
843 
844  OrthoExp.FwdTrans(tophys, coeffs);
845  break;
846  }
847  case LibUtilities::eQuadrilateral:
848  {
849  int np0 = e->GetBasis(0)->GetNumPoints();
850  int np1 = e->GetBasis(1)->GetNumPoints();
851  int np = max(np0, np1);
852 
853  LibUtilities::PointsKey pa(np + 1, e->GetPointsType(0));
854  LibUtilities::PointsKey pb(np + 1, e->GetPointsType(1));
855  Array<OneD, NekDouble> tophys((np + 1) * (np + 1));
856 
857  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, np, pa);
858  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, np, pb);
859  StdRegions::StdQuadExp OrthoExp(Ba, Bb);
860 
861  // interpolate points to new phys points!
862  LibUtilities::Interp2D(e->GetBasis(0)->GetBasisKey(),
863  e->GetBasis(1)->GetBasisKey(), phys, Ba, Bb,
864  tophys);
865 
866  OrthoExp.FwdTrans(phys, coeffs);
867  break;
868  }
869  default:
870  ASSERTL0(false, "Shape needs setting up");
871  break;
872  }
873 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::FieldUtils::Module::m_f
FieldSharedPtr m_f
Field object.
Definition: Module.h:234
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:106
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
Nektar::LibUtilities::eOrtho_B
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68

References ASSERTL0, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::StdRegions::StdExpansion::FwdTrans(), Nektar::LibUtilities::Interp2D(), and Nektar::FieldUtils::Module::m_f.

Referenced by v_Process().

◆ SetHomogeneousConnectivity()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::SetHomogeneousConnectivity ( void  )
protected

Definition at line 423 of file ProcessEquiSpacedOutput.cpp.

424 {
425  LocalRegions::ExpansionSharedPtr e;
426  int nel = m_f->m_exp[0]->GetPlane(0)->GetExpSize();
427  int nPlanes = m_f->m_exp[0]->GetZIDs().size();
428  int npts = m_f->m_fieldPts->GetNpoints();
429  int nptsPerPlane = npts / nPlanes;
430  int coordim = 3;
431 
432  if (m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
433  LibUtilities::eFourier &&
434  m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
435  LibUtilities::eFourierEvenlySpaced)
436  {
437  // Write points with extra plane
438  Array<OneD, Array<OneD, NekDouble>> pts;
439  m_f->m_fieldPts->GetPts(pts);
440  Array<OneD, Array<OneD, NekDouble>> newPts(pts.size());
441  for (int i = 0; i < pts.size(); i++)
442  {
443  newPts[i] = Array<OneD, NekDouble>(npts + nptsPerPlane);
444  // Copy old points
445  Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
446  // Get new plane
447  Array<OneD, NekDouble> extraPlane(nptsPerPlane, newPts[i] + npts);
448  if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
449  {
450  if (i == (coordim - 1))
451  {
452  // z-coordinate
453  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
454  pts[i], 1, extraPlane, 1);
455  }
456  else
457  {
458  Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
459  }
460  }
461  else
462  {
463  // Determine to and from rank for communication
464  int size = m_f->m_comm->GetColumnComm()->GetSize();
465  int rank = m_f->m_comm->GetColumnComm()->GetRank();
466  int fromRank = (rank + 1) % size;
467  int toRank = (rank == 0) ? size - 1 : rank - 1;
468  // Communicate using SendRecv
469  Array<OneD, NekDouble> send(nptsPerPlane, pts[i]);
470  m_f->m_comm->GetColumnComm()->SendRecv(toRank, send, fromRank,
471  extraPlane);
472  // Adjust z-coordinate of last process
473  if (i == (coordim - 1) && (rank == size - 1))
474  {
475  Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
476  extraPlane, 1, extraPlane, 1);
477  }
478  }
479  }
480  m_f->m_fieldPts->SetPts(newPts);
481  // Now extend plane connectivity
482  vector<Array<OneD, int>> oldConn;
483  Array<OneD, int> conn;
484  m_f->m_fieldPts->GetConnectivity(oldConn);
485  vector<Array<OneD, int>> newConn = oldConn;
486  int connPerPlane = oldConn.size() / nPlanes;
487  for (int i = 0; i < connPerPlane; i++)
488  {
489  conn = Array<OneD, int>(oldConn[i].size());
490  for (int j = 0; j < conn.size(); j++)
491  {
492  conn[j] = oldConn[i][j] + npts;
493  }
494  newConn.push_back(conn);
495  }
496  m_f->m_fieldPts->SetConnectivity(newConn);
497 
498  nPlanes++;
499  npts += nptsPerPlane;
500  }
501 
502  vector<Array<OneD, int>> oldConn;
503  vector<Array<OneD, int>> newConn;
504  Array<OneD, int> conn;
505  m_f->m_fieldPts->GetConnectivity(oldConn);
506 
507  // 2DH1D case
508  if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTriBlock)
509  {
510  for (int i = 0; i < nel; ++i)
511  {
512  int nLines = oldConn[i].size() / 2;
513  // Create array for new connectivity
514  // (2 triangles between each plane for each line)
515  conn = Array<OneD, int>(2 * 3 * nLines * (nPlanes - 1));
516  int cnt = 0;
517  for (int n = 0; n < nLines; n++)
518  {
519  // Define new connectivity with triangles
520  for (int p = 0; p < nPlanes - 1; p++)
521  {
522  conn[cnt++] = oldConn[i + p * nel][2 * n + 0];
523  conn[cnt++] = oldConn[i + p * nel][2 * n + 1];
524  conn[cnt++] = oldConn[i + (p + 1) * nel][2 * n + 0];
525 
526  conn[cnt++] = oldConn[i + p * nel][2 * n + 1];
527  conn[cnt++] = oldConn[i + (p + 1) * nel][2 * n + 0];
528  conn[cnt++] = oldConn[i + (p + 1) * nel][2 * n + 1];
529  }
530  }
531  newConn.push_back(conn);
532  }
533  }
534  else if (m_f->m_fieldPts->GetPtsType() == LibUtilities::ePtsTetBlock)
535  {
536  // Get maximum number of points per direction in the mesh
537  int maxN = 0;
538  for (int i = 0; i < nel; ++i)
539  {
540  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
541 
542  int np0 = e->GetBasis(0)->GetNumPoints();
543  int np1 = e->GetBasis(1)->GetNumPoints();
544  int np = max(np0, np1);
545  maxN = max(np, maxN);
546  }
547 
548  // Create a global numbering for points in plane 0
549  Array<OneD, int> vId(nptsPerPlane);
550  int cnt1 = 0;
551  int cnt2 = 0;
552  for (int i = 0; i < nel; ++i)
553  {
554  e = m_f->m_exp[0]->GetPlane(0)->GetExp(i);
555  int np0 = e->GetBasis(0)->GetNumPoints();
556  int np1 = e->GetBasis(1)->GetNumPoints();
557  int np = max(np0, np1);
558  switch (e->DetShapeType())
559  {
560  case LibUtilities::eTriangle:
561  {
562  int newpoints =
563  LibUtilities::StdTriData::getNumberOfCoefficients(np,
564  np);
565 
566  // Vertex numbering
567  vId[cnt1] = e->GetGeom()->GetVid(0);
568  vId[cnt1 + np - 1] = e->GetGeom()->GetVid(1);
569  vId[cnt1 + newpoints - 1] = e->GetGeom()->GetVid(2);
570 
571  // Edge numbering
572  StdRegions::Orientation edgeOrient;
573  int edge1 = 0;
574  int edge2 = 0;
575  for (int n = 1; n < np - 1; n++)
576  {
577  // Edge 0
578  edgeOrient = e->GetGeom()->GetEorient(0);
579  if (edgeOrient == StdRegions::eForwards)
580  {
581  vId[cnt1 + n] =
582  4 * nel + maxN * e->GetGeom()->GetEid(0) + n;
583  }
584  else
585  {
586  vId[cnt1 + np - 1 - n] =
587  4 * nel + maxN * e->GetGeom()->GetEid(0) + n;
588  }
589 
590  // Edge 1
591  edgeOrient = e->GetGeom()->GetEorient(1);
592  if (edgeOrient == StdRegions::eForwards)
593  {
594  edge1 += np - n;
595  vId[cnt1 + np - 1 + edge1] =
596  4 * nel + maxN * e->GetGeom()->GetEid(1) + n;
597  }
598  else
599  {
600  edge1 += n;
601  vId[cnt1 + newpoints - 1 - edge1] =
602  4 * nel + maxN * e->GetGeom()->GetEid(1) + n;
603  }
604 
605  // Edge 2
606  edgeOrient = e->GetGeom()->GetEorient(2);
607  if (edgeOrient == StdRegions::eForwards)
608  {
609  edge2 += n + 1;
610  vId[cnt1 + newpoints - 1 - edge2] =
611  4 * nel + maxN * e->GetGeom()->GetEid(2) + n;
612  }
613  else
614  {
615  edge2 += np + 1 - n;
616  vId[cnt1 + edge2] =
617  4 * nel + maxN * e->GetGeom()->GetEid(2) + n;
618  }
619  }
620 
621  // Interior numbering
622  edge2 = 0;
623  for (int n = 1; n < np - 1; n++)
624  {
625  edge2 += np + 1 - n;
626  for (int m = 1; m < np - n - 1; m++)
627  {
628  vId[cnt1 + edge2 + m] =
629  4 * nel + maxN * 4 * nel + cnt2;
630  cnt2++;
631  }
632  }
633  cnt1 += newpoints;
634  }
635  break;
636  case LibUtilities::eQuadrilateral:
637  {
638  int newpoints =
639  LibUtilities::StdQuadData::getNumberOfCoefficients(np,
640  np);
641  // Vertex numbering
642  vId[cnt1] = e->GetGeom()->GetVid(0);
643  vId[cnt1 + np - 1] = e->GetGeom()->GetVid(1);
644  vId[cnt1 + np * np - 1] = e->GetGeom()->GetVid(2);
645  vId[cnt1 + np * (np - 1)] = e->GetGeom()->GetVid(3);
646 
647  // Edge numbering
648  StdRegions::Orientation edgeOrient;
649  for (int n = 1; n < np - 1; n++)
650  {
651  // Edge 0
652  edgeOrient = e->GetGeom()->GetEorient(0);
653  if (edgeOrient == StdRegions::eForwards)
654  {
655  vId[cnt1 + n] =
656  4 * nel + maxN * e->GetGeom()->GetEid(0) + n;
657  }
658  else
659  {
660  vId[cnt1 + np - 1 - n] =
661  4 * nel + maxN * e->GetGeom()->GetEid(0) + n;
662  }
663 
664  // Edge 1
665  edgeOrient = e->GetGeom()->GetEorient(1);
666  if (edgeOrient == StdRegions::eForwards)
667  {
668  vId[cnt1 + np - 1 + n * np] =
669  4 * nel + maxN * e->GetGeom()->GetEid(1) + n;
670  }
671  else
672  {
673  vId[cnt1 + np * np - 1 - n * np] =
674  4 * nel + maxN * e->GetGeom()->GetEid(1) + n;
675  }
676 
677  // Edge 2
678  edgeOrient = e->GetGeom()->GetEorient(2);
679  if (edgeOrient == StdRegions::eForwards)
680  {
681  vId[cnt1 + np * np - 1 - n] =
682  4 * nel + maxN * e->GetGeom()->GetEid(2) + n;
683  }
684  else
685  {
686  vId[cnt1 + np * (np - 1) + n] =
687  4 * nel + maxN * e->GetGeom()->GetEid(2) + n;
688  }
689 
690  // Edge 3
691  edgeOrient = e->GetGeom()->GetEorient(3);
692  if (edgeOrient == StdRegions::eForwards)
693  {
694  vId[cnt1 + np * (np - 1) - n * np] =
695  4 * nel + maxN * e->GetGeom()->GetEid(3) + n;
696  }
697  else
698  {
699  vId[cnt1 + n * np] =
700  4 * nel + maxN * e->GetGeom()->GetEid(3) + n;
701  }
702  }
703 
704  // Interior numbering
705  for (int n = 1; n < np - 1; n++)
706  {
707  for (int m = 1; m < np - 1; m++)
708  {
709  vId[cnt1 + m * np + n] =
710  4 * nel + maxN * 4 * nel + cnt2;
711  cnt2++;
712  }
713  }
714  cnt1 += newpoints;
715  }
716  break;
717  default:
718  {
719  ASSERTL0(false, "Points not known");
720  }
721  }
722  }
723 
724  // Crete new connectivity using homogeneous information
725  for (int i = 0; i < nel; ++i)
726  {
727  int nTris = oldConn[i].size() / 3;
728  // Create array for new connectivity
729  // (3 tetrahedra between each plane for each triangle)
730  conn = Array<OneD, int>(4 * 3 * nTris * (nPlanes - 1));
731  cnt2 = 0;
732  for (int n = 0; n < nTris; n++)
733  {
734  // Get id of vertices in this triangle (on plane 0)
735  int vId0 = vId[oldConn[i][3 * n + 0]];
736  int vId1 = vId[oldConn[i][3 * n + 1]];
737  int vId2 = vId[oldConn[i][3 * n + 2]];
738 
739  // Determine ordering based on global Id of pts
740  Array<OneD, int> rot(3);
741  if ((vId0 < vId1) && (vId0 < vId2))
742  {
743  rot[0] = 0;
744  if (vId1 < vId2)
745  {
746  rot[1] = 1;
747  rot[2] = 2;
748  }
749  else
750  {
751  rot[1] = 2;
752  rot[2] = 1;
753  }
754  }
755  else if ((vId1 < vId0) && (vId1 < vId2))
756  {
757  rot[0] = 1;
758  if (vId0 < vId2)
759  {
760  rot[1] = 0;
761  rot[2] = 2;
762  }
763  else
764  {
765  rot[1] = 2;
766  rot[2] = 0;
767  }
768  }
769  else
770  {
771  rot[0] = 2;
772  if (vId0 < vId1)
773  {
774  rot[1] = 0;
775  rot[2] = 1;
776  }
777  else
778  {
779  rot[1] = 1;
780  rot[2] = 0;
781  }
782  }
783 
784  // Define new connectivity with tetrahedra
785  for (int p = 0; p < nPlanes - 1; p++)
786  {
787  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[0]];
788  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[1]];
789  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[2]];
790  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[2]];
791 
792  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[0]];
793  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[1]];
794  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[2]];
795  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[1]];
796 
797  conn[cnt2++] = oldConn[i + p * nel][3 * n + rot[0]];
798  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[1]];
799  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[2]];
800  conn[cnt2++] = oldConn[i + (p + 1) * nel][3 * n + rot[0]];
801  }
802  }
803  newConn.push_back(conn);
804  }
805  }
806  else
807  {
808  ASSERTL0(false, "Points type not supported.");
809  }
810 
811  m_f->m_fieldPts->SetConnectivity(newConn);
812 }
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:138
Nektar::LibUtilities::StdTriData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:114
Nektar::LibUtilities::eFourierEvenlySpaced
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
Nektar::LibUtilities::eFourier
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:384
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

References ASSERTL0, Nektar::StdRegions::eForwards, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::eFourierEvenlySpaced, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::FieldUtils::Module::m_f, CellMLToNektar.cellml_metadata::p, Vmath::Sadd(), and Vmath::Vcopy().

Referenced by v_Process().

◆ v_GetModuleDescription()

virtual std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModuleDescription ( )
inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 70 of file ProcessEquiSpacedOutput.h.

71  {
72  return "Interpolating fields to equispaced";
73  }

◆ v_GetModuleName()

virtual std::string Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModuleName ( )
inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 65 of file ProcessEquiSpacedOutput.h.

66  {
67  return "ProcessEquiSpacedOutput";
68  }

◆ v_GetModulePriority()

virtual ModulePriority Nektar::FieldUtils::ProcessEquiSpacedOutput::v_GetModulePriority ( )
inlineoverrideprotectedvirtual

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 75 of file ProcessEquiSpacedOutput.h.

76  {
77  return eConvertExpToPts;
78  }

References Nektar::FieldUtils::eConvertExpToPts.

◆ v_Process()

void Nektar::FieldUtils::ProcessEquiSpacedOutput::v_Process ( po::variables_map &  vm)
overrideprotectedvirtual

Write mesh to output file.

Reimplemented from Nektar::FieldUtils::Module.

Definition at line 74 of file ProcessEquiSpacedOutput.cpp.

75 {
76  m_f->SetUpExp(vm);
77 
78  int nel = m_f->m_exp[0]->GetExpSize();
79  if (!nel)
80  {
81  // Create empty PtsField
82  int nfields = m_f->m_variables.size();
83  int coordim = 3;
84 
85  Array<OneD, Array<OneD, NekDouble>> pts(nfields + coordim);
86  for (int i = 0; i < nfields + coordim; ++i)
87  {
88  pts[i] = Array<OneD, NekDouble>(0);
89  }
90  vector<Array<OneD, int>> ptsConn;
91 
92  m_f->m_fieldPts =
93  MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
94  coordim, m_f->m_variables, pts);
95  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
96  m_f->m_fieldPts->SetConnectivity(ptsConn);
97  return;
98  }
99 
100  int coordim = m_f->m_exp[0]->GetCoordim(0);
101  int shapedim = m_f->m_exp[0]->GetExp(0)->GetShapeDimension();
102  int npts = m_f->m_exp[0]->GetTotPoints();
103  Array<OneD, Array<OneD, NekDouble>> coords(3);
104 
105  // Check if we have a homogeneous expansion
106  bool homogeneous1D = false;
107  if (m_f->m_numHomogeneousDir == 1)
108  {
109  coordim++;
110  shapedim++;
111  homogeneous1D = true;
112  }
113  else if (m_f->m_numHomogeneousDir == 2)
114  {
115  ASSERTL0(false, "Homegeneous2D case not supported");
116  }
117 
118  // set up the number of points in each element
119  int newpoints = 0;
120  int newtotpoints = 0;
121 
122  Array<OneD, int> conn;
123  int prevNcoeffs = 0;
124  int prevNpoints = 0;
125  int cnt = 0;
126 
127  // identify face 1 connectivity for prisms
128  map<int, StdRegions::Orientation> face0orient;
129  set<int> prismorient;
130  LocalRegions::ExpansionSharedPtr e;
131 
132  // prepare PtsField
133  vector<int> ppe;
134  vector<Array<OneD, int>> ptsConn;
135  int nfields;
136 
137  for (int i = 0; i < nel; ++i)
138  {
139  e = m_f->m_exp[0]->GetExp(i);
140  if (e->DetShapeType() == LibUtilities::ePrism)
141  {
142  StdRegions::Orientation forient = e->GetTraceOrient(0);
143  int fid = e->GetGeom()->GetFid(0);
144  if (face0orient.count(fid))
145  { // face 1 meeting face 1 so reverse this id
146  prismorient.insert(i);
147  }
148  else
149  {
150  // just store if Dir 1 is fwd or bwd
151  if ((forient == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
152  (forient == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
153  (forient == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
154  (forient == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
155  {
156  face0orient[fid] = StdRegions::eBwd;
157  }
158  else
159  {
160  face0orient[fid] = StdRegions::eFwd;
161  }
162  }
163  }
164  }
165 
166  for (int i = 0; i < nel; ++i)
167  {
168  e = m_f->m_exp[0]->GetExp(i);
169  if (e->DetShapeType() == LibUtilities::ePrism)
170  {
171  int fid = e->GetGeom()->GetFid(2);
172  // check to see if face 2 meets face 1
173  if (face0orient.count(fid))
174  {
175  // check to see how face 2 is orientated
176  StdRegions::Orientation forient2 = e->GetTraceOrient(2);
177  StdRegions::Orientation forient0 = face0orient[fid];
178 
179  // If dir 1 or forient2 is bwd then check agains
180  // face 1 value
181  if ((forient2 == StdRegions::eDir1BwdDir1_Dir2FwdDir2) ||
182  (forient2 == StdRegions::eDir1BwdDir1_Dir2BwdDir2) ||
183  (forient2 == StdRegions::eDir1BwdDir2_Dir2FwdDir1) ||
184  (forient2 == StdRegions::eDir1BwdDir2_Dir2BwdDir1))
185  {
186  if (forient0 == StdRegions::eFwd)
187  {
188  prismorient.insert(i);
189  }
190  }
191  else
192  {
193  if (forient0 == StdRegions::eBwd)
194  {
195  prismorient.insert(i);
196  }
197  }
198  }
199  }
200  }
201 
202  for (int i = 0; i < nel; ++i)
203  {
204  e = m_f->m_exp[0]->GetExp(i);
205  if (m_config["tetonly"].m_beenSet && m_config["tetonly"].as<bool>())
206  {
207  if (m_f->m_exp[0]->GetExp(i)->DetShapeType() !=
208  LibUtilities::eTetrahedron)
209  {
210  continue;
211  }
212  }
213 
214  switch (e->DetShapeType())
215  {
216  case LibUtilities::eSegment:
217  {
218  int npoints0 = e->GetBasis(0)->GetNumPoints();
219 
220  newpoints =
221  LibUtilities::StdSegData::getNumberOfCoefficients(npoints0);
222  }
223  break;
224  case LibUtilities::eTriangle:
225  {
226  int np0 = e->GetBasis(0)->GetNumPoints();
227  int np1 = e->GetBasis(1)->GetNumPoints();
228  int np = max(np0, np1);
229  newpoints =
230  LibUtilities::StdTriData::getNumberOfCoefficients(np, np);
231  }
232  break;
233  case LibUtilities::eQuadrilateral:
234  {
235  int np0 = e->GetBasis(0)->GetNumPoints();
236  int np1 = e->GetBasis(1)->GetNumPoints();
237  int np = max(np0, np1);
238 
239  newpoints =
240  LibUtilities::StdQuadData::getNumberOfCoefficients(np, np);
241  }
242  break;
243  case LibUtilities::eTetrahedron:
244  {
245  int np0 = e->GetBasis(0)->GetNumPoints();
246  int np1 = e->GetBasis(1)->GetNumPoints();
247  int np2 = e->GetBasis(2)->GetNumPoints();
248  int np = max(np0, max(np1, np2));
249 
250  newpoints = LibUtilities::StdTetData::getNumberOfCoefficients(
251  np, np, np);
252  }
253  break;
254  case LibUtilities::ePrism:
255  {
256  int np0 = e->GetBasis(0)->GetNumPoints();
257  int np1 = e->GetBasis(1)->GetNumPoints();
258  int np2 = e->GetBasis(2)->GetNumPoints();
259  int np = max(np0, max(np1, np2));
260 
261  newpoints = LibUtilities::StdPrismData::getNumberOfCoefficients(
262  np, np, np);
263  }
264  break;
265  case LibUtilities::ePyramid:
266  {
267  int np0 = e->GetBasis(0)->GetNumPoints();
268  int np1 = e->GetBasis(1)->GetNumPoints();
269  int np2 = e->GetBasis(2)->GetNumPoints();
270  int np = max(np0, max(np1, np2));
271 
272  newpoints = LibUtilities::StdPyrData::getNumberOfCoefficients(
273  np, np, np);
274  }
275  break;
276  case LibUtilities::eHexahedron:
277  {
278  int np0 = e->GetBasis(0)->GetNumPoints();
279  int np1 = e->GetBasis(1)->GetNumPoints();
280  int np2 = e->GetBasis(2)->GetNumPoints();
281  int np = max(np0, max(np1, np2));
282 
283  newpoints = LibUtilities::StdHexData::getNumberOfCoefficients(
284  np, np, np);
285  }
286  break;
287  default:
288  {
289  ASSERTL0(false, "Points not known");
290  }
291  }
292 
293  ppe.push_back(newpoints);
294  newtotpoints += newpoints;
295 
296  if (e->DetShapeType() == LibUtilities::ePrism)
297  {
298  bool standard = true;
299 
300  if (prismorient.count(i))
301  {
302  standard = false; // reverse direction
303  }
304 
305  e->GetSimplexEquiSpacedConnectivity(conn, standard);
306  }
307  else
308  {
309 
310  if ((prevNcoeffs != e->GetNcoeffs()) ||
311  (prevNpoints != e->GetTotPoints()))
312  {
313  prevNcoeffs = e->GetNcoeffs();
314  prevNpoints = e->GetTotPoints();
315 
316  e->GetSimplexEquiSpacedConnectivity(conn);
317  }
318  }
319  Array<OneD, int> newconn(conn.size());
320  for (int j = 0; j < conn.size(); ++j)
321  {
322  newconn[j] = conn[j] + cnt;
323  }
324 
325  ptsConn.push_back(newconn);
326  cnt += newpoints;
327  }
328 
329  nfields = m_f->m_variables.size();
330 
331  Array<OneD, Array<OneD, NekDouble>> pts(nfields + coordim);
332 
333  for (int i = 0; i < nfields + coordim; ++i)
334  {
335  pts[i] = Array<OneD, NekDouble>(newtotpoints);
336  }
337 
338  // Interpolate coordinates
339  for (int i = 0; i < coordim; ++i)
340  {
341  coords[i] = Array<OneD, NekDouble>(npts);
342  }
343 
344  for (int i = coordim; i < 3; ++i)
345  {
346  coords[i] = NullNekDouble1DArray;
347  }
348 
349  m_f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
350 
351  Array<OneD, NekDouble> tmp;
352 
353  for (int n = 0; n < coordim; ++n)
354  {
355  cnt = 0;
356  int cnt1 = 0;
357  for (int i = 0; i < nel; ++i)
358  {
359  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
360  coords[n] + cnt, tmp = pts[n] + cnt1);
361  cnt1 += ppe[i];
362  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
363  }
364  }
365 
366  for (int n = 0; n < m_f->m_variables.size(); ++n)
367  {
368  cnt = 0;
369  int cnt1 = 0;
370 
371  if (m_config["modalenergy"].m_beenSet &&
372  m_config["modalenergy"].as<bool>())
373  {
374  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
375  for (int i = 0; i < nel; ++i)
376  {
377  GenOrthoModes(i, phys + cnt, tmp = pts[coordim + n] + cnt1);
378  cnt1 += ppe[i];
379  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
380  }
381  }
382  else
383  {
384  Array<OneD, const NekDouble> phys = m_f->m_exp[n]->GetPhys();
385  for (int i = 0; i < nel; ++i)
386  {
387  m_f->m_exp[0]->GetExp(i)->PhysInterpToSimplexEquiSpaced(
388  phys + cnt, tmp = pts[coordim + n] + cnt1);
389  cnt1 += ppe[i];
390  cnt += m_f->m_exp[0]->GetExp(i)->GetTotPoints();
391  }
392  }
393  }
394 
395  m_f->m_fieldPts = MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(
396  coordim, m_f->m_variables, pts);
397  if (shapedim == 1)
398  {
399  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsSegBlock);
400  }
401 
402  if (shapedim == 2)
403  {
404  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTriBlock);
405  }
406  else if (shapedim == 3)
407  {
408  m_f->m_fieldPts->SetPtsType(LibUtilities::ePtsTetBlock);
409  }
410  m_f->m_fieldPts->SetConnectivity(ptsConn);
411  if (homogeneous1D)
412  {
413  SetHomogeneousConnectivity();
414  }
415 
416  // Save points per element in the point field
417  m_f->m_fieldPts->SetPointsPerElement(ppe);
418 
419  // Clear m_exp
420  m_f->m_exp = vector<MultiRegions::ExpListSharedPtr>();
421 }
Nektar::FieldUtils::ProcessEquiSpacedOutput::GenOrthoModes
void GenOrthoModes(int n, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &coeffs)
Definition: ProcessEquiSpacedOutput.cpp:814
Nektar::LibUtilities::StdHexData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:158
Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:284
Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:237
Nektar::LibUtilities::StdTetData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:192
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:900

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::StdRegions::eBwd, Nektar::StdRegions::eDir1BwdDir1_Dir2BwdDir2, Nektar::StdRegions::eDir1BwdDir1_Dir2FwdDir2, Nektar::StdRegions::eDir1BwdDir2_Dir2BwdDir1, Nektar::StdRegions::eDir1BwdDir2_Dir2FwdDir1, Nektar::StdRegions::eFwd, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePtsSegBlock, Nektar::LibUtilities::ePtsTetBlock, Nektar::LibUtilities::ePtsTriBlock, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTetrahedron, Nektar::LibUtilities::eTriangle, GenOrthoModes(), Nektar::LibUtilities::StdSegData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTriData::getNumberOfCoefficients(), Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients(), Nektar::LibUtilities::StdHexData::getNumberOfCoefficients(), Nektar::LibUtilities::StdTetData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients(), Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients(), Nektar::FieldUtils::Module::m_config, Nektar::FieldUtils::Module::m_f, Nektar::NullNekDouble1DArray, and SetHomogeneousConnectivity().

Member Data Documentation

◆ className

ModuleKey Nektar::FieldUtils::ProcessEquiSpacedOutput::className
static
Initial value:
=
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "equispacedoutput"),
ProcessEquiSpacedOutput::create,
"Write data as equi-spaced output using simplices to represent the "
"data for connecting points")
Nektar::FieldUtils::ProcessEquiSpacedOutput::create
static std::shared_ptr< Module > create(FieldSharedPtr f)
Creates an instance of this class.
Definition: ProcessEquiSpacedOutput.h:52
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::FieldUtils::ModuleKey
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:317
Nektar::FieldUtils::GetModuleFactory
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:49

Definition at line 56 of file ProcessEquiSpacedOutput.h.

Generated on Sat Mar 25 2023 12:35:57 for Nektar++ by doxygen 1.9.1