Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Private Member Functions | Private Attributes | Friends | List of all members
Nektar::SolverUtils::FilterAeroForces Class Reference

#include <FilterAeroForces.h>

Inheritance diagram for Nektar::SolverUtils::FilterAeroForces:
Inheritance graph
[legend]
Collaboration diagram for Nektar::SolverUtils::FilterAeroForces:
Collaboration graph
[legend]

Public Member Functions

SOLVER_UTILS_EXPORT FilterAeroForces (const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
 
virtual SOLVER_UTILS_EXPORT ~FilterAeroForces ()
 
SOLVER_UTILS_EXPORT void GetForces (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)
 
- Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT Filter (const LibUtilities::SessionReaderSharedPtr &pSession)
 
virtual SOLVER_UTILS_EXPORT ~Filter ()
 
SOLVER_UTILS_EXPORT void Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT void Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT void Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT bool IsTimeDependent ()
 

Static Public Member Functions

static FilterSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of the class. More...
 

Protected Member Functions

virtual void v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
virtual void v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
virtual void v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
virtual bool v_IsTimeDependent ()
 

Private Member Functions

void CalculateForces (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
void CalculateForcesMapping (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 

Private Attributes

std::vector< unsigned int > m_boundaryRegionsIdList
 ID's of boundary regions where we want the forces. More...
 
std::vector< bool > m_boundaryRegionIsInList
 Determines if a given Boundary Region is in m_boundaryRegionsIdList. More...
 
unsigned int m_index
 
unsigned int m_outputFrequency
 
bool m_outputAllPlanes
 if using a homogeneous1D expansion, determine if should output all planes or just the average More...
 
bool m_isHomogeneous1D
 
std::string m_outputFile
 
std::ofstream m_outputStream
 
LibUtilities::BasisSharedPtr m_homogeneousBasis
 
std::string m_BoundaryString
 
Array< OneD, int > m_BCtoElmtID
 
Array< OneD, int > m_BCtoTraceID
 
int m_nPlanes
 number of planes for homogeneous1D expansion More...
 
Array< OneD, int > m_planesID
 
NekDouble m_startTime
 
Array< OneD, Array< OneD,
NekDouble > > 
m_directions
 
Array< OneD, Array< OneD,
NekDouble > > 
m_Fpplane
 
Array< OneD, Array< OneD,
NekDouble > > 
m_Fvplane
 
Array< OneD, Array< OneD,
NekDouble > > 
m_Ftplane
 
NekDouble m_lastTime
 
GlobalMapping::MappingSharedPtr m_mapping
 

Friends

class MemoryManager< FilterAeroForces >
 

Additional Inherited Members

- Public Types inherited from Nektar::SolverUtils::Filter
typedef std::map< std::string,
std::string > 
ParamMap
 
- Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr m_session
 

Detailed Description

Definition at line 47 of file FilterAeroForces.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::FilterAeroForces::FilterAeroForces ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::map< std::string, std::string > &  pParams 
)

Definition at line 60 of file FilterAeroForces.cpp.

References ASSERTL0, Nektar::LibUtilities::Equation::Evaluate(), m_BoundaryString, m_directions, m_isHomogeneous1D, m_outputAllPlanes, m_outputFile, m_outputFrequency, Nektar::SolverUtils::Filter::m_session, and m_startTime.

62  :
63  Filter(pSession)
64 {
65  ParamMap::const_iterator it;
66 
67  // OutputFile
68  it = pParams.find("OutputFile");
69  if (it == pParams.end())
70  {
71  m_outputFile = m_session->GetSessionName();
72  }
73  else
74  {
75  ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
76  m_outputFile = it->second;
77  }
78  if (!(m_outputFile.length() >= 4
79  && m_outputFile.substr(m_outputFile.length() - 4) == ".fce"))
80  {
81  m_outputFile += ".fce";
82  }
83 
84  // OutputFrequency
85  it = pParams.find("OutputFrequency");
86  if (it == pParams.end())
87  {
89  }
90  else
91  {
92  LibUtilities::Equation equ(m_session, it->second);
93  m_outputFrequency = floor(equ.Evaluate());
94  }
95 
96  // Time after which we need to calculate the forces
97  it = pParams.find("StartTime");
98  if (it == pParams.end())
99  {
100  m_startTime = 0;
101  }
102  else
103  {
104  LibUtilities::Equation equ(m_session, it->second);
105  m_startTime = equ.Evaluate();
106  }
107 
108  // For 3DH1D, OutputAllPlanes or just average forces?
109  m_session->MatchSolverInfo("Homogeneous", "1D",
110  m_isHomogeneous1D, false);
112  {
113  it = pParams.find("OutputAllPlanes");
114  if (it == pParams.end())
115  {
116  m_outputAllPlanes = false;
117  }
118  else
119  {
120  std::string sOption =
121  it->second.c_str();
122  m_outputAllPlanes = ( boost::iequals(sOption,"true")) ||
123  ( boost::iequals(sOption,"yes"));
124  }
125  }
126  else
127  {
128  m_outputAllPlanes = false;
129  }
130 
131  // Boundary (to calculate forces on)
132  it = pParams.find("Boundary");
133  ASSERTL0(it != pParams.end(), "Missing parameter 'Boundary");
134  ASSERTL0(it->second.length() > 0, "Empty parameter 'Boundary'.");
135  m_BoundaryString = it->second;
136 
137  //
138  // Directions (to project forces)
139  //
140 
141  // Allocate m_directions
142  m_directions = Array<OneD, Array<OneD, NekDouble> > (3);
143  //Initialise directions to default values (ex, ey, ez)
144  for (int i = 0; i < 3; ++i)
145  {
146  m_directions[i] = Array<OneD, NekDouble>(3, 0.0);
147  m_directions[i][i] = 1.0;
148  }
149  std::stringstream directionStream;
150  std::string directionString;
151  //Override with input from xml file (if defined)
152  for (int i = 0; i < 3; ++i)
153  {
154  std::stringstream tmp;
155  tmp << i+1;
156  std::string dir = "Direction" + tmp.str();
157  it = pParams.find(dir);
158  if ( it != pParams.end() )
159  {
160  ASSERTL0(!(it->second.empty()),
161  "Missing parameter '"+dir+"'.");
162  directionStream.str(it->second);
163  // Guarantee the stream is in its start position
164  // before extracting
165  directionStream.clear();
166  // normalisation factor
167  NekDouble norm = 0.0;
168  for (int j = 0; j < 3; j++)
169  {
170  directionStream >> directionString;
171  if (!directionString.empty())
172  {
173  LibUtilities::Equation equ(m_session, directionString);
174  m_directions[i][j] = equ.Evaluate();
175  norm += m_directions[i][j]*m_directions[i][j];
176  }
177  }
178  // Normalise direction
179  for( int j = 0; j < 3; j++)
180  {
181  m_directions[i][j] /= sqrt(norm);
182  }
183  }
184  }
185 
186 
187 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
SOLVER_UTILS_EXPORT Filter(const LibUtilities::SessionReaderSharedPtr &pSession)
Definition: Filter.cpp:52
Array< OneD, Array< OneD, NekDouble > > m_directions
double NekDouble
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...
Nektar::SolverUtils::FilterAeroForces::~FilterAeroForces ( )
virtual

Definition at line 193 of file FilterAeroForces.cpp.

194 {
195 
196 }

Member Function Documentation

void Nektar::SolverUtils::FilterAeroForces::CalculateForces ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
private

This function calculates the forces

Definition at line 558 of file FilterAeroForces.cpp.

References ASSERTL0, CalculateForcesMapping(), m_BCtoElmtID, m_BCtoTraceID, m_boundaryRegionIsInList, m_directions, m_Fpplane, m_Ftplane, m_Fvplane, m_isHomogeneous1D, m_lastTime, m_mapping, m_nPlanes, m_outputAllPlanes, m_planesID, Nektar::SolverUtils::Filter::m_session, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by GetForces(), and v_Update().

561 {
562  // Store time so we can check if result is up-to-date
563  m_lastTime = time;
564 
565  // If a mapping is defined, call specific function
566  // Note: CalculateForcesMapping should work without a mapping,
567  // but since it is not very efficient the way it is now,
568  // it is only used when actually required
569  if (m_mapping->IsDefined())
570  {
571  CalculateForcesMapping( pFields, time);
572  return;
573  }
574 
575  int i, j, k, n, cnt, elmtid, nq, offset, boundary, plane;
576  // Get number of quadrature points and dimensions
577  int physTot = pFields[0]->GetNpoints();
578  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
579  int nVel = expdim;
580  if( m_isHomogeneous1D )
581  {
582  nVel = nVel + 1;
583  }
584 
586 
587  // Fields used to calculate forces (a single plane for 3DH1D)
588  Array<OneD, MultiRegions::ExpListSharedPtr>
589  fields( pFields.num_elements() );
590 
591  // Arrays of variables in the element
592  Array<OneD, Array<OneD, NekDouble> > velocity(expdim);
593  Array<OneD, NekDouble> P(physTot);
594 
595  // Velocity gradient
596  Array<OneD, Array<OneD, NekDouble> > grad( expdim*expdim);
597 
598  // Values at the boundary
599  Array<OneD, NekDouble> Pb;
600  Array<OneD, Array<OneD, NekDouble> > gradb( expdim*expdim);
601 
602  // Communicators to exchange results
603  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
604  LibUtilities::CommSharedPtr rowComm = vComm->GetRowComm();
605  LibUtilities::CommSharedPtr colComm =
606  m_session->DefinesSolverInfo("HomoStrip") ?
607  vComm->GetColumnComm()->GetColumnComm():
608  vComm->GetColumnComm();
609 
610  // Arrays with forces in each plane
611  m_Fpplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
612  m_Fvplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
613  m_Ftplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
614  for( i = 0; i < expdim; i++)
615  {
616  m_Fpplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
617  m_Fvplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
618  m_Ftplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
619  }
620 
621  // Forces per element length in a boundary
622  Array<OneD, Array<OneD, NekDouble> > fp( expdim );
623  Array<OneD, Array<OneD, NekDouble> > fv( expdim );
624 
625  // Get viscosity
626  NekDouble rho = (m_session->DefinesParameter("rho"))
627  ? (m_session->GetParameter("rho"))
628  : 1;
629  NekDouble mu = rho*m_session->GetParameter("Kinvis");
630 
631  // Perform BwdTrans: when we only want the mean force in a 3DH1D
632  // we work in wavespace, otherwise we use physical space
633  for(i = 0; i < pFields.num_elements(); ++i)
634  {
636  {
637  pFields[i]->SetWaveSpace(false);
638  }
639  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
640  pFields[i]->UpdatePhys());
641  pFields[i]->SetPhysState(true);
642  }
643 
644  // Define boundary expansions
645  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
646  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
648  {
649  BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
650  BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
651  }
652  else
653  {
654  BndConds = pFields[0]->GetBndConditions();
655  BndExp = pFields[0]->GetBndCondExpansions();
656  }
657 
658  // For Homogeneous, calculate force on each 2D plane
659  // Otherwise, m_nPlanes = 1, and loop only runs once
660  for(plane = 0; plane < m_nPlanes; plane++ )
661  {
662  // Check if plane is in this proc
663  if( m_planesID[plane] != -1 )
664  {
665  // For Homogeneous, consider the 2D expansion
666  // on the current plane
668  {
669  for(n = 0; n < pFields.num_elements(); n++)
670  {
671  fields[n] = pFields[n]->GetPlane(m_planesID[plane]);
672  }
673  }
674  else
675  {
676  for(n = 0; n < pFields.num_elements(); n++)
677  {
678  fields[n] = pFields[n];
679  }
680  }
681 
682  //Loop all the Boundary Regions
683  for( cnt = n = 0; n < BndConds.num_elements(); n++)
684  {
685  if(m_boundaryRegionIsInList[n] == 1)
686  {
687  for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
688  {
689  elmtid = m_BCtoElmtID[cnt];
690  elmt = fields[0]->GetExp(elmtid);
691  nq = elmt->GetTotPoints();
692  offset = fields[0]->GetPhys_Offset(elmtid);
693 
694  // Extract fields on this element
695  for( j=0; j<expdim; j++)
696  {
697  velocity[j] = fields[j]->GetPhys() + offset;
698  }
699  P = fields[nVel]->GetPhys() + offset;
700 
701  // Compute the velocity gradients
702  for (j=0; j<expdim; j++)
703  {
704  for (k=0; k<expdim; k++)
705  {
706  grad[j*expdim+k] =
707  Array<OneD, NekDouble>(nq,0.0);
708  elmt->PhysDeriv(k,velocity[j],
709  grad[j*expdim+k]);
710  }
711  }
712 
713  // identify boundary of element
714  boundary = m_BCtoTraceID[cnt];
715 
716  // Dimension specific part for obtaining values
717  // at boundary and normal vector
718  Array<OneD, Array<OneD, NekDouble> > normals;
719  int nbc;
720  switch(expdim)
721  {
722  case 2:
723  {
724  // Get expansion on boundary
726  bc = BndExp[n]->GetExp(i)->
727  as<LocalRegions::Expansion1D> ();
728 
729  // Get number of points on the boundary
730  nbc = bc->GetTotPoints();
731 
732  // Get normals
733  normals = elmt->GetEdgeNormal(boundary);
734 
735  // Extract values at boundary
736  Pb = Array<OneD, NekDouble>(nbc,0.0);
737  elmt->GetEdgePhysVals(boundary,bc,P,Pb);
738  for(int j = 0; j < expdim*expdim; ++j)
739  {
740  gradb[j] = Array<OneD, NekDouble>
741  (nbc,0.0);
742  elmt->GetEdgePhysVals(boundary,
743  bc,grad[j],gradb[j]);
744  }
745  }
746  break;
747 
748  case 3:
749  {
750  // Get expansion on boundary
752  bc = BndExp[n]->GetExp(i)->
753  as<LocalRegions::Expansion2D> ();
754 
755  // Get number of points on the boundary
756  nbc = bc->GetTotPoints();
757 
758  // Get normals
759  normals = elmt->GetFaceNormal(boundary);
760 
761  // Extract values at boundary
762  Pb = Array<OneD, NekDouble>(nbc,0.0);
763  elmt->GetFacePhysVals(boundary,bc,P,Pb);
764  for(int j = 0; j < expdim*expdim; ++j)
765  {
766  gradb[j] = Array<OneD, NekDouble>
767  (nbc,0.0);
768  elmt->GetFacePhysVals(boundary,
769  bc,grad[j],gradb[j]);
770  }
771  }
772  break;
773 
774  default:
775  ASSERTL0(false,
776  "Expansion not supported by FilterForces");
777  break;
778  }
779 
780  // Calculate forces per unit length
781 
782  // Pressure component: fp[j] = p*n[j]
783  for ( j = 0; j < expdim; j++)
784  {
785  fp[j] = Array<OneD, NekDouble> (nbc,0.0);
786  Vmath::Vmul (nbc, Pb, 1,
787  normals[j], 1,
788  fp[j], 1);
789  }
790 
791  // Viscous component:
792  // fv[j] = -mu*{(grad[k,j]+grad[j,k]) *n[k]}
793  for ( j = 0; j < expdim; j++ )
794  {
795  fv[j] = Array<OneD, NekDouble> (nbc,0.0);
796  for ( k = 0; k < expdim; k++ )
797  {
798  Vmath::Vvtvp (nbc, gradb[k*expdim+j], 1,
799  normals[k], 1,
800  fv[j], 1,
801  fv[j], 1);
802  Vmath::Vvtvp (nbc, gradb[j*expdim+k], 1,
803  normals[k], 1,
804  fv[j], 1,
805  fv[j], 1);
806  }
807  Vmath::Smul(nbc, -mu, fv[j], 1, fv[j], 1);
808  }
809 
810  // Integrate to obtain force
811  for ( j = 0; j < expdim; j++)
812  {
813  m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
814  Integral(fp[j]);
815  m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
816  Integral(fv[j]);
817  }
818  }
819  }
820  else
821  {
822  cnt += BndExp[n]->GetExpSize();
823  }
824  }
825  }
826  }
827 
828  // Combine contributions from different processes
829  // this is split between row and col comm because of
830  // homostrips case, which only keeps its own strip
831  for( i = 0; i < expdim; i++)
832  {
833  rowComm->AllReduce(m_Fpplane[i], LibUtilities::ReduceSum);
834  colComm->AllReduce(m_Fpplane[i], LibUtilities::ReduceSum);
835 
836  rowComm->AllReduce(m_Fvplane[i], LibUtilities::ReduceSum);
837  colComm->AllReduce(m_Fvplane[i], LibUtilities::ReduceSum);
838  }
839 
840  // Project results to new directions
841  for(plane = 0; plane < m_nPlanes; plane++)
842  {
843  Array< OneD, NekDouble> tmpP(expdim, 0.0);
844  Array< OneD, NekDouble> tmpV(expdim, 0.0);
845  for( i = 0; i < expdim; i++)
846  {
847  for( j = 0; j < expdim; j++ )
848  {
849  tmpP[i] += m_Fpplane[j][plane]*m_directions[i][j];
850  tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
851  }
852  }
853  // Copy result
854  for( i = 0; i < expdim; i++)
855  {
856  m_Fpplane[i][plane] = tmpP[i];
857  m_Fvplane[i][plane] = tmpV[i];
858  }
859  }
860 
861  // Sum viscous and pressure components
862  for(plane = 0; plane < m_nPlanes; plane++)
863  {
864  for( i = 0; i < expdim; i++)
865  {
866  m_Ftplane[i][plane] = m_Fpplane[i][plane] + m_Fvplane[i][plane];
867  }
868  }
869 
870  // Put results back to wavespace, if necessary
872  {
873  for (i = 0; i < pFields.num_elements(); ++i)
874  {
875  pFields[i]->SetWaveSpace(true);
876  pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
877  pFields[i]->UpdatePhys());
878  }
879  }
880 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, Array< OneD, NekDouble > > m_directions
double NekDouble
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
int m_nPlanes
number of planes for homogeneous1D expansion
Array< OneD, Array< OneD, NekDouble > > m_Ftplane
Array< OneD, Array< OneD, NekDouble > > m_Fpplane
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, Array< OneD, NekDouble > > m_Fvplane
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...
void CalculateForcesMapping(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
void Nektar::SolverUtils::FilterAeroForces::CalculateForcesMapping ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
private

This function calculates the forces when we have a mapping defining a coordinate system transformation

Definition at line 886 of file FilterAeroForces.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::DirCartesianMap, m_BCtoElmtID, m_BCtoTraceID, m_boundaryRegionIsInList, m_directions, m_Fpplane, m_Ftplane, m_Fvplane, m_isHomogeneous1D, m_mapping, m_nPlanes, m_planesID, Nektar::SolverUtils::Filter::m_session, Nektar::LibUtilities::ReduceSum, Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Vvtvp().

Referenced by CalculateForces().

889 {
890  int i, j, k, n, cnt, elmtid, offset, boundary, plane;
891  // Get number of quadrature points and dimensions
892  int physTot = pFields[0]->GetNpoints();
893  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
894  int nVel = expdim;
895  if( m_isHomogeneous1D )
896  {
897  nVel = nVel + 1;
898  }
899 
901 
902  // Pressure stress tensor
903  // (global, in a plane, in element and boundary)
904  Array<OneD, MultiRegions::ExpListSharedPtr> P ( nVel*nVel);
905  Array<OneD, MultiRegions::ExpListSharedPtr> PPlane ( nVel*nVel);
906  Array<OneD, Array<OneD, NekDouble> > PElmt ( nVel*nVel);
907  Array<OneD, Array<OneD, NekDouble> > PBnd ( nVel*nVel);
908  // Velocity gradient
909  Array<OneD, MultiRegions::ExpListSharedPtr> grad ( nVel*nVel);
910  Array<OneD, MultiRegions::ExpListSharedPtr> gradPlane( nVel*nVel);
911  Array<OneD, Array<OneD, NekDouble> > gradElmt ( nVel*nVel);
912  Array<OneD, Array<OneD, NekDouble> > gradBnd ( nVel*nVel);
913 
914  // Transformation to cartesian system
915  Array<OneD, MultiRegions::ExpListSharedPtr> C ( nVel*nVel);
916  Array<OneD, MultiRegions::ExpListSharedPtr> CPlane( nVel*nVel);
917  Array<OneD, Array<OneD, NekDouble> > CElmt ( nVel*nVel);
918  Array<OneD, Array<OneD, NekDouble> > CBnd ( nVel*nVel);
919 
920  // Jacobian
923  Array<OneD, NekDouble> JacElmt;
924  Array<OneD, NekDouble> JacBnd;
925 
926  // Communicators to exchange results
927  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
928  LibUtilities::CommSharedPtr rowComm = vComm->GetRowComm();
929  LibUtilities::CommSharedPtr colComm =
930  m_session->DefinesSolverInfo("HomoStrip") ?
931  vComm->GetColumnComm()->GetColumnComm():
932  vComm->GetColumnComm();
933 
934  // Arrays with forces in each plane
935  m_Fpplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
936  m_Fvplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
937  m_Ftplane = Array<OneD, Array<OneD, NekDouble> > (expdim);
938  for( i = 0; i < expdim; i++)
939  {
940  m_Fpplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
941  m_Fvplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
942  m_Ftplane[i] = Array<OneD, NekDouble>(m_nPlanes,0.0);
943  }
944 
945  // Forces per element length in a boundary
946  Array<OneD, Array<OneD, NekDouble> > fp( nVel );
947  Array<OneD, Array<OneD, NekDouble> > fv( nVel );
948 
949  // Get viscosity
950  NekDouble rho = (m_session->DefinesParameter("rho"))
951  ? (m_session->GetParameter("rho"))
952  : 1;
953  NekDouble mu = rho*m_session->GetParameter("Kinvis");
954 
955  // Perform BwdTrans: for case with mapping, we can never work
956  // in wavespace
957  for(i = 0; i < pFields.num_elements(); ++i)
958  {
959  if (m_isHomogeneous1D)
960  {
961  pFields[i]->SetWaveSpace(false);
962  }
963  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
964  pFields[i]->UpdatePhys());
965  pFields[i]->SetPhysState(true);
966  }
967 
968  // Define boundary expansions
969  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
970  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
972  {
973  BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
974  BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
975  }
976  else
977  {
978  BndConds = pFields[0]->GetBndConditions();
979  BndExp = pFields[0]->GetBndCondExpansions();
980  }
981 
982  //
983  // Calculate pressure stress tensor, velocity gradient
984  // and get informations about the mapping
985 
986  // Initialise variables
987  switch (expdim)
988  {
989  case 2:
990  {
991  if (m_isHomogeneous1D)
992  {
994  Exp3DH1 = boost::dynamic_pointer_cast
995  <MultiRegions::ExpList3DHomogeneous1D>
996  (pFields[0]);
997  for(i = 0; i < nVel*nVel; i++)
998  {
1000  AllocateSharedPtr(*Exp3DH1);
1001 
1003  AllocateSharedPtr(*Exp3DH1);
1004 
1006  AllocateSharedPtr(*Exp3DH1);
1007  }
1009  AllocateSharedPtr(*Exp3DH1);
1010  }
1011  else
1012  {
1014  Exp2D = boost::dynamic_pointer_cast
1015  <MultiRegions::ExpList2D>
1016  (pFields[0]);
1017  for(i = 0; i < nVel*nVel; i++)
1018  {
1020  AllocateSharedPtr(*Exp2D);
1021 
1023  AllocateSharedPtr(*Exp2D);
1024 
1026  AllocateSharedPtr(*Exp2D);
1027  }
1029  AllocateSharedPtr(*Exp2D);
1030  }
1031  break;
1032  }
1033  case 3:
1034  {
1036  Exp3D = boost::dynamic_pointer_cast
1037  <MultiRegions::ExpList3D>
1038  (pFields[0]);
1039  for(i = 0; i < nVel*nVel; i++)
1040  {
1042  AllocateSharedPtr(*Exp3D);
1043 
1045  AllocateSharedPtr(*Exp3D);
1046 
1048  AllocateSharedPtr(*Exp3D);
1049  }
1051  AllocateSharedPtr(*Exp3D);
1052 
1053  break;
1054  }
1055  default:
1056  ASSERTL0(false,"Expansion dimension not supported by FilterAeroForces");
1057  break;
1058  }
1059 
1060 
1061  // Get g^ij
1062  Array<OneD, Array<OneD, NekDouble> > tmp( nVel*nVel );
1063  m_mapping->GetInvMetricTensor(tmp);
1064 
1065  // Calculate P^ij = g^ij*p
1066  for (i = 0; i < nVel*nVel; i++)
1067  {
1068  Vmath::Vmul(physTot, tmp[i], 1,
1069  pFields[nVel]->GetPhys(), 1,
1070  P[i]->UpdatePhys(), 1);
1071  }
1072 
1073  // Calculate covariant derivatives of U = u^i_,k
1074  Array<OneD, Array<OneD, NekDouble> > wk( nVel );
1075  for (i=0; i<nVel; i++)
1076  {
1077  wk[i] = Array<OneD, NekDouble>(physTot, 0.0);
1078  Vmath::Vcopy(physTot, pFields[i]->GetPhys(), 1,
1079  wk[i], 1);
1080  }
1081  m_mapping->ApplyChristoffelContravar(wk, tmp);
1082  for (i=0; i< nVel; i++)
1083  {
1084  for (k=0; k< nVel; k++)
1085  {
1086  pFields[0]->PhysDeriv(MultiRegions::DirCartesianMap[k],
1087  wk[i], grad[i*nVel+k]->UpdatePhys());
1088 
1089  Vmath::Vadd(physTot,tmp[i*nVel+k],1,
1090  grad[i*nVel+k]->UpdatePhys(), 1,
1091  grad[i*nVel+k]->UpdatePhys(), 1);
1092  }
1093  }
1094  // Raise index to obtain Grad^ij = g^jk u^i_,k
1095  for (i=0; i< nVel; i++)
1096  {
1097  for (k=0; k< nVel; k++)
1098  {
1099  Vmath::Vcopy(physTot, grad[i*nVel+k]->GetPhys(), 1,
1100  wk[k], 1);
1101  }
1102  m_mapping->RaiseIndex(wk, wk);
1103  for (j=0; j<nVel; j++)
1104  {
1105  Vmath::Vcopy(physTot, wk[j], 1,
1106  grad[i*nVel+j]->UpdatePhys(), 1);
1107  }
1108  }
1109 
1110  // Get Jacobian
1111  m_mapping->GetJacobian( Jac->UpdatePhys());
1112 
1113  // Get transformation to Cartesian system
1114  for (i=0; i< nVel; i++)
1115  {
1116  // Zero wk storage
1117  for (k=0; k< nVel; k++)
1118  {
1119  wk[k] = Array<OneD, NekDouble>(physTot, 0.0);
1120  }
1121  // Make wk[i] = 1
1122  wk[i] = Array<OneD, NekDouble>(physTot, 1.0);
1123  // Transform wk to Cartesian
1124  m_mapping->ContravarToCartesian(wk,wk);
1125  // Copy result to a column in C
1126  for (k=0; k< nVel; k++)
1127  {
1128  Vmath::Vcopy(physTot, wk[k], 1,
1129  C[k*nVel+i]->UpdatePhys(), 1);
1130  }
1131  }
1132 
1133  //
1134  // Calculate forces
1135  //
1136 
1137  // For Homogeneous, calculate force on each 2D plane
1138  // Otherwise, m_nPlanes = 1, and loop only runs once
1139  for(plane = 0; plane < m_nPlanes; plane++ )
1140  {
1141  // Check if plane is in this proc
1142  if( m_planesID[plane] != -1 )
1143  {
1144  // For Homogeneous, consider the 2D expansion
1145  // on the current plane
1146  if(m_isHomogeneous1D)
1147  {
1148  for(n = 0; n < nVel*nVel; n++)
1149  {
1150  PPlane[n] = P[n]->GetPlane(m_planesID[plane]);
1151  gradPlane[n] = grad[n]->GetPlane(m_planesID[plane]);
1152  CPlane[n] = C[n]->GetPlane(m_planesID[plane]);
1153  }
1154  JacPlane = Jac->GetPlane(m_planesID[plane]);
1155  }
1156  else
1157  {
1158  for(n = 0; n < nVel*nVel; n++)
1159  {
1160  PPlane[n] = P[n];
1161  gradPlane[n] = grad[n];
1162  CPlane[n] = C[n];
1163  }
1164  JacPlane = Jac;
1165  }
1166 
1167  //Loop all the Boundary Regions
1168  for( cnt = n = 0; n < BndConds.num_elements(); n++)
1169  {
1170  if(m_boundaryRegionIsInList[n] == 1)
1171  {
1172  for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
1173  {
1174  elmtid = m_BCtoElmtID[cnt];
1175  elmt = PPlane[0]->GetExp(elmtid);
1176  offset = PPlane[0]->GetPhys_Offset(elmtid);
1177 
1178  // Extract fields on this element
1179  for( j=0; j<nVel*nVel; j++)
1180  {
1181  PElmt[j] = PPlane[j]->GetPhys()
1182  + offset;
1183  gradElmt[j] = gradPlane[j]->GetPhys()
1184  + offset;
1185  CElmt[j] = CPlane[j]->GetPhys()
1186  + offset;
1187  }
1188  JacElmt = JacPlane->GetPhys() + offset;
1189 
1190  // identify boundary of element
1191  boundary = m_BCtoTraceID[cnt];
1192 
1193  // Dimension specific part for obtaining values
1194  // at boundary and normal vector
1195  Array<OneD, Array<OneD, NekDouble> > normals;
1196  int nbc;
1197  switch(expdim)
1198  {
1199  case 2:
1200  {
1201  // Get expansion on boundary
1203  bc = BndExp[n]->GetExp(i)->
1204  as<LocalRegions::Expansion1D> ();
1205 
1206  // Get number of points on the boundary
1207  nbc = bc->GetTotPoints();
1208 
1209  // Get normals
1210  normals = elmt->GetEdgeNormal(boundary);
1211 
1212  // Extract values at boundary
1213  for(int j = 0; j < nVel*nVel; ++j)
1214  {
1215  gradBnd[j] = Array<OneD, NekDouble>
1216  (nbc,0.0);
1217  elmt->GetEdgePhysVals(boundary,
1218  bc,gradElmt[j],
1219  gradBnd[j]);
1220 
1221  PBnd[j] = Array<OneD, NekDouble>
1222  (nbc,0.0);
1223  elmt->GetEdgePhysVals(boundary,
1224  bc,PElmt[j],
1225  PBnd[j]);
1226  CBnd[j] = Array<OneD, NekDouble>
1227  (nbc,0.0);
1228  elmt->GetEdgePhysVals(boundary,
1229  bc,CElmt[j],
1230  CBnd[j]);
1231  }
1232  JacBnd = Array<OneD, NekDouble>
1233  (nbc,0.0);
1234  elmt->GetEdgePhysVals(boundary,
1235  bc,JacElmt,
1236  JacBnd);
1237  }
1238  break;
1239 
1240  case 3:
1241  {
1242  // Get expansion on boundary
1244  bc = BndExp[n]->GetExp(i)->
1245  as<LocalRegions::Expansion2D> ();
1246 
1247  // Get number of points on the boundary
1248  nbc = bc->GetTotPoints();
1249 
1250  // Get normals
1251  normals = elmt->GetFaceNormal(boundary);
1252 
1253  // Extract values at boundary
1254  for(int j = 0; j < nVel*nVel; ++j)
1255  {
1256  gradBnd[j] = Array<OneD, NekDouble>
1257  (nbc,0.0);
1258  elmt->GetFacePhysVals(boundary,
1259  bc,gradElmt[j],
1260  gradBnd[j]);
1261 
1262  PBnd[j] = Array<OneD, NekDouble>
1263  (nbc,0.0);
1264  elmt->GetFacePhysVals(boundary,
1265  bc,PElmt[j],
1266  PBnd[j]);
1267 
1268  CBnd[j] = Array<OneD, NekDouble>
1269  (nbc,0.0);
1270  elmt->GetFacePhysVals(boundary,
1271  bc,CElmt[j],
1272  CBnd[j]);
1273  }
1274  JacBnd = Array<OneD, NekDouble>
1275  (nbc,0.0);
1276  elmt->GetFacePhysVals(boundary,
1277  bc,JacElmt,
1278  JacBnd);
1279  }
1280  break;
1281 
1282  default:
1283  ASSERTL0(false,
1284  "Expansion not supported by FilterForces");
1285  break;
1286  }
1287 
1288  // Calculate forces per unit length
1289 
1290  // Pressure component: fp[j] = P[j,k]*n[k]
1291  for ( j = 0; j < nVel; j++)
1292  {
1293  fp[j] = Array<OneD, NekDouble> (nbc,0.0);
1294  // Normals only has expdim elements
1295  for ( k = 0; k < expdim; k++)
1296  {
1297  Vmath::Vvtvp (nbc, PBnd[ j*nVel + k], 1,
1298  normals[k], 1,
1299  fp[j], 1,
1300  fp[j], 1);
1301  }
1302  }
1303 
1304  // Viscous component:
1305  // fv[j] = -mu*{(grad[k,j]+grad[j,k]) *n[k]}
1306  for ( j = 0; j < nVel; j++ )
1307  {
1308  fv[j] = Array<OneD, NekDouble> (nbc,0.0);
1309  for ( k = 0; k < expdim; k++ )
1310  {
1311  Vmath::Vvtvp (nbc,gradBnd[k*nVel+j],1,
1312  normals[k], 1,
1313  fv[j], 1,
1314  fv[j], 1);
1315  Vmath::Vvtvp (nbc,gradBnd[j*nVel+k],1,
1316  normals[k], 1,
1317  fv[j], 1,
1318  fv[j], 1);
1319  }
1320  Vmath::Smul(nbc, -mu, fv[j], 1, fv[j], 1);
1321  }
1322 
1323  // Multiply by Jacobian
1324  for ( k = 0; k < nVel; k++ )
1325  {
1326  Vmath::Vmul(nbc, JacBnd, 1, fp[k], 1,
1327  fp[k], 1);
1328  Vmath::Vmul(nbc, JacBnd, 1, fv[k], 1,
1329  fv[k], 1);
1330  }
1331 
1332  // Convert to cartesian system
1333  for ( k = 0; k < nVel; k++ )
1334  {
1335  wk[k] = Array<OneD, NekDouble>(physTot,0.0);
1336  for ( j = 0; j < nVel; j++ )
1337  {
1338  Vmath::Vvtvp(nbc, CBnd[k*nVel+j], 1,
1339  fp[j], 1,
1340  wk[k], 1,
1341  wk[k], 1);
1342  }
1343  }
1344  for ( k = 0; k < nVel; k++ )
1345  {
1346  Vmath::Vcopy(nbc, wk[k], 1, fp[k], 1);
1347  }
1348 
1349  for ( k = 0; k < nVel; k++ )
1350  {
1351  wk[k] = Array<OneD, NekDouble>(physTot,0.0);
1352  for ( j = 0; j < nVel; j++ )
1353  {
1354  Vmath::Vvtvp(nbc, CBnd[k*nVel+j], 1,
1355  fv[j], 1,
1356  wk[k], 1,
1357  wk[k], 1);
1358  }
1359  }
1360  for ( k = 0; k < nVel; k++ )
1361  {
1362  Vmath::Vcopy(nbc, wk[k], 1, fv[k], 1);
1363  }
1364 
1365  // Integrate to obtain force
1366  for ( j = 0; j < expdim; j++)
1367  {
1368  m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
1369  Integral(fp[j]);
1370  m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
1371  Integral(fv[j]);
1372  }
1373  }
1374  }
1375  else
1376  {
1377  cnt += BndExp[n]->GetExpSize();
1378  }
1379  }
1380  }
1381  }
1382 
1383  // Combine contributions from different processes
1384  // this is split between row and col comm because of
1385  // homostrips case, which only keeps its own strip
1386  for( i = 0; i < expdim; i++)
1387  {
1388  rowComm->AllReduce(m_Fpplane[i], LibUtilities::ReduceSum);
1389  colComm->AllReduce(m_Fpplane[i], LibUtilities::ReduceSum);
1390 
1391  rowComm->AllReduce(m_Fvplane[i], LibUtilities::ReduceSum);
1392  colComm->AllReduce(m_Fvplane[i], LibUtilities::ReduceSum);
1393  }
1394 
1395  // Project results to new directions
1396  for(plane = 0; plane < m_nPlanes; plane++)
1397  {
1398  Array< OneD, NekDouble> tmpP(expdim, 0.0);
1399  Array< OneD, NekDouble> tmpV(expdim, 0.0);
1400  for( i = 0; i < expdim; i++)
1401  {
1402  for( j = 0; j < expdim; j++ )
1403  {
1404  tmpP[i] += m_Fpplane[j][plane]*m_directions[i][j];
1405  tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
1406  }
1407  }
1408  // Copy result
1409  for( i = 0; i < expdim; i++)
1410  {
1411  m_Fpplane[i][plane] = tmpP[i];
1412  m_Fvplane[i][plane] = tmpV[i];
1413  }
1414  }
1415 
1416  // Sum viscous and pressure components
1417  for(plane = 0; plane < m_nPlanes; plane++)
1418  {
1419  for( i = 0; i < expdim; i++)
1420  {
1421  m_Ftplane[i][plane] = m_Fpplane[i][plane] + m_Fvplane[i][plane];
1422  }
1423  }
1424 
1425  // Put results back to wavespace, if necessary
1426  if( m_isHomogeneous1D)
1427  {
1428  for (i = 0; i < pFields.num_elements(); ++i)
1429  {
1430  pFields[i]->SetWaveSpace(true);
1431  pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
1432  pFields[i]->UpdatePhys());
1433  }
1434  }
1435 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, Array< OneD, NekDouble > > m_directions
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
int m_nPlanes
number of planes for homogeneous1D expansion
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:114
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
Array< OneD, Array< OneD, NekDouble > > m_Ftplane
Array< OneD, Array< OneD, NekDouble > > m_Fpplane
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, Array< OneD, NekDouble > > m_Fvplane
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
static FilterSharedPtr Nektar::SolverUtils::FilterAeroForces::create ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::map< std::string, std::string > &  pParams 
)
inlinestatic

Creates an instance of this class.

Definition at line 53 of file FilterAeroForces.h.

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

56  {
58  AllocateSharedPtr(pSession, pParams);
59  //p->InitObject();
60  return p;
61  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< Filter > FilterSharedPtr
A shared pointer to a Driver object.
Definition: Filter.h:50
void Nektar::SolverUtils::FilterAeroForces::GetForces ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
Array< OneD, NekDouble > &  Aeroforces,
const NekDouble time 
)

This function outputs the force on all planes of the current process, in the format required by ForcingMovingBody

Definition at line 515 of file FilterAeroForces.cpp.

References CalculateForces(), m_Ftplane, m_lastTime, and m_outputAllPlanes.

519 {
520  // Calculate forces if the result we have is not up-to-date
521  if(time > m_lastTime)
522  {
523  CalculateForces(pFields, time);
524  }
525  // Get information to write result
526  Array<OneD, unsigned int> ZIDs = pFields[0]->GetZIDs();
527  int local_planes = ZIDs.num_elements();
528  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
529 
530  // Copy results to Aeroforces
531  if (m_outputAllPlanes)
532  {
533  for(int plane = 0 ; plane < local_planes; plane++)
534  {
535  for(int dir =0; dir < expdim; dir++)
536  {
537  Aeroforces[plane + dir*local_planes] =
538  m_Ftplane[dir][ZIDs[plane]];
539  }
540  }
541  }
542  else
543  {
544  for(int plane = 0 ; plane < local_planes; plane++)
545  {
546  for(int dir =0; dir < expdim; dir++)
547  {
548  Aeroforces[plane + dir*local_planes] =
549  m_Ftplane[dir][0];
550  }
551  }
552  }
553 }
void CalculateForces(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Array< OneD, Array< OneD, NekDouble > > m_Ftplane
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...
void Nektar::SolverUtils::FilterAeroForces::v_Finalise ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
protectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 492 of file FilterAeroForces.cpp.

References m_outputStream.

495 {
496  if (pFields[0]->GetComm()->GetRank() == 0)
497  {
498  m_outputStream.close();
499  }
500 }
void Nektar::SolverUtils::FilterAeroForces::v_Initialise ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
protectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 201 of file FilterAeroForces.cpp.

References ASSERTL0, Nektar::StdRegions::find(), Nektar::ParseUtils::GenerateSeqVector(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::GlobalMapping::Mapping::Load(), m_BCtoElmtID, m_BCtoTraceID, m_boundaryRegionIsInList, m_boundaryRegionsIdList, m_BoundaryString, m_directions, m_index, m_isHomogeneous1D, m_lastTime, m_mapping, m_nPlanes, m_outputAllPlanes, m_outputFile, m_outputStream, m_planesID, Nektar::SolverUtils::Filter::m_session, and v_Update().

204 {
205  // Load mapping
207 
208  // Parse the boundary regions into a list.
209  std::string::size_type FirstInd =
210  m_BoundaryString.find_first_of('[') + 1;
211  std::string::size_type LastInd =
212  m_BoundaryString.find_last_of(']') - 1;
213 
214  ASSERTL0(FirstInd <= LastInd,
215  (std::string("Error reading boundary region definition:") +
216  m_BoundaryString).c_str());
217 
218  std::string IndString =
219  m_BoundaryString.substr(FirstInd, LastInd - FirstInd + 1);
220  bool parseGood = ParseUtils::GenerateSeqVector(IndString.c_str(),
222  ASSERTL0(parseGood && !m_boundaryRegionsIdList.empty(),
223  (std::string("Unable to read boundary regions index "
224  "range for FilterAeroForces: ") + IndString).c_str());
225 
226  // determine what boundary regions need to be considered
227  int cnt;
228  unsigned int numBoundaryRegions =
229  pFields[0]->GetBndConditions().num_elements();
231  numBoundaryRegions, 0);
232 
233  SpatialDomains::BoundaryConditions bcs(m_session,
234  pFields[0]->GetGraph());
236  bcs.GetBoundaryRegions();
237  SpatialDomains::BoundaryRegionCollection::const_iterator it;
238 
239  for (cnt = 0, it = bregions.begin(); it != bregions.end();
240  ++it, cnt++)
241  {
242  if ( std::find(m_boundaryRegionsIdList.begin(),
243  m_boundaryRegionsIdList.end(), it->first) !=
245  {
246  m_boundaryRegionIsInList[cnt] = 1;
247  }
248  }
249 
250  // Create map for element and edge/face of each boundary expansion
252  {
253  pFields[0]->GetPlane(0)->GetBoundaryToElmtMap
255  }
256  else
257  {
258  pFields[0]->GetBoundaryToElmtMap(m_BCtoElmtID,m_BCtoTraceID);
259  }
260 
261  // Define number of planes to calculate the forces
262  // in the Homogeneous direction ( if m_outputAllPlanes is false,
263  // consider only first plane in wave space)
264  // If flow has no Homogeneous direction, use 1 to make code general
265  if(m_isHomogeneous1D &&(m_outputAllPlanes || m_mapping->IsDefined()))
266  {
267  m_nPlanes = pFields[0]->GetHomogeneousBasis()->
268  GetZ().num_elements();
269  }
270  else
271  {
272  m_nPlanes = 1;
273  }
274 
275  // Create map for Planes ID for Homogeneous direction
276  // If flow has no Homogeneous direction, create trivial map
277  int j;
278  m_planesID = Array<OneD, int> (m_nPlanes,-1);
280  {
281  Array<OneD, const unsigned int> IDs = pFields[0]->GetZIDs();
282  //Loop through all planes
283  for(cnt = 0; cnt < m_nPlanes; cnt++)
284  {
285  for(j = 0; j < IDs.num_elements(); ++j)
286  {
287  //Look for current plane ID in this process
288  if(IDs[j] == cnt)
289  {
290  break;
291  }
292  }
293  // Assign ID to planesID
294  // If plane is not found in this process, leave it with -1
295  if(j != IDs.num_elements())
296  {
297  m_planesID[cnt] = j;
298  }
299  }
300  }
301  else
302  {
303  m_planesID[0] = 0;
304  }
305 
306  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
307 
308  // Write header
309  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
310  if (vComm->GetRank() == 0)
311  {
312  // Open output stream
313  m_outputStream.open(m_outputFile.c_str());
314  m_outputStream << "# Forces acting on bodies" << endl;
315  for( int i = 0; i < expdim; i++ )
316  {
317  m_outputStream << "#" << " Direction" << i+1 << " = (";
318  m_outputStream.width(8);
319  m_outputStream << setprecision(4) << m_directions[i][0];
320  m_outputStream.width(8);
321  m_outputStream << setprecision(4) << m_directions[i][1];
322  m_outputStream.width(8);
323  m_outputStream << setprecision(4) << m_directions[i][2];
324  m_outputStream << ")" << endl;
325  }
326  m_outputStream << "# Boundary regions: " << IndString.c_str() << endl;
327  m_outputStream << "#";
328  m_outputStream.width(7);
329  m_outputStream << "Time";
330  for( int i = 1; i <= expdim; i++ )
331  {
332  m_outputStream.width(8);
333  m_outputStream << "F" << i << "-press";
334  m_outputStream.width(9);
335  m_outputStream << "F" << i << "-visc";
336  m_outputStream.width(8);
337  m_outputStream << "F" << i << "-total";
338  }
339  if( m_outputAllPlanes )
340  {
341  m_outputStream.width(10);
342  m_outputStream << "Plane";
343  }
344  if (m_session->DefinesSolverInfo("HomoStrip"))
345  {
347  "Output forces on all planes not compatible with HomoStrips");
348  m_outputStream.width(10);
349  m_outputStream << "Strip";
350  }
351 
352  m_outputStream << endl;
353  }
354 
355  m_lastTime = -1;
356  m_index = 0;
357  v_Update(pFields, time);
358 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
virtual void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
std::vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, Array< OneD, NekDouble > > m_directions
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
int m_nPlanes
number of planes for homogeneous1D expansion
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:266
bool Nektar::SolverUtils::FilterAeroForces::v_IsTimeDependent ( )
protectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 506 of file FilterAeroForces.cpp.

507 {
508  return true;
509 }
void Nektar::SolverUtils::FilterAeroForces::v_Update ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
protectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 363 of file FilterAeroForces.cpp.

References CalculateForces(), m_Fpplane, m_Ftplane, m_Fvplane, m_index, m_nPlanes, m_outputAllPlanes, m_outputFrequency, m_outputStream, Nektar::SolverUtils::Filter::m_session, m_startTime, and Vmath::Vsum().

Referenced by v_Initialise().

366 {
367  // Only output every m_outputFrequency.
368  if ((m_index++) % m_outputFrequency || (time < m_startTime))
369  {
370  return;
371  }
372  // Calculate the forces
373  CalculateForces(pFields, time);
374 
375  // Calculate forces including all planes
376  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
377  Array<OneD, NekDouble> Fp(expdim,0.0);
378  Array<OneD, NekDouble> Fv(expdim,0.0);
379  Array<OneD, NekDouble> Ft(expdim,0.0);
380  for( int i = 0; i < expdim; i++)
381  {
382  Fp[i] = Vmath::Vsum(m_nPlanes, m_Fpplane[i], 1) / m_nPlanes;
383  Fv[i] = Vmath::Vsum(m_nPlanes, m_Fvplane[i], 1) / m_nPlanes;
384  Ft[i] = Fp[i] + Fv[i];
385  }
386 
387  // Communicators to exchange results
388  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
389 
390  //Write Results
391  if (vComm->GetRank() == 0)
392  {
393  // Write result in each plane
394  if( m_outputAllPlanes)
395  {
396  for( int plane = 0; plane < m_nPlanes; plane++)
397  {
398  // Write time
399  m_outputStream.width(8);
400  m_outputStream << setprecision(6) << time;
401  // Write forces
402  for( int i = 0; i < expdim; i++ )
403  {
404  m_outputStream.width(15);
405  m_outputStream << setprecision(8)
406  << m_Fpplane[i][plane];
407  m_outputStream.width(15);
408  m_outputStream << setprecision(8)
409  << m_Fvplane[i][plane];
410  m_outputStream.width(15);
411  m_outputStream << setprecision(8)
412  << m_Ftplane[i][plane];
413  }
414  m_outputStream.width(10);
415  m_outputStream << plane;
416  m_outputStream << endl;
417  }
418  }
419  // Output average (or total) force
420  m_outputStream.width(8);
421  m_outputStream << setprecision(6) << time;
422  for( int i = 0; i < expdim; i++)
423  {
424  m_outputStream.width(15);
425  m_outputStream << setprecision(8) << Fp[i];
426  m_outputStream.width(15);
427  m_outputStream << setprecision(8) << Fv[i];
428  m_outputStream.width(15);
429  m_outputStream << setprecision(8) << Ft[i];
430  }
431  if( m_outputAllPlanes)
432  {
433  m_outputStream.width(10);
434  m_outputStream << "average";
435  }
436 
437  if( m_session->DefinesSolverInfo("HomoStrip"))
438  {
439  // The result we already wrote is for strip 0
440  m_outputStream.width(10);
441  m_outputStream << 0;
442  m_outputStream << endl;
443 
444  // Now get result from other strips
445  int nstrips;
446  m_session->LoadParameter("Strip_Z", nstrips);
447  for(int i = 1; i<nstrips; i++)
448  {
449  vComm->GetColumnComm()->Recv(i, Fp);
450  vComm->GetColumnComm()->Recv(i, Fv);
451  vComm->GetColumnComm()->Recv(i, Ft);
452 
453  m_outputStream.width(8);
454  m_outputStream << setprecision(6) << time;
455  for( int j = 0; j < expdim; j++)
456  {
457  m_outputStream.width(15);
458  m_outputStream << setprecision(8) << Fp[j];
459  m_outputStream.width(15);
460  m_outputStream << setprecision(8) << Fv[j];
461  m_outputStream.width(15);
462  m_outputStream << setprecision(8) << Ft[j];
463  }
464  m_outputStream.width(10);
465  m_outputStream << i;
466  m_outputStream << endl;
467  }
468  }
469  else
470  {
471  m_outputStream << endl;
472  }
473  }
474  else
475  {
476  // In the homostrips case, we need to send result to root
477  if (m_session->DefinesSolverInfo("HomoStrip") &&
478  (vComm->GetRowComm()->GetRank() == 0) )
479  {
480  vComm->GetColumnComm()->Send(0, Fp);
481  vComm->GetColumnComm()->Send(0, Fv);
482  vComm->GetColumnComm()->Send(0, Ft);
483  }
484  }
485 
486 }
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
void CalculateForces(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:84
int m_nPlanes
number of planes for homogeneous1D expansion
Array< OneD, Array< OneD, NekDouble > > m_Ftplane
Array< OneD, Array< OneD, NekDouble > > m_Fpplane
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
Array< OneD, Array< OneD, NekDouble > > m_Fvplane
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...

Friends And Related Function Documentation

friend class MemoryManager< FilterAeroForces >
friend

Definition at line 50 of file FilterAeroForces.h.

Member Data Documentation

std::string Nektar::SolverUtils::FilterAeroForces::className
static
Initial value:

Name of the class.

Definition at line 64 of file FilterAeroForces.h.

Array<OneD, int> Nektar::SolverUtils::FilterAeroForces::m_BCtoElmtID
private

Definition at line 105 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Initialise().

Array<OneD, int> Nektar::SolverUtils::FilterAeroForces::m_BCtoTraceID
private

Definition at line 106 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Initialise().

std::vector<bool> Nektar::SolverUtils::FilterAeroForces::m_boundaryRegionIsInList
private

Determines if a given Boundary Region is in m_boundaryRegionsIdList.

Definition at line 94 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Initialise().

std::vector<unsigned int> Nektar::SolverUtils::FilterAeroForces::m_boundaryRegionsIdList
private

ID's of boundary regions where we want the forces.

Definition at line 91 of file FilterAeroForces.h.

Referenced by v_Initialise().

std::string Nektar::SolverUtils::FilterAeroForces::m_BoundaryString
private

Definition at line 104 of file FilterAeroForces.h.

Referenced by FilterAeroForces(), and v_Initialise().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::FilterAeroForces::m_directions
private
Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::FilterAeroForces::m_Fpplane
private

Definition at line 115 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Update().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::FilterAeroForces::m_Ftplane
private

Definition at line 117 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), GetForces(), and v_Update().

Array<OneD, Array<OneD, NekDouble> > Nektar::SolverUtils::FilterAeroForces::m_Fvplane
private

Definition at line 116 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Update().

LibUtilities::BasisSharedPtr Nektar::SolverUtils::FilterAeroForces::m_homogeneousBasis
private

Definition at line 103 of file FilterAeroForces.h.

unsigned int Nektar::SolverUtils::FilterAeroForces::m_index
private

Definition at line 95 of file FilterAeroForces.h.

Referenced by v_Initialise(), and v_Update().

bool Nektar::SolverUtils::FilterAeroForces::m_isHomogeneous1D
private
NekDouble Nektar::SolverUtils::FilterAeroForces::m_lastTime
private

Definition at line 119 of file FilterAeroForces.h.

Referenced by CalculateForces(), GetForces(), and v_Initialise().

GlobalMapping::MappingSharedPtr Nektar::SolverUtils::FilterAeroForces::m_mapping
private

Definition at line 120 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Initialise().

int Nektar::SolverUtils::FilterAeroForces::m_nPlanes
private

number of planes for homogeneous1D expansion

Definition at line 108 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), v_Initialise(), and v_Update().

bool Nektar::SolverUtils::FilterAeroForces::m_outputAllPlanes
private

if using a homogeneous1D expansion, determine if should output all planes or just the average

Definition at line 99 of file FilterAeroForces.h.

Referenced by CalculateForces(), FilterAeroForces(), GetForces(), v_Initialise(), and v_Update().

std::string Nektar::SolverUtils::FilterAeroForces::m_outputFile
private

Definition at line 101 of file FilterAeroForces.h.

Referenced by FilterAeroForces(), and v_Initialise().

unsigned int Nektar::SolverUtils::FilterAeroForces::m_outputFrequency
private

Definition at line 96 of file FilterAeroForces.h.

Referenced by FilterAeroForces(), and v_Update().

std::ofstream Nektar::SolverUtils::FilterAeroForces::m_outputStream
private

Definition at line 102 of file FilterAeroForces.h.

Referenced by v_Finalise(), v_Initialise(), and v_Update().

Array<OneD, int> Nektar::SolverUtils::FilterAeroForces::m_planesID
private

Definition at line 109 of file FilterAeroForces.h.

Referenced by CalculateForces(), CalculateForcesMapping(), and v_Initialise().

NekDouble Nektar::SolverUtils::FilterAeroForces::m_startTime
private

Definition at line 111 of file FilterAeroForces.h.

Referenced by FilterAeroForces(), and v_Update().