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