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