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