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

vector< unsigned int > m_boundaryRegionsIdList
 ID's of boundary regions where we want the forces. More...
 
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 58 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.

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

192 {
193 
194 }

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 556 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().

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

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

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

517 {
518  // Calculate forces if the result we have is not up-to-date
519  if(time > m_lastTime)
520  {
521  CalculateForces(pFields, time);
522  }
523  // Get information to write result
524  Array<OneD, unsigned int> ZIDs = pFields[0]->GetZIDs();
525  int local_planes = ZIDs.num_elements();
526  int expdim = pFields[0]->GetGraph()->GetMeshDimension();
527 
528  // Copy results to Aeroforces
529  if (m_outputAllPlanes)
530  {
531  for(int plane = 0 ; plane < local_planes; plane++)
532  {
533  for(int dir =0; dir < expdim; dir++)
534  {
535  Aeroforces[plane + dir*local_planes] =
536  m_Ftplane[dir][ZIDs[plane]];
537  }
538  }
539  }
540  else
541  {
542  for(int plane = 0 ; plane < local_planes; plane++)
543  {
544  for(int dir =0; dir < expdim; dir++)
545  {
546  Aeroforces[plane + dir*local_planes] =
547  m_Ftplane[dir][0];
548  }
549  }
550  }
551 }
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 490 of file FilterAeroForces.cpp.

References m_outputStream.

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

Implements Nektar::SolverUtils::Filter.

Definition at line 199 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().

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

Implements Nektar::SolverUtils::Filter.

Definition at line 504 of file FilterAeroForces.cpp.

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

Implements Nektar::SolverUtils::Filter.

Definition at line 361 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().

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

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().

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().