Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SubSteppingExtrapolate.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SubSteppingExtrapolate.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 SubSteppingExtrapolate.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 namespace Nektar
40 {
41  /**
42  * Registers the class with the Factory.
43  */
45  "SubStepping",
47  "SubStepping");
48 
51  Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
53  const Array<OneD, int> pVel,
54  const SolverUtils::AdvectionSharedPtr advObject)
55  : Extrapolate(pSession,pFields,pPressure,pVel,advObject)
56  {
57  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
58  m_session->LoadParameter("CFL", m_cflSafetyFactor, 0.5);
59  m_session->LoadParameter("MinSubSteps", minsubsteps,0);
60  }
61 
63  {
64  }
65 
67  int intMethod,
68  const LibUtilities::TimeIntegrationWrapperSharedPtr &IntegrationScheme)
69  {
70  int i;
71 
72  // Set to 1 for first step and it will then be increased in
73  // time advance routines
74  switch(intMethod)
75  {
78  {
80 
81  // Fields for linear interpolation
82  m_previousVelFields = Array<OneD, Array<OneD, NekDouble> >(2*m_fields.num_elements());
83  int ntotpts = m_fields[0]->GetTotPoints();
84  m_previousVelFields[0] = Array<OneD, NekDouble>(2*m_fields.num_elements()*ntotpts);
85  for(i = 1; i < 2*m_fields.num_elements(); ++i)
86  {
87  m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
88  }
89 
90  }
91  break;
93  {
95 
96  int nvel = m_velocity.num_elements();
97 
98  // Fields for quadratic interpolation
99  m_previousVelFields = Array<OneD, Array<OneD, NekDouble> >(3*nvel);
100 
101  int ntotpts = m_fields[0]->GetTotPoints();
102  m_previousVelFields[0] = Array<OneD, NekDouble>(3*nvel*ntotpts);
103  for(i = 1; i < 3*nvel; ++i)
104  {
105  m_previousVelFields[i] = m_previousVelFields[i-1] + ntotpts;
106  }
107 
108  }
109  break;
110  default:
111  ASSERTL0(0,"Integration method not suitable: Options include BackwardEuler or BDFImplicitOrder1");
112  break;
113  }
114 
115  m_intSteps = IntegrationScheme->GetIntegrationSteps();
116 
117  // set explicit time-integration class operators
120  }
121 
122  /**
123  * Explicit Advection terms used by SubStepAdvance time integration
124  */
126  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
127  Array<OneD, Array<OneD, NekDouble> > &outarray,
128  const NekDouble time)
129  {
130  int i;
131  int nVariables = inarray.num_elements();
132  int nQuadraturePts = inarray[0].num_elements();
133 
134  /// Get the number of coefficients
135  int ncoeffs = m_fields[0]->GetNcoeffs();
136 
137  /// Define an auxiliary variable to compute the RHS
138  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nVariables);
139  WeakAdv[0] = Array<OneD, NekDouble> (ncoeffs*nVariables);
140  for(i = 1; i < nVariables; ++i)
141  {
142  WeakAdv[i] = WeakAdv[i-1] + ncoeffs;
143  }
144 
145  Array<OneD, Array<OneD, NekDouble> > Velfields(m_velocity.num_elements());
146  Array<OneD, int> VelIds(m_velocity.num_elements());
147 
148  Velfields[0] = Array<OneD, NekDouble> (nQuadraturePts*m_velocity.num_elements());
149 
150  for(i = 1; i < m_velocity.num_elements(); ++i)
151  {
152  Velfields[i] = Velfields[i-1] + nQuadraturePts;
153  }
154 
155  SubStepExtrapolateField(fmod(time,m_timestep), Velfields);
156 
157  m_advObject->Advect(m_velocity.num_elements(), m_fields, Velfields, inarray, outarray, time);
158 
159  for(i = 0; i < nVariables; ++i)
160  {
161  m_fields[i]->IProductWRTBase(outarray[i],WeakAdv[i]);
162  // negation requried due to sign of DoAdvection term to be consistent
163  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
164  }
165 
166  AddAdvectionPenaltyFlux(Velfields, inarray, WeakAdv);
167 
168  /// Operations to compute the RHS
169  for(i = 0; i < nVariables; ++i)
170  {
171  // Negate the RHS
172  Vmath::Neg(ncoeffs, WeakAdv[i], 1);
173 
174  /// Multiply the flux by the inverse of the mass matrix
175  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
176 
177  /// Store in outarray the physical values of the RHS
178  m_fields[i]->BwdTrans(WeakAdv[i], outarray[i]);
179  }
180 
181  //add the force
182  /*if(m_session->DefinesFunction("BodyForce"))
183  {
184  if(m_SingleMode || m_HalfMode)
185  {
186  for(int i = 0; i < m_nConvectiveFields; ++i)
187  {
188  m_forces[i]->SetWaveSpace(true);
189  m_forces[i]->BwdTrans(m_forces[i]->GetCoeffs(),
190  m_forces[i]->UpdatePhys());
191  }
192  }
193 
194  int nqtot = m_fields[0]->GetTotPoints();
195  for(int i = 0; i < m_nConvectiveFields; ++i)
196  {
197  Vmath::Vadd(nqtot,outarray[i],1,(m_forces[i]->GetPhys()),1,outarray[i],1);
198  }
199  }*/
200  }
201 
202  /**
203  * Projection used by SubStepAdvance time integration
204  */
206  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
207  Array<OneD, Array<OneD, NekDouble> > &outarray,
208  const NekDouble time)
209  {
210  ASSERTL1(inarray.num_elements() == outarray.num_elements(),"Inarray and outarray of different sizes ");
211 
212  for(int i = 0; i < inarray.num_elements(); ++i)
213  {
214  Vmath::Vcopy(inarray[i].num_elements(),inarray[i],1,outarray[i],1);
215  }
216  }
217 
218 
219  /**
220  *
221  */
223  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
224  const NekDouble Aii_Dt,
225  NekDouble kinvis)
226  {
227  int nConvectiveFields =m_fields.num_elements()-1;
228  Array<OneD, Array<OneD, NekDouble> > velfields(nConvectiveFields);
229 
230  for(int i = 0; i < nConvectiveFields; ++i)
231  {
232  velfields[i] = m_fields[m_velocity[i]]->GetPhys();
233  }
234 
235  EvaluatePressureBCs(velfields,inarray,kinvis);
236 
237  AddDuDt(inarray,Aii_Dt);
238  }
239 
240  /**
241  *
242  */
244  const int nstep)
245  {
246  int i,n;
247  int nvel = m_velocity.num_elements();
248  int npts = m_fields[0]->GetTotPoints();
249 
250  // rotate fields
251  int nblocks = m_previousVelFields.num_elements()/nvel;
252  Array<OneD, NekDouble> save;
253 
254  for(n = 0; n < nvel; ++n)
255  {
256  save = m_previousVelFields[(nblocks-1)*nvel+n];
257 
258  for(i = nblocks-1; i > 0; --i)
259  {
260  m_previousVelFields[i*nvel+n] = m_previousVelFields[(i-1)*nvel+n];
261  }
262 
263  m_previousVelFields[n] = save;
264  }
265 
266  // Put previous field
267  for(i = 0; i < nvel; ++i)
268  {
269  m_fields[m_velocity[i]]->BwdTrans(m_fields[m_velocity[i]]->GetCoeffs(),
270  m_fields[m_velocity[i]]->UpdatePhys());
271  Vmath::Vcopy(npts,m_fields[m_velocity[i]]->GetPhys(),1,
272  m_previousVelFields[i],1);
273  }
274 
275  if(nstep == 0)// initialise all levels with first field
276  {
277  for(n = 0; n < nvel; ++n)
278  {
279  for(i = 1; i < nblocks; ++i)
280  {
281  Vmath::Vcopy(npts,m_fields[m_velocity[n]]->GetPhys(),1,
282  m_previousVelFields[i*nvel+n],1);
283 
284  }
285  }
286  }
287  }
288 
289  /**
290  *
291  */
293  const LibUtilities::TimeIntegrationSolutionSharedPtr &integrationSoln,
294  int nstep,
295  NekDouble time)
296  {
297  int n;
298  int nsubsteps;
299 
300  NekDouble dt;
301 
302  Array<OneD, Array<OneD, NekDouble> > fields, velfields;
303 
304  static int ncalls = 1;
305  int nint = min(ncalls++, m_intSteps);
306 
307  Array<OneD, NekDouble> CFL(m_fields[0]->GetExpSize(),
309  //this needs to change
310  m_comm = m_fields[0]->GetComm()->GetRowComm();
311 
312  // Get the proper time step with CFL control
313  dt = GetSubstepTimeStep();
314 
315  nsubsteps = (m_timestep > dt)? ((int)(m_timestep/dt)+1):1;
316  nsubsteps = max(minsubsteps, nsubsteps);
317 
318  dt = m_timestep/nsubsteps;
319 
320  if (m_infosteps && !((nstep+1)%m_infosteps) && m_comm->GetRank() == 0)
321  {
322  cout << "Sub-integrating using "<< nsubsteps
323  << " steps over Dt = " << m_timestep
324  << " (SubStep CFL=" << m_cflSafetyFactor << ")"<< endl;
325  }
326 
327  for (int m = 0; m < nint; ++m)
328  {
329  // We need to update the fields held by the m_integrationSoln
330  fields = integrationSoln->UpdateSolutionVector()[m];
331 
332  // Initialise NS solver which is set up to use a GLM method
333  // with calls to EvaluateAdvection_SetPressureBCs and
334  // SolveUnsteadyStokesSystem
336  SubIntegrationSoln = m_subStepIntegrationScheme->
337  InitializeScheme(dt, fields, time, m_subStepIntegrationOps);
338 
339  for(n = 0; n < nsubsteps; ++n)
340  {
341  fields = m_subStepIntegrationScheme->TimeIntegrate(
342  n, dt, SubIntegrationSoln,
344  }
345 
346  // Reset time integrated solution in m_integrationSoln
347  integrationSoln->SetSolVector(m,fields);
348  }
349  }
350 
351 
352  /**
353  *
354  */
356  {
357  int n_element = m_fields[0]->GetExpSize();
358 
359  const Array<OneD, int> ExpOrder=m_fields[0]->EvalBasisNumModesMaxPerExp();
360  Array<OneD, int> ExpOrderList (n_element, ExpOrder);
361 
362  const NekDouble cLambda = 0.2; // Spencer book pag. 317
363 
364  Array<OneD, NekDouble> tstep (n_element, 0.0);
365  Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
366  Array<OneD, Array<OneD, NekDouble> > velfields(m_velocity.num_elements());
367 
368  for(int i = 0; i < m_velocity.num_elements(); ++i)
369  {
370  velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
371  }
372  stdVelocity = GetMaxStdVelocity(velfields);
373 
374  for(int el = 0; el < n_element; ++el)
375  {
376  tstep[el] = m_cflSafetyFactor /
377  (stdVelocity[el] * cLambda *
378  (ExpOrder[el]-1) * (ExpOrder[el]-1));
379  }
380 
381  NekDouble TimeStep = Vmath::Vmin(n_element, tstep, 1);
382  m_comm->AllReduce(TimeStep,LibUtilities::ReduceMin);
383 
384  return TimeStep;
385  }
386 
387 
388 
389 
391  const Array<OneD, const Array<OneD, NekDouble> > &velfield,
392  const Array<OneD, const Array<OneD, NekDouble> > &physfield,
393  Array<OneD, Array<OneD, NekDouble> > &Outarray)
394  {
395  ASSERTL1(
396  physfield.num_elements() == Outarray.num_elements(),
397  "Physfield and outarray are of different dimensions");
398 
399  int i;
400 
401  /// Number of trace points
402  int nTracePts = m_fields[0]->GetTrace()->GetNpoints();
403 
404  /// Number of spatial dimensions
405  int nDimensions = m_curl_dim;
406 
407  Array<OneD, Array<OneD, NekDouble> > traceNormals(m_curl_dim);
408 
409  for(i = 0; i < nDimensions; ++i)
410  {
411  traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
412  }
413 
414  //trace normals
415  m_fields[0]->GetTrace()->GetNormals(traceNormals);
416 
417  /// Forward state array
418  Array<OneD, NekDouble> Fwd(3*nTracePts);
419 
420  /// Backward state array
421  Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
422 
423  /// upwind numerical flux state array
424  Array<OneD, NekDouble> numflux = Bwd + nTracePts;
425 
426  /// Normal velocity array
427  Array<OneD, NekDouble> Vn (nTracePts, 0.0);
428 
429  // Extract velocity field along the trace space and multiply by trace normals
430  for(i = 0; i < nDimensions; ++i)
431  {
432  m_fields[0]->ExtractTracePhys(m_fields[m_velocity[i]]->GetPhys(), Fwd);
433  Vmath::Vvtvp(nTracePts, traceNormals[i], 1, Fwd, 1, Vn, 1, Vn, 1);
434  }
435 
436  for(i = 0; i < physfield.num_elements(); ++i)
437  {
438  /// Extract forwards/backwards trace spaces
439  /// Note: Needs to have correct i value to get boundary conditions
440  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);
441 
442  /// Upwind between elements
443  m_fields[0]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux);
444 
445  /// Construct difference between numflux and Fwd,Bwd
446  Vmath::Vsub(nTracePts, numflux, 1, Fwd, 1, Fwd, 1);
447  Vmath::Vsub(nTracePts, numflux, 1, Bwd, 1, Bwd, 1);
448 
449  /// Calculate the numerical fluxes multipling Fwd, Bwd and
450  /// numflux by the normal advection velocity
451  Vmath::Vmul(nTracePts, Fwd, 1, Vn, 1, Fwd, 1);
452  Vmath::Vmul(nTracePts, Bwd, 1, Vn, 1, Bwd, 1);
453 
454  m_fields[0]->AddFwdBwdTraceIntegral(Fwd,Bwd,Outarray[i]);
455  }
456  }
457 
458  /**
459  * Extrapolate field using equally time spaced field un,un-1,un-2, (at
460  * dt intervals) to time n+t at order Ord
461  */
463  NekDouble toff,
464  Array< OneD, Array<OneD, NekDouble> > &ExtVel)
465  {
466  int npts = m_fields[0]->GetTotPoints();
467  int nvel = m_velocity.num_elements();
468  int i,j;
469  Array<OneD, NekDouble> l(4);
470 
471  int ord = m_intSteps;
472 
473  // calculate Lagrange interpolants
474  Vmath::Fill(4,1.0,l,1);
475 
476  for(i = 0; i <= ord; ++i)
477  {
478  for(j = 0; j <= ord; ++j)
479  {
480  if(i != j)
481  {
482  l[i] *= (j*m_timestep+toff);
483  l[i] /= (j*m_timestep-i*m_timestep);
484  }
485  }
486  }
487 
488  for(i = 0; i < nvel; ++i)
489  {
490  Vmath::Smul(npts,l[0],m_previousVelFields[i],1,ExtVel[i],1);
491 
492  for(j = 1; j <= ord; ++j)
493  {
494  Blas::Daxpy(npts,l[j],m_previousVelFields[j*nvel+i],1,
495  ExtVel[i],1);
496  }
497  }
498  }
499 
500  /**
501  *
502  */
504  int HBCdata,
505  NekDouble kinvis,
506  Array<OneD, NekDouble> &Q,
507  Array<OneD, const NekDouble> &Advection)
508  {
509  Vmath::Smul(HBCdata,-kinvis,Q,1,Q,1);
510  }
511 
512  /**
513  *
514  */
516  const Array<OneD, const Array<OneD, NekDouble> > &N,
517  NekDouble Aii_Dt)
518  {
519  switch(m_velocity.num_elements())
520  {
521  case 1:
522  ASSERTL0(
523  false,
524  "Velocity correction scheme not designed to have just one velocity component");
525  break;
526  case 2:
527  AddDuDt2D(N,Aii_Dt);
528  break;
529  case 3:
530  AddDuDt3D(N,Aii_Dt);
531  break;
532  }
533  }
534 
535  /**
536  *
537  */
539  const Array<OneD, const Array<OneD, NekDouble> > &N,
540  NekDouble Aii_Dt)
541  {
542  int i,n;
543  ASSERTL0(m_velocity.num_elements() == 2," Routine currently only set up for 2D");
544 
545  int pindex=2;
546 
547  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > PBndConds;
548  Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp,UBndExp,VBndExp;
549 
550  PBndConds = m_fields[pindex]->GetBndConditions();
551  PBndExp = m_fields[pindex]->GetBndCondExpansions();
552 
553  UBndExp = m_fields[m_velocity[0]]->GetBndCondExpansions();
554  VBndExp = m_fields[m_velocity[1]]->GetBndCondExpansions();
555 
558 
559 
560  int cnt,elmtid,nq,offset, boundary,ncoeffs;
561 
562  Array<OneD, NekDouble> N1(m_pressureBCsMaxPts), N2(m_pressureBCsMaxPts);
563  Array<OneD, NekDouble> ubc(m_pressureBCsMaxPts),vbc(m_pressureBCsMaxPts);
564  Array<OneD, NekDouble> Pvals(m_pressureBCsMaxPts);
565  Array<OneD, NekDouble> Nu,Nv,Ptmp;
566 
567  for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
568  {
569  SpatialDomains::BndUserDefinedType type = PBndConds[n]->GetUserDefined();
570 
571  if(type == SpatialDomains::eHigh)
572  {
573  for(i = 0; i < PBndExp[n]->GetExpSize(); ++i,cnt++)
574  {
575  // find element and edge of this expansion.
576  elmtid = m_pressureBCtoElmtID[cnt];
577  elmt = m_fields[0]->GetExp(elmtid);
578  offset = m_fields[0]->GetPhys_Offset(elmtid);
579 
580  Nu = N[0] + offset;
581  Nv = N[1] + offset;
582 
583  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (PBndExp[n]->GetExp(i));
584  nq = Pbc->GetTotPoints();
585  ncoeffs = Pbc->GetNcoeffs();
586  boundary = m_pressureBCtoTraceID[cnt];
587 
588  // Get velocity bc
589  UBndExp[n]->GetExp(i)->BwdTrans(UBndExp[n]->GetCoeffs() + UBndExp[n]->GetCoeff_Offset(i),ubc);
590  VBndExp[n]->GetExp(i)->BwdTrans(VBndExp[n]->GetCoeffs() + VBndExp[n]->GetCoeff_Offset(i),vbc);
591 
592 
593  // Get edge values and put into Nu,Nv
594  elmt->GetEdgePhysVals(boundary,Pbc,Nu,N1);
595  elmt->GetEdgePhysVals(boundary,Pbc,Nv,N2);
596 
597 
598  // Take different as Forward Euler but N1,N2
599  // actually contain the integration of the
600  // previous steps from the time integration
601  // scheme.
602  Vmath::Vsub(nq,ubc,1,N1,1,ubc,1);
603  Vmath::Vsub(nq,vbc,1,N2,1,vbc,1);
604 
605 
606  // Divide by aii_Dt to get correct Du/Dt. This is
607  // because all coefficients in the integration
608  // scheme are normalised so u^{n+1} has unit
609  // coefficient and N is already multiplied by
610  // local coefficient when taken from integration
611  // scheme
612  Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
613  Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
614 
615  // subtrace off du/dt derivative
616  Pbc->NormVectorIProductWRTBase(ubc,vbc,Pvals);
617 
618  Vmath::Vsub(ncoeffs,Ptmp = PBndExp[n]->UpdateCoeffs()
619  +PBndExp[n]->GetCoeff_Offset(i),1, Pvals,1,
620  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1);
621  }
622  }
623  // setting if just standard BC no High order
625  {
626  cnt += PBndExp[n]->GetExpSize();
627  }
628  else
629  {
630  ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
631  }
632  }
633  }
634 
635  /**
636  *
637  */
639  const Array<OneD, const Array<OneD, NekDouble> > &N,
640  NekDouble Aii_Dt)
641  {
642  int i,n;
643  ASSERTL0(m_velocity.num_elements() == 3," Routine currently only set up for 3D");
644  int pindex=3;
645 
646  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > PBndConds;
647  Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp,UBndExp,VBndExp,WBndExp;
648 
649  PBndConds = m_fields[pindex]->GetBndConditions();
650  PBndExp = m_fields[pindex]->GetBndCondExpansions();
651 
652  UBndExp = m_fields[m_velocity[0]]->GetBndCondExpansions();
653  VBndExp = m_fields[m_velocity[1]]->GetBndCondExpansions();
654  WBndExp = m_fields[m_velocity[2]]->GetBndCondExpansions();
655 
658 
659  int cnt,elmtid,nq,offset, boundary,ncoeffs;
660 
661  Array<OneD, NekDouble> N1(m_pressureBCsMaxPts), N2(m_pressureBCsMaxPts),
663  Array<OneD, NekDouble> ubc(m_pressureBCsMaxPts),vbc(m_pressureBCsMaxPts),
664  wbc(m_pressureBCsMaxPts);
665  Array<OneD, NekDouble> Pvals(m_pressureBCsMaxPts);
666  Array<OneD, NekDouble> Nu,Nv,Nw,Ptmp;
667 
668  for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
669  {
670  SpatialDomains::BndUserDefinedType type = PBndConds[n]->GetUserDefined();
671 
672  if(type == SpatialDomains::eHigh)
673  {
674  for(i = 0; i < PBndExp[n]->GetExpSize(); ++i,cnt++)
675  {
676  // find element and face of this expansion.
677  elmtid = m_pressureBCtoElmtID[cnt];
678  elmt = m_fields[0]->GetExp(elmtid);
679  offset = m_fields[0]->GetPhys_Offset(elmtid);
680 
681  Nu = N[0] + offset;
682  Nv = N[1] + offset;
683  Nw = N[2] + offset;
684 
685  Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (PBndExp[n]->GetExp(i));
686  nq = Pbc->GetTotPoints();
687  ncoeffs = Pbc->GetNcoeffs();
688  boundary = m_pressureBCtoTraceID[cnt];
689 
690  // Get velocity bc
691  UBndExp[n]->GetExp(i)->BwdTrans(UBndExp[n]->GetCoeffs() +
692  UBndExp[n]->GetCoeff_Offset(i),ubc);
693  VBndExp[n]->GetExp(i)->BwdTrans(VBndExp[n]->GetCoeffs() +
694  VBndExp[n]->GetCoeff_Offset(i),vbc);
695  WBndExp[n]->GetExp(i)->BwdTrans(WBndExp[n]->GetCoeffs() +
696  WBndExp[n]->GetCoeff_Offset(i),wbc);
697 
698  // Get edge values and put into N1,N2,N3
699  elmt->GetFacePhysVals(boundary,Pbc,Nu,N1);
700  elmt->GetFacePhysVals(boundary,Pbc,Nv,N2);
701  elmt->GetFacePhysVals(boundary,Pbc,Nw,N3);
702 
703 
704  // Take different as Forward Euler but N1,N2,N3
705  // actually contain the integration of the
706  // previous steps from the time integration
707  // scheme.
708  Vmath::Vsub(nq,ubc,1,N1,1,ubc,1);
709  Vmath::Vsub(nq,vbc,1,N2,1,vbc,1);
710  Vmath::Vsub(nq,wbc,1,N3,1,wbc,1);
711 
712  // Divide by aii_Dt to get correct Du/Dt. This is
713  // because all coefficients in the integration
714  // scheme are normalised so u^{n+1} has unit
715  // coefficient and N is already multiplied by
716  // local coefficient when taken from integration
717  // scheme
718  Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
719  Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
720  Blas::Dscal(nq,1.0/Aii_Dt,&wbc[0],1);
721 
722  // subtrace off du/dt derivative
723  Pbc->NormVectorIProductWRTBase(ubc,vbc,wbc,Pvals);
724 
725  Vmath::Vsub(
726  ncoeffs,
727  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),
728  1,Pvals,1,
729  Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),
730  1);
731  }
732  }
733  // setting if just standard BC no High order
735  {
736  cnt += PBndExp[n]->GetExpSize();
737  }
738  else
739  {
740  ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
741  }
742  }
743  }
744 }
745