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/AssemblyMapCG.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<AssemblyMapCG>
129  ::AllocateSharedPtr(m_session,m_ncoeffs,*this,
130  m_bndCondExpansions,
131  m_bndConditions,
132  false,
133  variable,
134  m_periodicVerts);
135  }
136 
137 
138  /**
139  * Constructs a 1D continuous field as a copy of an existing field
140  * including all boundary conditions.
141  * @param In Existing continuous field to duplicate.
142  */
143  ContField1D::ContField1D(const ContField1D &In):
144  DisContField1D(In),
145  m_locToGloMap(In.m_locToGloMap),
146  m_globalLinSysManager(
147  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
148  std::string("GlobalLinSys"))
149  {
150  }
151 
152  /**
153  * Constructs a 1D continuous field as a copy of an existing explist 1D field
154  * and adding all the boundary conditions.
155  * @param In Existing explist1D field .
156  */
157  ContField1D::ContField1D(const LibUtilities::SessionReaderSharedPtr &pSession, const ExpList1D & In):
158  DisContField1D(In),
159  m_locToGloMap(),
160  m_globalLinSysManager(
161  boost::bind(&ContField1D::GenGlobalLinSys, this, _1),
162  std::string("GlobalLinSys"))
163  {
164  m_locToGloMap = MemoryManager<AssemblyMapCG>
165  ::AllocateSharedPtr(pSession, m_ncoeffs, In);
166 
167  }
168 
169  /**
170  *
171  */
172  ContField1D::~ContField1D()
173  {
174  }
175 
176 
177  /**
178  * Given a function \f$f(x)\f$ defined at the quadrature
179  * points, this function determines the unknown global coefficients
180  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ employing a discrete
181  * Galerkin projection from physical space to coefficient
182  * space. The operation is evaluated by the function #GlobalSolve using
183  * the global mass matrix.
184  *
185  * The values of the function \f$f(x)\f$ evaluated at the
186  * quadrature points \f$x_i\f$ should be contained in the
187  * variable #m_phys of the ExpList object \a Sin. The resulting global
188  * coefficients \f$\hat{u}_g\f$ are stored in the array #m_coeffs.
189  *
190  * @param inarray Discrete \f$f(x)\f$.
191  * @param outarray Storage for result.
192  * @param coeffstate
193  */
194  void ContField1D::FwdTrans(const Array<OneD, const NekDouble> &inarray,
195  Array<OneD, NekDouble> &outarray,
196  CoeffState coeffstate)
197  {
198  // Inner product of forcing
199  Array<OneD,NekDouble> wsp(m_ncoeffs);
200  IProductWRTBase(inarray,wsp,eGlobal);
201 
202  // Solve the system
203  GlobalLinSysKey key(StdRegions::eMass, m_locToGloMap);
204 
205  GlobalSolve(key,wsp,outarray);
206  if(coeffstate == eLocal)
207  {
208  GlobalToLocal(outarray,outarray);
209  }
210  }
211 
212 
213  /**
214  * Given the coefficients of an expansion, this function evaluates the
215  * spectral/hp expansion \f$u^{\delta}(x)\f$ at the quadrature
216  * points \f$x_i\f$. This operation is evaluated locally by the
217  * function ExpList#BwdTrans.
218  *
219  * The coefficients of the expansion should be contained in the
220  * variable #m_coeffs of the ExpList object \a In. The resulting
221  * physical values at the quadrature points \f$u^{\delta}(x_i)\f$ are
222  * stored in the array #m_phys.
223  *
224  * @param In An ExpList, containing the local
225  * coefficients \f$\hat{u}_n^e\f$ in its array
226  * #m_coeffs.
227  */
228  void ContField1D::BwdTrans(
229  const Array<OneD, const NekDouble> &inarray,
230  Array<OneD, NekDouble> &outarray,
231  CoeffState coeffstate)
232  {
233  Array<OneD, NekDouble> tmpinarray;
234  if(coeffstate != eLocal)
235  {
236  tmpinarray = Array<OneD, NekDouble>(inarray);
237  GlobalToLocal(inarray,tmpinarray);
238  }
239  else
240  {
241  tmpinarray = inarray;
242  }
243 
244  BwdTrans_IterPerExp(tmpinarray,outarray);
245  }
246 
247 
248  /**
249  *
250  */
251  void ContField1D::MultiplyByInvMassMatrix(
252  const Array<OneD, const NekDouble> &inarray,
253  Array<OneD, NekDouble> &outarray,
254  CoeffState coeffstate)
255  {
256  GlobalLinSysKey key(StdRegions::eMass, m_locToGloMap);
257  if(coeffstate == eGlobal)
258  {
259  if(inarray.data() == outarray.data())
260  {
261  Array<OneD, NekDouble> tmp(inarray);
262  GlobalSolve(key,tmp,outarray);
263  }
264  else
265  {
266  GlobalSolve(key,inarray,outarray);
267  }
268  }
269  else
270  {
271  Array<OneD, NekDouble> globaltmp(m_ncoeffs,0.0);
272 
273  if(inarray.data() == outarray.data())
274  {
275  Array<OneD,NekDouble> tmp(inarray);
276  Assemble(tmp,outarray);
277  }
278  else
279  {
280  Assemble(inarray,outarray);
281  }
282 
283  GlobalSolve(key,outarray,globaltmp);
284  GlobalToLocal(globaltmp,outarray);
285  }
286  }
287 
288 
289  /**
290  * Given a linear system specified by the key \a key,
291  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}},\f]
292  * this function solves this linear system taking into account the
293  * boundary conditions specified in the data member
294  * #m_bndCondExpansions.
295  * Therefore, it adds an array \f$\boldsymbol{\hat{g}}\f$ which
296  * represents the non-zero surface integral resulting from the weak
297  * boundary conditions (e.g. Neumann boundary conditions) to the right
298  * hand side, that is,
299  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}+
300  * \boldsymbol{\hat{g}}.\f]
301  * Furthermore, it lifts the known degrees of freedom which are
302  * prescribed by the Dirichlet boundary conditions. As these known
303  * coefficients \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ are numbered
304  * first in the global coefficient array \f$\boldsymbol{\hat{u}}_g\f$,
305  * the linear system can be decomposed as,
306  * \f[\left[\begin{array}{cc}
307  * \boldsymbol{M}^{\mathcal{DD}}&\boldsymbol{M}^{\mathcal{DH}}\\
308  * \boldsymbol{M}^{\mathcal{HD}}&\boldsymbol{M}^{\mathcal{HH}}
309  * \end{array}\right]
310  * \left[\begin{array}{c}
311  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
312  * \boldsymbol{\hat{u}}^{\mathcal{H}}
313  * \end{array}\right]=
314  * \left[\begin{array}{c}
315  * \boldsymbol{\hat{f}}^{\mathcal{D}}\\
316  * \boldsymbol{\hat{f}}^{\mathcal{H}}
317  * \end{array}\right]+
318  * \left[\begin{array}{c}
319  * \boldsymbol{\hat{g}}^{\mathcal{D}}\\
320  * \boldsymbol{\hat{g}}^{\mathcal{H}}
321  * \end{array}\right]
322  * \f]
323  * which will then be solved for the unknown coefficients
324  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ as,
325  * \f[
326  * \boldsymbol{M}^{\mathcal{HH}}\boldsymbol{\hat{u}}^{\mathcal{H}}
327  * = \boldsymbol{\hat{f}}^{\mathcal{H}}
328  * +\boldsymbol{\hat{g}}^{\mathcal{H}}
329  * -\boldsymbol{M}^{\mathcal{HD}}\boldsymbol{\hat{u}}^{\mathcal{D}}\f]
330  *
331  * @param key Specifes the linear system to solve.
332  * @param rhs Forcing term \f$\boldsymbol{f}\f$.
333  * @param inout Solution vector \f$\boldsymbol{\hat{u}}\f$.
334  * @param dirForcing .
335  */
336  void ContField1D::GlobalSolve(const GlobalLinSysKey &key,
337  const Array<OneD, const NekDouble>& rhs,
338  Array<OneD, NekDouble>& inout,
339  const Array<OneD, const NekDouble>& dirForcing)
340  {
341  int NumDirBcs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
342  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
343 
344  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE
345  // IN THE SOLUTION ARRAY
346  v_ImposeDirichletConditions(inout);
347 
348  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
349  {
350  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
351  {
352  inout[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
353  = m_bndCondExpansions[i]->GetCoeff(0);
354  }
355  }
356 
357  // STEP 2: CALCULATE THE HOMOGENEOUS COEFFICIENTS
358  if(contNcoeffs - NumDirBcs > 0)
359  {
360  GlobalLinSysSharedPtr LinSys = GetGlobalLinSys(key);
361  LinSys->Solve(rhs,inout,m_locToGloMap,dirForcing);
362  }
363  }
364 
365 
366  /**
367  * The function searches the map #m_globalLinSys to see if the global
368  * matrix has been created before. If not, it calls the function
369  * #GenglobalLinSys to generate the requested global system.
370  *
371  * @param mkey Key specifying the linear system.
372  * @returns Pointer to the required linear system.
373  */
374  GlobalLinSysSharedPtr ContField1D::GetGlobalLinSys(
375  const GlobalLinSysKey &mkey)
376  {
377  return m_globalLinSysManager[mkey];
378  }
379 
380  GlobalLinSysSharedPtr ContField1D::GenGlobalLinSys(
381  const GlobalLinSysKey &mkey)
382  {
383  ASSERTL1(mkey.LocToGloMapIsDefined(),
384  "To use method must have a AssemblyMap "
385  "attached to key");
386  return ExpList::GenGlobalLinSys(mkey, m_locToGloMap);
387  }
388 
389 
390  /**
391  * The operation is evaluated locally (i.e. with respect to all local
392  * expansion modes) by the function ExpList#IProductWRTBase. The inner
393  * product with respect to the global expansion modes is than obtained
394  * by a global assembly operation.
395  *
396  * The values of the function \f$f(x)\f$ evaluated at the quadrature
397  * points \f$x_i\f$ should be contained in the variable #m_phys of the
398  * ExpList object \a in. The result is stored in the array
399  * #m_coeffs.
400  *
401  * @param In An ExpList, containing the discrete evaluation
402  * of \f$f(x)\f$ at the quadrature points in its
403  * array #m_phys.
404  */
405  void ContField1D::IProductWRTBase(
406  const Array<OneD, const NekDouble> &inarray,
407  Array<OneD, NekDouble> &outarray,
408  CoeffState coeffstate)
409  {
410  if(coeffstate == eGlobal)
411  {
412  Array<OneD, NekDouble> wsp(m_ncoeffs);
413  IProductWRTBase_IterPerExp(inarray,wsp);
414  Assemble(wsp,outarray);
415  }
416  else
417  {
418  IProductWRTBase_IterPerExp(inarray,outarray);
419  }
420  }
421 
422 
423  void ContField1D::v_FwdTrans(
424  const Array<OneD, const NekDouble> &inarray,
425  Array<OneD, NekDouble> &outarray,
426  CoeffState coeffstate)
427  {
428  FwdTrans(inarray,outarray,coeffstate);
429  }
430 
431  void ContField1D::v_MultiplyByInvMassMatrix(
432  const Array<OneD, const NekDouble> &inarray,
433  Array<OneD, NekDouble> &outarray,
434  CoeffState coeffstate)
435  {
436  MultiplyByInvMassMatrix(inarray,outarray,coeffstate);
437  }
438 
439  void ContField1D::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
440  {
441  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
442  {
443  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
444  {
445  outarray[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
446  = m_bndCondExpansions[i]->GetCoeff(0);
447  }
448  }
449  }
450 
451  /**
452  * This operation is evaluated as:
453  * \f{tabbing}
454  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
455  * > > Do \= $i=$ $0,N_m^e-1$ \\
456  * > > > $\boldsymbol{\hat{u}}_g[\mbox{map}[e][i]] =
457  * \mbox{sign}[e][i] \cdot \boldsymbol{\hat{u}}^{e}[i]$\\
458  * > > continue\\
459  * > continue
460  * \f}
461  * where \a map\f$[e][i]\f$ is the mapping array and \a
462  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
463  * correct modal connectivity between the different elements (both
464  * these arrays are contained in the data member #m_locToGloMap). This
465  * operation is equivalent to the gather operation
466  * \f$\boldsymbol{\hat{u}}_g=\mathcal{A}^{-1}\boldsymbol{\hat{u}}_l\f$,
467  * where \f$\mathcal{A}\f$ is the
468  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
469  *
470  */
471  void ContField1D::v_LocalToGlobal(void)
472  {
473  m_locToGloMap->LocalToGlobal(m_coeffs,m_coeffs);
474  }
475 
476  /**
477  * This operation is evaluated as:
478  * \f{tabbing}
479  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
480  * > > Do \= $i=$ $0,N_m^e-1$ \\
481  * > > > $\boldsymbol{\hat{u}}^{e}[i] = \mbox{sign}[e][i] \cdot
482  * \boldsymbol{\hat{u}}_g[\mbox{map}[e][i]]$ \\
483  * > > continue \\
484  * > continue
485  * \f}
486  * where \a map\f$[e][i]\f$ is the mapping array and
487  * \a sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
488  * correct modal connectivity between the different elements (both
489  * these arrays are contained in the data member #m_locToGloMap). This
490  * operation is equivalent to the scatter operation
491  * \f$\boldsymbol{\hat{u}}_l=\mathcal{A}\boldsymbol{\hat{u}}_g\f$, where
492  * \f$\mathcal{A}\f$ is the
493  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
494  *
495  */
496  void ContField1D::v_GlobalToLocal(void)
497  {
498  m_locToGloMap->GlobalToLocal(m_coeffs,m_coeffs);
499  }
500 
501  /**
502  * Consider the one dimensional Helmholtz equation,
503  * \f[\frac{d^2u}{dx^2}-\lambda u(x) = f(x),\f]
504  * supplemented with appropriate boundary conditions (which are
505  * contained in the data member #m_bndCondExpansions). Applying a
506  * \f$C^0\f$ continuous Galerkin discretisation, this equation leads to
507  * the following linear system:
508  * \f[\left( \boldsymbol{M}+\lambda\boldsymbol{L}\right)
509  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f]
510  * where \f$\boldsymbol{M}\f$ and \f$\boldsymbol{L}\f$ are the mass and
511  * Laplacian matrix respectively. This function solves the system above
512  * for the global coefficients \f$\boldsymbol{\hat{u}}\f$ by a call to
513  * the function #GlobalSolve.
514  *
515  * The values of the function \f$f(x)\f$ evaluated at the
516  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
517  * variable #m_phys of the ExpList object \a inarray. The resulting
518  * global coefficients \f$\boldsymbol{\hat{u}}_g\f$ are stored in the
519  * array #m_coeffs.
520  *
521  * @param inarray Input containing forcing function
522  * \f$\boldsymbol{f}\f$ at the quadrature points.
523  * @param outarray Output containing the coefficients
524  * \f$\boldsymbol{u}_g\f$
525  * @param lambda Parameter value.
526  * @param Sigma Coefficients of lambda.
527  * @param varcoeff Variable diffusivity coefficients.
528  * @param coeffstate
529  * @param dirForcing Dirichlet Forcing.
530  */
531  void ContField1D::v_HelmSolve(
532  const Array<OneD, const NekDouble> &inarray,
533  Array<OneD, NekDouble> &outarray,
534  const FlagList &flags,
535  const StdRegions::ConstFactorMap &factors,
536  const StdRegions::VarCoeffMap &varcoeff,
537  const Array<OneD, const NekDouble> &dirForcing)
538  {
539  // Inner product of forcing
540  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
541  Array<OneD,NekDouble> wsp(contNcoeffs);
542  IProductWRTBase(inarray,wsp,eGlobal);
543  // Note -1.0 term necessary to invert forcing function to
544  // be consistent with matrix definition
545  Vmath::Neg(contNcoeffs, wsp, 1);
546 
547  // Forcing function with weak boundary conditions
548  int i;
549  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
550  {
551  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
552  {
553  wsp[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(i)]
554  += m_bndCondExpansions[i]->GetCoeff(0);
555  }
556  }
557 
558  // Solve the system
559  GlobalLinSysKey key(StdRegions::eHelmholtz,
560  m_locToGloMap,factors,varcoeff);
561 
562  if(flags.isSet(eUseGlobal))
563  {
564  GlobalSolve(key,wsp,outarray,dirForcing);
565  }
566  else
567  {
568  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
569  GlobalSolve(key,wsp,tmp,dirForcing);
570  GlobalToLocal(tmp,outarray);
571  }
572  }
573 
574  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>&
575  ContField1D::v_GetBndConditions()
576  {
577  return GetBndConditions();
578  }
579 
580  void ContField1D::v_BwdTrans(
581  const Array<OneD, const NekDouble> &inarray,
582  Array<OneD, NekDouble> &outarray,
583  CoeffState coeffstate)
584  {
585  BwdTrans(inarray,outarray,coeffstate);
586  }
587 
588  void ContField1D::v_IProductWRTBase(
589  const Array<OneD, const NekDouble> &inarray,
590  Array<OneD, NekDouble> &outarray,
591  CoeffState coeffstate)
592  {
593  IProductWRTBase(inarray,outarray,coeffstate);
594  }
595 
596  /**
597  * This is equivalent to the operation:
598  * \f[\boldsymbol{M\hat{u}}_g\f]
599  * where \f$\boldsymbol{M}\f$ is the global matrix of type specified by
600  * \a mkey. After scattering the global array \a inarray to local
601  * level, this operation is evaluated locally by the function
602  * ExpList#GeneralMatrixOp. The global result is then obtained by a
603  * global assembly procedure.
604  *
605  * @param mkey This key uniquely defines the type matrix
606  * required for the operation.
607  * @param inarray The vector \f$\boldsymbol{\hat{u}}_g\f$ of size
608  * \f$N_{\mathrm{dof}}\f$.
609  * @param outarray The resulting vector of size
610  * \f$N_{\mathrm{dof}}\f$.
611  */
612  void ContField1D::v_GeneralMatrixOp(
613  const GlobalMatrixKey &gkey,
614  const Array<OneD,const NekDouble> &inarray,
615  Array<OneD, NekDouble> &outarray,
616  CoeffState coeffstate)
617  {
618  if(coeffstate == eGlobal)
619  {
620  Array<OneD,NekDouble> tmp1(2*m_ncoeffs);
621  Array<OneD,NekDouble> tmp2(tmp1+m_ncoeffs);
622  GlobalToLocal(inarray,tmp1);
623  GeneralMatrixOp_IterPerExp(gkey,tmp1,tmp2);
624  Assemble(tmp2,outarray);
625  }
626  else
627  {
628  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
629  }
630  }
631 
632  } // end of namespace
633 } //end of namespace
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2