Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ContField2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ContField2D.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 2D domain with boundary conditions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/ContField2D.h>
37 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h>
38 
39 namespace Nektar
40 {
41  namespace MultiRegions
42  {
43  /**
44  * @class ContField2D
45  * The class #ContField2D 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, contains the
54  * one-dimensional spectral/hp 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, this class is
60  * capable of lifting a known solution satisfying these boundary
61  * 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}}-1}
69  * \hat{u}_n^{\mathcal{H}} \Phi_n(\boldsymbol{x}_i).\f]
70  * This lifting is accomplished by ordering the known global degrees of
71  * freedom, prescribed by the Dirichlet boundary conditions, first in
72  * the global array
73  * \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 continuous fields.
79  * This class should be used when solving 2D problems using a standard
80  * Galerkin approach.
81  */
82 
83  /**
84  *
85  */
86  ContField2D::ContField2D():
87  DisContField2D(),
88  m_locToGloMap(),
89  m_globalMat(),
90  m_globalLinSysManager(
91  boost::bind(&ContField2D::GenGlobalLinSys, this, _1),
92  std::string("GlobalLinSys"))
93  {
94  }
95 
96 
97  /**
98  * Given a mesh \a graph2D, 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 graph2D 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  */
118  ContField2D::ContField2D(const LibUtilities::SessionReaderSharedPtr &pSession,
119  const SpatialDomains::MeshGraphSharedPtr &graph2D,
120  const std::string &variable,
121  const bool DeclareCoeffPhysArrays,
122  const bool CheckIfSingularSystem):
123  DisContField2D(pSession,graph2D,variable,false,DeclareCoeffPhysArrays),
124  m_globalMat(MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
125  m_globalLinSysManager(
126  boost::bind(&ContField2D::GenGlobalLinSys, this, _1),
127  std::string("GlobalLinSys"))
128  {
129  m_locToGloMap = MemoryManager<AssemblyMapCG>
130  ::AllocateSharedPtr(m_session,m_ncoeffs,*this,
131  m_bndCondExpansions,
132  m_bndConditions,
133  CheckIfSingularSystem,
134  variable,
135  m_periodicVerts,
136  m_periodicEdges);
137 
138  if (m_session->DefinesCmdLineArgument("verbose"))
139  {
140  m_locToGloMap->PrintStats(std::cout, variable);
141  }
142  }
143 
144 
145  /**
146  * Given a mesh \a graph2D, containing information about the domain and
147  * the spectral/hp element expansion, this constructor fills the list
148  * of local expansions #m_exp with the proper expansions, calculates
149  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
150  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory
151  * for the arrays #m_coeffs and #m_phys. Furthermore, it constructs the
152  * mapping array (contained in #m_locToGloMap) for the transformation
153  * between local elemental level and global level, it calculates the
154  * total number global expansion coefficients \f$\hat{u}_n\f$ and
155  * allocates memory for the array #m_coeffs. The constructor also
156  * discretises the boundary conditions, specified by the argument \a
157  * bcs, by expressing them in terms of the coefficient of the expansion
158  * on the boundary.
159  *
160  * @param In Existing ContField2D object used to provide the
161  * local to global mapping information and
162  * global solution type.
163  * @param graph2D A mesh, containing information about the domain
164  * and the spectral/hp element expansion.
165  * @param bcs The boundary conditions.
166  * @param bc_loc
167  */
168  ContField2D::ContField2D(const ContField2D &In,
169  const SpatialDomains::MeshGraphSharedPtr &graph2D,
170  const std::string &variable,
171  bool DeclareCoeffPhysArrays,
172  const bool CheckIfSingularSystem):
173  DisContField2D(In,graph2D,variable,false,DeclareCoeffPhysArrays),
174  m_globalMat (MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
175  m_globalLinSysManager(
176  boost::bind(&ContField2D::GenGlobalLinSys, this, _1),
177  std::string("GlobalLinSys"))
178  {
179  if(!SameTypeOfBoundaryConditions(In) || CheckIfSingularSystem)
180  {
181  m_locToGloMap = MemoryManager<AssemblyMapCG>
182  ::AllocateSharedPtr(m_session, m_ncoeffs,*this,
183  m_bndCondExpansions,
184  m_bndConditions,
185  CheckIfSingularSystem,
186  variable,
187  m_periodicVerts,
188  m_periodicEdges);
189 
190  if (m_session->DefinesCmdLineArgument("verbose"))
191  {
192  m_locToGloMap->PrintStats(std::cout, variable);
193  }
194  }
195  else
196  {
197  m_locToGloMap = In.m_locToGloMap;
198  }
199  }
200 
201 
202  /**
203  * Initialises the object as a copy of an existing ContField2D object.
204  * @param In Existing ContField2D object.
205  * @param DeclareCoeffPhysArrays bool to declare if \a m_phys
206  * and \a m_coeffs should be declared. Default is true
207  */
208  ContField2D::ContField2D(const ContField2D &In, bool DeclareCoeffPhysArrays):
209  DisContField2D(In,DeclareCoeffPhysArrays),
210  m_locToGloMap(In.m_locToGloMap),
211  m_globalMat(In.m_globalMat),
212  m_globalLinSysManager(In.m_globalLinSysManager)
213  {
214  }
215 
216 
217  /**
218  *
219  */
220  ContField2D::~ContField2D()
221  {
222  }
223 
224 
225  /**
226  * Given a function \f$f(\boldsymbol{x})\f$ defined at the quadrature
227  * points, this function determines the unknown global coefficients
228  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ employing a discrete
229  * Galerkin projection from physical space to coefficient
230  * space. The operation is evaluated by the function #GlobalSolve using
231  * the global mass matrix.
232  *
233  * The values of the function \f$f(\boldsymbol{x})\f$ evaluated at the
234  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
235  * variable #m_phys of the ExpList object \a Sin. The resulting global
236  * coefficients \f$\hat{u}_g\f$ are stored in the array #m_coeffs.
237  *
238  * @param Sin An ExpList, containing the discrete evaluation
239  * of \f$f(\boldsymbol{x})\f$ at the quadrature
240  * points in its array #m_phys.
241  */
242  void ContField2D::FwdTrans(const Array<OneD, const NekDouble> &inarray,
243  Array<OneD, NekDouble> &outarray,
244  CoeffState coeffstate)
245 
246  {
247  // Inner product of forcing
248  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
249  Array<OneD,NekDouble> wsp(contNcoeffs);
250  IProductWRTBase(inarray,wsp,eGlobal);
251 
252  // Solve the system
253  GlobalLinSysKey key(StdRegions::eMass, m_locToGloMap);
254 
255  if(coeffstate == eGlobal)
256  {
257  GlobalSolve(key,wsp,outarray);
258  }
259  else
260  {
261  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
262  GlobalSolve(key,wsp,tmp);
263  GlobalToLocal(tmp,outarray);
264  }
265  }
266 
267  /**
268  *
269  */
270  void ContField2D::v_SmoothField(Array<OneD,NekDouble> &field)
271  {
272  int gloNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
273  Array<OneD,NekDouble> tmp1(gloNcoeffs);
274  Array<OneD,NekDouble> tmp2(gloNcoeffs);
275 
276  IProductWRTBase(field,tmp1,eGlobal);
277  MultiplyByInvMassMatrix(tmp1,tmp2,eGlobal);
278  BwdTrans(tmp2,field,eGlobal);
279  }
280 
281 
282  /**
283  * Computes the matrix vector product
284  * @f$ \mathbf{y} = \mathbf{M}^{-1}\mathbf{x} @f$. If \a coeffstate == eGlobal
285  * is set then the elemental system is used directly. If not set, the
286  * global system is assembled, the system is solved, and mapped back to
287  * the local elemental system.
288  *
289  * @param inarray Input vector @f$\mathbf{x}@f$.
290  * @param outarray Output vector @f$\mathbf{y}@f$.
291  * @param coeffState Flag for using global system.
292  */
293  void ContField2D::MultiplyByInvMassMatrix(
294  const Array<OneD, const NekDouble> &inarray,
295  Array<OneD, NekDouble> &outarray,
296  CoeffState coeffstate)
297 
298  {
299  GlobalLinSysKey key(StdRegions::eMass,m_locToGloMap);
300  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
301 
302  if(coeffstate == eGlobal)
303  {
304  if(inarray.data() == outarray.data())
305  {
306  Array<OneD, NekDouble> tmp(contNcoeffs,0.0);
307  Vmath::Vcopy(contNcoeffs,inarray,1,tmp,1);
308  GlobalSolve(key,tmp,outarray);
309  }
310  else
311  {
312  GlobalSolve(key,inarray,outarray);
313  }
314  }
315  else
316  {
317  Array<OneD, NekDouble> globaltmp(contNcoeffs,0.0);
318 
319  if(inarray.data() == outarray.data())
320  {
321  Array<OneD,NekDouble> tmp(inarray.num_elements());
322  Vmath::Vcopy(inarray.num_elements(),inarray,1,tmp,1);
323  Assemble(tmp,outarray);
324  }
325  else
326  {
327  Assemble(inarray,outarray);
328  }
329 
330  GlobalSolve(key,outarray,globaltmp);
331  GlobalToLocal(globaltmp,outarray);
332  }
333  }
334 
335 
336  /**
337  * Consider the two dimensional Laplace equation,
338  * \f[\nabla\cdot\left(\boldsymbol{\sigma}\nabla
339  * u(\boldsymbol{x})\right) = f(\boldsymbol{x}),\f] supplemented with
340  * appropriate boundary conditions (which are contained in the data
341  * member #m_bndCondExpansions). In the equation above
342  * \f$\boldsymbol{\sigma}\f$ is the (symmetric positive definite)
343  * diffusion tensor:
344  * \f[ \sigma = \left[ \begin{array}{cc}
345  * \sigma_{00}(\boldsymbol{x},t) & \sigma_{01}(\boldsymbol{x},t) \\
346  * \sigma_{01}(\boldsymbol{x},t) & \sigma_{11}(\boldsymbol{x},t)
347  * \end{array} \right]. \f]
348  * Applying a \f$C^0\f$ continuous Galerkin discretisation, this
349  * equation leads to the following linear system:
350  * \f[\boldsymbol{L}
351  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f]
352  * where \f$\boldsymbol{L}\f$ is the Laplacian matrix. This function
353  * solves the system above for the global coefficients
354  * \f$\boldsymbol{\hat{u}}\f$ by a call to the function #GlobalSolve.
355  *
356  * The values of the function \f$f(\boldsymbol{x})\f$ evaluated at the
357  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
358  * variable #m_phys of the ExpList object \a Sin. The resulting global
359  * coefficients \f$\boldsymbol{\hat{u}}_g\f$ are stored in the array
360  * #m_coeffs.
361  *
362  * @param Sin An ExpList, containing the discrete evaluation
363  * of the forcing function \f$f(\boldsymbol{x})\f$
364  * at the quadrature points in its array #m_phys.
365  * @param variablecoeffs The (optional) parameter containing the
366  * coefficients evaluated at the quadrature
367  * points. It is an Array of (three) arrays which
368  * stores the laplacian coefficients in the
369  * following way
370  * \f[\mathrm{variablecoeffs} = \left[ \begin{array}{c}
371  * \left[\sigma_{00}(\boldsymbol{x_i},t)\right]_i \\
372  * \left[\sigma_{01}(\boldsymbol{x_i},t)\right]_i \\
373  * \left[\sigma_{11}(\boldsymbol{x_i},t)\right]_i
374  * \end{array}\right]
375  * \f]
376  * If this argument is not passed to the function, the following
377  * equation will be solved:
378  * \f[\nabla^2u(\boldsymbol{x}) = f(\boldsymbol{x}),\f]
379  *
380  * @param time The time-level at which the coefficients are
381  * evaluated
382  */
383  void ContField2D::LaplaceSolve(
384  const Array<OneD, const NekDouble> &inarray,
385  Array<OneD, NekDouble> &outarray,
386  const Array<OneD, const NekDouble> &dirForcing,
387  const Array<OneD, Array<OneD,NekDouble> >& variablecoeffs,
388  NekDouble time,
389  CoeffState coeffstate)
390  {
391  // Inner product of forcing
392  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
393  Array<OneD,NekDouble> wsp(contNcoeffs);
394  IProductWRTBase(inarray,wsp,eGlobal);
395  // Note -1.0 term necessary to invert forcing function to
396  // be consistent with matrix definition
397  Vmath::Neg(m_ncoeffs, wsp, 1);
398 
399  // Forcing function with weak boundary conditions
400  int i,j;
401  int bndcnt=0;
402  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
403  {
404  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
405  {
406  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
407  {
408  wsp[m_locToGloMap
409  ->GetBndCondCoeffsToGlobalCoeffsMap(bndcnt++)]
410  += (m_bndCondExpansions[i]->GetCoeffs())[j];
411  }
412  }
413  else
414  {
415  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
416  }
417  }
418 
419  StdRegions::VarCoeffMap varcoeffs;
420  varcoeffs[StdRegions::eVarCoeffD00] = variablecoeffs[0];
421  varcoeffs[StdRegions::eVarCoeffD11] = variablecoeffs[3];
422  varcoeffs[StdRegions::eVarCoeffD22] = variablecoeffs[5];
423  StdRegions::ConstFactorMap factors;
424  factors[StdRegions::eFactorTime] = time;
425 
426  // Solve the system
427  GlobalLinSysKey key(StdRegions::eLaplacian,m_locToGloMap,factors,
428  varcoeffs);
429 
430  if(coeffstate == eGlobal)
431  {
432  GlobalSolve(key,wsp,outarray,dirForcing);
433  }
434  else
435  {
436  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
437  GlobalSolve(key,wsp,tmp,dirForcing);
438  GlobalToLocal(tmp,outarray);
439  }
440  }
441 
442 
443  /**
444  * Constructs the GlobalLinearSysKey for the linear advection operator
445  * with the supplied parameters, and computes the eigenvectors and
446  * eigenvalues of the associated matrix.
447  * @param ax Advection parameter, x.
448  * @param ay Advection parameter, y.
449  * @param Real Computed eigenvalues, real component.
450  * @param Imag Computed eigenvalues, imag component.
451  * @param Evecs Computed eigenvectors.
452  */
453  void ContField2D::LinearAdvectionEigs(const NekDouble ax,
454  const NekDouble ay,
455  Array<OneD, NekDouble> &Real,
456  Array<OneD, NekDouble> &Imag,
457  Array<OneD, NekDouble> &Evecs)
458  {
459  // Solve the system
460  Array<OneD, Array<OneD, NekDouble> > vel(2);
461  Array<OneD, NekDouble> vel_x(m_npoints,ax);
462  Array<OneD, NekDouble> vel_y(m_npoints,ay);
463  vel[0] = vel_x;
464  vel[1] = vel_y;
465 
466  StdRegions::VarCoeffMap varcoeffs;
467  varcoeffs[StdRegions::eVarCoeffVelX] = Array<OneD, NekDouble>(m_npoints,ax);
468  varcoeffs[StdRegions::eVarCoeffVelY] = Array<OneD, NekDouble>(m_npoints,ay);
469  StdRegions::ConstFactorMap factors;
470  factors[StdRegions::eFactorTime] = 0.0;
471  GlobalLinSysKey key(StdRegions::eLinearAdvectionReaction,m_locToGloMap,
472  factors,varcoeffs);
473 
474  DNekMatSharedPtr Gmat = GenGlobalMatrixFull(key,m_locToGloMap);
475  Gmat->EigenSolve(Real,Imag,Evecs);
476  }
477 
478 
479 
480 
481  /**
482  * Given a linear system specified by the key \a key,
483  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}},\f]
484  * this function solves this linear system taking into account the
485  * boundary conditions specified in the data member
486  * #m_bndCondExpansions. Therefore, it adds an array
487  * \f$\boldsymbol{\hat{g}}\f$ which represents the non-zero surface
488  * integral resulting from the weak boundary conditions (e.g. Neumann
489  * boundary conditions) to the right hand side, that is,
490  * \f[\boldsymbol{M}\boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}+
491  * \boldsymbol{\hat{g}}.\f]
492  * Furthermore, it lifts the known degrees of freedom which are
493  * prescribed by the Dirichlet boundary conditions. As these known
494  * coefficients \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ are numbered
495  * first in the global coefficient array \f$\boldsymbol{\hat{u}}_g\f$,
496  * the linear system can be decomposed as,
497  * \f[\left[\begin{array}{cc}
498  * \boldsymbol{M}^{\mathcal{DD}}&\boldsymbol{M}^{\mathcal{DH}}\\
499  * \boldsymbol{M}^{\mathcal{HD}}&\boldsymbol{M}^{\mathcal{HH}}
500  * \end{array}\right]
501  * \left[\begin{array}{c}
502  * \boldsymbol{\hat{u}}^{\mathcal{D}}\\
503  * \boldsymbol{\hat{u}}^{\mathcal{H}}
504  * \end{array}\right]=
505  * \left[\begin{array}{c}
506  * \boldsymbol{\hat{f}}^{\mathcal{D}}\\
507  * \boldsymbol{\hat{f}}^{\mathcal{H}}
508  * \end{array}\right]+
509  * \left[\begin{array}{c}
510  * \boldsymbol{\hat{g}}^{\mathcal{D}}\\
511  * \boldsymbol{\hat{g}}^{\mathcal{H}}
512  * \end{array}\right]
513  * \f]
514  * which will then be solved for the unknown coefficients
515  * \f$\boldsymbol{\hat{u}}^{\mathcal{H}}\f$ as,
516  * \f[
517  * \boldsymbol{M}^{\mathcal{HH}}\boldsymbol{\hat{u}}^{\mathcal{H}}=
518  * \boldsymbol{\hat{f}}^{\mathcal{H}}+
519  * \boldsymbol{\hat{g}}^{\mathcal{H}}-
520  * \boldsymbol{M}^{\mathcal{HD}}\boldsymbol{\hat{u}}^{\mathcal{D}}\f]
521  *
522  * @param mkey This key uniquely defines the linear system to
523  * be solved.
524  * @param Sin An ExpList, containing the discrete evaluation
525  * of the forcing function \f$f(\boldsymbol{x})\f$
526  * at the quadrature points in its array #m_phys.
527  * @param ScaleForcing An optional parameter with which the forcing
528  * vector \f$\boldsymbol{\hat{f}}\f$ should be
529  * multiplied.
530  * @note inout contains initial guess and final output.
531  */
532  void ContField2D::GlobalSolve(
533  const GlobalLinSysKey &key,
534  const Array<OneD, const NekDouble>& rhs,
535  Array<OneD, NekDouble>& inout,
536  const Array<OneD, const NekDouble>& dirForcing)
537  {
538  int NumDirBcs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
539  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
540 
541  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE
542  // IN THE SOLUTION ARRAY
543  v_ImposeDirichletConditions(inout);
544 
545  // STEP 2: CALCULATE THE HOMOGENEOUS COEFFICIENTS
546  if(contNcoeffs - NumDirBcs > 0)
547  {
548  GlobalLinSysSharedPtr LinSys = GetGlobalLinSys(key);
549  LinSys->Solve(rhs,inout,m_locToGloMap,dirForcing);
550  }
551  }
552 
553 
554  /**
555  * Returns the global matrix associated with the given GlobalMatrixKey.
556  * If the global matrix has not yet been constructed on this field,
557  * it is first constructed using GenGlobalMatrix().
558  * @param mkey Global matrix key.
559  * @returns Assocated global matrix.
560  */
561  GlobalMatrixSharedPtr ContField2D::GetGlobalMatrix(
562  const GlobalMatrixKey &mkey)
563  {
564  ASSERTL1(mkey.LocToGloMapIsDefined(),
565  "To use method must have a AssemblyMap "
566  "attached to key");
567 
568  GlobalMatrixSharedPtr glo_matrix;
569  GlobalMatrixMap::iterator matrixIter = m_globalMat->find(mkey);
570 
571  if(matrixIter == m_globalMat->end())
572  {
573  glo_matrix = GenGlobalMatrix(mkey,m_locToGloMap);
574  (*m_globalMat)[mkey] = glo_matrix;
575  }
576  else
577  {
578  glo_matrix = matrixIter->second;
579  }
580 
581  return glo_matrix;
582  }
583 
584 
585  /**
586  * The function searches the map #m_globalLinSys to see if the
587  * global matrix has been created before. If not, it calls the function
588  * #GenGlobalLinSys to generate the requested global system.
589  *
590  * @param mkey This key uniquely defines the requested
591  * linear system.
592  */
593  GlobalLinSysSharedPtr ContField2D::GetGlobalLinSys(
594  const GlobalLinSysKey &mkey)
595  {
596  return m_globalLinSysManager[mkey];
597  }
598 
599  GlobalLinSysSharedPtr ContField2D::GenGlobalLinSys(
600  const GlobalLinSysKey &mkey)
601  {
602  ASSERTL1(mkey.LocToGloMapIsDefined(),
603  "To use method must have a AssemblyMap "
604  "attached to key");
605  return ExpList::GenGlobalLinSys(mkey, m_locToGloMap);
606  }
607 
608 
609  /**
610  *
611  */
612  void ContField2D::v_BwdTrans(
613  const Array<OneD, const NekDouble> &inarray,
614  Array<OneD, NekDouble> &outarray,
615  CoeffState coeffstate)
616  {
617  BwdTrans(inarray,outarray,coeffstate);
618  }
619 
620 
621  /**
622  *
623  */
624  void ContField2D::v_FwdTrans(
625  const Array<OneD, const NekDouble> &inarray,
626  Array<OneD, NekDouble> &outarray,
627  CoeffState coeffstate)
628  {
629  FwdTrans(inarray,outarray,coeffstate);
630  }
631 
632  void ContField2D::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
633  {
634  int i,j;
635  int bndcnt=0;
636  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
637 
638  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE IN THE SOLUTION
639  // ARRAY
640  NekDouble sign;
641  const Array<OneD,const int> &bndMap =
642  m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap();
643 
644  Array<OneD, NekDouble> tmp(
645  m_locToGloMap->GetNumGlobalBndCoeffs(), 0.0);
646 
647  // Fill in Dirichlet coefficients that are to be sent to
648  // other processors. This code block uses a
649  // tuple<int,int.NekDouble> which stores the local id of
650  // coefficent the global id of the data location and the
651  // inverse of the values of the data (arising from
652  // periodic boundary conditiosn)
653  map<int, vector<ExtraDirDof> > &extraDirDofs =
654  m_locToGloMap->GetExtraDirDofs();
655  map<int, vector<ExtraDirDof> >::iterator it;
656  for (it = extraDirDofs.begin(); it != extraDirDofs.end(); ++it)
657  {
658  for (i = 0; i < it->second.size(); ++i)
659  {
660  tmp[it->second.at(i).get<1>()] =
661  m_bndCondExpansions[it->first]->GetCoeffs()[
662  it->second.at(i).get<0>()]*it->second.at(i).get<2>();
663  }
664  }
665  m_locToGloMap->UniversalAssembleBnd(tmp);
666 
667  // Now fill in all other Dirichlet coefficients.
668  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
669  {
670  if(m_bndConditions[i]->GetBoundaryConditionType() ==
671  SpatialDomains::eDirichlet)
672  {
673  const Array<OneD,const NekDouble>& coeffs =
674  m_bndCondExpansions[i]->GetCoeffs();
675  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
676  {
677  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(
678  bndcnt);
679  tmp[bndMap[bndcnt++]] = sign * coeffs[j];
680  }
681  }
682  else
683  {
684  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
685  }
686  }
687 
688  Vmath::Vcopy(nDir, tmp, 1, outarray, 1);
689  }
690 
691  void ContField2D::v_FillBndCondFromField(void)
692  {
693  NekDouble sign;
694  int bndcnt = 0;
695  const Array<OneD,const int> &bndMap =
696  m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap();
697 
698  Array<OneD, NekDouble> tmp(m_locToGloMap->GetNumGlobalCoeffs());
699  LocalToGlobal(m_coeffs,tmp);
700 
701  // Now fill in all other Dirichlet coefficients.
702  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
703  {
704  Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[i]->UpdateCoeffs();
705 
706  for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
707  {
708  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
709  coeffs[j] = sign * tmp[bndMap[bndcnt++]];
710  }
711  }
712  }
713 
714  /**
715  * This operation is evaluated as:
716  * \f{tabbing}
717  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
718  * > > Do \= $i=$ $0,N_m^e-1$ \\
719  * > > > $\boldsymbol{\hat{u}}^{e}[i] = \mbox{sign}[e][i] \cdot
720  * \boldsymbol{\hat{u}}_g[\mbox{map}[e][i]]$ \\
721  * > > continue \\
722  * > continue
723  * \f}
724  * where \a map\f$[e][i]\f$ is the mapping array and \a
725  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
726  * correct modal connectivity between the different elements (both
727  * these arrays are contained in the data member #m_locToGloMap). This
728  * operation is equivalent to the scatter operation
729  * \f$\boldsymbol{\hat{u}}_l=\mathcal{A}\boldsymbol{\hat{u}}_g\f$,
730  * where \f$\mathcal{A}\f$ is the
731  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
732  *
733  * @note The array #m_coeffs should be filled with the global
734  * coefficients \f$\boldsymbol{\hat{u}}_g\f$ and that the resulting
735  * local coefficients \f$\boldsymbol{\hat{u}}_l\f$ will be stored in
736  * #m_coeffs.
737  */
738  void ContField2D::v_GlobalToLocal(void)
739  {
740  m_locToGloMap->GlobalToLocal(m_coeffs,m_coeffs);
741  }
742 
743 
744 
745  /**
746  * This operation is evaluated as:
747  * \f{tabbing}
748  * \hspace{1cm} \= Do \= $e=$ $1, N_{\mathrm{el}}$ \\
749  * > > Do \= $i=$ $0,N_m^e-1$ \\
750  * > > > $\boldsymbol{\hat{u}}_g[\mbox{map}[e][i]] =
751  * \mbox{sign}[e][i] \cdot \boldsymbol{\hat{u}}^{e}[i]$\\
752  * > > continue\\
753  * > continue
754  * \f}
755  * where \a map\f$[e][i]\f$ is the mapping array and \a
756  * sign\f$[e][i]\f$ is an array of similar dimensions ensuring the
757  * correct modal connectivity between the different elements (both
758  * these arrays are contained in the data member #m_locToGloMap). This
759  * operation is equivalent to the gather operation
760  * \f$\boldsymbol{\hat{u}}_g=\mathcal{A}^{-1}\boldsymbol{\hat{u}}_l\f$,
761  * where \f$\mathcal{A}\f$ is the
762  * \f$N_{\mathrm{eof}}\times N_{\mathrm{dof}}\f$ permutation matrix.
763  *
764  * @note The array #m_coeffs should be filled with the local
765  * coefficients \f$\boldsymbol{\hat{u}}_l\f$ and that the
766  * resulting global coefficients \f$\boldsymbol{\hat{u}}_g\f$
767  * will be stored in #m_coeffs.
768  */
769  void ContField2D::v_LocalToGlobal(void)
770  {
771  m_locToGloMap->LocalToGlobal(m_coeffs,m_coeffs);
772  }
773 
774  /**
775  *
776  */
777  void ContField2D::v_MultiplyByInvMassMatrix(
778  const Array<OneD, const NekDouble> &inarray,
779  Array<OneD, NekDouble> &outarray,
780  CoeffState coeffstate)
781  {
782  MultiplyByInvMassMatrix(inarray,outarray,coeffstate);
783  }
784 
785 
786  /**
787  * Consider the two dimensional Helmholtz equation,
788  * \f[\nabla^2u(\boldsymbol{x})-\lambda u(\boldsymbol{x})
789  * = f(\boldsymbol{x}),\f] supplemented with appropriate boundary
790  * conditions (which are contained in the data member
791  * #m_bndCondExpansions). Applying a \f$C^0\f$ continuous Galerkin
792  * discretisation, this equation leads to the following linear system:
793  * \f[\left(\boldsymbol{L}+\lambda\boldsymbol{M}\right)
794  * \boldsymbol{\hat{u}}_g=\boldsymbol{\hat{f}}\f] where
795  * \f$\boldsymbol{L}\f$ and \f$\boldsymbol{M}\f$ are the Laplacian and
796  * mass matrix respectively. This function solves the system above for
797  * the global coefficients \f$\boldsymbol{\hat{u}}\f$ by a call to the
798  * function #GlobalSolve. It is assumed #m_coeff contains an
799  * initial estimate for the solution.
800  *
801  * The values of the function \f$f(\boldsymbol{x})\f$
802  * evaluated at the quadrature points \f$\boldsymbol{x}_i\f$
803  * should be contained in the variable #m_phys of the ExpList
804  * object \a inarray. The resulting global coefficients
805  * \f$\boldsymbol{\hat{u}}_g\f$ are stored in the array
806  * #m_contCoeffs or #m_coeffs depending on whether
807  * \a coeffstate is eGlobal or eLocal
808  *
809  * @param inarray An ExpList, containing the discrete evaluation
810  * of the forcing function \f$f(\boldsymbol{x})\f$
811  * at the quadrature points in its array #m_phys.
812  * @param lambda The parameter \f$\lambda\f$ of the Helmholtz
813  * equation
814  */
815  void ContField2D::v_HelmSolve(
816  const Array<OneD, const NekDouble> &inarray,
817  Array<OneD, NekDouble> &outarray,
818  const FlagList &flags,
819  const StdRegions::ConstFactorMap &factors,
820  const StdRegions::VarCoeffMap &varcoeff,
821  const Array<OneD, const NekDouble> &dirForcing)
822  {
823  //----------------------------------
824  // Setup RHS Inner product
825  //----------------------------------
826  // Inner product of forcing
827  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
828  Array<OneD,NekDouble> wsp(contNcoeffs);
829  IProductWRTBase(inarray,wsp,eGlobal);
830  // Note -1.0 term necessary to invert forcing function to
831  // be consistent with matrix definition
832  Vmath::Neg(contNcoeffs, wsp, 1);
833 
834  // Fill weak boundary conditions
835  int i,j;
836  int bndcnt=0;
837  Array<OneD, NekDouble> gamma(contNcoeffs, 0.0);
838 
839  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
840  {
841  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
842  {
843  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
844  {
845  gamma[m_locToGloMap
846  ->GetBndCondCoeffsToGlobalCoeffsMap(bndcnt++)]
847  += (m_bndCondExpansions[i]->GetCoeffs())[j];
848  }
849  }
850  else
851  {
852  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
853  }
854  }
855 
856  m_locToGloMap->UniversalAssemble(gamma);
857 
858  // Add weak boundary conditions to forcing
859  Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
860 
861  GlobalLinSysKey key(StdRegions::eHelmholtz,m_locToGloMap,factors,varcoeff);
862 
863  if(flags.isSet(eUseGlobal))
864  {
865  GlobalSolve(key,wsp,outarray,dirForcing);
866  }
867  else
868  {
869  Array<OneD,NekDouble> tmp(contNcoeffs);
870  LocalToGlobal(outarray,tmp);
871  GlobalSolve(key,wsp,tmp,dirForcing);
872  GlobalToLocal(tmp,outarray);
873  }
874  }
875 
876 
877  /**
878  * This is equivalent to the operation:
879  * \f[\boldsymbol{M\hat{u}}_g\f]
880  * where \f$\boldsymbol{M}\f$ is the global matrix of type specified by
881  * \a mkey. After scattering the global array \a inarray to local
882  * level, this operation is evaluated locally by the function
883  * ExpList#GeneralMatrixOp. The global result is then obtained by a
884  * global assembly procedure.
885  *
886  * @param mkey This key uniquely defines the type matrix
887  * required for the operation.
888  * @param inarray The vector \f$\boldsymbol{\hat{u}}_g\f$ of size
889  * \f$N_{\mathrm{dof}}\f$.
890  * @param outarray The resulting vector of size
891  * \f$N_{\mathrm{dof}}\f$.
892  */
893  void ContField2D::v_GeneralMatrixOp(
894  const GlobalMatrixKey &gkey,
895  const Array<OneD,const NekDouble> &inarray,
896  Array<OneD, NekDouble> &outarray,
897  CoeffState coeffstate)
898  {
899  if(coeffstate == eGlobal)
900  {
901  bool doGlobalOp = m_globalOptParam->DoGlobalMatOp(
902  gkey.GetMatrixType());
903 
904  if(doGlobalOp)
905  {
906  GlobalMatrixSharedPtr mat = GetGlobalMatrix(gkey);
907  mat->Multiply(inarray,outarray);
908  m_locToGloMap->UniversalAssemble(outarray);
909  }
910  else
911  {
912  Array<OneD,NekDouble> tmp1(2*m_ncoeffs);
913  Array<OneD,NekDouble> tmp2(tmp1+m_ncoeffs);
914  GlobalToLocal(inarray,tmp1);
915  GeneralMatrixOp_IterPerExp(gkey,tmp1,tmp2);
916  Assemble(tmp2,outarray);
917  }
918  }
919  else
920  {
921  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
922  }
923  }
924 
925  /**
926  * First compute the inner product of forcing function with respect to
927  * base, and then solve the system with the linear advection operator.
928  * @param velocity Array of advection velocities in physical space
929  * @param inarray Forcing function.
930  * @param outarray Result.
931  * @param lambda reaction coefficient
932  * @param coeffstate State of Coefficients, Local or Global
933  * @param dirForcing Dirichlet Forcing.
934  */
935 
936  // could combine this with HelmholtzCG.
937  void ContField2D::v_LinearAdvectionDiffusionReactionSolve(const Array<OneD, Array<OneD, NekDouble> > &velocity,
938  const Array<OneD, const NekDouble> &inarray,
939  Array<OneD, NekDouble> &outarray,
940  const NekDouble lambda,
941  CoeffState coeffstate,
942  const Array<OneD, const NekDouble>& dirForcing)
943  {
944  // Inner product of forcing
945  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
946  Array<OneD,NekDouble> wsp(contNcoeffs);
947  IProductWRTBase(inarray,wsp,eGlobal);
948  // Note -1.0 term necessary to invert forcing function to
949  // be consistent with matrix definition
950  Vmath::Neg(contNcoeffs, wsp, 1);
951 
952  // Forcing function with weak boundary conditions
953  int i,j;
954  int bndcnt=0;
955  Array<OneD, NekDouble> gamma(contNcoeffs, 0.0);
956  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
957  {
958  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
959  {
960  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
961  {
962  gamma[m_locToGloMap
963  ->GetBndCondCoeffsToGlobalCoeffsMap(bndcnt++)]
964  += (m_bndCondExpansions[i]->GetCoeffs())[j];
965  }
966  }
967  else
968  {
969  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
970  }
971  }
972  m_locToGloMap->UniversalAssemble(wsp);
973  // Add weak boundary conditions to forcing
974  Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
975 
976  // Solve the system
977  StdRegions::ConstFactorMap factors;
978  factors[StdRegions::eFactorLambda] = lambda;
979  StdRegions::VarCoeffMap varcoeffs;
980  varcoeffs[StdRegions::eVarCoeffVelX] = velocity[0];
981  varcoeffs[StdRegions::eVarCoeffVelY] = velocity[1];
982  GlobalLinSysKey key(StdRegions::eLinearAdvectionDiffusionReaction,m_locToGloMap,factors,varcoeffs);
983 
984  if(coeffstate == eGlobal)
985  {
986  GlobalSolve(key,wsp,outarray,dirForcing);
987  }
988  else
989  {
990  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
991  GlobalSolve(key,wsp,tmp,dirForcing);
992  GlobalToLocal(tmp,outarray);
993  }
994  }
995 
996  /**
997  * First compute the inner product of forcing function with respect to
998  * base, and then solve the system with the linear advection operator.
999  * @param velocity Array of advection velocities in physical space
1000  * @param inarray Forcing function.
1001  * @param outarray Result.
1002  * @param lambda reaction coefficient
1003  * @param coeffstate State of Coefficients, Local or Global
1004  * @param dirForcing Dirichlet Forcing.
1005  */
1006  void ContField2D::v_LinearAdvectionReactionSolve(const Array<OneD, Array<OneD, NekDouble> > &velocity,
1007  const Array<OneD, const NekDouble> &inarray,
1008  Array<OneD, NekDouble> &outarray,
1009  const NekDouble lambda,
1010  CoeffState coeffstate,
1011  const Array<OneD, const NekDouble>& dirForcing)
1012  {
1013  // Inner product of forcing
1014  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
1015  Array<OneD,NekDouble> wsp(contNcoeffs);
1016  IProductWRTBase(inarray,wsp,eGlobal);
1017 
1018  // Solve the system
1019  StdRegions::ConstFactorMap factors;
1020  factors[StdRegions::eFactorLambda] = lambda;
1021  StdRegions::VarCoeffMap varcoeffs;
1022  varcoeffs[StdRegions::eVarCoeffVelX] = velocity[0];
1023  varcoeffs[StdRegions::eVarCoeffVelY] = velocity[1];
1024  GlobalLinSysKey key(StdRegions::eLinearAdvectionReaction,m_locToGloMap,factors,varcoeffs);
1025 
1026  if(coeffstate == eGlobal)
1027  {
1028  GlobalSolve(key,wsp,outarray,dirForcing);
1029  }
1030  else
1031  {
1032  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
1033  GlobalSolve(key,wsp,tmp,dirForcing);
1034  GlobalToLocal(tmp,outarray);
1035  }
1036  }
1037 
1038 
1039  /**
1040  *
1041  */
1042  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>&
1043  ContField2D::v_GetBndConditions()
1044  {
1045  return GetBndConditions();
1046  }
1047 
1048  } // end of namespace
1049 } //end of namespace
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2