Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NodalTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File NodalTriExp.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: NodalTriExp routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LocalRegions/NodalTriExp.h>
37 #include <LibUtilities/Foundations/Interp.h>
38 
39 
40 namespace Nektar
41 {
42  namespace LocalRegions
43  {
44  NodalTriExp::NodalTriExp(const LibUtilities::BasisKey &Ba,
45  const LibUtilities::BasisKey &Bb,
46  const LibUtilities::PointsType Ntype,
47  const SpatialDomains::TriGeomSharedPtr &geom):
48  StdExpansion (LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),2,Ba,Bb),
49  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(Ba.GetNumModes(),(Bb.GetNumModes())),Ba,Bb),
50  StdNodalTriExp(Ba,Bb,Ntype),
51  Expansion (geom),
52  Expansion2D (geom),
53  m_matrixManager(
54  boost::bind(&NodalTriExp::CreateMatrix, this, _1),
55  std::string("NodalTriExpMatrix")),
56  m_staticCondMatrixManager(
57  boost::bind(&NodalTriExp::CreateStaticCondMatrix, this, _1),
58  std::string("NodalTriExpStaticCondMatrix"))
59  {
60  }
61 
62  NodalTriExp::NodalTriExp(const NodalTriExp &T):
63  StdExpansion(T),
64  StdExpansion2D(T),
65  StdRegions::StdNodalTriExp(T),
66  Expansion (T),
67  Expansion2D (T),
68  m_matrixManager(T.m_matrixManager),
69  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
70  {
71  }
72 
73  NodalTriExp::~NodalTriExp()
74  {
75  }
76 
77  //----------------------------
78  // Integration Methods
79  //----------------------------
80 
81  /** \brief Integrate the physical point list \a inarray over region
82  and return the value
83 
84  Inputs:\n
85 
86  - \a inarray: definition of function to be returned at quadrature point
87  of expansion.
88 
89  Outputs:\n
90 
91  - returns \f$\int^1_{-1}\int^1_{-1} u(\xi_1, \xi_2) J[i,j] d
92  \xi_1 d \xi_2 \f$ where \f$inarray[i,j] = u(\xi_{1i},\xi_{2j})
93  \f$ and \f$ J[i,j] \f$ is the Jacobian evaluated at the
94  quadrature point.
95  */
96 
97 
98  NekDouble NodalTriExp::Integral(const Array<OneD, const NekDouble> &inarray)
99  {
100  int nquad0 = m_base[0]->GetNumPoints();
101  int nquad1 = m_base[1]->GetNumPoints();
102  Array<OneD, const NekDouble> jac = m_metricinfo->GetJac(GetPointsKeys());
103  NekDouble ival;
104  Array<OneD,NekDouble> tmp(nquad0*nquad1);
105 
106  // multiply inarray with Jacobian
107  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
108  {
109  Vmath::Vmul(nquad0*nquad1, jac, 1, inarray, 1,tmp, 1);
110  }
111  else
112  {
113  Vmath::Smul(nquad0*nquad1, jac[0], inarray, 1, tmp, 1);
114  }
115 
116  // call StdQuadExp version;
117  ival = StdNodalTriExp::v_Integral(tmp);
118  return ival;
119  }
120 
121 
122  void NodalTriExp::IProductWRTBase_SumFac(const Array<OneD, const NekDouble>& inarray,
123  Array<OneD, NekDouble> &outarray)
124  {
125  int nquad0 = m_base[0]->GetNumPoints();
126  int nquad1 = m_base[1]->GetNumPoints();
127  int order1 = m_base[1]->GetNumModes();
128 
129  Array<OneD,NekDouble> tmp(nquad0*nquad1+nquad0*order1);
130  Array<OneD,NekDouble> wsp(tmp+nquad0*nquad1);
131 
132  MultiplyByQuadratureMetric(inarray,tmp);
133  StdTriExp::IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),m_base[1]->GetBdata(),tmp,outarray,wsp);
134  NodalToModalTranspose(outarray,outarray);
135  }
136 
137  void NodalTriExp::IProductWRTBase_MatOp(const Array<OneD, const NekDouble>& inarray,
138  Array<OneD, NekDouble> &outarray)
139  {
140  int nq = GetTotPoints();
141  MatrixKey iprodmatkey(StdRegions::eIProductWRTBase,DetShapeType(),*this);
142  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
143 
144  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
145  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
146  }
147 
148  void NodalTriExp::IProductWRTDerivBase_SumFac(const int dir,
149  const Array<OneD, const NekDouble>& inarray,
150  Array<OneD, NekDouble> & outarray)
151  {
152  ASSERTL1((dir==0)||(dir==1)||(dir==2),"Invalid direction.");
153  ASSERTL1((dir==2)?(m_geom->GetCoordim()==3):true,"Invalid direction.");
154 
155  int i;
156  int nquad0 = m_base[0]->GetNumPoints();
157  int nquad1 = m_base[1]->GetNumPoints();
158  int nqtot = nquad0*nquad1;
159  int wspsize = max(nqtot,m_ncoeffs);
160 
161  const Array<TwoD, const NekDouble>& df =
162  m_metricinfo->GetDerivFactors(GetPointsKeys());
163 
164  Array<OneD, NekDouble> tmp0 (6*wspsize);
165  Array<OneD, NekDouble> tmp1 (tmp0 + wspsize);
166  Array<OneD, NekDouble> tmp2 (tmp0 + 2*wspsize);
167  Array<OneD, NekDouble> tmp3 (tmp0 + 3*wspsize);
168  Array<OneD, NekDouble> gfac0(tmp0 + 4*wspsize);
169  Array<OneD, NekDouble> gfac1(tmp0 + 5*wspsize);
170 
171  const Array<OneD, const NekDouble>& z0 = m_base[0]->GetZ();
172  const Array<OneD, const NekDouble>& z1 = m_base[1]->GetZ();
173 
174  // set up geometric factor: 2/(1-z1)
175  for(i = 0; i < nquad1; ++i)
176  {
177  gfac0[i] = 2.0/(1-z1[i]);
178  }
179  for(i = 0; i < nquad0; ++i)
180  {
181  gfac1[i] = 0.5*(1+z0[i]);
182  }
183 
184  for(i = 0; i < nquad1; ++i)
185  {
186  Vmath::Smul(nquad0,gfac0[i],&inarray[0]+i*nquad0,1,&tmp0[0]+i*nquad0,1);
187  }
188 
189  for(i = 0; i < nquad1; ++i)
190  {
191  Vmath::Vmul(nquad0,&gfac1[0],1,&tmp0[0]+i*nquad0,1,&tmp1[0]+i*nquad0,1);
192  }
193 
194  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
195  {
196  Vmath::Vmul(nqtot,&df[2*dir][0], 1,&tmp0[0], 1,&tmp0[0],1);
197  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&tmp1[0], 1,&tmp1[0],1);
198  Vmath::Vmul(nqtot,&df[2*dir+1][0],1,&inarray[0],1,&tmp2[0],1);
199  }
200  else
201  {
202  Vmath::Smul(nqtot, df[2*dir][0], tmp0, 1, tmp0, 1);
203  Vmath::Smul(nqtot, df[2*dir+1][0], tmp1, 1, tmp1, 1);
204  Vmath::Smul(nqtot, df[2*dir+1][0], inarray, 1, tmp2, 1);
205  }
206  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
207 
208  MultiplyByQuadratureMetric(tmp1,tmp1);
209  MultiplyByQuadratureMetric(tmp2,tmp2);
210 
211  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),m_base[1]->GetBdata() ,tmp1,tmp3 ,tmp0);
212  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata() ,m_base[1]->GetDbdata(),tmp2,outarray,tmp0);
213  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
214 
215  NodalToModalTranspose(outarray,outarray);
216  }
217 
218  void NodalTriExp::IProductWRTDerivBase_MatOp(const int dir,
219  const Array<OneD, const NekDouble>& inarray,
220  Array<OneD, NekDouble> &outarray)
221  {
222  int nq = GetTotPoints();
223  StdRegions::MatrixType mtype = StdRegions::eIProductWRTDerivBase0;
224 
225  switch(dir)
226  {
227  case 0:
228  {
229  mtype = StdRegions::eIProductWRTDerivBase0;
230  }
231  break;
232  case 1:
233  {
234  mtype = StdRegions::eIProductWRTDerivBase1;
235  }
236  break;
237  case 2:
238  {
239  mtype = StdRegions::eIProductWRTDerivBase2;
240  }
241  break;
242  default:
243  {
244  ASSERTL1(false,"input dir is out of range");
245  }
246  break;
247  }
248 
249  MatrixKey iprodmatkey(mtype,DetShapeType(),*this);
250  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
251 
252  Blas::Dgemv('N',m_ncoeffs,nq,iprodmat->Scale(),(iprodmat->GetOwnedMatrix())->GetPtr().get(),
253  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
254 
255  }
256 
257  ///////////////////////////////
258  /// Differentiation Methods
259  ///////////////////////////////
260 
261  /**
262  \brief Calculate the deritive of the physical points
263  **/
264  void NodalTriExp::PhysDeriv(const Array<OneD, const NekDouble> & inarray,
265  Array<OneD,NekDouble> &out_d0,
266  Array<OneD,NekDouble> &out_d1,
267  Array<OneD,NekDouble> &out_d2)
268  {
269  int nquad0 = m_base[0]->GetNumPoints();
270  int nquad1 = m_base[1]->GetNumPoints();
271  int nqtot = nquad0*nquad1;
272  const Array<TwoD, const NekDouble>& df
273  = m_metricinfo->GetDerivFactors(GetPointsKeys());
274 
275  Array<OneD,NekDouble> diff0(2*nqtot);
276  Array<OneD,NekDouble> diff1(diff0+nqtot);
277 
278  StdNodalTriExp::v_PhysDeriv(inarray, diff0, diff1);
279 
280  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
281  {
282  if(out_d0.num_elements())
283  {
284  Vmath::Vmul (nqtot,df[0],1,diff0,1, out_d0, 1);
285  Vmath::Vvtvp (nqtot,df[1],1,diff1,1, out_d0, 1, out_d0,1);
286  }
287 
288  if(out_d1.num_elements())
289  {
290  Vmath::Vmul (nqtot,df[2],1,diff0,1, out_d1, 1);
291  Vmath::Vvtvp (nqtot,df[3],1,diff1,1, out_d1, 1, out_d1,1);
292  }
293 
294  if(out_d2.num_elements())
295  {
296  Vmath::Vmul (nqtot,df[4],1,diff0,1, out_d2, 1);
297  Vmath::Vvtvp (nqtot,df[5],1,diff1,1, out_d2, 1, out_d2,1);
298  }
299  }
300  else // regular geometry
301  {
302  if(out_d0.num_elements())
303  {
304  Vmath::Smul (nqtot, df[0][0], diff0, 1, out_d0, 1);
305  Blas::Daxpy (nqtot, df[1][0], diff1, 1, out_d0, 1);
306  }
307 
308  if(out_d1.num_elements())
309  {
310  Vmath::Smul (nqtot, df[2][0], diff0, 1, out_d1, 1);
311  Blas::Daxpy (nqtot, df[3][0], diff1, 1, out_d1, 1);
312  }
313 
314  if(out_d2.num_elements())
315  {
316  Vmath::Smul (nqtot, df[4][0], diff0, 1, out_d2, 1);
317  Blas::Daxpy (nqtot, df[5][0], diff1, 1, out_d2, 1);
318  }
319  }
320  }
321 
322  /** \brief Forward transform from physical quadrature space
323  stored in \a inarray and evaluate the expansion coefficients and
324  store in \a (this)->m_coeffs
325 
326  Inputs:\n
327 
328  - \a inarray: array of physical quadrature points to be transformed
329 
330  Outputs:\n
331 
332  - (this)->_coeffs: updated array of expansion coefficients.
333 
334  */
335  void NodalTriExp::FwdTrans(const Array<OneD, const NekDouble> & inarray,
336  Array<OneD,NekDouble> &outarray)
337  {
338  IProductWRTBase(inarray,outarray);
339 
340  // get Mass matrix inverse
341  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(),*this,StdRegions::NullConstFactorMap,StdRegions::NullVarCoeffMap,
342  m_nodalPointsKey->GetPointsType());
343  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
344 
345  // copy inarray in case inarray == outarray
346  NekVector<NekDouble> in(m_ncoeffs,outarray,eCopy);
347  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
348 
349  out = (*matsys)*in;
350  }
351 
352  void NodalTriExp::GeneralMatrixOp_MatOp(const Array<OneD, const NekDouble> &inarray,
353  Array<OneD,NekDouble> &outarray,
354  const StdRegions::StdMatrixKey &mkey)
355  {
356  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
357 
358  if(inarray.get() == outarray.get())
359  {
360  Array<OneD,NekDouble> tmp(m_ncoeffs);
361  Vmath::Vcopy(m_ncoeffs,inarray.get(),1,tmp.get(),1);
362 
363  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
364  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
365  }
366  else
367  {
368  Blas::Dgemv('N',m_ncoeffs,m_ncoeffs,mat->Scale(),(mat->GetOwnedMatrix())->GetPtr().get(),
369  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
370  }
371  }
372 
373  void NodalTriExp::GetCoords(Array<OneD,NekDouble> &coords_0,
374  Array<OneD,NekDouble> &coords_1,
375  Array<OneD,NekDouble> &coords_2)
376  {
377  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
378  }
379 
380  // get the coordinates "coords" at the local coordinates "Lcoords"
381  void NodalTriExp::GetCoord(const Array<OneD, const NekDouble> &Lcoords,
382  Array<OneD,NekDouble> &coords)
383  {
384  int i;
385 
386  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 &&
387  Lcoords[1] >= -1.0 && Lcoords[1] <=1.0,
388  "Local coordinates are not in region [-1,1]");
389 
390  m_geom->FillGeom();
391 
392  for(i = 0; i < m_geom->GetCoordim(); ++i)
393  {
394  coords[i] = m_geom->GetCoord(i,Lcoords);
395  }
396  }
397 
398  DNekMatSharedPtr NodalTriExp::CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
399  {
400  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
401  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
402  LibUtilities::PointsType ntype = m_nodalPointsKey->GetPointsType();
403  StdRegions::StdNodalTriExpSharedPtr tmp = MemoryManager<StdNodalTriExp>::
404  AllocateSharedPtr(bkey0,bkey1,ntype);
405 
406  return tmp->GetStdMatrix(mkey);
407  }
408 
409  NekDouble NodalTriExp::PhysEvaluate(
410  const Array<OneD, const NekDouble> &coord,
411  const Array<OneD, const NekDouble> &physvals)
412 
413  {
414  Array<OneD,NekDouble> Lcoord = Array<OneD,NekDouble>(2);
415 
416  ASSERTL0(m_geom,"m_geom not defined");
417  m_geom->GetLocCoords(coord,Lcoord);
418 
419  return StdNodalTriExp::v_PhysEvaluate(Lcoord, physvals);
420  }
421 
422  DNekScalMatSharedPtr NodalTriExp::CreateMatrix(const MatrixKey &mkey)
423  {
424  DNekScalMatSharedPtr returnval;
425  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
426 
427  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
428 
429  StdRegions::MatrixType mtype = mkey.GetMatrixType();
430 
431  switch(mtype)
432  {
433  case StdRegions::eMass:
434  {
435  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
436  {
437  NekDouble one = 1.0;
438  DNekMatSharedPtr mat = GenMatrix(mkey);
439  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
440  }
441  else
442  {
443  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
444  DNekMatSharedPtr mat = GetStdMatrix(mkey);
445  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,mat);
446  }
447  }
448  break;
449  case StdRegions::eInvMass:
450  {
451  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
452  {
453  NekDouble one = 1.0;
454  StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetShapeType(),
455  *this);
456  DNekMatSharedPtr mat = GenMatrix(masskey);
457  mat->Invert();
458 
459  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
460  }
461  else
462  {
463  NekDouble fac = 1.0/(m_metricinfo->GetJac(ptsKeys))[0];
464  DNekMatSharedPtr mat = GetStdMatrix(mkey);
465  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac,mat);
466  }
467  }
468  break;
469  case StdRegions::eLaplacian:
470  {
471  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
472  {
473  NekDouble one = 1.0;
474  DNekMatSharedPtr mat = GenMatrix(mkey);
475 
476  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,mat);
477  }
478  else
479  {
480  ASSERTL1(m_geom->GetCoordim() == 2,"Standard Region Laplacian is only set up for Quads in two-dimensional");
481  MatrixKey lap00key(StdRegions::eLaplacian00,
482  mkey.GetShapeType(), *this);
483  MatrixKey lap01key(StdRegions::eLaplacian01,
484  mkey.GetShapeType(), *this);
485  MatrixKey lap11key(StdRegions::eLaplacian11,
486  mkey.GetShapeType(), *this);
487 
488  DNekMat &lap00 = *GetStdMatrix(lap00key);
489  DNekMat &lap01 = *GetStdMatrix(lap01key);
490  DNekMat &lap11 = *GetStdMatrix(lap11key);
491 
492  NekDouble jac = (m_metricinfo->GetJac(ptsKeys))[0];
493  Array<TwoD, const NekDouble> gmat =
494  m_metricinfo->GetGmat(ptsKeys);
495 
496  int rows = lap00.GetRows();
497  int cols = lap00.GetColumns();
498 
499  DNekMatSharedPtr lap = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
500 
501  (*lap) = gmat[0][0] * lap00 +
502  gmat[1][0] * (lap01 + Transpose(lap01)) +
503  gmat[3][0] * lap11;
504 
505  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(jac,lap);
506  }
507  }
508  break;
509  case StdRegions::eHelmholtz:
510  {
511  NekDouble factor = mkey.GetConstFactor(StdRegions::eFactorLambda);
512  MatrixKey masskey(StdRegions::eMass,
513  mkey.GetShapeType(), *this);
514  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
515  MatrixKey lapkey(StdRegions::eLaplacian,
516  mkey.GetShapeType(), *this);
517  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
518 
519  int rows = LapMat.GetRows();
520  int cols = LapMat.GetColumns();
521 
522  DNekMatSharedPtr helm = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols);
523 
524  NekDouble one = 1.0;
525  (*helm) = LapMat + factor*MassMat;
526 
527  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,helm);
528  }
529  break;
530  default:
531  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
532  break;
533  }
534 
535  return returnval;
536  }
537 
538  DNekScalBlkMatSharedPtr NodalTriExp::CreateStaticCondMatrix(const MatrixKey &mkey)
539  {
540  DNekScalBlkMatSharedPtr returnval;
541 
542  ASSERTL2(m_metricinfo->GetGtype() != SpatialDomains::eNoGeomType,"Geometric information is not set up");
543 
544  // set up block matrix system
545  unsigned int nbdry = NumBndryCoeffs();
546  unsigned int nint = (unsigned int)(m_ncoeffs - nbdry);
547  unsigned int exp_size[] = {nbdry,nint};
548  unsigned int nblks = 2;
549  returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks,nblks,exp_size,exp_size); //Really need a constructor which takes Arrays
550  NekDouble factor = 1.0;
551 
552  switch(mkey.GetMatrixType())
553  {
554  case StdRegions::eLaplacian:
555  case StdRegions::eHelmholtz: // special case since Helmholtz not defined in StdRegions
556 
557  // use Deformed case for both regular and deformed geometries
558  factor = 1.0;
559  goto UseLocRegionsMatrix;
560  break;
561  default:
562  if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
563  {
564  factor = 1.0;
565  goto UseLocRegionsMatrix;
566  }
567  else
568  {
569  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
570  factor = mat->Scale();
571  goto UseStdRegionsMatrix;
572  }
573  break;
574  UseStdRegionsMatrix:
575  {
576  NekDouble invfactor = 1.0/factor;
577  NekDouble one = 1.0;
578  DNekBlkMatSharedPtr mat = GetStdStaticCondMatrix(mkey);
579  DNekScalMatSharedPtr Atmp;
580  DNekMatSharedPtr Asubmat;
581 
582  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(0,0)));
583  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,Asubmat = mat->GetBlock(0,1)));
584  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,Asubmat = mat->GetBlock(1,0)));
585  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,Asubmat = mat->GetBlock(1,1)));
586  }
587  break;
588  UseLocRegionsMatrix:
589  {
590  int i,j;
591  NekDouble invfactor = 1.0/factor;
592  NekDouble one = 1.0;
593  DNekScalMat &mat = *GetLocMatrix(mkey);
594  DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
595  DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
596  DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
597  DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
598 
599  Array<OneD,unsigned int> bmap(nbdry);
600  Array<OneD,unsigned int> imap(nint);
601  GetBoundaryMap(bmap);
602  GetInteriorMap(imap);
603 
604  for(i = 0; i < nbdry; ++i)
605  {
606  for(j = 0; j < nbdry; ++j)
607  {
608  (*A)(i,j) = mat(bmap[i],bmap[j]);
609  }
610 
611  for(j = 0; j < nint; ++j)
612  {
613  (*B)(i,j) = mat(bmap[i],imap[j]);
614  }
615  }
616 
617  for(i = 0; i < nint; ++i)
618  {
619  for(j = 0; j < nbdry; ++j)
620  {
621  (*C)(i,j) = mat(imap[i],bmap[j]);
622  }
623 
624  for(j = 0; j < nint; ++j)
625  {
626  (*D)(i,j) = mat(imap[i],imap[j]);
627  }
628  }
629 
630  // Calculate static condensed system
631  if(nint)
632  {
633  D->Invert();
634  (*B) = (*B)*(*D);
635  (*A) = (*A) - (*B)*(*C);
636  }
637 
638  DNekScalMatSharedPtr Atmp;
639 
640  returnval->SetBlock(0,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,A));
641  returnval->SetBlock(0,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,B));
642  returnval->SetBlock(1,0,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(factor,C));
643  returnval->SetBlock(1,1,Atmp = MemoryManager<DNekScalMat>::AllocateSharedPtr(invfactor,D));
644 
645  }
646  }
647 
648  return returnval;
649  }
650 
651  DNekMatSharedPtr NodalTriExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
652  {
653  DNekMatSharedPtr returnval;
654 
655  switch(mkey.GetMatrixType())
656  {
657  case StdRegions::eHybridDGHelmholtz:
658  case StdRegions::eHybridDGLamToU:
659  case StdRegions::eHybridDGLamToQ0:
660  case StdRegions::eHybridDGLamToQ1:
661  case StdRegions::eHybridDGLamToQ2:
662  case StdRegions::eHybridDGHelmBndLam:
663  returnval = Expansion2D::v_GenMatrix(mkey);
664  break;
665  default:
666  returnval = StdNodalTriExp::v_GenMatrix(mkey);
667  break;
668  }
669  return returnval;
670  }
671 
672  void NodalTriExp::v_ComputeEdgeNormal(const int edge)
673  {
674  int i;
675  const SpatialDomains::GeomFactorsSharedPtr & geomFactors = GetGeom()->GetMetricInfo();
676  const SpatialDomains::GeomType type = geomFactors->GetGtype();
677 
678  LibUtilities::PointsKeyVector ptsKeys = GetPointsKeys();
679  const Array<TwoD, const NekDouble> & df = geomFactors->GetDerivFactors(ptsKeys);
680  const Array<OneD, const NekDouble> & jac = geomFactors->GetJac(ptsKeys);
681  int nqe = m_base[0]->GetNumPoints();
682  int dim = GetCoordim();
683 
684  m_edgeNormals[edge] = Array<OneD, Array<OneD, NekDouble> >(dim);
685  Array<OneD, Array<OneD, NekDouble> > &normal = m_edgeNormals[edge];
686  for (i = 0; i < dim; ++i)
687  {
688  normal[i] = Array<OneD, NekDouble>(nqe);
689  }
690 
691  // Regular geometry case
692  if((type == SpatialDomains::eRegular)||(type == SpatialDomains::eMovingRegular))
693  {
694  NekDouble fac;
695  // Set up normals
696  switch(edge)
697  {
698  case 0:
699  for(i = 0; i < GetCoordim(); ++i)
700  {
701  Vmath::Fill(nqe,-df[2*i+1][0],normal[i],1);
702  }
703  break;
704  case 1:
705  for(i = 0; i < GetCoordim(); ++i)
706  {
707  Vmath::Fill(nqe,df[2*i+1][0] + df[2*i][0],normal[i],1);
708  }
709  break;
710  case 2:
711  for(i = 0; i < GetCoordim(); ++i)
712  {
713  Vmath::Fill(nqe,-df[2*i][0],normal[i],1);
714  }
715  break;
716  default:
717  ASSERTL0(false,"Edge is out of range (edge < 3)");
718  }
719 
720  // normalise
721  fac = 0.0;
722  for(i =0 ; i < GetCoordim(); ++i)
723  {
724  fac += normal[i][0]*normal[i][0];
725  }
726  fac = 1.0/sqrt(fac);
727  for (i = 0; i < GetCoordim(); ++i)
728  {
729  Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
730  }
731  }
732  else // Set up deformed normals
733  {
734  int j;
735 
736  int nquad0 = ptsKeys[0].GetNumPoints();
737  int nquad1 = ptsKeys[1].GetNumPoints();
738 
739  LibUtilities::PointsKey from_key;
740 
741  Array<OneD,NekDouble> normals(GetCoordim()*max(nquad0,nquad1),0.0);
742  Array<OneD,NekDouble> edgejac(GetCoordim()*max(nquad0,nquad1),0.0);
743 
744  // Extract Jacobian along edges and recover local
745  // derivates (dx/dr) for polynomial interpolation by
746  // multiplying m_gmat by jacobian
747  switch(edge)
748  {
749  case 0:
750  for(j = 0; j < nquad0; ++j)
751  {
752  edgejac[j] = jac[j];
753  for(i = 0; i < GetCoordim(); ++i)
754  {
755  normals[i*nquad0+j] = -df[2*i+1][j]*edgejac[j];
756  }
757  }
758  from_key = ptsKeys[0];
759  break;
760  case 1:
761  for(j = 0; j < nquad1; ++j)
762  {
763  edgejac[j] = jac[nquad0*j+nquad0-1];
764  for(i = 0; i < GetCoordim(); ++i)
765  {
766  normals[i*nquad1+j] = (df[2*i][nquad0*j + nquad0-1] + df[2*i+1][nquad0*j + nquad0-1])*edgejac[j];
767  }
768  }
769  from_key = ptsKeys[1];
770  break;
771  case 2:
772  for(j = 0; j < nquad1; ++j)
773  {
774  edgejac[j] = jac[nquad0*j];
775  for(i = 0; i < GetCoordim(); ++i)
776  {
777  normals[i*nquad1+j] = -df[2*i][nquad0*j]*edgejac[j];
778  }
779  }
780  from_key = ptsKeys[1];
781  break;
782  default:
783  ASSERTL0(false,"edge is out of range (edge < 3)");
784 
785  }
786 
787  int nq = from_key.GetNumPoints();
788  Array<OneD,NekDouble> work(nqe,0.0);
789 
790  // interpolate Jacobian and invert
791  LibUtilities::Interp1D(from_key,jac,m_base[0]->GetPointsKey(),work);
792  Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
793 
794  // interpolate
795  for(i = 0; i < GetCoordim(); ++i)
796  {
797  LibUtilities::Interp1D(from_key,&normals[i*nq],m_base[0]->GetPointsKey(),&normal[i][0]);
798  Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
799  }
800 
801  //normalise normal vectors
802  Vmath::Zero(nqe,work,1);
803  for(i = 0; i < GetCoordim(); ++i)
804  {
805  Vmath::Vvtvp(nqe,normal[i],1, normal[i],1,work,1,work,1);
806  }
807 
808  Vmath::Vsqrt(nqe,work,1,work,1);
809  Vmath::Sdiv(nqe,1.0,work,1,work,1);
810 
811  for(i = 0; i < GetCoordim(); ++i)
812  {
813  Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
814  }
815 
816  // Reverse direction so that points are in
817  // anticlockwise direction if edge >=2
818  if(edge >= 2)
819  {
820  for(i = 0; i < GetCoordim(); ++i)
821  {
822  Vmath::Reverse(nqe,normal[i],1, normal[i],1);
823  }
824  }
825  }
826  }
827 
828  }//end of namespace
829 }//end of namespace
Generated on Sun Mar 15 2015 00:11:33 for Nektar++ by   doxygen 1.8.1.2