Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StdNodalTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdNodalTriExp.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: Nodal triangle routines built upon StdExpansion2D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <StdRegions/StdNodalTriExp.h>
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 
39 namespace Nektar
40 {
41  namespace StdRegions
42  {
43  StdNodalTriExp::StdNodalTriExp():
44  StdTriExp(),
45  m_nodalPointsKey()
46  {
47  }
48 
49  StdNodalTriExp::StdNodalTriExp(
50  const LibUtilities::BasisKey &Ba,
51  const LibUtilities::BasisKey &Bb,
52  LibUtilities::PointsType Ntype):
53  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(
54  Ba.GetNumModes(),
55  Bb.GetNumModes()),
56  2,Ba,Bb),
57  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
58  Ba.GetNumModes(),
59  Bb.GetNumModes()),
60  Ba,Bb),
61  StdTriExp (Ba,Bb),
62  m_nodalPointsKey()
63  {
64  ASSERTL0(m_base[0]->GetNumModes() == m_base[1]->GetNumModes(),
65  "Nodal basis initiated with different orders in the a "
66  "and b directions");
67  int nummodes = Ba.GetNumModes();
68  m_nodalPointsKey = MemoryManager<LibUtilities::PointsKey>::
69  AllocateSharedPtr(nummodes,Ntype);
70  }
71 
72  StdNodalTriExp::StdNodalTriExp(const StdNodalTriExp &T):
73  StdExpansion(T),
74  StdExpansion2D(T),
75  StdTriExp(T),
76  m_nodalPointsKey(T.m_nodalPointsKey)
77  {
78  }
79 
80  StdNodalTriExp::~StdNodalTriExp()
81  {
82  }
83 
84 
85  //-------------------------------
86  // Nodal basis specific routines
87  //-------------------------------
88 
89  void StdNodalTriExp::NodalToModal(
90  const Array<OneD, const NekDouble>& inarray,
91  Array<OneD, NekDouble>& outarray)
92  {
93  StdMatrixKey Nkey(eInvNBasisTrans, DetShapeType(), *this,
94  NullConstFactorMap, NullVarCoeffMap,
95  m_nodalPointsKey->GetPointsType());
96  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
97 
98  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eWrapper);
99  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
100  modal = (*inv_vdm) * nodal;
101  }
102 
103  // Operate with transpose of NodalToModal transformation
104  void StdNodalTriExp::NodalToModalTranspose(
105  const Array<OneD, const NekDouble>& inarray,
106  Array<OneD, NekDouble>& outarray)
107  {
108  StdMatrixKey Nkey(eInvNBasisTrans, DetShapeType(), *this,
109  NullConstFactorMap, NullVarCoeffMap,
110  m_nodalPointsKey->GetPointsType());
111  DNekMatSharedPtr inv_vdm = GetStdMatrix(Nkey);
112 
113  NekVector<NekDouble> nodal(m_ncoeffs,inarray,eCopy);
114  NekVector<NekDouble> modal(m_ncoeffs,outarray,eWrapper);
115  modal = Transpose(*inv_vdm) * nodal;
116  }
117 
118  void StdNodalTriExp::ModalToNodal(
119  const Array<OneD, const NekDouble>& inarray,
120  Array<OneD, NekDouble>& outarray)
121  {
122  StdMatrixKey Nkey(eNBasisTrans, DetShapeType(), *this,
123  NullConstFactorMap, NullVarCoeffMap,
124  m_nodalPointsKey->GetPointsType());
125  DNekMatSharedPtr vdm = GetStdMatrix(Nkey);
126 
127  // Multiply out matrix
128  NekVector<NekDouble> modal(m_ncoeffs,inarray,eWrapper);
129  NekVector<NekDouble> nodal(m_ncoeffs,outarray,eWrapper);
130  nodal = (*vdm)*modal;
131  }
132 
133  void StdNodalTriExp::GetNodalPoints(
134  Array<OneD, const NekDouble> &x,
135  Array<OneD, const NekDouble> &y)
136  {
137  LibUtilities::PointsManager()[*m_nodalPointsKey]->GetPoints(x,y);
138  }
139 
140  DNekMatSharedPtr StdNodalTriExp::GenNBasisTransMatrix()
141  {
142  int i,j;
143  Array<OneD, const NekDouble> r, s;
144  Array<OneD, NekDouble> c(2);
145  DNekMatSharedPtr Mat;
146 
147  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(
148  m_ncoeffs, m_ncoeffs);
149  GetNodalPoints(r,s);
150 
151  //Store the values of m_phys in a temporary array
152  int nqtot = GetTotPoints();
153  Array<OneD,NekDouble> phys(nqtot);
154 
155  for(i = 0; i < m_ncoeffs; ++i)
156  {
157  // fill physical space with mode i
158  StdTriExp::v_FillMode(i,phys);
159 
160  // interpolate mode i to the Nodal points 'j' and
161  // store in outarray
162  for(j = 0; j < m_ncoeffs; ++j)
163  {
164  c[0] = r[j];
165  c[1] = s[j];
166  (*Mat)(j,i) = StdTriExp::v_PhysEvaluate(c,phys);
167  }
168  }
169  return Mat;
170  }
171 
172 
173  //---------------------------------------
174  // Transforms
175  //---------------------------------------
176 
177  void StdNodalTriExp::v_BwdTrans(
178  const Array<OneD, const NekDouble>& inarray,
179  Array<OneD, NekDouble>& outarray)
180  {
181  v_BwdTrans_SumFac(inarray,outarray);
182  }
183 
184  void StdNodalTriExp::v_BwdTrans_SumFac(
185  const Array<OneD, const NekDouble>& inarray,
186  Array<OneD, NekDouble>& outarray)
187  {
188  Array<OneD, NekDouble> tmp(m_ncoeffs);
189  NodalToModal(inarray,tmp);
190  StdTriExp::v_BwdTrans_SumFac(tmp,outarray);
191  }
192 
193  void StdNodalTriExp::v_FwdTrans(
194  const Array<OneD, const NekDouble>& inarray,
195  Array<OneD, NekDouble>& outarray)
196  {
197  v_IProductWRTBase(inarray,outarray);
198 
199  // get Mass matrix inverse
200  StdMatrixKey masskey(eInvMass, DetShapeType(), *this,
201  NullConstFactorMap, NullVarCoeffMap,
202  m_nodalPointsKey->GetPointsType());
203  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
204 
205  // copy inarray in case inarray == outarray
206  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
207  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
208 
209  out = (*matsys)*in;
210  }
211 
212 
213  //---------------------------------------
214  // Inner product functions
215  //---------------------------------------
216 
217  void StdNodalTriExp::v_IProductWRTBase(
218  const Array<OneD, const NekDouble>& inarray,
219  Array<OneD, NekDouble>& outarray)
220  {
221  v_IProductWRTBase_SumFac(inarray,outarray);
222  }
223 
224  void StdNodalTriExp::v_IProductWRTBase_SumFac(
225  const Array<OneD, const NekDouble>& inarray,
226  Array<OneD, NekDouble>& outarray)
227  {
228  StdTriExp::v_IProductWRTBase_SumFac(inarray,outarray);
229  NodalToModalTranspose(outarray,outarray);
230  }
231 
232  void StdNodalTriExp::v_IProductWRTDerivBase(
233  const int dir,
234  const Array<OneD, const NekDouble>& inarray,
235  Array<OneD, NekDouble>& outarray)
236  {
237  v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
238  }
239 
240  void StdNodalTriExp::v_IProductWRTDerivBase_SumFac(
241  const int dir,
242  const Array<OneD, const NekDouble>& inarray,
243  Array<OneD, NekDouble>& outarray)
244  {
245  StdTriExp::v_IProductWRTDerivBase_SumFac(dir,inarray,outarray);
246  NodalToModalTranspose(outarray,outarray);
247  }
248 
249  //---------------------------------------
250  // Evaluation functions
251  //---------------------------------------
252 
253  void StdNodalTriExp::v_FillMode(
254  const int mode,
255  Array<OneD, NekDouble> &outarray)
256  {
257  ASSERTL2(mode >= m_ncoeffs,
258  "calling argument mode is larger than total expansion order");
259 
260  Vmath::Zero(m_ncoeffs, outarray, 1);
261  outarray[mode] = 1.0;
262  v_BwdTrans(outarray,outarray);
263  }
264 
265 
266  //---------------------------
267  // Helper functions
268  //---------------------------
269 
270  int StdNodalTriExp::v_NumBndryCoeffs() const
271  {
272  return 3 + (GetBasisNumModes(0)-2) + 2*(GetBasisNumModes(1)-2);
273  }
274 
275  //--------------------------
276  // Mappings
277  //--------------------------
278 
279  void StdNodalTriExp::v_GetEdgeToElementMap(
280  const int eid,
281  const Orientation edgeOrient,
282  Array<OneD, unsigned int> &maparray,
283  Array<OneD, int> &signarray)
284  {
285  ASSERTL0(eid >= 0 && eid <= 2,
286  "Local Edge ID must be between 0 and 2");
287 
288  const int nEdgeCoeffs = GetEdgeNcoeffs(eid);
289 
290  if (maparray.num_elements() != nEdgeCoeffs)
291  {
292  maparray = Array<OneD, unsigned int>(nEdgeCoeffs);
293  }
294 
295  if (signarray.num_elements() != nEdgeCoeffs)
296  {
297  signarray = Array<OneD, int>(nEdgeCoeffs,1);
298  }
299  else
300  {
301  fill(signarray.get(), signarray.get()+nEdgeCoeffs, 1);
302  }
303 
304  maparray[0] = eid;
305  maparray[1] = eid == 2 ? 0 : eid+1;
306  for (int i = 2; i < nEdgeCoeffs; i++)
307  {
308  maparray[i] = eid*(nEdgeCoeffs-2)+1+i;
309  }
310 
311  if (edgeOrient == eBackwards)
312  {
313  reverse(maparray.get(), maparray.get()+nEdgeCoeffs);
314  }
315  }
316 
317  int StdNodalTriExp::v_GetVertexMap(const int localVertexId,
318  bool useCoeffPacking)
319  {
320  ASSERTL0(localVertexId >= 0 && localVertexId <= 2,
321  "Local Vertex ID must be between 0 and 2");
322  return localVertexId;
323  }
324 
325  void StdNodalTriExp::v_GetEdgeInteriorMap(
326  const int eid,
327  const Orientation edgeOrient,
328  Array<OneD, unsigned int> &maparray,
329  Array<OneD, int> &signarray)
330  {
331  ASSERTL0(eid >= 0 && eid <= 2,
332  "Local Edge ID must be between 0 and 2");
333 
334  const int nEdgeIntCoeffs = GetEdgeNcoeffs(eid)-2;
335 
336  if (maparray.num_elements() != nEdgeIntCoeffs)
337  {
338  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
339  }
340 
341  if (signarray.num_elements() != nEdgeIntCoeffs)
342  {
343  signarray = Array<OneD, int>(nEdgeIntCoeffs,1);
344  }
345  else
346  {
347  fill(signarray.get(), signarray.get()+nEdgeIntCoeffs, 1);
348  }
349 
350  for (int i = 0; i < nEdgeIntCoeffs; i++)
351  {
352  maparray[i] = eid*nEdgeIntCoeffs+3+i;
353  }
354 
355  if (edgeOrient == eBackwards)
356  {
357  reverse(maparray.get(), maparray.get()+nEdgeIntCoeffs);
358  }
359  }
360 
361  void StdNodalTriExp::v_GetInteriorMap(
362  Array<OneD, unsigned int>& outarray)
363  {
364  unsigned int i;
365  if (outarray.num_elements() != GetNcoeffs()-NumBndryCoeffs())
366  {
367  outarray = Array<OneD, unsigned int>(
368  GetNcoeffs()-NumBndryCoeffs());
369  }
370 
371  for (i = NumBndryCoeffs(); i < GetNcoeffs(); i++)
372  {
373  outarray[i-NumBndryCoeffs()] = i;
374  }
375  }
376 
377  void StdNodalTriExp::v_GetBoundaryMap(
378  Array<OneD, unsigned int>& outarray)
379  {
380  unsigned int i;
381  if (outarray.num_elements()!=NumBndryCoeffs())
382  {
383  outarray = Array<OneD, unsigned int>(NumBndryCoeffs());
384  }
385 
386  for (i = 0; i < NumBndryCoeffs(); i++)
387  {
388  outarray[i] = i;
389  }
390  }
391 
392 
393  //---------------------------------------
394  // Wrapper functions
395  //---------------------------------------
396 
397  DNekMatSharedPtr StdNodalTriExp::v_GenMatrix(const StdMatrixKey &mkey)
398  {
399  DNekMatSharedPtr Mat;
400 
401  switch(mkey.GetMatrixType())
402  {
403  case eNBasisTrans:
404  Mat = GenNBasisTransMatrix();
405  break;
406  default:
407  Mat = StdExpansion::CreateGeneralMatrix(mkey);
408  break;
409  }
410 
411  return Mat;
412  }
413 
414  DNekMatSharedPtr StdNodalTriExp::v_CreateStdMatrix(
415  const StdMatrixKey &mkey)
416  {
417  return StdNodalTriExp::v_GenMatrix(mkey);
418  }
419 
420 
421  //---------------------------------------
422  // Operator evaluation functions
423  //---------------------------------------
424 
425  void StdNodalTriExp::v_MassMatrixOp(
426  const Array<OneD, const NekDouble> &inarray,
427  Array<OneD, NekDouble> &outarray,
428  const StdMatrixKey &mkey)
429  {
430  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
431  }
432 
433  void StdNodalTriExp::v_LaplacianMatrixOp(
434  const Array<OneD, const NekDouble> &inarray,
435  Array<OneD, NekDouble> &outarray,
436  const StdMatrixKey &mkey)
437  {
438  StdExpansion::LaplacianMatrixOp_MatFree_GenericImpl(
439  inarray,outarray,mkey);
440  }
441 
442  void StdNodalTriExp::v_LaplacianMatrixOp(
443  const int k1,
444  const int k2,
445  const Array<OneD, const NekDouble> &inarray,
446  Array<OneD, NekDouble> &outarray,
447  const StdMatrixKey &mkey)
448 
449  {
450  StdExpansion::LaplacianMatrixOp_MatFree(
451  k1,k2,inarray,outarray,mkey);
452  }
453 
454  void StdNodalTriExp::v_WeakDerivMatrixOp(
455  const int i,
456  const Array<OneD, const NekDouble> &inarray,
457  Array<OneD, NekDouble> &outarray,
458  const StdMatrixKey &mkey)
459  {
460  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
461  }
462 
463  void StdNodalTriExp::v_HelmholtzMatrixOp(
464  const Array<OneD, const NekDouble> &inarray,
465  Array<OneD, NekDouble> &outarray,
466  const StdMatrixKey &mkey)
467  {
468  StdExpansion::HelmholtzMatrixOp_MatFree_GenericImpl(
469  inarray,outarray,mkey);
470  }
471 
472 
473  //---------------------------------------
474  // Private helper functions
475  //---------------------------------------
476 
477  } // end of namespace
478 } // end of namespace
479 
Generated on Sun Mar 15 2015 00:11:35 for Nektar++ by   doxygen 1.8.1.2