Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AddModeTo2DFld.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <SpatialDomains/MeshGraph.h> // for FieldDefinitions, etc
4 #include <StdRegions/StdTriExp.h>
5 
6 using namespace Nektar;
7 
8 int main(int argc, char *argv[])
9 {
10  NekDouble scal1,scal2;
11 
12  if(argc != 6)
13  {
14  fprintf(stderr,"Usage: AddModeTo2DFld scal1 scal2 2Dfieldfile1 fieldfile2 outfield\n"
15  "\t produces scal1*2Dfieldfiel1 + scal2*fieldfile2 in outfield\n" );
16  exit(1);
17  }
18 
19  scal1 = boost::lexical_cast<double>(argv[argc-5]);
20  scal2 = boost::lexical_cast<double>(argv[argc-4]);
21 
22  //default meshgraph
24 
25  //----------------------------------------------
26  // Import fieldfile1.
27  string fieldfile1(argv[argc-3]);
28  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef1;
29  vector<vector<NekDouble> > fielddata1;
30  LibUtilities::Import(fieldfile1,fielddef1,fielddata1);
31  //----------------------------------------------
32 
33  //----------------------------------------------
34  // Import fieldfile2.
35  string fieldfile2(argv[argc-2]);
36  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef2;
37  vector<vector<NekDouble> > fielddata2;
38  LibUtilities::Import(fieldfile2,fielddef2,fielddata2);
39  //----------------------------------------------
40 
41  vector<vector<NekDouble> > combineddata;
42 
43  ASSERTL0(fielddata1.size() == fielddata2.size(),"Inner has different size");
44  //----------------------------------------------
45  // Add fielddata2 to fielddata1 using m_fields definition to align data.
46 
47  for(int i = 0; i < fielddata2.size(); ++i)
48  {
49  ASSERTL0(fielddef2[i]->m_numHomogeneousDir == 1,"Expected second fld to have one homogeneous direction");
50  ASSERTL0(fielddef2[i]->m_numModes[2] == 2,"Expected Fourier field to have 2 modes");
51 
52  int j;
53  int datalen1 = fielddata1[i].size()/fielddef1[i]->m_fields.size();
54  int datalen2 = fielddata2[i].size()/fielddef2[i]->m_fields.size();
55 
56  ASSERTL0(datalen1*2 == datalen2,"Data per fields is note compatible");
57 
58  // Determine the number of coefficients per element
59  int ncoeffs = 0;
60  switch(fielddef2[i]->m_shapeType)
61  {
63  ncoeffs = LibUtilities::StdTriData::getNumberOfCoefficients(fielddef2[i]->m_numModes[0], fielddef2[i]->m_numModes[1]);
64  break;
66  ncoeffs = fielddef2[i]->m_numModes[0]*fielddef2[i]->m_numModes[1];
67  break;
68  default:
69  ASSERTL0(false,"Shape not recognised");
70  break;
71  }
72 
73  // array for zero packing
74  Array<OneD,NekDouble> Zero(ncoeffs,0.0);
75 
76 
77  // scale first field
78  for(j = 0; j < fielddata1[i].size(); ++j)
79  {
80  fielddata1[i][j] *= scal1;
81 
82  }
83 
84  // scale second field
85  for(j = 0; j < fielddata2[i].size(); ++j)
86  {
87  fielddata2[i][j] *= scal2;
88 
89  }
90 
92 
93  vector<NekDouble> newdata;
94  vec_iter = fielddata2[i].begin();
95 
96  for(int k = 0; k < fielddef2[i]->m_fields.size(); ++k)
97  {
98 
99  // get location of 2D field information in order of field2 ordering
100  int offset = 0;
101  for(j = 0; j < fielddef1[i]->m_fields.size(); ++j)
102  {
103  if(fielddef1[i]->m_fields[j] == fielddef2[i]->m_fields[k])
104  {
105  break;
106  }
107  offset += datalen1;
108  }
109 
110  if(j != fielddef1[i]->m_fields.size())
111  {
112  for(int n = 0; n < fielddef2[i]->m_elementIDs.size(); ++n)
113  {
114  // Real zero component
115  newdata.insert(newdata.end(),&(fielddata1[i][offset+n*ncoeffs]),
116 &(fielddata1[i][offset+n*ncoeffs]) + ncoeffs);
117  // Imaginary zero component;
118  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
119 
120  // Put orginal mode in here.
121  newdata.insert(newdata.end(),vec_iter, vec_iter+2*ncoeffs);
122  vec_iter += 2*ncoeffs;
123  }
124  }
125  else
126  {
127 
128  for(int n = 0; n < fielddef2[i]->m_elementIDs.size(); ++n)
129  {
130  // Real & Imag zero component
131  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
132  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
133 
134  // Put orginal mode in here.
135  newdata.insert(newdata.end(),vec_iter, vec_iter+2*ncoeffs);
136  vec_iter += 2*ncoeffs;
137  }
138  }
139  }
140  combineddata.push_back(newdata);
141  fielddef2[i]->m_numModes[2] += 2;
142  fielddef2[i]->m_homogeneousZIDs.push_back(2);
143  fielddef2[i]->m_homogeneousZIDs.push_back(3);
144 
145  // check to see if any field in fielddef1[i]->m_fields is
146  // not defined in fielddef2[i]->m_fields
147  for(int k = 0; k < fielddef1[i]->m_fields.size(); ++k)
148  {
149  for(j = 0; j < fielddef2[i]->m_fields.size(); ++j)
150  {
151  if(fielddef1[i]->m_fields[k] == fielddef2[i]->m_fields[j])
152  {
153  break;
154  }
155  }
156 
157  if(j == fielddef2[i]->m_fields.size())
158  {
159  cout << "Warning: Field \'" << fielddef1[i]->m_fields[k] << "\' was not included in output " << endl;
160  }
161 
162  }
163 
164  }
165  //----------------------------------------------
166 
167  //-----------------------------------------------
168  // Write out datafile.
169  LibUtilities::Write(argv[argc-1], fielddef2, combineddata);
170  //-----------------------------------------------
171 
172  return 0;
173 }
174