Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TriExp.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 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Expasion for triangular elements.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 #include <LibUtilities/Foundations/InterpCoeff.h>
35 #include <LocalRegions/TriExp.h>
36 #include <LocalRegions/SegExp.h>
37 #include <LocalRegions/Expansion3D.h>
38 #include <StdRegions/StdNodalTriExp.h>
39 #include <LibUtilities/Foundations/Interp.h>
40 
41 
42 namespace Nektar
43 {
44  namespace LocalRegions
45  {
46  TriExp::TriExp(const LibUtilities::BasisKey &Ba,
47  const LibUtilities::BasisKey &Bb,
48  const SpatialDomains::TriGeomSharedPtr &geom):
49  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),2,Ba,Bb),
50  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),Ba,Bb),
51  StdTriExp(Ba,Bb),
52  Expansion (geom),
53  Expansion2D (geom),
54  m_matrixManager(
55  boost::bind(&TriExp::CreateMatrix, this, _1),
56  std::string("TriExpMatrix")),
57  m_staticCondMatrixManager(
58  boost::bind(&TriExp::CreateStaticCondMatrix, this, _1),
59  std::string("TriExpStaticCondMatrix"))
60  {
61  }
62 
63 
64  TriExp::TriExp(const TriExp &T):
65  StdExpansion(T),
66  StdExpansion2D(T),
67  StdTriExp(T),
68  Expansion(T),
69  Expansion2D(T),
70  m_matrixManager(T.m_matrixManager),
71  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
72  {
73  }
74 
75 
76  TriExp::~TriExp()
77  {
78  }
79 
80 
81 
82  NekDouble TriExp::v_Integral(const Array<OneD, const NekDouble> &inarray)
83  {
84  int nquad0 = m_base[0]->GetNumPoints();
85  int nquad1 = m_base[1]->GetNumPoints();
86  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
87  NekDouble ival;
88  Array<OneD,NekDouble> tmp(nquad0*nquad1);
89 
90  // multiply inarray with Jacobian
91  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
92  {
93  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
94  }
95  else
96  {
97  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
98  }
99 
100  // call StdQuadExp version;
101  ival = StdTriExp::v_Integral(tmp);
102  return ival;
103  }
104 
105 
106  void TriExp::v_PhysDeriv(const Array<OneD, const NekDouble> & inarray,
107  Array<OneD,NekDouble> &out_d0,
108  Array<OneD,NekDouble> &out_d1,
109  Array<OneD,NekDouble> &out_d2)
110  {
111  int nquad0 = m_base[0]->GetNumPoints();
112  int nquad1 = m_base[1]->GetNumPoints();
113  int nqtot = nquad0*nquad1;
114  const Array<TwoD, const NekDouble>& df
115  = m_metricinfo->GetDerivFactors(GetPointsKeys());
116 
117  Array<OneD,NekDouble> diff0(2*nqtot);
118  Array<OneD,NekDouble> diff1(diff0+nqtot);
119 
120  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
121 
122  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
123  {
124  if(out_d0.num_elements())
125  {
126  Vmath::Vmul (nqtot,df[0],1,diff0,1, out_d0, 1);
127  Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
128  }
129 
130  if(out_d1.num_elements())
131  {
132  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
133  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
134  }
135 
136  if(out_d2.num_elements())
137  {
138  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
139  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
140  }
141  }
142  else // regular geometry
143  {
144  if(out_d0.num_elements())
145  {
146  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
147  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
148  }
149 
150  if(out_d1.num_elements())
151  {
152  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
153  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
154  }
155 
156  if(out_d2.num_elements())
157  {
158  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
159  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
160  }
161  }
162  }
163 
164 
165  void TriExp::v_PhysDeriv(const int dir,
166  const Array<OneD, const NekDouble>& inarray,
167  Array<OneD, NekDouble> &outarray)
168  {
169  switch(dir)
170  {
171  case 0:
172  {
173  PhysDeriv(inarray, outarray, NullNekDouble1DArray, NullNekDouble1DArray);
174  }
175  break;
176  case 1:
177  {
178  PhysDeriv(inarray, NullNekDouble1DArray, outarray, NullNekDouble1DArray);
179  }
180  break;
181  case 2:
182  {
183  PhysDeriv(inarray, NullNekDouble1DArray, NullNekDouble1DArray, outarray);
184  }
185  break;
186  default:
187  {
188  ASSERTL1(false,"input dir is out of range");
189  }
190  break;
191  }
192  }
193 
194  void TriExp::v_PhysDirectionalDeriv(
195  const Array<OneD, const NekDouble> &inarray,
196  const Array<OneD, const NekDouble> &direction,
197  Array<OneD, NekDouble> &out)
198  {
199  if(! out.num_elements())
200  {
201  return;
202  }
203 
204  int nquad0 = m_base[0]->GetNumPoints();
205  int nquad1 = m_base[1]->GetNumPoints();
206  int nqtot = nquad0*nquad1;
207 
208  const Array<TwoD, const NekDouble>& df =
209  m_metricinfo->GetDerivFactors(GetPointsKeys());
210 
211  Array<OneD,NekDouble> diff0(2*nqtot);
212  Array<OneD,NekDouble> diff1(diff0+nqtot);
213 
214  // diff0 = du/d_xi, diff1 = du/d_eta
215  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
216 
217  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
218  {
219  Array<OneD, Array<OneD, NekDouble> > tangmat(2);
220 
221 
222  // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
223  // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
224  for (int i=0; i< 2; ++i)
225  {
226  tangmat[i] = Array<OneD, NekDouble>(nqtot,0.0);
227  for (int k=0; k<(m_geom->GetCoordim()); ++k)
228  {
229  Vmath::Vvtvp(nqtot,&df[2*k+i][0],1,&direction[k*nqtot],1,&tangmat[i][0],1,&tangmat[i][0],1);
230  }
231  }
232 
233  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
234  Vmath::Vmul (nqtot,&tangmat[0][0],1,&diff0[0],1, &out[0], 1);
235  Vmath::Vvtvp (nqtot,&tangmat[1][0],1,&diff1[0],1, &out[0], 1, &out[0],1);
236  }
237  else
238  {
239  ASSERTL1(m_metricinfo->GetGtype() == SpatialDomains::eDeformed,"Wrong route");
240  }
241  }
242 
243 
244  void TriExp::v_FwdTrans(const Array<OneD, const NekDouble> & inarray,
245  Array<OneD,NekDouble> &outarray)
246  {
247  IProductWRTBase(inarray,outarray);
248 
249  // get Mass matrix inverse
250  MatrixKey masskey(StdRegions::eInvMass,
251  DetShapeType(),*this);
252  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
253 
254  // copy inarray in case inarray == outarray
255  NekVector<NekDouble> in (m_ncoeffs,outarray,eCopy);
256  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
257 
258  out = (*matsys)*in;
259  }
260 
261 
262  void TriExp::v_FwdTrans_BndConstrained(const Array<OneD, const NekDouble>& inarray,
263  Array<OneD, NekDouble> &outarray)
264  {
265  int i,j;
266  int npoints[2] = {m_base[0]->GetNumPoints(),
267  m_base[1]->GetNumPoints()};
268  int nmodes[2] = {m_base[0]->GetNumModes(),
269  m_base[1]->GetNumModes()};
270 
271  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
272 
273  Array<OneD, NekDouble> physEdge[3];
274  Array<OneD, NekDouble> coeffEdge[3];
275  StdRegions::Orientation orient[3];
276  for(i = 0; i < 3; i++)
277  {
278  physEdge[i] = Array<OneD, NekDouble>(npoints[i!=0]);
279  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i!=0]);
280  orient[i] = GetEorient(i);
281  }
282 
283  for(i = 0; i < npoints[0]; i++)
284  {
285  physEdge[0][i] = inarray[i];
286  }
287 
288  for(i = 0; i < npoints[1]; i++)
289  {
290  physEdge[1][i] = inarray[npoints[0]-1+i*npoints[0]];
291  physEdge[2][i] = inarray[(npoints[1]-1)*npoints[0]-i*npoints[0]];
292  }
293 
294  for(i = 0; i < 3; i++)
295  {
296  if( orient[i] == StdRegions::eBackwards )
297  {
298  reverse( (physEdge[i]).get() , (physEdge[i]).get() + npoints[i!=0] );
299  }
300  }
301 
302  SegExpSharedPtr segexp[3];
303  for(i = 0; i < 3; i++)
304  {
305  segexp[i] = MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(m_base[i!=0]->GetBasisKey(),GetGeom2D()->GetEdge(i));
306  }
307 
308  Array<OneD, unsigned int> mapArray;
309  Array<OneD, int> signArray;
310  NekDouble sign;
311 
312  for(i = 0; i < 3; i++)
313  {
314  segexp[i!=0]->FwdTrans_BndConstrained(physEdge[i],coeffEdge[i]);
315 
316  GetEdgeToElementMap(i,orient[i],mapArray,signArray);
317  for(j=0; j < nmodes[i!=0]; j++)
318  {
319  sign = (NekDouble) signArray[j];
320  outarray[ mapArray[j] ] = sign * coeffEdge[i][j];
321  }
322  }
323 
324  int nBoundaryDofs = NumBndryCoeffs();
325  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
326 
327  if (nInteriorDofs > 0) {
328  Array<OneD, NekDouble> tmp0(m_ncoeffs);
329  Array<OneD, NekDouble> tmp1(m_ncoeffs);
330 
331  StdRegions::StdMatrixKey stdmasskey(StdRegions::eMass,DetShapeType(),*this);
332  MassMatrixOp(outarray,tmp0,stdmasskey);
333  IProductWRTBase(inarray,tmp1);
334 
335  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
336 
337  // get Mass matrix inverse (only of interior DOF)
338  // use block (1,1) of the static condensed system
339  // note: this block alreay contains the inverse matrix
340  MatrixKey masskey(StdRegions::eMass,DetShapeType(),*this);
341  DNekScalMatSharedPtr matsys = (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
342 
343  Array<OneD, NekDouble> rhs(nInteriorDofs);
344  Array<OneD, NekDouble> result(nInteriorDofs);
345 
346  GetInteriorMap(mapArray);
347 
348  for(i = 0; i < nInteriorDofs; i++)
349  {
350  rhs[i] = tmp1[ mapArray[i] ];
351  }
352 
353  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(), &((matsys->GetOwnedMatrix())->GetPtr())[0],
354  nInteriorDofs,rhs.get(),1,0.0,result.get(),1);
355 
356  for(i = 0; i < nInteriorDofs; i++)
357  {
358  outarray[ mapArray[i] ] = result[i];
359  }
360  }
361  }
362 
363 
364  void TriExp::v_IProductWRTBase(const Array<OneD, const NekDouble>& inarray,
365  Array<OneD, NekDouble> &outarray)
366  {
367  IProductWRTBase_SumFac(inarray,outarray);
368  }
369 
370 
371  void TriExp::v_IProductWRTDerivBase(const int dir,
372  const Array<OneD, const NekDouble>& inarray,
373  Array<OneD, NekDouble> & outarray)
374  {
375  IProductWRTDerivBase_SumFac(dir,inarray,outarray);
376  }
377 
378 
379  void TriExp::v_IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
380  Array<OneD, NekDouble> &outarray)
381  {
382  int nquad0 = m_base[0]->GetNumPoints();
383  int nquad1 = m_base[1]->GetNumPoints();
384  int order0 = m_base[0]->GetNumModes();
385 
386  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad1*order0);
387  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
388 
389  MultiplyByQuadratureMetric(inarray,tmp);
390  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),tmp,outarray,wsp);
391  }
392 
393 
394  void TriExp::v_IProductWRTBase_MatOp(const Array<OneD, const NekDouble>& inarray,
395  Array<OneD, NekDouble> &outarray)
396  {
397  int nq = GetTotPoints();
398  MatrixKey iprodmatkey(StdRegions::eIProductWRTBase,DetShapeType(),*this);
399  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
400 
401  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
402  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
403 
404  }
405 
406 
407  void TriExp::v_IProductWRTDerivBase_SumFac(const int dir,
408  const Array<OneD, const NekDouble>& inarray,
409  Array<OneD, NekDouble> & outarray)
410  {
411  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
412  ASSERTL1((dir==2)?(m_geom->GetCoordim()==3):true,"Invalid direction.");
413 
414  int i;
415  int nquad0 = m_base[0]->GetNumPoints();
416  int nquad1 = m_base[1]->GetNumPoints();
417  int nqtot = nquad0*nquad1;
418  int nmodes0 = m_base[0]->GetNumModes();
419  int wspsize = max(max(nqtot,m_ncoeffs),nquad1*nmodes0);
420 
421  const Array<TwoD, const NekDouble>& df =
422  m_metricinfo->GetDerivFactors(GetPointsKeys());
423 
424  Array<OneD, NekDouble> tmp0 (6*wspsize);
425  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
426  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
427  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
428  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
429  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
430 
431  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
432  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
433 
434  // set up geometric factor: 2/(1-z1)
435  for(i = 0; i < nquad1; ++i)
436  {
437  gfac0[i] = 2.0/(1-z1[i]);
438  }
439  for(i = 0; i < nquad0; ++i)
440  {
441  gfac1[i] = 0.5*(1+z0[i]);
442  }
443 
444  for(i = 0; i < nquad1; ++i)
445  {
446  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
447  }
448 
449  for(i = 0; i < nquad1; ++i)
450  {
451  Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
452  }
453 
454  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
455  {
456  Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
457  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
458  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
459  }
460  else
461  {
462  Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
463  Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
464  Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
465  }
466  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
467 
468  MultiplyByQuadratureMetric(tmp1,tmp1);
469  MultiplyByQuadratureMetric(tmp2,tmp2);
470 
471  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),m_base[1]->GetBdata() ,tmp1,tmp3 ,tmp0);
472  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata() ,m_base[1]->GetDbdata(),tmp2,outarray,tmp0);
473  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
474  }
475 
476 
477  void TriExp::v_IProductWRTDerivBase_MatOp(const int dir,
478  const Array<OneD, const NekDouble>& inarray,
479  Array<OneD, NekDouble> &outarray)
480  {
481  int nq = GetTotPoints();
482  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
483 
484  switch(dir)
485  {
486  case 0:
487  {
488  mtype = StdRegions::eIProductWRTDerivBase0;
489  }
490  break;
491  case 1:
492  {
493  mtype = StdRegions::eIProductWRTDerivBase1;
494  }
495  break;
496  case 2:
497  {
498  mtype = StdRegions::eIProductWRTDerivBase2;
499  }
500  break;
501  default:
502  {
503  ASSERTL1(false,"input dir is out of range");
504  }
505  break;
506  }
507 
508  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
509  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
510 
511  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
512  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
513 
514  }
515 
516  void TriExp::v_NormVectorIProductWRTBase(
517  const Array<OneD, const NekDouble> &Fx,
518  const Array<OneD, const NekDouble> &Fy,
519  const Array<OneD, const NekDouble> &Fz,
520  Array<OneD, NekDouble> &outarray)
521  {
522  int nq = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
523  Array<OneD, NekDouble > Fn(nq);
524 
525  const Array<OneD, const Array<OneD, NekDouble> > &normals =
526  GetLeftAdjacentElementExp()->GetFaceNormal(
527  GetLeftAdjacentElementFace());
528 
529  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
530  {
531  Vmath::Vvtvvtp(nq,&normals[0][0],1,&Fx[0],1,
532  &normals[1][0],1,&Fy[0],1,&Fn[0],1);
533  Vmath::Vvtvp (nq,&normals[2][0],1,&Fz[0],1,&Fn[0],1,&Fn[0],1);
534  }
535  else
536  {
537  Vmath::Svtsvtp(nq,normals[0][0],&Fx[0],1,
538  normals[1][0],&Fy[0],1,&Fn[0],1);
539  Vmath::Svtvp (nq,normals[2][0],&Fz[0],1,&Fn[0],1,&Fn[0],1);
540  }
541 
542  IProductWRTBase(Fn,outarray);
543  }
544 
545  void TriExp::v_GetCoord(const Array<OneD, const NekDouble> &Lcoords,
546  Array<OneD,NekDouble> &coords)
547  {
548  int i;
549 
550  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
551  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
552  "Local coordinates are not in region [-1,1]");
553 
554  m_geom->FillGeom();
555 
556  for(i = 0; i < m_geom->GetCoordim(); ++i)
557  {
558  coords[i] = m_geom->GetCoord(i,Lcoords);
559  }
560  }
561 
562  void TriExp::v_GetCoords(
563  Array<OneD, NekDouble> &coords_0,
564  Array<OneD, NekDouble> &coords_1,
565  Array<OneD, NekDouble> &coords_2)
566  {
567  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
568  }
569 
570 
571  /**
572  * Given the local cartesian coordinate \a Lcoord evaluate the
573  * value of physvals at this point by calling through to the
574  * StdExpansion method
575  */
576  NekDouble TriExp::v_StdPhysEvaluate(
577  const Array<OneD, const NekDouble> &Lcoord,
578  const Array<OneD, const NekDouble> &physvals)
579  {
580  // Evaluate point in local (eta) coordinates.
581  return StdTriExp::v_PhysEvaluate(Lcoord,physvals);
582  }
583 
584  NekDouble TriExp::v_PhysEvaluate(const Array<OneD, const NekDouble> &coord, const Array<OneD, const NekDouble> & physvals)
585  {
586  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
587 
588  ASSERTL0(m_geom,"m_geom not defined");
589  m_geom->GetLocCoords(coord,Lcoord);
590 
591  return StdTriExp::v_PhysEvaluate(Lcoord, physvals);
592  }
593 
594 
595  void TriExp::v_GetTracePhysVals(
596  const int edge,
597  const StdRegions::StdExpansionSharedPtr &EdgeExp,
598  const Array<OneD, const NekDouble> &inarray,
599  Array<OneD,NekDouble> &outarray,
600  StdRegions::Orientation orient)
601  {
602  v_GetEdgePhysVals(edge,EdgeExp,inarray,outarray);
603  }
604 
605  void TriExp::v_GetEdgePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
606  const Array<OneD, const NekDouble> &inarray,
607  Array<OneD,NekDouble> &outarray)
608  {
609  int nquad0 = m_base[0]->GetNumPoints();
610  int nquad1 = m_base[1]->GetNumPoints();
611 
612  // get points in Cartesian orientation
613  switch(edge)
614  {
615  case 0:
616  Vmath::Vcopy(nquad0,&(inarray[0]),1,&(outarray[0]),1);
617  break;
618  case 1:
619  Vmath::Vcopy(nquad1,&(inarray[0])+(nquad0-1),nquad0,&(outarray[0]),1);
620  break;
621  case 2:
622  Vmath::Vcopy(nquad1,&(inarray[0]),nquad0,&(outarray[0]),1);
623  break;
624  default:
625  ASSERTL0(false,"edge value (< 3) is out of range");
626  break;
627  }
628 
629  // Interpolate if required
630  if(m_base[edge?1:0]->GetPointsKey() != EdgeExp->GetBasis(0)->GetPointsKey())
631  {
632  Array<OneD,NekDouble> outtmp(max(nquad0,nquad1));
633 
634  outtmp = outarray;
635 
636  LibUtilities::Interp1D(m_base[edge?1:0]->GetPointsKey(),outtmp,
637  EdgeExp->GetBasis(0)->GetPointsKey(),outarray);
638  }
639 
640  //Reverse data if necessary
641  if(GetCartesianEorient(edge) == StdRegions::eBackwards)
642  {
643  Vmath::Reverse(EdgeExp->GetNumPoints(0),&outarray[0],1,
644  &outarray[0],1);
645  }
646 
647  }
648 
649 
650  void TriExp::v_GetEdgeInterpVals(
651  const int edge,const Array<OneD, const NekDouble> &inarray,
652  Array<OneD, NekDouble> &outarray)
653  {
654  ASSERTL0(false,
655  "Routine not implemented for triangular elements");
656  }
657 
658  void TriExp::v_GetEdgeQFactors(
659  const int edge,
660  Array<OneD, NekDouble> &outarray)
661  {
662  ASSERTL0(false,
663  "Routine not implemented for triangular elements");
664  }
665 
666 
667  void TriExp::v_ComputeEdgeNormal(const int edge)
668  {
669  int i;
670  const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
671  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
672  const SpatialDomains::GeomType type = geomFactors->GetGtype();
673  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
674  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
675  int nqe = m_base[0]->GetNumPoints();
676  int dim = GetCoordim();
677 
678  m_edgeNormals[edge] = Array<OneD, Array<OneD, NekDouble> >(dim);
679  Array<OneD, Array<OneD, NekDouble> > &normal = m_edgeNormals[edge];
680  for (i = 0; i < dim; ++i)
681  {
682  normal[i] = Array<OneD, NekDouble>(nqe);
683  }
684 
685  // Regular geometry case
686  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
687  {
688  NekDouble fac;
689  // Set up normals
690  switch(edge)
691  {
692  case 0:
693  for(i = 0; i < GetCoordim(); ++i)
694  {
695  Vmath::Fill(nqe,-df[2*i+1][0],normal[i],1);
696  }
697  break;
698  case 1:
699  for(i = 0; i < GetCoordim(); ++i)
700  {
701  Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
702  }
703  break;
704  case 2:
705  for(i = 0; i < GetCoordim(); ++i)
706  {
707  Vmath::Fill(nqe,-df[2*i][0],normal[i],1);
708  }
709  break;
710  default:
711  ASSERTL0(false,"Edge is out of range (edge < 3)");
712  }
713 
714  // normalise
715  fac = 0.0;
716  for(i =0 ; i < GetCoordim(); ++i)
717  {
718  fac += normal[i][0]*normal[i][0];
719  }
720  fac = 1.0/sqrt(fac);
721  for (i = 0; i < GetCoordim(); ++i)
722  {
723  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
724  }
725  }
726  else // Set up deformed normals
727  {
728  int j;
729 
730  int nquad0 = ptsKeys[0].GetNumPoints();
731  int nquad1 = ptsKeys[1].GetNumPoints();
732 
733  LibUtilities::PointsKey from_key;
734 
735  Array<OneD,NekDouble> normals(GetCoordim()*max(nquad0,nquad1),0.0);
736  Array<OneD,NekDouble> edgejac(GetCoordim()*max(nquad0,nquad1),0.0);
737 
738  // Extract Jacobian along edges and recover local
739  // derivates (dx/dr) for polynomial interpolation by
740  // multiplying m_gmat by jacobian
741  switch(edge)
742  {
743  case 0:
744  for(j = 0; j < nquad0; ++j)
745  {
746  edgejac[j] = jac[j];
747  for(i = 0; i < GetCoordim(); ++i)
748  {
749  normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
750  }
751  }
752  from_key = ptsKeys[0];
753  break;
754  case 1:
755  for(j = 0; j < nquad1; ++j)
756  {
757  edgejac[j] = jac[nquad0*j+nquad0-1];
758  for(i = 0; i < GetCoordim(); ++i)
759  {
760  normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
761  }
762  }
763  from_key = ptsKeys[1];
764  break;
765  case 2:
766  for(j = 0; j < nquad1; ++j)
767  {
768  edgejac[j] = jac[nquad0*j];
769  for(i = 0; i < GetCoordim(); ++i)
770  {
771  normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
772  }
773  }
774  from_key = ptsKeys[1];
775  break;
776  default:
777  ASSERTL0(false,"edge is out of range (edge < 3)");
778 
779  }
780 
781  int nq = from_key.GetNumPoints();
782  Array<OneD,NekDouble> work(nqe,0.0);
783 
784  // interpolate Jacobian and invert
785  LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
786  Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
787 
788  // interpolate
789  for(i = 0; i < GetCoordim(); ++i)
790  {
791  LibUtilities::Interp1D(from_key,&normals[i*nq],m_base[0]->GetPointsKey(),&normal[i][0]);
792  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
793  }
794 
795  //normalise normal vectors
796  Vmath::Zero(nqe,work,1);
797  for(i = 0; i < GetCoordim(); ++i)
798  {
799  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
800  }
801 
802  Vmath::Vsqrt(nqe,work,1,work,1);
803  Vmath::Sdiv(nqe,1.0,work,1,work,1);
804 
805  for(i = 0; i < GetCoordim(); ++i)
806  {
807  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
808  }
809 
810  // Reverse direction so that points are in
811  // anticlockwise direction if edge >=2
812  if(edge >= 2)
813  {
814  for(i = 0; i < GetCoordim(); ++i)
815  {
816  Vmath::Reverse(nqe,normal[i],1, normal[i],1);
817  }
818  }
819  }
820  if(GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
821  {
822  for(i = 0; i < GetCoordim(); ++i)
823  {
824  if(geomFactors->GetGtype() == SpatialDomains::eDeformed)
825  {
826  Vmath::Reverse(nqe, normal[i], 1, normal[i],1);
827  }
828  }
829  }
830  }
831 
832  int TriExp::v_GetCoordim()
833  {
834  return m_geom->GetCoordim();
835  }
836 
837 
838  void TriExp::v_ExtractDataToCoeffs(const NekDouble *data,
839  const std::vector<unsigned int > &nummodes, const int mode_offset, NekDouble * coeffs)
840  {
841  int data_order0 = nummodes[mode_offset];
842  int fillorder0 = min(m_base[0]->GetNumModes(),data_order0);
843  int data_order1 = nummodes[mode_offset+1];
844  int order1 = m_base[1]->GetNumModes();
845  int fillorder1 = min(order1,data_order1);
846 
847  switch(m_base[0]->GetBasisType())
848  {
849  case LibUtilities::eModified_A:
850  {
851  int i;
852  int cnt = 0;
853  int cnt1 = 0;
854 
855  ASSERTL1(m_base[1]->GetBasisType() == LibUtilities::eModified_B,
856  "Extraction routine not set up for this basis");
857 
858  Vmath::Zero(m_ncoeffs,coeffs,1);
859  for(i = 0; i < fillorder0; ++i)
860  {
861  Vmath::Vcopy(fillorder1-i,&data[cnt],1,&coeffs[cnt1],1);
862  cnt += data_order1-i;
863  cnt1 += order1-i;
864  }
865  }
866  break;
867  default:
868  ASSERTL0(false,"basis is either not set up or not hierarchicial");
869  }
870  }
871 
872 
873  StdRegions::Orientation TriExp::v_GetEorient(int edge)
874  {
875  return GetGeom2D()->GetEorient(edge);
876  }
877 
878 
879  StdRegions::Orientation TriExp::v_GetCartesianEorient(int edge)
880  {
881  return GetGeom2D()->GetCartesianEorient(edge);
882  }
883 
884 
885  const LibUtilities::BasisSharedPtr& TriExp::v_GetBasis(int dir) const
886  {
887  ASSERTL1(dir >= 0 &&dir <= 1,"input dir is out of range");
888  return m_base[dir];
889  }
890 
891 
892  int TriExp::v_GetNumPoints(const int dir) const
893  {
894  return GetNumPoints(dir);
895  }
896 
897 
898  DNekMatSharedPtr TriExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
899  {
900  DNekMatSharedPtr returnval;
901  switch(mkey.GetMatrixType())
902  {
903  case StdRegions::eHybridDGHelmholtz:
904  case StdRegions::eHybridDGLamToU:
905  case StdRegions::eHybridDGLamToQ0:
906  case StdRegions::eHybridDGLamToQ1:
907  case StdRegions::eHybridDGLamToQ2:
908  case StdRegions::eHybridDGHelmBndLam:
909  case StdRegions::eInvLaplacianWithUnityMean:
910  returnval = Expansion2D::v_GenMatrix(mkey);
911  break;
912  default:
913  returnval = StdTriExp::v_GenMatrix(mkey);
914  break;
915  }
916 
917  return returnval;
918  }
919 
920 
921  DNekMatSharedPtr TriExp::v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
922  {
923  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
924  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
925  StdRegions::StdTriExpSharedPtr tmp = MemoryManager<StdTriExp>::
926  AllocateSharedPtr(bkey0,bkey1);
927 
928  return tmp->GetStdMatrix(mkey);
929  }
930 
931 
932  DNekScalMatSharedPtr TriExp::CreateMatrix(const MatrixKey &mkey)
933  {
934  DNekScalMatSharedPtr returnval;
935  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
936 
937  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
938 
939  switch(mkey.GetMatrixType())
940  {
941  case StdRegions::eMass:
942  {
943  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
944  (mkey.GetNVarCoeff()))
945  {
946  NekDouble one = 1.0;
947  DNekMatSharedPtr mat = GenMatrix(mkey);
948  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
949  }
950  else
951  {
952  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
953  DNekMatSharedPtr mat = GetStdMatrix(mkey);
954  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
955  }
956  }
957  break;
958  case StdRegions::eInvMass:
959  {
960  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
961  {
962  NekDouble one = 1.0;
963  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),
964  *this);
965  DNekMatSharedPtr mat = GenMatrix(masskey);
966  mat->Invert();
967 
968  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
969  }
970  else
971  {
972  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
973  DNekMatSharedPtr mat = GetStdMatrix(mkey);
974  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
975 
976  }
977  }
978  break;
979  case StdRegions::eWeakDeriv0:
980  case StdRegions::eWeakDeriv1:
981  case StdRegions::eWeakDeriv2:
982  {
983  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed || mkey.GetNVarCoeff())
984  {
985  NekDouble one = 1.0;
986  DNekMatSharedPtr mat = GenMatrix(mkey);
987 
988  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
989  }
990  else
991  {
992  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
993  Array<TwoD, const NekDouble> df = m_metricinfo->GetDerivFactors(ptsKeys);
994  int dir = 0;
995  switch(mkey.GetMatrixType())
996  {
997  case StdRegions::eWeakDeriv0:
998  dir = 0;
999  break;
1000  case StdRegions::eWeakDeriv1:
1001  dir = 1;
1002  break;
1003  case StdRegions::eWeakDeriv2:
1004  dir = 2;
1005  break;
1006  default:
1007  break;
1008  }
1009 
1010  MatrixKey deriv0key(StdRegions::eWeakDeriv0,
1011  mkey.GetShapeType(), *this);
1012  MatrixKey deriv1key(StdRegions::eWeakDeriv1,
1013  mkey.GetShapeType(), *this);
1014 
1015  DNekMat &deriv0 = *GetStdMatrix(deriv0key);
1016  DNekMat &deriv1 = *GetStdMatrix(deriv1key);
1017 
1018  int rows = deriv0.GetRows();
1019  int cols = deriv1.GetColumns();
1020 
1021  DNekMatSharedPtr WeakDeriv = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1022  (*WeakDeriv) = df[2*dir][0]*deriv0 + df[2*dir+1][0]*deriv1;
1023 
1024  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,WeakDeriv);
1025  }
1026  }
1027  break;
1028  case StdRegions::eLaplacian:
1029  {
1030  if( (m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
1031  (mkey.GetNVarCoeff() > 0)||(mkey.ConstFactorExists(StdRegions::eFactorSVVCutoffRatio)))
1032  {
1033  NekDouble one = 1.0;
1034  DNekMatSharedPtr mat = GenMatrix(mkey);
1035 
1036  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1037  }
1038  else
1039  {
1040  MatrixKey lap00key(StdRegions::eLaplacian00,
1041  mkey.GetShapeType(), *this);
1042  MatrixKey lap01key(StdRegions::eLaplacian01,
1043  mkey.GetShapeType(), *this);
1044  MatrixKey lap11key(StdRegions::eLaplacian11,
1045  mkey.GetShapeType(), *this);
1046 
1047  DNekMat &lap00 = *GetStdMatrix(lap00key);
1048  DNekMat &lap01 = *GetStdMatrix(lap01key);
1049  DNekMat &lap11 = *GetStdMatrix(lap11key);
1050 
1051  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1052  Array<TwoD, const NekDouble> gmat =
1053  m_metricinfo->GetGmat(ptsKeys);
1054 
1055  int rows = lap00.GetRows();
1056  int cols = lap00.GetColumns();
1057 
1058  DNekMatSharedPtr lap = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1059 
1060  (*lap) = gmat[0][0] * lap00 +
1061  gmat[1][0] * (lap01 + Transpose(lap01)) +
1062  gmat[3][0] * lap11;
1063 
1064  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,lap);
1065  }
1066  }
1067  break;
1068  case StdRegions::eInvLaplacianWithUnityMean:
1069  {
1070  DNekMatSharedPtr mat = GenMatrix(mkey);
1071  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0,mat);
1072  }
1073  break;
1074  case StdRegions::eHelmholtz:
1075  {
1076  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
1077 
1078  MatrixKey masskey(mkey, StdRegions::eMass);
1079  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1080 
1081  MatrixKey lapkey(mkey, StdRegions::eLaplacian);
1082  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1083 
1084  int rows = LapMat.GetRows();
1085  int cols = LapMat.GetColumns();
1086 
1087  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1088 
1089  NekDouble one = 1.0;
1090  (*helm) = LapMat + factor*MassMat;
1091 
1092  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
1093  }
1094  break;
1095  case StdRegions::eIProductWRTBase:
1096  {
1097  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1098  {
1099  NekDouble one = 1.0;
1100  DNekMatSharedPtr mat = GenMatrix(mkey);
1101  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1102  }
1103  else
1104  {
1105  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1106  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1107  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1108  }
1109  }
1110  break;
1111  case StdRegions::eIProductWRTDerivBase0:
1112  case StdRegions::eIProductWRTDerivBase1:
1113  case StdRegions::eIProductWRTDerivBase2:
1114  {
1115  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1116  {
1117  NekDouble one = 1.0;
1118  DNekMatSharedPtr mat = GenMatrix(mkey);
1119  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1120  }
1121  else
1122  {
1123  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
1124 
1125  const Array<TwoD, const NekDouble>& df = m_metricinfo->GetDerivFactors(ptsKeys);
1126  int dir = 0;
1127 
1128  switch(mkey.GetMatrixType())
1129  {
1130  case StdRegions::eIProductWRTDerivBase0:
1131  dir = 0;
1132  break;
1133  case StdRegions::eIProductWRTDerivBase1:
1134  dir = 1;
1135  break;
1136  case StdRegions::eIProductWRTDerivBase2:
1137  dir = 2;
1138  break;
1139  default:
1140  break;
1141  }
1142 
1143  MatrixKey iProdDeriv0Key(StdRegions::eIProductWRTDerivBase0,
1144  mkey.GetShapeType(), *this);
1145  MatrixKey iProdDeriv1Key(StdRegions::eIProductWRTDerivBase1,
1146  mkey.GetShapeType(), *this);
1147 
1148  DNekMat &stdiprod0 = *GetStdMatrix(iProdDeriv0Key);
1149  DNekMat &stdiprod1 = *GetStdMatrix(iProdDeriv0Key);
1150 
1151  int rows = stdiprod0.GetRows();
1152  int cols = stdiprod1.GetColumns();
1153 
1154  DNekMatSharedPtr mat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
1155  (*mat) = df[2*dir][0]*stdiprod0 + df[2*dir+1][0]*stdiprod1;
1156 
1157  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
1158  }
1159  }
1160  break;
1161 
1162  case StdRegions::eInvHybridDGHelmholtz:
1163  {
1164  NekDouble one = 1.0;
1165 
1166  MatrixKey hkey(StdRegions::eHybridDGHelmholtz, DetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1167 
1168  DNekMatSharedPtr mat = GenMatrix(hkey);
1169 
1170  mat->Invert();
1171  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1172  }
1173  break;
1174  case StdRegions::ePreconLinearSpace:
1175  {
1176  NekDouble one = 1.0;
1177  MatrixKey helmkey(StdRegions::eHelmholtz, mkey.GetShapeType(), *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1178  DNekScalBlkMatSharedPtr helmStatCond = GetLocStaticCondMatrix(helmkey);
1179  DNekScalMatSharedPtr A =helmStatCond->GetBlock(0,0);
1180  DNekMatSharedPtr R=BuildVertexMatrix(A);
1181 
1182  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,R);
1183  }
1184  break;
1185  default:
1186  {
1187  NekDouble one = 1.0;
1188  DNekMatSharedPtr mat = GenMatrix(mkey);
1189 
1190  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
1191  }
1192  break;
1193  }
1194 
1195  return returnval;
1196  }
1197 
1198 
1199  DNekScalBlkMatSharedPtr TriExp::CreateStaticCondMatrix(const MatrixKey &mkey)
1200  {
1201  DNekScalBlkMatSharedPtr returnval;
1202  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
1203 
1204  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
1205 
1206  // set up block matrix system
1207  unsigned int nbdry = NumBndryCoeffs();
1208  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
1209  unsigned int exp_size[] = {nbdry,nint};
1210  unsigned int nblks = 2;
1211  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
1212  NekDouble factor = 1.0;
1213 
1214  switch(mkey.GetMatrixType())
1215  {
1216  // this can only use stdregions statically condensed system for mass matrix
1217  case StdRegions::eMass:
1218  if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||(mkey.GetNVarCoeff()))
1219  {
1220  factor = 1.0;
1221  goto UseLocRegionsMatrix;
1222  }
1223  else
1224  {
1225  factor = (m_metricinfo->GetJac(ptsKeys))[0];
1226  goto UseStdRegionsMatrix;
1227  }
1228  break;
1229  default: // use Deformed case for both regular and deformed geometries
1230  factor = 1.0;
1231  goto UseLocRegionsMatrix;
1232  break;
1233  UseStdRegionsMatrix:
1234  {
1235  NekDouble invfactor = 1.0/factor;
1236  NekDouble one = 1.0;
1237  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
1238  DNekScalMatSharedPtr Atmp;
1239  DNekMatSharedPtr Asubmat;
1240 
1241  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
1242  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
1243  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
1244  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
1245  }
1246  break;
1247 
1248  UseLocRegionsMatrix:
1249  {
1250  int i,j;
1251  NekDouble invfactor = 1.0/factor;
1252  NekDouble one = 1.0;
1253 
1254  DNekScalMat &mat = *GetLocMatrix(mkey);
1255 
1256  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
1257  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
1258  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
1259  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
1260 
1261  Array<OneD,unsigned int> bmap(nbdry);
1262  Array<OneD,unsigned int> imap(nint);
1263  GetBoundaryMap(bmap);
1264  GetInteriorMap(imap);
1265 
1266  for(i = 0; i < nbdry; ++i)
1267  {
1268  for(j = 0; j < nbdry; ++j)
1269  {
1270  (*A)(i,j) = mat(bmap[i],bmap[j]);
1271  }
1272 
1273  for(j = 0; j < nint; ++j)
1274  {
1275  (*B)(i,j) = mat(bmap[i],imap[j]);
1276  }
1277  }
1278 
1279  for(i = 0; i < nint; ++i)
1280  {
1281  for(j = 0; j < nbdry; ++j)
1282  {
1283  (*C)(i,j) = mat(imap[i],bmap[j]);
1284  }
1285 
1286  for(j = 0; j < nint; ++j)
1287  {
1288  (*D)(i,j) = mat(imap[i],imap[j]);
1289  }
1290  }
1291 
1292  // Calculate static condensed system
1293  if(nint)
1294  {
1295  D->Invert();
1296  (*B) = (*B)*(*D);
1297  (*A) = (*A) - (*B)*(*C);
1298  }
1299 
1300  DNekScalMatSharedPtr Atmp;
1301 
1302  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
1303  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
1304  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
1305  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
1306 
1307  }
1308  }
1309 
1310  return returnval;
1311  }
1312 
1313 
1314  DNekScalMatSharedPtr TriExp::v_GetLocMatrix(const MatrixKey &mkey)
1315  {
1316  return m_matrixManager[mkey];
1317  }
1318 
1319 
1320  DNekScalBlkMatSharedPtr TriExp::v_GetLocStaticCondMatrix(const MatrixKey &mkey)
1321  {
1322  return m_staticCondMatrixManager[mkey];
1323  }
1324 
1325  void TriExp::v_DropLocStaticCondMatrix(const MatrixKey &mkey)
1326  {
1327  m_staticCondMatrixManager.DeleteObject(mkey);
1328  }
1329 
1330 
1331 
1332  void TriExp::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
1333  Array<OneD,NekDouble> &outarray,
1334  const StdRegions::StdMatrixKey &mkey)
1335  {
1336  StdExpansion::MassMatrixOp_MatFree(inarray,outarray,mkey);
1337  }
1338 
1339 
1340  void TriExp::v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
1341  Array<OneD,NekDouble> &outarray,
1342  const StdRegions::StdMatrixKey &mkey)
1343  {
1344  TriExp::LaplacianMatrixOp_MatFree(inarray,outarray,mkey);
1345  }
1346 
1347 
1348  void TriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1349  const Array<OneD, const NekDouble> &inarray,
1350  Array<OneD,NekDouble> &outarray,
1351  const StdRegions::StdMatrixKey &mkey)
1352  {
1353  StdExpansion::LaplacianMatrixOp_MatFree(k1,k2,inarray,outarray,mkey);
1354  }
1355 
1356 
1357  void TriExp::v_WeakDerivMatrixOp(const int i,
1358  const Array<OneD, const NekDouble> &inarray,
1359  Array<OneD,NekDouble> &outarray,
1360  const StdRegions::StdMatrixKey &mkey)
1361  {
1362  StdExpansion::WeakDerivMatrixOp_MatFree(i,inarray,outarray,mkey);
1363  }
1364 
1365 
1366  void TriExp::v_WeakDirectionalDerivMatrixOp(const Array<OneD, const NekDouble> &inarray,
1367  Array<OneD,NekDouble> &outarray,
1368  const StdRegions::StdMatrixKey &mkey)
1369  {
1370  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray,outarray,mkey);
1371  }
1372 
1373 
1374  void TriExp::v_MassLevelCurvatureMatrixOp(const Array<OneD, const NekDouble> &inarray,
1375  Array<OneD,NekDouble> &outarray,
1376  const StdRegions::StdMatrixKey &mkey)
1377  {
1378  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray,outarray,mkey);
1379  }
1380 
1381 
1382  void TriExp::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
1383  Array<OneD,NekDouble> &outarray,
1384  const StdRegions::StdMatrixKey &mkey)
1385  {
1386  TriExp::HelmholtzMatrixOp_MatFree(inarray,outarray,mkey);
1387  }
1388 
1389 
1390  void TriExp::v_GeneralMatrixOp_MatOp(const Array<OneD, const NekDouble> &inarray,
1391  Array<OneD,NekDouble> &outarray,
1392  const StdRegions::StdMatrixKey &mkey)
1393  {
1394  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1395 
1396  if(inarray.get() == outarray.get())
1397  {
1398  Array<OneD,NekDouble> tmp(m_ncoeffs);
1399  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
1400 
1401  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1402  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1403  }
1404  else
1405  {
1406  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
1407  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1408  }
1409  }
1410 
1411 
1412  void TriExp::v_LaplacianMatrixOp_MatFree_Kernel(
1413  const Array<OneD, const NekDouble> &inarray,
1414  Array<OneD, NekDouble> &outarray,
1415  Array<OneD, NekDouble> &wsp)
1416  {
1417  if (m_metrics.count(MetricLaplacian00) == 0)
1418  {
1419  ComputeLaplacianMetric();
1420  }
1421 
1422  int nquad0 = m_base[0]->GetNumPoints();
1423  int nquad1 = m_base[1]->GetNumPoints();
1424  int nqtot = nquad0*nquad1;
1425  int nmodes0 = m_base[0]->GetNumModes();
1426  int nmodes1 = m_base[1]->GetNumModes();
1427  int wspsize = max(max(max(nqtot,m_ncoeffs),nquad1*nmodes0),nquad0*nmodes1);
1428 
1429  ASSERTL1(wsp.num_elements() >= 3*wspsize,
1430  "Workspace is of insufficient size.");
1431 
1432  const Array<OneD, const NekDouble>& base0 = m_base[0]->GetBdata();
1433  const Array<OneD, const NekDouble>& base1 = m_base[1]->GetBdata();
1434  const Array<OneD, const NekDouble>& dbase0 = m_base[0]->GetDbdata();
1435  const Array<OneD, const NekDouble>& dbase1 = m_base[1]->GetDbdata();
1436  const Array<OneD, const NekDouble>& metric00 = m_metrics[MetricLaplacian00];
1437  const Array<OneD, const NekDouble>& metric01 = m_metrics[MetricLaplacian01];
1438  const Array<OneD, const NekDouble>& metric11 = m_metrics[MetricLaplacian11];
1439 
1440  // Allocate temporary storage
1441  Array<OneD,NekDouble> wsp0(wsp);
1442  Array<OneD,NekDouble> wsp1(wsp+wspsize);
1443  Array<OneD,NekDouble> wsp2(wsp+2*wspsize);
1444 
1445  StdExpansion2D::PhysTensorDeriv(inarray,wsp1,wsp2);
1446 
1447  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1448  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1449  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1450  // especially for this purpose
1451  Vmath::Vvtvvtp(nqtot,&metric00[0],1,&wsp1[0],1,&metric01[0],1,&wsp2[0],1,&wsp0[0],1);
1452  Vmath::Vvtvvtp(nqtot,&metric01[0],1,&wsp1[0],1,&metric11[0],1,&wsp2[0],1,&wsp2[0],1);
1453 
1454  // outarray = m = (D_xi1 * B)^T * k
1455  // wsp1 = n = (D_xi2 * B)^T * l
1456  IProductWRTBase_SumFacKernel(dbase0,base1,wsp0,outarray,wsp1);
1457  IProductWRTBase_SumFacKernel(base0,dbase1,wsp2,wsp1, wsp0);
1458 
1459  // outarray = outarray + wsp1
1460  // = L * u_hat
1461  Vmath::Vadd(m_ncoeffs,wsp1.get(),1,outarray.get(),1,outarray.get(),1);
1462  }
1463 
1464 
1465  void TriExp::v_ComputeLaplacianMetric()
1466  {
1467  if (m_metrics.count(MetricQuadrature) == 0)
1468  {
1469  ComputeQuadratureMetric();
1470  }
1471 
1472  unsigned int i, j;
1473  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1474  const unsigned int nqtot = GetTotPoints();
1475  const unsigned int dim = 2;
1476  const MetricType m[3][3] = { {MetricLaplacian00, MetricLaplacian01, MetricLaplacian02},
1477  {MetricLaplacian01, MetricLaplacian11, MetricLaplacian12},
1478  {MetricLaplacian02, MetricLaplacian12, MetricLaplacian22}
1479  };
1480 
1481  Array<OneD, NekDouble> dEta_dXi[2] = {Array<OneD, NekDouble>(nqtot,1.0),
1482  Array<OneD, NekDouble>(nqtot,1.0)};
1483 
1484  for (i = 0; i < dim; ++i)
1485  {
1486  for (j = i; j < dim; ++j)
1487  {
1488  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1489  }
1490  }
1491 
1492  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
1493  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
1494  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1495  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1496  const Array<TwoD, const NekDouble>& df =
1497  m_metricinfo->GetDerivFactors(GetPointsKeys());
1498 
1499  for(i = 0; i < nquad1; i++)
1500  {
1501  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[0][0]+i*nquad0,1);
1502  Blas::Dscal(nquad0,2.0/(1-z1[i]),&dEta_dXi[1][0]+i*nquad0,1);
1503  }
1504  for(i = 0; i < nquad0; i++)
1505  {
1506  Blas::Dscal(nquad1,0.5*(1+z0[i]),&dEta_dXi[1][0]+i,nquad0);
1507  }
1508 
1509  Array<OneD, NekDouble> tmp(nqtot);
1510  if((type == SpatialDomains::eRegular ||
1511  type == SpatialDomains::eMovingRegular))
1512  {
1513  Vmath::Smul (nqtot,df[0][0],&dEta_dXi[0][0],1,&tmp[0],1);
1514  Vmath::Svtvp(nqtot,df[1][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1515 
1516  Vmath::Vmul (nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[MetricLaplacian00][0],1);
1517  Vmath::Smul (nqtot,df[1][0],&tmp[0],1,&m_metrics[MetricLaplacian01][0],1);
1518 
1519 
1520  Vmath::Smul (nqtot,df[2][0],&dEta_dXi[0][0],1,&tmp[0],1);
1521  Vmath::Svtvp(nqtot,df[3][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1522 
1523  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[MetricLaplacian00][0],1,&m_metrics[MetricLaplacian00][0],1);
1524  Vmath::Svtvp(nqtot,df[3][0],&tmp[0],1,&m_metrics[MetricLaplacian01][0],1,&m_metrics[MetricLaplacian01][0],1);
1525 
1526  if(GetCoordim() == 3)
1527  {
1528  Vmath::Smul (nqtot,df[4][0],&dEta_dXi[0][0],1,&tmp[0],1);
1529  Vmath::Svtvp(nqtot,df[5][0],&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1530 
1531  Vmath::Vvtvp(nqtot,&tmp[0],1, &tmp[0],1,&m_metrics[MetricLaplacian00][0],1,&m_metrics[MetricLaplacian00][0],1);
1532  Vmath::Svtvp(nqtot,df[5][0],&tmp[0],1,&m_metrics[MetricLaplacian01][0],1,&m_metrics[MetricLaplacian01][0],1);
1533  }
1534 
1535  NekDouble g2 = df[1][0]*df[1][0] + df[3][0]*df[3][0];
1536  if(GetCoordim() == 3)
1537  {
1538  g2 += df[5][0]*df[5][0];
1539  }
1540  Vmath::Fill(nqtot,g2,&m_metrics[MetricLaplacian11][0],1);
1541  }
1542  else
1543  {
1544 
1545  Vmath::Vmul (nqtot,&df[0][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1546  Vmath::Vvtvp(nqtot,&df[1][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1547 
1548  Vmath::Vmul (nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[MetricLaplacian00][0],1);
1549  Vmath::Vmul (nqtot,&df[1][0],1,&tmp[0], 1,&m_metrics[MetricLaplacian01][0],1);
1550  Vmath::Vmul (nqtot,&df[1][0],1,&df[1][0],1,&m_metrics[MetricLaplacian11][0],1);
1551 
1552 
1553  Vmath::Vmul (nqtot,&df[2][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1554  Vmath::Vvtvp(nqtot,&df[3][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1555 
1556  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[MetricLaplacian00][0],1,&m_metrics[MetricLaplacian00][0],1);
1557  Vmath::Vvtvp(nqtot,&df[3][0],1,&tmp[0], 1,&m_metrics[MetricLaplacian01][0],1,&m_metrics[MetricLaplacian01][0],1);
1558  Vmath::Vvtvp(nqtot,&df[3][0],1,&df[3][0],1,&m_metrics[MetricLaplacian11][0],1,&m_metrics[MetricLaplacian11][0],1);
1559 
1560  if(GetCoordim() == 3)
1561  {
1562  Vmath::Vmul (nqtot,&df[4][0],1,&dEta_dXi[0][0],1,&tmp[0],1);
1563  Vmath::Vvtvp(nqtot,&df[5][0],1,&dEta_dXi[1][0],1,&tmp[0],1,&tmp[0],1);
1564 
1565  Vmath::Vvtvp(nqtot,&tmp[0], 1,&tmp[0], 1,&m_metrics[MetricLaplacian00][0],1,&m_metrics[MetricLaplacian00][0],1);
1566  Vmath::Vvtvp(nqtot,&df[5][0],1,&tmp[0], 1,&m_metrics[MetricLaplacian01][0],1,&m_metrics[MetricLaplacian01][0],1);
1567  Vmath::Vvtvp(nqtot,&df[5][0],1,&df[5][0],1,&m_metrics[MetricLaplacian11][0],1,&m_metrics[MetricLaplacian11][0],1);
1568  }
1569  }
1570 
1571  for (unsigned int i = 0; i < dim; ++i)
1572  {
1573  for (unsigned int j = i; j < dim; ++j)
1574  {
1575  MultiplyByQuadratureMetric(m_metrics[m[i][j]],
1576  m_metrics[m[i][j]]);
1577 
1578  }
1579  }
1580  }
1581 
1582  /**
1583  * Function is used to compute exactly the advective numerical flux on
1584  * theinterface of two elements with different expansions, hence an
1585  * appropriate number of Gauss points has to be used. The number of
1586  * Gauss points has to be equal to the number used by the highest
1587  * polynomial degree of the two adjacent elements. Furthermore, this
1588  * function is used to compute the sensor value in each element.
1589  *
1590  * @param numMin Is the reduced polynomial order
1591  * @param inarray Input array of coefficients
1592  * @param dumpVar Output array of reduced coefficients.
1593  */
1594  void TriExp::v_ReduceOrderCoeffs(
1595  int numMin,
1596  const Array<OneD, const NekDouble> &inarray,
1597  Array<OneD, NekDouble> &outarray)
1598  {
1599  int n_coeffs = inarray.num_elements();
1600  int nquad0 = m_base[0]->GetNumPoints();
1601  int nquad1 = m_base[1]->GetNumPoints();
1602  int nqtot = nquad0*nquad1;
1603  int nmodes0 = m_base[0]->GetNumModes();
1604  int nmodes1 = m_base[1]->GetNumModes();
1605  int numMin2 = nmodes0, i;
1606 
1607  Array<OneD, NekDouble> coeff(n_coeffs,0.0);
1608  Array<OneD, NekDouble> phys_tmp(nqtot,0.0);
1609  Array<OneD, NekDouble> tmp, tmp2;
1610 
1611  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1612  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1613 
1614  LibUtilities::BasisKey b0(
1615  m_base[0]->GetBasisType(), nmodes0, Pkey0);
1616  LibUtilities::BasisKey b1(
1617  m_base[1]->GetBasisType(), nmodes1, Pkey1);
1618  LibUtilities::BasisKey bortho0(
1619  LibUtilities::eOrtho_A, nmodes0, Pkey0);
1620  LibUtilities::BasisKey bortho1(
1621  LibUtilities::eOrtho_B, nmodes1, Pkey1);
1622 
1623  // Check if it is also possible to use the same InterCoeff routine
1624  // which is also used for Quadrilateral and Hexagonal shaped
1625  // elements
1626 
1627  // For now, set up the used basis on the standard element to
1628  // calculate the phys values, set up the orthogonal basis to do a
1629  // forward transform, to obtain the coefficients in orthogonal
1630  // coefficient space
1631  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1632  StdRegions::StdTriExpSharedPtr m_TriExp;
1633 
1634  m_TriExp = MemoryManager<StdRegions::StdTriExp>
1635  ::AllocateSharedPtr(b0, b1);
1636  m_OrthoTriExp = MemoryManager<StdRegions::StdTriExp>
1637  ::AllocateSharedPtr(bortho0, bortho1);
1638 
1639  m_TriExp ->BwdTrans(inarray,phys_tmp);
1640  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1641 
1642  for (i = 0; i < n_coeffs; i++)
1643  {
1644  if (i == numMin)
1645  {
1646  coeff[i] = 0.0;
1647  numMin += numMin2 - 1;
1648  numMin2 -= 1.0;
1649  }
1650  }
1651 
1652  m_OrthoTriExp->BwdTrans(coeff,phys_tmp);
1653  m_TriExp ->FwdTrans(phys_tmp, outarray);
1654  }
1655  }
1656 }
1657 
Generated on Sun Mar 15 2015 00:11:35 for Nektar++ by   doxygen 1.8.1.2