Nektar++
ContField.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ContField.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: Field definition for a continuous domain with boundary
32 // conditions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ContField.h>
38 #include <tuple>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44 namespace MultiRegions
45 {
46 /**
47  * @class ContField
48  * The class #ContField is
49  * able to incorporate the boundary conditions imposed to the problem
50  * to be solved. Therefore, the class is equipped with three additional
51  * data members:
52  * - #m_bndCondExpansions
53  * - #m_bndTypes
54  * - #m_bndCondEquations
55  *
56  * The first data structure, #m_bndCondExpansions, contains the
57  * one-dimensional spectral/hp expansion on the boundary, #m_bndTypes
58  * stores information about the type of boundary condition on the
59  * different parts of the boundary while #m_bndCondEquations holds the
60  * equation of the imposed boundary conditions.
61  *
62  * Furthermore, in case of Dirichlet boundary conditions, this class is
63  * capable of lifting a known solution satisfying these boundary
64  * conditions. If we denote the unknown solution by
65  * \f$u^{\mathcal{H}}(\boldsymbol{x})\f$ and the known Dirichlet
66  * boundary conditions by \f$u^{\mathcal{D}}(\boldsymbol{x})\f$, the
67  * expansion then can be decomposed as
68  * \f[ u^{\delta}(\boldsymbol{x}_i)=u^{\mathcal{D}}(\boldsymbol{x}_i)+
69  * u^{\mathcal{H}}(\boldsymbol{x}_i)=\sum_{n=0}^{N^{\mathcal{D}}-1}
70  * \hat{u}_n^{\mathcal{D}}\Phi_n(\boldsymbol{x}_i)+
71  * \sum_{n={N^{\mathcal{D}}}}^{N_{\mathrm{dof}}-1}
72  * \hat{u}_n^{\mathcal{H}} \Phi_n(\boldsymbol{x}_i).\f]
73  * This lifting is accomplished by ordering the known global degrees of
74  * freedom, prescribed by the Dirichlet boundary conditions, first in
75  * the global array
76  * \f$\boldsymbol{\hat{u}}\f$, that is,
77  * \f[\boldsymbol{\hat{u}}=\left[ \begin{array}{c}
78  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
79  * \boldsymbol{\hat{u}}^{\mathcal{H}}
80  * \end{array} \right].\f]
81  * Such kind of expansions are also referred to as continuous fields.
82  * This class should be used when solving 2D problems using a standard
83  * Galerkin approach.
84  */
85 
86 /**
87  *
88  */
89 ContField::ContField()
90  : DisContField(), m_locToGloMap(), m_globalMat(),
91  m_globalLinSysManager(
92  std::bind(&ContField::GenGlobalLinSys, this, std::placeholders::_1),
93  std::string("GlobalLinSys"))
94 {
95 }
96 
97 /**
98  * Given a mesh \a graph, containing information about the domain and
99  * the spectral/hp element expansion, this constructor fills the list
100  * of local expansions #m_exp with the proper expansions, calculates
101  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
102  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory
103  * for the arrays #m_coeffs and #m_phys. Furthermore, it constructs the
104  * mapping array (contained in #m_locToGloMap) for the transformation
105  * between local elemental level and global level, it calculates the
106  * total number global expansion coefficients \f$\hat{u}_n\f$ and
107  * allocates memory for the array #m_contCoeffs. The constructor also
108  * discretises the boundary conditions, specified by the argument \a
109  * bcs, by expressing them in terms of the coefficient of the expansion
110  * on the boundary.
111  *
112  * @param graph A mesh, containing information about the domain
113  * and the spectral/hp element expansion.
114  * @param bcs The boundary conditions.
115  * @param variable An optional parameter to indicate for which
116  * variable the field should be constructed.
117  */
120  const std::string &variable,
121  const bool DeclareCoeffPhysArrays,
122  const bool CheckIfSingularSystem,
123  const Collections::ImplementationType ImpType)
124  : DisContField(pSession, graph, variable, false, DeclareCoeffPhysArrays,
125  ImpType),
126  m_globalMat(MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
127  m_globalLinSysManager(
128  std::bind(&ContField::GenGlobalLinSys, this, std::placeholders::_1),
129  std::string("GlobalLinSys"))
130 {
133  CheckIfSingularSystem, variable, m_periodicVerts, m_periodicEdges,
135 
136  if (m_session->DefinesCmdLineArgument("verbose"))
137  {
138  m_locToGloMap->PrintStats(std::cout, variable);
139  }
140 }
141 
142 /**
143  * Given a mesh \a graph, containing information about the domain and
144  * the spectral/hp element expansion, this constructor fills the list
145  * of local expansions #m_exp with the proper expansions, calculates
146  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
147  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory
148  * for the arrays #m_coeffs and #m_phys. Furthermore, it constructs the
149  * mapping array (contained in #m_locToGloMap) for the transformation
150  * between local elemental level and global level, it calculates the
151  * total number global expansion coefficients \f$\hat{u}_n\f$ and
152  * allocates memory for the array #m_coeffs. The constructor also
153  * discretises the boundary conditions, specified by the argument \a
154  * bcs, by expressing them in terms of the coefficient of the expansion
155  * on the boundary.
156  *
157  * @param In Existing ContField object used to provide the
158  * local to global mapping information and
159  * global solution type.
160  * @param graph A mesh, containing information about the domain
161  * and the spectral/hp element expansion.
162  * @param bcs The boundary conditions.
163  * @param bc_loc
164  */
167  const std::string &variable, bool DeclareCoeffPhysArrays,
168  const bool CheckIfSingularSystem)
169  : DisContField(In, graph, variable, false, DeclareCoeffPhysArrays),
170  m_globalMat(MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
171  m_globalLinSysManager(
172  std::bind(&ContField::GenGlobalLinSys, this, std::placeholders::_1),
173  std::string("GlobalLinSys")),
174  m_GJPData(In.m_GJPData)
175 {
176  if (!SameTypeOfBoundaryConditions(In) || CheckIfSingularSystem)
177  {
180  CheckIfSingularSystem, variable, m_periodicVerts, m_periodicEdges,
182 
183  if (m_session->DefinesCmdLineArgument("verbose"))
184  {
185  m_locToGloMap->PrintStats(std::cout, variable);
186  }
187  }
188  else
189  {
191  }
192 }
193 
194 /**
195  * Initialises the object as a copy of an existing ContField object.
196  * @param In Existing ContField object.
197  * @param DeclareCoeffPhysArrays bool to declare if \a m_phys
198  * and \a m_coeffs should be declared. Default is true
199  */
200 ContField::ContField(const ContField &In, bool DeclareCoeffPhysArrays)
201  : DisContField(In, DeclareCoeffPhysArrays), m_locToGloMap(In.m_locToGloMap),
202  m_globalMat(In.m_globalMat),
203  m_globalLinSysManager(In.m_globalLinSysManager), m_GJPData(In.m_GJPData)
204 {
205 }
206 
207 /**
208  * Constructs a continuous field as a copy of an existing
209  * explist field and adding all the boundary conditions.
210  *
211  * @param In Existing explist1D field .
212  */
214  const ExpList &In)
215  : DisContField(In), m_locToGloMap(),
216  m_globalLinSysManager(
217  std::bind(&ContField::GenGlobalLinSys, this, std::placeholders::_1),
218  std::string("GlobalLinSys"))
219 {
221  pSession, m_ncoeffs, In);
222 }
223 
224 /**
225  *
226  */
228 {
229 }
230 
231 /**
232  * Given a function \f$f(\boldsymbol{x})\f$ defined at the quadrature
233  * points, this function determines the unknown global coefficients
234  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ employing a discrete
235  * Galerkin projection from physical space to coefficient
236  * space. The operation is evaluated by the function #GlobalSolve using
237  * the global mass matrix.
238  *
239  * The values of the function \f$f(\boldsymbol{x})\f$ evaluated at the
240  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
241  * variable #m_phys of the ExpList object \a Sin. The resulting global
242  * coefficients \f$\hat{u}_g\f$ are stored in the array #m_coeffs.
243  *
244  * @param Sin An ExpList, containing the discrete evaluation
245  * of \f$f(\boldsymbol{x})\f$ at the quadrature
246  * points in its array #m_phys.
247  */
249  Array<OneD, NekDouble> &outarray)
250 
251 {
252  // Inner product of forcing
254  IProductWRTBase(inarray, wsp);
255 
256  // Solve the system
258 
259  GlobalSolve(key, wsp, outarray);
260 }
261 
262 /**
263  *
264  */
266 {
267  int Ncoeffs = m_locToGloMap->GetNumLocalCoeffs();
268  Array<OneD, NekDouble> tmp1(Ncoeffs);
269  Array<OneD, NekDouble> tmp2(Ncoeffs);
270 
271  IProductWRTBase(field, tmp1);
272  MultiplyByInvMassMatrix(tmp1, tmp2);
273  BwdTrans(tmp2, field);
274 }
275 
276 /**
277  * Computes the matrix vector product
278  * @f$ \mathbf{y} = \mathbf{M}^{-1}\mathbf{x} @f$.
279  *
280  * @param inarray Input vector @f$\mathbf{x}@f$.
281  * @param outarray Output vector @f$\mathbf{y}@f$.
282  */
284  const Array<OneD, const NekDouble> &inarray,
285  Array<OneD, NekDouble> &outarray)
286 
287 {
288 
290  GlobalSolve(key, inarray, outarray);
291 }
292 
293 /**
294  * Consider the two dimensional Laplace equation,
295  * \f[\nabla\cdot\left(\boldsymbol{\sigma}\nabla
296  * u(\boldsymbol{x})\right) = f(\boldsymbol{x}),\f] supplemented with
297  * appropriate boundary conditions (which are contained in the data
298  * member #m_bndCondExpansions). In the equation above
299  * \f$\boldsymbol{\sigma}\f$ is the (symmetric positive definite)
300  * diffusion tensor:
301  * \f[ \sigma = \left[ \begin{array}{cc}
302  * \sigma_{00}(\boldsymbol{x},t) & \sigma_{01}(\boldsymbol{x},t) \\
303  * \sigma_{01}(\boldsymbol{x},t) & \sigma_{11}(\boldsymbol{x},t)
304  * \end{array} \right]. \f]
305  * Applying a \f$C^0\f$ continuous Galerkin discretisation, this
306  * equation leads to the following linear system:
307  * \f[\boldsymbol{L}
308  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f]
309  * where \f$\boldsymbol{L}\f$ is the Laplacian matrix. This function
310  * solves the system above for the global coefficients
311  * \f$\boldsymbol{\hat{u}}\f$ by a call to the function #GlobalSolve.
312  *
313  * The values of the function \f$f(\boldsymbol{x})\f$ evaluated at the
314  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
315  * variable #m_phys of the ExpList object \a Sin. The resulting global
316  * coefficients \f$\boldsymbol{\hat{u}}_g\f$ are stored in the array
317  * #m_coeffs.
318  *
319  * @param Sin An ExpList, containing the discrete evaluation
320  * of the forcing function \f$f(\boldsymbol{x})\f$
321  * at the quadrature points in its array #m_phys.
322  * @param variablecoeffs The (optional) parameter containing the
323  * coefficients evaluated at the quadrature
324  * points. It is an Array of (three) arrays which
325  * stores the laplacian coefficients in the
326  * following way
327  * \f[\mathrm{variablecoeffs} = \left[ \begin{array}{c}
328  * \left[\sigma_{00}(\boldsymbol{x_i},t)\right]_i \\
329  * \left[\sigma_{01}(\boldsymbol{x_i},t)\right]_i \\
330  * \left[\sigma_{11}(\boldsymbol{x_i},t)\right]_i
331  * \end{array}\right]
332  * \f]
333  * If this argument is not passed to the function, the following
334  * equation will be solved:
335  * \f[\nabla^2u(\boldsymbol{x}) = f(\boldsymbol{x}),\f]
336  *
337  * @param time The time-level at which the coefficients are
338  * evaluated
339  */
341  const Array<OneD, const NekDouble> &inarray,
342  Array<OneD, NekDouble> &outarray,
343  const Array<OneD, const NekDouble> &dirForcing,
344  const Array<OneD, Array<OneD, NekDouble>> &variablecoeffs, NekDouble time)
345 {
346  // Inner product of forcing
348  IProductWRTBase(inarray, wsp);
349 
350  // Note -1.0 term necessary to invert forcing function to
351  // be consistent with matrix definition
352  Vmath::Neg(m_ncoeffs, wsp, 1);
353 
354  // Forcing function with weak boundary conditions
355  int i, j;
356  int bndcnt = 0;
358  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
359  const Array<OneD, const int> map =
360  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
361 
362  // Add weak boundary conditions to forcing
363  for (i = 0; i < m_bndCondExpansions.size(); ++i)
364  {
365  if (m_bndConditions[i]->GetBoundaryConditionType() ==
367  m_bndConditions[i]->GetBoundaryConditionType() ==
369  {
370  const Array<OneD, NekDouble> bndcoeff =
372 
373  if (m_locToGloMap->GetSignChange())
374  {
375  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
376  {
377  wsp[map[bndcnt + j]] += sign[bndcnt + j] * bndcoeff[j];
378  }
379  }
380  else
381  {
382  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
383  {
384  wsp[map[bndcnt + j]] += bndcoeff[bndcnt + j];
385  }
386  }
387  }
388 
389  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
390  }
391 
392  StdRegions::VarCoeffMap varcoeffs;
393  varcoeffs[StdRegions::eVarCoeffD00] = variablecoeffs[0];
394  varcoeffs[StdRegions::eVarCoeffD01] = variablecoeffs[1];
395  varcoeffs[StdRegions::eVarCoeffD11] = variablecoeffs[3];
396  varcoeffs[StdRegions::eVarCoeffD22] = variablecoeffs[5];
398  factors[StdRegions::eFactorTime] = time;
399 
400  // Solve the system
402  varcoeffs);
403 
404  GlobalSolve(key, wsp, outarray, dirForcing);
405 }
406 
407 /**
408  * Constructs the GlobalLinearSysKey for the linear advection operator
409  * with the supplied parameters, and computes the eigenvectors and
410  * eigenvalues of the associated matrix.
411  * @param ax Advection parameter, x.
412  * @param ay Advection parameter, y.
413  * @param Real Computed eigenvalues, real component.
414  * @param Imag Computed eigenvalues, imag component.
415  * @param Evecs Computed eigenvectors.
416  */
420  Array<OneD, NekDouble> &Evecs)
421 {
422  // Solve the system
426  vel[0] = vel_x;
427  vel[1] = vel_y;
428 
429  StdRegions::VarCoeffMap varcoeffs;
430  varcoeffs[StdRegions::eVarCoeffVelX] =
432  varcoeffs[StdRegions::eVarCoeffVelY] =
435  factors[StdRegions::eFactorTime] = 0.0;
437  factors, varcoeffs);
438 
440  Gmat->EigenSolve(Real, Imag, Evecs);
441 }
442 
443 /**
444  * Given a linear system specified by the key \a key,
445  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}},\f]
446  * this function solves this linear system taking into account the
447  * boundary conditions specified in the data member
448  * #m_bndCondExpansions. Therefore, it adds an array
449  * \f$\boldsymbol{\hat{g}}\f$ which represents the non-zero surface
450  * integral resulting from the weak boundary conditions (e.g. Neumann
451  * boundary conditions) to the right hand side, that is,
452  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}+
453  * \boldsymbol{\hat{g}}.\f]
454  * Furthermore, it lifts the known degrees of freedom which are
455  * prescribed by the Dirichlet boundary conditions. As these known
456  * coefficients \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ are numbered
457  * first in the global coefficient array \f$\boldsymbol{\hat{u}}_g\f$,
458  * the linear system can be decomposed as,
459  * \f[\left[\begin{array}{cc}
460  * \boldsymbol{M}^{\mathcal{DD}}&\boldsymbol{M}^{\mathcal{DH}}\\
461  * \boldsymbol{M}^{\mathcal{HD}}&\boldsymbol{M}^{\mathcal{HH}}
462  * \end{array}\right]
463  * \left[\begin{array}{c}
464  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
465  * \boldsymbol{\hat{u}}^{\mathcal{H}}
466  * \end{array}\right]=
467  * \left[\begin{array}{c}
468  * \boldsymbol{\hat{f}}^{\mathcal{D}}\\
469  * \boldsymbol{\hat{f}}^{\mathcal{H}}
470  * \end{array}\right]+
471  * \left[\begin{array}{c}
472  * \boldsymbol{\hat{g}}^{\mathcal{D}}\\
473  * \boldsymbol{\hat{g}}^{\mathcal{H}}
474  * \end{array}\right]
475  * \f]
476  * which will then be solved for the unknown coefficients
477  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ as,
478  * \f[
479  * \boldsymbol{M}^{\mathcal{HH}}\boldsymbol{\hat{u}}^{\mathcal{H}}=
480  * \boldsymbol{\hat{f}}^{\mathcal{H}}+
481  * \boldsymbol{\hat{g}}^{\mathcal{H}}-
482  * \boldsymbol{M}^{\mathcal{HD}}\boldsymbol{\hat{u}}^{\mathcal{D}}\f]
483  *
484  * @param mkey This key uniquely defines the linear system to
485  * be solved.
486  * @param locrhs contains the forcing term in local coefficient space
487  * @note inout contains initial guess and final output in local coeffs.
488  */
490  const Array<OneD, const NekDouble> &locrhs,
491  Array<OneD, NekDouble> &inout,
492  const Array<OneD, const NekDouble> &dirForcing)
493 {
494  int NumDirBcs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
495  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
496 
497  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE
498  // IN THE SOLUTION ARRAY
500 
501  // STEP 2: CALCULATE THE HOMOGENEOUS COEFFICIENTS
502  if (contNcoeffs - NumDirBcs > 0)
503  {
505  LinSys->Solve(locrhs, inout, m_locToGloMap, dirForcing);
506  }
507 }
508 
509 /**
510  * Returns the global matrix associated with the given GlobalMatrixKey.
511  * If the global matrix has not yet been constructed on this field,
512  * it is first constructed using GenGlobalMatrix().
513  * @param mkey Global matrix key.
514  * @returns Assocated global matrix.
515  */
517 {
519  "To use method must have a AssemblyMap "
520  "attached to key");
521 
522  GlobalMatrixSharedPtr glo_matrix;
523  auto matrixIter = m_globalMat->find(mkey);
524 
525  if (matrixIter == m_globalMat->end())
526  {
527  glo_matrix = GenGlobalMatrix(mkey, m_locToGloMap);
528  (*m_globalMat)[mkey] = glo_matrix;
529  }
530  else
531  {
532  glo_matrix = matrixIter->second;
533  }
534 
535  return glo_matrix;
536 }
537 
538 /**
539  * The function searches the map #m_globalLinSys to see if the
540  * global matrix has been created before. If not, it calls the function
541  * #GenGlobalLinSys to generate the requested global system.
542  *
543  * @param mkey This key uniquely defines the requested
544  * linear system.
545  */
547 {
548  return m_globalLinSysManager[mkey];
549 }
550 
552 {
554  "To use method must have a AssemblyMap "
555  "attached to key");
557 }
558 
559 /**
560  *
561  */
563  Array<OneD, NekDouble> &outarray)
564 {
565  FwdTrans(inarray, outarray);
566 }
567 
569 {
570  int i, j;
571  int bndcnt = 0;
572 
574  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
575  const Array<OneD, const int> map =
576  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
577 
578  for (i = 0; i < m_bndCondExpansions.size(); ++i)
579  {
580  if (m_bndConditions[i]->GetBoundaryConditionType() ==
582  {
583  const Array<OneD, NekDouble> bndcoeff =
585 
586  if (m_locToGloMap->GetSignChange())
587  {
588  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
589  {
590  outarray[map[bndcnt + j]] = sign[bndcnt + j] * bndcoeff[j];
591  }
592  }
593  else
594  {
595  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
596  {
597  outarray[map[bndcnt + j]] = bndcoeff[j];
598  }
599  }
600  }
601  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
602  }
603 
604  // communicate local Dirichlet coeffs that are just
605  // touching a dirichlet boundary on another partition
606  set<int> &ParallelDirBndSign = m_locToGloMap->GetParallelDirBndSign();
607 
608  for (auto &it : ParallelDirBndSign)
609  {
610  outarray[it] *= -1;
611  }
612 
613  m_locToGloMap->UniversalAbsMaxBnd(outarray);
614 
615  for (auto &it : ParallelDirBndSign)
616  {
617  outarray[it] *= -1;
618  }
619 
620  set<ExtraDirDof> &copyLocalDirDofs = m_locToGloMap->GetCopyLocalDirDofs();
621  for (auto &it : copyLocalDirDofs)
622  {
623  outarray[std::get<0>(it)] = outarray[std::get<1>(it)] * std::get<2>(it);
624  }
625 }
626 
628 {
629  int bndcnt = 0;
630 
632  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
633  const Array<OneD, const int> bndmap =
634  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
635 
636  for (int i = 0; i < m_bndCondExpansions.size(); ++i)
637  {
638  Array<OneD, NekDouble> &coeffs = m_bndCondExpansions[i]->UpdateCoeffs();
639 
640  if (m_locToGloMap->GetSignChange())
641  {
642  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
643  {
644  coeffs[j] = sign[bndcnt + j] * m_coeffs[bndmap[bndcnt + j]];
645  }
646  }
647  else
648  {
649  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
650  {
651  coeffs[j] = m_coeffs[bndmap[bndcnt + j]];
652  }
653  }
654 
655  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
656  }
657 }
658 
660 {
661  int bndcnt = 0;
662 
663  ASSERTL1(nreg < m_bndCondExpansions.size(),
664  "nreg is out or range since this many boundary "
665  "regions to not exist");
666 
668  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
669  const Array<OneD, const int> bndmap =
670  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
671 
672  // Now fill in all other Dirichlet coefficients.
673  Array<OneD, NekDouble> &coeffs = m_bndCondExpansions[nreg]->UpdateCoeffs();
674 
675  for (int j = 0; j < nreg; ++j)
676  {
677  bndcnt += m_bndCondExpansions[j]->GetNcoeffs();
678  }
679 
680  if (m_locToGloMap->GetSignChange())
681  {
682  for (int j = 0; j < (m_bndCondExpansions[nreg])->GetNcoeffs(); ++j)
683  {
684  coeffs[j] = sign[bndcnt + j] * m_coeffs[bndmap[bndcnt + j]];
685  }
686  }
687  else
688  {
689  for (int j = 0; j < (m_bndCondExpansions[nreg])->GetNcoeffs(); ++j)
690  {
691  coeffs[j] = m_coeffs[bndmap[bndcnt + j]];
692  }
693  }
694 }
695 
696 /**
697  * This operation is evaluated as:
698  * \f{tabbing}
699  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
700  * > > Do \= $i=$ $0,N_m^e-1$ \\
701  * > > > $\boldsymbol{\hat{u}}^{e}[i] = \mbox{sign}[e][i] \cdot
702  * \boldsymbol{\hat{u}}_g[\mbox{map}[e][i]]$ \\
703  * > > continue \\
704  * > continue
705  * \f}
706  * where \a map\f$[e][i]\f$ is the mapping array and \a
707  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
708  * correct modal connectivity between the different elements (both
709  * these arrays are contained in the data member #m_locToGloMap). This
710  * operation is equivalent to the scatter operation
711  * \f$\boldsymbol{\hat{u}}_l=\mathcal{A}\boldsymbol{\hat{u}}_g\f$,
712  * where \f$\mathcal{A}\f$ is the
713  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
714  *
715  * @note The array #m_coeffs should be filled with the global
716  * coefficients \f$\boldsymbol{\hat{u}}_g\f$ and that the resulting
717  * local coefficients \f$\boldsymbol{\hat{u}}_l\f$ will be stored in
718  * #m_coeffs.
719  */
721  Array<OneD, NekDouble> &outarray)
722 {
723  m_locToGloMap->GlobalToLocal(inarray, outarray);
724 }
725 
727 {
728  m_locToGloMap->GlobalToLocal(m_coeffs, m_coeffs);
729 }
730 
731 /**
732  * This operation is evaluated as:
733  * \f{tabbing}
734  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
735  * > > Do \= $i=$ $0,N_m^e-1$ \\
736  * > > > $\boldsymbol{\hat{u}}_g[\mbox{map}[e][i]] =
737  * \mbox{sign}[e][i] \cdot \boldsymbol{\hat{u}}^{e}[i]$\\
738  * > > continue\\
739  * > continue
740  * \f}
741  * where \a map\f$[e][i]\f$ is the mapping array and \a
742  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
743  * correct modal connectivity between the different elements (both
744  * these arrays are contained in the data member #m_locToGloMap). This
745  * operation is equivalent to the gather operation
746  * \f$\boldsymbol{\hat{u}}_g=\mathcal{A}^{-1}\boldsymbol{\hat{u}}_l\f$,
747  * where \f$\mathcal{A}\f$ is the
748  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
749  *
750  * @note The array #m_coeffs should be filled with the local
751  * coefficients \f$\boldsymbol{\hat{u}}_l\f$ and that
752  * the resulting global coefficients
753  * \f$\boldsymbol{\hat{u}}_g\f$ will be stored in
754  * #m_coeffs. Also if useComm is set to false then no
755  * communication call will be made to check if all
756  * values are consistent over processors
757  */
758 
760  Array<OneD, NekDouble> &outarray, bool useComm)
761 {
762  m_locToGloMap->LocalToGlobal(inarray, outarray, useComm);
763 }
764 
765 void ContField::v_LocalToGlobal(bool useComm)
766 
767 {
768  m_locToGloMap->LocalToGlobal(m_coeffs, m_coeffs, useComm);
769 }
770 
771 /**
772  *
773  */
775  const Array<OneD, const NekDouble> &inarray,
776  Array<OneD, NekDouble> &outarray)
777 {
778  MultiplyByInvMassMatrix(inarray, outarray);
779 }
780 
781 /**
782  * Consider the two dimensional Helmholtz equation,
783  * \f[\nabla^2u(\boldsymbol{x})-\lambda u(\boldsymbol{x})
784  * = f(\boldsymbol{x}),\f] supplemented with appropriate boundary
785  * conditions (which are contained in the data member
786  * #m_bndCondExpansions). Applying a \f$C^0\f$ continuous Galerkin
787  * discretisation, this equation leads to the following linear system:
788  * \f[\left(\boldsymbol{L}+\lambda\boldsymbol{M}\right)
789  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f] where
790  * \f$\boldsymbol{L}\f$ and \f$\boldsymbol{M}\f$ are the Laplacian and
791  * mass matrix respectively. This function solves the system above for
792  * the global coefficients \f$\boldsymbol{\hat{u}}\f$ by a call to the
793  * function #GlobalSolve. It is assumed #m_coeff contains an
794  * initial estimate for the solution.
795  *
796  * The values of the function \f$f(\boldsymbol{x})\f$
797  * evaluated at the quadrature points \f$\boldsymbol{x}_i\f$
798  * should be contained in the variable #m_phys of the ExpList
799  * object \a inarray.
800  *
801  * @param inarray An ExpList, containing the discrete evaluation
802  * of the forcing function \f$f(\boldsymbol{x})\f$
803  * at the quadrature points in its array #m_phys.
804  * @param factors The parameter \f$\lambda\f$ of the Helmholtz
805  * equation is specified through the factors map
806  */
808  Array<OneD, NekDouble> &outarray,
809  const StdRegions::ConstFactorMap &factors,
810  const StdRegions::VarCoeffMap &pvarcoeff,
811  const MultiRegions::VarFactorsMap &varfactors,
812  const Array<OneD, const NekDouble> &dirForcing,
813  const bool PhysSpaceForcing)
814 
815 {
816  int i, j;
817 
818  //----------------------------------
819  // Setup RHS Inner product
820  //----------------------------------
821  // Inner product of forcing
823  if (PhysSpaceForcing)
824  {
825  IProductWRTBase(inarray, wsp);
826  // Note -1.0 term necessary to invert forcing function to
827  // be consistent with matrix definition
828  Vmath::Neg(m_ncoeffs, wsp, 1);
829  }
830  else
831  {
832  Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, wsp, 1);
833  }
834 
835  int bndcnt = 0;
837  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
838  const Array<OneD, const int> map =
839  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
840  // Add weak boundary conditions to forcing
841  for (i = 0; i < m_bndCondExpansions.size(); ++i)
842  {
843  if (m_bndConditions[i]->GetBoundaryConditionType() ==
845  m_bndConditions[i]->GetBoundaryConditionType() ==
847  {
848  const Array<OneD, NekDouble> bndcoeff =
850 
851  if (m_locToGloMap->GetSignChange())
852  {
853  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
854  {
855  wsp[map[bndcnt + j]] += sign[bndcnt + j] * bndcoeff[j];
856  }
857  }
858  else
859  {
860  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
861  {
862  wsp[map[bndcnt + j]] += bndcoeff[j];
863  }
864  }
865  }
866  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
867  }
868 
870 
871  StdRegions::VarCoeffMap varcoeff(pvarcoeff);
872  if (factors.count(StdRegions::eFactorGJP))
873  {
874  // initialize if required
875  if (!m_GJPData)
876  {
878  GetSharedThisPtr());
879  }
880 
881  if (m_GJPData->IsSemiImplicit())
882  {
884  }
885 
886  // to set up forcing need initial guess in physical space
888  BwdTrans(outarray, phys);
889  NekDouble scale = -1.0 * factors.find(StdRegions::eFactorGJP)->second;
890 
891  m_GJPData->Apply(
892  phys, wsp,
893  pvarcoeff.count(StdRegions::eVarCoeffGJPNormVel)
894  ? pvarcoeff.find(StdRegions::eVarCoeffGJPNormVel)->second
896  scale);
897 
898  varcoeff.erase(StdRegions::eVarCoeffGJPNormVel);
899  }
900 
901  GlobalLinSysKey key(mtype, m_locToGloMap, factors, varcoeff, varfactors);
902 
903  GlobalSolve(key, wsp, outarray, dirForcing);
904 }
905 
906 /**
907  * First compute the inner product of forcing function with respect to
908  * base, and then solve the system with the linear advection operator.
909  * @param velocity Array of advection velocities in physical space
910  * @param inarray Forcing function.
911  * @param outarray Result.
912  * @param lambda reaction coefficient
913  * @param dirForcing Dirichlet Forcing.
914  */
915 
916 // could combine this with HelmholtzCG.
918  const Array<OneD, Array<OneD, NekDouble>> &velocity,
919  const Array<OneD, const NekDouble> &inarray,
920  Array<OneD, NekDouble> &outarray, const NekDouble lambda,
921  const Array<OneD, const NekDouble> &dirForcing)
922 {
923  // Inner product of forcing
925  IProductWRTBase(inarray, wsp);
926 
927  // Note -1.0 term necessary to invert forcing function to
928  // be consistent with matrix definition
929  Vmath::Neg(m_ncoeffs, wsp, 1);
930 
931  // Forcing function with weak boundary conditions
932  int i, j;
933  int bndcnt = 0;
935  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsSign();
936  const Array<OneD, const int> map =
937  m_locToGloMap->GetBndCondCoeffsToLocalCoeffsMap();
938  // Add weak boundary conditions to forcing
939  for (i = 0; i < m_bndCondExpansions.size(); ++i)
940  {
941  if (m_bndConditions[i]->GetBoundaryConditionType() ==
943  m_bndConditions[i]->GetBoundaryConditionType() ==
945  {
946  const Array<OneD, NekDouble> bndcoeff =
948 
949  if (m_locToGloMap->GetSignChange())
950  {
951  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
952  {
953  wsp[map[bndcnt + j]] += sign[bndcnt + j] * bndcoeff[j];
954  }
955  }
956  else
957  {
958  for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
959  {
960  wsp[map[bndcnt + j]] += bndcoeff[bndcnt + j];
961  }
962  }
963  }
964 
965  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
966  }
967 
968  // Solve the system
970  factors[StdRegions::eFactorLambda] = lambda;
971  StdRegions::VarCoeffMap varcoeffs;
972  varcoeffs[StdRegions::eVarCoeffVelX] = velocity[0];
973  varcoeffs[StdRegions::eVarCoeffVelY] = velocity[1];
974  if (m_expType == e3D)
975  {
976  varcoeffs[StdRegions::eVarCoeffVelZ] = velocity[2];
977  }
978 
980  m_locToGloMap, factors, varcoeffs);
981 
982  GlobalSolve(key, wsp, outarray, dirForcing);
983 }
984 
985 /**
986  * First compute the inner product of forcing function with respect to
987  * base, and then solve the system with the linear advection operator.
988  * @param velocity Array of advection velocities in physical space
989  * @param inarray Forcing function.
990  * @param outarray Result.
991  * @param lambda reaction coefficient
992  * @param dirForcing Dirichlet Forcing.
993  */
995  const Array<OneD, Array<OneD, NekDouble>> &velocity,
996  const Array<OneD, const NekDouble> &inarray,
997  Array<OneD, NekDouble> &outarray, const NekDouble lambda,
998  const Array<OneD, const NekDouble> &dirForcing)
999 {
1000  // Inner product of forcing
1002  IProductWRTBase(inarray, wsp);
1003 
1004  // Solve the system
1006  factors[StdRegions::eFactorLambda] = lambda;
1007  StdRegions::VarCoeffMap varcoeffs;
1008  varcoeffs[StdRegions::eVarCoeffVelX] = velocity[0];
1009  varcoeffs[StdRegions::eVarCoeffVelY] = velocity[1];
1011  factors, varcoeffs);
1012 
1013  GlobalSolve(key, wsp, outarray, dirForcing);
1014 }
1015 
1016 /**
1017  *
1018  */
1021 {
1022  return GetBndConditions();
1023 }
1024 
1025 /**
1026  * Reset the GlobalLinSys Manager
1027  */
1029 {
1030  m_globalLinSysManager.ClearManager("GlobalLinSys");
1031 }
1032 
1033 } // namespace MultiRegions
1034 } // namespace Nektar
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
This class is the abstraction of a global continuous two- dimensional spectral/hp element expansion w...
Definition: ContField.h:57
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Template method virtual forwarder for MultiplyByInvMassMatrix().
Definition: ContField.cpp:774
virtual ~ContField()
The default destructor.
Definition: ContField.cpp:227
LibUtilities::NekManager< GlobalLinSysKey, GlobalLinSys > m_globalLinSysManager
A manager which collects all the global linear systems being assembled, such that they should be cons...
Definition: ContField.h:183
void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ContField.cpp:994
GlobalLinSysSharedPtr GenGlobalLinSys(const GlobalLinSysKey &mkey)
Definition: ContField.cpp:551
AssemblyMapCGSharedPtr m_locToGloMap
(A shared pointer to) the object which contains all the required information for the transformation f...
Definition: ContField.h:172
virtual void v_LinearAdvectionDiffusionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ContField.cpp:917
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
Template method virtual forwarded for SmoothField().
Definition: ContField.cpp:265
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool useComm)
Gathers the global coefficients from the local coefficients .
Definition: ContField.cpp:759
const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & GetBndConditions()
Returns the boundary conditions.
Definition: ContField.h:354
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Template method virtual forwarder for FwdTrans().
Definition: ContField.cpp:562
void MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Multiply a solution by the inverse mass matrix.
Definition: ContField.cpp:283
void LaplaceSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray, const Array< OneD, Array< OneD, NekDouble >> &variablecoeffs=NullNekDoubleArrayOfArray, NekDouble time=0.0)
Solves the two-dimensional Laplace equation, subject to the boundary conditions specified.
Definition: ContField.cpp:340
virtual void v_GlobalToLocal(void)
Definition: ContField.cpp:726
GlobalMatrixMapShPtr m_globalMat
(A shared pointer to) a list which collects all the global matrices being assembled,...
Definition: ContField.h:177
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Solves the two-dimensional Helmholtz equation, subject to the boundary conditions specified.
Definition: ContField.cpp:807
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
Impose the Dirichlet Boundary Conditions on outarray.
Definition: ContField.cpp:568
virtual void v_ClearGlobalLinSysManager(void)
Definition: ContField.cpp:1028
GlobalLinSysSharedPtr GetGlobalLinSys(const GlobalLinSysKey &mkey)
Returns the linear system specified by the key mkey.
Definition: ContField.cpp:546
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Performs the global forward transformation of a function , subject to the boundary conditions specifi...
Definition: ContField.cpp:248
GlobalMatrixSharedPtr GetGlobalMatrix(const GlobalMatrixKey &mkey)
Returns the global matrix specified by mkey.
Definition: ContField.cpp:516
virtual void v_FillBndCondFromField()
Definition: ContField.cpp:627
ContField()
The default constructor.
Definition: ContField.cpp:89
virtual const Array< OneD, const SpatialDomains ::BoundaryConditionShPtr > & v_GetBndConditions()
Template method virtual forwarder for GetBndConditions().
Definition: ContField.cpp:1020
void LinearAdvectionEigs(const NekDouble ax, const NekDouble ay, Array< OneD, NekDouble > &Real, Array< OneD, NekDouble > &Imag, Array< OneD, NekDouble > &Evecs=NullNekDouble1DArray)
Compute the eigenvalues of the linear advection operator.
Definition: ContField.cpp:417
void GlobalSolve(const GlobalLinSysKey &key, const Array< OneD, const NekDouble > &rhs, Array< OneD, NekDouble > &inout, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Solves the linear system specified by the key key.
Definition: ContField.cpp:489
GJPStabilisationSharedPtr m_GJPData
Data for Gradient Jump Penalisation (GJP) stabilisaiton.
Definition: ContField.h:186
This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansi...
Definition: DisContField.h:56
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
Definition: DisContField.h:173
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
Definition: DisContField.h:178
bool SameTypeOfBoundaryConditions(const DisContField &In)
Check to see if expansion has the same BCs as In.
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
Definition: DisContField.h:149
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
Definition: DisContField.h:168
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:143
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1158
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1787
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1641
std::shared_ptr< GlobalMatrix > GenGlobalMatrix(const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Generates a global matrix from the given key and map.
Definition: ExpList.cpp:2126
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:2100
std::shared_ptr< GlobalLinSys > GenGlobalLinSys(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
This operation constructs the global linear system of type mkey.
Definition: ExpList.cpp:2413
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1136
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Definition: ExpList.cpp:2260
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: ExpList.h:1034
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
Definition: ExpList.h:1849
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1119
Describes a matrix with ordering defined by a local to global map.
bool LocToGloMapIsDefined() const
Returns true if a local to global map is defined.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
Definition: GlobalMatrix.h:86
std::map< GlobalMatrixKey, GlobalMatrixSharedPtr > GlobalMatrixMap
Mapping from global matrix keys to global matrices.
Definition: GlobalMatrix.h:88
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248