Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FldAddFld.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 
4 #include <LibUtilities/BasicUtils/ErrorUtil.hpp> // for ASSERTL0
5 #include <SpatialDomains/MeshGraph.h> // for FieldDefinitions, etc
6 
7 using namespace Nektar;
8 
9 int main(int argc, char *argv[])
10 {
11  NekDouble scal1,scal2;
12 
13  if(argc != 6)
14  {
15  fprintf(stderr,"Usage: FldAddFld scal1 scal2 fieldfile1 fieldfile2 outfield\n"
16  "\t produces scal1*fieldfiel1 + scal2*fieldfile2 in outfield\n" );
17  exit(1);
18  }
19 
20  scal1 = boost::lexical_cast<double>(argv[argc-5]);
21  scal2 = boost::lexical_cast<double>(argv[argc-4]);
22 
23  //default meshgraph
25 
26  //----------------------------------------------
27  // Import fieldfile1.
28  string fieldfile1(argv[argc-3]);
29  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef1;
30  vector<vector<NekDouble> > fielddata1;
31  LibUtilities::Import(fieldfile1,fielddef1,fielddata1);
32  //----------------------------------------------
33 
34  //----------------------------------------------
35  // Import fieldfile2.
36  string fieldfile2(argv[argc-2]);
37  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef2;
38  vector<vector<NekDouble> > fielddata2;
39  LibUtilities::Import(fieldfile2,fielddef2,fielddata2);
40  //----------------------------------------------
41 
42 
43  ASSERTL0(fielddata1.size() == fielddata2.size(),"Inner has different size");
44 
45  //----------------------------------------------
46  // Add fielddata2 to fielddata1 using m_fields definition to align data.
47 
48 
49  for(int i = 0; i < fielddata1.size(); ++i)
50  {
51  int j;
52  int datalen1 = fielddata1[i].size()/fielddef1[i]->m_fields.size();
53  int datalen2 = fielddata2[i].size()/fielddef2[i]->m_fields.size();
54 
55  ASSERTL0(datalen1 == datalen2,"Data per field is of different length");
56 
57  for(int k = 0; k < fielddef1[i]->m_fields.size(); ++k)
58  {
59  int offset = 0;
60  for(j = 0; j < fielddef2[i]->m_fields.size(); ++j)
61  {
62  if(fielddef1[i]->m_fields[k] == fielddef2[i]->m_fields[j])
63  {
64  break;
65  }
66  offset += datalen1;
67  }
68 
69  if(j == fielddef2[i]->m_fields.size())
70  {
71  for(j = 0; j < datalen1; ++j)
72  {
73  fielddata1[i][datalen1*k+j] *= scal1;
74  }
75  }
76  else // add fields
77  {
78  for(j = 0; j < datalen1; ++j)
79  {
80  fielddata1[i][datalen1*k+j] *= scal1;
81  fielddata1[i][datalen1*k+j] += scal2*fielddata2[i][offset + j];
82  }
83  }
84 
85  }
86 
87  // now check to see if any field in fielddef2[i]->m_fields is
88  // not defined in fielddef1[i]->m_fields
89  for(int k = 0; k < fielddef2[i]->m_fields.size(); ++k)
90  {
91  for(j = 0; j < fielddef1[i]->m_fields.size(); ++j)
92  {
93  if(fielddef2[i]->m_fields[k] == fielddef1[i]->m_fields[j])
94  {
95  break;
96  }
97  }
98 
99  if(j == fielddef1[i]->m_fields.size())
100  {
101  for(j = 0; j < datalen2; ++j)
102  {
103  fielddata2[i][datalen2*k+j] *= scal2;
104  }
105 
106  // add this field to fielddata1
107  fielddef1[i]->m_fields.push_back(fielddef2[i]->m_fields[k]);
108  fielddata1[i].insert(fielddata1[i].end(),&(fielddata2[i][k*datalen2]),
109  &(fielddata2[i][k*datalen2])+datalen1);
110  }
111 
112  }
113 
114  }
115  //----------------------------------------------
116 
117  //-----------------------------------------------
118  // Write out datafile.
119  LibUtilities::Write(argv[argc-1], fielddef1, fielddata1);
120  //-----------------------------------------------
121 
122  return 0;
123 }
124