Nektar++
MMFSystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MMFSystem.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Base class for MMF systems.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
36 #include <SolverUtils/MMFSystem.h>
37 
38 namespace Nektar
39 {
40 namespace SolverUtils
41 {
42 
43 MMFSystem::MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession,
44  const SpatialDomains::MeshGraphSharedPtr &pGraph)
45  : UnsteadySystem(pSession, pGraph)
46 {
47 }
48 
49 MMFSystem::~MMFSystem()
50 {
51 }
52 
53 void MMFSystem::MMFInitObject(
54  const Array<OneD, const Array<OneD, NekDouble>> &Anisotropy,
55  const int TangentXelem)
56 {
57  m_pi = 3.14159265358979323846;
58  m_shapedim = m_fields[0]->GetShapeDimension();
59 
60  ASSERTL0(m_spacedim == 3, "Space Dimension should be 3.");
61 
62  // Factor for Numerical Flux
63  m_session->LoadParameter("alpha", m_alpha, 1.0);
64 
65  // Factor for Numerical Flux
66  m_session->LoadParameter("Incfreq", m_Incfreq, 1.0);
67 
68  // SmoothFactor
69  m_session->LoadParameter("SmoothFactor", m_SmoothFactor, 1);
70  m_session->LoadParameter("SFinit", m_SFinit, 0.0);
71 
72  // Define SurfaceType
73  if (m_session->DefinesSolverInfo("SURFACETYPE"))
74  {
75  std::string SurfaceTypeStr = m_session->GetSolverInfo("SURFACETYPE");
76  int i;
77  for (i = 0; i < (int)SIZE_SurfaceType; ++i)
78  {
79  if (boost::iequals(SurfaceTypeMap[i], SurfaceTypeStr))
80  {
81  m_surfaceType = (SurfaceType)i;
82  break;
83  }
84  }
85  }
86  else
87  {
88  m_surfaceType = (SurfaceType)0;
89  }
90 
91  // if discontinuous Galerkin determine numerical flux to use
92  for (int i = 0; i < (int)SIZE_UpwindType; ++i)
93  {
94  bool match;
95  m_session->MatchSolverInfo("UPWINDTYPE", UpwindTypeMap[i], match,
96  false);
97  if (match)
98  {
99  m_upwindType = (UpwindType)i;
100  break;
101  }
102  }
103 
104  // SetUpMovingFrames: To generate m_movingframes
105  SetUpMovingFrames(Anisotropy, TangentXelem);
106 
107  // Derive movingfrmaes at interfaces
108  // Get: m_MFtraceFwd and m_MFtraceBwd
109  ComputeMFtrace();
110 
111  // Check Movingframes and surfraceNormal
112  // Get: m_ncdotMFFwd,m_ncdotMFBwd,m_nperpcdotMFFwd,m_nperpcdotMFBwd
113  ComputencdotMF();
114 
115  // Compute nabla cdot m_movingframes
116  // Get: m_DivMF, m_CurlMF
117  ComputeDivCurlMF();
118 }
119 
120 void MMFSystem::SetUpMovingFrames(
121  const Array<OneD, const Array<OneD, NekDouble>> &Anisotropy,
122  const int TangentXelem)
123 {
124 
125  int nq = m_fields[0]->GetNpoints();
126  int MFdim = 3;
127 
128  // Construct The Moving Frames
129  m_movingframes = Array<OneD, Array<OneD, NekDouble>>(MFdim);
130 
131  for (int j = 0; j < MFdim; ++j)
132  {
133  m_movingframes[j] = Array<OneD, NekDouble>(m_spacedim * nq, 0.0);
134  }
135 
136  // Read MMF Geom Info
137  std::string conn = "TangentX";
138  m_MMFfactors = Array<OneD, NekDouble>(4);
139 
140  m_session->LoadSolverInfo("MMFDir", conn, "LOCAL");
141 
142  // (x-x_0)^2/a^2 + (y-y_0)^2/b^2 = 1
143  // factors[0] = a
144  // factors[1] = b
145  // factors[2] = x_0
146  // factors[3] = y_0
147  m_session->LoadParameter("MMFCircAxisX", m_MMFfactors[0], 1.0);
148  m_session->LoadParameter("MMFCircAxisY", m_MMFfactors[1], 1.0);
149  m_session->LoadParameter("MMFCircCentreX", m_MMFfactors[2], 0.0);
150  m_session->LoadParameter("MMFCircCentreY", m_MMFfactors[3], 0.0);
151 
152  if (conn == "TangentX")
153  m_MMFdir = SpatialDomains::eTangentX;
154  if (conn == "TangentY")
155  m_MMFdir = SpatialDomains::eTangentY;
156  if (conn == "TangentXY")
157  m_MMFdir = SpatialDomains::eTangentXY;
158  if (conn == "TangentZ")
159  m_MMFdir = SpatialDomains::eTangentZ;
160  if (conn == "TangentCircular")
161  m_MMFdir = SpatialDomains::eTangentCircular;
162  if (conn == "TangentIrregular")
163  m_MMFdir = SpatialDomains::eTangentIrregular;
164  if (conn == "TangentNonconvex")
165  m_MMFdir = SpatialDomains::eTangentNonconvex;
166  if (conn == "LOCAL")
167  m_MMFdir = SpatialDomains::eLOCAL;
168 
169  /*for (int j=0; j<m_shapedim; j++)
170  {
171  for (int k=0; k<m_spacedim; k++)
172  {
173  std::cout << "before ehsan " << nq << "\t"<< m_shapedim << "\t" <<
174  m_spacedim << "\t" <<"m_movingframes " << m_movingframes[j][k*nq] <<
175  std::endl;
176  }
177  }*/
178  // Get Tangetn vectors from GeomFactors2D, Orthonormalized = true
179  m_fields[0]->GetMovingFrames(m_MMFdir, m_MMFfactors, m_movingframes);
180  /* for (int j=0; j<m_shapedim; j++)
181  {
182  for (int k=0; k<m_spacedim; k++)
183  {
184  std::cout << "ehsan " << nq << "\t"<< m_shapedim << "\t" <<
185  m_spacedim << "\t" <<"m_movingframes " << m_movingframes[j][k*nq] <<
186  std::endl;
187  }
188  }*/
189  // Align the tangentX direction after TangentXelem
190  if (TangentXelem > 0)
191  {
192  Array<OneD, NekDouble> tmp(nq);
193  m_fields[0]->GenerateElementVector(TangentXelem, 1.0, 0.0, tmp);
194  for (int j = 0; j < m_shapedim; j++)
195  {
196  for (int k = 0; k < m_spacedim; k++)
197  {
198  Vmath::Vmul(nq, &tmp[0], 1, &m_movingframes[j][k * nq], 1,
199  &m_movingframes[j][k * nq], 1);
200  }
201  }
202 
203  m_fields[0]->GenerateElementVector(TangentXelem, 0.0, 1.0, tmp);
204  Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[0][0 * nq], 1,
205  &m_movingframes[0][0 * nq], 1);
206  Vmath::Vadd(nq, &tmp[0], 1, &m_movingframes[1][1 * nq], 1,
207  &m_movingframes[1][1 * nq], 1);
208 
209  int indxtmp = Vmath::Imax(nq, tmp, 1);
210  std::cout << "*** MF in PML Region is aligned as MF1 = ( "
211  << m_movingframes[0][indxtmp] << " , "
212  << m_movingframes[0][nq + indxtmp] << " , "
213  << m_movingframes[0][2 * nq + indxtmp] << " ) "
214  << ", MF2 = ( " << m_movingframes[1][indxtmp] << " , "
215  << m_movingframes[1][nq + indxtmp] << " , "
216  << m_movingframes[1][2 * nq + indxtmp] << " ) " << std::endl;
217  }
218 
219  // Multiply Anisotropy to movingframes
220  for (int j = 0; j < m_shapedim; ++j)
221  {
222  for (int k = 0; k < m_spacedim; ++k)
223  {
224  Vmath::Vmul(nq, &Anisotropy[j][0], 1, &m_movingframes[j][k * nq], 1,
225  &m_movingframes[j][k * nq], 1);
226  }
227  }
228 
229  // Test the moving frames
230  CheckMovingFrames(m_movingframes);
231 }
232 
233 void MMFSystem::CheckMovingFrames(
234  const Array<OneD, const Array<OneD, NekDouble>> &movingframes)
235 {
236  NekDouble t1x, t1y, t1z, t2x, t2y, t2z, t3x, t3y, t3z;
237  NekDouble dot12 = 0.0, dot23 = 0.0, dot31 = 0.0;
238  NekDouble Tol = 0.0001;
239 
240  int nq = m_fields[0]->GetNpoints();
241 
242  for (int i = 0; i < nq; ++i)
243  {
244  t1x = movingframes[0][i];
245  t1y = movingframes[0][i + nq];
246  t1z = movingframes[0][i + 2 * nq];
247 
248  t2x = movingframes[1][i];
249  t2y = movingframes[1][i + nq];
250  t2z = movingframes[1][i + 2 * nq];
251 
252  t3x = movingframes[2][i];
253  t3y = movingframes[2][i + nq];
254  t3z = movingframes[2][i + 2 * nq];
255 
256  dot12 = t1x * t2x + t1y * t2y + t1z * t2z;
257  dot23 = t2x * t3x + t2y * t3y + t2z * t3z;
258  dot31 = t3x * t1x + t3y * t1y + t3z * t1z;
259  }
260 
261  std::cout << "======================================================"
262  << std::endl;
263  std::cout << "======================================================"
264  << std::endl;
265  std::cout << "*** The first moving frame is alinged along"
266  << SpatialDomains::GeomMMFMap[m_MMFdir] << std::endl;
267 
268  Array<OneD, NekDouble> tmpx(nq), tmpy(nq), tmpz(nq);
269 
270  Vmath::Vcopy(nq, &movingframes[0][0], 1, &tmpx[0], 1);
271  Vmath::Vcopy(nq, &movingframes[0][nq], 1, &tmpy[0], 1);
272  Vmath::Vcopy(nq, &movingframes[0][2 * nq], 1, &tmpz[0], 1);
273  std::cout << nq << " , "
274  << "*** Avg MF1 = ( " << AvgAbsInt(tmpx) << " , "
275  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
276  << std::endl;
277 
278  Vmath::Vcopy(nq, &movingframes[1][0], 1, &tmpx[0], 1);
279  Vmath::Vcopy(nq, &movingframes[1][nq], 1, &tmpy[0], 1);
280  Vmath::Vcopy(nq, &movingframes[1][2 * nq], 1, &tmpz[0], 1);
281  std::cout << "*** Avg MF2 = ( " << AvgAbsInt(tmpx) << " , "
282  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
283  << std::endl;
284 
285  if (m_shapedim == 3)
286  {
287  Vmath::Vcopy(nq, &movingframes[2][0], 1, &tmpx[0], 1);
288  Vmath::Vcopy(nq, &movingframes[2][nq], 1, &tmpy[0], 1);
289  Vmath::Vcopy(nq, &movingframes[2][2 * nq], 1, &tmpz[0], 1);
290  std::cout << "*** Avg MF3 = ( " << AvgAbsInt(tmpx) << " , "
291  << AvgAbsInt(tmpy) << " , " << AvgAbsInt(tmpz) << " ) "
292  << std::endl;
293  }
294 
295  if ((fabs(dot12) + fabs(dot23) + fabs(dot31)) < Tol)
296  {
297  std::cout << "*** Moving frames are Orthogonal" << std::endl;
298  }
299 
300  else
301  {
302  std::cout << "*** Moving frames are NOT Orthogonal" << std::endl;
303  }
304 
305  Array<OneD, NekDouble> tmp;
306  for (int j = 0; j < m_shapedim; ++j)
307  {
308  tmp = Array<OneD, NekDouble>(nq, 0.0);
309  for (int k = 0; k < m_spacedim; ++k)
310  {
311  Vmath::Vvtvp(nq, &movingframes[j][k * nq], 1,
312  &movingframes[j][k * nq], 1, &tmp[0], 1, &tmp[0], 1);
313  }
314  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
315  std::cout << "*** Avg. Magnitude of MF" << j << " = " << AvgAbsInt(tmp)
316  << std::endl;
317  }
318 }
319 
320 void MMFSystem::ComputencdotMF()
321 {
322  int nq = m_fields[0]->GetNpoints();
323  int nTracePointsTot = GetTraceNpoints();
324 
325  // Compute MFjFwd and MFjBwd
326  Array<OneD, NekDouble> tmp(nq);
327 
328  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceFwd(m_shapedim);
329  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtraceBwd(m_shapedim);
330  Array<OneD, Array<OneD, NekDouble>> SurfaceNormalFwd;
331  Array<OneD, Array<OneD, NekDouble>> SurfaceNormalBwd;
332 
333  for (int j = 0; j < m_shapedim; ++j)
334  {
335  MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
336  MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
337 
338  SurfaceNormalFwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
339  SurfaceNormalBwd = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
340 
341  for (int k = 0; k < m_spacedim; ++k)
342  {
343  MFtraceFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
344  MFtraceBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
345 
346  SurfaceNormalFwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
347  SurfaceNormalBwd[k] = Array<OneD, NekDouble>(nTracePointsTot);
348 
349  Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
350 
351  m_fields[0]->GetFwdBwdTracePhys(tmp, MFtraceFwd[j][k],
352  MFtraceBwd[j][k]);
353 
354  CopyBoundaryTrace(MFtraceFwd[j][k], MFtraceBwd[j][k], eFwdEQBwd);
355  }
356  }
357 
358  VectorCrossProd(MFtraceFwd[0], MFtraceFwd[1], SurfaceNormalFwd);
359  VectorCrossProd(MFtraceBwd[0], MFtraceBwd[1], SurfaceNormalBwd);
360 
361  // Compute n \times e^i
362  m_ncdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
363  m_ncdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
364  for (int j = 0; j < m_shapedim; ++j)
365  {
366  m_ncdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
367  m_ncdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
368 
369  VectorDotProd(m_traceNormals, MFtraceFwd[j], m_ncdotMFFwd[j]);
370  VectorDotProd(m_traceNormals, MFtraceBwd[j], m_ncdotMFBwd[j]);
371  }
372 
373  // Compute n^{\perp} \times e^i
374  Array<OneD, Array<OneD, NekDouble>> SurfaceNormal(m_spacedim);
375  Array<OneD, Array<OneD, NekDouble>> Tracevector(m_spacedim);
376  for (int k = 0; k < m_spacedim; k++)
377  {
378  SurfaceNormal[k] = Array<OneD, NekDouble>(nTracePointsTot);
379  Tracevector[k] = Array<OneD, NekDouble>(nTracePointsTot);
380 
381  Vmath::Vcopy(nTracePointsTot, &m_MFtraceFwd[2][k][0], 1,
382  &SurfaceNormal[k][0], 1);
383  }
384 
385  VectorCrossProd(m_traceNormals, SurfaceNormal, Tracevector);
386 
387  m_nperpcdotMFFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
388  m_nperpcdotMFBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
389  for (int j = 0; j < m_shapedim; ++j)
390  {
391  m_nperpcdotMFFwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
392  m_nperpcdotMFBwd[j] = Array<OneD, NekDouble>(nTracePointsTot, 0.0);
393 
394  VectorDotProd(Tracevector, MFtraceFwd[j], m_nperpcdotMFFwd[j]);
395  VectorDotProd(Tracevector, MFtraceBwd[j], m_nperpcdotMFBwd[j]);
396  }
397 
398  if (m_shapedim == 2)
399  {
400  std::cout << "*** m_ncdotMFFwd = ( " << RootMeanSquare(m_ncdotMFFwd[0])
401  << " , " << RootMeanSquare(m_ncdotMFFwd[1]) << " ) "
402  << std::endl;
403  std::cout << "*** m_ncdotMFBwd = ( " << RootMeanSquare(m_ncdotMFBwd[0])
404  << " , " << RootMeanSquare(m_ncdotMFBwd[1]) << " ) "
405  << std::endl;
406 
407  std::cout << "*** m_nperpcdotMFFwd = ( "
408  << RootMeanSquare(m_nperpcdotMFFwd[0]) << " , "
409  << RootMeanSquare(m_nperpcdotMFFwd[1]) << " ) " << std::endl;
410  std::cout << "*** m_nperpcdotMFBwd = ( "
411  << RootMeanSquare(m_nperpcdotMFBwd[0]) << " , "
412  << RootMeanSquare(m_nperpcdotMFBwd[1]) << " ) " << std::endl;
413  }
414 
415  else if (m_shapedim == 3)
416  {
417  std::cout << "*** m_ncdotMFFwd = ( "
418  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[0], 1) << " , "
419  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[1], 1) << " , "
420  << Vmath::Vsum(nTracePointsTot, m_ncdotMFFwd[2], 1) << " ) "
421  << std::endl;
422  std::cout << "*** m_ncdotMFBwd = ( "
423  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[0], 1) << " , "
424  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[1], 1) << " , "
425  << Vmath::Vsum(nTracePointsTot, m_ncdotMFBwd[2], 1) << " ) "
426  << std::endl;
427 
428  std::cout << "*** m_nperpcdotMFFwd = ( "
429  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[0], 1)
430  << " , "
431  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[1], 1)
432  << " , "
433  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFFwd[2], 1)
434  << " ) " << std::endl;
435  std::cout << "*** m_nperpcdotMFBwd = ( "
436  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[0], 1)
437  << " , "
438  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[1], 1)
439  << " , "
440  << Vmath::Vsum(nTracePointsTot, m_nperpcdotMFBwd[2], 1)
441  << " ) " << std::endl;
442  }
443 }
444 
445 void MMFSystem::ComputeDivCurlMF()
446 {
447  int nq = m_fields[0]->GetNpoints();
448  int MMFdim = 3;
449 
450  Array<OneD, NekDouble> tmp(nq);
451  Array<OneD, NekDouble> Dtmp(nq);
452 
453  m_DivMF = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
454  for (int j = 0; j < MMFdim; ++j)
455  {
456  m_DivMF[j] = Array<OneD, NekDouble>(nq, 0.0);
457  for (int k = 0; k < m_spacedim; ++k)
458  {
459  Vmath::Vcopy(nq, &m_movingframes[j][k * nq], 1, &tmp[0], 1);
460 
461  m_fields[0]->PhysDeriv(k, tmp, Dtmp);
462  Vmath::Vadd(nq, &Dtmp[0], 1, &m_DivMF[j][0], 1, &m_DivMF[j][0], 1);
463  }
464  }
465 
466  std::cout << "*** Divergence of MF1 = " << AvgInt(m_DivMF[0])
467  << ", MF2 = " << AvgInt(m_DivMF[1])
468  << ", MF3 = " << AvgInt(m_DivMF[2]) << std::endl;
469 
470  // Compute Curl of MF: CurlMF[i][j] = (\nabla \times e^i) cdot e^j
471  m_CurlMF = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MMFdim);
472  for (int i = 0; i < MMFdim; ++i)
473  {
474  m_CurlMF[i] = Array<OneD, Array<OneD, NekDouble>>(MMFdim);
475 
476  for (int j = 0; j < MMFdim; ++j)
477  {
478  m_CurlMF[i][j] = Array<OneD, NekDouble>(nq, 0.0);
479  }
480  }
481 
482  Array<OneD, Array<OneD, NekDouble>> MFtmp(m_spacedim);
483  Array<OneD, Array<OneD, NekDouble>> CurlMFtmp(m_spacedim);
484 
485  for (int i = 0; i < m_spacedim; ++i)
486  {
487  MFtmp[i] = Array<OneD, NekDouble>(nq);
488  CurlMFtmp[i] = Array<OneD, NekDouble>(nq);
489  }
490 
491  for (int dir = 0; dir < MMFdim; dir++)
492  {
493  for (int i = 0; i < m_spacedim; ++i)
494  {
495  Vmath::Vcopy(nq, &m_movingframes[dir][i * nq], 1, &MFtmp[i][0], 1);
496  }
497 
498  ComputeCurl(MFtmp, CurlMFtmp);
499 
500  for (int j = 0; j < MMFdim; ++j)
501  {
502  for (int i = 0; i < m_spacedim; ++i)
503  {
504  Vmath::Vvtvp(nq, &m_movingframes[j][i * nq], 1,
505  &CurlMFtmp[i][0], 1, &m_CurlMF[dir][j][0], 1,
506  &m_CurlMF[dir][j][0], 1);
507  }
508  }
509  }
510 
511  std::cout << "*** Curl of MF1 = ( " << AvgInt(m_CurlMF[0][0]) << " , "
512  << AvgInt(m_CurlMF[0][1]) << " , " << AvgInt(m_CurlMF[0][2])
513  << " ) " << std::endl;
514  std::cout << "*** Curl of MF2 = ( " << AvgInt(m_CurlMF[1][0]) << " , "
515  << AvgInt(m_CurlMF[1][1]) << " , " << AvgInt(m_CurlMF[1][2])
516  << " ) " << std::endl;
517  std::cout << "*** Curl of MF3 = ( " << AvgInt(m_CurlMF[2][0]) << " , "
518  << AvgInt(m_CurlMF[2][1]) << " , " << AvgInt(m_CurlMF[2][2])
519  << " ) " << std::endl;
520 }
521 
522 void MMFSystem::ComputeMFtrace()
523 {
524  int MFdim = 3;
525 
526  int nq = m_fields[0]->GetNpoints();
527  int nTraceNumPoints = GetTraceTotPoints();
528 
529  Array<OneD, NekDouble> tmp(nq);
530  Array<OneD, NekDouble> Fwdtmp(nq);
531  Array<OneD, NekDouble> Bwdtmp(nq);
532 
533  m_MFtraceFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
534  m_MFtraceBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
535 
536  for (int j = 0; j < MFdim; ++j)
537  {
538  m_MFtraceFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
539  m_MFtraceBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
540  }
541 
542  // m_MFtraceFwd[0] = e^1_{Fwd}, m_MFtraceFwd[1] = e^2_{Fwd}
543  for (int j = 0; j < MFdim; ++j)
544  {
545  for (int i = 0; i < m_spacedim; ++i)
546  {
547  m_MFtraceFwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
548  m_MFtraceBwd[j][i] = Array<OneD, NekDouble>(nTraceNumPoints);
549 
550  Vmath::Vcopy(nq, &m_movingframes[j][i * nq], 1, &tmp[0], 1);
551 
552  m_fields[0]->GetFwdBwdTracePhys(tmp, Fwdtmp, Bwdtmp);
553 
554  CopyBoundaryTrace(Fwdtmp, Bwdtmp, eFwdEQBwd);
555 
556  Vmath::Vcopy(nTraceNumPoints, &Fwdtmp[0], 1, &m_MFtraceFwd[j][i][0],
557  1);
558  Vmath::Vcopy(nTraceNumPoints, &Bwdtmp[0], 1, &m_MFtraceBwd[j][i][0],
559  1);
560  }
561  }
562 
563  std::cout << "*** MFtraceFwd = ( " << VectorAvgMagnitude(m_MFtraceFwd[0])
564  << " , " << VectorAvgMagnitude(m_MFtraceFwd[1]) << " , "
565  << VectorAvgMagnitude(m_MFtraceFwd[2]) << " ) " << std::endl;
566  std::cout << "*** MFtraceBwd = ( " << VectorAvgMagnitude(m_MFtraceBwd[0])
567  << " , " << VectorAvgMagnitude(m_MFtraceBwd[1]) << " , "
568  << VectorAvgMagnitude(m_MFtraceBwd[2]) << " ) " << std::endl;
569 }
570 
571 void MMFSystem::DeriveCrossProductMF(
572  Array<OneD, Array<OneD, NekDouble>> &CrossProductMF)
573 {
574  int MFdim = 3;
575  int nq = GetTotPoints();
576 
577  CrossProductMF = Array<OneD, Array<OneD, NekDouble>>(MFdim);
578 
579  Array<OneD, Array<OneD, NekDouble>> MF1tmp(m_spacedim);
580  Array<OneD, Array<OneD, NekDouble>> MF2tmp(m_spacedim);
581  Array<OneD, Array<OneD, NekDouble>> MF3tmp(m_spacedim);
582  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> MFtmpCurl(MFdim);
583  for (int j = 0; j < MFdim; ++j)
584  {
585  CrossProductMF[j] = Array<OneD, NekDouble>(nq * m_spacedim);
586  MFtmpCurl[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
587  for (int k = 0; k < m_spacedim; ++k)
588  {
589  MFtmpCurl[j][k] = Array<OneD, NekDouble>(nq);
590  }
591  }
592 
593  for (int k = 0; k < m_spacedim; ++k)
594  {
595  MF1tmp[k] = Array<OneD, NekDouble>(nq);
596  MF2tmp[k] = Array<OneD, NekDouble>(nq);
597  MF3tmp[k] = Array<OneD, NekDouble>(nq);
598 
599  Vmath::Vcopy(nq, &m_movingframes[0][k * nq], 1, &MF1tmp[k][0], 1);
600  Vmath::Vcopy(nq, &m_movingframes[1][k * nq], 1, &MF2tmp[k][0], 1);
601  Vmath::Vcopy(nq, &m_movingframes[2][k * nq], 1, &MF3tmp[k][0], 1);
602  }
603 
604  VectorCrossProd(MF3tmp, MF1tmp, MFtmpCurl[0]);
605  VectorCrossProd(MF2tmp, MF3tmp, MFtmpCurl[1]);
606  VectorCrossProd(MF1tmp, MF2tmp, MFtmpCurl[2]);
607 
608  for (int j = 0; j < MFdim; ++j)
609  {
610  for (int k = 0; k < m_spacedim; ++k)
611  {
612  Vmath::Vcopy(nq, &MFtmpCurl[j][k][0], 1, &CrossProductMF[j][k * nq],
613  1);
614  }
615  }
616 }
617 
618 void MMFSystem::ComputeNtimesMF()
619 {
620  int MFdim = 3;
621  int nTracePointsTot = GetTraceNpoints();
622 
623  m_ntimesMFFwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
624  m_ntimesMFBwd = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
625  m_ntimes_ntimesMFFwd =
626  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
627  m_ntimes_ntimesMFBwd =
628  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
629  for (int j = 0; j < MFdim; ++j)
630  {
631  m_ntimesMFFwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
632  m_ntimesMFBwd[j] = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
633  m_ntimes_ntimesMFFwd[j] =
634  Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
635  m_ntimes_ntimesMFBwd[j] =
636  Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
637 
638  for (int k = 0; k < m_spacedim; ++k)
639  {
640  m_ntimesMFFwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
641  m_ntimesMFBwd[j][k] = Array<OneD, NekDouble>(nTracePointsTot);
642  m_ntimes_ntimesMFFwd[j][k] =
643  Array<OneD, NekDouble>(nTracePointsTot);
644  m_ntimes_ntimesMFBwd[j][k] =
645  Array<OneD, NekDouble>(nTracePointsTot);
646  }
647 
648  VectorCrossProd(m_traceNormals, m_MFtraceFwd[j], m_ntimesMFFwd[j]);
649  VectorCrossProd(m_traceNormals, m_MFtraceBwd[j], m_ntimesMFBwd[j]);
650  VectorCrossProd(m_traceNormals, m_ntimesMFFwd[j],
651  m_ntimes_ntimesMFFwd[j]);
652  VectorCrossProd(m_traceNormals, m_ntimesMFBwd[j],
653  m_ntimes_ntimesMFBwd[j]);
654  }
655 
656  std::cout << "*** m_ntimesMFFwd = ( "
657  << VectorAvgMagnitude(m_ntimesMFFwd[0]) << " , "
658  << VectorAvgMagnitude(m_ntimesMFFwd[1]) << " , "
659  << VectorAvgMagnitude(m_ntimesMFFwd[2]) << " ) " << std::endl;
660  std::cout << "*** m_ntimesMFBwd = ( "
661  << VectorAvgMagnitude(m_ntimesMFBwd[0]) << " , "
662  << VectorAvgMagnitude(m_ntimesMFBwd[1]) << " , "
663  << VectorAvgMagnitude(m_ntimesMFBwd[2]) << " ) " << std::endl;
664  std::cout << "*** m_ntimes_ntimesMFFwd = ( "
665  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[0]) << " , "
666  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[1]) << " , "
667  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[2]) << " ) "
668  << std::endl;
669  std::cout << "*** m_ntimes_ntimesMFBwd = ( "
670  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[0]) << " , "
671  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[1]) << " , "
672  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[2]) << " ) "
673  << std::endl;
674 }
675 
676 void MMFSystem::VectorDotProd(
677  const Array<OneD, const Array<OneD, NekDouble>> &v1,
678  const Array<OneD, const Array<OneD, NekDouble>> &v2,
679  Array<OneD, NekDouble> &v3)
680 {
681  int coordim = v1.size();
682  int nq = v1[0].size();
683 
684  v3 = Array<OneD, NekDouble>(nq, 0.0);
685  for (int i = 0; i < coordim; ++i)
686  {
687  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &v3[0], 1, &v3[0], 1);
688  }
689 }
690 
691 /**
692  * Computes the vector cross-product in 3D of \a v1 and \a v2, storing
693  * the result in \a v3.
694  * @param v1 First input vector.
695  * @param v2 Second input vector.
696  * @param v3 Output vector computed to be orthogonal to
697  * both \a v1 and \a v2.
698  */
699 void MMFSystem::VectorCrossProd(
700  const Array<OneD, const Array<OneD, NekDouble>> &v1,
701  const Array<OneD, const Array<OneD, NekDouble>> &v2,
702  Array<OneD, Array<OneD, NekDouble>> &v3)
703 {
704  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
705  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
706  ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
707 
708  int nq = v1[0].size();
709  Array<OneD, NekDouble> temp(nq);
710 
711  Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
712  Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
713 
714  Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
715  Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
716 
717  Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
718  Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
719 }
720 
721 void MMFSystem::VectorCrossProd(const Array<OneD, NekDouble> &v1,
722  const Array<OneD, NekDouble> &v2,
723  Array<OneD, NekDouble> &v3)
724 {
725  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
726  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
727  ASSERTL0(v3.size() == 3, "Output vector has dimension not equal to 3.");
728 
729  v3[0] = v1[1] * v2[2] - v1[2] * v2[1];
730  v3[1] = v1[2] * v2[0] - v1[0] * v2[2];
731  v3[2] = v1[0] * v2[1] - v1[1] * v2[0];
732 }
733 
734 void MMFSystem::ComputeCurl(
735  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
736  Array<OneD, Array<OneD, NekDouble>> &outarray)
737 
738 {
739  int nq = inarray[0].size();
740 
741  Array<OneD, NekDouble> tmpx, tmpy, tmpz;
742  Array<OneD, NekDouble> Dtmpzdx, Dtmpydx, Dtmpxdy, Dtmpzdy, Dtmpxdz, Dtmpydz;
743 
744  tmpx = Array<OneD, NekDouble>(nq);
745  tmpy = Array<OneD, NekDouble>(nq);
746  tmpz = Array<OneD, NekDouble>(nq);
747 
748  Dtmpzdx = Array<OneD, NekDouble>(nq);
749  Dtmpydx = Array<OneD, NekDouble>(nq);
750  Dtmpxdy = Array<OneD, NekDouble>(nq);
751  Dtmpzdy = Array<OneD, NekDouble>(nq);
752  Dtmpxdz = Array<OneD, NekDouble>(nq);
753  Dtmpydz = Array<OneD, NekDouble>(nq);
754 
755  for (int k = 0; k < m_spacedim; ++k)
756  {
757  Vmath::Vcopy(nq, &inarray[0][0], 1, &tmpx[0], 1);
758  Vmath::Vcopy(nq, &inarray[1][0], 1, &tmpy[0], 1);
759  Vmath::Vcopy(nq, &inarray[2][0], 1, &tmpz[0], 1);
760 
761  m_fields[0]->PhysDeriv(0, tmpz, Dtmpzdx);
762  m_fields[0]->PhysDeriv(0, tmpy, Dtmpydx);
763  m_fields[0]->PhysDeriv(1, tmpx, Dtmpxdy);
764  m_fields[0]->PhysDeriv(1, tmpz, Dtmpzdy);
765  m_fields[0]->PhysDeriv(2, tmpx, Dtmpxdz);
766  m_fields[0]->PhysDeriv(2, tmpy, Dtmpydz);
767 
768  Vmath::Vsub(nq, &Dtmpzdy[0], 1, &Dtmpydz[0], 1, &outarray[0][0], 1);
769  Vmath::Vsub(nq, &Dtmpxdz[0], 1, &Dtmpzdx[0], 1, &outarray[1][0], 1);
770  Vmath::Vsub(nq, &Dtmpydx[0], 1, &Dtmpxdy[0], 1, &outarray[2][0], 1);
771  }
772 }
773 
774 Array<OneD, NekDouble> MMFSystem::CartesianToMovingframes(
775  const Array<OneD, const Array<OneD, NekDouble>> &uvec, unsigned int field)
776 {
777  int nq = m_fields[0]->GetNpoints();
778 
779  Array<OneD, NekDouble> outarray(nq, 0.0);
780 
781  // new u0 = ( [u v] \cdot e^1 )/ |e^1|^2
782  Vmath::Vmul(nq, &m_movingframes[field][0], 1, &uvec[0][0], 1, &outarray[0],
783  1);
784  Vmath::Vvtvp(nq, &m_movingframes[field][nq], 1, &uvec[1][0], 1,
785  &outarray[0], 1, &outarray[0], 1);
786  Vmath::Vvtvp(nq, &m_movingframes[field][2 * nq], 1, &uvec[2][0], 1,
787  &outarray[0], 1, &outarray[0], 1);
788 
789  return outarray;
790 }
791 
792 // x = r \cos \theta \cos \varphi
793 // y = r \cos \theta \sin \varphi
794 // z = r \sin \theta
795 void MMFSystem::CartesianToSpherical(const NekDouble x0j, const NekDouble x1j,
796  const NekDouble x2j, NekDouble &sin_varphi,
797  NekDouble &cos_varphi,
798  NekDouble &sin_theta, NekDouble &cos_theta)
799 {
800  NekDouble radius;
801  NekDouble radxy;
802  NekDouble Tol = 0.0000001;
803 
804  NekDouble m_Xscale = 1.0;
805  NekDouble m_Yscale = 1.0;
806  NekDouble m_Zscale = 1.0;
807 
808  radius = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
809  x1j * x1j / (m_Yscale * m_Yscale) +
810  x2j * x2j / (m_Zscale * m_Zscale));
811  radxy = sqrt(x0j * x0j / (m_Xscale * m_Xscale) +
812  x1j * x1j / (m_Yscale * m_Yscale));
813 
814  if (radxy > Tol)
815  {
816  sin_varphi = x1j / (radxy * m_Yscale);
817  cos_varphi = x0j / (radxy * m_Xscale);
818  }
819 
820  else
821  {
822  sin_varphi = 0.0;
823  if (x2j > 0)
824  {
825  cos_varphi = 1.0;
826  }
827 
828  else
829  {
830  cos_varphi = -1.0;
831  }
832  }
833 
834  sin_theta = x2j / (radius * m_Zscale);
835  cos_theta = radxy / radius;
836 }
837 
838 void MMFSystem::CopyBoundaryTrace(const Array<OneD, const NekDouble> &Fwd,
839  Array<OneD, NekDouble> &Bwd,
840  const BoundaryCopyType BDCopyType,
841  const int var, const std::string BDtype)
842 {
843  int id1, id2, npts, nptselem, cnt = 0, bdrycnt = 0;
844  Array<OneD, NekDouble> Dirichlet, x0, x1, x2;
845 
846  // loop over Boundary Regions
847  for (int n = 0; n < m_fields[var]->GetBndConditions().size(); ++n)
848  {
849  nptselem = m_fields[var]->GetBndCondExpansions()[n]->GetNpoints();
850 
851  Dirichlet = Array<OneD, NekDouble>(nptselem);
852  x0 = Array<OneD, NekDouble>(nptselem);
853  x1 = Array<OneD, NekDouble>(nptselem);
854  x2 = Array<OneD, NekDouble>(nptselem);
855 
856  if (BDCopyType == eDirichlet)
857  {
858  m_fields[var]->GetBndCondExpansions()[n]->GetCoords(x0, x1, x2);
859  LibUtilities::EquationSharedPtr ifunc =
860  m_session->GetFunction("BoundaryConditions", 0);
861  ifunc->Evaluate(x0, x1, x2, 0.0, Dirichlet);
862  }
863 
864  for (int e = 0;
865  e < m_fields[var]->GetBndCondExpansions()[n]->GetExpSize(); ++e)
866  {
867  npts = m_fields[var]
868  ->GetBndCondExpansions()[n]
869  ->GetExp(e)
870  ->GetNumPoints(0);
871  id1 = m_fields[var]->GetBndCondExpansions()[n]->GetPhys_Offset(e);
872  id2 = m_fields[var]->GetTrace()->GetPhys_Offset(
873  m_fields[var]->GetTraceMap()->GetBndCondIDToGlobalTraceID(cnt +
874  e));
875 
876  if (m_fields[var]->GetBndConditions()[n]->GetUserDefined() ==
877  BDtype ||
878  BDtype == "NoUserDefined")
879  {
880  switch (BDCopyType)
881  {
882  case eDirichlet:
883  {
884  Vmath::Vcopy(npts, &Dirichlet[id1], 1, &Bwd[id2], 1);
885  bdrycnt++;
886  }
887  break;
888 
889  case eFwdEQBwd:
890  {
891  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
892  bdrycnt++;
893  }
894  break;
895 
896  case eFwdEQNegBwd:
897  {
898  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
899  Vmath::Neg(npts, &Bwd[id2], 1);
900  bdrycnt++;
901  }
902  break;
903 
904  default:
905  break;
906  }
907  }
908  }
909 
910  cnt += m_fields[var]->GetBndCondExpansions()[n]->GetExpSize();
911  }
912 }
913 
914 // Compute (e^[dir] \cdot ( n \times e^3)) * (imFwd * Fwd^3 + imBwd * Bwd^3 )
915 void MMFSystem::ComputeNtimesFz(const int dir,
916  const Array<OneD, Array<OneD, NekDouble>> &Fwd,
917  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
918  const Array<OneD, const NekDouble> &imFwd,
919  const Array<OneD, const NekDouble> &imBwd,
920  Array<OneD, NekDouble> &outarrayFwd,
921  Array<OneD, NekDouble> &outarrayBwd)
922 
923 {
924  int nTraceNumPoints = GetTraceTotPoints();
925 
926  NekDouble tmpFwd, tmpBwd;
927  NekDouble Aver, ntimesz;
928 
929  for (int i = 0; i < nTraceNumPoints; ++i)
930  {
931  Aver = 0.5 * (imFwd[i] + imBwd[i]);
932 
933  tmpFwd = 0.0;
934  tmpBwd = 0.0;
935  for (int k = 0; k < m_spacedim; ++k)
936  {
937  ntimesz = 0.5 * (imFwd[i] * Fwd[2][i] + imBwd[i] * Bwd[2][i]);
938 
939  tmpFwd +=
940  m_MFtraceFwd[dir][k][i] * m_ntimesMFFwd[2][k][i] * ntimesz;
941  tmpBwd +=
942  m_MFtraceBwd[dir][k][i] * m_ntimesMFBwd[2][k][i] * ntimesz;
943  }
944 
945  outarrayFwd[i] = tmpFwd / Aver;
946  outarrayBwd[i] = tmpBwd / Aver;
947  }
948 }
949 
950 // Compute e^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd EBwd ) ) / 2{{Y_i}}
951 void MMFSystem::ComputeNtimesF12(const Array<OneD, Array<OneD, NekDouble>> &Fwd,
952  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
953  const Array<OneD, const NekDouble> &im1Fwd,
954  const Array<OneD, const NekDouble> &im1Bwd,
955  const Array<OneD, const NekDouble> &im2Fwd,
956  const Array<OneD, const NekDouble> &im2Bwd,
957  Array<OneD, NekDouble> &outarrayFwd,
958  Array<OneD, NekDouble> &outarrayBwd)
959 {
960  int nTraceNumPoints = GetTraceTotPoints();
961 
962  NekDouble tmpFwd, tmpBwd, Aver1, Aver2, HFwdk, HBwdk;
963 
964  Array<OneD, NekDouble> z1HAver(m_spacedim);
965  Array<OneD, NekDouble> z2HAver(m_spacedim);
966  Array<OneD, NekDouble> n1e1(m_spacedim);
967  Array<OneD, NekDouble> n2e2(m_spacedim);
968 
969  Array<OneD, NekDouble> n1e1_times_z1HAver(m_spacedim);
970  Array<OneD, NekDouble> n2e2_times_z2HAver(m_spacedim);
971 
972  for (int i = 0; i < nTraceNumPoints; ++i)
973  {
974  Aver1 = 0.5 * (im1Fwd[i] + im1Bwd[i]);
975  Aver2 = 0.5 * (im2Fwd[i] + im2Bwd[i]);
976 
977  for (int k = 0; k < m_spacedim; k++)
978  {
979  // Compute \vec{HFwd} and \vec{HBwd}
980  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
981  Fwd[1][i] * m_MFtraceFwd[1][k][i];
982  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
983  Bwd[1][i] * m_MFtraceBwd[1][k][i];
984 
985  // Compute z_i {{ \vec{H} }}
986  z1HAver[k] = 0.5 * (im1Fwd[i] * HFwdk + im1Bwd[i] * HBwdk);
987  z2HAver[k] = 0.5 * (im2Fwd[i] * HFwdk + im2Bwd[i] * HBwdk);
988 
989  // Choose e^i for the one in anisotropy region
990  n1e1[k] = m_ncdotMFFwd[0][i] * m_MFtraceFwd[0][k][i];
991  n2e2[k] = m_ncdotMFFwd[1][i] * m_MFtraceFwd[1][k][i];
992  }
993 
994  // Compute n1e1 \times z1HAver and n2e2 \times z2HAver
995  VectorCrossProd(n1e1, z1HAver, n1e1_times_z1HAver);
996  VectorCrossProd(n2e2, z2HAver, n2e2_times_z2HAver);
997 
998  // e^3 \cdot ( n1e1 \times z1HAver + n2e2 \times z2HAver)
999  tmpFwd = 0.0;
1000  tmpBwd = 0.0;
1001  for (int k = 0; k < m_spacedim; k++)
1002  {
1003  tmpFwd += m_MFtraceFwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1004  n2e2_times_z2HAver[k] / Aver2);
1005  tmpBwd += m_MFtraceBwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1006  n2e2_times_z2HAver[k] / Aver2);
1007  }
1008 
1009  outarrayFwd[i] = tmpFwd;
1010  outarrayBwd[i] = tmpBwd;
1011  }
1012 }
1013 
1014 void MMFSystem::ComputeNtimestimesdFz(
1015  const int dir, const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1016  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1017  const Array<OneD, const NekDouble> &imFwd,
1018  const Array<OneD, const NekDouble> &imBwd,
1019  Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1020 {
1021  int nTraceNumPoints = GetTraceTotPoints();
1022 
1023  Array<OneD, NekDouble> dH(m_spacedim);
1024  Array<OneD, NekDouble> nFwd(m_spacedim);
1025 
1026  Array<OneD, NekDouble> eiFwd(m_spacedim);
1027  Array<OneD, NekDouble> eiBwd(m_spacedim);
1028 
1029  Array<OneD, NekDouble> eitimesdHFwd(m_spacedim);
1030  Array<OneD, NekDouble> eitimesdHBwd(m_spacedim);
1031 
1032  Array<OneD, NekDouble> ntimeseitimesdHFwd(m_spacedim);
1033  Array<OneD, NekDouble> ntimeseitimesdHBwd(m_spacedim);
1034 
1035  NekDouble Aver, HFwdk, HBwdk, tmpFwd, tmpBwd;
1036  for (int i = 0; i < nTraceNumPoints; ++i)
1037  {
1038  Aver = 0.5 * (imFwd[i] + imBwd[i]);
1039 
1040  // Get [H]
1041  for (int k = 0; k < m_spacedim; k++)
1042  {
1043  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
1044  Fwd[1][i] * m_MFtraceFwd[1][k][i];
1045  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
1046  Bwd[1][i] * m_MFtraceBwd[1][k][i];
1047  dH[k] = HFwdk - HBwdk;
1048 
1049  eiFwd[k] = m_MFtraceFwd[dir][k][i];
1050  eiBwd[k] = m_MFtraceBwd[dir][k][i];
1051 
1052  nFwd[k] = m_traceNormals[k][i];
1053  }
1054 
1055  // MFtraceFwd (MFtraceBwd) \times [H]
1056  // VectorCrossProd(eiFwd, dH, eitimesdHFwd);
1057  // VectorCrossProd(eiBwd, dH, eitimesdHBwd);
1058  VectorCrossProd(nFwd, dH, eitimesdHFwd);
1059  VectorCrossProd(nFwd, dH, eitimesdHBwd);
1060 
1061  // n times eitimesdH
1062  VectorCrossProd(nFwd, eitimesdHFwd, ntimeseitimesdHFwd);
1063  VectorCrossProd(nFwd, eitimesdHBwd, ntimeseitimesdHBwd);
1064 
1065  // MFtraceFwd \cdot ntimeseitimesdH
1066  tmpFwd = 0.0;
1067  tmpBwd = 0.0;
1068  for (int k = 0; k < m_spacedim; k++)
1069  {
1070  tmpFwd += eiFwd[k] * ntimeseitimesdHFwd[k];
1071  tmpBwd += eiBwd[k] * ntimeseitimesdHBwd[k];
1072  }
1073 
1074  outarrayFwd[i] = 0.5 * m_alpha * tmpFwd / Aver;
1075  outarrayBwd[i] = 0.5 * m_alpha * tmpBwd / Aver;
1076  }
1077 }
1078 
1079 // Compute - \alpha [E3] / ( {{im1}} + {{im2}} )
1080 void MMFSystem::ComputeNtimestimesdF12(
1081  const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1082  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1083  const Array<OneD, const NekDouble> &im1Fwd,
1084  const Array<OneD, const NekDouble> &im1Bwd,
1085  const Array<OneD, const NekDouble> &im2Fwd,
1086  const Array<OneD, const NekDouble> &im2Bwd,
1087  Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1088 {
1089  int nTraceNumPoints = GetTraceTotPoints();
1090 
1091  Array<OneD, NekDouble> directFwd(nTraceNumPoints);
1092  Array<OneD, NekDouble> directBwd(nTraceNumPoints);
1093 
1094  NekDouble Aver1, Aver2;
1095  for (int i = 0; i < nTraceNumPoints; ++i)
1096  {
1097  Aver1 = im1Fwd[i] + im1Bwd[i];
1098  Aver2 = im2Fwd[i] + im2Bwd[i];
1099 
1100  outarrayFwd[i] =
1101  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1102  outarrayBwd[i] =
1103  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1104  }
1105 }
1106 
1107 // dim = dimension, pol = polarization, 0 = TM, 1 = TE.
1108 void MMFSystem::ComputeZimYim(Array<OneD, Array<OneD, NekDouble>> &epsvec,
1109  Array<OneD, Array<OneD, NekDouble>> &muvec)
1110 
1111 {
1112  int nTraceNumPoints = GetTraceNpoints();
1113 
1114  switch (m_TestMaxwellType)
1115  {
1116  case eMaxwell1D:
1117  case eScatField1D:
1118  {
1119  Array<OneD, NekDouble> Fwdeps(nTraceNumPoints, 1.0);
1120  Array<OneD, NekDouble> Bwdeps(nTraceNumPoints, 1.0);
1121  Array<OneD, NekDouble> Fwdmu(nTraceNumPoints, 1.0);
1122  Array<OneD, NekDouble> Bwdmu(nTraceNumPoints, 1.0);
1123  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps, Bwdeps);
1124  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdeps, Bwdeps);
1125 
1126  CopyBoundaryTrace(Fwdeps, Bwdeps, eFwdEQBwd, 0);
1127  CopyBoundaryTrace(Fwdmu, Bwdmu, eFwdEQBwd, 1);
1128 
1129  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1130  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1131  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1132  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1133 
1134  m_ZimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1135  m_ZimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1136  m_YimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1137  m_YimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1138 
1139  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd / epsBwd)
1140  for (int i = 0; i < nTraceNumPoints; ++i)
1141  {
1142  m_ZimFwd[0][i] = sqrt(Fwdmu[i] / Fwdeps[i]);
1143  m_ZimBwd[0][i] = sqrt(Bwdmu[i] / Bwdeps[i]);
1144 
1145  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1146  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1147  }
1148 
1149  std::cout << "*** ZimFwd = " << RootMeanSquare(m_ZimFwd[0])
1150  << ", ZimBwd = " << RootMeanSquare(m_ZimBwd[0])
1151  << ", YimFwd = " << RootMeanSquare(m_YimFwd[0])
1152  << ", YimBwd = " << RootMeanSquare(m_YimBwd[0])
1153  << std::endl;
1154  }
1155  break;
1156 
1157  case eTestMaxwell2DPEC:
1158  case eTestMaxwell2DPECAVGFLUX:
1159  case eTestMaxwell2DPMC:
1160  case eScatField2D:
1161  case eTotField2D:
1162  case eMaxwellSphere:
1163  case eELF2DSurface:
1164  {
1165  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1166  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1167  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1168  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1169 
1170  for (int j = 0; j < m_shapedim; ++j)
1171  {
1172  m_ZimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1173  m_ZimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1174  m_YimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1175  m_YimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1176  }
1177 
1178  switch (m_PolType)
1179  {
1180  case eTransMagnetic:
1181  {
1182  Array<OneD, NekDouble> Fwdmu1(nTraceNumPoints);
1183  Array<OneD, NekDouble> Bwdmu1(nTraceNumPoints);
1184  Array<OneD, NekDouble> Fwdmu2(nTraceNumPoints);
1185  Array<OneD, NekDouble> Bwdmu2(nTraceNumPoints);
1186  Array<OneD, NekDouble> Fwdeps3(nTraceNumPoints);
1187  Array<OneD, NekDouble> Bwdeps3(nTraceNumPoints);
1188 
1189  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdmu1, Bwdmu1);
1190  m_fields[0]->GetFwdBwdTracePhys(muvec[1], Fwdmu2, Bwdmu2);
1191  m_fields[0]->GetFwdBwdTracePhys(epsvec[2], Fwdeps3,
1192  Bwdeps3);
1193 
1194  CopyBoundaryTrace(Fwdmu1, Bwdmu1, eFwdEQBwd, 0);
1195  CopyBoundaryTrace(Fwdmu2, Bwdmu2, eFwdEQBwd, 1);
1196  CopyBoundaryTrace(Fwdeps3, Bwdeps3, eFwdEQBwd, 2);
1197 
1198  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd /
1199  // epsBwd)
1200  for (int i = 0; i < nTraceNumPoints; ++i)
1201  {
1202  m_ZimFwd[0][i] = sqrt(Fwdmu2[i] / Fwdeps3[i]);
1203  m_ZimBwd[0][i] = sqrt(Bwdmu2[i] / Bwdeps3[i]);
1204 
1205  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1206  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1207 
1208  m_ZimFwd[1][i] = sqrt(Fwdmu1[i] / Fwdeps3[i]);
1209  m_ZimBwd[1][i] = sqrt(Bwdmu1[i] / Bwdeps3[i]);
1210 
1211  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1212  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1213  }
1214  }
1215  break; // eTransMagnetic
1216 
1217  case eTransElectric:
1218  {
1219  Array<OneD, NekDouble> Fwdeps1(nTraceNumPoints);
1220  Array<OneD, NekDouble> Bwdeps1(nTraceNumPoints);
1221  Array<OneD, NekDouble> Fwdeps2(nTraceNumPoints);
1222  Array<OneD, NekDouble> Bwdeps2(nTraceNumPoints);
1223  Array<OneD, NekDouble> Fwdmu3(nTraceNumPoints);
1224  Array<OneD, NekDouble> Bwdmu3(nTraceNumPoints);
1225 
1226  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps1,
1227  Bwdeps1);
1228  m_fields[0]->GetFwdBwdTracePhys(epsvec[1], Fwdeps2,
1229  Bwdeps2);
1230  m_fields[0]->GetFwdBwdTracePhys(muvec[2], Fwdmu3, Bwdmu3);
1231 
1232  CopyBoundaryTrace(Fwdeps1, Bwdeps1, eFwdEQBwd, 0);
1233  CopyBoundaryTrace(Fwdeps2, Bwdeps2, eFwdEQBwd, 1);
1234  CopyBoundaryTrace(Fwdmu3, Bwdmu3, eFwdEQBwd, 2);
1235 
1236  for (int i = 0; i < nTraceNumPoints; ++i)
1237  {
1238  m_ZimFwd[0][i] = sqrt(Fwdmu3[i] / Fwdeps2[i]);
1239  m_ZimBwd[0][i] = sqrt(Bwdmu3[i] / Bwdeps2[i]);
1240 
1241  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1242  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1243 
1244  m_ZimFwd[1][i] = sqrt(Fwdmu3[i] / Fwdeps1[i]);
1245  m_ZimBwd[1][i] = sqrt(Bwdmu3[i] / Bwdeps1[i]);
1246 
1247  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1248  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1249  }
1250  }
1251  break; // eTransELectric
1252 
1253  default:
1254  break;
1255  } // PolType
1256 
1257  std::cout << "*** ZimFwd0 = [ "
1258  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[0], 1) << " , "
1259  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[0], 1)
1260  << " ], ZimBwd0 = [ "
1261  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[0], 1) << " , "
1262  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[0], 1) << " ] "
1263  << std::endl;
1264  std::cout << "*** ZimFwd1 = [ "
1265  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[1], 1) << " , "
1266  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[1], 1)
1267  << " ], ZimBwd1 = [ "
1268  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[1], 1) << " , "
1269  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[1], 1) << " ] "
1270  << std::endl;
1271  }
1272  break; // eMaxwell2D
1273 
1274  default:
1275  break;
1276  } // TestMaxwellType
1277 }
1278 
1279 // m_dedxi_cdot_e[m][j][k][] = de^m / d \xi^j \cdot e^n
1280 void MMFSystem::Computedemdxicdote()
1281 {
1282  int MFdim = 3;
1283  int nq = GetTotPoints();
1284 
1285  // Initialization
1286  m_dedxi_cdot_e =
1287  Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(MFdim);
1288  for (int indm = 0; indm < MFdim; ++indm)
1289  {
1290  m_dedxi_cdot_e[indm] =
1291  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
1292  for (int indj = 0; indj < MFdim; ++indj)
1293  {
1294  m_dedxi_cdot_e[indm][indj] =
1295  Array<OneD, Array<OneD, NekDouble>>(MFdim);
1296  for (int indn = 0; indn < MFdim; ++indn)
1297  {
1298  m_dedxi_cdot_e[indm][indj][indn] =
1299  Array<OneD, NekDouble>(nq, 0.0);
1300  }
1301  }
1302  }
1303 
1304  Array<OneD, NekDouble> tmp(nq);
1305  Array<OneD, NekDouble> tmpderiv(nq);
1306  Array<OneD, NekDouble> dedt;
1307  for (int indm = 0; indm < MFdim; ++indm)
1308  {
1309  for (int indj = 0; indj < MFdim; ++indj)
1310  {
1311  for (int indn = 0; indn < MFdim; ++indn)
1312  {
1313  dedt = Array<OneD, NekDouble>(nq, 0.0);
1314  for (int k = 0; k < m_spacedim; ++k)
1315  {
1316  // Compute d e^m / d \xi_j cdot e^n
1317  Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0],
1318  1);
1319  m_fields[0]->PhysDirectionalDeriv(m_movingframes[indj], tmp,
1320  tmpderiv);
1321 
1322  Vmath::Vvtvp(nq, &tmpderiv[0], 1,
1323  &m_movingframes[indn][k * nq], 1, &dedt[0], 1,
1324  &dedt[0], 1);
1325  }
1326 
1327  Vmath::Vcopy(nq, &dedt[0], 1,
1328  &m_dedxi_cdot_e[indm][indj][indn][0], 1);
1329  }
1330  }
1331  }
1332 
1333  int indx = 0;
1334  std::cout << "*** m_dedxi_cdot_e[0]/dxi1 = ( "
1335  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1336  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1337  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1338  << std::endl;
1339  std::cout << "*** m_dedxi_cdot_e[0]/dxi2 = ( "
1340  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1341  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1342  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1343  << std::endl;
1344 
1345  indx = 1;
1346  std::cout << "*** m_dedxi_cdot_e[1]/dxi1 = ( "
1347  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1348  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1349  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1350  << std::endl;
1351  std::cout << "*** m_dedxi_cdot_e[1]/dxi2 = ( "
1352  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1353  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1354  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1355  << std::endl;
1356 
1357  indx = 2;
1358  std::cout << "*** m_dedxi_cdot_e[2]/dxi1 = ( "
1359  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1360  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1361  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1362  << std::endl;
1363  std::cout << "*** m_dedxi_cdot_e[2]/dxi2 = ( "
1364  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1365  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1366  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1367  << std::endl;
1368 }
1369 
1370 void MMFSystem::AdddedtMaxwell(
1371  const Array<OneD, const Array<OneD, NekDouble>> &physarray,
1372  Array<OneD, Array<OneD, NekDouble>> &outarray)
1373 {
1374  int nq = GetTotPoints();
1375 
1376  // m_dedxi_cdot_e[m][j][n][] = de^m / d \xi^j \cdot e^n
1377  Array<OneD, NekDouble> dedtej(nq);
1378  Array<OneD, NekDouble> de1dtcdotej(nq);
1379  Array<OneD, NekDouble> de2dtcdotej(nq);
1380 
1381  Array<OneD, NekDouble> normH(nq);
1382  Array<OneD, NekDouble> NormedH1(nq, 0.0);
1383  Array<OneD, NekDouble> NormednegH2(nq, 0.0);
1384 
1385  Vmath::Vmul(nq, &physarray[0][0], 1, &physarray[0][0], 1, &normH[0], 1);
1386  Vmath::Vvtvp(nq, &physarray[1][0], 1, &physarray[1][0], 1, &normH[0], 1,
1387  &normH[0], 1);
1388  Vmath::Vsqrt(nq, normH, 1, normH, 1);
1389 
1390  NekDouble Tol = 0.001;
1391  for (int i = 0; i < nq; ++i)
1392  {
1393  if (normH[i] > Tol)
1394  {
1395  NormedH1[i] = physarray[0][i] / normH[i];
1396  NormednegH2[i] = -1.0 * physarray[1][i] / normH[i];
1397  }
1398  }
1399 
1400  for (int j = 0; j < m_shapedim; ++j)
1401  {
1402  // Compute de1 / dt \cdot ej = (-H2 de^1/d\xi1 \cdot e^j + H1 de^1/d\xi2
1403  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1404  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[0][0][j][0], 1,
1405  &de1dtcdotej[0], 1);
1406  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[0][1][j][0], 1,
1407  &de1dtcdotej[0], 1, &de1dtcdotej[0], 1);
1408 
1409  // Compute de2 / dt \cdot ej = (-H2 de2/d\xi1 \cdot e^j + H1 de2/d\xi2
1410  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1411  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[1][0][j][0], 1,
1412  &de2dtcdotej[0], 1);
1413  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[1][1][j][0], 1,
1414  &de2dtcdotej[0], 1, &de2dtcdotej[0], 1);
1415 
1416  // Add dedt component: (H1 (de1/dt) + H2 (de2/dt) ) \cdot ej
1417  Vmath::Vmul(nq, &physarray[0][0], 1, &de1dtcdotej[0], 1, &dedtej[0], 1);
1418  Vmath::Vvtvp(nq, &physarray[1][0], 1, &de2dtcdotej[0], 1, &dedtej[0], 1,
1419  &dedtej[0], 1);
1420 
1421  Vmath::Neg(nq, dedtej, 1);
1422 
1423  switch (m_PolType)
1424  {
1425  case SolverUtils::eTransMagnetic:
1426  {
1427  if (j == 0)
1428  {
1429  Vmath::Vmul(nq, m_muvec[0], 1, dedtej, 1, dedtej, 1);
1430  }
1431 
1432  else if (j == 1)
1433  {
1434  Vmath::Vmul(nq, m_muvec[1], 1, dedtej, 1, dedtej, 1);
1435  }
1436  }
1437  break;
1438 
1439  case SolverUtils::eTransElectric:
1440  {
1441  if (j == 0)
1442  {
1443  Vmath::Vmul(nq, m_epsvec[0], 1, dedtej, 1, dedtej, 1);
1444  }
1445 
1446  else if (j == 1)
1447  {
1448  Vmath::Vmul(nq, m_epsvec[1], 1, dedtej, 1, dedtej, 1);
1449  }
1450  }
1451  break;
1452 
1453  default:
1454  break;
1455  }
1456 
1457  Vmath::Vadd(nq, &dedtej[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
1458  }
1459 }
1460 
1461 void MMFSystem::UpwindMaxwellFlux1D(
1462  Array<OneD, Array<OneD, NekDouble>> &physfield,
1463  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1464  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1465 {
1466  int i;
1467  int nTraceNumPoints = GetTraceTotPoints();
1468  int nvar = 2;
1469 
1470  // get temporary arrays
1471  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1472  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1473 
1474  for (i = 0; i < nvar; ++i)
1475  {
1476  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1477  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1478  }
1479 
1480  // get the physical values at the trace from the dependent variables
1481  for (i = 0; i < nvar; ++i)
1482  {
1483  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1484  }
1485 
1486  // E = 0 at the boundaries
1487  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1488 
1489  // d H / d n = 0 at the boundaries
1490  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1491 
1492  Array<OneD, NekDouble> dE(nTraceNumPoints);
1493  Array<OneD, NekDouble> dH(nTraceNumPoints);
1494 
1495  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1496  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1497 
1498  NekDouble nx, AverZ, AverY, AverZH, AverYE;
1499  for (i = 0; i < nTraceNumPoints; ++i)
1500  {
1501  nx = m_traceNormals[0][i];
1502  AverZ = m_ZimFwd[0][i] + m_ZimBwd[0][i];
1503  AverY = m_YimFwd[0][i] + m_YimBwd[0][i];
1504  AverZH = m_ZimFwd[0][i] * Fwd[1][i] + m_ZimBwd[0][i] * Bwd[1][i];
1505  AverYE = m_YimFwd[0][i] * Fwd[0][i] + m_YimBwd[0][i] * Bwd[0][i];
1506 
1507  numfluxFwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1508  numfluxFwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1509 
1510  numfluxBwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1511  numfluxBwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1512  }
1513 }
1514 
1515 void MMFSystem::LaxFriedrichMaxwellFlux1D(
1516  Array<OneD, Array<OneD, NekDouble>> &physfield,
1517  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1518  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1519 {
1520  int i;
1521  int nTraceNumPoints = GetTraceTotPoints();
1522  int nvar = 2;
1523 
1524  // get temporary arrays
1525  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1526  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1527 
1528  for (i = 0; i < nvar; ++i)
1529  {
1530  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1531  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1532  }
1533 
1534  // get the physical values at the trace from the dependent variables
1535  for (i = 0; i < nvar; ++i)
1536  {
1537  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1538  }
1539 
1540  // E = 0 at the boundaries
1541  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1542 
1543  // d H / d n = 0 at the boundaries
1544  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1545 
1546  Array<OneD, NekDouble> dE(nTraceNumPoints);
1547  Array<OneD, NekDouble> dH(nTraceNumPoints);
1548 
1549  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1550  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1551 
1552  NekDouble nx;
1553  for (i = 0; i < nTraceNumPoints; ++i)
1554  {
1555  nx = m_traceNormals[0][i];
1556  numfluxFwd[0][i] =
1557  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1558  numfluxFwd[1][i] =
1559  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1560 
1561  numfluxBwd[0][i] =
1562  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1563  numfluxBwd[1][i] =
1564  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1565  }
1566 }
1567 
1568 void MMFSystem::AverageMaxwellFlux1D(
1569  Array<OneD, Array<OneD, NekDouble>> &physfield,
1570  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1571  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1572 {
1573  int i;
1574  int nTraceNumPoints = GetTraceTotPoints();
1575  int nvar = 2;
1576 
1577  // get temporary arrays
1578  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1579  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1580 
1581  for (i = 0; i < nvar; ++i)
1582  {
1583  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1584  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1585  }
1586 
1587  // get the physical values at the trace from the dependent variables
1588  for (i = 0; i < nvar; ++i)
1589  {
1590  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1591  }
1592 
1593  // E = 0 at the boundaries
1594  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1595 
1596  // d H / d n = 0 at the boundaries
1597  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1598 
1599  for (i = 0; i < nTraceNumPoints; ++i)
1600  {
1601  numfluxFwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1602  numfluxFwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1603 
1604  numfluxBwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1605  numfluxBwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1606  }
1607 }
1608 
1609 void MMFSystem::GetMaxwellFluxVector(
1610  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1611  Array<OneD, Array<OneD, NekDouble>> &flux)
1612 {
1613  switch (m_TestMaxwellType)
1614  {
1615  case eMaxwell1D:
1616  case eScatField1D:
1617  {
1618  GetMaxwellFlux1D(var, physfield, flux);
1619  }
1620  break;
1621 
1622  case eTestMaxwell2DPEC:
1623  case eTestMaxwell2DPECAVGFLUX:
1624  case eTestMaxwell2DPMC:
1625  case eScatField2D:
1626  case eTotField2D:
1627  case eMaxwellSphere:
1628  case eELF2DSurface:
1629  {
1630  GetMaxwellFlux2D(var, physfield, flux);
1631  }
1632  break;
1633 
1634  default:
1635  break;
1636  }
1637 }
1638 
1639 void MMFSystem::GetMaxwellFlux1D(
1640  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1641  Array<OneD, Array<OneD, NekDouble>> &flux)
1642 {
1643  int nq = m_fields[0]->GetTotPoints();
1644 
1645  switch (var)
1646  {
1647  case 0:
1648  {
1649  // H in flux 0
1650  Vmath::Vcopy(nq, physfield[1], 1, flux[0], 1);
1651 
1652  // E in flux 1
1653  Vmath::Zero(nq, flux[1], 1);
1654  }
1655  break;
1656 
1657  case 1:
1658  {
1659  // E in flux 0
1660  Vmath::Vcopy(nq, physfield[0], 1, flux[0], 1);
1661 
1662  // H in flux 1
1663  Vmath::Zero(nq, flux[1], 1);
1664  }
1665  break;
1666  //----------------------------------------------------
1667 
1668  default:
1669  break;
1670  }
1671 }
1672 
1673 void MMFSystem::GetMaxwellFlux2D(
1674  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1675  Array<OneD, Array<OneD, NekDouble>> &flux)
1676 {
1677  int nq = m_fields[0]->GetTotPoints();
1678 
1679  NekDouble sign = 1.0;
1680  switch (m_PolType)
1681  {
1682  // TransMagnetic
1683  case 0:
1684  {
1685  sign = -1.0;
1686  }
1687  break;
1688 
1689  // TransElectric
1690  case 1:
1691  {
1692  sign = 1.0;
1693  }
1694  break;
1695 
1696  default:
1697  break;
1698  }
1699 
1700  switch (var)
1701  {
1702  case 0:
1703  {
1704  // -Ez in flux 1
1705  Vmath::Smul(nq, sign, physfield[2], 1, flux[0], 1);
1706  Vmath::Zero(nq, flux[1], 1);
1707  }
1708  break;
1709 
1710  case 1:
1711  {
1712  // Ez in flux 0
1713  Vmath::Zero(nq, flux[0], 1);
1714  Vmath::Smul(nq, -sign, physfield[2], 1, flux[1], 1);
1715  }
1716  break;
1717 
1718  case 2:
1719  {
1720  Vmath::Smul(nq, sign, physfield[0], 1, flux[0], 1);
1721  Vmath::Smul(nq, -sign, physfield[1], 1, flux[1], 1);
1722  }
1723  break;
1724 
1725  default:
1726  ASSERTL0(false, "GetFluxVector2D: illegal vector index");
1727  }
1728 }
1729 
1730 void MMFSystem::NumericalMaxwellFlux(
1731  Array<OneD, Array<OneD, NekDouble>> &physfield,
1732  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1733  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1734 {
1735 
1736  switch (m_TestMaxwellType)
1737  {
1738  case eMaxwell1D:
1739  case eScatField1D:
1740  {
1741  switch (m_upwindType)
1742  {
1743  case eAverage:
1744  {
1745  AverageMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1746  }
1747  break;
1748 
1749  case eLaxFriedrich:
1750  {
1751  LaxFriedrichMaxwellFlux1D(physfield, numfluxFwd,
1752  numfluxBwd);
1753  }
1754  break;
1755 
1756  case eUpwind:
1757  {
1758  UpwindMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1759  }
1760  break;
1761 
1762  default:
1763  {
1764  ASSERTL0(false,
1765  "populate switch statement for upwind flux");
1766  }
1767  break; // upwindType
1768  }
1769  }
1770  break; // eMaxwell1D
1771 
1772  case eTestMaxwell2DPEC:
1773  case eTestMaxwell2DPECAVGFLUX:
1774  case eTestMaxwell2DPMC:
1775  case eScatField2D:
1776  case eTotField2D:
1777  case eMaxwellSphere:
1778  case eELF2DSurface:
1779  {
1780  switch (m_PolType)
1781  {
1782  case eTransMagnetic:
1783  {
1784  NumericalMaxwellFluxTM(physfield, numfluxFwd, numfluxBwd,
1785  time);
1786  }
1787  break;
1788 
1789  case eTransElectric:
1790  {
1791  NumericalMaxwellFluxTE(physfield, numfluxFwd, numfluxBwd,
1792  time);
1793  }
1794  break;
1795 
1796  default:
1797  break;
1798  }
1799  }
1800  break; // eMaxwell2D
1801 
1802  default:
1803  break;
1804  } // m_TestMaxwellType
1805 }
1806 
1807 void MMFSystem::NumericalMaxwellFluxTM(
1808  Array<OneD, Array<OneD, NekDouble>> &physfield,
1809  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1810  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1811 {
1812  int nq = m_fields[0]->GetNpoints();
1813  int nTraceNumPoints = GetTraceTotPoints();
1814  int nvar = physfield.size();
1815 
1816  // get temporary arrays
1817  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1818  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1819  for (int i = 0; i < nvar; ++i)
1820  {
1821  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1822  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1823 
1824  // get the physical values at the trace
1825  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1826  }
1827 
1828  // E^|| = 0 at the PEC boundaries vs. H^|| = 0 at PMC boundaries
1829  Array<OneD, NekDouble> IncField(nq, 0.0);
1830  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1831  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1832  for (int i = 0; i < m_spacedim; ++i)
1833  {
1834  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1835 
1836  IncField = GetIncidentField(i, time);
1837  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1838 
1839  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1840  &IncFieldFwd[i][0], 1);
1841  }
1842 
1843  // Total Field Formulation
1844  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC_Forces");
1845  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC_Forces");
1846  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1847  "PEC_Forces");
1848 
1849  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1850  "PMC_Forces");
1851  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1852  "PMC_Forces");
1853  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC_Forces");
1854 
1855  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PEC");
1856  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PEC");
1857  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PEC");
1858 
1859  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PMC");
1860  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PMC");
1861  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PMC");
1862 
1863  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1864  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1865  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1866  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1867  Array<OneD, NekDouble> e3Fwd_cdot_dEe3(nTraceNumPoints, 0.0);
1868  Array<OneD, NekDouble> e3Bwd_cdot_dEe3(nTraceNumPoints, 0.0);
1869 
1870  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (YimFwd *
1871  // EFwd^3 + YimBwd * EBwd^3 )
1872  ComputeNtimesFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], numfluxFwd[0],
1873  numfluxBwd[0]);
1874  ComputeNtimesFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1], numfluxFwd[1],
1875  numfluxBwd[1]);
1876 
1877  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( ZimFwd HFwd + ZimBwd
1878  // HBwd ) ) / 2 {{Z_i}}
1879  ComputeNtimesF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1880  m_ZimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1881 
1882  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1883  // [H] / 2 {{YimFwd}}
1884  ComputeNtimestimesdFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0],
1885  e1Fwd_cdot_ncrossdH, e1Bwd_cdot_ncrossdH);
1886  ComputeNtimestimesdFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1],
1887  e2Fwd_cdot_ncrossdH, e2Bwd_cdot_ncrossdH);
1888 
1889  // Compute \alpha [E3] * ( 1/2{{Zim1}} + 1/2{{Zim2}} )
1890  ComputeNtimestimesdF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1891  m_ZimBwd[1], e3Fwd_cdot_dEe3, e3Bwd_cdot_dEe3);
1892 
1893  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdH,
1894  1, numfluxFwd[0], 1);
1895  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdH,
1896  1, numfluxFwd[1], 1);
1897  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[2], 1, e3Fwd_cdot_dEe3, 1,
1898  numfluxFwd[2], 1);
1899 
1900  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdH,
1901  1, numfluxBwd[0], 1);
1902  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdH,
1903  1, numfluxBwd[1], 1);
1904  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[2], 1, e3Bwd_cdot_dEe3, 1,
1905  numfluxBwd[2], 1);
1906 }
1907 
1908 void MMFSystem::NumericalMaxwellFluxTE(
1909  Array<OneD, Array<OneD, NekDouble>> &physfield,
1910  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1911  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1912 
1913 {
1914  int nq = m_fields[0]->GetNpoints();
1915  int nTraceNumPoints = GetTraceTotPoints();
1916  int nvar = physfield.size();
1917 
1918  // Get temporary arrays
1919  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1920  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1921  for (int i = 0; i < nvar; ++i)
1922  {
1923  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1924  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1925 
1926  // get the physical values at the trace
1927  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1928  }
1929 
1930  // E = 0 at the PEC boundaries:
1931  Array<OneD, NekDouble> IncField(nq, 0.0);
1932  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1933  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1934  for (int i = 0; i < m_spacedim; ++i)
1935  {
1936  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1937 
1938  IncField = GetIncidentField(i, time);
1939  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1940 
1941  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1942  &IncFieldFwd[i][0], 1);
1943  }
1944 
1945  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1946  "PEC_Forces");
1947  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1948  "PEC_Forces");
1949  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC_Forces");
1950 
1951  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC_Forces");
1952  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC_Forces");
1953  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1954  "PMC_Forces");
1955 
1956  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0, "PEC");
1957  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1, "PEC");
1958  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2, "PEC");
1959 
1960  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0, "PMC");
1961  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1, "PMC");
1962  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2, "PMC");
1963 
1964  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdE(nTraceNumPoints);
1965  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdE(nTraceNumPoints);
1966  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdE(nTraceNumPoints);
1967  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdE(nTraceNumPoints);
1968  Array<OneD, NekDouble> e3Fwd_cdot_dHe3(nTraceNumPoints);
1969  Array<OneD, NekDouble> e3Bwd_cdot_dHe3(nTraceNumPoints);
1970 
1971  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (ZimFwd *
1972  // HFwd^3 + ZimBwd * HBwd^3 )
1973  // Compute numfluxBwd[dir] = (eBwd^[dir] \cdot n \times e^3) * (ZimFwd *
1974  // HFwd^3 + ZimBwd * HBwd^3 )
1975  ComputeNtimesFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], numfluxFwd[0],
1976  numfluxBwd[0]);
1977  ComputeNtimesFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1], numfluxFwd[1],
1978  numfluxBwd[1]);
1979 
1980  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd
1981  // EBwd ) ) / 2 {{Y_i}}
1982  ComputeNtimesF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1983  m_YimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1984 
1985  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1986  // [E] / 2 {{ZimFwd}}
1987  ComputeNtimestimesdFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0],
1988  e1Fwd_cdot_ncrossdE, e1Bwd_cdot_ncrossdE);
1989  ComputeNtimestimesdFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1],
1990  e2Fwd_cdot_ncrossdE, e2Bwd_cdot_ncrossdE);
1991 
1992  // Compute - \alpha [H3] * ( 1/2{{Yim1}} + 1/2{{Yim2}} )
1993  ComputeNtimestimesdF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
1994  m_YimBwd[1], e3Fwd_cdot_dHe3, e3Bwd_cdot_dHe3);
1995 
1996  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdE, 1,
1997  numfluxFwd[0], 1);
1998  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdE, 1,
1999  numfluxFwd[1], 1);
2000  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[2], 1, e3Fwd_cdot_dHe3, 1,
2001  numfluxFwd[2], 1);
2002 
2003  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdE, 1,
2004  numfluxBwd[0], 1);
2005  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdE, 1,
2006  numfluxBwd[1], 1);
2007  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[2], 1, e3Bwd_cdot_dHe3, 1,
2008  numfluxBwd[2], 1);
2009 }
2010 
2011 Array<OneD, NekDouble> MMFSystem::GetIncidentField(const int var,
2012  const NekDouble time)
2013 {
2014  int nq = m_fields[0]->GetNpoints();
2015 
2016  Array<OneD, NekDouble> x0(nq);
2017  Array<OneD, NekDouble> x1(nq);
2018  Array<OneD, NekDouble> x2(nq);
2019 
2020  m_fields[0]->GetCoords(x0, x1, x2);
2021 
2022  // GetSmoothFactor such that wave propages from the left to the object.
2023  // a = 0.1, ta = 1, f = 1.0./(1.0 + exp( -0.5.*(time-ta)/a ));
2024  Array<OneD, NekDouble> SmoothFactor(nq, 1.0);
2025  switch (m_SmoothFactor)
2026  {
2027  case 0:
2028  {
2029  for (int i = 0; i < nq; i++)
2030  {
2031  SmoothFactor[i] = 1.0 / (1.0 + exp(-1.0 * (time - 1.0) / 0.1));
2032  }
2033  }
2034  break;
2035 
2036  case 1:
2037  {
2038  NekDouble xmin = Vmath::Vmin(nq, x0, 1);
2039  NekDouble xp;
2040  for (int i = 0; i < nq; i++)
2041  {
2042  xp = x0[i] - xmin - time - m_SFinit;
2043  if (xp > 0.0)
2044  {
2045  SmoothFactor[i] =
2046  2.0 / (1.0 + exp(0.5 * (sqrt(xp * xp) - 0.1)));
2047  }
2048 
2049  else
2050  {
2051  SmoothFactor[i] = 1.0;
2052  }
2053  }
2054  }
2055  break;
2056 
2057  default:
2058  break;
2059  }
2060 
2061  // Generate a factor for smoothly increasing wave
2062  Array<OneD, NekDouble> F1(nq);
2063  Array<OneD, NekDouble> F2(nq);
2064  Array<OneD, NekDouble> F3(nq);
2065 
2066  Array<OneD, NekDouble> dF1dt(nq);
2067  Array<OneD, NekDouble> dF2dt(nq);
2068  Array<OneD, NekDouble> dF3dt(nq);
2069 
2070  Array<OneD, NekDouble> F1int(nq);
2071  Array<OneD, NekDouble> F2int(nq);
2072  Array<OneD, NekDouble> F3int(nq);
2073 
2074  Array<OneD, NekDouble> outarray(nq);
2075 
2076  switch (m_IncType)
2077  {
2078  case ePlaneWave:
2079  {
2080  NekDouble cs, sn;
2081  NekDouble e1y, e2y;
2082  switch (m_PolType)
2083  {
2084  case eTransMagnetic:
2085  {
2086  // H = 0 \hat{x} - e^{ikx} \hat{y}
2087  // E = e^{ikx} \hat{z}
2088  // outarray1 = Hr1inc
2089  // outarray2 = Hr2inc
2090  // outarray3 = Ezrinc
2091  for (int i = 0; i < nq; i++)
2092  {
2093  e1y = m_movingframes[0][nq + i];
2094  e2y = m_movingframes[1][nq + i];
2095 
2096  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2097  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2098 
2099  F1[i] = -1.0 * cs * e1y;
2100  F2[i] = -1.0 * cs * e2y;
2101  F3[i] = cs;
2102 
2103  dF1dt[i] = -1.0 * m_Incfreq * sn * e1y;
2104  dF2dt[i] = -1.0 * m_Incfreq * sn * e2y;
2105  dF3dt[i] = 1.0 * m_Incfreq * sn;
2106 
2107  F1int[i] = (1.0 / m_Incfreq) * sn * e1y;
2108  F2int[i] = (1.0 / m_Incfreq) * sn * e2y;
2109  F3int[i] = (-1.0 / m_Incfreq) * sn;
2110  }
2111  }
2112  break;
2113 
2114  case eTransElectric:
2115  {
2116  // E = 0 \hat{x} + e^{ikx} \hat{y}
2117  // H = e^{ikx} \hat{z}
2118  // outarray1 = Er1inc
2119  // outarray2 = Er2inc
2120  // outarray3 = Hzrinc
2121  // outarray1 = Ei1inc
2122  // outarray2 = Ei2inc
2123  // outarray3 = Hziinc
2124  for (int i = 0; i < nq; i++)
2125  {
2126  e1y = m_movingframes[0][nq + i];
2127  e2y = m_movingframes[1][nq + i];
2128 
2129  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2130  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2131 
2132  F1[i] = cs * e1y;
2133  F2[i] = cs * e2y;
2134  F3[i] = cs;
2135 
2136  dF1dt[i] = m_Incfreq * sn * e1y;
2137  dF2dt[i] = m_Incfreq * sn * e2y;
2138  dF3dt[i] = m_Incfreq * sn;
2139 
2140  F1int[i] = (-1.0 / m_Incfreq) * sn * e1y;
2141  F2int[i] = (-1.0 / m_Incfreq) * sn * e2y;
2142  F3int[i] = (-1.0 / m_Incfreq) * sn;
2143  }
2144  }
2145  break;
2146 
2147  default:
2148  break;
2149  }
2150  }
2151  break;
2152 
2153  case ePlaneWaveImag:
2154  {
2155  NekDouble cs, sn;
2156  NekDouble e1y, e2y;
2157  switch (m_PolType)
2158  {
2159  case eTransMagnetic:
2160  {
2161  // H = 0 \hat{x} - e^{ikx} \hat{y}
2162  // E = e^{ikx} \hat{z}
2163  // outarray1 = Hr1inc
2164  // outarray2 = Hr2inc
2165  // outarray3 = Ezrinc
2166  for (int i = 0; i < nq; i++)
2167  {
2168  e1y = m_movingframes[0][nq + i];
2169  e2y = m_movingframes[1][nq + i];
2170 
2171  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2172  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2173 
2174  F1[i] = -1.0 * sn * e1y;
2175  F2[i] = -1.0 * sn * e2y;
2176  F3[i] = sn;
2177 
2178  dF1dt[i] = m_Incfreq * cs * e1y;
2179  dF2dt[i] = m_Incfreq * cs * e2y;
2180  dF3dt[i] = -1.0 * m_Incfreq * cs;
2181 
2182  F1int[i] = (-1.0 / m_Incfreq) * cs * e1y;
2183  F2int[i] = (-1.0 / m_Incfreq) * cs * e2y;
2184  F3int[i] = (1.0 / m_Incfreq) * cs;
2185  }
2186  }
2187  break;
2188 
2189  case eTransElectric:
2190  {
2191  // E = 0 \hat{x} + e^{ikx} \hat{y}
2192  // H = e^{ikx} \hat{z}
2193  // outarray1 = Er1inc
2194  // outarray2 = Er2inc
2195  // outarray3 = Hzrinc
2196  // outarray1 = Ei1inc
2197  // outarray2 = Ei2inc
2198  // outarray3 = Hziinc
2199  for (int i = 0; i < nq; i++)
2200  {
2201  e1y = m_movingframes[0][nq + i];
2202  e2y = m_movingframes[1][nq + i];
2203 
2204  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2205  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2206 
2207  F1[i] = sn * e1y;
2208  F2[i] = sn * e2y;
2209  F3[i] = sn;
2210 
2211  dF1dt[i] = -1.0 * m_Incfreq * cs * e1y;
2212  dF2dt[i] = -1.0 * m_Incfreq * cs * e2y;
2213  dF3dt[i] = -1.0 * m_Incfreq * cs;
2214 
2215  F1int[i] = (1.0 / m_Incfreq) * cs * e1y;
2216  F2int[i] = (1.0 / m_Incfreq) * cs * e2y;
2217  F3int[i] = (1.0 / m_Incfreq) * cs;
2218  }
2219  }
2220  break;
2221 
2222  default:
2223  break;
2224  }
2225  }
2226  break;
2227 
2228  default:
2229  break;
2230  }
2231 
2232  switch (var)
2233  {
2234  case 0:
2235  {
2236  outarray = F1;
2237  }
2238  break;
2239 
2240  case 1:
2241  {
2242  outarray = F2;
2243  }
2244  break;
2245 
2246  case 2:
2247  {
2248  outarray = F3;
2249  }
2250  break;
2251 
2252  case 10:
2253  {
2254  outarray = dF1dt;
2255  }
2256  break;
2257 
2258  case 11:
2259  {
2260  outarray = dF2dt;
2261  }
2262  break;
2263 
2264  case 12:
2265  {
2266  outarray = dF3dt;
2267  }
2268  break;
2269 
2270  case 20:
2271  {
2272  outarray = F1int;
2273  }
2274  break;
2275 
2276  case 21:
2277  {
2278  outarray = F2int;
2279  }
2280  break;
2281 
2282  case 22:
2283  {
2284  outarray = F3int;
2285  }
2286  break;
2287 
2288  default:
2289  {
2290  Vmath::Zero(nq, outarray, 1);
2291  }
2292  break;
2293  }
2294 
2295  return outarray;
2296 }
2297 
2298 NekDouble MMFSystem::AvgInt(const Array<OneD, const NekDouble> &inarray)
2299 {
2300  int nq = m_fields[0]->GetNpoints();
2301  Array<OneD, NekDouble> Ones(nq, 1.0);
2302 
2303  if (inarray.size() != nq)
2304  {
2305  ASSERTL0(false, "AvgInt Error: Vector size is not correct");
2306  }
2307 
2308  NekDouble jac = m_fields[0]->Integral(Ones);
2309 
2310  return (m_fields[0]->Integral(inarray)) / jac;
2311 }
2312 
2313 NekDouble MMFSystem::AvgAbsInt(const Array<OneD, const NekDouble> &inarray)
2314 {
2315  int nq = m_fields[0]->GetNpoints();
2316  Array<OneD, NekDouble> Ones(nq, 1.0);
2317  Array<OneD, NekDouble> tmp(nq);
2318 
2319  if (inarray.size() != nq)
2320  {
2321  ASSERTL0(false, "AvgAbsInt Error: Vector size is not correct");
2322  }
2323 
2324  NekDouble jac = m_fields[0]->Integral(Ones);
2325 
2326  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2327  return (m_fields[0]->Integral(tmp)) / jac;
2328 }
2329 
2330 NekDouble MMFSystem::AbsIntegral(const Array<OneD, const NekDouble> &inarray)
2331 {
2332  int nq = m_fields[0]->GetNpoints();
2333  Array<OneD, NekDouble> tmp(nq);
2334 
2335  if (inarray.size() != nq)
2336  {
2337  ASSERTL0(false, "AbsIntegral Error: Vector size is not correct");
2338  }
2339 
2340  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2341  return m_fields[0]->Integral(tmp);
2342 }
2343 
2344 NekDouble MMFSystem::RootMeanSquare(const Array<OneD, const NekDouble> &inarray)
2345 {
2346  int Ntot = inarray.size();
2347 
2348  NekDouble reval = 0.0;
2349  for (int i = 0; i < Ntot; ++i)
2350  {
2351  reval = reval + inarray[i] * inarray[i];
2352  }
2353  reval = sqrt(reval / Ntot);
2354 
2355  return reval;
2356 }
2357 
2358 NekDouble MMFSystem::VectorAvgMagnitude(
2359  const Array<OneD, const Array<OneD, NekDouble>> &inarray)
2360 {
2361  int nq = inarray[0].size();
2362 
2363  Array<OneD, NekDouble> tmp(nq, 0.0);
2364  for (int k = 0; k < m_spacedim; k++)
2365  {
2366  Vmath::Vvtvp(nq, &inarray[k][0], 1, &inarray[k][0], 1, &tmp[0], 1,
2367  &tmp[0], 1);
2368  }
2369  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
2370 
2371  return RootMeanSquare(tmp);
2372 }
2373 
2374 void MMFSystem::BubbleSort(Array<OneD, NekDouble> &refarray,
2375  Array<OneD, NekDouble> &sortarray)
2376 {
2377  int nq = refarray.size();
2378 
2379  bool swapped = true;
2380  int j = 0;
2381  NekDouble tmp;
2382 
2383  while (swapped)
2384  {
2385  swapped = false;
2386  j++;
2387  for (int i = 0; i < nq - j; i++)
2388  {
2389  if (refarray[i] > refarray[i + 1])
2390  {
2391  tmp = refarray[i];
2392  refarray[i] = refarray[i + 1];
2393  refarray[i + 1] = tmp;
2394 
2395  tmp = sortarray[i];
2396  sortarray[i] = sortarray[i + 1];
2397  sortarray[i + 1] = tmp;
2398 
2399  swapped = true;
2400  }
2401  }
2402  }
2403 }
2404 
2405 void MMFSystem::GramSchumitz(
2406  const Array<OneD, const Array<OneD, NekDouble>> &v1,
2407  const Array<OneD, const Array<OneD, NekDouble>> &v2,
2408  Array<OneD, Array<OneD, NekDouble>> &outarray, bool KeepTheMagnitude)
2409 {
2410 
2411  int nq = v1[0].size();
2412  Array<OneD, NekDouble> tmp(nq, 0.0);
2413  Array<OneD, NekDouble> mag(nq, 0.0);
2414 
2415  for (int i = 0; i < m_spacedim; ++i)
2416  {
2417  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2418  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &tmp[0], 1, &tmp[0], 1);
2419  Vmath::Vvtvp(nq, &v1[i][0], 1, &v1[i][0], 1, &mag[0], 1, &mag[0], 1);
2420  }
2421  Vmath::Vdiv(nq, &tmp[0], 1, &mag[0], 1, &tmp[0], 1);
2422  Vmath::Neg(nq, &tmp[0], 1);
2423 
2424  // outarray = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
2425 
2426  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2427  for (int i = 0; i < m_spacedim; ++i)
2428  {
2429  outarray[i] = Array<OneD, NekDouble>(nq, 0.0);
2430  Vmath::Vvtvp(nq, &tmp[0], 1, &v1[i][0], 1, &v2[i][0], 1,
2431  &outarray[i][0], 1);
2432  }
2433 
2434  if (KeepTheMagnitude)
2435  {
2436  Array<OneD, NekDouble> magorig(nq, 0.0);
2437  Array<OneD, NekDouble> magnew(nq, 0.0);
2438 
2439  for (int i = 0; i < m_spacedim; ++i)
2440  {
2441  Vmath::Vmul(nq, &v2[0][0], 1, &v2[0][0], 1, &magorig[0], 1);
2442  Vmath::Vvtvp(nq, &v2[1][0], 1, &v2[1][0], 1, &magorig[0], 1,
2443  &magorig[0], 1);
2444  Vmath::Vvtvp(nq, &v2[2][0], 1, &v2[2][0], 1, &magorig[0], 1,
2445  &magorig[0], 1);
2446 
2447  Vmath::Vmul(nq, &outarray[0][0], 1, &outarray[0][0], 1, &magnew[0],
2448  1);
2449  Vmath::Vvtvp(nq, &outarray[1][0], 1, &outarray[1][0], 1, &magnew[0],
2450  1, &magnew[0], 1);
2451  Vmath::Vvtvp(nq, &outarray[2][0], 1, &outarray[2][0], 1, &magnew[0],
2452  1, &magnew[0], 1);
2453  }
2454 
2455  for (int i = 0; i < m_spacedim; ++i)
2456  {
2457  for (int j = 0; j < nq; ++j)
2458  {
2459  if (fabs(magnew[j]) > 0.000000001)
2460  {
2461  outarray[i][j] =
2462  outarray[i][j] * sqrt(magorig[j] / magnew[j]);
2463  }
2464  }
2465  }
2466  }
2467 }
2468 
2469 void MMFSystem::v_GenerateSummary(SummaryList &s)
2470 {
2471  int nq = m_fields[0]->GetNpoints();
2472  UnsteadySystem::v_GenerateSummary(s);
2473 
2474  AddSummaryItem(s, "Total grids", nq);
2475  AddSummaryItem(s, "Shape Dimension", m_shapedim);
2476  AddSummaryItem(s, "Surface", SurfaceTypeMap[m_surfaceType]);
2477  if (m_surfaceType == eSphere)
2478  {
2479  NekDouble MeshError;
2480 
2481  Array<OneD, NekDouble> x(nq);
2482  Array<OneD, NekDouble> y(nq);
2483  Array<OneD, NekDouble> z(nq);
2484  Array<OneD, NekDouble> rad(nq, 0.0);
2485 
2486  m_fields[0]->GetCoords(x, y, z);
2487 
2488  Vmath::Vvtvp(nq, x, 1, x, 1, rad, 1, rad, 1);
2489  Vmath::Vvtvp(nq, y, 1, y, 1, rad, 1, rad, 1);
2490  Vmath::Vvtvp(nq, z, 1, z, 1, rad, 1, rad, 1);
2491  Vmath::Vsqrt(nq, rad, 1, rad, 1);
2492 
2493  Vmath::Sadd(nq, -1.0, rad, 1, rad, 1);
2494  Vmath::Vabs(nq, rad, 1, rad, 1);
2495 
2496  MeshError = m_fields[0]->Integral(rad);
2497  SolverUtils::AddSummaryItem(s, "Mesh Error", MeshError);
2498  }
2499 
2500  AddSummaryItem(s, "MMFdir", SpatialDomains::GeomMMFMap[m_MMFdir]);
2501  AddSummaryItem(s, "SmoothFactor", m_SmoothFactor);
2502  AddSummaryItem(s, "SFinit", m_SFinit);
2503 
2504  if (fabs(m_alpha - 1.0) > 0.001)
2505  {
2506  AddSummaryItem(s, "Alpha", m_alpha);
2507  }
2508 
2509  if (m_Incfreq > 0.0)
2510  {
2511  AddSummaryItem(s, "Incfreq", m_Incfreq);
2512  }
2513 
2514  if (m_MMFdir == 4)
2515  {
2516  AddSummaryItem(s, "MMFCircAxisX", m_MMFfactors[0]);
2517  AddSummaryItem(s, "MMFCircAxisY", m_MMFfactors[1]);
2518  AddSummaryItem(s, "MMFCircCentreX", m_MMFfactors[2]);
2519  AddSummaryItem(s, "MMFCircCentreY", m_MMFfactors[3]);
2520  }
2521 }
2522 } // namespace SolverUtils
2523 } // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:49
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:412
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:748
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:743
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:375
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:753
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:368
Nektar::SolverUtils::EquationSystem::m_traceNormals
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
Definition: EquationSystem.h:439
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:768
Nektar::SolverUtils::MMFSystem::GetMaxwellFlux2D
SOLVER_UTILS_EXPORT void GetMaxwellFlux2D(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1673
Nektar::SolverUtils::MMFSystem::CartesianToMovingframes
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > CartesianToMovingframes(const Array< OneD, const Array< OneD, NekDouble >> &uvec, unsigned int field)
Definition: MMFSystem.cpp:774
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFlux
SOLVER_UTILS_EXPORT void NumericalMaxwellFlux(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time=0.0)
Definition: MMFSystem.cpp:1730
Nektar::SolverUtils::MMFSystem::m_dedxi_cdot_e
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
Definition: MMFSystem.h:220
Nektar::SolverUtils::MMFSystem::m_ZimFwd
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
Definition: MMFSystem.h:207
Nektar::SolverUtils::MMFSystem::GramSchumitz
SOLVER_UTILS_EXPORT void GramSchumitz(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &outarray, bool KeepTheMagnitude=true)
Definition: MMFSystem.cpp:2405
Nektar::SolverUtils::MMFSystem::BubbleSort
SOLVER_UTILS_EXPORT void BubbleSort(Array< OneD, NekDouble > &refarray, Array< OneD, NekDouble > &sortarray)
Definition: MMFSystem.cpp:2374
Nektar::SolverUtils::MMFSystem::ComputeZimYim
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
Definition: MMFSystem.cpp:1108
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
Definition: MMFSystem.h:205
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:195
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:192
Nektar::SolverUtils::MMFSystem::VectorCrossProd
SOLVER_UTILS_EXPORT void VectorCrossProd(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, Array< OneD, NekDouble >> &v3)
Definition: MMFSystem.cpp:699
Nektar::SolverUtils::MMFSystem::GetMaxwellFluxVector
SOLVER_UTILS_EXPORT void GetMaxwellFluxVector(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1609
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2469
Nektar::SolverUtils::MMFSystem::AbsIntegral
SOLVER_UTILS_EXPORT NekDouble AbsIntegral(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2330
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
Definition: MMFSystem.h:204
Nektar::SolverUtils::MMFSystem::ComputeNtimestimesdFz
SOLVER_UTILS_EXPORT void ComputeNtimestimesdFz(const int dir, const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:1014
Nektar::SolverUtils::MMFSystem::AvgInt
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2298
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTM
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTM(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time)
Definition: MMFSystem.cpp:1807
Nektar::SolverUtils::MMFSystem::MMFSystem
SOLVER_UTILS_EXPORT MMFSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: MMFSystem.cpp:43
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFBwd
Definition: MMFSystem.h:193
Nektar::SolverUtils::MMFSystem::VectorDotProd
SOLVER_UTILS_EXPORT void VectorDotProd(const Array< OneD, const Array< OneD, NekDouble >> &v1, const Array< OneD, const Array< OneD, NekDouble >> &v2, Array< OneD, NekDouble > &v3)
Definition: MMFSystem.cpp:676
Nektar::SolverUtils::MMFSystem::CheckMovingFrames
void CheckMovingFrames(const Array< OneD, const Array< OneD, NekDouble >> &movingframes)
Definition: MMFSystem.cpp:233
Nektar::SolverUtils::MMFSystem::m_YimBwd
Array< OneD, Array< OneD, NekDouble > > m_YimBwd
Definition: MMFSystem.h:210
Nektar::SolverUtils::MMFSystem::DeriveCrossProductMF
SOLVER_UTILS_EXPORT void DeriveCrossProductMF(Array< OneD, Array< OneD, NekDouble >> &CrossProductMF)
Definition: MMFSystem.cpp:571
Nektar::SolverUtils::MMFSystem::m_muvec
Array< OneD, Array< OneD, NekDouble > > m_muvec
Definition: MMFSystem.h:213
Nektar::SolverUtils::MMFSystem::m_MMFdir
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:222
Nektar::SolverUtils::MMFSystem::AvgAbsInt
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2313
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:186
Nektar::SolverUtils::MMFSystem::m_MMFfactors
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:160
Nektar::SolverUtils::MMFSystem::m_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
Definition: MMFSystem.h:203
Nektar::SolverUtils::MMFSystem::m_MFtraceFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:199
Nektar::SolverUtils::MMFSystem::VectorAvgMagnitude
SOLVER_UTILS_EXPORT NekDouble VectorAvgMagnitude(const Array< OneD, const Array< OneD, NekDouble >> &inarray)
Definition: MMFSystem.cpp:2358
Nektar::SolverUtils::MMFSystem::m_MFtraceBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:200
Nektar::SolverUtils::MMFSystem::ComputeNtimestimesdF12
SOLVER_UTILS_EXPORT void ComputeNtimestimesdF12(const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:1080
Nektar::SolverUtils::MMFSystem::SetUpMovingFrames
void SetUpMovingFrames(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem)
Definition: MMFSystem.cpp:120
Nektar::SolverUtils::MMFSystem::GetIncidentField
SOLVER_UTILS_EXPORT Array< OneD, NekDouble > GetIncidentField(const int var, const NekDouble time)
Definition: MMFSystem.cpp:2011
Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:202
Nektar::SolverUtils::MMFSystem::NumericalMaxwellFluxTE
SOLVER_UTILS_EXPORT void NumericalMaxwellFluxTE(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd, const NekDouble time)
Definition: MMFSystem.cpp:1908
Nektar::SolverUtils::MMFSystem::AverageMaxwellFlux1D
SOLVER_UTILS_EXPORT void AverageMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1568
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
Nektar::SolverUtils::MMFSystem::~MMFSystem
virtual SOLVER_UTILS_EXPORT ~MMFSystem()
Definition: MMFSystem.cpp:49
Nektar::SolverUtils::MMFSystem::ComputeDivCurlMF
SOLVER_UTILS_EXPORT void ComputeDivCurlMF()
Definition: MMFSystem.cpp:445
Nektar::SolverUtils::MMFSystem::GetMaxwellFlux1D
SOLVER_UTILS_EXPORT void GetMaxwellFlux1D(const int var, const Array< OneD, const Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &flux)
Definition: MMFSystem.cpp:1639
Nektar::SolverUtils::MMFSystem::LaxFriedrichMaxwellFlux1D
SOLVER_UTILS_EXPORT void LaxFriedrichMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1515
Nektar::SolverUtils::MMFSystem::ComputeNtimesF12
SOLVER_UTILS_EXPORT void ComputeNtimesF12(const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &im1Fwd, const Array< OneD, const NekDouble > &im1Bwd, const Array< OneD, const NekDouble > &im2Fwd, const Array< OneD, const NekDouble > &im2Bwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:951
Nektar::SolverUtils::MMFSystem::m_ZimBwd
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
Definition: MMFSystem.h:208
Nektar::SolverUtils::MMFSystem::AdddedtMaxwell
SOLVER_UTILS_EXPORT void AdddedtMaxwell(const Array< OneD, const Array< OneD, NekDouble >> &physarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSystem.cpp:1370
Nektar::SolverUtils::MMFSystem::CopyBoundaryTrace
SOLVER_UTILS_EXPORT void CopyBoundaryTrace(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, const BoundaryCopyType BDCopyType, const int var=0, const std::string btype="NoUserDefined")
Definition: MMFSystem.cpp:838
Nektar::SolverUtils::MMFSystem::m_ncdotMFBwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFBwd
Definition: MMFSystem.h:190
Nektar::SolverUtils::MMFSystem::ComputeMFtrace
SOLVER_UTILS_EXPORT void ComputeMFtrace()
Definition: MMFSystem.cpp:522
Nektar::SolverUtils::MMFSystem::UpwindMaxwellFlux1D
SOLVER_UTILS_EXPORT void UpwindMaxwellFlux1D(Array< OneD, Array< OneD, NekDouble >> &physfield, Array< OneD, Array< OneD, NekDouble >> &numfluxFwd, Array< OneD, Array< OneD, NekDouble >> &numfluxBwd)
Definition: MMFSystem.cpp:1461
Nektar::SolverUtils::MMFSystem::ComputeNtimesFz
SOLVER_UTILS_EXPORT void ComputeNtimesFz(const int dir, const Array< OneD, Array< OneD, NekDouble >> &Fwd, const Array< OneD, Array< OneD, NekDouble >> &Bwd, const Array< OneD, const NekDouble > &imFwd, const Array< OneD, const NekDouble > &imBwd, Array< OneD, NekDouble > &outarrayFwd, Array< OneD, NekDouble > &outarrayBwd)
Definition: MMFSystem.cpp:915
Nektar::SolverUtils::MMFSystem::m_YimFwd
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
Definition: MMFSystem.h:209
Nektar::SolverUtils::MMFSystem::m_epsvec
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:212
Nektar::SolverUtils::MMFSystem::ComputencdotMF
SOLVER_UTILS_EXPORT void ComputencdotMF()
Definition: MMFSystem.cpp:320
Nektar::SolverUtils::MMFSystem::CartesianToSpherical
SOLVER_UTILS_EXPORT void CartesianToSpherical(const NekDouble x0j, const NekDouble x1j, const NekDouble x2j, NekDouble &sin_varphi, NekDouble &cos_varphi, NekDouble &sin_theta, NekDouble &cos_theta)
Definition: MMFSystem.cpp:795
Nektar::SolverUtils::MMFSystem::ComputeCurl
SOLVER_UTILS_EXPORT void ComputeCurl(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: MMFSystem.cpp:734
Nektar::SolverUtils::MMFSystem::ComputeNtimesMF
SOLVER_UTILS_EXPORT void ComputeNtimesMF()
Definition: MMFSystem.cpp:618
Nektar::SolverUtils::MMFSystem::RootMeanSquare
SOLVER_UTILS_EXPORT NekDouble RootMeanSquare(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2344
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:189
Nektar::SolverUtils::MMFSystem::m_TestMaxwellType
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:156
Nektar::SolverUtils::MMFSystem::m_CurlMF
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:196
Nektar::SolverUtils::MMFSystem::Computedemdxicdote
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Definition: MMFSystem.cpp:1280
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:49
Nektar::SolverUtils::UnsteadySystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: UnsteadySystem.cpp:610
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter/Interpreter.cpp:86
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LibUtilities::EquationSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:129
Nektar::SolverUtils::UpwindTypeMap
const char *const UpwindTypeMap[]
Definition: MMFSystem.h:89
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
Nektar::SolverUtils::SIZE_UpwindType
@ SIZE_UpwindType
Length of enum list.
Definition: MMFSystem.h:86
Nektar::SolverUtils::eAverage
@ eAverage
averaged (or centred) flux
Definition: MMFSystem.h:80
Nektar::SolverUtils::eLaxFriedrich
@ eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:81
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
Nektar::SolverUtils::SurfaceTypeMap
const char *const SurfaceTypeMap[]
Definition: MMFSystem.h:57
Nektar::SpatialDomains::GeomMMFMap
const char *const GeomMMFMap[]
Session file names associated with tangent principle directions.
Definition: SpatialDomains.hpp:89
Nektar::SpatialDomains::eLOCAL
@ eLOCAL
No Principal direction.
Definition: SpatialDomains.hpp:82
Nektar::SpatialDomains::eTangentIrregular
@ eTangentIrregular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:80
Nektar::SpatialDomains::eTangentX
@ eTangentX
X coordinate direction.
Definition: SpatialDomains.hpp:75
Nektar::SpatialDomains::eTangentCircular
@ eTangentCircular
Circular around the centre of domain.
Definition: SpatialDomains.hpp:79
Nektar::SpatialDomains::eTangentNonconvex
@ eTangentNonconvex
Circular around the centre of domain.
Definition: SpatialDomains.hpp:81
Nektar::SpatialDomains::eTangentXY
@ eTangentXY
XY direction.
Definition: SpatialDomains.hpp:77
Nektar::SpatialDomains::eTangentZ
@ eTangentZ
Z coordinate direction.
Definition: SpatialDomains.hpp:78
Nektar::SpatialDomains::eTangentY
@ eTangentY
Y coordinate direction.
Definition: SpatialDomains.hpp:76
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:553
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:1050
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
Vmath::Vvtvm
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
Definition: Vmath.cpp:598
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:918
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:284
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:384
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:945
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Sat Mar 25 2023 12:35:51 for Nektar++ by doxygen 1.9.1