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