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