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