Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
AddExprToField.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
Include dependency graph for AddExprToField.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 15 of file AddExprToField.cpp.

References ASSERTL0, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::Import(), and Nektar::LibUtilities::Write().

16 {
17  int i,j;
18 
19  if(argc != 3)
20  {
21  fprintf(stderr,"Usage: AddExprToField meshfile fieldfile\n");
22  exit(1);
23  }
24 
25  LibUtilities::SessionReaderSharedPtr vSession
26  = LibUtilities::SessionReader::CreateInstance(argc, argv);
27 
28 
29  //----------------------------------------------
30  // Read in mesh from input file
31  string meshfile(argv[argc-2]);
32  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
33  //----------------------------------------------
34 
35  //----------------------------------------------
36  // Import field file.
37  string fieldfile(argv[argc-1]);
38  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
39  vector<vector<NekDouble> > fielddata;
40  LibUtilities::Import(fieldfile,fielddef,fielddata);
41  bool useFFT = false;
42  bool deal = false;
43  //----------------------------------------------
44 
45  //----------------------------------------------
46  // Define Expansion
47  int expdim = graphShPt->GetMeshDimension();
48  int nfields = fielddef[0]->m_fields.size();
49  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
50 
51  switch(expdim)
52  {
53  case 1:
54  {
55  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 2,"Quasi-3D approach is only set up for 1 or 2 homogeneous directions");
56 
57  if(fielddef[0]->m_numHomogeneousDir == 1)
58  {
59  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
60 
61  // Define Homogeneous expansion
62  //int nplanes = fielddef[0]->m_numModes[1];
63  int nplanes;
64  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[1]);
65 
66  // choose points to be at evenly spaced points at
67  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
68  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[1],nplanes,Pkey);
69  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
70 
71  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,ly,useFFT,deal,graphShPt);
72  Exp[0] = Exp2DH1;
73 
74  for(i = 1; i < nfields; ++i)
75  {
76  Exp[i] = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>::AllocateSharedPtr(*Exp2DH1);
77  }
78  }
79  else if(fielddef[0]->m_numHomogeneousDir == 2)
80  {
81  MultiRegions::ExpList3DHomogeneous2DSharedPtr Exp3DH2;
82 
83  // Define Homogeneous expansion
84  //int nylines = fielddef[0]->m_numModes[1];
85  //int nzlines = fielddef[0]->m_numModes[2];
86 
87  int nylines;
88  int nzlines;
89  vSession->LoadParameter("HomModesY",nylines,fielddef[0]->m_numModes[1]);
90  vSession->LoadParameter("HomModesZ",nzlines,fielddef[0]->m_numModes[2]);
91 
92  // choose points to be at evenly spaced points at
93  const LibUtilities::PointsKey PkeyY(nylines+1,LibUtilities::ePolyEvenlySpaced);
94  const LibUtilities::BasisKey BkeyY(fielddef[0]->m_basis[1],nylines,PkeyY);
95 
96  const LibUtilities::PointsKey PkeyZ(nzlines+1,LibUtilities::ePolyEvenlySpaced);
97  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
98 
99  NekDouble ly = fielddef[0]->m_homogeneousLengths[0];
100  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
101 
102  Exp3DH2 = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(vSession,BkeyY,BkeyZ,ly,lz,useFFT,deal,graphShPt);
103  Exp[0] = Exp3DH2;
104 
105  for(i = 1; i < nfields; ++i)
106  {
107  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous2D>::AllocateSharedPtr(*Exp3DH2);
108  }
109  }
110  else
111  {
112  MultiRegions::ExpList1DSharedPtr Exp1D;
113  Exp1D = MemoryManager<MultiRegions::ExpList1D>
114  ::AllocateSharedPtr(vSession,graphShPt);
115  Exp[0] = Exp1D;
116  for(i = 1; i < nfields ; ++i)
117  {
118  Exp[i] = MemoryManager<MultiRegions::ExpList1D>
119  ::AllocateSharedPtr(*Exp1D);
120  }
121  }
122  }
123  break;
124  case 2:
125  {
126  ASSERTL0(fielddef[0]->m_numHomogeneousDir <= 1,"NumHomogeneousDir is only set up for 1");
127 
128  if(fielddef[0]->m_numHomogeneousDir == 1)
129  {
130  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
131 
132  // Define Homogeneous expansion
133  //int nplanes = fielddef[0]->m_numModes[2];
134 
135  int nplanes;
136  vSession->LoadParameter("HomModesZ",nplanes,fielddef[0]->m_numModes[2]);
137 
138  // choose points to be at evenly spaced points at
139  // nplanes + 1 points
140  const LibUtilities::PointsKey Pkey(nplanes+1,LibUtilities::ePolyEvenlySpaced);
141  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
142  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
143 
144  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,useFFT,deal,graphShPt);
145  Exp[0] = Exp3DH1;
146 
147  for(i = 1; i < nfields; ++i)
148  {
149  Exp[i] = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(*Exp3DH1);
150  }
151  }
152  else
153  {
154  MultiRegions::ExpList2DSharedPtr Exp2D;
155  Exp2D = MemoryManager<MultiRegions::ExpList2D>
156  ::AllocateSharedPtr(vSession,graphShPt);
157  Exp[0] = Exp2D;
158 
159  for(i = 1; i < nfields; ++i)
160  {
161  Exp[i] = MemoryManager<MultiRegions::ExpList2D>
162  ::AllocateSharedPtr(*Exp2D);
163  }
164  }
165  }
166  break;
167  case 3:
168  {
169  MultiRegions::ExpList3DSharedPtr Exp3D;
170  Exp3D = MemoryManager<MultiRegions::ExpList3D>
171  ::AllocateSharedPtr(vSession,graphShPt);
172  Exp[0] = Exp3D;
173 
174  for(i = 1; i < nfields; ++i)
175  {
176  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
177  ::AllocateSharedPtr(*Exp3D);
178  }
179  }
180  break;
181  default:
182  ASSERTL0(false,"Expansion dimension not recognised");
183  break;
184  }
185  //----------------------------------------------
186 
187  //----------------------------------------------
188  // Copy data from field file
189  for(j = 0; j < nfields; ++j)
190  {
191  for(int i = 0; i < fielddata.size(); ++i)
192  {
193  Exp[j]->ExtractDataToCoeffs(fielddef [i],
194  fielddata[i],
195  fielddef [i]->m_fields[j],
196  Exp[j]->UpdateCoeffs());
197  }
198  Exp[j]->BwdTrans_IterPerExp(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
199  }
200  //----------------------------------------------
201 
202  //----------------------------------------------
203  // Add expression to field
204  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
205  int nq = Exp[0]->GetNpoints();
206  Array<OneD, NekDouble> x(nq);
207  Array<OneD, NekDouble> y(nq);
208  Exp[0]->GetCoords(x,y);
209  //Array<OneD, NekDouble> tmp(nq, 1000);
210  NekDouble pi = 3.14159265;
211  NekDouble tmp;
212  //the Exp.num_elements()==1!!!
213  ASSERTL0(Exp.num_elements()==1, "the field is not a streak");
214 cout<<"before Exp[0][1]="<<Exp[0]->GetPhys()[9]<<endl;
215  for (int i = 0; i < nq; ++i)
216  {
217  //sin(pi*y)*cos(x)
218  tmp = 0.01*sin(pi*y[i])*cos(x[i]);
219 //cout<<"add[0][1]="<<tmp<<endl;
220  //Vmath::Vadd(nq, Exp[0]->GetPhys(),1,tmp,1,Exp[0]->UpdatePhys(),1);
221  Exp[0]->UpdatePhys()[i] = Exp[0]->GetPhys()[i] +tmp;
222  }
223  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
224  //-----------------------------------------------
225  // Write solution to file with additional computed fields
226  string fldfilename(argv[2]);
227  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
228  string endfile("_add.fld");
229  out += endfile;
230  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
231  = Exp[0]->GetFieldDefinitions();
232  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
233  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(Exp.num_elements());
234 
235  for(int j = 0; j < nfields ; ++j)
236  {
237  Exp[j]->FwdTrans_IterPerExp(Exp[j]->GetPhys(),Exp[j]->UpdateCoeffs());
238  cout<<" Exp[0][0]="<<Exp[0]->GetPhys()[9]<<endl;
239  fieldcoeffs[j] = Exp[j]->UpdateCoeffs();
240  for(int i = 0; i < FieldDef.size(); ++i)
241  {
242  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
243  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
244  }
245  }
246  LibUtilities::Write(out, FieldDef, FieldData);
247  //-----------------------------------------------
248 
249  return 0;
250 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::MultiRegions::ExpList3DHomogeneous2DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous2D > ExpList3DHomogeneous2DSharedPtr
Shared pointer to an ExpList3DHomogeneous2D object.
Definition: ExpList3DHomogeneous2D.h:50
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
Definition: FieldIO.cpp:115
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::MultiRegions::ExpList2DSharedPtr
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
Nektar::MultiRegions::ExpList2DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
Definition: ExpList2DHomogeneous1D.h:50
Nektar::MultiRegions::ExpList3DSharedPtr
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:114
Nektar::MultiRegions::ExpList3DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
Definition: ExpList3DHomogeneous1D.h:50
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:81
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Generated on Fri Sep 9 2016 12:03:21 for Nektar++ by   doxygen 1.8.8