Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
NodalTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalTriExp.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: NodalTriExp routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43  namespace LocalRegions
44  {
45  NodalTriExp::NodalTriExp(const LibUtilities::BasisKey &Ba,
46  const LibUtilities::BasisKey &Bb,
47  const LibUtilities::PointsType Ntype,
49  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),2,Ba,Bb),
50  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),Ba,Bb),
51  StdNodalTriExp(Ba,Bb,Ntype),
52  Expansion (geom),
53  Expansion2D (geom),
54  m_matrixManager(
55  boost::bind(&NodalTriExp::CreateMatrix, this, _1),
56  std::string("NodalTriExpMatrix")),
57  m_staticCondMatrixManager(
58  boost::bind(&NodalTriExp::CreateStaticCondMatrix, this, _1),
59  std::string("NodalTriExpStaticCondMatrix"))
60  {
61  }
62 
64  StdExpansion(T),
65  StdExpansion2D(T),
66  StdRegions::StdNodalTriExp(T),
67  Expansion (T),
68  Expansion2D (T),
69  m_matrixManager(T.m_matrixManager),
70  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
71  {
72  }
73 
75  {
76  }
77 
78  //----------------------------
79  // Integration Methods
80  //----------------------------
81 
82  /** \brief Integrate the physical point list \a inarray over region
83  and return the value
84 
85  Inputs:\n
86 
87  - \a inarray: definition of function to be returned at quadrature point
88  of expansion.
89 
90  Outputs:\n
91 
92  - returns \f$\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2) J[i,j] d
93  \xi_1 d \xi_2 \f$ where \f$inarray[i,j] = u(\xi_{1i},\xi_{2j})
94  \f$ and \f$ J[i,j] \f$ is the Jacobian evaluated at the
95  quadrature point.
96  */
97 
98 
100  {
101  int nquad0 = m_base[0]->GetNumPoints();
102  int nquad1 = m_base[1]->GetNumPoints();
104  NekDouble ival;
105  Array<OneD,NekDouble> tmp(nquad0*nquad1);
106 
107  // multiply inarray with Jacobian
108  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
109  {
110  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
111  }
112  else
113  {
114  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
115  }
116 
117  // call StdQuadExp version;
118  ival = StdNodalTriExp::v_Integral(tmp);
119  return ival;
120  }
121 
122 
124  Array<OneD, NekDouble> &outarray,
125  bool multiplybyweights)
126  {
127  int nquad0 = m_base[0]->GetNumPoints();
128  int nquad1 = m_base[1]->GetNumPoints();
129  int order1 = m_base[1]->GetNumModes();
130 
131  if(multiplybyweights)
132  {
133  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad0*order1);
134  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
135 
136  MultiplyByQuadratureMetric(inarray,tmp);
137  StdTriExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),tmp,outarray,wsp);
138  NodalToModalTranspose(outarray,outarray);
139  }
140  else
141  {
142  Array<OneD,NekDouble> wsp(nquad0*order1);
143 
144  StdTriExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),inarray,outarray,wsp);
145  NodalToModalTranspose(outarray,outarray);
146  }
147  }
148 
150  Array<OneD, NekDouble> &outarray)
151  {
152  int nq = GetTotPoints();
154  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
155 
156  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
157  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
158  }
159 
161  const Array<OneD, const NekDouble>& inarray,
162  Array<OneD, NekDouble> & outarray)
163  {
164  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
165  ASSERTL1((dir==2)?(m_geom->GetCoordim()==3):true,"Invalid direction.");
166 
167  int i;
168  int nquad0 = m_base[0]->GetNumPoints();
169  int nquad1 = m_base[1]->GetNumPoints();
170  int nqtot = nquad0*nquad1;
171  int wspsize = max(nqtot,m_ncoeffs);
172 
173  const Array<TwoD, const NekDouble>& df =
174  m_metricinfo->GetDerivFactors(GetPointsKeys());
175 
176  Array<OneD, NekDouble> tmp0 (6*wspsize);
177  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
178  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
179  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
180  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
181  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
182 
183  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
184  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
185 
186  // set up geometric factor: 2/(1-z1)
187  for(i = 0; i < nquad1; ++i)
188  {
189  gfac0[i] = 2.0/(1-z1[i]);
190  }
191  for(i = 0; i < nquad0; ++i)
192  {
193  gfac1[i] = 0.5*(1+z0[i]);
194  }
195 
196  for(i = 0; i < nquad1; ++i)
197  {
198  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
199  }
200 
201  for(i = 0; i < nquad1; ++i)
202  {
203  Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
204  }
205 
206  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
207  {
208  Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
209  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
210  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
211  }
212  else
213  {
214  Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
215  Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
216  Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
217  }
218  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
219 
220  MultiplyByQuadratureMetric(tmp1,tmp1);
221  MultiplyByQuadratureMetric(tmp2,tmp2);
222 
223  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),m_base[1]->GetBdata() ,tmp1,tmp3 ,tmp0);
224  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata() ,m_base[1]->GetDbdata(),tmp2,outarray,tmp0);
225  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
226 
227  NodalToModalTranspose(outarray,outarray);
228  }
229 
231  const Array<OneD, const NekDouble>& inarray,
232  Array<OneD, NekDouble> &outarray)
233  {
234  int nq = GetTotPoints();
236 
237  switch(dir)
238  {
239  case 0:
240  {
242  }
243  break;
244  case 1:
245  {
247  }
248  break;
249  case 2:
250  {
252  }
253  break;
254  default:
255  {
256  ASSERTL1(false,"input dir is out of range");
257  }
258  break;
259  }
260 
261  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
262  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
263 
264  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
265  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
266 
267  }
268 
269  ///////////////////////////////
270  /// Differentiation Methods
271  ///////////////////////////////
272 
273  /**
274  \brief Calculate the deritive of the physical points
275  **/
277  Array<OneD,NekDouble> &out_d0,
278  Array<OneD,NekDouble> &out_d1,
279  Array<OneD,NekDouble> &out_d2)
280  {
281  int nquad0 = m_base[0]->GetNumPoints();
282  int nquad1 = m_base[1]->GetNumPoints();
283  int nqtot = nquad0*nquad1;
285  = m_metricinfo->GetDerivFactors(GetPointsKeys());
286 
287  Array<OneD,NekDouble> diff0(2*nqtot);
288  Array<OneD,NekDouble> diff1(diff0+nqtot);
289 
290  StdNodalTriExp::v_PhysDeriv(inarray, diff0, diff1);
291 
292  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
293  {
294  if(out_d0.num_elements())
295  {
296  Vmath::Vmul (nqtot,df[0],1,diff0,1, out_d0, 1);
297  Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
298  }
299 
300  if(out_d1.num_elements())
301  {
302  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
303  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
304  }
305 
306  if(out_d2.num_elements())
307  {
308  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
309  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
310  }
311  }
312  else // regular geometry
313  {
314  if(out_d0.num_elements())
315  {
316  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
317  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
318  }
319 
320  if(out_d1.num_elements())
321  {
322  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
323  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
324  }
325 
326  if(out_d2.num_elements())
327  {
328  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
329  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
330  }
331  }
332  }
333 
334  /** \brief Forward transform from physical quadrature space
335  stored in \a inarray and evaluate the expansion coefficients and
336  store in \a (this)->m_coeffs
337 
338  Inputs:\n
339 
340  - \a inarray: array of physical quadrature points to be transformed
341 
342  Outputs:\n
343 
344  - (this)->_coeffs: updated array of expansion coefficients.
345 
346  */
348  Array<OneD,NekDouble> &outarray)
349  {
350  IProductWRTBase(inarray,outarray);
351 
352  // get Mass matrix inverse
354  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
355 
356  // copy inarray in case inarray == outarray
357  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
359 
360  out = (*matsys)*in;
361  }
362 
364  Array<OneD,NekDouble> &outarray,
365  const StdRegions::StdMatrixKey &mkey)
366  {
368 
369  if(inarray.get() == outarray.get())
370  {
372  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
373 
374  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
375  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
376  }
377  else
378  {
379  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
380  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
381  }
382  }
383 
385  Array<OneD,NekDouble> &coords_1,
386  Array<OneD,NekDouble> &coords_2)
387  {
388  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
389  }
390 
391  // get the coordinates "coords" at the local coordinates "Lcoords"
393  Array<OneD,NekDouble> &coords)
394  {
395  int i;
396 
397  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
398  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
399  "Local coordinates are not in region [-1,1]");
400 
401  m_geom->FillGeom();
402 
403  for(i = 0; i < m_geom->GetCoordim(); ++i)
404  {
405  coords[i] = m_geom->GetCoord(i,Lcoords);
406  }
407  }
408 
410  {
411  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
412  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
415  AllocateSharedPtr(bkey0,bkey1,ntype);
416 
417  return tmp->GetStdMatrix(mkey);
418  }
419 
421  const Array<OneD, const NekDouble> &coord,
422  const Array<OneD, const NekDouble> &physvals)
423 
424  {
426 
427  ASSERTL0(m_geom,"m_geom not defined");
428  m_geom->GetLocCoords(coord,Lcoord);
429 
430  return StdNodalTriExp::v_PhysEvaluate(Lcoord, physvals);
431  }
432 
434  {
435  DNekScalMatSharedPtr returnval;
437 
438  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
439 
440  StdRegions::MatrixType mtype = mkey.GetMatrixType();
441 
442  switch(mtype)
443  {
444  case StdRegions::eMass:
445  {
446  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
447  {
448  NekDouble one = 1.0;
449  DNekMatSharedPtr mat = GenMatrix(mkey);
451  }
452  else
453  {
454  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
455  DNekMatSharedPtr mat = GetStdMatrix(mkey);
457  }
458  }
459  break;
461  {
462  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
463  {
464  NekDouble one = 1.0;
466  *this);
467  DNekMatSharedPtr mat = GenMatrix(masskey);
468  mat->Invert();
469 
471  }
472  else
473  {
474  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
475  DNekMatSharedPtr mat = GetStdMatrix(mkey);
477  }
478  }
479  break;
481  {
482  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
483  {
484  NekDouble one = 1.0;
485  DNekMatSharedPtr mat = GenMatrix(mkey);
486 
488  }
489  else
490  {
491  ASSERTL1(m_geom->GetCoordim() == 2,"Standard Region Laplacian is only set up for Quads in two-dimensional");
493  mkey.GetShapeType(), *this);
495  mkey.GetShapeType(), *this);
497  mkey.GetShapeType(), *this);
498 
499  DNekMat &lap00 = *GetStdMatrix(lap00key);
500  DNekMat &lap01 = *GetStdMatrix(lap01key);
501  DNekMat &lap11 = *GetStdMatrix(lap11key);
502 
503  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
505  m_metricinfo->GetGmat(ptsKeys);
506 
507  int rows = lap00.GetRows();
508  int cols = lap00.GetColumns();
509 
511 
512  (*lap) = gmat[0][0] * lap00 +
513  gmat[1][0] * (lap01 + Transpose(lap01)) +
514  gmat[3][0] * lap11;
515 
517  }
518  }
519  break;
521  {
524  mkey.GetShapeType(), *this);
525  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
527  mkey.GetShapeType(), *this);
528  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
529 
530  int rows = LapMat.GetRows();
531  int cols = LapMat.GetColumns();
532 
534 
535  NekDouble one = 1.0;
536  (*helm) = LapMat + factor*MassMat;
537 
538  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
539  }
540  break;
541  default:
542  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
543  break;
544  }
545 
546  return returnval;
547  }
548 
550  {
551  DNekScalBlkMatSharedPtr returnval;
552 
553  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
554 
555  // set up block matrix system
556  unsigned int nbdry = NumBndryCoeffs();
557  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
558  unsigned int exp_size[] = {nbdry,nint};
559  unsigned int nblks = 2;
560  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
561  NekDouble factor = 1.0;
562 
563  switch(mkey.GetMatrixType())
564  {
566  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
567 
568  // use Deformed case for both regular and deformed geometries
569  factor = 1.0;
570  goto UseLocRegionsMatrix;
571  break;
572  default:
573  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
574  {
575  factor = 1.0;
576  goto UseLocRegionsMatrix;
577  }
578  else
579  {
581  factor = mat->Scale();
582  goto UseStdRegionsMatrix;
583  }
584  break;
585  UseStdRegionsMatrix:
586  {
587  NekDouble invfactor = 1.0/factor;
588  NekDouble one = 1.0;
591  DNekMatSharedPtr Asubmat;
592 
593  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
594  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
595  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
596  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
597  }
598  break;
599  UseLocRegionsMatrix:
600  {
601  int i,j;
602  NekDouble invfactor = 1.0/factor;
603  NekDouble one = 1.0;
604  DNekScalMat &mat = *GetLocMatrix(mkey);
609 
610  Array<OneD,unsigned int> bmap(nbdry);
611  Array<OneD,unsigned int> imap(nint);
612  GetBoundaryMap(bmap);
613  GetInteriorMap(imap);
614 
615  for(i = 0; i < nbdry; ++i)
616  {
617  for(j = 0; j < nbdry; ++j)
618  {
619  (*A)(i,j) = mat(bmap[i],bmap[j]);
620  }
621 
622  for(j = 0; j < nint; ++j)
623  {
624  (*B)(i,j) = mat(bmap[i],imap[j]);
625  }
626  }
627 
628  for(i = 0; i < nint; ++i)
629  {
630  for(j = 0; j < nbdry; ++j)
631  {
632  (*C)(i,j) = mat(imap[i],bmap[j]);
633  }
634 
635  for(j = 0; j < nint; ++j)
636  {
637  (*D)(i,j) = mat(imap[i],imap[j]);
638  }
639  }
640 
641  // Calculate static condensed system
642  if(nint)
643  {
644  D->Invert();
645  (*B) = (*B)*(*D);
646  (*A) = (*A) - (*B)*(*C);
647  }
648 
650 
651  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
652  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
653  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
654  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
655 
656  }
657  }
658 
659  return returnval;
660  }
661 
662 
664  {
665 
667  ::AllocateSharedPtr(m_base[0]->GetBasisKey(),
668  m_base[1]->GetBasisKey(),
670  }
671 
673  {
675  2, m_base[0]->GetPointsKey());
677  2, m_base[1]->GetPointsKey());
678 
681  }
682 
684  {
685  DNekMatSharedPtr returnval;
686 
687  switch(mkey.GetMatrixType())
688  {
695  returnval = Expansion2D::v_GenMatrix(mkey);
696  break;
697  default:
698  returnval = StdNodalTriExp::v_GenMatrix(mkey);
699  break;
700  }
701  return returnval;
702  }
703 
704  void NodalTriExp::v_ComputeEdgeNormal(const int edge)
705  {
706  int i;
707  const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
708  const SpatialDomains::GeomType type = geomFactors->GetGtype();
709 
711  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
712  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
713  int nqe = m_base[0]->GetNumPoints();
714  int dim = GetCoordim();
715 
718  for (i = 0; i < dim; ++i)
719  {
720  normal[i] = Array<OneD, NekDouble>(nqe);
721  }
722 
723  // Regular geometry case
725  {
726  NekDouble fac;
727  // Set up normals
728  switch(edge)
729  {
730  case 0:
731  for(i = 0; i < GetCoordim(); ++i)
732  {
733  Vmath::Fill(nqe,-df[2*i+1][0],normal[i],1);
734  }
735  break;
736  case 1:
737  for(i = 0; i < GetCoordim(); ++i)
738  {
739  Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
740  }
741  break;
742  case 2:
743  for(i = 0; i < GetCoordim(); ++i)
744  {
745  Vmath::Fill(nqe,-df[2*i][0],normal[i],1);
746  }
747  break;
748  default:
749  ASSERTL0(false,"Edge is out of range (edge < 3)");
750  }
751 
752  // normalise
753  fac = 0.0;
754  for(i =0 ; i < GetCoordim(); ++i)
755  {
756  fac += normal[i][0]*normal[i][0];
757  }
758  fac = 1.0/sqrt(fac);
759  for (i = 0; i < GetCoordim(); ++i)
760  {
761  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
762  }
763  }
764  else // Set up deformed normals
765  {
766  int j;
767 
768  int nquad0 = ptsKeys[0].GetNumPoints();
769  int nquad1 = ptsKeys[1].GetNumPoints();
770 
771  LibUtilities::PointsKey from_key;
772 
773  Array<OneD,NekDouble> normals(GetCoordim()*max(nquad0,nquad1),0.0);
774  Array<OneD,NekDouble> edgejac(GetCoordim()*max(nquad0,nquad1),0.0);
775 
776  // Extract Jacobian along edges and recover local
777  // derivates (dx/dr) for polynomial interpolation by
778  // multiplying m_gmat by jacobian
779  switch(edge)
780  {
781  case 0:
782  for(j = 0; j < nquad0; ++j)
783  {
784  edgejac[j] = jac[j];
785  for(i = 0; i < GetCoordim(); ++i)
786  {
787  normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
788  }
789  }
790  from_key = ptsKeys[0];
791  break;
792  case 1:
793  for(j = 0; j < nquad1; ++j)
794  {
795  edgejac[j] = jac[nquad0*j+nquad0-1];
796  for(i = 0; i < GetCoordim(); ++i)
797  {
798  normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
799  }
800  }
801  from_key = ptsKeys[1];
802  break;
803  case 2:
804  for(j = 0; j < nquad1; ++j)
805  {
806  edgejac[j] = jac[nquad0*j];
807  for(i = 0; i < GetCoordim(); ++i)
808  {
809  normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
810  }
811  }
812  from_key = ptsKeys[1];
813  break;
814  default:
815  ASSERTL0(false,"edge is out of range (edge < 3)");
816 
817  }
818 
819  int nq = from_key.GetNumPoints();
820  Array<OneD,NekDouble> work(nqe,0.0);
821 
822  // interpolate Jacobian and invert
823  LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
824  Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
825 
826  // interpolate
827  for(i = 0; i < GetCoordim(); ++i)
828  {
829  LibUtilities::Interp1D(from_key,&normals[i*nq],m_base[0]->GetPointsKey(),&normal[i][0]);
830  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
831  }
832 
833  //normalise normal vectors
834  Vmath::Zero(nqe,work,1);
835  for(i = 0; i < GetCoordim(); ++i)
836  {
837  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
838  }
839 
840  Vmath::Vsqrt(nqe,work,1,work,1);
841  Vmath::Sdiv(nqe,1.0,work,1,work,1);
842 
843  for(i = 0; i < GetCoordim(); ++i)
844  {
845  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
846  }
847 
848  // Reverse direction so that points are in
849  // anticlockwise direction if edge >=2
850  if(edge >= 2)
851  {
852  for(i = 0; i < GetCoordim(); ++i)
853  {
854  Vmath::Reverse(nqe,normal[i],1, normal[i],1);
855  }
856  }
857  }
858  }
859 
860  }//end of namespace
861 }//end of namespace
const LibUtilities::PointsKeyVector GetPointsKeys() const
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
NodalTriExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::PointsType Ntype, const SpatialDomains::TriGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
Definition: NodalTriExp.cpp:45
void v_ComputeEdgeNormal(const int edge)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:947
void GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
void IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Differentiation Methods.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
std::map< int, StdRegions::NormalVector > m_edgeNormals
Definition: Expansion2D.h:135
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:131
STL namespace.
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:271
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:130
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
boost::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
void NodalToModalTranspose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: NodalTriExp.h:187
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:711
LibUtilities::PointsKey m_nodalPointsKey
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1085
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
void IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:224
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Defines a specification for a set of points.
Definition: Points.h:58
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
double NekDouble
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
NekDouble Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region.
Definition: NodalTriExp.cpp:99
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:85
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:161
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Geometry is straight-sided with constant geometric factors.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:54
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
unsigned int GetNumPoints() const
Definition: Points.h:106
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
GeomType
Indicates the type of element geometry.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Geometry is curved or has non-constant factors.
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Inner product of inarray over region with respect to the expansion basis (this)->_Base[0] and return ...
Definition: NodalTriExp.h:83
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:814
Describes the specification for a Basis.
Definition: Basis.h:50
PointsType GetPointsType() const
Definition: Points.h:111
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
DNekMatSharedPtr CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:228
static ConstFactorMap NullConstFactorMap
Definition: StdRegions.hpp:253