Nektar++
SegExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File SegExp.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: SegExp routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LocalRegions/SegExp.h>
39 
40 
41 namespace Nektar
42 {
43  namespace LocalRegions
44  {
45 
46  /**
47  * @class SegExp
48  * Defines a Segment local expansion.
49  */
50 
51  /// Constructor using BasisKey class for quadrature points and
52  /// order definition.
53  /**
54  * @param Ba Basis key of segment expansion.
55  * @param geom Description of geometry.
56  */
59  StdExpansion(Ba.GetNumModes(), 1, Ba),
60  StdExpansion1D(Ba.GetNumModes(), Ba),
61  StdRegions::StdSegExp(Ba),
62  Expansion(geom),
63  Expansion1D(geom),
64  m_matrixManager(
65  boost::bind(&SegExp::CreateMatrix, this, _1),
66  std::string("SegExpMatrix")),
67  m_staticCondMatrixManager(
68  boost::bind(&SegExp::CreateStaticCondMatrix, this, _1),
69  std::string("SegExpStaticCondMatrix"))
70  {
71  }
72 
73 
74  /// Copy Constructor
75  /**
76  * @param S Existing segment to duplicate.
77  */
79  StdExpansion(S),
80  StdExpansion1D(S),
81  StdRegions::StdSegExp(S),
82  Expansion(S),
83  Expansion1D(S),
84  m_matrixManager(S.m_matrixManager),
85  m_staticCondMatrixManager(S.m_staticCondMatrixManager)
86  {
87  }
88 
89 
90  /**
91  *
92  */
94  {
95  }
96 
97 
98  //----------------------------
99  // Integration Methods
100  //----------------------------
101 
102  /** \brief Integrate the physical point list \a inarray over region
103  and return the value
104 
105  Inputs:\n
106 
107  - \a inarray: definition of function to be returned at
108  quadrature point of expansion.
109 
110  Outputs:\n
111 
112  - returns \f$\int^1_{-1} u(\xi_1)d \xi_1 \f$ where \f$inarray[i]
113  = u(\xi_{1i}) \f$
114  */
115 
117  const Array<OneD, const NekDouble>& inarray)
118  {
119  int nquad0 = m_base[0]->GetNumPoints();
121  NekDouble ival;
122  Array<OneD,NekDouble> tmp(nquad0);
123 
124  // multiply inarray with Jacobian
125  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
126  {
127  Vmath::Vmul(nquad0, jac, 1, inarray, 1, tmp,1);
128  }
129  else
130  {
131  Vmath::Smul(nquad0, jac[0], inarray, 1, tmp, 1);
132  }
133 
134  // call StdSegExp version;
135  ival = StdSegExp::v_Integral(tmp);
136  //ival = StdSegExp::Integral(tmp);
137  return ival;
138  }
139 
140  //-----------------------------
141  // Differentiation Methods
142  //-----------------------------
143 
144  /** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the
145  physical quadrature points given by \a inarray and return in \a
146  outarray.
147 
148  This is a wrapper around StdExpansion1D::Tensor_Deriv
149 
150  Input:\n
151 
152  - \a n: number of derivatives to be evaluated where \f$ n \leq dim\f$
153 
154  - \a inarray: array of function evaluated at the quadrature points
155 
156  Output: \n
157 
158  - \a outarray: array of the derivatives \f$
159  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dx,
160  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dy,
161  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dz,
162  \f$ depending on value of \a dim
163  */
165  const Array<OneD, const NekDouble>& inarray,
166  Array<OneD,NekDouble> &out_d0,
167  Array<OneD,NekDouble> &out_d1,
168  Array<OneD,NekDouble> &out_d2)
169  {
170  int nquad0 = m_base[0]->GetNumPoints();
172  m_metricinfo->GetDerivFactors(GetPointsKeys());
173  Array<OneD,NekDouble> diff(nquad0);
174 
175  //StdExpansion1D::PhysTensorDeriv(inarray,diff);
176  PhysTensorDeriv(inarray,diff);
177  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
178  {
179  if(out_d0.num_elements())
180  {
181  Vmath::Vmul(nquad0,&gmat[0][0],1,&diff[0],1,
182  &out_d0[0],1);
183  }
184 
185  if(out_d1.num_elements())
186  {
187  Vmath::Vmul(nquad0,&gmat[1][0],1,&diff[0],1,
188  &out_d1[0],1);
189  }
190 
191  if(out_d2.num_elements())
192  {
193  Vmath::Vmul(nquad0,&gmat[2][0],1,&diff[0],1,
194  &out_d2[0],1);
195  }
196  }
197  else
198  {
199  if(out_d0.num_elements())
200  {
201  Vmath::Smul(nquad0, gmat[0][0], diff, 1,
202  out_d0, 1);
203  }
204 
205  if(out_d1.num_elements())
206  {
207  Vmath::Smul(nquad0, gmat[1][0], diff, 1,
208  out_d1, 1);
209  }
210 
211  if(out_d2.num_elements())
212  {
213  Vmath::Smul(nquad0, gmat[2][0], diff, 1,
214  out_d2, 1);
215  }
216  }
217  }
218 
219  /**
220  *\brief Evaluate the derivative along a line:
221  * \f$ d/ds=\frac{spacedim}{||tangent||}d/d{\xi} \f$.
222  * The derivative is calculated performing
223  *the product \f$ du/d{s}=\nabla u \cdot tangent \f$.
224  *\param inarray function to derive
225  *\param out_ds result of the derivative operation
226  **/
228  const Array<OneD, const NekDouble>& inarray,
229  Array<OneD,NekDouble> &out_ds)
230  {
231  int nquad0 = m_base[0]->GetNumPoints();
232  int coordim = m_geom->GetCoordim();
233  Array<OneD, NekDouble> diff (nquad0);
234  //this operation is needed if you put out_ds==inarray
235  Vmath::Zero(nquad0,out_ds,1);
236  switch(coordim)
237  {
238  case 2:
239  //diff= dU/de
240  Array<OneD,NekDouble> diff(nquad0);
241 
242  PhysTensorDeriv(inarray,diff);
243 
244  //get dS/de= (Jac)^-1
246  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
247  {
248  //calculate the derivative as (dU/de)*(Jac)^-1
249  Vmath::Vdiv(nquad0,diff,1,Jac ,1,out_ds,1);
250  }
251  else
252  {
253  NekDouble invJac = 1/Jac[0];
254  Vmath::Smul(nquad0, invJac,diff,1,out_ds,1);
255  }
256  }
257  }
258 
259  /**
260  *\brief Evaluate the derivative normal to a line:
261  * \f$ d/dn=\frac{spacedim}{||normal||}d/d{\xi} \f$.
262  * The derivative is calculated performing
263  *the product \f$ du/d{s}=\nabla u \cdot normal \f$.
264  *\param inarray function to derive
265  *\param out_dn result of the derivative operation
266  **/
268  const Array<OneD, const NekDouble>& inarray,
269  Array<OneD, NekDouble>& out_dn)
270  {
271  int nquad0 = m_base[0]->GetNumPoints();
273  m_metricinfo->GetDerivFactors(GetPointsKeys());
274  int coordim = m_geom->GetCoordim();
275  Array<OneD, NekDouble> out_dn_tmp(nquad0,0.0);
276  switch(coordim)
277  {
278  case 2:
279 
280  Array<OneD, NekDouble> inarray_d0(nquad0);
281  Array<OneD, NekDouble> inarray_d1(nquad0);
282 
283  v_PhysDeriv(inarray,inarray_d0,inarray_d1);
285  normals = Array<OneD, Array<OneD, NekDouble> >(coordim);
286  cout<<"der_n"<<endl;
287  for(int k=0; k<coordim; ++k)
288  {
289  normals[k]= Array<OneD, NekDouble>(nquad0);
290  }
291 // @TODO: this routine no longer makes sense, since normals are not unique on
292 // an edge
293 // normals = GetMetricInfo()->GetNormal();
294  for(int i=0; i<nquad0; i++)
295  {
296 cout<<"nx= "<<normals[0][i]<<" ny="<<normals[1][i]<<endl;
297  }
299  "normal vectors do not exist: check if a"
300  "boundary region is defined as I ");
301  // \nabla u \cdot normal
302  Vmath::Vmul(nquad0,normals[0],1,inarray_d0,1,out_dn_tmp,1);
303  Vmath::Vadd(nquad0,out_dn_tmp,1,out_dn,1,out_dn,1);
304  Vmath::Zero(nquad0,out_dn_tmp,1);
305  Vmath::Vmul(nquad0,normals[1],1,inarray_d1,1,out_dn_tmp,1);
306  Vmath::Vadd(nquad0,out_dn_tmp,1,out_dn,1,out_dn,1);
307 
308  for(int i=0; i<nquad0; i++)
309  {
310 cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
311  }
312 
313 
314  }
315 
316  }
317  void SegExp::v_PhysDeriv(const int dir,
318  const Array<OneD, const NekDouble>& inarray,
319  Array<OneD, NekDouble> &outarray)
320  {
321  switch(dir)
322  {
323  case 0:
324  {
325  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
327  }
328  break;
329  case 1:
330  {
331  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
333  }
334  break;
335  case 2:
336  {
338  NullNekDouble1DArray, outarray);
339  }
340  break;
341  default:
342  {
343  ASSERTL1(false,"input dir is out of range");
344  }
345  break;
346  }
347  }
348 
349 
350  //-----------------------------
351  // Transforms
352  //-----------------------------
353 
354  /** \brief Forward transform from physical quadrature space
355  stored in \a inarray and evaluate the expansion coefficients and
356  store in \a outarray
357 
358  Perform a forward transform using a Galerkin projection by
359  taking the inner product of the physical points and multiplying
360  by the inverse of the mass matrix using the Solve method of the
361  standard matrix container holding the local mass matrix, i.e.
362  \f$ {\bf \hat{u}} = {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] =
363  \int^1_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
364 
365  Inputs:\n
366 
367  - \a inarray: array of physical quadrature points to be transformed
368 
369  Outputs:\n
370 
371  - \a outarray: updated array of expansion coefficients.
372 
373  */
374  // need to sort out family of matrices
376  const Array<OneD, const NekDouble>& inarray,
377  Array<OneD,NekDouble> &outarray)
378  {
379  if (m_base[0]->Collocation())
380  {
381  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
382  }
383  else
384  {
385  v_IProductWRTBase(inarray,outarray);
386 
387  // get Mass matrix inverse
388  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
389  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
390 
391  // copy inarray in case inarray == outarray
392  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
394 
395  out = (*matsys)*in;
396  }
397  }
398 
400  const Array<OneD, const NekDouble>& inarray,
401  Array<OneD, NekDouble> &outarray)
402  {
403  if(m_base[0]->Collocation())
404  {
405  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
406  }
407  else
408  {
409  int nInteriorDofs = m_ncoeffs-2;
410  int offset;
411 
412  switch (m_base[0]->GetBasisType())
413  {
415  {
416  offset = 1;
417  }
418  break;
420  {
421  nInteriorDofs = m_ncoeffs;
422  offset = 0;
423  }
424  break;
427  {
428  offset = 2;
429  }
430  break;
431  default:
432  ASSERTL0(false,"This type of FwdTrans is not defined"
433  "for this expansion type");
434  }
435 
436  fill(outarray.get(), outarray.get()+m_ncoeffs, 0.0 );
437 
439  {
440 
441  outarray[GetVertexMap(0)] = inarray[0];
442  outarray[GetVertexMap(1)] =
443  inarray[m_base[0]->GetNumPoints()-1];
444 
445  if (m_ncoeffs>2)
446  {
447  // ideally, we would like to have tmp0 to be replaced
448  // by outarray (currently MassMatrixOp does not allow
449  // aliasing)
452 
453  StdRegions::StdMatrixKey stdmasskey(
455  MassMatrixOp(outarray,tmp0,stdmasskey);
456  v_IProductWRTBase(inarray,tmp1);
457 
458  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
459 
460  // get Mass matrix inverse (only of interior DOF)
461  MatrixKey masskey(
463  DNekScalMatSharedPtr matsys =
464  (m_staticCondMatrixManager[masskey])->GetBlock(1,1);
465 
466  Blas::Dgemv('N',nInteriorDofs,nInteriorDofs,
467  matsys->Scale(),
468  &((matsys->GetOwnedMatrix())->GetPtr())[0],
469  nInteriorDofs,tmp1.get()+offset,1,0.0,
470  outarray.get()+offset,1);
471  }
472  }
473  else
474  {
475  SegExp::v_FwdTrans(inarray, outarray);
476 
477  }
478  }
479  }
480 
481 
482  //-----------------------------
483  // Inner product functions
484  //-----------------------------
485 
486  /** \brief Inner product of \a inarray over region with respect to
487  the expansion basis (this)->_Base[0] and return in \a outarray
488 
489  Wrapper call to SegExp::IProduct_WRT_B
490 
491  Input:\n
492 
493  - \a inarray: array of function evaluated at the physical
494  collocation points
495 
496  Output:\n
497 
498  - \a outarray: array of inner product with respect to each
499  basis over region
500  */
502  const Array<OneD, const NekDouble>& inarray,
503  Array<OneD, NekDouble> &outarray)
504  {
505  v_IProductWRTBase(m_base[0]->GetBdata(),inarray,outarray,1);
506  }
507 
508 
509  /**
510  \brief Inner product of \a inarray over region with respect to
511  expansion basis \a base and return in \a outarray
512 
513  Calculate \f$ I[p] = \int^{1}_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1
514  = \sum_{i=0}^{nq-1} \phi_p(\xi_{1i}) u(\xi_{1i}) w_i \f$ where
515  \f$ outarray[p] = I[p], inarray[i] = u(\xi_{1i}), base[p*nq+i] =
516  \phi_p(\xi_{1i}) \f$.
517 
518  Inputs: \n
519 
520  - \a base: an array definiing the local basis for the inner
521  product usually passed from Basis->get_bdata() or
522  Basis->get_Dbdata()
523  - \a inarray: physical point array of function to be integrated
524  \f$ u(\xi_1) \f$
525  - \a coll_check: Flag to identify when a Basis->collocation()
526  call should be performed to see if this is a GLL_Lagrange basis
527  with a collocation property. (should be set to 0 if taking the
528  inner product with respect to the derivative of basis)
529 
530  Output: \n
531 
532  - \a outarray: array of coefficients representing the inner
533  product of function with ever mode in the exapnsion
534 
535  **/
537  const Array<OneD, const NekDouble>& base,
538  const Array<OneD, const NekDouble>& inarray,
539  Array<OneD, NekDouble> &outarray,
540  int coll_check)
541  {
542  int nquad0 = m_base[0]->GetNumPoints();
544  Array<OneD,NekDouble> tmp(nquad0);
545 
546 
547  // multiply inarray with Jacobian
548  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
549  {
550  Vmath::Vmul(nquad0, jac, 1, inarray, 1, tmp, 1);
551  }
552  else
553  {
554  Vmath::Smul(nquad0, jac[0], inarray, 1, tmp, 1);
555  }
556  StdSegExp::v_IProductWRTBase(base,tmp,outarray,coll_check);
557  }
558 
559 
561  const int dir,
562  const Array<OneD, const NekDouble>& inarray,
563  Array<OneD, NekDouble> & outarray)
564  {
565  int nquad = m_base[0]->GetNumPoints();
566  const Array<TwoD, const NekDouble>& gmat =
567  m_metricinfo->GetDerivFactors(GetPointsKeys());
568 
569  Array<OneD, NekDouble> tmp1(nquad);
570 
571  switch(dir)
572  {
573  case 0:
574  {
575  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
576  {
577  Vmath::Vmul(nquad,gmat[0],1,inarray,1,tmp1,1);
578  }
579  else
580  {
581  Vmath::Smul(nquad, gmat[0][0], inarray, 1, tmp1, 1);
582  }
583  }
584  break;
585  case 1:
586  {
587  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
588  {
589  Vmath::Vmul(nquad,gmat[1],1,inarray,1,tmp1,1);
590  }
591  else
592  {
593  Vmath::Smul(nquad, gmat[1][0], inarray, 1, tmp1, 1);
594  }
595  }
596  break;
597  case 2:
598  {
599  ASSERTL1(m_geom->GetCoordim() == 3,"input dir is out of range");
600  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
601  {
602  Vmath::Vmul(nquad,gmat[2],1,inarray,1,tmp1,1);
603  }
604  else
605  {
606  Vmath::Smul(nquad, gmat[2][0], inarray, 1, tmp1, 1);
607  }
608  }
609  break;
610  default:
611  {
612  ASSERTL1(false,"input dir is out of range");
613  }
614  break;
615  }
616  v_IProductWRTBase(m_base[0]->GetDbdata(),tmp1,outarray,1);
617  }
618 
622  Array<OneD, NekDouble> &outarray)
623  {
624  int nq = m_base[0]->GetNumPoints();
626 // cout << "I am segment " << GetGeom()->GetGlobalID() << endl;
627 // cout << "I want edge " << GetLeftAdjacentElementEdge() << endl;
628 // @TODO: This routine no longer makes sense as a normal is not unique to an edge
630  &normals =
633  Vmath::Vmul (nq, &Fx[0], 1, &normals[0][0], 1, &Fn[0], 1);
634  Vmath::Vvtvp(nq, &Fy[0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
635 
636  v_IProductWRTBase(Fn,outarray);
637  }
638 
639  //-----------------------------
640  // Evaluation functions
641  //-----------------------------
642 
643 
644  /**
645  * Given the local cartesian coordinate \a Lcoord evaluate the
646  * value of physvals at this point by calling through to the
647  * StdExpansion method
648  */
650  const Array<OneD, const NekDouble> &Lcoord,
651  const Array<OneD, const NekDouble> &physvals)
652  {
653  // Evaluate point in local (eta) coordinates.
654  return StdSegExp::v_PhysEvaluate(Lcoord,physvals);
655  }
656 
658  const Array<OneD, const NekDouble>& coord,
659  const Array<OneD, const NekDouble> &physvals)
660  {
662 
663  ASSERTL0(m_geom,"m_geom not defined");
664  m_geom->GetLocCoords(coord,Lcoord);
665 
666  return StdSegExp::v_PhysEvaluate(Lcoord, physvals);
667  }
668 
669 
671  const Array<OneD, const NekDouble>& Lcoords,
672  Array<OneD,NekDouble> &coords)
673  {
674  int i;
675 
676  ASSERTL1(Lcoords[0] >= -1.0&& Lcoords[0] <= 1.0,
677  "Local coordinates are not in region [-1,1]");
678 
679  m_geom->FillGeom();
680  for(i = 0; i < m_geom->GetCoordim(); ++i)
681  {
682  coords[i] = m_geom->GetCoord(i,Lcoords);
683  }
684  }
685 
687  Array<OneD, NekDouble> &coords_0,
688  Array<OneD, NekDouble> &coords_1,
689  Array<OneD, NekDouble> &coords_2)
690  {
691  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
692  }
693 
694  // Get vertex value from the 1D Phys space.
696  const int vertex,
697  const Array<OneD, const NekDouble> &inarray,
698  NekDouble &outarray)
699  {
700  int nquad = m_base[0]->GetNumPoints();
701 
703  {
704  switch (vertex)
705  {
706  case 0:
707  outarray = inarray[0];
708  break;
709  case 1:
710  outarray = inarray[nquad - 1];
711  break;
712  }
713  }
714  else
715  {
717  factors[StdRegions::eFactorGaussVertex] = vertex;
718 
721  DetShapeType(),*this,factors);
722 
723  DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
724 
725  outarray = Blas::Ddot(nquad, mat_gauss->GetOwnedMatrix()
726  ->GetPtr().get(), 1, &inarray[0], 1);
727  }
728  }
729 
730 
731  //-----------------------------
732  // Helper functions
733  //-----------------------------
734 
736  Array<OneD, NekDouble> &coeffs,
738  {
739  v_SetCoeffsToOrientation(dir,coeffs,coeffs);
740  }
741 
745  Array<OneD, NekDouble> &outarray)
746  {
747 
748  if (dir == StdRegions::eBackwards)
749  {
750  if (&inarray[0] != &outarray[0])
751  {
752  Array<OneD,NekDouble> intmp (inarray);
753  ReverseCoeffsAndSign(intmp,outarray);
754  }
755  else
756  {
757  ReverseCoeffsAndSign(inarray,outarray);
758  }
759  }
760  }
761 
763  {
764  return m_geom->GetPorient(point);
765  }
766 
767 
769  {
771  ::AllocateSharedPtr(m_base[0]->GetBasisKey());
772  }
773 
775  {
776  return m_geom->GetCoordim();
777  }
778 
780  {
781  NEKERROR(ErrorUtil::efatal, "Got to SegExp");
782  return NullNekDouble1DArray;
783  }
784 
786  {
787  return m_metricinfo->GetGtype();
788  }
789 
790  int SegExp::v_GetNumPoints(const int dir) const
791  {
792  return GetNumPoints(dir);
793  }
794 
795  int SegExp::v_GetNcoeffs(void) const
796  {
797  return m_ncoeffs;
798  }
799 
801  {
802  return GetBasis(dir);
803  }
804 
805 
807  {
808  return 2;
809  }
810 
811 
813  {
814  return 2;
815  }
816 
817 
818  /// Unpack data from input file assuming it comes from
819  // the same expansion type
821  const NekDouble *data,
822  const std::vector<unsigned int > &nummodes,
823  const int mode_offset,
824  NekDouble *coeffs)
825  {
826  switch(m_base[0]->GetBasisType())
827  {
829  {
830  int fillorder = min((int) nummodes[mode_offset],m_ncoeffs);
831 
832  Vmath::Zero(m_ncoeffs,coeffs,1);
833  Vmath::Vcopy(fillorder,&data[0],1,&coeffs[0],1);
834  }
835  break;
837  {
838  // Assume that input is also Gll_Lagrange
839  // but no way to check;
841  nummodes[mode_offset],
844  p0,data, m_base[0]->GetPointsKey(), coeffs);
845  }
846  break;
848  {
849  // Assume that input is also Gauss_Lagrange
850  // but no way to check;
852  nummodes[mode_offset],
855  p0,data, m_base[0]->GetPointsKey(), coeffs);
856  }
857  break;
858  default:
859  ASSERTL0(false,
860  "basis is either not set up or not hierarchicial");
861  }
862  }
863 
864  void SegExp::v_ComputeVertexNormal(const int vertex)
865  {
866  int i;
867  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
868  GetGeom()->GetMetricInfo();
869  SpatialDomains::GeomType type = geomFactors->GetGtype();
870  const Array<TwoD, const NekDouble> &gmat =
871  geomFactors->GetDerivFactors(GetPointsKeys());
872  int nqe = m_base[0]->GetNumPoints();
873  int vCoordDim = GetCoordim();
874 
875  m_vertexNormals[vertex] =
878  m_vertexNormals[vertex];
879  for (i = 0; i < vCoordDim; ++i)
880  {
881  normal[i] = Array<OneD, NekDouble>(nqe);
882  }
883 
884  // Regular geometry case
885  if ((type == SpatialDomains::eRegular) ||
887  {
888  NekDouble vert;
889  // Set up normals
890  switch (vertex)
891  {
892  case 0:
893  for(i = 0; i < vCoordDim; ++i)
894  {
895  Vmath::Fill(nqe, -gmat[i][0], normal[i], 1);
896  }
897  break;
898  case 1:
899  for(i = 0; i < vCoordDim; ++i)
900  {
901  Vmath::Fill(nqe, gmat[i][0], normal[i], 1);
902  }
903  break;
904  default:
905  ASSERTL0(false,
906  "point is out of range (point < 2)");
907  }
908 
909  // normalise
910  vert = 0.0;
911  for (i =0 ; i < vCoordDim; ++i)
912  {
913  vert += normal[i][0]*normal[i][0];
914  }
915  vert = 1.0/sqrt(vert);
916  for (i = 0; i < vCoordDim; ++i)
917  {
918  Vmath::Smul(nqe, vert, normal[i], 1, normal[i], 1);
919  }
920  }
921  }
922 
923  //-----------------------------
924  // Operator creation functions
925  //-----------------------------
926 
927 
929  const Array<OneD, const NekDouble> &inarray,
930  Array<OneD, NekDouble> &outarray,
931  const StdRegions::StdMatrixKey &mkey)
932  {
933  int nquad = m_base[0]->GetNumPoints();
934  const Array<TwoD, const NekDouble>& gmat =
935  m_metricinfo->GetDerivFactors(GetPointsKeys());
936 
937  Array<OneD,NekDouble> physValues(nquad);
938  Array<OneD,NekDouble> dPhysValuesdx(nquad);
939 
940  BwdTrans(inarray,physValues);
941 
942  // Laplacian matrix operation
943  switch (m_geom->GetCoordim())
944  {
945  case 1:
946  {
947  PhysDeriv(physValues,dPhysValuesdx);
948 
949  // multiply with the proper geometric factors
950  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
951  {
952  Vmath::Vmul(nquad,
953  &gmat[0][0],1,dPhysValuesdx.get(),1,
954  dPhysValuesdx.get(),1);
955  }
956  else
957  {
958  Blas::Dscal(nquad,
959  gmat[0][0], dPhysValuesdx.get(), 1);
960  }
961  }
962  break;
963  case 2:
964  {
965  Array<OneD,NekDouble> dPhysValuesdy(nquad);
966 
967  PhysDeriv(physValues,dPhysValuesdx,dPhysValuesdy);
968 
969  // multiply with the proper geometric factors
970  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
971  {
972  Vmath::Vmul (nquad,
973  &gmat[0][0],1,dPhysValuesdx.get(),1,
974  dPhysValuesdx.get(),1);
975  Vmath::Vvtvp(nquad,
976  &gmat[1][0],1,dPhysValuesdy.get(),1,
977  dPhysValuesdx.get(),1,
978  dPhysValuesdx.get(),1);
979  }
980  else
981  {
982  Blas::Dscal(nquad,
983  gmat[0][0], dPhysValuesdx.get(), 1);
984  Blas::Daxpy(nquad,
985  gmat[1][0], dPhysValuesdy.get(), 1,
986  dPhysValuesdx.get(), 1);
987  }
988  }
989  break;
990  case 3:
991  {
992  Array<OneD,NekDouble> dPhysValuesdy(nquad);
993  Array<OneD,NekDouble> dPhysValuesdz(nquad);
994 
995  PhysDeriv(physValues,dPhysValuesdx,
996  dPhysValuesdy,dPhysValuesdz);
997 
998  // multiply with the proper geometric factors
999  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1000  {
1001  Vmath::Vmul (nquad,
1002  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1003  dPhysValuesdx.get(),1);
1004  Vmath::Vvtvp(nquad,
1005  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1006  dPhysValuesdx.get(),1,
1007  dPhysValuesdx.get(),1);
1008  Vmath::Vvtvp(nquad,
1009  &gmat[2][0], 1, dPhysValuesdz.get(), 1,
1010  dPhysValuesdx.get(),1,
1011  dPhysValuesdx.get(),1);
1012  }
1013  else
1014  {
1015  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1016  Blas::Daxpy(nquad,
1017  gmat[1][0], dPhysValuesdy.get(), 1,
1018  dPhysValuesdx.get(), 1);
1019  Blas::Daxpy(nquad,
1020  gmat[2][0], dPhysValuesdz.get(), 1,
1021  dPhysValuesdx.get(), 1);
1022  }
1023  }
1024  break;
1025  default:
1026  ASSERTL0(false,"Wrong number of dimensions");
1027  break;
1028  }
1029 
1030  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
1031  }
1032 
1033 
1035  const Array<OneD, const NekDouble> &inarray,
1036  Array<OneD, NekDouble> &outarray,
1037  const StdRegions::StdMatrixKey &mkey)
1038  {
1039  int nquad = m_base[0]->GetNumPoints();
1040  const Array<TwoD, const NekDouble>& gmat =
1041  m_metricinfo->GetDerivFactors(GetPointsKeys());
1042  const NekDouble lambda =
1044 
1045  Array<OneD,NekDouble> physValues(nquad);
1046  Array<OneD,NekDouble> dPhysValuesdx(nquad);
1048 
1049  BwdTrans(inarray, physValues);
1050 
1051  // mass matrix operation
1052  v_IProductWRTBase((m_base[0]->GetBdata()),physValues,wsp,1);
1053 
1054  // Laplacian matrix operation
1055  switch (m_geom->GetCoordim())
1056  {
1057  case 1:
1058  {
1059  PhysDeriv(physValues,dPhysValuesdx);
1060 
1061  // multiply with the proper geometric factors
1062  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1063  {
1064  Vmath::Vmul(nquad,
1065  &gmat[0][0],1,dPhysValuesdx.get(),1,
1066  dPhysValuesdx.get(),1);
1067  }
1068  else
1069  {
1070  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1071  }
1072  }
1073  break;
1074  case 2:
1075  {
1076  Array<OneD,NekDouble> dPhysValuesdy(nquad);
1077 
1078  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy);
1079 
1080  // multiply with the proper geometric factors
1081  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1082  {
1083  Vmath::Vmul (nquad,
1084  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1085  dPhysValuesdx.get(), 1);
1086  Vmath::Vvtvp(nquad,
1087  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1088  dPhysValuesdx.get(), 1,
1089  dPhysValuesdx.get(), 1);
1090  }
1091  else
1092  {
1093  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1094  Blas::Daxpy(nquad,
1095  gmat[1][0], dPhysValuesdy.get(), 1,
1096  dPhysValuesdx.get(), 1);
1097  }
1098  }
1099  break;
1100  case 3:
1101  {
1102  Array<OneD,NekDouble> dPhysValuesdy(nquad);
1103  Array<OneD,NekDouble> dPhysValuesdz(nquad);
1104 
1105  PhysDeriv(physValues, dPhysValuesdx,
1106  dPhysValuesdy, dPhysValuesdz);
1107 
1108  // multiply with the proper geometric factors
1109  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1110  {
1111  Vmath::Vmul (nquad,
1112  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1113  dPhysValuesdx.get(), 1);
1114  Vmath::Vvtvp(nquad,
1115  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1116  dPhysValuesdx.get(), 1,
1117  dPhysValuesdx.get(), 1);
1118  Vmath::Vvtvp(nquad,
1119  &gmat[2][0], 1, dPhysValuesdz.get(), 1,
1120  dPhysValuesdx.get(), 1,
1121  dPhysValuesdx.get(), 1);
1122  }
1123  else
1124  {
1125  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1126  Blas::Daxpy(nquad,
1127  gmat[1][0], dPhysValuesdy.get(), 1,
1128  dPhysValuesdx.get(), 1);
1129  Blas::Daxpy(nquad,
1130  gmat[2][0], dPhysValuesdz.get(),
1131  1, dPhysValuesdx.get(), 1);
1132  }
1133  }
1134  break;
1135  default:
1136  ASSERTL0(false,"Wrong number of dimensions");
1137  break;
1138  }
1139 
1140  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
1141  Blas::Daxpy(m_ncoeffs, lambda, wsp.get(), 1, outarray.get(), 1);
1142  }
1143 
1144  //-----------------------------
1145  // Matrix creation functions
1146  //-----------------------------
1147 
1148 
1150  const MatrixKey &mkey)
1151  {
1152  return m_staticCondMatrixManager[mkey];
1153  }
1154 
1156  {
1157  m_staticCondMatrixManager.DeleteObject(mkey);
1158  }
1159 
1161  {
1162  return m_matrixManager[mkey];
1163  }
1164 
1165 
1167  const StdRegions::StdMatrixKey &mkey)
1168  {
1169  LibUtilities::BasisKey bkey = m_base[0]->GetBasisKey();
1172 
1173  return tmp->GetStdMatrix(mkey);
1174  }
1175 
1176 
1178  {
1179  DNekScalMatSharedPtr returnval;
1180  NekDouble fac;
1182 
1183  ASSERTL2(m_metricinfo->GetGtype() !=
1185  "Geometric information is not set up");
1186 
1187  switch (mkey.GetMatrixType())
1188  {
1189  case StdRegions::eMass:
1190  {
1191  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1192  || (mkey.GetNVarCoeff()))
1193  {
1194  fac = 1.0;
1195  goto UseLocRegionsMatrix;
1196  }
1197  else
1198  {
1199  fac = (m_metricinfo->GetJac(ptsKeys))[0];
1200  goto UseStdRegionsMatrix;
1201  }
1202  }
1203  break;
1204  case StdRegions::eInvMass:
1205  {
1206  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1207  || (mkey.GetNVarCoeff()))
1208  {
1209  NekDouble one = 1.0;
1210  StdRegions::StdMatrixKey masskey(
1211  StdRegions::eMass,DetShapeType(), *this);
1212  DNekMatSharedPtr mat = GenMatrix(masskey);
1213  mat->Invert();
1214 
1215  returnval = MemoryManager<DNekScalMat>::
1216  AllocateSharedPtr(one,mat);
1217  }
1218  else
1219  {
1220  fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1221  goto UseStdRegionsMatrix;
1222  }
1223  }
1224  break;
1228  {
1229  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1230  mkey.GetNVarCoeff())
1231  {
1232  fac = 1.0;
1233  goto UseLocRegionsMatrix;
1234  }
1235  else
1236  {
1237  int dir = 0;
1238  switch(mkey.GetMatrixType())
1239  {
1241  dir = 0;
1242  break;
1244  ASSERTL1(m_geom->GetCoordim() >= 2,
1245  "Cannot call eWeakDeriv2 in a "
1246  "coordinate system which is not at "
1247  "least two-dimensional");
1248  dir = 1;
1249  break;
1251  ASSERTL1(m_geom->GetCoordim() == 3,
1252  "Cannot call eWeakDeriv2 in a "
1253  "coordinate system which is not "
1254  "three-dimensional");
1255  dir = 2;
1256  break;
1257  default:
1258  break;
1259  }
1260 
1262  mkey.GetShapeType(), *this);
1263 
1264  DNekMatSharedPtr WeakDerivStd = GetStdMatrix(deriv0key);
1265  fac = m_metricinfo->GetDerivFactors(ptsKeys)[dir][0]*
1266  m_metricinfo->GetJac(ptsKeys)[0];
1267 
1268  returnval = MemoryManager<DNekScalMat>::
1269  AllocateSharedPtr(fac,WeakDerivStd);
1270  }
1271  }
1272  break;
1274  {
1275  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1276  {
1277  fac = 1.0;
1278  goto UseLocRegionsMatrix;
1279  }
1280  else
1281  {
1282  int coordim = m_geom->GetCoordim();
1283  fac = 0.0;
1284  for (int i = 0; i < coordim; ++i)
1285  {
1286  fac += m_metricinfo->GetDerivFactors(ptsKeys)[i][0]*
1287  m_metricinfo->GetDerivFactors(ptsKeys)[i][0];
1288  }
1289  fac *= m_metricinfo->GetJac(ptsKeys)[0];
1290  goto UseStdRegionsMatrix;
1291  }
1292  }
1293  break;
1295  {
1296  NekDouble factor =
1298  MatrixKey masskey(StdRegions::eMass,
1299  mkey.GetShapeType(), *this);
1300  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1302  *this, mkey.GetConstFactors(),
1303  mkey.GetVarCoeffs());
1304  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1305 
1306  int rows = LapMat.GetRows();
1307  int cols = LapMat.GetColumns();
1308 
1309  DNekMatSharedPtr helm =
1311 
1312  NekDouble one = 1.0;
1313  (*helm) = LapMat + factor*MassMat;
1314 
1315  returnval =
1317  }
1318  break;
1323  {
1324  NekDouble one = 1.0;
1325 
1326  DNekMatSharedPtr mat = GenMatrix(mkey);
1327  returnval =
1329  }
1330  break;
1332  {
1333  NekDouble one = 1.0;
1334 
1335 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1336 // DetShapeType(),*this,
1337 // mkey.GetConstant(0),
1338 // mkey.GetConstant(1));
1340  DetShapeType(),
1341  *this, mkey.GetConstFactors(),
1342  mkey.GetVarCoeffs());
1343  DNekMatSharedPtr mat = GenMatrix(hkey);
1344 
1345  mat->Invert();
1346  returnval =
1348  }
1349  break;
1351  {
1352  DNekMatSharedPtr m_Ix;
1353  Array<OneD, NekDouble> coords(1, 0.0);
1355  int vertex = (int)factors[StdRegions::eFactorGaussVertex];
1356 
1357  coords[0] = (vertex == 0) ? -1.0 : 1.0;
1358 
1359  m_Ix = m_base[0]->GetI(coords);
1360  returnval =
1362  }
1363  break;
1364 
1365  UseLocRegionsMatrix:
1366  {
1367  DNekMatSharedPtr mat = GenMatrix(mkey);
1368  returnval =
1370  }
1371  break;
1372  UseStdRegionsMatrix:
1373  {
1374  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1375  returnval =
1377  }
1378  break;
1379  default:
1380  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
1381  break;
1382  }
1383 
1384  return returnval;
1385  }
1386 
1388  const StdRegions::StdMatrixKey &mkey)
1389  {
1390  DNekMatSharedPtr returnval;
1391 
1392  switch (mkey.GetMatrixType())
1393  {
1400  returnval = Expansion1D::v_GenMatrix(mkey);
1401  break;
1402  default:
1403  returnval = StdSegExp::v_GenMatrix(mkey);
1404  break;
1405  }
1406 
1407  return returnval;
1408  }
1409 
1410 
1412  const MatrixKey &mkey)
1413  {
1414  DNekScalBlkMatSharedPtr returnval;
1415 
1416  ASSERTL2(m_metricinfo->GetGtype() !=
1418  "Geometric information is not set up");
1419 
1420  // set up block matrix system
1421  int nbdry = NumBndryCoeffs();
1422  int nint = m_ncoeffs - nbdry;
1423  Array<OneD, unsigned int> exp_size(2);
1424  exp_size[0] = nbdry;
1425  exp_size[1] = nint;
1426 
1427  /// \todo Really need a constructor which takes Arrays
1429  AllocateSharedPtr(exp_size,exp_size);
1430  NekDouble factor = 1.0;
1431 
1432  switch (mkey.GetMatrixType())
1433  {
1435  case StdRegions::eHelmholtz: // special case since Helmholtz
1436  // not defined in StdRegions
1437 
1438  // use Deformed case for both regular and deformed geometries
1439  factor = 1.0;
1440  goto UseLocRegionsMatrix;
1441  break;
1442  default:
1443  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1444  {
1445  factor = 1.0;
1446  goto UseLocRegionsMatrix;
1447  }
1448  else
1449  {
1450  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1451  factor = mat->Scale();
1452  goto UseStdRegionsMatrix;
1453  }
1454  break;
1455  UseStdRegionsMatrix:
1456  {
1457  NekDouble invfactor = 1.0/factor;
1458  NekDouble one = 1.0;
1460  DNekScalMatSharedPtr Atmp;
1461  DNekMatSharedPtr Asubmat;
1462 
1463  returnval->SetBlock(0,0,Atmp =
1465  factor,Asubmat = mat->GetBlock(0,0)));
1466  returnval->SetBlock(0,1,Atmp =
1468  one,Asubmat = mat->GetBlock(0,1)));
1469  returnval->SetBlock(1,0,Atmp =
1471  factor,Asubmat = mat->GetBlock(1,0)));
1472  returnval->SetBlock(1,1,Atmp =
1474  invfactor,Asubmat = mat->GetBlock(1,1)));
1475  }
1476  break;
1477  UseLocRegionsMatrix:
1478  {
1479  int i,j;
1480  NekDouble invfactor = 1.0/factor;
1481  NekDouble one = 1.0;
1482  DNekScalMat &mat = *GetLocMatrix(mkey);
1483  DNekMatSharedPtr A =
1485  DNekMatSharedPtr B =
1487  DNekMatSharedPtr C =
1489  DNekMatSharedPtr D =
1491 
1492  Array<OneD,unsigned int> bmap(nbdry);
1493  Array<OneD,unsigned int> imap(nint);
1494  GetBoundaryMap(bmap);
1495  GetInteriorMap(imap);
1496 
1497  for (i = 0; i < nbdry; ++i)
1498  {
1499  for (j = 0; j < nbdry; ++j)
1500  {
1501  (*A)(i,j) = mat(bmap[i],bmap[j]);
1502  }
1503 
1504  for (j = 0; j < nint; ++j)
1505  {
1506  (*B)(i,j) = mat(bmap[i],imap[j]);
1507  }
1508  }
1509 
1510  for (i = 0; i < nint; ++i)
1511  {
1512  for (j = 0; j < nbdry; ++j)
1513  {
1514  (*C)(i,j) = mat(imap[i],bmap[j]);
1515  }
1516 
1517  for (j = 0; j < nint; ++j)
1518  {
1519  (*D)(i,j) = mat(imap[i],imap[j]);
1520  }
1521  }
1522 
1523  // Calculate static condensed system
1524  if (nint)
1525  {
1526  D->Invert();
1527  (*B) = (*B)*(*D);
1528  (*A) = (*A) - (*B)*(*C);
1529  }
1530 
1531  DNekScalMatSharedPtr Atmp;
1532 
1533  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
1534  AllocateSharedPtr(factor,A));
1535  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
1536  AllocateSharedPtr(one,B));
1537  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
1538  AllocateSharedPtr(factor,C));
1539  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
1540  AllocateSharedPtr(invfactor,D));
1541  }
1542  }
1543 
1544 
1545  return returnval;
1546  }
1547 
1548 
1549  //-----------------------------
1550  // Private methods
1551  //-----------------------------
1552 
1553 
1554  /// Reverse the coefficients in a boundary interior expansion
1555  /// this routine is of use when we need the segment
1556  /// coefficients corresponding to a expansion in the reverse
1557  /// coordinate direction
1559  const Array<OneD,NekDouble> &inarray,
1560  Array<OneD,NekDouble> &outarray)
1561  {
1562 
1563  int m;
1564  NekDouble sgn = 1;
1565 
1566  ASSERTL1(&inarray[0] != &outarray[0],
1567  "inarray and outarray can not be the same");
1568  switch(GetBasisType(0))
1569  {
1571  //Swap vertices
1572  outarray[0] = inarray[1];
1573  outarray[1] = inarray[0];
1574  // negate odd modes
1575  for(m = 2; m < m_ncoeffs; ++m)
1576  {
1577  outarray[m] = sgn*inarray[m];
1578  sgn = -sgn;
1579  }
1580  break;
1583  for(m = 0; m < m_ncoeffs; ++m)
1584  {
1585  outarray[m_ncoeffs-1-m] = inarray[m];
1586  }
1587  break;
1588  default:
1589  ASSERTL0(false,"This basis is not allowed in this method");
1590  break;
1591  }
1592  }
1593 
1594  /* \brief Mass inversion product from \a inarray to \a outarray
1595  *
1596  * Multiply by the inverse of the mass matrix
1597  * \f$ {\bf \hat{u}} = {\bf M}^{-1} {\bf I} \f$
1598  *
1599  **/
1601  const Array<OneD, const NekDouble>& inarray,
1602  Array<OneD,NekDouble> &outarray)
1603  {
1604  // get Mass matrix inverse
1606  DetShapeType(),*this);
1607  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
1608 
1609  NekVector<NekDouble> in(m_ncoeffs,inarray,eCopy);
1610  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
1611 
1612  out = (*matsys)*in;
1613  }
1614 
1615 
1616  } // end of namespace
1617 }//end of namespace
1618 
const LibUtilities::PointsKeyVector GetPointsKeys() const
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1149
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
const NormalVector & GetEdgeNormal(const int edge) const
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:118
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:142
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: SegExp.cpp:800
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: SegExp.cpp:670
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:168
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1155
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:931
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: SegExp.h:237
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion1D.cpp:43
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: SegExp.cpp:649
std::map< int, NormalVector > m_vertexNormals
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: SegExp.cpp:1034
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
virtual void v_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...
Definition: SegExp.cpp:375
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:229
virtual void v_PhysDeriv_n(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dn)
Evaluate the derivative normal to a line: . The derivative is calculated performing the product ...
Definition: SegExp.cpp:267
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
Principle Modified Functions .
Definition: BasisType.h:49
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs)
Unpack data from input file assuming it comes from.
Definition: SegExp.cpp:820
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:126
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
STL namespace.
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
virtual int v_GetNcoeffs(void) const
Definition: SegExp.cpp:795
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: SegExp.cpp:768
void Vdiv(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:227
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:125
virtual void v_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: SegExp.cpp:501
virtual void v_PhysDeriv_s(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_ds)
Evaluate the derivative along a line: . The derivative is calculated performing the product ...
Definition: SegExp.cpp:227
boost::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:47
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
virtual int v_GetNumPoints(const int dir) const
Definition: SegExp.cpp:790
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
virtual void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:619
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion1D.h:116
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:689
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:800
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 ReverseCoeffsAndSign(const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Reverse the coefficients in a boundary interior expansion this routine is of use when we need the seg...
Definition: SegExp.cpp:1558
virtual SpatialDomains::GeomType v_MetricInfoType()
Definition: SegExp.cpp:785
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: SegExp.cpp:1387
virtual StdRegions::Orientation v_GetPorient(int point)
Definition: SegExp.cpp:762
Principle Modified Functions .
Definition: BasisType.h:50
virtual void v_ComputeVertexNormal(const int vertex)
Definition: SegExp.cpp:864
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
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
virtual const Array< OneD, const NekDouble > & v_GetPhysNormals(void)
Definition: SegExp.cpp:779
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1411
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:399
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: SegExp.cpp:928
T Ddot(int n, const Array< OneD, const T > &w, const int incw, const Array< OneD, const T > &x, const int incx, const Array< OneD, const int > &y, const int incy)
Definition: VmathArray.hpp:434
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1177
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1160
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
virtual void v_GetVertexPhysVals(const int vertex, const Array< OneD, const NekDouble > &inarray, NekDouble &outarray)
Definition: SegExp.cpp:695
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:83
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: SegExp.cpp:686
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:148
virtual void v_SetCoeffsToOrientation(Array< OneD, NekDouble > &coeffs, StdRegions::Orientation dir)
Definition: SegExp.cpp:735
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
virtual int v_NumDGBndryCoeffs() const
Definition: SegExp.cpp:812
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Evaluate the derivative at the physical quadrature points given by inarray and return in outarray...
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:1600
#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
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Evaluate the derivative at the physical quadrature points given by inarray and return in outarray...
Definition: SegExp.cpp:164
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:509
virtual int v_GetCoordim()
Definition: SegExp.cpp:774
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrate the physical point list inarray over region and return the value.
Definition: SegExp.cpp:116
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: SegExp.cpp:1166
GeomType
Indicates the type of element geometry.
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayofArray
boost::shared_ptr< Basis > BasisSharedPtr
Lagrange for SEM basis .
Definition: BasisType.h:53
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:560
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:216
#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 GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:790
Describes the specification for a Basis.
Definition: Basis.h:50
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
Definition: SegExp.cpp:657
virtual int v_NumBndryCoeffs() const
Definition: SegExp.cpp:806
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
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
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: SegExp.h:235