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  int i = 0;
48  int j = 0;
49  int k = 0;
50  int n = 0;
51 
52  for(i = 0; i < fielddata2.size(); ++i)
53  {
54  ASSERTL0(fielddef2[i]->m_numHomogeneousDir == 1,"Expected second fld to have one homogeneous direction");
55  ASSERTL0(fielddef2[i]->m_numModes[2] == 2,"Expected Fourier field to have 2 modes");
56 
57  int datalen1 = fielddata1[i].size()/fielddef1[i]->m_fields.size();
58  int datalen2 = fielddata2[i].size()/fielddef2[i]->m_fields.size();
59 
60  ASSERTL0(datalen1*2 == datalen2,"Data per fields is note compatible");
61 
62  // Determine the number of coefficients per element
63  int ncoeffs = 0;
64  switch(fielddef2[i]->m_shapeType)
65  {
67  ncoeffs = LibUtilities::StdTriData::getNumberOfCoefficients(fielddef2[i]->m_numModes[0], fielddef2[i]->m_numModes[1]);
68  break;
70  ncoeffs = fielddef2[i]->m_numModes[0]*fielddef2[i]->m_numModes[1];
71  break;
72  default:
73  ASSERTL0(false,"Shape not recognised");
74  break;
75  }
76 
77  // array for zero packing
78  Array<OneD,NekDouble> Zero(ncoeffs,0.0);
79 
80  // scale first and second fields
81  Vmath::Smul(fielddata1[i].size(), scal1, &fielddata1[i][0], 1,
82  &fielddata1[i][0], 1);
83  Vmath::Smul(fielddata2[i].size(), scal2, &fielddata2[i][0], 1,
84  &fielddata2[i][0], 1);
85 
87 
88  vector<NekDouble> newdata;
89  vec_iter = fielddata2[i].begin();
90 
91  for(k = 0; k < fielddef2[i]->m_fields.size(); ++k)
92  {
93  // get location of 2D field information in order of field2 ordering
94  int offset = 0;
95  for(j = 0; j < fielddef1[i]->m_fields.size(); ++j)
96  {
97  if(fielddef1[i]->m_fields[j] == fielddef2[i]->m_fields[k])
98  {
99  break;
100  }
101  offset += datalen1;
102  }
103 
104  if(j != fielddef1[i]->m_fields.size())
105  {
106  for(n = 0; n < fielddef2[i]->m_elementIDs.size(); ++n)
107  {
108  // Real zero component
109  newdata.insert(newdata.end(),
110  &(fielddata1[i][offset+n*ncoeffs]),
111  &(fielddata1[i][offset+n*ncoeffs])
112  + ncoeffs);
113 
114  // Imaginary zero component;
115  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
116 
117  // Put orginal mode in here.
118  newdata.insert(newdata.end(),vec_iter, vec_iter+2*ncoeffs);
119  vec_iter += 2*ncoeffs;
120  }
121  }
122  else
123  {
124 
125  for(n = 0; n < fielddef2[i]->m_elementIDs.size(); ++n)
126  {
127  // Real & Imag zero component
128  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
129  newdata.insert(newdata.end(),&Zero[0],&Zero[0] + ncoeffs);
130 
131  // Put orginal mode in here.
132  newdata.insert(newdata.end(),vec_iter, vec_iter+2*ncoeffs);
133  vec_iter += 2*ncoeffs;
134  }
135  }
136  }
137  combineddata.push_back(newdata);
138  fielddef2[i]->m_numModes[2] += 2;
139  fielddef2[i]->m_homogeneousZIDs.push_back(2);
140  fielddef2[i]->m_homogeneousZIDs.push_back(3);
141 
142  // check to see if any field in fielddef1[i]->m_fields is
143  // not defined in fielddef2[i]->m_fields
144  for(k = 0; k < fielddef1[i]->m_fields.size(); ++k)
145  {
146  for(j = 0; j < fielddef2[i]->m_fields.size(); ++j)
147  {
148  if(fielddef1[i]->m_fields[k] == fielddef2[i]->m_fields[j])
149  {
150  break;
151  }
152  }
153 
154  if(j == fielddef2[i]->m_fields.size())
155  {
156  cout << "Warning: Field \'" << fielddef1[i]->m_fields[k]
157  << "\' was not included in output " << endl;
158  }
159 
160  }
161 
162  }
163  //----------------------------------------------
164 
165  //-----------------------------------------------
166  // Write out datafile.
167  LibUtilities::Write(argv[argc-1], fielddef2, combineddata);
168  //-----------------------------------------------
169 
170  return 0;
171 }
172 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
double NekDouble
int main(int argc, char *argv[])
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Base class for a spectral/hp element mesh.
Definition: MeshGraph.h:186
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359