Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessQCriterion.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessQCriterion.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Computes Q Criterion field.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessQCriterion.h"
41 
42 #include <LibUtilities/BasicUtils/SharedArray.hpp>
43 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
44 
45 namespace Nektar
46 {
47  namespace Utilities
48  {
49  ModuleKey ProcessQCriterion::className =
50  GetModuleFactory().RegisterCreatorFunction(
51  ModuleKey(eProcessModule, "QCriterion"),
52  ProcessQCriterion::create, "Computes Q-Criterion.");
53 
54  ProcessQCriterion::ProcessQCriterion(FieldSharedPtr f) :
55  ProcessModule(f)
56  {
57  }
58 
59  ProcessQCriterion::~ProcessQCriterion()
60  {
61  }
62 
63  void ProcessQCriterion::Process(po::variables_map &vm)
64  {
65  if (m_f->m_verbose)
66  {
67  cout << "ProcessQCriterion: Calculating Q Criterion..." << endl;
68  }
69 
70  int i, j;
71  int expdim = m_f->m_graph->GetMeshDimension();
72  int spacedim = expdim;
73  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
74  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
75  {
76  spacedim = 3;
77  }
78  int nfields = m_f->m_fielddef[0]->m_fields.size();
79  if (spacedim == 1 || spacedim == 2)
80  {
81  cerr << "\n Error: ProcessQCriterion must be computed for a 3D"
82  " (or quasi-3D) case. \n" << endl;
83  }
84 
85  //For calculating Q-Criterion only 1 field must be added
86  int addfields = 1;
87 
88  int npoints = m_f->m_exp[0]->GetNpoints();
89 
90  Array<OneD, Array<OneD, NekDouble> > grad(nfields * nfields);
91 
92  Array<OneD, Array<OneD, NekDouble> > omega(nfields * nfields);
93  Array<OneD, Array<OneD, NekDouble> > S (nfields * nfields);
94 
95  Array<OneD, Array<OneD, NekDouble> > outfield (addfields);
96  Array<OneD, Array<OneD, NekDouble> > outfield1(addfields);
97  Array<OneD, Array<OneD, NekDouble> > outfield2(addfields);
98  Array<OneD, Array<OneD, NekDouble> > outfield3(addfields);
99 
100  m_f->m_exp.resize(nfields+addfields);
101 
102  for (i = 0; i < nfields*nfields; ++i)
103  {
104  grad[i] = Array<OneD, NekDouble>(npoints);
105  }
106 
107  for (i = 0; i < addfields; ++i)
108  {
109  //Will store the Q-Criterion
110  outfield[i] = Array<OneD, NekDouble>(npoints);
111  outfield1[i] = Array<OneD, NekDouble>(npoints);
112  outfield2[i] = Array<OneD, NekDouble>(npoints);
113  outfield3[i] = Array<OneD, NekDouble>(npoints);
114 
115  omega[i] = Array<OneD, NekDouble>(npoints);
116  S[i] = Array<OneD, NekDouble>(npoints);
117  }
118 
119  for (i = 0; i < nfields; ++i)
120  {
121  m_f->m_exp[i]->PhysDeriv(m_f->m_exp[i]->GetPhys(),
122  grad[i*nfields],
123  grad[i*nfields+1],
124  grad[i*nfields+2]);
125  }
126 
127  // W_x = Wy - Vz
128  Vmath::Vsub(npoints, grad[2 * nfields + 1], 1,
129  grad[1 * nfields + 2], 1,
130  outfield1[0], 1);
131  // W_x^2
132  Vmath::Vmul(npoints, outfield1[0], 1,
133  outfield1[0], 1,
134  outfield1[0], 1);
135 
136  // W_y = Uz - Wx
137  Vmath::Vsub(npoints, grad[0 * nfields + 2], 1,
138  grad[2 * nfields + 0], 1,
139  outfield2[0], 1);
140  // W_y^2
141  Vmath::Vmul(npoints, outfield2[0], 1,
142  outfield2[0], 1,
143  outfield2[0], 1);
144 
145  // W_z = Vx - Uy
146  Vmath::Vsub(npoints, grad[1 * nfields + 0], 1,
147  grad[0 * nfields + 1], 1,
148  outfield3[0], 1);
149  // W_z^2
150  Vmath::Vmul(npoints, outfield3[0], 1,
151  outfield3[0], 1,
152  outfield3[0], 1);
153 
154  // add fields omega = 0.5*(W_x^2 + W_y^2 + W_z^2)
155 
156  NekDouble fac = 0.5;
157  Vmath::Vadd(npoints, &outfield1[0][0], 1,
158  &outfield2[0][0], 1,
159  &omega[0][0], 1);
160  Vmath::Vadd(npoints, &omega[0][0], 1,
161  &outfield3[0][0], 1,
162  &omega[0][0], 1);
163 
164  for (int k = 0; k < addfields; ++k)
165  {
166  Vmath::Smul(npoints, fac, &omega[k][0], 1, &omega[k][0], 1);
167  }
168 
169  Vmath::Zero(npoints, &outfield1[0][0], 1);
170  Vmath::Zero(npoints, &outfield2[0][0], 1);
171  Vmath::Zero(npoints, &outfield3[0][0], 1);
172 
173  Vmath::Vmul(npoints, grad[0 * nfields + 0], 1,
174  grad[0 * nfields + 0], 1,
175  outfield1[0], 1);
176  Vmath::Vmul(npoints, grad[1 * nfields + 1], 1,
177  grad[1 * nfields + 1], 1,
178  outfield2[0], 1);
179  Vmath::Vmul(npoints, grad[2 * nfields + 2], 1,
180  grad[2 * nfields + 2], 1,
181  outfield3[0], 1);
182 
183  Vmath::Vadd(npoints, &outfield1[0][0], 1,
184  &outfield2[0][0], 1,
185  &S[0][0], 1);
186  Vmath::Vadd(npoints, &S[0][0], 1,
187  &outfield3[0][0], 1,
188  &S[0][0], 1);
189 
190  // W_y + V_z
191  Vmath::Vadd(npoints, grad[2 * nfields + 1], 1,
192  grad[1 * nfields + 2], 1,
193  outfield1[0], 1);
194  Vmath::Vmul(npoints, &outfield1[0][0], 1,
195  &outfield1[0][0], 1,
196  &outfield1[0][0], 1);
197 
198  // U_z + W_x
199  Vmath::Vadd(npoints, grad[0 * nfields + 2], 1,
200  grad[2 * nfields + 0], 1,
201  outfield2[0], 1);
202  Vmath::Vmul(npoints, &outfield2[0][0], 1,
203  &outfield2[0][0], 1,
204  &outfield2[0][0], 1);
205 
206  // V_x + U_y
207  Vmath::Vadd(npoints, grad[1 * nfields + 0], 1,
208  grad[0 * nfields + 1], 1,
209  outfield3[0], 1);
210  Vmath::Vmul(npoints, &outfield3[0][0], 1,
211  &outfield3[0][0], 1,
212  &outfield3[0][0], 1);
213 
214  Vmath::Vadd(npoints, &outfield1[0][0], 1,
215  &outfield2[0][0], 1,
216  &outfield2[0][0], 1);
217  Vmath::Vadd(npoints, &outfield2[0][0], 1,
218  &outfield3[0][0], 1,
219  &outfield3[0][0], 1);
220 
221  for (int k = 0; k < addfields; ++k)
222  {
223  Vmath::Smul(npoints, fac, &outfield3[k][0], 1,
224  &outfield3[k][0], 1);
225  }
226 
227  Vmath::Vadd(npoints, &outfield3[0][0], 1, &S[0][0], 1, &S[0][0], 1);
228  Vmath::Vsub(npoints, omega[0], 1, S[0], 1, outfield[0], 1);
229 
230  for (int k = 0; k < addfields; ++k)
231  {
232  Vmath::Smul(npoints, fac, &outfield[k][0], 1,
233  &outfield[k][0], 1);
234  }
235 
236 
237  for (i = 0; i < addfields; ++i)
238  {
239  m_f->m_exp[nfields + i] = m_f->AppendExpList(m_f->m_fielddef[0]->m_numHomogeneousDir);
240  m_f->m_exp[nfields + i]->UpdatePhys() = outfield[i];
241  m_f->m_exp[nfields + i]->FwdTrans(outfield[i],
242  m_f->m_exp[nfields + i]->UpdateCoeffs());
243  }
244 
245  vector<string> outname;
246  outname.push_back("Q");
247 
248  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
249  = m_f->m_exp[0]->GetFieldDefinitions();
250  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
251 
252  for (j = 0; j < nfields + addfields; ++j)
253  {
254  for (i = 0; i < FieldDef.size(); ++i)
255  {
256  if (j >= nfields)
257  {
258  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
259  }
260  else
261  {
262  FieldDef[i]->m_fields.push_back(
263  m_f->m_fielddef[0]->m_fields[j]);
264  }
265  m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
266  }
267  }
268 
269  m_f->m_fielddef = FieldDef;
270  m_f->m_data = FieldData;
271  }
272  }
273 }
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2