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