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 = ( " << VectorAvgMagnitude(m_ntimesMFFwd[0])
657  << " , " << VectorAvgMagnitude(m_ntimesMFFwd[1]) << " , "
658  << VectorAvgMagnitude(m_ntimesMFFwd[2]) << " ) " << std::endl;
659  std::cout << "*** m_ntimesMFBwd = ( " << VectorAvgMagnitude(m_ntimesMFBwd[0])
660  << " , " << VectorAvgMagnitude(m_ntimesMFBwd[1]) << " , "
661  << VectorAvgMagnitude(m_ntimesMFBwd[2]) << " ) " << std::endl;
662  std::cout << "*** m_ntimes_ntimesMFFwd = ( "
663  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[0]) << " , "
664  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[1]) << " , "
665  << VectorAvgMagnitude(m_ntimes_ntimesMFFwd[2]) << " ) "
666  << std::endl;
667  std::cout << "*** m_ntimes_ntimesMFBwd = ( "
668  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[0]) << " , "
669  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[1]) << " , "
670  << VectorAvgMagnitude(m_ntimes_ntimesMFBwd[2]) << " ) "
671  << std::endl;
672 }
673 
674 void MMFSystem::VectorDotProd(
675  const Array<OneD, const Array<OneD, NekDouble>> &v1,
676  const Array<OneD, const Array<OneD, NekDouble>> &v2,
677  Array<OneD, NekDouble> &v3)
678 {
679  int coordim = v1.size();
680  int nq = v1[0].size();
681 
682  v3 = Array<OneD, NekDouble>(nq, 0.0);
683  for (int i = 0; i < coordim; ++i)
684  {
685  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &v3[0], 1, &v3[0], 1);
686  }
687 }
688 
689 /**
690  * Computes the vector cross-product in 3D of \a v1 and \a v2, storing
691  * the result in \a v3.
692  * @param v1 First input vector.
693  * @param v2 Second input vector.
694  * @param v3 Output vector computed to be orthogonal to
695  * both \a v1 and \a v2.
696  */
697 void MMFSystem::VectorCrossProd(
698  const Array<OneD, const Array<OneD, NekDouble>> &v1,
699  const Array<OneD, const Array<OneD, NekDouble>> &v2,
700  Array<OneD, Array<OneD, NekDouble>> &v3)
701 {
702  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
703  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
704  ASSERTL0(v3.size() == 3,
705  "Output vector has dimension not equal to 3.");
706 
707  int nq = v1[0].size();
708  Array<OneD, NekDouble> temp(nq);
709 
710  Vmath::Vmul(nq, v1[2], 1, v2[1], 1, temp, 1);
711  Vmath::Vvtvm(nq, v1[1], 1, v2[2], 1, temp, 1, v3[0], 1);
712 
713  Vmath::Vmul(nq, v1[0], 1, v2[2], 1, temp, 1);
714  Vmath::Vvtvm(nq, v1[2], 1, v2[0], 1, temp, 1, v3[1], 1);
715 
716  Vmath::Vmul(nq, v1[1], 1, v2[0], 1, temp, 1);
717  Vmath::Vvtvm(nq, v1[0], 1, v2[1], 1, temp, 1, v3[2], 1);
718 }
719 
720 void MMFSystem::VectorCrossProd(const Array<OneD, NekDouble> &v1,
721  const Array<OneD, NekDouble> &v2,
722  Array<OneD, NekDouble> &v3)
723 {
724  ASSERTL0(v1.size() == 3, "Input 1 has dimension not equal to 3.");
725  ASSERTL0(v2.size() == 3, "Input 2 has dimension not equal to 3.");
726  ASSERTL0(v3.size() == 3,
727  "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(
839  const Array<OneD, const NekDouble> &Fwd, Array<OneD, NekDouble> &Bwd,
840  const BoundaryCopyType BDCopyType, const int var,
841  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(
874  cnt + e));
875 
876  if (m_fields[var]->GetBndConditions()[n]->GetUserDefined() ==
877  BDtype || BDtype == "NoUserDefined")
878  {
879  switch (BDCopyType)
880  {
881  case eDirichlet:
882  {
883  Vmath::Vcopy(npts, &Dirichlet[id1], 1, &Bwd[id2], 1);
884  bdrycnt++;
885  }
886  break;
887 
888  case eFwdEQBwd:
889  {
890  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
891  bdrycnt++;
892  }
893  break;
894 
895  case eFwdEQNegBwd:
896  {
897  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
898  Vmath::Neg(npts, &Bwd[id2], 1);
899  bdrycnt++;
900  }
901  break;
902 
903  default:
904  break;
905  }
906  }
907  }
908 
909  cnt += m_fields[var]->GetBndCondExpansions()[n]->GetExpSize();
910  }
911 }
912 
913 // Compute (e^[dir] \cdot ( n \times e^3)) * (imFwd * Fwd^3 + imBwd * Bwd^3 )
914 void MMFSystem::ComputeNtimesFz(const int dir,
915  const Array<OneD, Array<OneD, NekDouble>> &Fwd,
916  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
917  const Array<OneD, const NekDouble> &imFwd,
918  const Array<OneD, const NekDouble> &imBwd,
919  Array<OneD, NekDouble> &outarrayFwd,
920  Array<OneD, NekDouble> &outarrayBwd)
921 
922 {
923  int nTraceNumPoints = GetTraceTotPoints();
924 
925  NekDouble tmpFwd, tmpBwd;
926  NekDouble Aver, ntimesz;
927 
928  for (int i = 0; i < nTraceNumPoints; ++i)
929  {
930  Aver = 0.5 * (imFwd[i] + imBwd[i]);
931 
932  tmpFwd = 0.0;
933  tmpBwd = 0.0;
934  for (int k = 0; k < m_spacedim; ++k)
935  {
936  ntimesz = 0.5 * (imFwd[i] * Fwd[2][i] + imBwd[i] * Bwd[2][i]);
937 
938  tmpFwd +=
939  m_MFtraceFwd[dir][k][i] * m_ntimesMFFwd[2][k][i] * ntimesz;
940  tmpBwd +=
941  m_MFtraceBwd[dir][k][i] * m_ntimesMFBwd[2][k][i] * ntimesz;
942  }
943 
944  outarrayFwd[i] = tmpFwd / Aver;
945  outarrayBwd[i] = tmpBwd / Aver;
946  }
947 }
948 
949 // Compute e^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd EBwd ) ) / 2{{Y_i}}
950 void MMFSystem::ComputeNtimesF12(const Array<OneD, Array<OneD, NekDouble>> &Fwd,
951  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
952  const Array<OneD, const NekDouble> &im1Fwd,
953  const Array<OneD, const NekDouble> &im1Bwd,
954  const Array<OneD, const NekDouble> &im2Fwd,
955  const Array<OneD, const NekDouble> &im2Bwd,
956  Array<OneD, NekDouble> &outarrayFwd,
957  Array<OneD, NekDouble> &outarrayBwd)
958 {
959  int nTraceNumPoints = GetTraceTotPoints();
960 
961  NekDouble tmpFwd, tmpBwd, Aver1, Aver2, HFwdk, HBwdk;
962 
963  Array<OneD, NekDouble> z1HAver(m_spacedim);
964  Array<OneD, NekDouble> z2HAver(m_spacedim);
965  Array<OneD, NekDouble> n1e1(m_spacedim);
966  Array<OneD, NekDouble> n2e2(m_spacedim);
967 
968  Array<OneD, NekDouble> n1e1_times_z1HAver(m_spacedim);
969  Array<OneD, NekDouble> n2e2_times_z2HAver(m_spacedim);
970 
971  for (int i = 0; i < nTraceNumPoints; ++i)
972  {
973  Aver1 = 0.5 * (im1Fwd[i] + im1Bwd[i]);
974  Aver2 = 0.5 * (im2Fwd[i] + im2Bwd[i]);
975 
976  for (int k = 0; k < m_spacedim; k++)
977  {
978  // Compute \vec{HFwd} and \vec{HBwd}
979  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
980  Fwd[1][i] * m_MFtraceFwd[1][k][i];
981  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
982  Bwd[1][i] * m_MFtraceBwd[1][k][i];
983 
984  // Compute z_i {{ \vec{H} }}
985  z1HAver[k] = 0.5 * (im1Fwd[i] * HFwdk + im1Bwd[i] * HBwdk);
986  z2HAver[k] = 0.5 * (im2Fwd[i] * HFwdk + im2Bwd[i] * HBwdk);
987 
988  // Choose e^i for the one in anisotropy region
989  n1e1[k] = m_ncdotMFFwd[0][i] * m_MFtraceFwd[0][k][i];
990  n2e2[k] = m_ncdotMFFwd[1][i] * m_MFtraceFwd[1][k][i];
991  }
992 
993  // Compute n1e1 \times z1HAver and n2e2 \times z2HAver
994  VectorCrossProd(n1e1, z1HAver, n1e1_times_z1HAver);
995  VectorCrossProd(n2e2, z2HAver, n2e2_times_z2HAver);
996 
997  // e^3 \cdot ( n1e1 \times z1HAver + n2e2 \times z2HAver)
998  tmpFwd = 0.0;
999  tmpBwd = 0.0;
1000  for (int k = 0; k < m_spacedim; k++)
1001  {
1002  tmpFwd += m_MFtraceFwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1003  n2e2_times_z2HAver[k] / Aver2);
1004  tmpBwd += m_MFtraceBwd[2][k][i] * (n1e1_times_z1HAver[k] / Aver1 +
1005  n2e2_times_z2HAver[k] / Aver2);
1006  }
1007 
1008  outarrayFwd[i] = tmpFwd;
1009  outarrayBwd[i] = tmpBwd;
1010  }
1011 }
1012 
1013 void MMFSystem::ComputeNtimestimesdFz(
1014  const int dir, const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1015  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1016  const Array<OneD, const NekDouble> &imFwd,
1017  const Array<OneD, const NekDouble> &imBwd,
1018  Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1019 {
1020  int nTraceNumPoints = GetTraceTotPoints();
1021 
1022  Array<OneD, NekDouble> dH(m_spacedim);
1023  Array<OneD, NekDouble> nFwd(m_spacedim);
1024 
1025  Array<OneD, NekDouble> eiFwd(m_spacedim);
1026  Array<OneD, NekDouble> eiBwd(m_spacedim);
1027 
1028  Array<OneD, NekDouble> eitimesdHFwd(m_spacedim);
1029  Array<OneD, NekDouble> eitimesdHBwd(m_spacedim);
1030 
1031  Array<OneD, NekDouble> ntimeseitimesdHFwd(m_spacedim);
1032  Array<OneD, NekDouble> ntimeseitimesdHBwd(m_spacedim);
1033 
1034  NekDouble Aver, HFwdk, HBwdk, tmpFwd, tmpBwd;
1035  for (int i = 0; i < nTraceNumPoints; ++i)
1036  {
1037  Aver = 0.5 * (imFwd[i] + imBwd[i]);
1038 
1039  // Get [H]
1040  for (int k = 0; k < m_spacedim; k++)
1041  {
1042  HFwdk = Fwd[0][i] * m_MFtraceFwd[0][k][i] +
1043  Fwd[1][i] * m_MFtraceFwd[1][k][i];
1044  HBwdk = Bwd[0][i] * m_MFtraceBwd[0][k][i] +
1045  Bwd[1][i] * m_MFtraceBwd[1][k][i];
1046  dH[k] = HFwdk - HBwdk;
1047 
1048  eiFwd[k] = m_MFtraceFwd[dir][k][i];
1049  eiBwd[k] = m_MFtraceBwd[dir][k][i];
1050 
1051  nFwd[k] = m_traceNormals[k][i];
1052  }
1053 
1054  // MFtraceFwd (MFtraceBwd) \times [H]
1055  // VectorCrossProd(eiFwd, dH, eitimesdHFwd);
1056  // VectorCrossProd(eiBwd, dH, eitimesdHBwd);
1057  VectorCrossProd(nFwd, dH, eitimesdHFwd);
1058  VectorCrossProd(nFwd, dH, eitimesdHBwd);
1059 
1060  // n times eitimesdH
1061  VectorCrossProd(nFwd, eitimesdHFwd, ntimeseitimesdHFwd);
1062  VectorCrossProd(nFwd, eitimesdHBwd, ntimeseitimesdHBwd);
1063 
1064  // MFtraceFwd \cdot ntimeseitimesdH
1065  tmpFwd = 0.0;
1066  tmpBwd = 0.0;
1067  for (int k = 0; k < m_spacedim; k++)
1068  {
1069  tmpFwd += eiFwd[k] * ntimeseitimesdHFwd[k];
1070  tmpBwd += eiBwd[k] * ntimeseitimesdHBwd[k];
1071  }
1072 
1073  outarrayFwd[i] = 0.5 * m_alpha * tmpFwd / Aver;
1074  outarrayBwd[i] = 0.5 * m_alpha * tmpBwd / Aver;
1075  }
1076 }
1077 
1078 // Compute - \alpha [E3] / ( {{im1}} + {{im2}} )
1079 void MMFSystem::ComputeNtimestimesdF12(
1080  const Array<OneD, Array<OneD, NekDouble>> &Fwd,
1081  const Array<OneD, Array<OneD, NekDouble>> &Bwd,
1082  const Array<OneD, const NekDouble> &im1Fwd,
1083  const Array<OneD, const NekDouble> &im1Bwd,
1084  const Array<OneD, const NekDouble> &im2Fwd,
1085  const Array<OneD, const NekDouble> &im2Bwd,
1086  Array<OneD, NekDouble> &outarrayFwd, Array<OneD, NekDouble> &outarrayBwd)
1087 {
1088  int nTraceNumPoints = GetTraceTotPoints();
1089 
1090  Array<OneD, NekDouble> directFwd(nTraceNumPoints);
1091  Array<OneD, NekDouble> directBwd(nTraceNumPoints);
1092 
1093  NekDouble Aver1, Aver2;
1094  for (int i = 0; i < nTraceNumPoints; ++i)
1095  {
1096  Aver1 = im1Fwd[i] + im1Bwd[i];
1097  Aver2 = im2Fwd[i] + im2Bwd[i];
1098 
1099  outarrayFwd[i] =
1100  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1101  outarrayBwd[i] =
1102  -m_alpha * (Fwd[2][i] - Bwd[2][i]) * (1.0 / Aver1 + 1.0 / Aver2);
1103  }
1104 }
1105 
1106 // dim = dimension, pol = polarization, 0 = TM, 1 = TE.
1107 void MMFSystem::ComputeZimYim(Array<OneD, Array<OneD, NekDouble>> &epsvec,
1108  Array<OneD, Array<OneD, NekDouble>> &muvec)
1109 
1110 {
1111  int nTraceNumPoints = GetTraceNpoints();
1112 
1113  switch (m_TestMaxwellType)
1114  {
1115  case eMaxwell1D:
1116  case eScatField1D:
1117  {
1118  Array<OneD, NekDouble> Fwdeps(nTraceNumPoints, 1.0);
1119  Array<OneD, NekDouble> Bwdeps(nTraceNumPoints, 1.0);
1120  Array<OneD, NekDouble> Fwdmu(nTraceNumPoints, 1.0);
1121  Array<OneD, NekDouble> Bwdmu(nTraceNumPoints, 1.0);
1122  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps, Bwdeps);
1123  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdeps, Bwdeps);
1124 
1125  CopyBoundaryTrace(Fwdeps, Bwdeps, eFwdEQBwd, 0);
1126  CopyBoundaryTrace(Fwdmu, Bwdmu, eFwdEQBwd, 1);
1127 
1128  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1129  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1130  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(1);
1131  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(1);
1132 
1133  m_ZimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1134  m_ZimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1135  m_YimFwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1136  m_YimBwd[0] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1137 
1138  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd / epsBwd)
1139  for (int i = 0; i < nTraceNumPoints; ++i)
1140  {
1141  m_ZimFwd[0][i] = sqrt(Fwdmu[i] / Fwdeps[i]);
1142  m_ZimBwd[0][i] = sqrt(Bwdmu[i] / Bwdeps[i]);
1143 
1144  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1145  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1146  }
1147 
1148  std::cout << "*** ZimFwd = " << RootMeanSquare(m_ZimFwd[0])
1149  << ", ZimBwd = " << RootMeanSquare(m_ZimBwd[0])
1150  << ", YimFwd = " << RootMeanSquare(m_YimFwd[0])
1151  << ", YimBwd = " << RootMeanSquare(m_YimBwd[0])
1152  << std::endl;
1153  }
1154  break;
1155 
1156  case eTestMaxwell2DPEC:
1157  case eTestMaxwell2DPECAVGFLUX:
1158  case eTestMaxwell2DPMC:
1159  case eScatField2D:
1160  case eTotField2D:
1161  case eMaxwellSphere:
1162  case eELF2DSurface:
1163  {
1164  m_ZimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1165  m_ZimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1166  m_YimFwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1167  m_YimBwd = Array<OneD, Array<OneD, NekDouble>>(m_shapedim);
1168 
1169  for (int j = 0; j < m_shapedim; ++j)
1170  {
1171  m_ZimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1172  m_ZimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1173  m_YimFwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1174  m_YimBwd[j] = Array<OneD, NekDouble>(nTraceNumPoints, 1.0);
1175  }
1176 
1177  switch (m_PolType)
1178  {
1179  case eTransMagnetic:
1180  {
1181  Array<OneD, NekDouble> Fwdmu1(nTraceNumPoints);
1182  Array<OneD, NekDouble> Bwdmu1(nTraceNumPoints);
1183  Array<OneD, NekDouble> Fwdmu2(nTraceNumPoints);
1184  Array<OneD, NekDouble> Bwdmu2(nTraceNumPoints);
1185  Array<OneD, NekDouble> Fwdeps3(nTraceNumPoints);
1186  Array<OneD, NekDouble> Bwdeps3(nTraceNumPoints);
1187 
1188  m_fields[0]->GetFwdBwdTracePhys(muvec[0], Fwdmu1, Bwdmu1);
1189  m_fields[0]->GetFwdBwdTracePhys(muvec[1], Fwdmu2, Bwdmu2);
1190  m_fields[0]->GetFwdBwdTracePhys(epsvec[2], Fwdeps3,
1191  Bwdeps3);
1192 
1193  CopyBoundaryTrace(Fwdmu1, Bwdmu1, eFwdEQBwd, 0);
1194  CopyBoundaryTrace(Fwdmu2, Bwdmu2, eFwdEQBwd, 1);
1195  CopyBoundaryTrace(Fwdeps3, Bwdeps3, eFwdEQBwd, 2);
1196 
1197  // ZimFwd = sqrt( muFwd / epsFwd), ZimBwd = sqrt( muBwd /
1198  // epsBwd)
1199  for (int i = 0; i < nTraceNumPoints; ++i)
1200  {
1201  m_ZimFwd[0][i] = sqrt(Fwdmu2[i] / Fwdeps3[i]);
1202  m_ZimBwd[0][i] = sqrt(Bwdmu2[i] / Bwdeps3[i]);
1203 
1204  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1205  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1206 
1207  m_ZimFwd[1][i] = sqrt(Fwdmu1[i] / Fwdeps3[i]);
1208  m_ZimBwd[1][i] = sqrt(Bwdmu1[i] / Bwdeps3[i]);
1209 
1210  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1211  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1212  }
1213  }
1214  break; // eTransMagnetic
1215 
1216  case eTransElectric:
1217  {
1218  Array<OneD, NekDouble> Fwdeps1(nTraceNumPoints);
1219  Array<OneD, NekDouble> Bwdeps1(nTraceNumPoints);
1220  Array<OneD, NekDouble> Fwdeps2(nTraceNumPoints);
1221  Array<OneD, NekDouble> Bwdeps2(nTraceNumPoints);
1222  Array<OneD, NekDouble> Fwdmu3(nTraceNumPoints);
1223  Array<OneD, NekDouble> Bwdmu3(nTraceNumPoints);
1224 
1225  m_fields[0]->GetFwdBwdTracePhys(epsvec[0], Fwdeps1,
1226  Bwdeps1);
1227  m_fields[0]->GetFwdBwdTracePhys(epsvec[1], Fwdeps2,
1228  Bwdeps2);
1229  m_fields[0]->GetFwdBwdTracePhys(muvec[2], Fwdmu3, Bwdmu3);
1230 
1231  CopyBoundaryTrace(Fwdeps1, Bwdeps1, eFwdEQBwd, 0);
1232  CopyBoundaryTrace(Fwdeps2, Bwdeps2, eFwdEQBwd, 1);
1233  CopyBoundaryTrace(Fwdmu3, Bwdmu3, eFwdEQBwd, 2);
1234 
1235  for (int i = 0; i < nTraceNumPoints; ++i)
1236  {
1237  m_ZimFwd[0][i] = sqrt(Fwdmu3[i] / Fwdeps2[i]);
1238  m_ZimBwd[0][i] = sqrt(Bwdmu3[i] / Bwdeps2[i]);
1239 
1240  m_YimFwd[0][i] = 1.0 / m_ZimFwd[0][i];
1241  m_YimBwd[0][i] = 1.0 / m_ZimBwd[0][i];
1242 
1243  m_ZimFwd[1][i] = sqrt(Fwdmu3[i] / Fwdeps1[i]);
1244  m_ZimBwd[1][i] = sqrt(Bwdmu3[i] / Bwdeps1[i]);
1245 
1246  m_YimFwd[1][i] = 1.0 / m_ZimFwd[1][i];
1247  m_YimBwd[1][i] = 1.0 / m_ZimBwd[1][i];
1248  }
1249  }
1250  break; // eTransELectric
1251 
1252  default:
1253  break;
1254  } // PolType
1255 
1256  std::cout << "*** ZimFwd0 = [ "
1257  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[0], 1) << " , "
1258  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[0], 1)
1259  << " ], ZimBwd0 = [ "
1260  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[0], 1) << " , "
1261  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[0], 1) << " ] "
1262  << std::endl;
1263  std::cout << "*** ZimFwd1 = [ "
1264  << Vmath::Vmin(nTraceNumPoints, m_ZimFwd[1], 1) << " , "
1265  << Vmath::Vmax(nTraceNumPoints, m_ZimFwd[1], 1)
1266  << " ], ZimBwd1 = [ "
1267  << Vmath::Vmin(nTraceNumPoints, m_ZimBwd[1], 1) << " , "
1268  << Vmath::Vmax(nTraceNumPoints, m_ZimBwd[1], 1) << " ] "
1269  << std::endl;
1270  }
1271  break; // eMaxwell2D
1272 
1273  default:
1274  break;
1275  } // TestMaxwellType
1276 }
1277 
1278 // m_dedxi_cdot_e[m][j][k][] = de^m / d \xi^j \cdot e^n
1279 void MMFSystem::Computedemdxicdote()
1280 {
1281  int MFdim = 3;
1282  int nq = GetTotPoints();
1283 
1284  // Initialization
1285  m_dedxi_cdot_e =
1286  Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble>>>>(MFdim);
1287  for (int indm = 0; indm < MFdim; ++indm)
1288  {
1289  m_dedxi_cdot_e[indm] =
1290  Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(MFdim);
1291  for (int indj = 0; indj < MFdim; ++indj)
1292  {
1293  m_dedxi_cdot_e[indm][indj] =
1294  Array<OneD, Array<OneD, NekDouble>>(MFdim);
1295  for (int indn = 0; indn < MFdim; ++indn)
1296  {
1297  m_dedxi_cdot_e[indm][indj][indn] =
1298  Array<OneD, NekDouble>(nq, 0.0);
1299  }
1300  }
1301  }
1302 
1303  Array<OneD, NekDouble> tmp(nq);
1304  Array<OneD, NekDouble> tmpderiv(nq);
1305  Array<OneD, NekDouble> dedt;
1306  for (int indm = 0; indm < MFdim; ++indm)
1307  {
1308  for (int indj = 0; indj < MFdim; ++indj)
1309  {
1310  for (int indn = 0; indn < MFdim; ++indn)
1311  {
1312  dedt = Array<OneD, NekDouble>(nq, 0.0);
1313  for (int k = 0; k < m_spacedim; ++k)
1314  {
1315  // Compute d e^m / d \xi_j cdot e^n
1316  Vmath::Vcopy(nq, &m_movingframes[indm][k * nq], 1, &tmp[0],
1317  1);
1318  m_fields[0]->PhysDirectionalDeriv(m_movingframes[indj], tmp,
1319  tmpderiv);
1320 
1321  Vmath::Vvtvp(nq, &tmpderiv[0], 1,
1322  &m_movingframes[indn][k * nq], 1, &dedt[0], 1,
1323  &dedt[0], 1);
1324  }
1325 
1326  Vmath::Vcopy(nq, &dedt[0], 1,
1327  &m_dedxi_cdot_e[indm][indj][indn][0], 1);
1328  }
1329  }
1330  }
1331 
1332  int indx = 0;
1333  std::cout << "*** m_dedxi_cdot_e[0]/dxi1 = ( "
1334  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1335  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1336  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1337  << std::endl;
1338  std::cout << "*** m_dedxi_cdot_e[0]/dxi2 = ( "
1339  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1340  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1341  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1342  << std::endl;
1343 
1344  indx = 1;
1345  std::cout << "*** m_dedxi_cdot_e[1]/dxi1 = ( "
1346  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1347  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1348  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1349  << std::endl;
1350  std::cout << "*** m_dedxi_cdot_e[1]/dxi2 = ( "
1351  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1352  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1353  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1354  << std::endl;
1355 
1356  indx = 2;
1357  std::cout << "*** m_dedxi_cdot_e[2]/dxi1 = ( "
1358  << RootMeanSquare(m_dedxi_cdot_e[indx][0][0]) << " , "
1359  << RootMeanSquare(m_dedxi_cdot_e[indx][0][1]) << " , "
1360  << RootMeanSquare(m_dedxi_cdot_e[indx][0][2]) << " )_1, "
1361  << std::endl;
1362  std::cout << "*** m_dedxi_cdot_e[2]/dxi2 = ( "
1363  << RootMeanSquare(m_dedxi_cdot_e[indx][1][0]) << " , "
1364  << RootMeanSquare(m_dedxi_cdot_e[indx][1][1]) << " , "
1365  << RootMeanSquare(m_dedxi_cdot_e[indx][1][2]) << " )_2 "
1366  << std::endl;
1367 }
1368 
1369 void MMFSystem::AdddedtMaxwell(
1370  const Array<OneD, const Array<OneD, NekDouble>> &physarray,
1371  Array<OneD, Array<OneD, NekDouble>> &outarray)
1372 {
1373  int nq = GetTotPoints();
1374 
1375  // m_dedxi_cdot_e[m][j][n][] = de^m / d \xi^j \cdot e^n
1376  Array<OneD, NekDouble> dedtej(nq);
1377  Array<OneD, NekDouble> de1dtcdotej(nq);
1378  Array<OneD, NekDouble> de2dtcdotej(nq);
1379 
1380  Array<OneD, NekDouble> normH(nq);
1381  Array<OneD, NekDouble> NormedH1(nq, 0.0);
1382  Array<OneD, NekDouble> NormednegH2(nq, 0.0);
1383 
1384  Vmath::Vmul(nq, &physarray[0][0], 1, &physarray[0][0], 1, &normH[0], 1);
1385  Vmath::Vvtvp(nq, &physarray[1][0], 1, &physarray[1][0], 1, &normH[0], 1,
1386  &normH[0], 1);
1387  Vmath::Vsqrt(nq, normH, 1, normH, 1);
1388 
1389  NekDouble Tol = 0.001;
1390  for (int i = 0; i < nq; ++i)
1391  {
1392  if (normH[i] > Tol)
1393  {
1394  NormedH1[i] = physarray[0][i] / normH[i];
1395  NormednegH2[i] = -1.0 * physarray[1][i] / normH[i];
1396  }
1397  }
1398 
1399  for (int j = 0; j < m_shapedim; ++j)
1400  {
1401  // Compute de1 / dt \cdot ej = (-H2 de^1/d\xi1 \cdot e^j + H1 de^1/d\xi2
1402  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1403  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[0][0][j][0], 1,
1404  &de1dtcdotej[0], 1);
1405  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[0][1][j][0], 1,
1406  &de1dtcdotej[0], 1, &de1dtcdotej[0], 1);
1407 
1408  // Compute de2 / dt \cdot ej = (-H2 de2/d\xi1 \cdot e^j + H1 de2/d\xi2
1409  // \cdot e^j) / sqrt{ H1^2 + H2^2 }
1410  Vmath::Vmul(nq, &NormednegH2[0], 1, &m_dedxi_cdot_e[1][0][j][0], 1,
1411  &de2dtcdotej[0], 1);
1412  Vmath::Vvtvp(nq, &NormedH1[0], 1, &m_dedxi_cdot_e[1][1][j][0], 1,
1413  &de2dtcdotej[0], 1, &de2dtcdotej[0], 1);
1414 
1415  // Add dedt component: (H1 (de1/dt) + H2 (de2/dt) ) \cdot ej
1416  Vmath::Vmul(nq, &physarray[0][0], 1, &de1dtcdotej[0], 1, &dedtej[0], 1);
1417  Vmath::Vvtvp(nq, &physarray[1][0], 1, &de2dtcdotej[0], 1, &dedtej[0], 1,
1418  &dedtej[0], 1);
1419 
1420  Vmath::Neg(nq, dedtej, 1);
1421 
1422  switch (m_PolType)
1423  {
1424  case SolverUtils::eTransMagnetic:
1425  {
1426  if (j == 0)
1427  {
1428  Vmath::Vmul(nq, m_muvec[0], 1, dedtej, 1, dedtej, 1);
1429  }
1430 
1431  else if (j == 1)
1432  {
1433  Vmath::Vmul(nq, m_muvec[1], 1, dedtej, 1, dedtej, 1);
1434  }
1435  }
1436  break;
1437 
1438  case SolverUtils::eTransElectric:
1439  {
1440  if (j == 0)
1441  {
1442  Vmath::Vmul(nq, m_epsvec[0], 1, dedtej, 1, dedtej, 1);
1443  }
1444 
1445  else if (j == 1)
1446  {
1447  Vmath::Vmul(nq, m_epsvec[1], 1, dedtej, 1, dedtej, 1);
1448  }
1449  }
1450  break;
1451 
1452  default:
1453  break;
1454  }
1455 
1456  Vmath::Vadd(nq, &dedtej[0], 1, &outarray[j][0], 1, &outarray[j][0], 1);
1457  }
1458 }
1459 
1460 void MMFSystem::UpwindMaxwellFlux1D(
1461  Array<OneD, Array<OneD, NekDouble>> &physfield,
1462  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1463  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1464 {
1465  int i;
1466  int nTraceNumPoints = GetTraceTotPoints();
1467  int nvar = 2;
1468 
1469  // get temporary arrays
1470  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1471  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1472 
1473  for (i = 0; i < nvar; ++i)
1474  {
1475  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1476  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1477  }
1478 
1479  // get the physical values at the trace from the dependent variables
1480  for (i = 0; i < nvar; ++i)
1481  {
1482  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1483  }
1484 
1485  // E = 0 at the boundaries
1486  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1487  "PEC");
1488 
1489  // d H / d n = 0 at the boundaries
1490  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1491  "PEC");
1492 
1493  Array<OneD, NekDouble> dE(nTraceNumPoints);
1494  Array<OneD, NekDouble> dH(nTraceNumPoints);
1495 
1496  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1497  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1498 
1499  NekDouble nx, AverZ, AverY, AverZH, AverYE;
1500  for (i = 0; i < nTraceNumPoints; ++i)
1501  {
1502  nx = m_traceNormals[0][i];
1503  AverZ = m_ZimFwd[0][i] + m_ZimBwd[0][i];
1504  AverY = m_YimFwd[0][i] + m_YimBwd[0][i];
1505  AverZH = m_ZimFwd[0][i] * Fwd[1][i] + m_ZimBwd[0][i] * Bwd[1][i];
1506  AverYE = m_YimFwd[0][i] * Fwd[0][i] + m_YimBwd[0][i] * Bwd[0][i];
1507 
1508  numfluxFwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1509  numfluxFwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1510 
1511  numfluxBwd[0][i] = nx / AverZ * (AverZH - nx * dE[i]);
1512  numfluxBwd[1][i] = nx / AverY * (AverYE - nx * dH[i]);
1513  }
1514 }
1515 
1516 void MMFSystem::LaxFriedrichMaxwellFlux1D(
1517  Array<OneD, Array<OneD, NekDouble>> &physfield,
1518  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1519  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1520 {
1521  int i;
1522  int nTraceNumPoints = GetTraceTotPoints();
1523  int nvar = 2;
1524 
1525  // get temporary arrays
1526  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1527  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1528 
1529  for (i = 0; i < nvar; ++i)
1530  {
1531  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1532  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1533  }
1534 
1535  // get the physical values at the trace from the dependent variables
1536  for (i = 0; i < nvar; ++i)
1537  {
1538  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1539  }
1540 
1541  // E = 0 at the boundaries
1542  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1543  "PEC");
1544 
1545  // d H / d n = 0 at the boundaries
1546  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1547  "PEC");
1548 
1549  Array<OneD, NekDouble> dE(nTraceNumPoints);
1550  Array<OneD, NekDouble> dH(nTraceNumPoints);
1551 
1552  Vmath::Vsub(nTraceNumPoints, &Fwd[0][0], 1, &Bwd[0][0], 1, &dE[0], 1);
1553  Vmath::Vsub(nTraceNumPoints, &Fwd[1][0], 1, &Bwd[1][0], 1, &dH[0], 1);
1554 
1555  NekDouble nx;
1556  for (i = 0; i < nTraceNumPoints; ++i)
1557  {
1558  nx = m_traceNormals[0][i];
1559  numfluxFwd[0][i] =
1560  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1561  numfluxFwd[1][i] =
1562  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1563 
1564  numfluxBwd[0][i] =
1565  0.5 * nx * ((Fwd[1][i] + Bwd[1][i]) - nx * m_YimFwd[0][i] * dE[i]);
1566  numfluxBwd[1][i] =
1567  0.5 * nx * ((Fwd[0][i] + Bwd[0][i]) - nx * m_ZimFwd[0][i] * dH[i]);
1568  }
1569 }
1570 
1571 void MMFSystem::AverageMaxwellFlux1D(
1572  Array<OneD, Array<OneD, NekDouble>> &physfield,
1573  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1574  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd)
1575 {
1576  int i;
1577  int nTraceNumPoints = GetTraceTotPoints();
1578  int nvar = 2;
1579 
1580  // get temporary arrays
1581  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1582  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1583 
1584  for (i = 0; i < nvar; ++i)
1585  {
1586  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1587  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints);
1588  }
1589 
1590  // get the physical values at the trace from the dependent variables
1591  for (i = 0; i < nvar; ++i)
1592  {
1593  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1594  }
1595 
1596  // E = 0 at the boundaries
1597  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1598  "PEC");
1599 
1600  // d H / d n = 0 at the boundaries
1601  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1602  "PEC");
1603 
1604  for (i = 0; i < nTraceNumPoints; ++i)
1605  {
1606  numfluxFwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1607  numfluxFwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1608 
1609  numfluxBwd[0][i] = 0.5 * m_traceNormals[0][i] * (Fwd[1][i] + Bwd[1][i]);
1610  numfluxBwd[1][i] = 0.5 * m_traceNormals[0][i] * (Fwd[0][i] + Bwd[0][i]);
1611  }
1612 }
1613 
1614 void MMFSystem::GetMaxwellFluxVector(
1615  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1616  Array<OneD, Array<OneD, NekDouble>> &flux)
1617 {
1618  switch (m_TestMaxwellType)
1619  {
1620  case eMaxwell1D:
1621  case eScatField1D:
1622  {
1623  GetMaxwellFlux1D(var, physfield, flux);
1624  }
1625  break;
1626 
1627  case eTestMaxwell2DPEC:
1628  case eTestMaxwell2DPECAVGFLUX:
1629  case eTestMaxwell2DPMC:
1630  case eScatField2D:
1631  case eTotField2D:
1632  case eMaxwellSphere:
1633  case eELF2DSurface:
1634  {
1635  GetMaxwellFlux2D(var, physfield, flux);
1636  }
1637  break;
1638 
1639  default:
1640  break;
1641  }
1642 }
1643 
1644 void MMFSystem::GetMaxwellFlux1D(
1645  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1646  Array<OneD, Array<OneD, NekDouble>> &flux)
1647 {
1648  int nq = m_fields[0]->GetTotPoints();
1649 
1650  switch (var)
1651  {
1652  case 0:
1653  {
1654  // H in flux 0
1655  Vmath::Vcopy(nq, physfield[1], 1, flux[0], 1);
1656 
1657  // E in flux 1
1658  Vmath::Zero(nq, flux[1], 1);
1659  }
1660  break;
1661 
1662  case 1:
1663  {
1664  // E in flux 0
1665  Vmath::Vcopy(nq, physfield[0], 1, flux[0], 1);
1666 
1667  // H in flux 1
1668  Vmath::Zero(nq, flux[1], 1);
1669  }
1670  break;
1671  //----------------------------------------------------
1672 
1673  default:
1674  break;
1675  }
1676 }
1677 
1678 void MMFSystem::GetMaxwellFlux2D(
1679  const int var, const Array<OneD, const Array<OneD, NekDouble>> &physfield,
1680  Array<OneD, Array<OneD, NekDouble>> &flux)
1681 {
1682  int nq = m_fields[0]->GetTotPoints();
1683 
1684  NekDouble sign = 1.0;
1685  switch (m_PolType)
1686  {
1687  // TransMagnetic
1688  case 0:
1689  {
1690  sign = -1.0;
1691  }
1692  break;
1693 
1694  // TransElectric
1695  case 1:
1696  {
1697  sign = 1.0;
1698  }
1699  break;
1700 
1701  default:
1702  break;
1703  }
1704 
1705  switch (var)
1706  {
1707  case 0:
1708  {
1709  // -Ez in flux 1
1710  Vmath::Smul(nq, sign, physfield[2], 1, flux[0], 1);
1711  Vmath::Zero(nq, flux[1], 1);
1712  }
1713  break;
1714 
1715  case 1:
1716  {
1717  // Ez in flux 0
1718  Vmath::Zero(nq, flux[0], 1);
1719  Vmath::Smul(nq, -sign, physfield[2], 1, flux[1], 1);
1720  }
1721  break;
1722 
1723  case 2:
1724  {
1725  Vmath::Smul(nq, sign, physfield[0], 1, flux[0], 1);
1726  Vmath::Smul(nq, -sign, physfield[1], 1, flux[1], 1);
1727  }
1728  break;
1729 
1730  default:
1731  ASSERTL0(false, "GetFluxVector2D: illegal vector index");
1732  }
1733 }
1734 
1735 void MMFSystem::NumericalMaxwellFlux(
1736  Array<OneD, Array<OneD, NekDouble>> &physfield,
1737  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1738  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1739 {
1740 
1741  switch (m_TestMaxwellType)
1742  {
1743  case eMaxwell1D:
1744  case eScatField1D:
1745  {
1746  switch (m_upwindType)
1747  {
1748  case eAverage:
1749  {
1750  AverageMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1751  }
1752  break;
1753 
1754  case eLaxFriedrich:
1755  {
1756  LaxFriedrichMaxwellFlux1D(physfield, numfluxFwd,
1757  numfluxBwd);
1758  }
1759  break;
1760 
1761  case eUpwind:
1762  {
1763  UpwindMaxwellFlux1D(physfield, numfluxFwd, numfluxBwd);
1764  }
1765  break;
1766 
1767  default:
1768  {
1769  ASSERTL0(false,
1770  "populate switch statement for upwind flux");
1771  }
1772  break; // upwindType
1773  }
1774  }
1775  break; // eMaxwell1D
1776 
1777  case eTestMaxwell2DPEC:
1778  case eTestMaxwell2DPECAVGFLUX:
1779  case eTestMaxwell2DPMC:
1780  case eScatField2D:
1781  case eTotField2D:
1782  case eMaxwellSphere:
1783  case eELF2DSurface:
1784  {
1785  switch (m_PolType)
1786  {
1787  case eTransMagnetic:
1788  {
1789  NumericalMaxwellFluxTM(physfield, numfluxFwd, numfluxBwd,
1790  time);
1791  }
1792  break;
1793 
1794  case eTransElectric:
1795  {
1796  NumericalMaxwellFluxTE(physfield, numfluxFwd, numfluxBwd,
1797  time);
1798  }
1799  break;
1800 
1801  default:
1802  break;
1803  }
1804  }
1805  break; // eMaxwell2D
1806 
1807  default:
1808  break;
1809  } // m_TestMaxwellType
1810 }
1811 
1812 void MMFSystem::NumericalMaxwellFluxTM(
1813  Array<OneD, Array<OneD, NekDouble>> &physfield,
1814  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1815  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1816 {
1817  int nq = m_fields[0]->GetNpoints();
1818  int nTraceNumPoints = GetTraceTotPoints();
1819  int nvar = physfield.size();
1820 
1821  // get temporary arrays
1822  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1823  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1824  for (int i = 0; i < nvar; ++i)
1825  {
1826  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1827  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1828 
1829  // get the physical values at the trace
1830  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1831  }
1832 
1833  // E^|| = 0 at the PEC boundaries vs. H^|| = 0 at PMC boundaries
1834  Array<OneD, NekDouble> IncField(nq, 0.0);
1835  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1836  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1837  for (int i = 0; i < m_spacedim; ++i)
1838  {
1839  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1840 
1841  IncField = GetIncidentField(i, time);
1842  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1843 
1844  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1845  &IncFieldFwd[i][0], 1);
1846  }
1847 
1848  // Total Field Formulation
1849  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0,
1850  "PEC_Forces");
1851  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1852  "PEC_Forces");
1853  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1854  "PEC_Forces");
1855 
1856  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1857  "PMC_Forces");
1858  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1859  "PMC_Forces");
1860  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2,
1861  "PMC_Forces");
1862 
1863  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0,
1864  "PEC");
1865  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1866  "PEC");
1867  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1868  "PEC");
1869 
1870  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1871  "PMC");
1872  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1873  "PMC");
1874  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2,
1875  "PMC");
1876 
1877  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1878  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1879  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1880  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdH(nTraceNumPoints, 0.0);
1881  Array<OneD, NekDouble> e3Fwd_cdot_dEe3(nTraceNumPoints, 0.0);
1882  Array<OneD, NekDouble> e3Bwd_cdot_dEe3(nTraceNumPoints, 0.0);
1883 
1884  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (YimFwd *
1885  // EFwd^3 + YimBwd * EBwd^3 )
1886  ComputeNtimesFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], numfluxFwd[0],
1887  numfluxBwd[0]);
1888  ComputeNtimesFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1], numfluxFwd[1],
1889  numfluxBwd[1]);
1890 
1891  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( ZimFwd HFwd + ZimBwd
1892  // HBwd ) ) / 2 {{Z_i}}
1893  ComputeNtimesF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1894  m_ZimBwd[1], numfluxFwd[2], numfluxBwd[2]);
1895 
1896  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
1897  // [H] / 2 {{YimFwd}}
1898  ComputeNtimestimesdFz(0, Fwd, Bwd, m_YimFwd[0], m_YimBwd[0],
1899  e1Fwd_cdot_ncrossdH, e1Bwd_cdot_ncrossdH);
1900  ComputeNtimestimesdFz(1, Fwd, Bwd, m_YimFwd[1], m_YimBwd[1],
1901  e2Fwd_cdot_ncrossdH, e2Bwd_cdot_ncrossdH);
1902 
1903  // Compute \alpha [E3] * ( 1/2{{Zim1}} + 1/2{{Zim2}} )
1904  ComputeNtimestimesdF12(Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], m_ZimFwd[1],
1905  m_ZimBwd[1], e3Fwd_cdot_dEe3, e3Bwd_cdot_dEe3);
1906 
1907  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdH,
1908  1, numfluxFwd[0], 1);
1909  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdH,
1910  1, numfluxFwd[1], 1);
1911  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[2], 1, e3Fwd_cdot_dEe3, 1,
1912  numfluxFwd[2], 1);
1913 
1914  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdH,
1915  1, numfluxBwd[0], 1);
1916  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdH,
1917  1, numfluxBwd[1], 1);
1918  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[2], 1, e3Bwd_cdot_dEe3, 1,
1919  numfluxBwd[2], 1);
1920 }
1921 
1922 void MMFSystem::NumericalMaxwellFluxTE(
1923  Array<OneD, Array<OneD, NekDouble>> &physfield,
1924  Array<OneD, Array<OneD, NekDouble>> &numfluxFwd,
1925  Array<OneD, Array<OneD, NekDouble>> &numfluxBwd, const NekDouble time)
1926 
1927 {
1928  int nq = m_fields[0]->GetNpoints();
1929  int nTraceNumPoints = GetTraceTotPoints();
1930  int nvar = physfield.size();
1931 
1932  // Get temporary arrays
1933  Array<OneD, Array<OneD, NekDouble>> Fwd(nvar);
1934  Array<OneD, Array<OneD, NekDouble>> Bwd(nvar);
1935  for (int i = 0; i < nvar; ++i)
1936  {
1937  Fwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1938  Bwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1939 
1940  // get the physical values at the trace
1941  m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd[i], Bwd[i]);
1942  }
1943 
1944  // E = 0 at the PEC boundaries:
1945  Array<OneD, NekDouble> IncField(nq, 0.0);
1946  Array<OneD, NekDouble> IncFieldBwd(nTraceNumPoints, 0.0);
1947  Array<OneD, Array<OneD, NekDouble>> IncFieldFwd(m_spacedim);
1948  for (int i = 0; i < m_spacedim; ++i)
1949  {
1950  IncFieldFwd[i] = Array<OneD, NekDouble>(nTraceNumPoints, 0.0);
1951 
1952  IncField = GetIncidentField(i, time);
1953  m_fields[0]->GetFwdBwdTracePhys(IncField, IncFieldFwd[i], IncFieldBwd);
1954 
1955  Vmath::Svtvp(nTraceNumPoints, 2.0, &IncFieldFwd[i][0], 1, &Fwd[i][0], 1,
1956  &IncFieldFwd[i][0], 1);
1957  }
1958 
1959  CopyBoundaryTrace(IncFieldFwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1960  "PEC_Forces");
1961  CopyBoundaryTrace(IncFieldFwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1962  "PEC_Forces");
1963  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2,
1964  "PEC_Forces");
1965 
1966  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0,
1967  "PMC_Forces");
1968  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1969  "PMC_Forces");
1970  CopyBoundaryTrace(IncFieldFwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1971  "PMC_Forces");
1972 
1973  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQNegBwd, 0,
1974  "PEC");
1975  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQNegBwd, 1,
1976  "PEC");
1977  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQBwd, 2,
1978  "PEC");
1979 
1980  CopyBoundaryTrace(Fwd[0], Bwd[0], SolverUtils::eFwdEQBwd, 0,
1981  "PMC");
1982  CopyBoundaryTrace(Fwd[1], Bwd[1], SolverUtils::eFwdEQBwd, 1,
1983  "PMC");
1984  CopyBoundaryTrace(Fwd[2], Bwd[2], SolverUtils::eFwdEQNegBwd, 2,
1985  "PMC");
1986 
1987  Array<OneD, NekDouble> e1Fwd_cdot_ncrossdE(nTraceNumPoints);
1988  Array<OneD, NekDouble> e1Bwd_cdot_ncrossdE(nTraceNumPoints);
1989  Array<OneD, NekDouble> e2Fwd_cdot_ncrossdE(nTraceNumPoints);
1990  Array<OneD, NekDouble> e2Bwd_cdot_ncrossdE(nTraceNumPoints);
1991  Array<OneD, NekDouble> e3Fwd_cdot_dHe3(nTraceNumPoints);
1992  Array<OneD, NekDouble> e3Bwd_cdot_dHe3(nTraceNumPoints);
1993 
1994  // Compute numfluxFwd[dir] = (eFwd^[dir] \cdot n \times e^3) * (ZimFwd *
1995  // HFwd^3 + ZimBwd * HBwd^3 )
1996  // Compute numfluxBwd[dir] = (eBwd^[dir] \cdot n \times e^3) * (ZimFwd *
1997  // HFwd^3 + ZimBwd * HBwd^3 )
1998  ComputeNtimesFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0], numfluxFwd[0],
1999  numfluxBwd[0]);
2000  ComputeNtimesFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1], numfluxFwd[1],
2001  numfluxBwd[1]);
2002 
2003  // Compute numfluxFwd[2] = eFwd^3 \cdot ( n1e1 \times ( imFwd EFwd + imBwd
2004  // EBwd ) ) / 2 {{Y_i}}
2005  ComputeNtimesF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
2006  m_YimBwd[1], numfluxFwd[2], numfluxBwd[2]);
2007 
2008  // Compute e1Fwd_cdot_ncrossdE = eFwd[dir] \cdot \alpha n \times n \times
2009  // [E] / 2 {{ZimFwd}}
2010  ComputeNtimestimesdFz(0, Fwd, Bwd, m_ZimFwd[0], m_ZimBwd[0],
2011  e1Fwd_cdot_ncrossdE, e1Bwd_cdot_ncrossdE);
2012  ComputeNtimestimesdFz(1, Fwd, Bwd, m_ZimFwd[1], m_ZimBwd[1],
2013  e2Fwd_cdot_ncrossdE, e2Bwd_cdot_ncrossdE);
2014 
2015  // Compute - \alpha [H3] * ( 1/2{{Yim1}} + 1/2{{Yim2}} )
2016  ComputeNtimestimesdF12(Fwd, Bwd, m_YimFwd[0], m_YimBwd[0], m_YimFwd[1],
2017  m_YimBwd[1], e3Fwd_cdot_dHe3, e3Bwd_cdot_dHe3);
2018 
2019  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[0], 1, e1Fwd_cdot_ncrossdE, 1,
2020  numfluxFwd[0], 1);
2021  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxFwd[1], 1, e2Fwd_cdot_ncrossdE, 1,
2022  numfluxFwd[1], 1);
2023  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxFwd[2], 1, e3Fwd_cdot_dHe3, 1,
2024  numfluxFwd[2], 1);
2025 
2026  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[0], 1, e1Bwd_cdot_ncrossdE, 1,
2027  numfluxBwd[0], 1);
2028  Vmath::Svtvp(nTraceNumPoints, 1.0, numfluxBwd[1], 1, e2Bwd_cdot_ncrossdE, 1,
2029  numfluxBwd[1], 1);
2030  Vmath::Svtvp(nTraceNumPoints, -1.0, numfluxBwd[2], 1, e3Bwd_cdot_dHe3, 1,
2031  numfluxBwd[2], 1);
2032 }
2033 
2034 Array<OneD, NekDouble> MMFSystem::GetIncidentField(const int var,
2035  const NekDouble time)
2036 {
2037  int nq = m_fields[0]->GetNpoints();
2038 
2039  Array<OneD, NekDouble> x0(nq);
2040  Array<OneD, NekDouble> x1(nq);
2041  Array<OneD, NekDouble> x2(nq);
2042 
2043  m_fields[0]->GetCoords(x0, x1, x2);
2044 
2045  // GetSmoothFactor such that wave propages from the left to the object.
2046  // a = 0.1, ta = 1, f = 1.0./(1.0 + exp( -0.5.*(time-ta)/a ));
2047  Array<OneD, NekDouble> SmoothFactor(nq, 1.0);
2048  switch (m_SmoothFactor)
2049  {
2050  case 0:
2051  {
2052  for (int i = 0; i < nq; i++)
2053  {
2054  SmoothFactor[i] = 1.0 / (1.0 + exp(-1.0 * (time - 1.0) / 0.1));
2055  }
2056  }
2057  break;
2058 
2059  case 1:
2060  {
2061  NekDouble xmin = Vmath::Vmin(nq, x0, 1);
2062  NekDouble xp;
2063  for (int i = 0; i < nq; i++)
2064  {
2065  xp = x0[i] - xmin - time - m_SFinit;
2066  if (xp > 0.0)
2067  {
2068  SmoothFactor[i] =
2069  2.0 / (1.0 + exp(0.5 * (sqrt(xp * xp) - 0.1)));
2070  }
2071 
2072  else
2073  {
2074  SmoothFactor[i] = 1.0;
2075  }
2076  }
2077  }
2078  break;
2079 
2080  default:
2081  break;
2082  }
2083 
2084  // Generate a factor for smoothly increasing wave
2085  Array<OneD, NekDouble> F1(nq);
2086  Array<OneD, NekDouble> F2(nq);
2087  Array<OneD, NekDouble> F3(nq);
2088 
2089  Array<OneD, NekDouble> dF1dt(nq);
2090  Array<OneD, NekDouble> dF2dt(nq);
2091  Array<OneD, NekDouble> dF3dt(nq);
2092 
2093  Array<OneD, NekDouble> F1int(nq);
2094  Array<OneD, NekDouble> F2int(nq);
2095  Array<OneD, NekDouble> F3int(nq);
2096 
2097  Array<OneD, NekDouble> outarray(nq);
2098 
2099  switch (m_IncType)
2100  {
2101  case ePlaneWave:
2102  {
2103  NekDouble cs, sn;
2104  NekDouble e1y, e2y;
2105  switch (m_PolType)
2106  {
2107  case eTransMagnetic:
2108  {
2109  // H = 0 \hat{x} - e^{ikx} \hat{y}
2110  // E = e^{ikx} \hat{z}
2111  // outarray1 = Hr1inc
2112  // outarray2 = Hr2inc
2113  // outarray3 = Ezrinc
2114  for (int i = 0; i < nq; i++)
2115  {
2116  e1y = m_movingframes[0][nq + i];
2117  e2y = m_movingframes[1][nq + i];
2118 
2119  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2120  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2121 
2122  F1[i] = -1.0 * cs * e1y;
2123  F2[i] = -1.0 * cs * e2y;
2124  F3[i] = cs;
2125 
2126  dF1dt[i] = -1.0 * m_Incfreq * sn * e1y;
2127  dF2dt[i] = -1.0 * m_Incfreq * sn * e2y;
2128  dF3dt[i] = 1.0 * m_Incfreq * sn;
2129 
2130  F1int[i] = (1.0 / m_Incfreq) * sn * e1y;
2131  F2int[i] = (1.0 / m_Incfreq) * sn * e2y;
2132  F3int[i] = (-1.0 / m_Incfreq) * sn;
2133  }
2134  }
2135  break;
2136 
2137  case eTransElectric:
2138  {
2139  // E = 0 \hat{x} + e^{ikx} \hat{y}
2140  // H = e^{ikx} \hat{z}
2141  // outarray1 = Er1inc
2142  // outarray2 = Er2inc
2143  // outarray3 = Hzrinc
2144  // outarray1 = Ei1inc
2145  // outarray2 = Ei2inc
2146  // outarray3 = Hziinc
2147  for (int i = 0; i < nq; i++)
2148  {
2149  e1y = m_movingframes[0][nq + i];
2150  e2y = m_movingframes[1][nq + i];
2151 
2152  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2153  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2154 
2155  F1[i] = cs * e1y;
2156  F2[i] = cs * e2y;
2157  F3[i] = cs;
2158 
2159  dF1dt[i] = m_Incfreq * sn * e1y;
2160  dF2dt[i] = m_Incfreq * sn * e2y;
2161  dF3dt[i] = m_Incfreq * sn;
2162 
2163  F1int[i] = (-1.0 / m_Incfreq) * sn * e1y;
2164  F2int[i] = (-1.0 / m_Incfreq) * sn * e2y;
2165  F3int[i] = (-1.0 / m_Incfreq) * sn;
2166  }
2167  }
2168  break;
2169 
2170  default:
2171  break;
2172  }
2173  }
2174  break;
2175 
2176  case ePlaneWaveImag:
2177  {
2178  NekDouble cs, sn;
2179  NekDouble e1y, e2y;
2180  switch (m_PolType)
2181  {
2182  case eTransMagnetic:
2183  {
2184  // H = 0 \hat{x} - e^{ikx} \hat{y}
2185  // E = e^{ikx} \hat{z}
2186  // outarray1 = Hr1inc
2187  // outarray2 = Hr2inc
2188  // outarray3 = Ezrinc
2189  for (int i = 0; i < nq; i++)
2190  {
2191  e1y = m_movingframes[0][nq + i];
2192  e2y = m_movingframes[1][nq + i];
2193 
2194  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2195  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2196 
2197  F1[i] = -1.0 * sn * e1y;
2198  F2[i] = -1.0 * sn * e2y;
2199  F3[i] = sn;
2200 
2201  dF1dt[i] = m_Incfreq * cs * e1y;
2202  dF2dt[i] = m_Incfreq * cs * e2y;
2203  dF3dt[i] = -1.0 * m_Incfreq * cs;
2204 
2205  F1int[i] = (-1.0 / m_Incfreq) * cs * e1y;
2206  F2int[i] = (-1.0 / m_Incfreq) * cs * e2y;
2207  F3int[i] = (1.0 / m_Incfreq) * cs;
2208  }
2209  }
2210  break;
2211 
2212  case eTransElectric:
2213  {
2214  // E = 0 \hat{x} + e^{ikx} \hat{y}
2215  // H = e^{ikx} \hat{z}
2216  // outarray1 = Er1inc
2217  // outarray2 = Er2inc
2218  // outarray3 = Hzrinc
2219  // outarray1 = Ei1inc
2220  // outarray2 = Ei2inc
2221  // outarray3 = Hziinc
2222  for (int i = 0; i < nq; i++)
2223  {
2224  e1y = m_movingframes[0][nq + i];
2225  e2y = m_movingframes[1][nq + i];
2226 
2227  cs = SmoothFactor[i] * cos(m_Incfreq * (x0[i] - time));
2228  sn = SmoothFactor[i] * sin(m_Incfreq * (x0[i] - time));
2229 
2230  F1[i] = sn * e1y;
2231  F2[i] = sn * e2y;
2232  F3[i] = sn;
2233 
2234  dF1dt[i] = -1.0 * m_Incfreq * cs * e1y;
2235  dF2dt[i] = -1.0 * m_Incfreq * cs * e2y;
2236  dF3dt[i] = -1.0 * m_Incfreq * cs;
2237 
2238  F1int[i] = (1.0 / m_Incfreq) * cs * e1y;
2239  F2int[i] = (1.0 / m_Incfreq) * cs * e2y;
2240  F3int[i] = (1.0 / m_Incfreq) * cs;
2241  }
2242  }
2243  break;
2244 
2245  default:
2246  break;
2247  }
2248  }
2249  break;
2250 
2251  default:
2252  break;
2253  }
2254 
2255  switch (var)
2256  {
2257  case 0:
2258  {
2259  outarray = F1;
2260  }
2261  break;
2262 
2263  case 1:
2264  {
2265  outarray = F2;
2266  }
2267  break;
2268 
2269  case 2:
2270  {
2271  outarray = F3;
2272  }
2273  break;
2274 
2275  case 10:
2276  {
2277  outarray = dF1dt;
2278  }
2279  break;
2280 
2281  case 11:
2282  {
2283  outarray = dF2dt;
2284  }
2285  break;
2286 
2287  case 12:
2288  {
2289  outarray = dF3dt;
2290  }
2291  break;
2292 
2293  case 20:
2294  {
2295  outarray = F1int;
2296  }
2297  break;
2298 
2299  case 21:
2300  {
2301  outarray = F2int;
2302  }
2303  break;
2304 
2305  case 22:
2306  {
2307  outarray = F3int;
2308  }
2309  break;
2310 
2311  default:
2312  {
2313  Vmath::Zero(nq, outarray, 1);
2314  }
2315  break;
2316  }
2317 
2318  return outarray;
2319 }
2320 
2321 NekDouble MMFSystem::AvgInt(const Array<OneD, const NekDouble> &inarray)
2322 {
2323  int nq = m_fields[0]->GetNpoints();
2324  Array<OneD, NekDouble> Ones(nq, 1.0);
2325 
2326  if (inarray.size() != nq)
2327  {
2328  ASSERTL0(false, "AvgInt Error: Vector size is not correct");
2329  }
2330 
2331  NekDouble jac = m_fields[0]->Integral(Ones);
2332 
2333  return (m_fields[0]->Integral(inarray)) / jac;
2334 }
2335 
2336 NekDouble MMFSystem::AvgAbsInt(const Array<OneD, const NekDouble> &inarray)
2337 {
2338  int nq = m_fields[0]->GetNpoints();
2339  Array<OneD, NekDouble> Ones(nq, 1.0);
2340  Array<OneD, NekDouble> tmp(nq);
2341 
2342  if (inarray.size() != nq)
2343  {
2344  ASSERTL0(false, "AvgAbsInt Error: Vector size is not correct");
2345  }
2346 
2347  NekDouble jac = m_fields[0]->Integral(Ones);
2348 
2349  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2350  return (m_fields[0]->Integral(tmp)) / jac;
2351 }
2352 
2353 NekDouble MMFSystem::AbsIntegral(const Array<OneD, const NekDouble> &inarray)
2354 {
2355  int nq = m_fields[0]->GetNpoints();
2356  Array<OneD, NekDouble> tmp(nq);
2357 
2358  if (inarray.size() != nq)
2359  {
2360  ASSERTL0(false, "AbsIntegral Error: Vector size is not correct");
2361  }
2362 
2363  Vmath::Vabs(nq, inarray, 1, tmp, 1);
2364  return m_fields[0]->Integral(tmp);
2365 }
2366 
2367 NekDouble MMFSystem::RootMeanSquare(const Array<OneD, const NekDouble> &inarray)
2368 {
2369  int Ntot = inarray.size();
2370 
2371  NekDouble reval = 0.0;
2372  for (int i = 0; i < Ntot; ++i)
2373  {
2374  reval = reval + inarray[i] * inarray[i];
2375  }
2376  reval = sqrt(reval / Ntot);
2377 
2378  return reval;
2379 }
2380 
2381 NekDouble MMFSystem::VectorAvgMagnitude(
2382  const Array<OneD, const Array<OneD, NekDouble>> &inarray)
2383 {
2384  int nq = inarray[0].size();
2385 
2386  Array<OneD, NekDouble> tmp(nq, 0.0);
2387  for (int k = 0; k < m_spacedim; k++)
2388  {
2389  Vmath::Vvtvp(nq, &inarray[k][0], 1, &inarray[k][0], 1, &tmp[0], 1,
2390  &tmp[0], 1);
2391  }
2392  Vmath::Vsqrt(nq, tmp, 1, tmp, 1);
2393 
2394  return RootMeanSquare(tmp);
2395 }
2396 
2397 void MMFSystem::BubbleSort(Array<OneD, NekDouble> &refarray,
2398  Array<OneD, NekDouble> &sortarray)
2399 {
2400  int nq = refarray.size();
2401 
2402  bool swapped = true;
2403  int j = 0;
2404  NekDouble tmp;
2405 
2406  while (swapped)
2407  {
2408  swapped = false;
2409  j++;
2410  for (int i = 0; i < nq - j; i++)
2411  {
2412  if (refarray[i] > refarray[i + 1])
2413  {
2414  tmp = refarray[i];
2415  refarray[i] = refarray[i + 1];
2416  refarray[i + 1] = tmp;
2417 
2418  tmp = sortarray[i];
2419  sortarray[i] = sortarray[i + 1];
2420  sortarray[i + 1] = tmp;
2421 
2422  swapped = true;
2423  }
2424  }
2425  }
2426 }
2427 
2428 void MMFSystem::GramSchumitz(
2429  const Array<OneD, const Array<OneD, NekDouble>> &v1,
2430  const Array<OneD, const Array<OneD, NekDouble>> &v2,
2431  Array<OneD, Array<OneD, NekDouble>> &outarray, bool KeepTheMagnitude)
2432 {
2433 
2434  int nq = v1[0].size();
2435  Array<OneD, NekDouble> tmp(nq, 0.0);
2436  Array<OneD, NekDouble> mag(nq, 0.0);
2437 
2438  for (int i = 0; i < m_spacedim; ++i)
2439  {
2440  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2441  Vmath::Vvtvp(nq, &v1[i][0], 1, &v2[i][0], 1, &tmp[0], 1, &tmp[0], 1);
2442  Vmath::Vvtvp(nq, &v1[i][0], 1, &v1[i][0], 1, &mag[0], 1, &mag[0], 1);
2443  }
2444  Vmath::Vdiv(nq, &tmp[0], 1, &mag[0], 1, &tmp[0], 1);
2445  Vmath::Neg(nq, &tmp[0], 1);
2446 
2447  // outarray = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
2448 
2449  // u2 = v2 - < u1 , v2 > ( u1 / < u1, u1 > )
2450  for (int i = 0; i < m_spacedim; ++i)
2451  {
2452  outarray[i] = Array<OneD, NekDouble>(nq, 0.0);
2453  Vmath::Vvtvp(nq, &tmp[0], 1, &v1[i][0], 1, &v2[i][0], 1,
2454  &outarray[i][0], 1);
2455  }
2456 
2457  if (KeepTheMagnitude)
2458  {
2459  Array<OneD, NekDouble> magorig(nq, 0.0);
2460  Array<OneD, NekDouble> magnew(nq, 0.0);
2461 
2462  for (int i = 0; i < m_spacedim; ++i)
2463  {
2464  Vmath::Vmul(nq, &v2[0][0], 1, &v2[0][0], 1, &magorig[0], 1);
2465  Vmath::Vvtvp(nq, &v2[1][0], 1, &v2[1][0], 1, &magorig[0], 1,
2466  &magorig[0], 1);
2467  Vmath::Vvtvp(nq, &v2[2][0], 1, &v2[2][0], 1, &magorig[0], 1,
2468  &magorig[0], 1);
2469 
2470  Vmath::Vmul(nq, &outarray[0][0], 1, &outarray[0][0], 1, &magnew[0],
2471  1);
2472  Vmath::Vvtvp(nq, &outarray[1][0], 1, &outarray[1][0], 1, &magnew[0],
2473  1, &magnew[0], 1);
2474  Vmath::Vvtvp(nq, &outarray[2][0], 1, &outarray[2][0], 1, &magnew[0],
2475  1, &magnew[0], 1);
2476  }
2477 
2478  for (int i = 0; i < m_spacedim; ++i)
2479  {
2480  for (int j = 0; j < nq; ++j)
2481  {
2482  if (fabs(magnew[j]) > 0.000000001)
2483  {
2484  outarray[i][j] =
2485  outarray[i][j] * sqrt(magorig[j] / magnew[j]);
2486  }
2487  }
2488  }
2489  }
2490 }
2491 
2492 void MMFSystem::v_GenerateSummary(SummaryList &s)
2493 {
2494  int nq = m_fields[0]->GetNpoints();
2495  UnsteadySystem::v_GenerateSummary(s);
2496 
2497  AddSummaryItem(s, "Total grids", nq);
2498  AddSummaryItem(s, "Shape Dimension", m_shapedim);
2499  AddSummaryItem(s, "Surface", SurfaceTypeMap[m_surfaceType]);
2500  if (m_surfaceType == eSphere)
2501  {
2502  NekDouble MeshError;
2503 
2504  Array<OneD, NekDouble> x(nq);
2505  Array<OneD, NekDouble> y(nq);
2506  Array<OneD, NekDouble> z(nq);
2507  Array<OneD, NekDouble> rad(nq, 0.0);
2508 
2509  m_fields[0]->GetCoords(x, y, z);
2510 
2511  Vmath::Vvtvp(nq, x, 1, x, 1, rad, 1, rad, 1);
2512  Vmath::Vvtvp(nq, y, 1, y, 1, rad, 1, rad, 1);
2513  Vmath::Vvtvp(nq, z, 1, z, 1, rad, 1, rad, 1);
2514  Vmath::Vsqrt(nq, rad, 1, rad, 1);
2515 
2516  Vmath::Sadd(nq, -1.0, rad, 1, rad, 1);
2517  Vmath::Vabs(nq, rad, 1, rad, 1);
2518 
2519  MeshError = m_fields[0]->Integral(rad);
2520  SolverUtils::AddSummaryItem(s, "Mesh Error", MeshError);
2521  }
2522 
2523  AddSummaryItem(s, "MMFdir", SpatialDomains::GeomMMFMap[m_MMFdir]);
2524  AddSummaryItem(s, "SmoothFactor", m_SmoothFactor);
2525  AddSummaryItem(s, "SFinit", m_SFinit);
2526 
2527  if (fabs(m_alpha - 1.0) > 0.001)
2528  {
2529  AddSummaryItem(s, "Alpha", m_alpha);
2530  }
2531 
2532  if (m_Incfreq > 0.0)
2533  {
2534  AddSummaryItem(s, "Incfreq", m_Incfreq);
2535  }
2536 
2537  if (m_MMFdir == 4)
2538  {
2539  AddSummaryItem(s, "MMFCircAxisX", m_MMFfactors[0]);
2540  AddSummaryItem(s, "MMFCircAxisY", m_MMFfactors[1]);
2541  AddSummaryItem(s, "MMFCircCentreX", m_MMFfactors[2]);
2542  AddSummaryItem(s, "MMFCircCentreY", m_MMFfactors[3]);
2543  }
2544 }
2545 }
2546 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:376
Nektar::SolverUtils::EquationSystem::GetTraceNpoints
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Definition: EquationSystem.h:693
Nektar::SolverUtils::EquationSystem::GetTraceTotPoints
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Definition: EquationSystem.h:688
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:343
Nektar::SolverUtils::EquationSystem::GetExpSize
SOLVER_UTILS_EXPORT int GetExpSize()
Definition: EquationSystem.h:698
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:337
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:400
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:713
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:1678
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:1735
Nektar::SolverUtils::MMFSystem::m_dedxi_cdot_e
Array< OneD, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > > m_dedxi_cdot_e
Definition: MMFSystem.h:214
Nektar::SolverUtils::MMFSystem::m_ZimFwd
Array< OneD, Array< OneD, NekDouble > > m_ZimFwd
Definition: MMFSystem.h:201
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:2428
Nektar::SolverUtils::MMFSystem::BubbleSort
SOLVER_UTILS_EXPORT void BubbleSort(Array< OneD, NekDouble > &refarray, Array< OneD, NekDouble > &sortarray)
Definition: MMFSystem.cpp:2397
Nektar::SolverUtils::MMFSystem::ComputeZimYim
SOLVER_UTILS_EXPORT void ComputeZimYim(Array< OneD, Array< OneD, NekDouble >> &epsvec, Array< OneD, Array< OneD, NekDouble >> &muvec)
Definition: MMFSystem.cpp:1107
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFBwd
Definition: MMFSystem.h:199
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:189
Nektar::SolverUtils::MMFSystem::m_nperpcdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_nperpcdotMFFwd
Definition: MMFSystem.h:186
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:697
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:1614
Nektar::SolverUtils::MMFSystem::AbsIntegral
SOLVER_UTILS_EXPORT NekDouble AbsIntegral(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2353
Nektar::SolverUtils::MMFSystem::m_ntimes_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimes_ntimesMFFwd
Definition: MMFSystem.h:198
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:1013
Nektar::SolverUtils::MMFSystem::AvgInt
SOLVER_UTILS_EXPORT NekDouble AvgInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2321
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:1812
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:187
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:674
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:204
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:207
Nektar::SolverUtils::MMFSystem::m_MMFdir
SpatialDomains::GeomMMF m_MMFdir
Definition: MMFSystem.h:216
Nektar::SolverUtils::MMFSystem::AvgAbsInt
SOLVER_UTILS_EXPORT NekDouble AvgAbsInt(const Array< OneD, const NekDouble > &inarray)
Definition: MMFSystem.cpp:2336
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:180
Nektar::SolverUtils::MMFSystem::m_MMFfactors
Array< OneD, NekDouble > m_MMFfactors
Definition: MMFSystem.h:154
Nektar::SolverUtils::MMFSystem::m_ntimesMFBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFBwd
Definition: MMFSystem.h:197
Nektar::SolverUtils::MMFSystem::m_MFtraceFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceFwd
Definition: MMFSystem.h:193
Nektar::SolverUtils::MMFSystem::VectorAvgMagnitude
SOLVER_UTILS_EXPORT NekDouble VectorAvgMagnitude(const Array< OneD, const Array< OneD, NekDouble >> &inarray)
Definition: MMFSystem.cpp:2381
Nektar::SolverUtils::MMFSystem::m_MFtraceBwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_MFtraceBwd
Definition: MMFSystem.h:194
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:1079
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:2034
Nektar::SolverUtils::MMFSystem::m_ntimesMFFwd
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_ntimesMFFwd
Definition: MMFSystem.h:196
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:1922
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:1571
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:1644
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:1516
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:950
Nektar::SolverUtils::MMFSystem::m_ZimBwd
Array< OneD, Array< OneD, NekDouble > > m_ZimBwd
Definition: MMFSystem.h:202
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:1369
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:184
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:1460
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:914
Nektar::SolverUtils::MMFSystem::m_YimFwd
Array< OneD, Array< OneD, NekDouble > > m_YimFwd
Definition: MMFSystem.h:203
Nektar::SolverUtils::MMFSystem::m_epsvec
Array< OneD, Array< OneD, NekDouble > > m_epsvec
Definition: MMFSystem.h:206
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:2367
Nektar::SolverUtils::MMFSystem::m_ncdotMFFwd
Array< OneD, Array< OneD, NekDouble > > m_ncdotMFFwd
Definition: MMFSystem.h:183
Nektar::SolverUtils::MMFSystem::m_TestMaxwellType
TestMaxwellType m_TestMaxwellType
Definition: MMFSystem.h:150
Nektar::SolverUtils::MMFSystem::m_CurlMF
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_CurlMF
Definition: MMFSystem.h:190
Nektar::SolverUtils::MMFSystem::Computedemdxicdote
SOLVER_UTILS_EXPORT void Computedemdxicdote()
Definition: MMFSystem.cpp:1279
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2492
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)
Print a summary of time stepping parameters.
Definition: UnsteadySystem.cpp:609
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter.cpp:86
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::LibUtilities::EquationSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:131
Nektar::SolverUtils::UpwindTypeMap
const char *const UpwindTypeMap[]
Definition: MMFSystem.h:86
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
Nektar::SolverUtils::SIZE_UpwindType
@ SIZE_UpwindType
Length of enum list.
Definition: MMFSystem.h:83
Nektar::SolverUtils::eAverage
@ eAverage
averaged (or centred) flux
Definition: MMFSystem.h:77
Nektar::SolverUtils::eLaxFriedrich
@ eLaxFriedrich
Lax-Friedrich flux.
Definition: MMFSystem.h:78
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:47
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:174
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
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:475
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:192
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:565
Vmath::Vabs
void Vabs(int n, const T *x, const int incx, T *y, const int incy)
vabs: y = |x|
Definition: Vmath.cpp:493
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
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:992
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:513
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
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:322
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 plus vector): z = w*x - y
Definition: Vmath.cpp:541
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:225
Vmath::Imax
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:866
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:257
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:341
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:892
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
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:372
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Generated on Wed Nov 24 2021 10:33:27 for Nektar++ by doxygen 1.9.1