Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ContField1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ContField1D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Field definition for 1D domain with boundary conditions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/ContField1D.h>
37 #include <MultiRegions/AssemblyMap/AssemblyMapCG1D.h>
38 
39 namespace Nektar
40 {
41  namespace MultiRegions
42  {
43  /**
44  * @class ContField1D
45  * As opposed to the class #ContExpList1D, the class #ContField1D is
46  * able to incorporate the boundary conditions imposed to the problem
47  * to be solved. Therefore, the class is equipped with three additional
48  * data members:
49  * - #m_bndCondExpansions
50  * - #m_bndTypes
51  * - #m_bndCondEquations
52  *
53  * The first data structure, #m_bndCondExpansions,
54  * contains the point Expansion on the boundary, #m_bndTypes
55  * stores information about the type of boundary condition on the
56  * different parts of the boundary while #m_bndCondEquations holds the
57  * equation of the imposed boundary conditions.
58  *
59  * Furthermore, in case of Dirichlet boundary conditions,
60  * this class is capable of lifting a known solution satisfying these
61  * boundary conditions. If we denote the unknown solution by
62  * \f$u^{\mathcal{H}}(\boldsymbol{x})\f$ and the known Dirichlet
63  * boundary conditions by \f$u^{\mathcal{D}}(\boldsymbol{x})\f$, the
64  * expansion then can be decomposed as
65  * \f[ u^{\delta}(\boldsymbol{x}_i)=u^{\mathcal{D}}(\boldsymbol{x}_i)+
66  * u^{\mathcal{H}}(\boldsymbol{x}_i)=\sum_{n=0}^{N^{\mathcal{D}}-1}
67  * \hat{u}_n^{\mathcal{D}}\Phi_n(\boldsymbol{x}_i)+
68  * \sum_{n={N^{\mathcal{D}}}}^{N_{\mathrm{dof}}
69  * -1}\hat{u}_n^{\mathcal{H}}
70  * \Phi_n(\boldsymbol{x}_i).\f]
71  * This lifting is accomplished by ordering the known global degrees of
72  * freedom, prescribed by the Dirichlet boundary conditions, first in
73  * the global array \f$\boldsymbol{\hat{u}}\f$, that is,
74  * \f[\boldsymbol{\hat{u}}=\left[ \begin{array}{c}
75  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
76  * \boldsymbol{\hat{u}}^{\mathcal{H}}
77  * \end{array} \right].\f]
78  * Such kind of expansions are also referred to as continuoous fields.
79  * This class should be used when solving 2D problems using a standard
80  * Galerkin approach.
81  */
82 
83  /**
84  * Constructs an empty 1D continuous field.
85  */
86  ContField1D::ContField1D():
87  DisContField1D(),
88  m_locToGloMap(),
89  m_globalLinSysManager(
90  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
91  std::string("GlobalLinSys"))
92  {
93  }
94 
95 
96  /**
97  * Given a mesh \a graph1D, containing information about the domain and
98  * the spectral/hp element expansion, this constructor fills the list
99  * of local expansions #m_exp with the proper expansions, calculates
100  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
101  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates
102  * memory for the arrays #m_coeffs and #m_phys. Furthermore, it
103  * constructs the mapping array (contained in #m_locToGloMap) for the
104  * transformation between local elemental level and global level, it
105  * calculates the total number global expansion coefficients
106  * \f$\hat{u}_n\f$.
107  * The constructor also discretises the boundary conditions, specified
108  * by the argument \a bcs, by expressing them in terms of the
109  * coefficient of the expansion on the boundary.
110  *
111  * @param graph1D A 1D mesh, containing information about the
112  * domain and the spectral/hp element expansion.
113  * @param bcs The boundary conditions.
114  * @param variable An optional parameter to indicate for which
115  * variable the field should be constructed.
116  */
117  ContField1D::ContField1D(const LibUtilities::SessionReaderSharedPtr &pSession,
118  const SpatialDomains::MeshGraphSharedPtr &graph1D,
119  const std::string &variable):
120  DisContField1D(pSession,graph1D,variable,false),
121  m_locToGloMap(),
122  m_globalLinSysManager(
123  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
124  std::string("GlobalLinSys"))
125  {
126  SpatialDomains::BoundaryConditions bcs(pSession, graph1D);
127 
128  m_locToGloMap = MemoryManager<AssemblyMapCG1D>
129  ::AllocateSharedPtr(m_session,m_ncoeffs,*this,
130  m_bndCondExpansions,
131  m_bndConditions,
132  m_periodicVerts,
133  variable);
134  }
135 
136 
137  /**
138  * Constructs a 1D continuous field as a copy of an existing field
139  * including all boundary conditions.
140  * @param In Existing continuous field to duplicate.
141  */
142  ContField1D::ContField1D(const ContField1D &In):
143  DisContField1D(In),
144  m_locToGloMap(In.m_locToGloMap),
145  m_globalLinSysManager(
146  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
147  std::string("GlobalLinSys"))
148  {
149  }
150 
151  /**
152  * Constructs a 1D continuous field as a copy of an existing explist 1D field
153  * and adding all the boundary conditions.
154  * @param In Existing explist1D field .
155  */
156  ContField1D::ContField1D(const LibUtilities::SessionReaderSharedPtr &pSession, const ExpList1D & In):
157  DisContField1D(In),
158  m_locToGloMap(),
159  m_globalLinSysManager(
160  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
161  std::string("GlobalLinSys"))
162  {
163  m_locToGloMap = MemoryManager<AssemblyMapCG1D>
164  ::AllocateSharedPtr(pSession,m_ncoeffs, In);
165 
166  }
167 
168  /**
169  *
170  */
171  ContField1D::~ContField1D()
172  {
173  }
174 
175 
176  /**
177  * Given a function \f$f(x)\f$ defined at the quadrature
178  * points, this function determines the unknown global coefficients
179  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ employing a discrete
180  * Galerkin projection from physical space to coefficient
181  * space. The operation is evaluated by the function #GlobalSolve using
182  * the global mass matrix.
183  *
184  * The values of the function \f$f(x)\f$ evaluated at the
185  * quadrature points \f$x_i\f$ should be contained in the
186  * variable #m_phys of the ExpList object \a Sin. The resulting global
187  * coefficients \f$\hat{u}_g\f$ are stored in the array #m_coeffs.
188  *
189  * @param inarray Discrete \f$f(x)\f$.
190  * @param outarray Storage for result.
191  * @param coeffstate
192  */
193  void ContField1D::FwdTrans(const Array<OneD, const NekDouble> &inarray,
194  Array<OneD, NekDouble> &outarray,
195  CoeffState coeffstate)
196  {
197  // Inner product of forcing
198  Array<OneD,NekDouble> wsp(m_ncoeffs);
199  IProductWRTBase(inarray,wsp,eGlobal);
200 
201  // Solve the system
202  GlobalLinSysKey key(StdRegions::eMass, m_locToGloMap);
203 
204  GlobalSolve(key,wsp,outarray);
205  if(coeffstate == eLocal)
206  {
207  GlobalToLocal(outarray,outarray);
208  }
209  }
210 
211 
212  /**
213  * Given the coefficients of an expansion, this function evaluates the
214  * spectral/hp expansion \f$u^{\delta}(x)\f$ at the quadrature
215  * points \f$x_i\f$. This operation is evaluated locally by the
216  * function ExpList#BwdTrans.
217  *
218  * The coefficients of the expansion should be contained in the
219  * variable #m_coeffs of the ExpList object \a In. The resulting
220  * physical values at the quadrature points \f$u^{\delta}(x_i)\f$ are
221  * stored in the array #m_phys.
222  *
223  * @param In An ExpList, containing the local
224  * coefficients \f$\hat{u}_n^e\f$ in its array
225  * #m_coeffs.
226  */
227  void ContField1D::BwdTrans(
228  const Array<OneD, const NekDouble> &inarray,
229  Array<OneD, NekDouble> &outarray,
230  CoeffState coeffstate)
231  {
232  Array<OneD, NekDouble> tmpinarray;
233  if(coeffstate != eLocal)
234  {
235  tmpinarray = Array<OneD, NekDouble>(inarray);
236  GlobalToLocal(inarray,tmpinarray);
237  }
238  else
239  {
240  tmpinarray = inarray;
241  }
242 
243  BwdTrans_IterPerExp(tmpinarray,outarray);
244  }
245 
246 
247  /**
248  *
249  */
250  void ContField1D::MultiplyByInvMassMatrix(
251  const Array<OneD, const NekDouble> &inarray,
252  Array<OneD, NekDouble> &outarray,
253  CoeffState coeffstate)
254  {
255  GlobalLinSysKey key(StdRegions::eMass, m_locToGloMap);
256  if(coeffstate == eGlobal)
257  {
258  if(inarray.data() == outarray.data())
259  {
260  Array<OneD, NekDouble> tmp(inarray);
261  GlobalSolve(key,tmp,outarray);
262  }
263  else
264  {
265  GlobalSolve(key,inarray,outarray);
266  }
267  }
268  else
269  {
270  Array<OneD, NekDouble> globaltmp(m_ncoeffs,0.0);
271 
272  if(inarray.data() == outarray.data())
273  {
274  Array<OneD,NekDouble> tmp(inarray);
275  Assemble(tmp,outarray);
276  }
277  else
278  {
279  Assemble(inarray,outarray);
280  }
281 
282  GlobalSolve(key,outarray,globaltmp);
283  GlobalToLocal(globaltmp,outarray);
284  }
285  }
286 
287 
288  /**
289  * Given a linear system specified by the key \a key,
290  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}},\f]
291  * this function solves this linear system taking into account the
292  * boundary conditions specified in the data member
293  * #m_bndCondExpansions.
294  * Therefore, it adds an array \f$\boldsymbol{\hat{g}}\f$ which
295  * represents the non-zero surface integral resulting from the weak
296  * boundary conditions (e.g. Neumann boundary conditions) to the right
297  * hand side, that is,
298  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}+
299  * \boldsymbol{\hat{g}}.\f]
300  * Furthermore, it lifts the known degrees of freedom which are
301  * prescribed by the Dirichlet boundary conditions. As these known
302  * coefficients \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ are numbered
303  * first in the global coefficient array \f$\boldsymbol{\hat{u}}_g\f$,
304  * the linear system can be decomposed as,
305  * \f[\left[\begin{array}{cc}
306  * \boldsymbol{M}^{\mathcal{DD}}&\boldsymbol{M}^{\mathcal{DH}}\\
307  * \boldsymbol{M}^{\mathcal{HD}}&\boldsymbol{M}^{\mathcal{HH}}
308  * \end{array}\right]
309  * \left[\begin{array}{c}
310  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
311  * \boldsymbol{\hat{u}}^{\mathcal{H}}
312  * \end{array}\right]=
313  * \left[\begin{array}{c}
314  * \boldsymbol{\hat{f}}^{\mathcal{D}}\\
315  * \boldsymbol{\hat{f}}^{\mathcal{H}}
316  * \end{array}\right]+
317  * \left[\begin{array}{c}
318  * \boldsymbol{\hat{g}}^{\mathcal{D}}\\
319  * \boldsymbol{\hat{g}}^{\mathcal{H}}
320  * \end{array}\right]
321  * \f]
322  * which will then be solved for the unknown coefficients
323  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ as,
324  * \f[
325  * \boldsymbol{M}^{\mathcal{HH}}\boldsymbol{\hat{u}}^{\mathcal{H}}
326  * = \boldsymbol{\hat{f}}^{\mathcal{H}}
327  * +\boldsymbol{\hat{g}}^{\mathcal{H}}
328  * -\boldsymbol{M}^{\mathcal{HD}}\boldsymbol{\hat{u}}^{\mathcal{D}}\f]
329  *
330  * @param key Specifes the linear system to solve.
331  * @param rhs Forcing term \f$\boldsymbol{f}\f$.
332  * @param inout Solution vector \f$\boldsymbol{\hat{u}}\f$.
333  * @param dirForcing .
334  */
335  void ContField1D::GlobalSolve(const GlobalLinSysKey &key,
336  const Array<OneD, const NekDouble>& rhs,
337  Array<OneD, NekDouble>& inout,
338  const Array<OneD, const NekDouble>& dirForcing)
339  {
340  int NumDirBcs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
341  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
342 
343  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE
344  // IN THE SOLUTION ARRAY
345  v_ImposeDirichletConditions(inout);
346 
347  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
348  {
349  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
350  {
351  inout[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
352  = m_bndCondExpansions[i]->GetCoeff(0);
353  }
354  }
355 
356  // STEP 2: CALCULATE THE HOMOGENEOUS COEFFICIENTS
357  if(contNcoeffs - NumDirBcs > 0)
358  {
359  GlobalLinSysSharedPtr LinSys = GetGlobalLinSys(key);
360  LinSys->Solve(rhs,inout,m_locToGloMap,dirForcing);
361  }
362  }
363 
364 
365  /**
366  * The function searches the map #m_globalLinSys to see if the global
367  * matrix has been created before. If not, it calls the function
368  * #GenglobalLinSys to generate the requested global system.
369  *
370  * @param mkey Key specifying the linear system.
371  * @returns Pointer to the required linear system.
372  */
373  GlobalLinSysSharedPtr ContField1D::GetGlobalLinSys(
374  const GlobalLinSysKey &mkey)
375  {
376  return m_globalLinSysManager[mkey];
377  }
378 
379  GlobalLinSysSharedPtr ContField1D::GenGlobalLinSys(
380  const GlobalLinSysKey &mkey)
381  {
382  ASSERTL1(mkey.LocToGloMapIsDefined(),
383  "To use method must have a AssemblyMap "
384  "attached to key");
385  return ExpList::GenGlobalLinSys(mkey, m_locToGloMap);
386  }
387 
388 
389  /**
390  * The operation is evaluated locally (i.e. with respect to all local
391  * expansion modes) by the function ExpList#IProductWRTBase. The inner
392  * product with respect to the global expansion modes is than obtained
393  * by a global assembly operation.
394  *
395  * The values of the function \f$f(x)\f$ evaluated at the quadrature
396  * points \f$x_i\f$ should be contained in the variable #m_phys of the
397  * ExpList object \a in. The result is stored in the array
398  * #m_coeffs.
399  *
400  * @param In An ExpList, containing the discrete evaluation
401  * of \f$f(x)\f$ at the quadrature points in its
402  * array #m_phys.
403  */
404  void ContField1D::IProductWRTBase(
405  const Array<OneD, const NekDouble> &inarray,
406  Array<OneD, NekDouble> &outarray,
407  CoeffState coeffstate)
408  {
409  if(coeffstate == eGlobal)
410  {
411  Array<OneD, NekDouble> wsp(m_ncoeffs);
412  IProductWRTBase_IterPerExp(inarray,wsp);
413  Assemble(wsp,outarray);
414  }
415  else
416  {
417  IProductWRTBase_IterPerExp(inarray,outarray);
418  }
419  }
420 
421 
422  void ContField1D::v_FwdTrans(
423  const Array<OneD, const NekDouble> &inarray,
424  Array<OneD, NekDouble> &outarray,
425  CoeffState coeffstate)
426  {
427  FwdTrans(inarray,outarray,coeffstate);
428  }
429 
430  void ContField1D::v_MultiplyByInvMassMatrix(
431  const Array<OneD, const NekDouble> &inarray,
432  Array<OneD, NekDouble> &outarray,
433  CoeffState coeffstate)
434  {
435  MultiplyByInvMassMatrix(inarray,outarray,coeffstate);
436  }
437 
438  void ContField1D::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
439  {
440  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
441  {
442  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
443  {
444  outarray[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
445  = m_bndCondExpansions[i]->GetCoeff(0);
446  }
447  }
448  }
449 
450  /**
451  * This operation is evaluated as:
452  * \f{tabbing}
453  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
454  * > > Do \= $i=$ $0,N_m^e-1$ \\
455  * > > > $\boldsymbol{\hat{u}}_g[\mbox{map}[e][i]] =
456  * \mbox{sign}[e][i] \cdot \boldsymbol{\hat{u}}^{e}[i]$\\
457  * > > continue\\
458  * > continue
459  * \f}
460  * where \a map\f$[e][i]\f$ is the mapping array and \a
461  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
462  * correct modal connectivity between the different elements (both
463  * these arrays are contained in the data member #m_locToGloMap). This
464  * operation is equivalent to the gather operation
465  * \f$\boldsymbol{\hat{u}}_g=\mathcal{A}^{-1}\boldsymbol{\hat{u}}_l\f$,
466  * where \f$\mathcal{A}\f$ is the
467  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
468  *
469  */
470  void ContField1D::v_LocalToGlobal(void)
471  {
472  m_locToGloMap->LocalToGlobal(m_coeffs,m_coeffs);
473  }
474 
475  /**
476  * This operation is evaluated as:
477  * \f{tabbing}
478  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
479  * > > Do \= $i=$ $0,N_m^e-1$ \\
480  * > > > $\boldsymbol{\hat{u}}^{e}[i] = \mbox{sign}[e][i] \cdot
481  * \boldsymbol{\hat{u}}_g[\mbox{map}[e][i]]$ \\
482  * > > continue \\
483  * > continue
484  * \f}
485  * where \a map\f$[e][i]\f$ is the mapping array and
486  * \a sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
487  * correct modal connectivity between the different elements (both
488  * these arrays are contained in the data member #m_locToGloMap). This
489  * operation is equivalent to the scatter operation
490  * \f$\boldsymbol{\hat{u}}_l=\mathcal{A}\boldsymbol{\hat{u}}_g\f$, where
491  * \f$\mathcal{A}\f$ is the
492  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
493  *
494  */
495  void ContField1D::v_GlobalToLocal(void)
496  {
497  m_locToGloMap->GlobalToLocal(m_coeffs,m_coeffs);
498  }
499 
500  /**
501  * Consider the one dimensional Helmholtz equation,
502  * \f[\frac{d^2u}{dx^2}-\lambda u(x) = f(x),\f]
503  * supplemented with appropriate boundary conditions (which are
504  * contained in the data member #m_bndCondExpansions). Applying a
505  * \f$C^0\f$ continuous Galerkin discretisation, this equation leads to
506  * the following linear system:
507  * \f[\left( \boldsymbol{M}+\lambda\boldsymbol{L}\right)
508  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f]
509  * where \f$\boldsymbol{M}\f$ and \f$\boldsymbol{L}\f$ are the mass and
510  * Laplacian matrix respectively. This function solves the system above
511  * for the global coefficients \f$\boldsymbol{\hat{u}}\f$ by a call to
512  * the function #GlobalSolve.
513  *
514  * The values of the function \f$f(x)\f$ evaluated at the
515  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
516  * variable #m_phys of the ExpList object \a inarray. The resulting
517  * global coefficients \f$\boldsymbol{\hat{u}}_g\f$ are stored in the
518  * array #m_coeffs.
519  *
520  * @param inarray Input containing forcing function
521  * \f$\boldsymbol{f}\f$ at the quadrature points.
522  * @param outarray Output containing the coefficients
523  * \f$\boldsymbol{u}_g\f$
524  * @param lambda Parameter value.
525  * @param Sigma Coefficients of lambda.
526  * @param varcoeff Variable diffusivity coefficients.
527  * @param coeffstate
528  * @param dirForcing Dirichlet Forcing.
529  */
530  void ContField1D::v_HelmSolve(
531  const Array<OneD, const NekDouble> &inarray,
532  Array<OneD, NekDouble> &outarray,
533  const FlagList &flags,
534  const StdRegions::ConstFactorMap &factors,
535  const StdRegions::VarCoeffMap &varcoeff,
536  const Array<OneD, const NekDouble> &dirForcing)
537  {
538  // Inner product of forcing
539  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
540  Array<OneD,NekDouble> wsp(contNcoeffs);
541  IProductWRTBase(inarray,wsp,eGlobal);
542  // Note -1.0 term necessary to invert forcing function to
543  // be consistent with matrix definition
544  Vmath::Neg(contNcoeffs, wsp, 1);
545 
546  // Forcing function with weak boundary conditions
547  int i;
548  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
549  {
550  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
551  {
552  wsp[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
553  += m_bndCondExpansions[i]->GetCoeff(0);
554  }
555  }
556 
557  // Solve the system
558  GlobalLinSysKey key(StdRegions::eHelmholtz,
559  m_locToGloMap,factors,varcoeff);
560 
561  if(flags.isSet(eUseGlobal))
562  {
563  GlobalSolve(key,wsp,outarray,dirForcing);
564  }
565  else
566  {
567  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
568  GlobalSolve(key,wsp,tmp,dirForcing);
569  GlobalToLocal(tmp,outarray);
570  }
571  }
572 
573  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>&
574  ContField1D::v_GetBndConditions()
575  {
576  return GetBndConditions();
577  }
578 
579  void ContField1D::v_BwdTrans(
580  const Array<OneD, const NekDouble> &inarray,
581  Array<OneD, NekDouble> &outarray,
582  CoeffState coeffstate)
583  {
584  BwdTrans(inarray,outarray,coeffstate);
585  }
586 
587  void ContField1D::v_IProductWRTBase(
588  const Array<OneD, const NekDouble> &inarray,
589  Array<OneD, NekDouble> &outarray,
590  CoeffState coeffstate)
591  {
592  IProductWRTBase(inarray,outarray,coeffstate);
593  }
594 
595  /**
596  * This is equivalent to the operation:
597  * \f[\boldsymbol{M\hat{u}}_g\f]
598  * where \f$\boldsymbol{M}\f$ is the global matrix of type specified by
599  * \a mkey. After scattering the global array \a inarray to local
600  * level, this operation is evaluated locally by the function
601  * ExpList#GeneralMatrixOp. The global result is then obtained by a
602  * global assembly procedure.
603  *
604  * @param mkey This key uniquely defines the type matrix
605  * required for the operation.
606  * @param inarray The vector \f$\boldsymbol{\hat{u}}_g\f$ of size
607  * \f$N_{\mathrm{dof}}\f$.
608  * @param outarray The resulting vector of size
609  * \f$N_{\mathrm{dof}}\f$.
610  */
611  void ContField1D::v_GeneralMatrixOp(
612  const GlobalMatrixKey &gkey,
613  const Array<OneD,const NekDouble> &inarray,
614  Array<OneD, NekDouble> &outarray,
615  CoeffState coeffstate)
616  {
617  if(coeffstate == eGlobal)
618  {
619  Array<OneD,NekDouble> tmp1(2*m_ncoeffs);
620  Array<OneD,NekDouble> tmp2(tmp1+m_ncoeffs);
621  GlobalToLocal(inarray,tmp1);
622  GeneralMatrixOp_IterPerExp(gkey,tmp1,tmp2);
623  Assemble(tmp2,outarray);
624  }
625  else
626  {
627  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
628  }
629  }
630 
631  } // end of namespace
632 } //end of namespace
Generated on Sun Mar 15 2015 00:11:30 for Nektar++ by   doxygen 1.8.1.2