Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  {
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  {
798  2, m_base[0]->GetPointsKey());
799 
801  ::AllocateSharedPtr( bkey0);
802  }
803 
805  {
806  return m_geom->GetCoordim();
807  }
808 
810  {
811  NEKERROR(ErrorUtil::efatal, "Got to SegExp");
812  return NullNekDouble1DArray;
813  }
814 
816  {
817  return m_metricinfo->GetGtype();
818  }
819 
820  int SegExp::v_GetNumPoints(const int dir) const
821  {
822  return GetNumPoints(dir);
823  }
824 
825  int SegExp::v_GetNcoeffs(void) const
826  {
827  return m_ncoeffs;
828  }
829 
831  {
832  return GetBasis(dir);
833  }
834 
835 
837  {
838  return 2;
839  }
840 
841 
843  {
844  return 2;
845  }
846 
847 
848  /// Unpack data from input file assuming it comes from
849  // the same expansion type
851  const NekDouble *data,
852  const std::vector<unsigned int > &nummodes,
853  const int mode_offset,
854  NekDouble *coeffs,
855  std::vector<LibUtilities::BasisType> &fromType)
856  {
857  switch(m_base[0]->GetBasisType())
858  {
860  {
861  int fillorder = min((int) nummodes[mode_offset],m_ncoeffs);
862 
863  Vmath::Zero(m_ncoeffs,coeffs,1);
864  Vmath::Vcopy(fillorder,&data[0],1,&coeffs[0],1);
865  }
866  break;
868  {
869  // Assume that input is also Gll_Lagrange
870  // but no way to check;
872  nummodes[mode_offset],
875  p0,data, m_base[0]->GetPointsKey(), coeffs);
876  }
877  break;
879  {
880  // Assume that input is also Gauss_Lagrange
881  // but no way to check;
883  nummodes[mode_offset],
886  p0,data, m_base[0]->GetPointsKey(), coeffs);
887  }
888  break;
889  default:
890  ASSERTL0(false,
891  "basis is either not set up or not hierarchicial");
892  }
893  }
894 
896  {
897  int i;
898  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
899  GetGeom()->GetMetricInfo();
900  SpatialDomains::GeomType type = geomFactors->GetGtype();
901  const Array<TwoD, const NekDouble> &gmat =
902  geomFactors->GetDerivFactors(GetPointsKeys());
903  int nqe = 1;
904  int vCoordDim = GetCoordim();
905 
910  for (i = 0; i < vCoordDim; ++i)
911  {
912  normal[i] = Array<OneD, NekDouble>(nqe);
913  }
914 
915  // Regular geometry case
916  if ((type == SpatialDomains::eRegular) ||
918  {
919  NekDouble vert;
920  // Set up normals
921  switch (vertex)
922  {
923  case 0:
924  for(i = 0; i < vCoordDim; ++i)
925  {
926  Vmath::Fill(nqe, -gmat[i][0], normal[i], 1);
927  }
928  break;
929  case 1:
930  for(i = 0; i < vCoordDim; ++i)
931  {
932  Vmath::Fill(nqe, gmat[i][0], normal[i], 1);
933  }
934  break;
935  default:
936  ASSERTL0(false,
937  "point is out of range (point < 2)");
938  }
939 
940  // normalise
941  vert = 0.0;
942  for (i =0 ; i < vCoordDim; ++i)
943  {
944  vert += normal[i][0]*normal[i][0];
945  }
946  vert = 1.0/sqrt(vert);
947  for (i = 0; i < vCoordDim; ++i)
948  {
949  Vmath::Smul(nqe, vert, normal[i], 1, normal[i], 1);
950  }
951  }
952  }
953 
954  //-----------------------------
955  // Operator creation functions
956  //-----------------------------
957 
958 
960  const Array<OneD, const NekDouble> &inarray,
961  Array<OneD, NekDouble> &outarray,
962  const StdRegions::StdMatrixKey &mkey)
963  {
964  int nquad = m_base[0]->GetNumPoints();
965  const Array<TwoD, const NekDouble>& gmat =
966  m_metricinfo->GetDerivFactors(GetPointsKeys());
967 
968  Array<OneD,NekDouble> physValues(nquad);
969  Array<OneD,NekDouble> dPhysValuesdx(nquad);
970 
971  BwdTrans(inarray,physValues);
972 
973  // Laplacian matrix operation
974  switch (m_geom->GetCoordim())
975  {
976  case 1:
977  {
978  PhysDeriv(physValues,dPhysValuesdx);
979 
980  // multiply with the proper geometric factors
981  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
982  {
983  Vmath::Vmul(nquad,
984  &gmat[0][0],1,dPhysValuesdx.get(),1,
985  dPhysValuesdx.get(),1);
986  }
987  else
988  {
989  Blas::Dscal(nquad,
990  gmat[0][0], dPhysValuesdx.get(), 1);
991  }
992  }
993  break;
994  case 2:
995  {
996  Array<OneD,NekDouble> dPhysValuesdy(nquad);
997 
998  PhysDeriv(physValues,dPhysValuesdx,dPhysValuesdy);
999 
1000  // multiply with the proper geometric factors
1001  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1002  {
1003  Vmath::Vmul (nquad,
1004  &gmat[0][0],1,dPhysValuesdx.get(),1,
1005  dPhysValuesdx.get(),1);
1006  Vmath::Vvtvp(nquad,
1007  &gmat[1][0],1,dPhysValuesdy.get(),1,
1008  dPhysValuesdx.get(),1,
1009  dPhysValuesdx.get(),1);
1010  }
1011  else
1012  {
1013  Blas::Dscal(nquad,
1014  gmat[0][0], dPhysValuesdx.get(), 1);
1015  Blas::Daxpy(nquad,
1016  gmat[1][0], dPhysValuesdy.get(), 1,
1017  dPhysValuesdx.get(), 1);
1018  }
1019  }
1020  break;
1021  case 3:
1022  {
1023  Array<OneD,NekDouble> dPhysValuesdy(nquad);
1024  Array<OneD,NekDouble> dPhysValuesdz(nquad);
1025 
1026  PhysDeriv(physValues,dPhysValuesdx,
1027  dPhysValuesdy,dPhysValuesdz);
1028 
1029  // multiply with the proper geometric factors
1030  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1031  {
1032  Vmath::Vmul (nquad,
1033  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1034  dPhysValuesdx.get(),1);
1035  Vmath::Vvtvp(nquad,
1036  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1037  dPhysValuesdx.get(),1,
1038  dPhysValuesdx.get(),1);
1039  Vmath::Vvtvp(nquad,
1040  &gmat[2][0], 1, dPhysValuesdz.get(), 1,
1041  dPhysValuesdx.get(),1,
1042  dPhysValuesdx.get(),1);
1043  }
1044  else
1045  {
1046  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1047  Blas::Daxpy(nquad,
1048  gmat[1][0], dPhysValuesdy.get(), 1,
1049  dPhysValuesdx.get(), 1);
1050  Blas::Daxpy(nquad,
1051  gmat[2][0], dPhysValuesdz.get(), 1,
1052  dPhysValuesdx.get(), 1);
1053  }
1054  }
1055  break;
1056  default:
1057  ASSERTL0(false,"Wrong number of dimensions");
1058  break;
1059  }
1060 
1061  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
1062  }
1063 
1064 
1066  const Array<OneD, const NekDouble> &inarray,
1067  Array<OneD, NekDouble> &outarray,
1068  const StdRegions::StdMatrixKey &mkey)
1069  {
1070  int nquad = m_base[0]->GetNumPoints();
1071  const Array<TwoD, const NekDouble>& gmat =
1072  m_metricinfo->GetDerivFactors(GetPointsKeys());
1073  const NekDouble lambda =
1075 
1076  Array<OneD,NekDouble> physValues(nquad);
1077  Array<OneD,NekDouble> dPhysValuesdx(nquad);
1079 
1080  BwdTrans(inarray, physValues);
1081 
1082  // mass matrix operation
1083  v_IProductWRTBase((m_base[0]->GetBdata()),physValues,wsp,1);
1084 
1085  // Laplacian matrix operation
1086  switch (m_geom->GetCoordim())
1087  {
1088  case 1:
1089  {
1090  PhysDeriv(physValues,dPhysValuesdx);
1091 
1092  // multiply with the proper geometric factors
1093  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1094  {
1095  Vmath::Vmul(nquad,
1096  &gmat[0][0],1,dPhysValuesdx.get(),1,
1097  dPhysValuesdx.get(),1);
1098  }
1099  else
1100  {
1101  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1102  }
1103  }
1104  break;
1105  case 2:
1106  {
1107  Array<OneD,NekDouble> dPhysValuesdy(nquad);
1108 
1109  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy);
1110 
1111  // multiply with the proper geometric factors
1112  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1113  {
1114  Vmath::Vmul (nquad,
1115  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1116  dPhysValuesdx.get(), 1);
1117  Vmath::Vvtvp(nquad,
1118  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1119  dPhysValuesdx.get(), 1,
1120  dPhysValuesdx.get(), 1);
1121  }
1122  else
1123  {
1124  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1125  Blas::Daxpy(nquad,
1126  gmat[1][0], dPhysValuesdy.get(), 1,
1127  dPhysValuesdx.get(), 1);
1128  }
1129  }
1130  break;
1131  case 3:
1132  {
1133  Array<OneD,NekDouble> dPhysValuesdy(nquad);
1134  Array<OneD,NekDouble> dPhysValuesdz(nquad);
1135 
1136  PhysDeriv(physValues, dPhysValuesdx,
1137  dPhysValuesdy, dPhysValuesdz);
1138 
1139  // multiply with the proper geometric factors
1140  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1141  {
1142  Vmath::Vmul (nquad,
1143  &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1144  dPhysValuesdx.get(), 1);
1145  Vmath::Vvtvp(nquad,
1146  &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1147  dPhysValuesdx.get(), 1,
1148  dPhysValuesdx.get(), 1);
1149  Vmath::Vvtvp(nquad,
1150  &gmat[2][0], 1, dPhysValuesdz.get(), 1,
1151  dPhysValuesdx.get(), 1,
1152  dPhysValuesdx.get(), 1);
1153  }
1154  else
1155  {
1156  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1157  Blas::Daxpy(nquad,
1158  gmat[1][0], dPhysValuesdy.get(), 1,
1159  dPhysValuesdx.get(), 1);
1160  Blas::Daxpy(nquad,
1161  gmat[2][0], dPhysValuesdz.get(),
1162  1, dPhysValuesdx.get(), 1);
1163  }
1164  }
1165  break;
1166  default:
1167  ASSERTL0(false,"Wrong number of dimensions");
1168  break;
1169  }
1170 
1171  v_IProductWRTBase(m_base[0]->GetDbdata(),dPhysValuesdx,outarray,1);
1172  Blas::Daxpy(m_ncoeffs, lambda, wsp.get(), 1, outarray.get(), 1);
1173  }
1174 
1175  //-----------------------------
1176  // Matrix creation functions
1177  //-----------------------------
1178 
1179 
1181  const MatrixKey &mkey)
1182  {
1183  return m_staticCondMatrixManager[mkey];
1184  }
1185 
1187  {
1188  m_staticCondMatrixManager.DeleteObject(mkey);
1189  }
1190 
1192  {
1193  return m_matrixManager[mkey];
1194  }
1195 
1196 
1198  const StdRegions::StdMatrixKey &mkey)
1199  {
1200  LibUtilities::BasisKey bkey = m_base[0]->GetBasisKey();
1203 
1204  return tmp->GetStdMatrix(mkey);
1205  }
1206 
1207 
1209  {
1210  DNekScalMatSharedPtr returnval;
1211  NekDouble fac;
1213 
1214  ASSERTL2(m_metricinfo->GetGtype() !=
1216  "Geometric information is not set up");
1217 
1218  switch (mkey.GetMatrixType())
1219  {
1220  case StdRegions::eMass:
1221  {
1222  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1223  || (mkey.GetNVarCoeff()))
1224  {
1225  fac = 1.0;
1226  goto UseLocRegionsMatrix;
1227  }
1228  else
1229  {
1230  fac = (m_metricinfo->GetJac(ptsKeys))[0];
1231  goto UseStdRegionsMatrix;
1232  }
1233  }
1234  break;
1235  case StdRegions::eInvMass:
1236  {
1237  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1238  || (mkey.GetNVarCoeff()))
1239  {
1240  NekDouble one = 1.0;
1241  StdRegions::StdMatrixKey masskey(
1242  StdRegions::eMass,DetShapeType(), *this);
1243  DNekMatSharedPtr mat = GenMatrix(masskey);
1244  mat->Invert();
1245 
1246  returnval = MemoryManager<DNekScalMat>::
1247  AllocateSharedPtr(one,mat);
1248  }
1249  else
1250  {
1251  fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
1252  goto UseStdRegionsMatrix;
1253  }
1254  }
1255  break;
1259  {
1260  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1261  mkey.GetNVarCoeff())
1262  {
1263  fac = 1.0;
1264  goto UseLocRegionsMatrix;
1265  }
1266  else
1267  {
1268  int dir = 0;
1269  switch(mkey.GetMatrixType())
1270  {
1272  dir = 0;
1273  break;
1275  ASSERTL1(m_geom->GetCoordim() >= 2,
1276  "Cannot call eWeakDeriv2 in a "
1277  "coordinate system which is not at "
1278  "least two-dimensional");
1279  dir = 1;
1280  break;
1282  ASSERTL1(m_geom->GetCoordim() == 3,
1283  "Cannot call eWeakDeriv2 in a "
1284  "coordinate system which is not "
1285  "three-dimensional");
1286  dir = 2;
1287  break;
1288  default:
1289  break;
1290  }
1291 
1293  mkey.GetShapeType(), *this);
1294 
1295  DNekMatSharedPtr WeakDerivStd = GetStdMatrix(deriv0key);
1296  fac = m_metricinfo->GetDerivFactors(ptsKeys)[dir][0]*
1297  m_metricinfo->GetJac(ptsKeys)[0];
1298 
1299  returnval = MemoryManager<DNekScalMat>::
1300  AllocateSharedPtr(fac,WeakDerivStd);
1301  }
1302  }
1303  break;
1305  {
1306  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1307  {
1308  fac = 1.0;
1309  goto UseLocRegionsMatrix;
1310  }
1311  else
1312  {
1313  int coordim = m_geom->GetCoordim();
1314  fac = 0.0;
1315  for (int i = 0; i < coordim; ++i)
1316  {
1317  fac += m_metricinfo->GetDerivFactors(ptsKeys)[i][0]*
1318  m_metricinfo->GetDerivFactors(ptsKeys)[i][0];
1319  }
1320  fac *= m_metricinfo->GetJac(ptsKeys)[0];
1321  goto UseStdRegionsMatrix;
1322  }
1323  }
1324  break;
1326  {
1327  NekDouble factor =
1329  MatrixKey masskey(StdRegions::eMass,
1330  mkey.GetShapeType(), *this);
1331  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1333  *this, mkey.GetConstFactors(),
1334  mkey.GetVarCoeffs());
1335  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1336 
1337  int rows = LapMat.GetRows();
1338  int cols = LapMat.GetColumns();
1339 
1340  DNekMatSharedPtr helm =
1342 
1343  NekDouble one = 1.0;
1344  (*helm) = LapMat + factor*MassMat;
1345 
1346  returnval =
1348  }
1349  break;
1354  {
1355  NekDouble one = 1.0;
1356 
1357  DNekMatSharedPtr mat = GenMatrix(mkey);
1358  returnval =
1360  }
1361  break;
1363  {
1364  NekDouble one = 1.0;
1365 
1366 // StdRegions::StdMatrixKey hkey(StdRegions::eHybridDGHelmholtz,
1367 // DetShapeType(),*this,
1368 // mkey.GetConstant(0),
1369 // mkey.GetConstant(1));
1371  DetShapeType(),
1372  *this, mkey.GetConstFactors(),
1373  mkey.GetVarCoeffs());
1374  DNekMatSharedPtr mat = GenMatrix(hkey);
1375 
1376  mat->Invert();
1377  returnval =
1379  }
1380  break;
1382  {
1383  DNekMatSharedPtr m_Ix;
1384  Array<OneD, NekDouble> coords(1, 0.0);
1386  int vertex = (int)factors[StdRegions::eFactorGaussVertex];
1387 
1388  coords[0] = (vertex == 0) ? -1.0 : 1.0;
1389 
1390  m_Ix = m_base[0]->GetI(coords);
1391  returnval =
1393  }
1394  break;
1395 
1396  UseLocRegionsMatrix:
1397  {
1398  DNekMatSharedPtr mat = GenMatrix(mkey);
1399  returnval =
1401  }
1402  break;
1403  UseStdRegionsMatrix:
1404  {
1405  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1406  returnval =
1408  }
1409  break;
1410  default:
1411  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
1412  break;
1413  }
1414 
1415  return returnval;
1416  }
1417 
1419  const StdRegions::StdMatrixKey &mkey)
1420  {
1421  DNekMatSharedPtr returnval;
1422 
1423  switch (mkey.GetMatrixType())
1424  {
1431  returnval = Expansion1D::v_GenMatrix(mkey);
1432  break;
1433  default:
1434  returnval = StdSegExp::v_GenMatrix(mkey);
1435  break;
1436  }
1437 
1438  return returnval;
1439  }
1440 
1441 
1443  const MatrixKey &mkey)
1444  {
1445  DNekScalBlkMatSharedPtr returnval;
1446 
1447  ASSERTL2(m_metricinfo->GetGtype() !=
1449  "Geometric information is not set up");
1450 
1451  // set up block matrix system
1452  int nbdry = NumBndryCoeffs();
1453  int nint = m_ncoeffs - nbdry;
1454  Array<OneD, unsigned int> exp_size(2);
1455  exp_size[0] = nbdry;
1456  exp_size[1] = nint;
1457 
1458  /// \todo Really need a constructor which takes Arrays
1460  AllocateSharedPtr(exp_size,exp_size);
1461  NekDouble factor = 1.0;
1462 
1463  switch (mkey.GetMatrixType())
1464  {
1466  case StdRegions::eHelmholtz: // special case since Helmholtz
1467  // not defined in StdRegions
1468 
1469  // use Deformed case for both regular and deformed geometries
1470  factor = 1.0;
1471  goto UseLocRegionsMatrix;
1472  break;
1473  default:
1474  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1475  {
1476  factor = 1.0;
1477  goto UseLocRegionsMatrix;
1478  }
1479  else
1480  {
1481  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1482  factor = mat->Scale();
1483  goto UseStdRegionsMatrix;
1484  }
1485  break;
1486  UseStdRegionsMatrix:
1487  {
1488  NekDouble invfactor = 1.0/factor;
1489  NekDouble one = 1.0;
1491  DNekScalMatSharedPtr Atmp;
1492  DNekMatSharedPtr Asubmat;
1493 
1494  returnval->SetBlock(0,0,Atmp =
1496  factor,Asubmat = mat->GetBlock(0,0)));
1497  returnval->SetBlock(0,1,Atmp =
1499  one,Asubmat = mat->GetBlock(0,1)));
1500  returnval->SetBlock(1,0,Atmp =
1502  factor,Asubmat = mat->GetBlock(1,0)));
1503  returnval->SetBlock(1,1,Atmp =
1505  invfactor,Asubmat = mat->GetBlock(1,1)));
1506  }
1507  break;
1508  UseLocRegionsMatrix:
1509  {
1510  int i,j;
1511  NekDouble invfactor = 1.0/factor;
1512  NekDouble one = 1.0;
1513  DNekScalMat &mat = *GetLocMatrix(mkey);
1516  DNekMatSharedPtr B =
1518  DNekMatSharedPtr C =
1520  DNekMatSharedPtr D =
1522 
1523  Array<OneD,unsigned int> bmap(nbdry);
1524  Array<OneD,unsigned int> imap(nint);
1525  GetBoundaryMap(bmap);
1526  GetInteriorMap(imap);
1527 
1528  for (i = 0; i < nbdry; ++i)
1529  {
1530  for (j = 0; j < nbdry; ++j)
1531  {
1532  (*A)(i,j) = mat(bmap[i],bmap[j]);
1533  }
1534 
1535  for (j = 0; j < nint; ++j)
1536  {
1537  (*B)(i,j) = mat(bmap[i],imap[j]);
1538  }
1539  }
1540 
1541  for (i = 0; i < nint; ++i)
1542  {
1543  for (j = 0; j < nbdry; ++j)
1544  {
1545  (*C)(i,j) = mat(imap[i],bmap[j]);
1546  }
1547 
1548  for (j = 0; j < nint; ++j)
1549  {
1550  (*D)(i,j) = mat(imap[i],imap[j]);
1551  }
1552  }
1553 
1554  // Calculate static condensed system
1555  if (nint)
1556  {
1557  D->Invert();
1558  (*B) = (*B)*(*D);
1559  (*A) = (*A) - (*B)*(*C);
1560  }
1561 
1562  DNekScalMatSharedPtr Atmp;
1563 
1564  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::
1565  AllocateSharedPtr(factor,A));
1566  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::
1567  AllocateSharedPtr(one,B));
1568  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::
1569  AllocateSharedPtr(factor,C));
1570  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::
1571  AllocateSharedPtr(invfactor,D));
1572  }
1573  }
1574 
1575 
1576  return returnval;
1577  }
1578 
1579 
1580  //-----------------------------
1581  // Private methods
1582  //-----------------------------
1583 
1584 
1585  /// Reverse the coefficients in a boundary interior expansion
1586  /// this routine is of use when we need the segment
1587  /// coefficients corresponding to a expansion in the reverse
1588  /// coordinate direction
1590  const Array<OneD,NekDouble> &inarray,
1591  Array<OneD,NekDouble> &outarray)
1592  {
1593 
1594  int m;
1595  NekDouble sgn = 1;
1596 
1597  ASSERTL1(&inarray[0] != &outarray[0],
1598  "inarray and outarray can not be the same");
1599  switch(GetBasisType(0))
1600  {
1602  //Swap vertices
1603  outarray[0] = inarray[1];
1604  outarray[1] = inarray[0];
1605  // negate odd modes
1606  for(m = 2; m < m_ncoeffs; ++m)
1607  {
1608  outarray[m] = sgn*inarray[m];
1609  sgn = -sgn;
1610  }
1611  break;
1614  for(m = 0; m < m_ncoeffs; ++m)
1615  {
1616  outarray[m_ncoeffs-1-m] = inarray[m];
1617  }
1618  break;
1619  default:
1620  ASSERTL0(false,"This basis is not allowed in this method");
1621  break;
1622  }
1623  }
1624 
1625  /* \brief Mass inversion product from \a inarray to \a outarray
1626  *
1627  * Multiply by the inverse of the mass matrix
1628  * \f$ {\bf \hat{u}} = {\bf M}^{-1} {\bf I} \f$
1629  *
1630  **/
1632  const Array<OneD, const NekDouble>& inarray,
1633  Array<OneD,NekDouble> &outarray)
1634  {
1635  // get Mass matrix inverse
1637  DetShapeType(),*this);
1638  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
1639 
1640  NekVector<NekDouble> in(m_ncoeffs,inarray,eCopy);
1641  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
1642 
1643  out = (*matsys)*in;
1644  }
1645 
1646 
1647  } // end of namespace
1648 }//end of namespace
1649 
const LibUtilities::PointsKeyVector GetPointsKeys() const
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1180
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:122
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: SegExp.cpp:795
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:191
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: SegExp.cpp:830
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:242
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:82
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1186
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:976
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:252
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:1065
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:442
Principle Modified Functions .
Definition: BasisType.h:49
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:733
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:131
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.
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
Unpack data from input file assuming it comes from.
Definition: SegExp.cpp:850
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:87
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
virtual int v_GetNcoeffs(void) const
Definition: SegExp.cpp:825
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:241
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:130
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:820
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:706
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:49
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion1D.h:123
DNekBlkMatSharedPtr GetStdStaticCondMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:711
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:824
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
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:1589
virtual SpatialDomains::GeomType v_MetricInfoType()
Definition: SegExp.cpp:815
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: SegExp.cpp:1418
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:895
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:224
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:819
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:809
DNekScalBlkMatSharedPtr CreateStaticCondMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1442
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:959
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:436
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1208
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:165
bg::model::point< double, 3, bg::cs::cartesian > point
Definition: BLMesh.cpp:54
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1191
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:343
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:161
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:842
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:1631
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
Geometry is straight-sided with constant geometric factors.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:54
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:531
virtual int v_GetCoordim()
Definition: SegExp.cpp:804
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:1197
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:373
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:228
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Geometry is curved or has non-constant factors.
void GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:814
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:836
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: SegExp.h:250