Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Extrapolate.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Extrapolate.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: Abstract base class for Extrapolate.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <IncNavierStokesSolver/EquationSystems/Extrapolate.h>
37 #include <LibUtilities/Communication/Comm.h>
38 
39 namespace Nektar
40 {
41  NekDouble Extrapolate::StifflyStable_Betaq_Coeffs[3][3] = {
42  { 1.0, 0.0, 0.0},{ 2.0, -1.0, 0.0},{ 3.0, -3.0, 1.0}};
43  NekDouble Extrapolate::StifflyStable_Alpha_Coeffs[3][3] = {
44  { 1.0, 0.0, 0.0},{ 2.0, -0.5, 0.0},{ 3.0, -1.5, 1.0/3.0}};
45  NekDouble Extrapolate::StifflyStable_Gamma0_Coeffs[3] = {
46  1.0, 1.5, 11.0/6.0};
47 
48  ExtrapolateFactory& GetExtrapolateFactory()
49  {
50  typedef Loki::SingletonHolder<ExtrapolateFactory,
51  Loki::CreateUsingNew,
52  Loki::NoDestroy > Type;
53  return Type::Instance();
54  }
55 
56  Extrapolate::Extrapolate(
57  const LibUtilities::SessionReaderSharedPtr pSession,
58  Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
59  MultiRegions::ExpListSharedPtr pPressure,
60  const Array<OneD, int> pVel,
61  const AdvectionTermSharedPtr advObject)
62  : m_session(pSession),
63  m_fields(pFields),
64  m_pressure(pPressure),
65  m_velocity(pVel),
66  m_advObject(advObject)
67  {
68  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
69  m_comm = m_session->GetComm();
70  }
71 
72  Extrapolate::~Extrapolate()
73  {
74  }
75 
76 
77  std::string Extrapolate::def =
78  LibUtilities::SessionReader::RegisterDefaultSolverInfo(
79  "StandardExtrapolate", "StandardExtrapolate");
80 
81  /**
82  * Function to extrapolate the new pressure boundary condition.
83  * Based on the velocity field and on the advection term.
84  * Acceleration term is also computed.
85  * This routine is a general one for 2d and 3D application and it can be called
86  * directly from velocity correction scheme. Specialisation on dimensionality is
87  * redirected to the CalcPressureBCs method.
88  */
89  void Extrapolate::EvaluatePressureBCs(
90  const Array<OneD, const Array<OneD, NekDouble> > &fields,
91  const Array<OneD, const Array<OneD, NekDouble> > &N,
92  NekDouble kinvis)
93  {
94  if(m_HBCdata.num_elements()>0)
95  {
96  Array<OneD, NekDouble> tmp;
97  Array<OneD, NekDouble> accelerationTerm;
98 
99  m_pressureCalls++;
100  int n,cnt;
101  int nint = min(m_pressureCalls,m_intSteps);
102  int nlevels = m_pressureHBCs.num_elements();
103 
104  int acc_order = 0;
105 
106  accelerationTerm =
107  Array<OneD, NekDouble>(m_acceleration[0].num_elements(), 0.0);
108 
109  // Rotate HOPBCs storage
110  RollOver(m_pressureHBCs);
111 
112  // Rotate acceleration term
113  RollOver(m_acceleration);
114 
115  // Calculate Neumann BCs at current level
116  CalcNeumannPressureBCs(fields, N, kinvis);
117 
118  // Copy High order values into storage array
119  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
120  {
121  // High order boundary condition;
122  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHigh)
123  {
124  int nq = m_PBndExp[n]->GetNcoeffs();
125  Vmath::Vcopy(nq, &(m_PBndExp[n]->GetCoeffs()[0]), 1,
126  &(m_pressureHBCs[0])[cnt], 1);
127  cnt += nq;
128  }
129  }
130 
131  //Calculate acceleration term at level n based on previous steps
132  if (m_pressureCalls > 2)
133  {
134  acc_order = min(m_pressureCalls-2,m_intSteps);
135  Vmath::Smul(cnt, StifflyStable_Gamma0_Coeffs[acc_order-1],
136  m_acceleration[0], 1,
137  accelerationTerm, 1);
138 
139  for(int i = 0; i < acc_order; i++)
140  {
141  Vmath::Svtvp(cnt,
142  -1*StifflyStable_Alpha_Coeffs[acc_order-1][i],
143  m_acceleration[i+1], 1,
144  accelerationTerm, 1,
145  accelerationTerm, 1);
146  }
147  }
148 
149  // Adding acceleration term to HOPBCs
150  Vmath::Svtvp(cnt, -1.0/m_timestep,
151  accelerationTerm, 1,
152  m_pressureHBCs[0], 1,
153  m_pressureHBCs[0], 1);
154 
155  // Extrapolate to n+1
156  Vmath::Smul(cnt, StifflyStable_Betaq_Coeffs[nint-1][nint-1],
157  m_pressureHBCs[nint-1], 1,
158  m_pressureHBCs[nlevels-1], 1);
159 
160  for(n = 0; n < nint-1; ++n)
161  {
162  Vmath::Svtvp(cnt, StifflyStable_Betaq_Coeffs[nint-1][n],
163  m_pressureHBCs[n],1,m_pressureHBCs[nlevels-1],1,
164  m_pressureHBCs[nlevels-1],1);
165  }
166 
167  // Copy values of [dP/dn]^{n+1} in the pressure bcs storage.
168  // m_pressureHBCS[nlevels-1] will be cancelled at next time step
169  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
170  {
171  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHigh)
172  {
173  int nq = m_PBndExp[n]->GetNcoeffs();
174  Vmath::Vcopy(nq, &(m_pressureHBCs[nlevels-1])[cnt], 1,
175  &(m_PBndExp[n]->UpdateCoeffs()[0]), 1);
176  cnt += nq;
177  }
178  }
179 
180  }
181 
182  CalcOutflowBCs(fields, N, kinvis);
183  }
184 
185 
186  /**
187  * Unified routine for calculation high-oder terms
188  */
189  void Extrapolate::CalcNeumannPressureBCs(
190  const Array<OneD, const Array<OneD, NekDouble> > &fields,
191  const Array<OneD, const Array<OneD, NekDouble> > &N,
192  NekDouble kinvis)
193  {
194  Array<OneD, NekDouble> Pvals;
195  Array<OneD, NekDouble> Uvals;
196  StdRegions::StdExpansionSharedPtr Pbc;
197  StdRegions::StdExpansionSharedPtr elmt;
198 
199  Array<OneD, Array<OneD, const NekDouble> > Velocity(m_curl_dim);
200  Array<OneD, Array<OneD, const NekDouble> > Advection(m_bnd_dim);
201 
202  Array<OneD, Array<OneD, NekDouble> > BndValues(m_bnd_dim);
203  Array<OneD, Array<OneD, NekDouble> > Q(m_bnd_dim);
204 
205  for(int i = 0; i < m_bnd_dim; i++)
206  {
207  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
208  Q[i] = Array<OneD, NekDouble> (m_pressureBCsElmtMaxPts,0.0);
209  }
210 
211  for(int j = 0 ; j < m_HBCdata.num_elements() ; j++)
212  {
213  /// Casting the boundary expansion to the specific case
214  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
215  (m_PBndExp[m_HBCdata[j].m_bndryElmtID]
216  ->GetExp(m_HBCdata[j].m_bndElmtOffset));
217 
218  /// Picking up the element where the HOPBc is located
219  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
220 
221  /// Assigning
222  for(int i = 0; i < m_bnd_dim; i++)
223  {
224  Velocity[i] = fields[i] + m_HBCdata[j].m_physOffset;
225  Advection[i] = N[i] + m_HBCdata[j].m_physOffset;
226  }
227 
228  // for the 3DH1D case we need to grab the conjugate mode
229  if(m_pressure->GetExpType() == MultiRegions::e3DH1D)
230  {
231  Velocity[2] = fields[2] + m_HBCdata[j].m_assPhysOffset;
232  }
233 
234  /// Calculating the curl-curl and storing it in Q
235  CurlCurl(Velocity,Q,j);
236 
237  // Mounting advection component into the high-order condition
238  for(int i = 0; i < m_bnd_dim; i++)
239  {
240  MountHOPBCs(m_HBCdata[j].m_ptsInElmt,kinvis,Q[i],Advection[i]);
241  }
242 
243  Pvals = m_PBndExp[m_HBCdata[j].m_bndryElmtID]->UpdateCoeffs()
244  + m_PBndExp[m_HBCdata[j].m_bndryElmtID]
245  ->GetCoeff_Offset(m_HBCdata[j].m_bndElmtOffset);
246  Uvals = (m_acceleration[0]) + m_HBCdata[j].m_coeffOffset;
247 
248  // Getting values on the edge and filling the pressure boundary
249  // expansion and the acceleration term. Multiplication by the
250  // normal is required
251  switch(m_pressure->GetExpType())
252  {
253  case MultiRegions::e2D:
254  case MultiRegions::e3DH1D:
255  {
256  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
257  Q[0], BndValues[0]);
258  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
259  Q[1], BndValues[1]);
260  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
261  Pvals);
262 
263  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
264  Velocity[0], BndValues[0]);
265  elmt->GetEdgePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
266  Velocity[1], BndValues[1]);
267  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
268  Uvals);
269  }
270  break;
271  case MultiRegions::e3DH2D:
272  {
273  if(m_HBCdata[j].m_elmtTraceID == 0)
274  {
275  (m_PBndExp[m_HBCdata[j].m_bndryElmtID]->UpdateCoeffs()
276  + m_PBndExp[m_HBCdata[j].m_bndryElmtID]
277  ->GetCoeff_Offset(
278  m_HBCdata[j].m_bndElmtOffset))[0]
279  = -1.0*Q[0][0];
280  }
281  else if (m_HBCdata[j].m_elmtTraceID == 1)
282  {
283  (m_PBndExp[m_HBCdata[j].m_bndryElmtID]->UpdateCoeffs()
284  + m_PBndExp[m_HBCdata[j].m_bndryElmtID]
285  ->GetCoeff_Offset(
286  m_HBCdata[j].m_bndElmtOffset))[0]
287  = Q[0][m_HBCdata[j].m_ptsInElmt-1];
288  }
289  else
290  {
291  ASSERTL0(false,
292  "In the 3D homogeneous 2D approach BCs edge "
293  "ID can be just 0 or 1 ");
294  }
295  }
296  break;
297  case MultiRegions::e3D:
298  {
299  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
300  Q[0], BndValues[0]);
301  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
302  Q[1], BndValues[1]);
303  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
304  Q[2], BndValues[2]);
305  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
306  BndValues[2], Pvals);
307 
308  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
309  Velocity[0], BndValues[0]);
310  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
311  Velocity[1], BndValues[1]);
312  elmt->GetFacePhysVals(m_HBCdata[j].m_elmtTraceID, Pbc,
313  Velocity[2], BndValues[2]);
314  Pbc->NormVectorIProductWRTBase(BndValues[0], BndValues[1],
315  BndValues[2], Uvals);
316  }
317  break;
318  default:
319  ASSERTL0(0,"Dimension not supported");
320  break;
321  }
322  }
323  }
324 
325 
326 
327  void Extrapolate::CalcOutflowBCs(
328  const Array<OneD, const Array<OneD, NekDouble> > &fields,
329  const Array<OneD, const Array<OneD, NekDouble> > &N,
330  NekDouble kinvis)
331  {
332 
333  static bool init = true;
334  static bool noHOBC = false;
335 
336  if(noHOBC == true)
337  {
338  return;
339  }
340 
341  if(init) // set up storage for boundary velocity at outflow
342  {
343  init = false;
344  int totbndpts = 0;
345  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
346  {
347  if(m_PBndConds[n]->GetUserDefined()
348  == SpatialDomains::eHighOutflow)
349  {
350  totbndpts += m_PBndExp[n]->GetTotPoints();
351  }
352  }
353 
354  if(totbndpts == 0)
355  {
356  noHOBC = true;
357  return;
358  }
359 
360  m_outflowVel = Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (m_bnd_dim);
361  for(int i = 0; i < m_bnd_dim; ++i)
362  {
363  m_outflowVel[i] = Array<OneD, Array<OneD, NekDouble> >(m_curl_dim);
364  for(int j = 0; j < m_curl_dim; ++j)
365  {
366  // currently just set up for 2nd order extrapolation
367  m_outflowVel[i][j] = Array<OneD, NekDouble>(totbndpts,0.0);
368  }
369  }
370  }
371 
372  StdRegions::StdExpansionSharedPtr Bc,Pbc;
373  Array<OneD, Array<OneD, const SpatialDomains::BoundaryConditionShPtr> >
374  UBndConds(m_curl_dim);
375  Array<OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
376  UBndExp(m_curl_dim);
377 
378  for (int i = 0; i < m_curl_dim; ++i)
379  {
380  UBndConds[i] = m_fields[m_velocity[i]]->GetBndConditions();
381  UBndExp[i] = m_fields[m_velocity[i]]->GetBndCondExpansions();
382  }
383 
384  Array<OneD, Array<OneD, NekDouble> > BndValues(m_bnd_dim);
385  Array<OneD, Array<OneD, NekDouble> > BndElmt (m_bnd_dim);
386  Array<OneD, Array<OneD, NekDouble> > nGradu(m_bnd_dim);
387  Array<OneD, NekDouble > gradtmp (m_pressureBCsElmtMaxPts),
388  fgradtmp(m_pressureBCsElmtMaxPts);
389 
390  nGradu[0] = Array<OneD, NekDouble>(m_bnd_dim*m_pressureBCsMaxPts);
391  for(int i = 0; i < m_bnd_dim; ++i)
392  {
393  BndElmt[i] = Array<OneD, NekDouble> (m_pressureBCsElmtMaxPts,
394  0.0);
395  BndValues[i] = Array<OneD, NekDouble> (m_pressureBCsMaxPts,0.0);
396  nGradu[i] = nGradu[0] + i*m_pressureBCsMaxPts;
397  RollOver(m_outflowVel[i]);
398  }
399 
400  int nbc,cnt,cnt_start;
401  int veloffset = 0;
402  int nint = min(m_pressureCalls,m_intSteps);
403 
404  StdRegions::StdExpansionSharedPtr elmt;
405  Array<OneD, NekDouble> PBCvals, UBCvals;
406  Array<OneD, Array<OneD, NekDouble> > ubc(m_curl_dim);
407  Array<OneD, Array<OneD, NekDouble> > normals;
408 
409  cnt = 0;
410  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
411  {
412  // Do outflow boundary conditions if they exist
413  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHighOutflow)
414  {
415  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i,cnt++)
416  {
417  cnt = max(cnt,m_PBndExp[n]->GetTotPoints());
418  }
419  }
420  }
421 
422  for(int i =0; i < m_curl_dim; ++i)
423  {
424  ubc[i] = Array<OneD, NekDouble>(cnt);
425  }
426 
427  NekDouble U0,delta;
428  m_session->LoadParameter("U0_HighOrderBC",U0,1.0);
429  m_session->LoadParameter("Delta_HighOrderBC",delta,1/20.0);
430 
431  cnt = 0;
432  for(int n = 0; n < m_PBndConds.num_elements(); ++n)
433  {
434  cnt_start = cnt;
435 
436  // Do outflow boundary conditions if they exist
437  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHighOutflow)
438  {
439  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i,cnt++)
440  {
441  // find element and edge of this expansion.
442  Bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion>
443  (m_PBndExp[n]->GetExp(i));
444 
445  int elmtid = m_pressureBCtoElmtID[cnt];
446  elmt = m_fields[0]->GetExp(elmtid);
447  int offset = m_fields[0]->GetPhys_Offset(elmtid);
448 
449  // Determine extrapolated U,V values
450  int nq = elmt->GetTotPoints();
451 
452  // currently just using first order approximation here.
453  // previously have obtained value from m_integrationSoln
454  for(int j = 0; j < m_bnd_dim; ++j)
455  {
456  Vmath::Vcopy(nq, &fields[m_velocity[j]][offset], 1,
457  &BndElmt[j][0], 1);
458  }
459 
460  int nbc = m_PBndExp[n]->GetExp(i)->GetTotPoints();
461  int boundary = m_pressureBCtoTraceID[cnt];
462 
463  Array<OneD, NekDouble> veltmp, ptmp(nbc,0.0),
464  normDotu(nbc,0.0), utot(nbc,0.0),
465  divU(nbc,0.0);
466 
467  normals=elmt->GetSurfaceNormal(boundary);
468  Vmath::Zero(m_bnd_dim*m_pressureBCsMaxPts,nGradu[0],1);
469 
470  for (int j = 0; j < m_bnd_dim; j++)
471  {
472  // Calculate Grad u = du/dx, du/dy, du/dz, etc.
473  for (int k = 0; k< m_bnd_dim; k++)
474  {
475  elmt->PhysDeriv(MultiRegions::DirCartesianMap[k],
476  BndElmt[j], gradtmp);
477  elmt->GetTracePhysVals(boundary, Bc, gradtmp,
478  fgradtmp);
479  Vmath::Vvtvp(nbc,normals[k], 1, fgradtmp, 1,
480  nGradu[j], 1, nGradu[j],1);
481  if(j == k)
482  {
483  Vmath::Vadd(nbc,fgradtmp, 1, divU, 1, divU, 1);
484  }
485  }
486  }
487 
488  // extract velocity and store
489  for(int j = 0; j < m_bnd_dim; ++j)
490  {
491  elmt->GetTracePhysVals(boundary,Bc,BndElmt[j],
492  veltmp = m_outflowVel[j][0] + veloffset);
493  }
494 
495  // extrapolate velocity
496  if(nint <= 1)
497  {
498  for(int j = 0; j < m_bnd_dim; ++j)
499  {
500  Vmath::Vcopy(nbc,
501  veltmp = m_outflowVel[j][0] +veloffset, 1,
502  BndValues[j], 1);
503  }
504  }
505  else // only set up for 2nd order extrapolation
506  {
507  for(int j = 0; j < m_bnd_dim; ++j)
508  {
509  Vmath::Smul(nbc, 2.0,
510  veltmp = m_outflowVel[j][0] + veloffset, 1,
511  BndValues[j], 1);
512  Vmath::Svtvp(nbc, -1.0,
513  veltmp = m_outflowVel[j][1] + veloffset, 1,
514  BndValues[j], 1,
515  BndValues[j], 1);
516  }
517  }
518 
519  // Set up |u|^2, n.u, div(u), and (n.grad(u) . n) for
520  // pressure condition
521  for(int j = 0; j < m_bnd_dim; ++j)
522  {
523  Vmath::Vvtvp(nbc, BndValues[j], 1, BndValues[j], 1,
524  utot, 1, utot, 1);
525  Vmath::Vvtvp(nbc, normals[j], 1, BndValues[j], 1,
526  normDotu, 1, normDotu, 1);
527  Vmath::Vvtvp(nbc, normals[j], 1, nGradu[j], 1,
528  ptmp, 1, ptmp, 1);
529  }
530 
531  PBCvals = m_PBndExp[n]->GetPhys() +
532  m_PBndExp[n]->GetPhys_Offset(i);
533 
534  for(int k = 0; k < nbc; ++k)
535  {
536  NekDouble fac = 0.5*(1.0-tanh(normDotu[k]/(U0*delta)));
537 
538  // Set up Dirichlet pressure condition and
539  // store in ptmp (PBCvals contains a
540  // function from the input file )
541  ptmp[k] = kinvis * ptmp[k] - 0.5 * utot[k] * fac
542  - PBCvals[k];
543  }
544 
545  int u_offset = UBndExp[0][n]->GetPhys_Offset(i);
546 
547  for(int j = 0; j < m_bnd_dim; ++j)
548  {
549  UBCvals = UBndExp[j][n]->GetPhys()
550  + UBndExp[j][n]->GetPhys_Offset(i);
551 
552  for(int k = 0; k < nbc; ++k)
553  {
554  NekDouble fac = 0.5 * (1.0 - tanh(normDotu[k]
555  / (U0 * delta)));
556  ubc[j][k + u_offset] = (1.0 / kinvis)
557  * (UBCvals[k] + 0.5 * utot[k] * fac
558  * normals[j][k]);
559  }
560  }
561 
562  // set up pressure boundary condition
563  PBCvals = m_PBndExp[n]->UpdateCoeffs()
564  + m_PBndExp[n]->GetCoeff_Offset(i);
565  m_PBndExp[n]->GetExp(i)->FwdTrans(ptmp,PBCvals);
566 
567  veloffset += nbc;
568  }
569 
570  // Now set up Velocity conditions.
571  for(int j = 0; j < m_bnd_dim; j++)
572  {
573  if(UBndConds[j][n]->GetUserDefined()
574  == SpatialDomains::eHighOutflow)
575  {
576  cnt = cnt_start;
577 
578  for(int i = 0; i < UBndExp[0][n]->GetExpSize();
579  ++i, cnt++)
580  {
581  Pbc = StdRegions::StdExpansionSharedPtr
582  (m_PBndExp[n]->GetExp(i));
583  Bc = StdRegions::StdExpansionSharedPtr
584  (UBndExp[0][n]->GetExp(i));
585 
586  nbc = UBndExp[0][n]->GetExp(i)->GetTotPoints();
587  int boundary = m_pressureBCtoTraceID[cnt];
588 
589  Array<OneD, NekDouble> pb(nbc), ub(nbc);
590  int elmtid = m_pressureBCtoElmtID[cnt];
591 
592  elmt = m_fields[0]->GetExp(elmtid);
593 
594  normals = elmt->GetSurfaceNormal(boundary);
595 
596  // Get p from projected boundary condition
597  PBCvals = m_PBndExp[n]->UpdateCoeffs()
598  + m_PBndExp[n]->GetCoeff_Offset(i);
599  Pbc->BwdTrans(PBCvals,pb);
600 
601  int u_offset = UBndExp[j][n]->GetPhys_Offset(i);
602 
603  for(int k = 0; k < nbc; ++k)
604  {
605  ub[k] = ubc[j][k + u_offset]
606  + pb[k] * normals[j][k] / kinvis;
607  }
608 
609  UBCvals = UBndExp[j][n]->UpdateCoeffs()
610  + UBndExp[j][n]->GetCoeff_Offset(i);
611  Bc->IProductWRTBase(ub,UBCvals);
612  }
613  }
614  }
615  }
616  else
617  {
618  cnt += m_PBndExp[n]->GetExpSize();
619  }
620  }
621  }
622 
623 
624  /**
625  * Curl Curl routine - dimension dependent
626  */
627  void Extrapolate::CurlCurl(
628  Array<OneD, Array<OneD, const NekDouble> > &Vel,
629  Array<OneD, Array<OneD, NekDouble> > &Q,
630  const int j)
631  {
632  StdRegions::StdExpansionSharedPtr elmt
633  = m_fields[0]->GetExp(m_HBCdata[j].m_globalElmtID);
634 
635  Array<OneD,NekDouble> Vx(m_pressureBCsElmtMaxPts);
636  Array<OneD,NekDouble> Uy(m_pressureBCsElmtMaxPts);
637 
638  switch(m_fields[0]->GetExpType())
639  {
640  case MultiRegions::e2D:
641  {
642  Array<OneD,NekDouble> Dummy(m_pressureBCsElmtMaxPts);
643 
644  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
645  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vel[0], Uy);
646 
647  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1, Dummy, 1);
648 
649  elmt->PhysDeriv(Dummy,Q[1],Q[0]);
650 
651  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, -1.0, Q[1], 1, Q[1], 1);
652  }
653  break;
654 
655  case MultiRegions::e3DH1D:
656  {
657  Array<OneD,NekDouble> Wz(m_pressureBCsElmtMaxPts);
658 
659  Array<OneD,NekDouble> Dummy1(m_pressureBCsElmtMaxPts);
660  Array<OneD,NekDouble> Dummy2(m_pressureBCsElmtMaxPts);
661 
662  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
663  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vel[0], Uy);
664  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
665  Vel[2], 1, Wz, 1);
666 
667  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Vx, Dummy1);
668  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Uy, Dummy2);
669  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Dummy1, 1, Dummy2, 1,
670  Q[0], 1);
671  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
672  Vel[0], 1, Dummy1, 1);
673  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Wz, Dummy2);
674  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Q[0], 1, Dummy1, 1,
675  Q[0], 1);
676  Vmath::Vadd(m_HBCdata[j].m_ptsInElmt, Q[0], 1, Dummy2, 1,
677  Q[0], 1);
678 
679  elmt->PhysDeriv(MultiRegions::DirCartesianMap[1], Wz, Dummy1);
680  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
681  Vel[1], 1, Dummy2, 1);
682  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Dummy1, 1, Dummy2, 1,
683  Q[1], 1);
684  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vx, Dummy1);
685  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Uy, Dummy2);
686  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Q[1], 1, Dummy1, 1,
687  Q[1], 1);
688  Vmath::Vadd(m_HBCdata[j].m_ptsInElmt, Q[1], 1, Dummy2, 1,
689  Q[1], 1);
690  }
691  break;
692 
693  case MultiRegions::e3DH2D:
694  {
695  Array<OneD,NekDouble> Wx(m_pressureBCsElmtMaxPts);
696  Array<OneD,NekDouble> Wz(m_pressureBCsElmtMaxPts);
697  Array<OneD,NekDouble> Uz(m_pressureBCsElmtMaxPts);
698  Array<OneD,NekDouble> qz(m_pressureBCsElmtMaxPts);
699  Array<OneD,NekDouble> qy(m_pressureBCsElmtMaxPts);
700 
701  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[2], Wx);
702  elmt->PhysDeriv(MultiRegions::DirCartesianMap[0], Vel[1], Vx);
703 
704  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
705  Vel[0], 1, Uy, 1);
706 
707  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
708  Vel[0], 1, Uz, 1);
709 
710  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wz, 1, Wx, 1,
711  qy, 1);
712  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1,
713  qz, 1);
714 
715  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_negWavenumberSq[j],
716  qz, 1, Uy, 1);
717 
718  Vmath::Smul(m_HBCdata[j].m_ptsInElmt, m_wavenumber[j],
719  qy, 1, Uz, 1);
720 
721  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uy, 1, Uz, 1,
722  Q[0], 1);
723  }
724  break;
725 
726  case MultiRegions::e3D:
727  {
728  Array<OneD,NekDouble> Dummy(m_pressureBCsElmtMaxPts);
729  Array<OneD,NekDouble> Vz(m_pressureBCsElmtMaxPts);
730  Array<OneD,NekDouble> Uz(m_pressureBCsElmtMaxPts);
731  Array<OneD,NekDouble> Wx(m_pressureBCsElmtMaxPts);
732  Array<OneD,NekDouble> Wy(m_pressureBCsElmtMaxPts);
733 
734  elmt->PhysDeriv(Vel[0], Dummy, Uy, Uz);
735  elmt->PhysDeriv(Vel[1], Vx, Dummy, Vz);
736  elmt->PhysDeriv(Vel[2], Wx, Wy, Dummy);
737 
738  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wy, 1, Vz, 1, Q[0], 1);
739  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uz, 1, Wx, 1, Q[1], 1);
740  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Uy, 1, Q[2], 1);
741 
742  elmt->PhysDeriv(Q[0], Dummy, Wy, Vx);
743  elmt->PhysDeriv(Q[1], Wx, Dummy, Uz);
744  elmt->PhysDeriv(Q[2], Vz, Uy, Dummy);
745 
746  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Uy, 1, Uz, 1, Q[0], 1);
747  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Vx, 1, Vz, 1, Q[1], 1);
748  Vmath::Vsub(m_HBCdata[j].m_ptsInElmt, Wx, 1, Wy, 1, Q[2], 1);
749  }
750  break;
751  default:
752  ASSERTL0(0,"Dimension not supported");
753  break;
754  }
755  }
756 
757 
758  /**
759  * Function to roll time-level storages to the next step layout.
760  * The stored data associated with the oldest time-level
761  * (not required anymore) are moved to the top, where they will
762  * be overwritten as the solution process progresses.
763  */
764  void Extrapolate::RollOver(Array<OneD, Array<OneD, NekDouble> > &input)
765  {
766  int nlevels = input.num_elements();
767 
768  Array<OneD, NekDouble> tmp;
769 
770  tmp = input[nlevels-1];
771 
772  for(int n = nlevels-1; n > 0; --n)
773  {
774  input[n] = input[n-1];
775  }
776 
777  input[0] = tmp;
778  }
779 
780 
781  /**
782  * Map to directly locate HOPBCs position and offsets in all scenarios
783  */
784  void Extrapolate::GenerateHOPBCMap()
785  {
786  m_PBndConds = m_pressure->GetBndConditions();
787  m_PBndExp = m_pressure->GetBndCondExpansions();
788 
789  // Set up mapping from pressure boundary condition to pressure element
790  // details.
791  m_pressure->GetBoundaryToElmtMap(m_pressureBCtoElmtID,
792  m_pressureBCtoTraceID);
793 
794  // find the maximum values of points for pressure BC evaluation
795  m_pressureBCsMaxPts = 0;
796  m_pressureBCsElmtMaxPts = 0;
797  int cnt, n;
798  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
799  {
800  for(int i = 0; i < m_PBndExp[n]->GetExpSize(); ++i)
801  {
802  m_pressureBCsMaxPts = max(m_pressureBCsMaxPts,
803  m_PBndExp[n]->GetExp(i)->GetTotPoints());
804  m_pressureBCsElmtMaxPts = max(m_pressureBCsElmtMaxPts,
805  m_pressure->GetExp(m_pressureBCtoElmtID[cnt++])
806  ->GetTotPoints());
807  }
808  }
809 
810  // Storage array for high order pressure BCs
811  m_pressureHBCs = Array<OneD, Array<OneD, NekDouble> > (m_intSteps);
812  m_acceleration = Array<OneD, Array<OneD, NekDouble> > (m_intSteps + 1);
813 
814  int HBCnumber = 0;
815  for(cnt = n = 0; n < m_PBndConds.num_elements(); ++n)
816  {
817  // High order boundary condition;
818  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHigh)
819  {
820  cnt += m_PBndExp[n]->GetNcoeffs();
821  HBCnumber += m_PBndExp[n]->GetExpSize();
822  }
823  }
824 
825  int checkHBC = HBCnumber;
826  m_comm->AllReduce(checkHBC,LibUtilities::ReduceSum);
827  //ASSERTL0(checkHBC > 0 ,"At least one high-order pressure boundary "
828  // "condition is required for scheme "
829  // "consistency");
830 
831  m_acceleration[0] = Array<OneD, NekDouble>(cnt, 0.0);
832  for(n = 0; n < m_intSteps; ++n)
833  {
834  m_pressureHBCs[n] = Array<OneD, NekDouble>(cnt, 0.0);
835  m_acceleration[n+1] = Array<OneD, NekDouble>(cnt, 0.0);
836  }
837 
838  m_pressureCalls = 0;
839 
840  switch(m_pressure->GetExpType())
841  {
842  case MultiRegions::e2D:
843  {
844  m_curl_dim = 2;
845  m_bnd_dim = 2;
846  }
847  break;
848  case MultiRegions::e3DH1D:
849  {
850  m_curl_dim = 3;
851  m_bnd_dim = 2;
852  }
853  break;
854  case MultiRegions::e3DH2D:
855  {
856  m_curl_dim = 3;
857  m_bnd_dim = 1;
858  }
859  break;
860  case MultiRegions::e3D:
861  {
862  m_curl_dim = 3;
863  m_bnd_dim = 3;
864  }
865  break;
866  default:
867  ASSERTL0(0,"Dimension not supported");
868  break;
869  }
870 
871 
872  m_HBCdata = Array<OneD, HBCInfo>(HBCnumber);
873  StdRegions::StdExpansionSharedPtr elmt;
874 
875  switch(m_pressure->GetExpType())
876  {
877  case MultiRegions::e2D:
878  case MultiRegions::e3D:
879  {
880  int coeff_count = 0;
881  int exp_size;
882  int j=0;
883  int cnt = 0;
884  for(int n = 0 ; n < m_PBndConds.num_elements(); ++n)
885  {
886  exp_size = m_PBndExp[n]->GetExpSize();
887 
888  if(m_PBndConds[n]->GetUserDefined()
889  == SpatialDomains::eHigh)
890  {
891  for(int i = 0; i < exp_size; ++i,cnt++)
892  {
893  m_HBCdata[j].m_globalElmtID =
894  m_pressureBCtoElmtID[cnt];
895  elmt = m_pressure->GetExp(
896  m_HBCdata[j].m_globalElmtID);
897  m_HBCdata[j].m_ptsInElmt =
898  elmt->GetTotPoints();
899  m_HBCdata[j].m_physOffset =
900  m_pressure->GetPhys_Offset(
901  m_HBCdata[j].m_globalElmtID);
902  m_HBCdata[j].m_bndElmtOffset = i;
903  m_HBCdata[j].m_elmtTraceID =
904  m_pressureBCtoTraceID[cnt];
905  m_HBCdata[j].m_bndryElmtID = n;
906  m_HBCdata[j].m_coeffOffset = coeff_count;
907  coeff_count += elmt->GetEdgeNcoeffs(
908  m_HBCdata[j].m_elmtTraceID);
909  j = j+1;
910  }
911  }
912  else // setting if just standard BC no High order
913  {
914  cnt += exp_size;
915  }
916  }
917  }
918  break;
919 
920  case MultiRegions::e3DH1D:
921  {
922  Array<OneD, unsigned int> planes;
923  planes = m_pressure->GetZIDs();
924  int num_planes = planes.num_elements();
925  int num_elm_per_plane = (m_pressure->GetExpSize())/num_planes;
926 
927  m_wavenumber = Array<OneD, NekDouble>(HBCnumber);
928  m_negWavenumberSq = Array<OneD, NekDouble>(HBCnumber);
929 
930  int coeff_count = 0;
931  int exp_size, exp_size_per_plane;
932  int j=0;
933  int K;
934  NekDouble sign = -1.0;
935  int cnt = 0;
936 
937  m_session->MatchSolverInfo("ModeType", "SingleMode",
938  m_SingleMode, false);
939  m_session->MatchSolverInfo("ModeType", "HalfMode",
940  m_HalfMode, false);
941  m_session->MatchSolverInfo("ModeType", "MultipleModes",
942  m_MultipleModes, false);
943  m_session->LoadParameter("LZ", m_LhomZ);
944 
945  // Stability Analysis flags
946  if(m_session->DefinesSolverInfo("ModeType"))
947  {
948  if(m_SingleMode)
949  {
950  m_npointsZ = 2;
951  }
952  else if(m_HalfMode)
953  {
954  m_npointsZ = 1;
955  }
956  else if(m_MultipleModes)
957  {
958  m_npointsZ = m_session->GetParameter("HomModesZ");
959  }
960  else
961  {
962  ASSERTL0(false, "SolverInfo ModeType not valid");
963  }
964  }
965  else
966  {
967  m_npointsZ = m_session->GetParameter("HomModesZ");
968  }
969 
970  for(int k = 0; k < num_planes; k++)
971  {
972  K = planes[k]/2;
973  for(int n = 0 ; n < m_PBndConds.num_elements(); ++n)
974  {
975  exp_size = m_PBndExp[n]->GetExpSize();
976  exp_size_per_plane = exp_size/num_planes;
977 
978  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHigh)
979  {
980  for(int i = 0; i < exp_size_per_plane; ++i,cnt++)
981  {
982  m_HBCdata[j].m_globalElmtID = m_pressureBCtoElmtID[cnt];
983  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
984  m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
985  m_HBCdata[j].m_physOffset = m_pressure->GetPhys_Offset(m_HBCdata[j].m_globalElmtID);
986  m_HBCdata[j].m_bndElmtOffset = i+k*exp_size_per_plane;
987  m_HBCdata[j].m_elmtTraceID = m_pressureBCtoTraceID[cnt];
988  m_HBCdata[j].m_bndryElmtID = n;
989  m_HBCdata[j].m_coeffOffset = coeff_count;
990  coeff_count += elmt->GetEdgeNcoeffs(m_HBCdata[j].m_elmtTraceID);
991 
992  if(m_SingleMode)
993  {
994  m_wavenumber[j] = -2*M_PI/m_LhomZ;
995  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
996  }
997  else if(m_HalfMode || m_MultipleModes)
998  {
999  m_wavenumber[j] = 2*M_PI/m_LhomZ;
1000  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
1001  }
1002  else
1003  {
1004  m_wavenumber[j] = 2*M_PI*sign*(NekDouble(K))/m_LhomZ;
1005  m_negWavenumberSq[j] = -1.0*m_wavenumber[j]*m_wavenumber[j];
1006  }
1007 
1008  int assElmtID;
1009 
1010  if(k%2==0)
1011  {
1012  if(m_HalfMode)
1013  {
1014  assElmtID = m_HBCdata[j].m_globalElmtID;
1015 
1016  }
1017  else
1018  {
1019  assElmtID = m_HBCdata[j].m_globalElmtID + num_elm_per_plane;
1020  }
1021  }
1022  else
1023  {
1024  assElmtID = m_HBCdata[j].m_globalElmtID - num_elm_per_plane;
1025  }
1026 
1027  m_HBCdata[j].m_assPhysOffset = m_pressure->GetPhys_Offset(assElmtID);
1028 
1029  j = j+1;
1030  }
1031  }
1032  else // setting if just standard BC no High order
1033  {
1034  cnt += exp_size_per_plane;
1035  }
1036  }
1037  sign = -1.0*sign;
1038  }
1039  }
1040  break;
1041 
1042  case MultiRegions::e3DH2D:
1043  {
1044  int cnt = 0;
1045  int exp_size, exp_size_per_line;
1046  int j=0;
1047 
1048  for(int k1 = 0; k1 < m_npointsZ; k1++)
1049  {
1050  for(int k2 = 0; k2 < m_npointsY; k2++)
1051  {
1052  for(int n = 0 ; n < m_PBndConds.num_elements(); ++n)
1053  {
1054  exp_size = m_PBndExp[n]->GetExpSize();
1055 
1056  exp_size_per_line = exp_size/(m_npointsZ*m_npointsY);
1057 
1058  if(m_PBndConds[n]->GetUserDefined() == SpatialDomains::eHigh)
1059  {
1060  for(int i = 0; i < exp_size_per_line; ++i,cnt++)
1061  {
1062  // find element and edge of this expansion.
1063  m_HBCdata[j].m_globalElmtID = m_pressureBCtoElmtID[cnt];
1064  elmt = m_pressure->GetExp(m_HBCdata[j].m_globalElmtID);
1065  m_HBCdata[j].m_ptsInElmt = elmt->GetTotPoints();
1066  m_HBCdata[j].m_physOffset = m_pressure->GetPhys_Offset(m_HBCdata[j].m_globalElmtID);
1067  m_HBCdata[j].m_bndElmtOffset = i+(k1*m_npointsY+k2)*exp_size_per_line;
1068  m_HBCdata[j].m_elmtTraceID = m_pressureBCtoTraceID[cnt];
1069  m_HBCdata[j].m_bndryElmtID = n;
1070  }
1071  }
1072  else
1073  {
1074  cnt += exp_size_per_line;
1075  }
1076  }
1077  }
1078  }
1079  }
1080  break;
1081  default:
1082  ASSERTL0(0,"Dimension not supported");
1083  break;
1084  }
1085  }
1086 
1087  /**
1088  *
1089  */
1090  Array<OneD, NekDouble> Extrapolate::GetMaxStdVelocity(
1091  const Array<OneD, Array<OneD,NekDouble> > inarray)
1092  {
1093  // Checking if the problem is 2D
1094  ASSERTL0(m_curl_dim >= 2, "Method not implemented for 1D");
1095 
1096  int n_points_0 = m_fields[0]->GetExp(0)->GetTotPoints();
1097  int n_element = m_fields[0]->GetExpSize();
1098  int nvel = inarray.num_elements();
1099  int cnt;
1100 
1101  NekDouble pntVelocity;
1102 
1103  // Getting the standard velocity vector on the 2D normal space
1104  Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
1105  Array<OneD, NekDouble> maxV(n_element, 0.0);
1106  LibUtilities::PointsKeyVector ptsKeys;
1107 
1108  for (int i = 0; i < nvel; ++i)
1109  {
1110  stdVelocity[i] = Array<OneD, NekDouble>(n_points_0);
1111  }
1112 
1113  if (nvel == 2)
1114  {
1115  cnt = 0.0;
1116  for (int el = 0; el < n_element; ++el)
1117  {
1118  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
1119  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
1120 
1121  // reset local space if necessary
1122  if(n_points != n_points_0)
1123  {
1124  for (int j = 0; j < nvel; ++j)
1125  {
1126  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
1127  }
1128  n_points_0 = n_points;
1129  }
1130 
1131  Array<TwoD, const NekDouble> gmat =
1132  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1133 
1134  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1135  == SpatialDomains::eDeformed)
1136  {
1137  for (int i = 0; i < n_points; i++)
1138  {
1139  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1140  + gmat[2][i]*inarray[1][i+cnt];
1141 
1142  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1143  + gmat[3][i]*inarray[1][i+cnt];
1144  }
1145  }
1146  else
1147  {
1148  for (int i = 0; i < n_points; i++)
1149  {
1150  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1151  + gmat[2][0]*inarray[1][i+cnt];
1152 
1153  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1154  + gmat[3][0]*inarray[1][i+cnt];
1155  }
1156  }
1157 
1158  cnt += n_points;
1159 
1160 
1161  for (int i = 0; i < n_points; i++)
1162  {
1163  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1164  + stdVelocity[1][i]*stdVelocity[1][i];
1165 
1166  if (pntVelocity>maxV[el])
1167  {
1168  maxV[el] = pntVelocity;
1169  }
1170  }
1171  maxV[el] = sqrt(maxV[el]);
1172  }
1173  }
1174  else
1175  {
1176  cnt = 0;
1177  for (int el = 0; el < n_element; ++el)
1178  {
1179 
1180  int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
1181  ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();
1182 
1183  // reset local space if necessary
1184  if(n_points != n_points_0)
1185  {
1186  for (int j = 0; j < nvel; ++j)
1187  {
1188  stdVelocity[j] = Array<OneD, NekDouble>(n_points);
1189  }
1190  n_points_0 = n_points;
1191  }
1192 
1193  Array<TwoD, const NekDouble> gmat =
1194  m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetDerivFactors(ptsKeys);
1195 
1196  if (m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()->GetGtype()
1197  == SpatialDomains::eDeformed)
1198  {
1199  for (int i = 0; i < n_points; i++)
1200  {
1201  stdVelocity[0][i] = gmat[0][i]*inarray[0][i+cnt]
1202  + gmat[3][i]*inarray[1][i+cnt]
1203  + gmat[6][i]*inarray[2][i+cnt];
1204 
1205  stdVelocity[1][i] = gmat[1][i]*inarray[0][i+cnt]
1206  + gmat[4][i]*inarray[1][i+cnt]
1207  + gmat[7][i]*inarray[2][i+cnt];
1208 
1209  stdVelocity[2][i] = gmat[2][i]*inarray[0][i+cnt]
1210  + gmat[5][i]*inarray[1][i+cnt]
1211  + gmat[8][i]*inarray[2][i+cnt];
1212  }
1213  }
1214  else
1215  {
1216  for (int i = 0; i < n_points; i++)
1217  {
1218  stdVelocity[0][i] = gmat[0][0]*inarray[0][i+cnt]
1219  + gmat[3][0]*inarray[1][i+cnt]
1220  + gmat[6][0]*inarray[2][i+cnt];
1221 
1222  stdVelocity[1][i] = gmat[1][0]*inarray[0][i+cnt]
1223  + gmat[4][0]*inarray[1][i+cnt]
1224  + gmat[7][0]*inarray[2][i+cnt];
1225 
1226  stdVelocity[2][i] = gmat[2][0]*inarray[0][i+cnt]
1227  + gmat[5][0]*inarray[1][i+cnt]
1228  + gmat[8][0]*inarray[2][i+cnt];
1229  }
1230  }
1231 
1232  cnt += n_points;
1233 
1234  for (int i = 0; i < n_points; i++)
1235  {
1236  pntVelocity = stdVelocity[0][i]*stdVelocity[0][i]
1237  + stdVelocity[1][i]*stdVelocity[1][i]
1238  + stdVelocity[2][i]*stdVelocity[2][i];
1239 
1240  if (pntVelocity > maxV[el])
1241  {
1242  maxV[el] = pntVelocity;
1243  }
1244  }
1245 
1246  maxV[el] = sqrt(maxV[el]);
1247  //cout << maxV[el]*maxV[el] << endl;
1248  }
1249  }
1250 
1251  return maxV;
1252  }
1253 
1254 }
1255 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2