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  bool useFFT = false;
178  bool deal = false;
179  ///Parameter for homogeneous expansions
180  enum HomogeneousType
181  {
182  eHomogeneous1D,
183  eHomogeneous2D,
184  eHomogeneous3D,
185  eNotHomogeneous
186  };
187 
188  enum HomogeneousType HomogeneousType = eNotHomogeneous;
189 
190  if(session->DefinesSolverInfo("HOMOGENEOUS"))
191  {
192  std::string HomoStr = session->GetSolverInfo("HOMOGENEOUS");
193  //m_spacedim = 3;
194 
195  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
196  (HomoStr == "1D")||(HomoStr == "Homo1D"))
197  {
198  HomogeneousType = eHomogeneous1D;
199  npointsZ = session->GetParameter("HomModesZ");
200  LhomZ = session->GetParameter("LZ");
201  }
202 
203  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
204  (HomoStr == "2D")||(HomoStr == "Homo2D"))
205  {
206  HomogeneousType = eHomogeneous2D;
207  npointsY = session->GetParameter("HomModesY");
208  LhomY = session->GetParameter("LY");
209  npointsZ = session->GetParameter("HomModesZ");
210  LhomZ = session->GetParameter("LZ");
211  }
212 
213  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
214  (HomoStr == "3D")||(HomoStr == "Homo3D"))
215  {
216  HomogeneousType = eHomogeneous3D;
217  //npointsX = session->GetParameter("HomModesX");
218  //LhomX = session->GetParameter("LX");
219  npointsY = session->GetParameter("HomModesY");
220  LhomY = session->GetParameter("LY");
221  npointsZ = session->GetParameter("HomModesZ");
222  LhomZ = session->GetParameter("LZ");
223  }
224 
225  if(session->DefinesSolverInfo("USEFFT"))
226  {
227  useFFT = true;
228  }
229  }
230 
231  int i;
232  int expdim = mesh->GetMeshDimension();
233  //Exp= Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
234  // I can always have 3 variables in a 2D mesh (oech vel component i a function which can depend on 1-3 var)
235  // Continuous Galerkin projection
236 
237  switch(expdim)
238  {
239  case 1:
240  {
241  if(HomogeneousType == eHomogeneous2D)
242  {
243  const LibUtilities::PointsKey PkeyY(npointsY,LibUtilities::eFourierEvenlySpaced);
244  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,npointsY,PkeyY);
245  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
246  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
247 
248  for(i = 0 ; i < nvariables; i++)
249  {
250  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
251  ::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
252  }
253  }
254  else
255  {
256  for(i = 0 ; i < nvariables; i++)
257  {
258  Exp[i] = MemoryManager<MultiRegions::ContField1D>
259  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
260  }
261  }
262 
263  break;
264  }
265  case 2:
266  {
267  if(HomogeneousType == eHomogeneous1D)
268  {
269  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
270  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
271  for(i = 0 ; i < nvariables; i++)
272  {
273  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
274  ::AllocateSharedPtr(session,BkeyZ,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
275  }
276  }
277  else
278  {
279  i = 0;
280  MultiRegions::ContField2DSharedPtr firstfield;
281  firstfield = MemoryManager<MultiRegions::ContField2D>
282  ::AllocateSharedPtr(session,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
283 
284  Exp[0] = firstfield;
285  for(i = 1 ; i < nvariables; i++)
286  {
287  Exp[i] = MemoryManager<MultiRegions::ContField2D>
288  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
289  }
290  }
291 
292  break;
293  }
294  case 3:
295  {
296  if(HomogeneousType == eHomogeneous3D)
297  {
298  ASSERTL0(false,"3D fully periodic problems not implemented yet");
299  }
300  else
301  {
302  i = 0;
303  MultiRegions::ContField3DSharedPtr firstfield =
304  MemoryManager<MultiRegions::ContField3D>
305  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
306 
307  Exp[0] = firstfield;
308  for(i = 1 ; i < nvariables; i++)
309  {
310  Exp[i] = MemoryManager<MultiRegions::ContField3D>
311  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i));
312  }
313  }
314  break;
315  }
316  default:
317  ASSERTL0(false,"Expansion dimension not recognised");
318  break;
319  }
320  }
321 
Generated on Wed Mar 4 2015 16:33:17 for Nektar++ by   doxygen 1.8.1.2