Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SplitFld.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <LibUtilities/BasicUtils/SessionReader.h>
4 #include <MultiRegions/ExpList.h>
5 #include <MultiRegions/ContField1D.h>
6 #include <MultiRegions/ContField2D.h>
7 #include <MultiRegions/ContField3D.h>
8 #include <MultiRegions/ContField3DHomogeneous1D.h>
9 #include <MultiRegions/ContField3DHomogeneous2D.h>
10 
11 
12 using namespace Nektar;
13 
14 int main(int argc, char *argv[])
15 {
16 
17  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
18  SpatialDomains::BoundaryConditionsSharedPtr& boundaryConditions,
19  LibUtilities::SessionReaderSharedPtr &session,
20  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables);
21  Array<OneD, int> GetReflectionIndex(Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,
22  int Ireg);
23 
24  if(argc != 3)
25  {
26  fprintf(stderr,"Usage: SplitFld meshfile fieldfile\n");
27  exit(1);
28  }
29 
30  LibUtilities::SessionReaderSharedPtr vSession
31  = LibUtilities::SessionReader::CreateInstance(argc, argv);
32 
33 
34  //----------------------------------------------
35 
36  // Read in mesh from input file
37  string meshfile(argv[argc-2]);
38  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
39  //----------------------------------------------
40 
41  // Also read and store the boundary conditions
42  SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
43  boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
44  ::AllocateSharedPtr(vSession, graphShPt);
45  SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
46  //----------------------------------------------
47 
48  //----------------------------------------------
49  // Import field file.
50  string fieldfile(argv[argc-1]);
51  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
52  vector<vector<NekDouble> > fielddata;
53  LibUtilities::Import(fieldfile,fielddef,fielddata);
54  //----------------------------------------------
55 
56  // Define Expansion
57  int nfields;
58  nfields = fielddef[0]->m_fields.size();
59  Array<OneD, MultiRegions::ExpListSharedPtr> Exp;
60  Exp = Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
61 
62  std::string solvtype = vSession->GetSolverInfo("SOLVERTYPE");
63  if(solvtype == "CoupledLinearisedNS" && vSession->DefinesSolverInfo("HOMOGENEOUS") )
64  {
65 
66  SetFields(graphShPt,boundaryConditions,vSession,Exp,nfields-1);
67  //decomment
68  //nfields = nfields-1;
69 //start
70  int lastfield = nfields-1;
71  cout<<"Set pressure: "<<lastfield<<endl;
72  int nplanes = fielddef[0]->m_numModes[2];
73  const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
74  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
75  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
76  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
77  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt,fielddef[0]->m_fields[0]);
78  Exp[lastfield] = Exp3DH1;
79 //end
80  }
81  else
82  {
83 
84  SetFields(graphShPt,boundaryConditions,vSession,Exp,nfields);
85  }
86  //----------------------------------------------
87 
88  //----------------------------------------------
89  // Copy data from field file
90  for(int j = 0; j < nfields; ++j)
91  {
92  for(int i = 0; i < fielddef.size(); ++i)
93  {
94  Exp[j]->ExtractDataToCoeffs(fielddef [i],
95  fielddata[i],
96  fielddef [i]->m_fields[j],
97  Exp[j]->UpdateCoeffs());
98  }
99  Exp[j]->BwdTrans_IterPerExp(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
100  }
101 
102  //----------------------------------------------
103 
104 
105  // Write solution to file with additional computed fields
106  string fldfilename(argv[2]);
107  string out = fldfilename.substr(0, fldfilename.find_last_of("."));
108  string endfile("split.fld");
109 
110 
111  //Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
112 
113  string outfile;
114  string var;
115 
116 
117  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(Exp.num_elements());
118  //NB in case of homo fields you CANNOT use the BwdTrans
119  //because the Im comp is set to 0
120  for(int i = 0; i < Exp.num_elements(); ++i)
121  {
122  fieldcoeffs[i] = Exp[i]->UpdateCoeffs();
123  }
124 
125 
126  // copy Data into FieldData and set variable
127  /*
128  for(int g=0; g<Exp[0]->GetPlane(1)->GetNcoeffs(); g++)
129  {
130  cout<<"g="<<g<<" coeff f0="<<Exp[lastfield]->GetPlane(0)->GetCoeff(g)<<" f1="<<Exp[lastfield]->GetPlane(1)->GetCoeff(g)<<endl;
131  }
132  */
133 
134  for(int j =0; j<nfields; j++)
135  {
136  outfile = out;
137  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
138  = Exp[j]->GetFieldDefinitions();
139  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
140  //fieldcoeffs[j] = fields[j]->UpdateCoeffs();
141  for(int i = 0; i < FieldDef.size(); ++i)
142  {
143  var = fielddef[i]->m_fields[j];
144  // Could do a search here to find correct variable
145  FieldDef[i]->m_fields.push_back(var);
146  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[j]);
147 
148  }
149  outfile += "_"+var+"_"+endfile;
150  LibUtilities::Write(outfile,FieldDef,FieldData);
151 
152  }
153 
154  //-----------------------------------------------
155 
156  return 0;
157 }
158 
159 
160 
161 
162  // Define Expansion
163  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
164  SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
165  LibUtilities::SessionReaderSharedPtr &session,
166  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables)
167  {
168  // Setting parameteres for homogenous problems
169  NekDouble LhomX; ///< physical length in X direction (if homogeneous)
170  NekDouble LhomY; ///< physical length in Y direction (if homogeneous)
171  NekDouble LhomZ; ///< physical length in Z direction (if homogeneous)
172 
173  bool DeclareCoeffPhysArrays = true;
174  int npointsX; ///< number of points in X direction (if homogeneous)
175  int npointsY; ///< number of points in Y direction (if homogeneous)
176  int npointsZ; ///< number of points in Z direction (if homogeneous)
177  int HomoDirec = 0;
178  bool useFFT = false;
179  bool deal = false;
180  ///Parameter for homogeneous expansions
181  enum HomogeneousType
182  {
183  eHomogeneous1D,
184  eHomogeneous2D,
185  eHomogeneous3D,
186  eNotHomogeneous
187  };
188 
189  enum HomogeneousType HomogeneousType = eNotHomogeneous;
190 
191  if(session->DefinesSolverInfo("HOMOGENEOUS"))
192  {
193  std::string HomoStr = session->GetSolverInfo("HOMOGENEOUS");
194  //m_spacedim = 3;
195 
196  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
197  (HomoStr == "1D")||(HomoStr == "Homo1D"))
198  {
199  HomogeneousType = eHomogeneous1D;
200  npointsZ = session->GetParameter("HomModesZ");
201  LhomZ = session->GetParameter("LZ");
202  HomoDirec = 1;
203  }
204 
205  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
206  (HomoStr == "2D")||(HomoStr == "Homo2D"))
207  {
208  HomogeneousType = eHomogeneous2D;
209  npointsY = session->GetParameter("HomModesY");
210  LhomY = session->GetParameter("LY");
211  npointsZ = session->GetParameter("HomModesZ");
212  LhomZ = session->GetParameter("LZ");
213  HomoDirec = 2;
214  }
215 
216  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
217  (HomoStr == "3D")||(HomoStr == "Homo3D"))
218  {
219  HomogeneousType = eHomogeneous3D;
220  npointsX = session->GetParameter("HomModesX");
221  LhomX = session->GetParameter("LX");
222  npointsY = session->GetParameter("HomModesY");
223  LhomY = session->GetParameter("LY");
224  npointsZ = session->GetParameter("HomModesZ");
225  LhomZ = session->GetParameter("LZ");
226  HomoDirec = 3;
227  }
228 
229  if(session->DefinesSolverInfo("USEFFT"))
230  {
231  useFFT = true;
232  }
233  }
234 
235  int i;
236  int expdim = mesh->GetMeshDimension();
237  //Exp= Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
238  // I can always have 3 variables in a 2D mesh (oech vel component i a function which can depend on 1-3 var)
239  // Continuous Galerkin projection
240 
241  switch(expdim)
242  {
243  case 1:
244  {
245  if(HomogeneousType == eHomogeneous2D)
246  {
247  const LibUtilities::PointsKey PkeyY(npointsY,LibUtilities::eFourierEvenlySpaced);
248  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,npointsY,PkeyY);
249  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
250  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
251 
252  for(i = 0 ; i < nvariables; i++)
253  {
254  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
255  ::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
256  }
257  }
258  else
259  {
260  for(i = 0 ; i < nvariables; i++)
261  {
262  Exp[i] = MemoryManager<MultiRegions::ContField1D>
263  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
264  }
265  }
266 
267  break;
268  }
269  case 2:
270  {
271  if(HomogeneousType == eHomogeneous1D)
272  {
273  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
274  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
275  for(i = 0 ; i < nvariables; i++)
276  {
277  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
278  ::AllocateSharedPtr(session,BkeyZ,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
279  }
280  }
281  else
282  {
283  i = 0;
284  MultiRegions::ContField2DSharedPtr firstfield;
285  firstfield = MemoryManager<MultiRegions::ContField2D>
286  ::AllocateSharedPtr(session,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
287 
288  Exp[0] = firstfield;
289  for(i = 1 ; i < nvariables; i++)
290  {
291  Exp[i] = MemoryManager<MultiRegions::ContField2D>
292  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
293  }
294  }
295 
296  break;
297  }
298  case 3:
299  {
300  if(HomogeneousType == eHomogeneous3D)
301  {
302  ASSERTL0(false,"3D fully periodic problems not implemented yet");
303  }
304  else
305  {
306  i = 0;
307  MultiRegions::ContField3DSharedPtr firstfield =
308  MemoryManager<MultiRegions::ContField3D>
309  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
310 
311  Exp[0] = firstfield;
312  for(i = 1 ; i < nvariables; i++)
313  {
314  Exp[i] = MemoryManager<MultiRegions::ContField3D>
315  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i));
316  }
317  }
318  break;
319  }
320  default:
321  ASSERTL0(false,"Expansion dimension not recognised");
322  break;
323  }
324  }
325 
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2