Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ContField3D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ContField3D.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 3D domain with boundary conditions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44 
47  m_locToGloMap(),
48  m_globalMat(),
49  m_globalLinSysManager(
50  boost::bind(&ContField3D::GenGlobalLinSys, this, _1),
51  std::string("GlobalLinSys"))
52  {
53  }
54 
55 
56  /**
57  * Given a mesh \a graph2D, containing information about the domain and
58  * the spectral/hp element expansion, this constructor fills the list
59  * of local expansions #m_exp with the proper expansions, calculates
60  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
61  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory
62  * for the arrays #m_coeffs and #m_phys. Furthermore, it constructs the
63  * mapping array (contained in #m_locToGloMap) for the transformation
64  * between local elemental level and global level, it calculates the
65  * total number global expansion coefficients \f$\hat{u}_n\f$ and
66  * allocates memory for the array #m_coeffs. The constructor also
67  * discretises the boundary conditions, specified by the argument \a
68  * bcs, by expressing them in terms of the coefficient of the expansion
69  * on the boundary.
70  *
71  * @param pSession Session information.
72  * @param graph3D A mesh, containing information about the domain
73  * and the spectral/hp element expansion.
74  * @param variable The variable associated with this field.
75  */
78  const std::string &variable,
79  const bool CheckIfSingularSystem):
80  DisContField3D(pSession,graph3D,variable,false),
81  m_globalMat(MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
82  m_globalLinSysManager(
83  boost::bind(&ContField3D::GenGlobalLinSys, this, _1),
84  std::string("GlobalLinSys"))
85  {
89  CheckIfSingularSystem, variable);
90 
91  if (m_session->DefinesCmdLineArgument("verbose"))
92  {
93  m_locToGloMap->PrintStats(std::cout, variable);
94  }
95  }
96 
97 
98  /**
99  * Given a mesh \a graph3D, containing information about the domain and
100  * the spectral/hp element expansion, this constructor fills the list
101  * of local expansions #m_exp with the proper expansions, calculates
102  * the total number of quadrature points \f$\boldsymbol{x}_i\f$ and
103  * local expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory
104  * for the arrays #m_coeffs and #m_phys. Furthermore, it constructs the
105  * mapping array (contained in #m_locToGloMap) for the transformation
106  * between local elemental level and global level, it calculates the
107  * total number global expansion coefficients \f$\hat{u}_n\f$ and
108  * allocates memory for the array #m_coeffs. The constructor also
109  * discretises the boundary conditions, specified by the argument \a
110  * bcs, by expressing them in terms of the coefficient of the expansion
111  * on the boundary.
112  *
113  * @param In Existing ContField2D object used to provide the
114  * local to global mapping information and
115  * global solution type.
116  * @param graph3D A mesh, containing information about the domain
117  * and the spectral/hp element expansion.
118  * @param bcs The boundary conditions.
119  * @param bc_loc
120  */
122  const SpatialDomains::MeshGraphSharedPtr &graph3D,
123  const std::string &variable,
124  const bool CheckIfSingularSystem):
125  DisContField3D(In,graph3D,variable,false),
126  m_globalMat (MemoryManager<GlobalMatrixMap>::AllocateSharedPtr()),
127  m_globalLinSysManager(boost::bind(&ContField3D::GenGlobalLinSys, this, _1),
128  std::string("GlobalLinSys"))
129 
130  {
131  if(!SameTypeOfBoundaryConditions(In) || CheckIfSingularSystem)
132  {
137  CheckIfSingularSystem,variable);
138 
139  if (m_session->DefinesCmdLineArgument("verbose"))
140  {
141  m_locToGloMap->PrintStats(std::cout, variable);
142  }
143  }
144  else
145  {
147  }
148  }
149 
150 
152  DisContField3D(In),
153  m_locToGloMap(In.m_locToGloMap),
154  m_globalMat(In.m_globalMat),
155  m_globalLinSysManager(In.m_globalLinSysManager)
156  {
157  }
158 
159 
161  {
162  }
163 
164 
165  /**
166  * Given the coefficients of an expansion, this function evaluates the
167  * spectral/hp expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
168  * quadrature points \f$\boldsymbol{x}_i\f$. This operation is
169  * evaluated locally by the function ExpList#BwdTrans.
170  *
171  * The coefficients of the expansion should be contained in the variable
172  * #inarray of the ExpList object \a In. The resulting physical values
173  * at the quadrature points \f$u^{\delta}(\boldsymbol{x}_i)\f$ are
174  * stored in the array #outarray.
175  *
176  * @param In An ExpList, containing the local coefficients
177  * \f$\hat{u}_n^e\f$ in its array #m_coeffs.
178  */
180  const Array<OneD, const NekDouble> &inarray,
181  Array<OneD, NekDouble> &outarray,
182  CoeffState coeffstate)
183  {
184  if(coeffstate == eGlobal)
185  {
186  bool doGlobalOp = m_globalOptParam->DoGlobalMatOp(
188 
189  if(doGlobalOp)
190  {
193  mat->Multiply(inarray,outarray);
194  }
195  else
196  {
197  Array<OneD, NekDouble> wsp(m_ncoeffs);
198  GlobalToLocal(inarray,wsp);
199  BwdTrans_IterPerExp(wsp,outarray);
200  }
201  }
202  else
203  {
204  BwdTrans_IterPerExp(inarray,outarray);
205  }
206  }
207 
208 
209  /**
210  * The operation is evaluated locally (i.e. with respect to all local
211  * expansion modes) by the function ExpList#IProductWRTBase. The inner
212  * product with respect to the global expansion modes is than obtained
213  * by a global assembly operation.
214  *
215  * The values of the function \f$f(\boldsymbol{x})\f$ evaluated at the
216  * quadrature points \f$\boldsymbol{x}_i\f$ should be contained in the
217  * variable #m_phys of the ExpList object \a in. The result is stored
218  * in the array #m_coeffs.
219  *
220  * @param In An ExpList, containing the discrete evaluation
221  * of \f$f(\boldsymbol{x})\f$ at the quadrature
222  * points in its array #m_phys.
223  */
224  void ContField3D::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
225  Array<OneD, NekDouble> &outarray,
226  CoeffState coeffstate)
227  {
228  if(coeffstate == eGlobal)
229  {
230  bool doGlobalOp = m_globalOptParam->DoGlobalMatOp(
232 
233  if(doGlobalOp)
234  {
236  m_locToGloMap);
238  mat->Multiply(inarray,outarray);
239  m_locToGloMap->UniversalAssemble(outarray);
240  }
241  else
242  {
243  Array<OneD, NekDouble> wsp(m_ncoeffs);
244  IProductWRTBase_IterPerExp(inarray,wsp);
245  Assemble(wsp,outarray);
246  }
247  }
248  else
249  {
250  IProductWRTBase_IterPerExp(inarray,outarray);
251  }
252  }
253 
254 
255  void ContField3D::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
256  Array<OneD, NekDouble> &outarray,
257  CoeffState coeffstate)
258  {
259  // Inner product of forcing
260  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
261  Array<OneD,NekDouble> wsp(contNcoeffs);
262  IProductWRTBase(inarray,wsp,eGlobal);
263 
264  // Solve the system
266 
267  if(coeffstate == eGlobal)
268  {
269  GlobalSolve(key,wsp,outarray);
270  }
271  else
272  {
273  Array<OneD,NekDouble> tmp(contNcoeffs,0.0);
274  GlobalSolve(key,wsp,tmp);
275  GlobalToLocal(tmp,outarray);
276  }
277  }
278 
279 
281  const Array<OneD, const NekDouble> &inarray,
282  Array<OneD, NekDouble> &outarray,
283  CoeffState coeffstate)
284 
285  {
286  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
288 
289  if(coeffstate == eGlobal)
290  {
291  if(inarray.data() == outarray.data())
292  {
293  Array<OneD, NekDouble> tmp(contNcoeffs,0.0);
294  Vmath::Vcopy(contNcoeffs,inarray,1,tmp,1);
295  GlobalSolve(key,tmp,outarray);
296  }
297  else
298  {
299  GlobalSolve(key,inarray,outarray);
300  }
301  }
302  else
303  {
304  Array<OneD, NekDouble> globaltmp(contNcoeffs,0.0);
305 
306  if(inarray.data() == outarray.data())
307  {
308  Array<OneD,NekDouble> tmp(inarray.num_elements());
309  Vmath::Vcopy(inarray.num_elements(),inarray,1,tmp,1);
310  Assemble(tmp,outarray);
311  }
312  else
313  {
314  Assemble(inarray,outarray);
315  }
316 
317  GlobalSolve(key,outarray,globaltmp);
318  GlobalToLocal(globaltmp,outarray);
319  }
320  }
321 
322 
324  Array<OneD, NekDouble> &inout,
325  Array<OneD, NekDouble> &outarray)
326  {
327  int bndcnt=0;
328  const Array<OneD,const int>& map = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap();
329  NekDouble sign;
330 
331  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
332  {
333  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
334  {
335  const Array<OneD,const NekDouble>& coeffs = m_bndCondExpansions[i]->GetCoeffs();
336  for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
337  {
338  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
339  inout[map[bndcnt++]] = sign * coeffs[j];
340  }
341  }
342  else
343  {
344  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
345  }
346  }
347  GeneralMatrixOp(key,inout,outarray,eGlobal);
348  }
349 
350 
351 
352  // Note inout contains initial guess and final output.
354  const GlobalLinSysKey &key,
355  const Array<OneD, const NekDouble> &rhs,
356  Array<OneD, NekDouble> &inout,
357  const Array<OneD, const NekDouble> &dirForcing)
358  {
359  int NumDirBcs = m_locToGloMap->GetNumGlobalDirBndCoeffs();
360  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
361 
362  // STEP 1: SET THE DIRICHLET DOFS TO THE RIGHT VALUE
363  // IN THE SOLUTION ARRAY
365 
366 
367  // STEP 2: CALCULATE THE HOMOGENEOUS COEFFICIENTS
368  if(contNcoeffs - NumDirBcs > 0)
369  {
371  LinSys->Solve(rhs,inout,m_locToGloMap,dirForcing);
372  }
373  }
374 
376  {
377  return m_globalLinSysManager[mkey];
378  }
379 
380 
382  {
384  "To use method must have a AssemblyMap "
385  "attached to key");
387  }
388 
389 
390  /**
391  * Returns the global matrix associated with the given GlobalMatrixKey.
392  * If the global matrix has not yet been constructed on this field,
393  * it is first constructed using GenGlobalMatrix().
394  * @param mkey Global matrix key.
395  * @returns Assocated global matrix.
396  */
398  {
400  "To use method must have a AssemblyMap "
401  "attached to key");
402 
403  GlobalMatrixSharedPtr glo_matrix;
404  GlobalMatrixMap::iterator matrixIter = m_globalMat->find(mkey);
405 
406  if(matrixIter == m_globalMat->end())
407  {
408  glo_matrix = GenGlobalMatrix(mkey,m_locToGloMap);
409  (*m_globalMat)[mkey] = glo_matrix;
410  }
411  else
412  {
413  glo_matrix = matrixIter->second;
414  }
415 
416  return glo_matrix;
417  }
418 
419 
420  void ContField3D::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
421  {
422  int i,j;
423  int bndcnt = 0;
424  int nDir = m_locToGloMap->GetNumGlobalDirBndCoeffs();
425 
426  NekDouble sign;
427  const Array<OneD,const int> &bndMap =
428  m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap();
429 
430  Array<OneD, NekDouble> tmp(
431  m_locToGloMap->GetNumGlobalBndCoeffs(), 0.0);
432 
433  // Fill in Dirichlet coefficients that are to be sent to other
434  // processors. This code block uses a tuple<int,int.NekDouble> which
435  // stores the local id of coefficent the global id of the data
436  // location and the inverse of the values of the data (arising from
437  // periodic boundary conditiosn)
438  map<int, vector<ExtraDirDof> > &extraDirDofs =
439  m_locToGloMap->GetExtraDirDofs();
440  map<int, vector<ExtraDirDof> >::iterator it;
441  for (it = extraDirDofs.begin(); it != extraDirDofs.end(); ++it)
442  {
443  for (i = 0; i < it->second.size(); ++i)
444  {
445  tmp[it->second.at(i).get<1>()] =
446  m_bndCondExpansions[it->first]->GetCoeffs()[
447  it->second.at(i).get<0>()]*it->second.at(i).get<2>();
448  }
449  }
450 
451  m_locToGloMap->UniversalAssembleBnd(tmp);
452 
453  // Now fill in all other Dirichlet coefficients.
454  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
455  {
456  if(m_bndConditions[i]->GetBoundaryConditionType() ==
458  {
459  const Array<OneD,const NekDouble>& coeffs =
460  m_bndCondExpansions[i]->GetCoeffs();
461  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
462  {
463  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(
464  bndcnt);
465  tmp[bndMap[bndcnt++]] = sign * coeffs[j];
466  }
467  }
468  else
469  {
470  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
471  }
472  }
473 
474  Vmath::Vcopy(nDir, tmp, 1, outarray, 1);
475  }
476 
478  {
479  NekDouble sign;
480  int bndcnt = 0;
481  const Array<OneD,const int> &bndMap =
482  m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap();
483 
484  Array<OneD, NekDouble> tmp(m_locToGloMap->GetNumGlobalCoeffs());
485  LocalToGlobal(m_coeffs,tmp);
486 
487  // Now fill in all other Dirichlet coefficients.
488  for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
489  {
490  Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[i]->UpdateCoeffs();
491 
492  for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
493  {
494  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
495  coeffs[j] = sign * tmp[bndMap[bndcnt++]];
496  }
497  }
498  }
499 
501  {
502  m_locToGloMap->LocalToGlobal(m_coeffs, m_coeffs);
503  }
504 
506  {
507  m_locToGloMap->GlobalToLocal(m_coeffs, m_coeffs);
508  }
509 
510 
511 
513  const Array<OneD, const NekDouble> &inarray,
514  Array<OneD, NekDouble> &outarray,
515  const FlagList &flags,
516  const StdRegions::ConstFactorMap &factors,
517  const StdRegions::VarCoeffMap &varcoeff,
518  const Array<OneD, const NekDouble> &dirForcing)
519  {
520  // Inner product of forcing
521  int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
522  Array<OneD,NekDouble> wsp(contNcoeffs);
523  IProductWRTBase(inarray,wsp,eGlobal);
524 
525  // Note -1.0 term necessary to invert forcing function to
526  // be consistent with matrix definition
527  Vmath::Neg(contNcoeffs, wsp, 1);
528 
529  // Forcing function with weak boundary conditions
530  int i,j;
531  int bndcnt = 0;
532  NekDouble sign;
533  Array<OneD, NekDouble> gamma(contNcoeffs, 0.0);
534  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
535  {
536  if(m_bndConditions[i]->GetBoundaryConditionType() != SpatialDomains::eDirichlet)
537  {
538  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
539  {
540  sign = m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
541  gamma[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(bndcnt++)] +=
542  sign * (m_bndCondExpansions[i]->GetCoeffs())[j];
543  }
544  }
545  else
546  {
547  bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
548  }
549  }
550  m_locToGloMap->UniversalAssemble(gamma);
551 
552  // Add weak boundary conditions to forcing
553  Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
554 
555  // Solve the system
556  GlobalLinSysKey key(StdRegions::eHelmholtz, m_locToGloMap, factors,varcoeff);
557 
558  if(flags.isSet(eUseGlobal))
559  {
560  GlobalSolve(key,wsp,outarray,dirForcing);
561  }
562  else
563  {
564  Array<OneD,NekDouble> tmp(contNcoeffs);
565  LocalToGlobal(outarray,tmp);
566  GlobalSolve(key,wsp,tmp,dirForcing);
567  GlobalToLocal(tmp,outarray);
568  }
569  }
570 
572  const GlobalMatrixKey &gkey,
573  const Array<OneD,const NekDouble> &inarray,
574  Array<OneD, NekDouble> &outarray,
575  CoeffState coeffstate)
576  {
577  if(coeffstate == eGlobal)
578  {
579  bool doGlobalOp = m_globalOptParam->DoGlobalMatOp(gkey.GetMatrixType());
580 
581  if(doGlobalOp)
582  {
584  mat->Multiply(inarray,outarray);
585  m_locToGloMap->UniversalAssemble(outarray);
586  }
587  else
588  {
589  Array<OneD,NekDouble> tmp1(2*m_ncoeffs);
590  Array<OneD,NekDouble> tmp2(tmp1+m_ncoeffs);
591  GlobalToLocal(inarray,tmp1);
592  GeneralMatrixOp_IterPerExp(gkey,tmp1,tmp2);
593  Assemble(tmp2,outarray);
594  }
595  }
596  else
597  {
598  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
599  }
600  }
601 
603  {
605  "To use method must have a AssemblyMap "
606  "attached to key");
607 
608  GlobalMatrixMap::iterator matrixIter = m_globalMat->find(gkey);
609 
610  if(matrixIter == m_globalMat->end())
611  {
612  return 0;
613  }
614  else
615  {
616  return matrixIter->second->GetNumNonZeroEntries();
617  }
618 
619  return 0;
620  }
621 
622  } //end of namespace
623 } //end of namespace