Nektar++
ExpList2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList2D.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Expansion list 2D definition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 
37 #include <boost/core/ignore_unused.hpp>
38 
39 #include <LocalRegions/TriExp.h>
40 #include <LocalRegions/QuadExp.h>
41 #include <LocalRegions/NodalTriExp.h>
42 #include <LocalRegions/Expansion3D.h>
43 #include <MultiRegions/ExpList2D.h>
44 #include <LibUtilities/Foundations/Interp.h>
45 #include <LibUtilities/Foundations/PhysGalerkinProject.h>
46 
47 using namespace std;
48 
49 namespace Nektar
50 {
51  namespace MultiRegions
52  {
53  /**
54  * @class ExpList2D
55  *
56  * This multi-elemental expansion, which does not exhibit any coupling
57  * between the expansion on the separate elements, can be formulated
58  * as,
59  * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
60  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
61  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
62  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
63  * This class inherits all its variables and member functions from the
64  * base class #ExpList.
65  */
66 
67  /**
68  *
69  */
70  ExpList2D::ExpList2D():
71  ExpList()
72  {
73  SetExpType(e2D);
74  }
75 
76 
77  /**
78  *
79  */
80  ExpList2D::~ExpList2D()
81  {
82  }
83 
84 
85  /**
86  * @param In ExpList2D object to copy.
87  */
88  ExpList2D::ExpList2D(
89  const ExpList2D &In,
90  const bool DeclareCoeffPhysArrays):
91  ExpList(In,DeclareCoeffPhysArrays)
92  {
93  SetExpType(e2D);
94  }
95 
96  /**
97  * @param In ExpList2D object to copy.
98  * @param eIDs Id of elements that should be copied.
99  */
100  ExpList2D::ExpList2D(
101  const ExpList2D &In,
102  const std::vector<unsigned int> &eIDs,
103  const bool DeclareCoeffPhysArrays,
104  const Collections::ImplementationType ImpType):
105  ExpList(In,eIDs,DeclareCoeffPhysArrays)
106  {
107  SetExpType(e2D);
108 
109  // Setup Default optimisation information.
110  int nel = GetExpSize();
111  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
112  ::AllocateSharedPtr(nel);
113 
114  // set up offset arrays.
115  SetCoeffPhysOffsets();
116 
117  ReadGlobalOptimizationParameters();
118  CreateCollections(ImpType);
119  }
120 
121 
122  /**
123  * Given a mesh \a graph2D, containing information about the domain and
124  * the spectral/hp element expansion, this constructor fills the list
125  * of local expansions \texttt{m_exp} with the proper expansions,
126  * calculates the total number of quadrature points
127  * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
128  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays #m_coeffs
129  * and #m_phys.
130  *
131  * @param graph2D A mesh, containing information about the domain
132  * and the spectral/hp element expansion.
133  */
134  ExpList2D::ExpList2D(
135  const LibUtilities::SessionReaderSharedPtr &pSession,
136  const SpatialDomains::MeshGraphSharedPtr &graph2D,
137  const bool DeclareCoeffPhysArrays,
138  const std::string &var,
139  const Collections::ImplementationType ImpType):
140  ExpList(pSession,graph2D)
141  {
142  SetExpType(e2D);
143 
144  int elmtid=0;
145  LocalRegions::TriExpSharedPtr tri;
146  LocalRegions::NodalTriExpSharedPtr Ntri;
147  LibUtilities::PointsType TriNb;
148  LocalRegions::QuadExpSharedPtr quad;
149  SpatialDomains::Composite comp;
150 
151  const SpatialDomains::ExpansionMap &expansions
152  = graph2D->GetExpansions(var);
153 
154  for (auto &expIt : expansions)
155  {
156  SpatialDomains::TriGeomSharedPtr TriangleGeom;
157  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
158 
159  if ((TriangleGeom = std::dynamic_pointer_cast<SpatialDomains
160  ::TriGeom>(expIt.second->m_geomShPtr)))
161  {
162  LibUtilities::BasisKey TriBa
163  = expIt.second->m_basisKeyVector[0];
164  LibUtilities::BasisKey TriBb
165  = expIt.second->m_basisKeyVector[1];
166 
167  // This is not elegantly implemented needs re-thinking.
168  if (TriBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
169  {
170  LibUtilities::BasisKey newBa(LibUtilities::eOrtho_A,
171  TriBa.GetNumModes(),
172  TriBa.GetPointsKey());
173 
174  TriNb = LibUtilities::eNodalTriElec;
175  Ntri = MemoryManager<LocalRegions::NodalTriExp>
176  ::AllocateSharedPtr(newBa,TriBb,TriNb,
177  TriangleGeom);
178  Ntri->SetElmtId(elmtid++);
179  (*m_exp).push_back(Ntri);
180  }
181  else
182  {
183  tri = MemoryManager<LocalRegions::TriExp>
184  ::AllocateSharedPtr(TriBa,TriBb,
185  TriangleGeom);
186  tri->SetElmtId(elmtid++);
187  (*m_exp).push_back(tri);
188  }
189  m_ncoeffs += (TriBa.GetNumModes()*(TriBa.GetNumModes()+1))/2
190  + TriBa.GetNumModes()*(TriBb.GetNumModes()
191  -TriBa.GetNumModes());
192  m_npoints += TriBa.GetNumPoints()*TriBb.GetNumPoints();
193  }
194  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
195  SpatialDomains::QuadGeom>(expIt.second->m_geomShPtr)))
196  {
197  LibUtilities::BasisKey QuadBa
198  = expIt.second->m_basisKeyVector[0];
199  LibUtilities::BasisKey QuadBb
200  = expIt.second->m_basisKeyVector[1];
201 
202  quad = MemoryManager<LocalRegions::QuadExp>
203  ::AllocateSharedPtr(QuadBa,QuadBb,
204  QuadrilateralGeom);
205  quad->SetElmtId(elmtid++);
206  (*m_exp).push_back(quad);
207 
208  m_ncoeffs += QuadBa.GetNumModes()*QuadBb.GetNumModes();
209  m_npoints += QuadBa.GetNumPoints()*QuadBb.GetNumPoints();
210  }
211  else
212  {
213  ASSERTL0(false, "dynamic cast to a proper Geometry2D "
214  "failed");
215  }
216 
217  }
218 
219  // set up element numbering
220  for(int i = 0; i < (*m_exp).size(); ++i)
221  {
222  (*m_exp)[i]->SetElmtId(i);
223  }
224 
225  // Setup Default optimisation information.
226  int nel = GetExpSize();
227  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
228  ::AllocateSharedPtr(nel);
229 
230 
231  // set up offset arrays.
232  SetCoeffPhysOffsets();
233 
234  if (DeclareCoeffPhysArrays)
235  {
236  // Set up m_coeffs, m_phys.
237  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
238  m_phys = Array<OneD, NekDouble>(m_npoints);
239  }
240 
241  ReadGlobalOptimizationParameters();
242  CreateCollections(ImpType);
243  }
244 
245 
246  /**
247  * Given an expansion vector \a expansions, containing
248  * information about the domain and the spectral/hp element
249  * expansion, this constructor fills the list of local
250  * expansions \texttt{m_exp} with the proper expansions,
251  * calculates the total number of quadrature points
252  * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
253  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays
254  * #m_coeffs and #m_phys.
255  *
256  * @param expansions A vector containing information about the
257  * domain and the spectral/hp element
258  * expansion.
259  */
260  ExpList2D::ExpList2D(
261  const LibUtilities::SessionReaderSharedPtr &pSession,
262  const SpatialDomains::ExpansionMap &expansions,
263  const bool DeclareCoeffPhysArrays,
264  const Collections::ImplementationType ImpType):
265  ExpList(pSession)
266  {
267  SetExpType(e2D);
268 
269  int elmtid=0;
270  LocalRegions::TriExpSharedPtr tri;
271  LocalRegions::NodalTriExpSharedPtr Ntri;
272  LibUtilities::PointsType TriNb;
273  LocalRegions::QuadExpSharedPtr quad;
274  SpatialDomains::Composite comp;
275 
276  for (auto &expIt : expansions)
277  {
278  SpatialDomains::TriGeomSharedPtr TriangleGeom;
279  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
280 
281  if ((TriangleGeom = std::dynamic_pointer_cast<SpatialDomains
282  ::TriGeom>(expIt.second->m_geomShPtr)))
283  {
284  LibUtilities::BasisKey TriBa
285  = expIt.second->m_basisKeyVector[0];
286  LibUtilities::BasisKey TriBb
287  = expIt.second->m_basisKeyVector[1];
288 
289  // This is not elegantly implemented needs re-thinking.
290  if (TriBa.GetBasisType() == LibUtilities::eGLL_Lagrange)
291  {
292  LibUtilities::BasisKey newBa(LibUtilities::eOrtho_A,
293  TriBa.GetNumModes(),
294  TriBa.GetPointsKey());
295 
296  TriNb = LibUtilities::eNodalTriElec;
297  Ntri = MemoryManager<LocalRegions::NodalTriExp>
298  ::AllocateSharedPtr(newBa,TriBb,TriNb,
299  TriangleGeom);
300  Ntri->SetElmtId(elmtid++);
301  (*m_exp).push_back(Ntri);
302  }
303  else
304  {
305  tri = MemoryManager<LocalRegions::TriExp>
306  ::AllocateSharedPtr(TriBa,TriBb,
307  TriangleGeom);
308  tri->SetElmtId(elmtid++);
309  (*m_exp).push_back(tri);
310  }
311  m_ncoeffs += (TriBa.GetNumModes()*(TriBa.GetNumModes()+1))/2
312  + TriBa.GetNumModes()*(TriBb.GetNumModes()
313  -TriBa.GetNumModes());
314  m_npoints += TriBa.GetNumPoints()*TriBb.GetNumPoints();
315  }
316  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
317  SpatialDomains::QuadGeom>(expIt.second->m_geomShPtr)))
318  {
319  LibUtilities::BasisKey QuadBa
320  = expIt.second->m_basisKeyVector[0];
321  LibUtilities::BasisKey QuadBb
322  = expIt.second->m_basisKeyVector[1];
323 
324  quad = MemoryManager<LocalRegions::QuadExp>
325  ::AllocateSharedPtr(QuadBa,QuadBb,
326  QuadrilateralGeom);
327  quad->SetElmtId(elmtid++);
328  (*m_exp).push_back(quad);
329 
330  m_ncoeffs += QuadBa.GetNumModes()*QuadBb.GetNumModes();
331  m_npoints += QuadBa.GetNumPoints()*QuadBb.GetNumPoints();
332  }
333  else
334  {
335  ASSERTL0(false, "dynamic cast to a proper Geometry2D "
336  "failed");
337  }
338 
339  }
340 
341  // Setup Default optimisation information.
342  int nel = GetExpSize();
343  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
344  ::AllocateSharedPtr(nel);
345 
346 
347  // set up offset arrays.
348  SetCoeffPhysOffsets();
349 
350  if (DeclareCoeffPhysArrays)
351  {
352  // Set up m_coeffs, m_phys.
353  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
354  m_phys = Array<OneD, NekDouble>(m_npoints);
355  }
356 
357  ReadGlobalOptimizationParameters();
358  CreateCollections(ImpType);
359  }
360 
361 
362  /**
363  * Given a mesh \a graph2D, containing information about the domain and
364  * the a list of basiskeys, this constructor fills the list
365  * of local expansions \texttt{m_exp} with the proper expansions,
366  * calculates the total number of quadrature points
367  * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
368  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays #m_coeffs
369  * and #m_phys.
370  *
371  * @param TriBa A BasisKey, containing the definition of the
372  * basis in the first coordinate direction for
373  * triangular elements
374  * @param TriBb A BasisKey, containing the definition of the
375  * basis in the second coordinate direction for
376  * triangular elements
377  * @param QuadBa A BasisKey, containing the definition of the
378  * basis in the first coordinate direction for
379  * quadrilateral elements
380  * @param QuadBb A BasisKey, containing the definition of the
381  * basis in the second coordinate direction for
382  * quadrilateral elements
383  * @param graph2D A mesh, containing information about the domain
384  * and the spectral/hp element expansion.
385  * @param TriNb The PointsType of possible nodal points
386  */
387  ExpList2D::ExpList2D(
388  const LibUtilities::SessionReaderSharedPtr &pSession,
389  const LibUtilities::BasisKey &TriBa,
390  const LibUtilities::BasisKey &TriBb,
391  const LibUtilities::BasisKey &QuadBa,
392  const LibUtilities::BasisKey &QuadBb,
393  const SpatialDomains::MeshGraphSharedPtr &graph2D,
394  const LibUtilities::PointsType TriNb,
395  const Collections::ImplementationType ImpType):
396  ExpList(pSession, graph2D)
397  {
398  SetExpType(e2D);
399 
400  int elmtid=0;
401  LocalRegions::TriExpSharedPtr tri;
402  LocalRegions::NodalTriExpSharedPtr Ntri;
403  LocalRegions::QuadExpSharedPtr quad;
404  SpatialDomains::Composite comp;
405 
406  const SpatialDomains::ExpansionMap &expansions =
407  graph2D->GetExpansions();
408  m_ncoeffs = 0;
409  m_npoints = 0;
410 
411  m_physState = false;
412 
413  for (auto &expIt : expansions)
414  {
415  SpatialDomains::TriGeomSharedPtr TriangleGeom;
416  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
417 
418  if ((TriangleGeom = std::dynamic_pointer_cast<SpatialDomains::
419  TriGeom>(expIt.second->m_geomShPtr)))
420  {
421  if (TriNb < LibUtilities::SIZE_PointsType)
422  {
423  Ntri = MemoryManager<LocalRegions::NodalTriExp>::
424  AllocateSharedPtr(TriBa, TriBb, TriNb,
425  TriangleGeom);
426  Ntri->SetElmtId(elmtid++);
427  (*m_exp).push_back(Ntri);
428  }
429  else
430  {
431  tri = MemoryManager<LocalRegions::TriExp>::
432  AllocateSharedPtr(TriBa, TriBb, TriangleGeom);
433  tri->SetElmtId(elmtid++);
434  (*m_exp).push_back(tri);
435  }
436 
437  m_ncoeffs += (TriBa.GetNumModes()*(TriBa.GetNumModes()+1))/2
438  + TriBa.GetNumModes() * (TriBb.GetNumModes() -
439  TriBa.GetNumModes());
440  m_npoints += TriBa.GetNumPoints()*TriBb.GetNumPoints();
441  }
442  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
443  SpatialDomains::QuadGeom>(expIt.second->m_geomShPtr)))
444  {
445  quad = MemoryManager<LocalRegions::QuadExp>::
446  AllocateSharedPtr(QuadBa, QuadBb, QuadrilateralGeom);
447  quad->SetElmtId(elmtid++);
448  (*m_exp).push_back(quad);
449 
450  m_ncoeffs += QuadBa.GetNumModes()*QuadBb.GetNumModes();
451  m_npoints += QuadBa.GetNumPoints()*QuadBb.GetNumPoints();
452  }
453  else
454  {
455  ASSERTL0(false,
456  "dynamic cast to a proper Geometry2D failed");
457  }
458 
459  }
460 
461  // Setup Default optimisation information.
462  int nel = GetExpSize();
463  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
464  ::AllocateSharedPtr(nel);
465 
466  // Set up m_coeffs, m_phys and offset arrays.
467  SetCoeffPhysOffsets();
468  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
469  m_phys = Array<OneD, NekDouble>(m_npoints);
470 
471  ReadGlobalOptimizationParameters();
472  CreateCollections(ImpType);
473  }
474 
475  /**
476  * Specialized constructor for trace expansions. Store
477  * expansions for the trace space used in DisContField3D
478  *
479  * @param bndConstraint Array of ExpList2D objects each containing a
480  * 2D spectral/hp element expansion on a single
481  * boundary region.
482  * @param bndCond Array of BoundaryCondition objects which contain
483  * information about the boundary conditions on the
484  * different boundary regions.
485  * @param locexp Complete domain expansion list.
486  * @param graph3D 3D mesh corresponding to the expansion list.
487  * @param periodicFaces List of periodic faces.
488  * @param DeclareCoeffPhysArrays If true, set up m_coeffs,
489  * m_phys arrays
490  **/
491  ExpList2D::ExpList2D(
492  const LibUtilities::SessionReaderSharedPtr &pSession,
493  const Array<OneD,const ExpListSharedPtr> &bndConstraint,
494  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr> &bndCond,
495  const LocalRegions::ExpansionVector &locexp,
496  const SpatialDomains::MeshGraphSharedPtr &graph3D,
497  const PeriodicMap &periodicFaces,
498  const bool DeclareCoeffPhysArrays,
499  const std::string variable,
500  const Collections::ImplementationType ImpType):
501  ExpList(pSession, graph3D)
502  {
503  boost::ignore_unused(periodicFaces, variable);
504 
505  SetExpType(e2D);
506 
507  int i, j, id, elmtid=0;
508  set<int> facesDone;
509 
510  SpatialDomains::Geometry2DSharedPtr FaceGeom;
511  SpatialDomains::QuadGeomSharedPtr FaceQuadGeom;
512  SpatialDomains::TriGeomSharedPtr FaceTriGeom;
513  LocalRegions::QuadExpSharedPtr FaceQuadExp;
514  LocalRegions::TriExpSharedPtr FaceTriExp;
515  LocalRegions::Expansion2DSharedPtr exp2D;
516  LocalRegions::Expansion3DSharedPtr exp3D;
517 
518  // First loop over boundary conditions to renumber
519  // Dirichlet boundaries
520  for (i = 0; i < bndCond.num_elements(); ++i)
521  {
522  if (bndCond[i]->GetBoundaryConditionType()
523  == SpatialDomains::eDirichlet)
524  {
525  for (j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
526  {
527  LibUtilities::BasisKey bkey0 = bndConstraint[i]
528  ->GetExp(j)->GetBasis(0)->GetBasisKey();
529  LibUtilities::BasisKey bkey1 = bndConstraint[i]
530  ->GetExp(j)->GetBasis(1)->GetBasisKey();
531  exp2D = bndConstraint[i]->GetExp(j)
532  ->as<LocalRegions::Expansion2D>();
533  FaceGeom = exp2D->GetGeom2D();
534 
535  //if face is a quad
536  if((FaceQuadGeom = std::dynamic_pointer_cast<
537  SpatialDomains::QuadGeom>(FaceGeom)))
538  {
539  FaceQuadExp = MemoryManager<LocalRegions::QuadExp>
540  ::AllocateSharedPtr(bkey0, bkey1, FaceQuadGeom);
541  facesDone.insert(FaceQuadGeom->GetGlobalID());
542  FaceQuadExp->SetElmtId(elmtid++);
543  (*m_exp).push_back(FaceQuadExp);
544  }
545  //if face is a triangle
546  else if((FaceTriGeom = std::dynamic_pointer_cast<
547  SpatialDomains::TriGeom>(FaceGeom)))
548  {
549  FaceTriExp = MemoryManager<LocalRegions::TriExp>
550  ::AllocateSharedPtr(bkey0, bkey1, FaceTriGeom);
551  facesDone.insert(FaceTriGeom->GetGlobalID());
552  FaceTriExp->SetElmtId(elmtid++);
553  (*m_exp).push_back(FaceTriExp);
554  }
555  else
556  {
557  ASSERTL0(false,"dynamic cast to a proper face geometry failed");
558  }
559  }
560  }
561  }
562 
563  map<int, pair<SpatialDomains::Geometry2DSharedPtr,
564  pair<LibUtilities::BasisKey,
565  LibUtilities::BasisKey> > > faceOrders;
566 
567  for(i = 0; i < locexp.size(); ++i)
568  {
569  exp3D = locexp[i]->as<LocalRegions::Expansion3D>();
570  for (j = 0; j < exp3D->GetNfaces(); ++j)
571  {
572  FaceGeom = exp3D->GetGeom3D()->GetFace(j);
573  id = FaceGeom->GetGlobalID();
574 
575  if(facesDone.count(id) != 0)
576  {
577  continue;
578  }
579  auto it = faceOrders.find(id);
580 
581  if (it == faceOrders.end())
582  {
583  LibUtilities::BasisKey face_dir0
584  = locexp[i]->DetFaceBasisKey(j,0);
585  LibUtilities::BasisKey face_dir1
586  = locexp[i]->DetFaceBasisKey(j,1);
587 
588  faceOrders.insert(
589  std::make_pair(
590  id, std::make_pair(
591  FaceGeom,
592  std::make_pair(face_dir0, face_dir1))));
593  }
594  else // variable modes/points
595  {
596  LibUtilities::BasisKey face0 =
597  locexp[i]->DetFaceBasisKey(j,0);
598  LibUtilities::BasisKey face1 =
599  locexp[i]->DetFaceBasisKey(j,1);
600  LibUtilities::BasisKey existing0 =
601  it->second.second.first;
602  LibUtilities::BasisKey existing1 =
603  it->second.second.second;
604 
605  int np11 = face0 .GetNumPoints();
606  int np12 = face1 .GetNumPoints();
607  int np21 = existing0.GetNumPoints();
608  int np22 = existing1.GetNumPoints();
609  int nm11 = face0 .GetNumModes ();
610  int nm12 = face1 .GetNumModes ();
611  int nm21 = existing0.GetNumModes ();
612  int nm22 = existing1.GetNumModes ();
613 
614  if ((np22 >= np12 || np21 >= np11) &&
615  (nm22 >= nm12 || nm21 >= nm11))
616  {
617  continue;
618  }
619  else if((np22 < np12 || np21 < np11) &&
620  (nm22 < nm12 || nm21 < nm11))
621  {
622  it->second.second.first = face0;
623  it->second.second.second = face1;
624  }
625  else
626  {
627  ASSERTL0(false,
628  "inappropriate number of points/modes (max "
629  "num of points is not set with max order)");
630  }
631  }
632  }
633  }
634 
635  LibUtilities::CommSharedPtr vComm = pSession->GetComm();
636  int nproc = vComm->GetSize(); // number of processors
637  int facepr = vComm->GetRank(); // ID processor
638 
639  if (nproc > 1)
640  {
641  int fCnt = 0;
642 
643  // Count the number of faces on each partition
644  for(i = 0; i < locexp.size(); ++i)
645  {
646  fCnt += locexp[i]->GetNfaces();
647  }
648 
649  // Set up the offset and the array that will contain the list of
650  // face IDs, then reduce this across processors.
651  Array<OneD, int> faceCnt(nproc,0);
652  faceCnt[facepr] = fCnt;
653  vComm->AllReduce(faceCnt, LibUtilities::ReduceSum);
654 
655  int totFaceCnt = Vmath::Vsum(nproc, faceCnt, 1);
656  Array<OneD, int> fTotOffsets(nproc,0);
657 
658  for (i = 1; i < nproc; ++i)
659  {
660  fTotOffsets[i] = fTotOffsets[i-1] + faceCnt[i-1];
661  }
662 
663  // Local list of the edges per element
664 
665  Array<OneD, int> FacesTotID (totFaceCnt, 0);
666  Array<OneD, int> FacesTotNm0 (totFaceCnt, 0);
667  Array<OneD, int> FacesTotNm1 (totFaceCnt, 0);
668  Array<OneD, int> FacesTotPnts0(totFaceCnt, 0);
669  Array<OneD, int> FacesTotPnts1(totFaceCnt, 0);
670 
671  int cntr = fTotOffsets[facepr];
672 
673  for(i = 0; i < locexp.size(); ++i)
674  {
675  exp3D = locexp[i]->as<LocalRegions::Expansion3D>();
676 
677  int nfaces = locexp[i]->GetNfaces();
678 
679  for(j = 0; j < nfaces; ++j, ++cntr)
680  {
681  LibUtilities::BasisKey face_dir0
682  = locexp[i]->DetFaceBasisKey(j,0);
683  LibUtilities::BasisKey face_dir1
684  = locexp[i]->DetFaceBasisKey(j,1);
685 
686  FacesTotID[cntr] = exp3D->GetGeom3D()->GetFid(j);
687  FacesTotNm0[cntr] = face_dir0.GetNumModes ();
688  FacesTotNm1[cntr] = face_dir1.GetNumModes ();
689  FacesTotPnts0[cntr] = face_dir0.GetNumPoints();
690  FacesTotPnts1[cntr] = face_dir1.GetNumPoints();
691  }
692  }
693 
694  vComm->AllReduce(FacesTotID, LibUtilities::ReduceSum);
695  vComm->AllReduce(FacesTotNm0, LibUtilities::ReduceSum);
696  vComm->AllReduce(FacesTotNm1, LibUtilities::ReduceSum);
697  vComm->AllReduce(FacesTotPnts0, LibUtilities::ReduceSum);
698  vComm->AllReduce(FacesTotPnts1, LibUtilities::ReduceSum);
699 
700  for (i = 0; i < totFaceCnt; ++i)
701  {
702  auto it = faceOrders.find(FacesTotID[i]);
703 
704  if (it == faceOrders.end())
705  {
706  continue;
707  }
708 
709  LibUtilities::BasisKey existing0 =
710  it->second.second.first;
711  LibUtilities::BasisKey existing1 =
712  it->second.second.second;
713  LibUtilities::BasisKey face0(
714  existing0.GetBasisType(), FacesTotNm0[i],
715  LibUtilities::PointsKey(FacesTotPnts0[i],
716  existing0.GetPointsType()));
717  LibUtilities::BasisKey face1(
718  existing1.GetBasisType(), FacesTotNm1[i],
719  LibUtilities::PointsKey(FacesTotPnts1[i],
720  existing1.GetPointsType()));
721 
722  int np11 = face0 .GetNumPoints();
723  int np12 = face1 .GetNumPoints();
724  int np21 = existing0.GetNumPoints();
725  int np22 = existing1.GetNumPoints();
726  int nm11 = face0 .GetNumModes ();
727  int nm12 = face1 .GetNumModes ();
728  int nm21 = existing0.GetNumModes ();
729  int nm22 = existing1.GetNumModes ();
730 
731  if ((np22 >= np12 || np21 >= np11) &&
732  (nm22 >= nm12 || nm21 >= nm11))
733  {
734  continue;
735  }
736  else if((np22 < np12 || np21 < np11) &&
737  (nm22 < nm12 || nm21 < nm11))
738  {
739  it->second.second.first = face0;
740  it->second.second.second = face1;
741  }
742  else
743  {
744  ASSERTL0(false,
745  "inappropriate number of points/modes (max "
746  "num of points is not set with max order)");
747  }
748  }
749  }
750 
751  for (auto &it : faceOrders)
752  {
753  FaceGeom = it.second.first;
754 
755  if ((FaceQuadGeom = std::dynamic_pointer_cast<
756  SpatialDomains::QuadGeom>(FaceGeom)))
757  {
758  FaceQuadExp = MemoryManager<LocalRegions::QuadExp>
759  ::AllocateSharedPtr(it.second.second.first,
760  it.second.second.second,
761  FaceQuadGeom);
762  FaceQuadExp->SetElmtId(elmtid++);
763  (*m_exp).push_back(FaceQuadExp);
764  }
765  else if ((FaceTriGeom = std::dynamic_pointer_cast<
766  SpatialDomains::TriGeom>(FaceGeom)))
767  {
768  FaceTriExp = MemoryManager<LocalRegions::TriExp>
769  ::AllocateSharedPtr(it.second.second.first,
770  it.second.second.second,
771  FaceTriGeom);
772  FaceTriExp->SetElmtId(elmtid++);
773  (*m_exp).push_back(FaceTriExp);
774  }
775  }
776 
777  // Setup Default optimisation information.
778  int nel = GetExpSize();
779 
780  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
781  ::AllocateSharedPtr(nel);
782 
783  // Set up offset information and array sizes
784  SetCoeffPhysOffsets();
785 
786  // Set up m_coeffs, m_phys.
787  if(DeclareCoeffPhysArrays)
788  {
789  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
790  m_phys = Array<OneD, NekDouble>(m_npoints);
791  }
792 
793  CreateCollections(ImpType);
794  }
795 
796  /**
797  * Fills the list of local expansions with the segments from the 3D
798  * mesh specified by \a domain. This CompositeMap contains a list of
799  * Composites which define the Neumann boundary.
800  * @see ExpList2D#ExpList2D(SpatialDomains::MeshGraph2D&)
801  * for details.
802  * @param domain A domain, comprising of one or more composite
803  * regions.
804  * @param graph3D A mesh, containing information about the domain
805  * and the spectral/hp element expansions.
806  */
807  ExpList2D::ExpList2D(
808  const LibUtilities::SessionReaderSharedPtr &pSession,
809  const SpatialDomains::CompositeMap &domain,
810  const SpatialDomains::MeshGraphSharedPtr &graph3D,
811  const std::string variable,
812  const LibUtilities::CommSharedPtr comm,
813  const Collections::ImplementationType ImpType)
814  : ExpList(pSession, graph3D)
815  {
816  SetExpType(e2D);
817 
818  if (comm)
819  {
820  m_comm = comm;
821  }
822 
823  int j, elmtid=0;
824  int nel = 0;
825 
826  SpatialDomains::Composite comp;
827  SpatialDomains::TriGeomSharedPtr TriangleGeom;
828  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
829 
830  LocalRegions::TriExpSharedPtr tri;
831  LocalRegions::NodalTriExpSharedPtr Ntri;
832  LibUtilities::PointsType TriNb;
833  LocalRegions::QuadExpSharedPtr quad;
834 
835  for (auto &compIt : domain)
836  {
837  nel += compIt.second->m_geomVec.size();
838  }
839 
840  for (auto &compIt : domain)
841  {
842  for (j = 0; j < compIt.second->m_geomVec.size(); ++j)
843  {
844  if ((TriangleGeom = std::dynamic_pointer_cast<
845  SpatialDomains::TriGeom>(compIt.second->m_geomVec[j])))
846  {
847  LibUtilities::BasisKey TriBa
848  = graph3D->GetFaceBasisKey(TriangleGeom, 0, variable);
849  LibUtilities::BasisKey TriBb
850  = graph3D->GetFaceBasisKey(TriangleGeom,1,variable);
851 
852  if (graph3D->GetExpansions().begin()->second->
853  m_basisKeyVector[0].GetBasisType() ==
854  LibUtilities::eGLL_Lagrange)
855  {
856  ASSERTL0(false,"This method needs sorting");
857  TriNb = LibUtilities::eNodalTriElec;
858 
859  Ntri = MemoryManager<LocalRegions::NodalTriExp>
860  ::AllocateSharedPtr(TriBa,TriBb,TriNb,
861  TriangleGeom);
862  Ntri->SetElmtId(elmtid++);
863  (*m_exp).push_back(Ntri);
864  }
865  else
866  {
867  tri = MemoryManager<LocalRegions::TriExp>
868  ::AllocateSharedPtr(TriBa, TriBb,
869  TriangleGeom);
870  tri->SetElmtId(elmtid++);
871  (*m_exp).push_back(tri);
872  }
873 
874  m_ncoeffs
875  += (TriBa.GetNumModes()*(TriBa.GetNumModes()+1))/2
876  + TriBa.GetNumModes()*(TriBb.GetNumModes()
877  -TriBa.GetNumModes());
878  m_npoints += TriBa.GetNumPoints()*TriBb.GetNumPoints();
879  }
880  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
881  SpatialDomains::QuadGeom>(compIt.second->m_geomVec[j])))
882  {
883  LibUtilities::BasisKey QuadBa
884  = graph3D->GetFaceBasisKey(QuadrilateralGeom, 0,
885  variable);
886  LibUtilities::BasisKey QuadBb
887  = graph3D->GetFaceBasisKey(QuadrilateralGeom, 1,
888  variable);
889 
890  quad = MemoryManager<LocalRegions::QuadExp>
891  ::AllocateSharedPtr(QuadBa, QuadBb,
892  QuadrilateralGeom);
893  quad->SetElmtId(elmtid++);
894  (*m_exp).push_back(quad);
895 
896  m_ncoeffs += QuadBa.GetNumModes()*QuadBb.GetNumModes();
897  m_npoints += QuadBa.GetNumPoints()
898  * QuadBb.GetNumPoints();
899  }
900  else
901  {
902  ASSERTL0(false,
903  "dynamic cast to a proper Geometry2D failed");
904  }
905  }
906 
907  }
908 
909  // Setup Default optimisation information.
910  nel = GetExpSize();
911  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
912  ::AllocateSharedPtr(nel);
913 
914  // Set up m_coeffs, m_phys and offset arrays.
915  SetCoeffPhysOffsets();
916  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
917  m_phys = Array<OneD, NekDouble>(m_npoints);
918 
919  ReadGlobalOptimizationParameters();
920  CreateCollections(ImpType);
921  }
922 
923  /**
924  * One-dimensional upwind.
925  * @param Vn Velocity field.
926  * @param Fwd Left state.
927  * @param Bwd Right state.
928  * @param Upwind Output vector.
929  */
930  void ExpList2D::v_Upwind(
931  const Array<OneD, const NekDouble> &Vn,
932  const Array<OneD, const NekDouble> &Fwd,
933  const Array<OneD, const NekDouble> &Bwd,
934  Array<OneD, NekDouble> &Upwind)
935  {
936  int i,j,f_npoints,offset;
937 
938  // Process each expansion.
939  for(i = 0; i < m_exp->size(); ++i)
940  {
941  // Get the number of points and the data offset.
942  f_npoints = (*m_exp)[i]->GetNumPoints(0)*
943  (*m_exp)[i]->GetNumPoints(1);
944  offset = m_phys_offset[i];
945 
946  // Process each point in the expansion.
947  for(j = 0; j < f_npoints; ++j)
948  {
949  // Upwind based on one-dimensional velocity.
950  if(Vn[offset + j] > 0.0)
951  {
952  Upwind[offset + j] = Fwd[offset + j];
953  }
954  else
955  {
956  Upwind[offset + j] = Bwd[offset + j];
957  }
958  }
959  }
960  }
961 
962  /**
963  * @brief Helper function to re-align face to a given orientation.
964  */
965  void AlignFace(const StdRegions::Orientation orient,
966  const int nquad1,
967  const int nquad2,
968  const Array<OneD, const NekDouble> &in,
969  Array<OneD, NekDouble> &out)
970  {
971  // Copy transpose.
972  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
973  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
974  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
975  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
976  {
977  for (int i = 0; i < nquad2; ++i)
978  {
979  for (int j = 0; j < nquad1; ++j)
980  {
981  out[i*nquad1 + j] = in[j*nquad2 + i];
982  }
983  }
984  }
985  else
986  {
987  for (int i = 0; i < nquad2; ++i)
988  {
989  for (int j = 0; j < nquad1; ++j)
990  {
991  out[i*nquad1 + j] = in[i*nquad1 + j];
992  }
993  }
994  }
995 
996  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
997  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
998  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
999  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1000  {
1001  // Reverse x direction
1002  for (int i = 0; i < nquad2; ++i)
1003  {
1004  for (int j = 0; j < nquad1/2; ++j)
1005  {
1006  swap(out[i*nquad1 + j],
1007  out[i*nquad1 + nquad1-j-1]);
1008  }
1009  }
1010  }
1011 
1012  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
1013  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
1014  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
1015  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1016  {
1017  // Reverse y direction
1018  for (int j = 0; j < nquad1; ++j)
1019  {
1020  for (int i = 0; i < nquad2/2; ++i)
1021  {
1022  swap(out[i*nquad1 + j],
1023  out[(nquad2-i-1)*nquad1 + j]);
1024  }
1025  }
1026  }
1027  }
1028 
1029  /**
1030  * @brief For each local element, copy the normals stored in the element
1031  * list into the array \a normals.
1032  *
1033  * @param normals Multidimensional array in which to copy normals
1034  * to. Must have dimension equal to or larger than
1035  * the spatial dimension of the elements.
1036  */
1037  void ExpList2D::v_GetNormals(
1038  Array<OneD, Array<OneD, NekDouble> > &normals)
1039  {
1040  Array<OneD, NekDouble> tmp;
1041  int i, j;
1042  const int coordim = GetCoordim(0);
1043 
1044  ASSERTL1(normals.num_elements() >= coordim,
1045  "Output vector does not have sufficient dimensions to "
1046  "match coordim");
1047 
1048  // Process each expansion.
1049  for (i = 0; i < m_exp->size(); ++i)
1050  {
1051  LocalRegions::Expansion2DSharedPtr traceExp = (*m_exp)[i]->as<
1052  LocalRegions::Expansion2D>();
1053  LocalRegions::Expansion3DSharedPtr exp3D =
1054  traceExp->GetLeftAdjacentElementExp();
1055 
1056  // Get the number of points and normals for this expansion.
1057  int faceNum = traceExp->GetLeftAdjacentElementFace();
1058  int offset = m_phys_offset[i];
1059 
1060  const Array<OneD, const Array<OneD, NekDouble> > &locNormals
1061  = exp3D->GetFaceNormal(faceNum);
1062 
1063  // Project normals from 3D element onto the same orientation as
1064  // the trace expansion.
1065  StdRegions::Orientation orient = exp3D->GetForient(faceNum);
1066 
1067 
1068  int fromid0,fromid1;
1069 
1070  if(orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
1071  {
1072  fromid0 = 0;
1073  fromid1 = 1;
1074  }
1075  else
1076  {
1077  fromid0 = 1;
1078  fromid1 = 0;
1079  }
1080 
1081  LibUtilities::BasisKey faceBasis0
1082  = exp3D->DetFaceBasisKey(faceNum, fromid0);
1083  LibUtilities::BasisKey faceBasis1
1084  = exp3D->DetFaceBasisKey(faceNum, fromid1);
1085  LibUtilities::BasisKey traceBasis0
1086  = traceExp->GetBasis(0)->GetBasisKey();
1087  LibUtilities::BasisKey traceBasis1
1088  = traceExp->GetBasis(1)->GetBasisKey();
1089 
1090  const int faceNq0 = faceBasis0.GetNumPoints();
1091  const int faceNq1 = faceBasis1.GetNumPoints();
1092 
1093  for (j = 0; j < coordim; ++j)
1094  {
1095  Array<OneD, NekDouble> traceNormals(faceNq0 * faceNq1);
1096  AlignFace(orient, faceNq0, faceNq1,
1097  locNormals[j], traceNormals);
1098  LibUtilities::Interp2D(
1099  faceBasis0.GetPointsKey(),
1100  faceBasis1.GetPointsKey(),
1101  traceNormals,
1102  traceBasis0.GetPointsKey(),
1103  traceBasis1.GetPointsKey(),
1104  tmp = normals[j]+offset);
1105  }
1106  }
1107  }
1108 
1109  /**
1110  *
1111  */
1112  void ExpList2D::v_SetUpPhysNormals()
1113  {
1114  int i, j;
1115  for (i = 0; i < m_exp->size(); ++i)
1116  {
1117  for (j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
1118  {
1119  (*m_exp)[i]->ComputeEdgeNormal(j);
1120  }
1121  }
1122  }
1123 
1124 
1125  void ExpList2D::v_ReadGlobalOptimizationParameters()
1126  {
1127  Array<OneD, int> NumShape(2,0);
1128 
1129  for(int i = 0; i < GetExpSize(); ++i)
1130  {
1131  if((*m_exp)[i]->DetShapeType() == LibUtilities::eTriangle)
1132  {
1133  NumShape[0] += 1;
1134  }
1135  else // Quadrilateral element
1136  {
1137  NumShape[1] += 1;
1138  }
1139  }
1140 
1141  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
1142  ::AllocateSharedPtr(m_session,2,NumShape);
1143  }
1144 
1145  void ExpList2D::v_WriteVtkPieceHeader(
1146  std::ostream &outfile,
1147  int expansion,
1148  int istrip)
1149  {
1150  boost::ignore_unused(istrip);
1151 
1152  int i,j;
1153  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1154  int nquad1 = (*m_exp)[expansion]->GetNumPoints(1);
1155  int ntot = nquad0*nquad1;
1156  int ntotminus = (nquad0-1)*(nquad1-1);
1157 
1158  Array<OneD,NekDouble> coords[3];
1159  coords[0] = Array<OneD,NekDouble>(ntot,0.0);
1160  coords[1] = Array<OneD,NekDouble>(ntot,0.0);
1161  coords[2] = Array<OneD,NekDouble>(ntot,0.0);
1162  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1163 
1164  outfile << " <Piece NumberOfPoints=\""
1165  << ntot << "\" NumberOfCells=\""
1166  << ntotminus << "\">" << endl;
1167  outfile << " <Points>" << endl;
1168  outfile << " <DataArray type=\"Float64\" "
1169  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1170  outfile << " ";
1171  for (i = 0; i < ntot; ++i)
1172  {
1173  for (j = 0; j < 3; ++j)
1174  {
1175  outfile << setprecision(8) << scientific
1176  << (float)coords[j][i] << " ";
1177  }
1178  outfile << endl;
1179  }
1180  outfile << endl;
1181  outfile << " </DataArray>" << endl;
1182  outfile << " </Points>" << endl;
1183  outfile << " <Cells>" << endl;
1184  outfile << " <DataArray type=\"Int32\" "
1185  << "Name=\"connectivity\" format=\"ascii\">" << endl;
1186  for (i = 0; i < nquad0-1; ++i)
1187  {
1188  for (j = 0; j < nquad1-1; ++j)
1189  {
1190  outfile << j*nquad0 + i << " "
1191  << j*nquad0 + i + 1 << " "
1192  << (j+1)*nquad0 + i + 1 << " "
1193  << (j+1)*nquad0 + i << endl;
1194  }
1195  }
1196  outfile << endl;
1197  outfile << " </DataArray>" << endl;
1198  outfile << " <DataArray type=\"Int32\" "
1199  << "Name=\"offsets\" format=\"ascii\">" << endl;
1200  for (i = 0; i < ntotminus; ++i)
1201  {
1202  outfile << i*4+4 << " ";
1203  }
1204  outfile << endl;
1205  outfile << " </DataArray>" << endl;
1206  outfile << " <DataArray type=\"UInt8\" "
1207  << "Name=\"types\" format=\"ascii\">" << endl;
1208  for (i = 0; i < ntotminus; ++i)
1209  {
1210  outfile << "9 ";
1211  }
1212  outfile << endl;
1213  outfile << " </DataArray>" << endl;
1214  outfile << " </Cells>" << endl;
1215  outfile << " <PointData>" << endl;
1216  }
1217 
1218 
1219  void ExpList2D::v_PhysInterp1DScaled(
1220  const NekDouble scale,
1221  const Array<OneD, NekDouble> &inarray,
1222  Array<OneD, NekDouble> &outarray)
1223  {
1224  int cnt,cnt1;
1225 
1226  cnt = cnt1 = 0;
1227  for(int i = 0; i < GetExpSize(); ++i)
1228  {
1229  // get new points key
1230  int pt0 = (*m_exp)[i]->GetNumPoints(0);
1231  int pt1 = (*m_exp)[i]->GetNumPoints(1);
1232  int npt0 = (int) pt0*scale;
1233  int npt1 = (int) pt1*scale;
1234 
1235  LibUtilities::PointsKey newPointsKey0(npt0,
1236  (*m_exp)[i]->GetPointsType(0));
1237  LibUtilities::PointsKey newPointsKey1(npt1,
1238  (*m_exp)[i]->GetPointsType(1));
1239 
1240  // Interpolate points;
1241  LibUtilities::Interp2D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
1242  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
1243  &inarray[cnt],newPointsKey0,
1244  newPointsKey1,&outarray[cnt1]);
1245 
1246  cnt += pt0*pt1;
1247  cnt1 += npt0*npt1;
1248  }
1249  }
1250 
1251  void ExpList2D::v_PhysGalerkinProjection1DScaled(
1252  const NekDouble scale,
1253  const Array<OneD, NekDouble> &inarray,
1254  Array<OneD, NekDouble> &outarray)
1255  {
1256  int cnt,cnt1;
1257 
1258  cnt = cnt1 = 0;
1259  for(int i = 0; i < GetExpSize(); ++i)
1260  {
1261  // get new points key
1262  int pt0 = (*m_exp)[i]->GetNumPoints(0);
1263  int pt1 = (*m_exp)[i]->GetNumPoints(1);
1264  int npt0 = (int) pt0*scale;
1265  int npt1 = (int) pt1*scale;
1266 
1267  LibUtilities::PointsKey newPointsKey0(npt0,
1268  (*m_exp)[i]->GetPointsType(0));
1269  LibUtilities::PointsKey newPointsKey1(npt1,
1270  (*m_exp)[i]->GetPointsType(1));
1271 
1272  // Project points;
1273  LibUtilities::PhysGalerkinProject2D(
1274  newPointsKey0,
1275  newPointsKey1,
1276  &inarray[cnt],
1277  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
1278  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
1279  &outarray[cnt1]);
1280 
1281  cnt += npt0*npt1;
1282  cnt1 += pt0*pt1;
1283  }
1284 
1285  }
1286 
1287 
1288  } //end of namespace
1289 } //end of namespace
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:133
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:144
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:150
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::LibUtilities::BasisKey::GetPointsType
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:156
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::LocalRegions::Expansion2D::GetLeftAdjacentElementExp
Expansion3DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion2D.h:245
Nektar::LocalRegions::Expansion2D::GetGeom2D
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:291
Nektar::LocalRegions::Expansion3D::GetGeom3D
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
Definition: Expansion3D.h:196
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:161
Nektar::MultiRegions::ExpList2D
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
Definition: ExpList2D.h:58
Nektar::MultiRegions::ExpList2D::~ExpList2D
virtual ~ExpList2D()
Destructor.
Definition: ExpList2D.cpp:80
Nektar::MultiRegions::ExpList2D::v_WriteVtkPieceHeader
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Definition: ExpList2D.cpp:1145
Nektar::MultiRegions::ExpList2D::v_PhysInterp1DScaled
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList2D.cpp:1219
Nektar::MultiRegions::ExpList2D::ExpList2D
ExpList2D()
Default constructor.
Definition: ExpList2D.cpp:70
Nektar::MultiRegions::ExpList2D::v_ReadGlobalOptimizationParameters
virtual void v_ReadGlobalOptimizationParameters()
Definition: ExpList2D.cpp:1125
Nektar::MultiRegions::ExpList2D::v_GetNormals
void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
For each local element, copy the normals stored in the element list into the array normals.
Definition: ExpList2D.cpp:1037
Nektar::MultiRegions::ExpList2D::v_PhysGalerkinProjection1DScaled
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList2D.cpp:1251
Nektar::MultiRegions::ExpList2D::v_Upwind
void v_Upwind(const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Upwind the Fwd and Bwd states based on the one- dimensional normal velocity field given by Vn.
Definition: ExpList2D.cpp:930
Nektar::MultiRegions::ExpList2D::v_SetUpPhysNormals
virtual void v_SetUpPhysNormals()
Set up the normals on each expansion.
Definition: ExpList2D.cpp:1112
Nektar::MultiRegions::ExpList
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:104
Nektar::MultiRegions::ExpList::m_coeffs
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1052
Nektar::MultiRegions::ExpList::CreateCollections
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
Definition: ExpList.cpp:3295
Nektar::MultiRegions::ExpList::Upwind
void Upwind(const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.h:2249
Nektar::MultiRegions::ExpList::GetExpSize
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2170
Nektar::MultiRegions::ExpList::m_physState
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1078
Nektar::MultiRegions::ExpList::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1020
Nektar::MultiRegions::ExpList::m_exp
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1090
Nektar::MultiRegions::ExpList::m_ncoeffs
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1030
Nektar::MultiRegions::ExpList::SetCoeffPhysOffsets
void SetCoeffPhysOffsets()
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
Definition: ExpList.cpp:249
Nektar::MultiRegions::ExpList::m_globalOptParam
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
Definition: ExpList.h:1106
Nektar::MultiRegions::ExpList::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1023
Nektar::MultiRegions::ExpList::m_phys_offset
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1104
Nektar::MultiRegions::ExpList::m_phys
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1069
Nektar::MultiRegions::ExpList::SetExpType
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:279
Nektar::MultiRegions::ExpList::GetCoordim
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:2013
Nektar::StdRegions::StdExpansion::GetNfaces
int GetNfaces() const
This function returns the number of faces of the expansion domain.
Definition: StdExpansion.h:437
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:112
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:115
Nektar::LibUtilities::SIZE_PointsType
@ SIZE_PointsType
Length of enum list.
Definition: PointsType.h:81
Nektar::LibUtilities::eNodalTriElec
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:69
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::LibUtilities::PhysGalerkinProject2D
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: PhysGalerkinProject.cpp:115
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LocalRegions::QuadExpSharedPtr
std::shared_ptr< QuadExp > QuadExpSharedPtr
Definition: QuadExp.h:293
Nektar::LocalRegions::Expansion2DSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
Nektar::LocalRegions::NodalTriExpSharedPtr
std::shared_ptr< NodalTriExp > NodalTriExpSharedPtr
Definition: NodalTriExp.h:383
Nektar::LocalRegions::TriExpSharedPtr
std::shared_ptr< TriExp > TriExpSharedPtr
Definition: TriExp.h:302
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:67
Nektar::LocalRegions::Expansion3DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:49
Nektar::MultiRegions::AlignFace
void AlignFace(const StdRegions::Orientation orient, const int nquad1, const int nquad2, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out)
Helper function to re-align face to a given orientation.
Definition: ExpList2D.cpp:965
Nektar::MultiRegions::PeriodicMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Definition: MultiRegions.hpp:200
Nektar::SpatialDomains::ExpansionMap
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:152
Nektar::SpatialDomains::QuadGeomSharedPtr
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
Nektar::SpatialDomains::CompositeMap
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:137
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:163
Nektar::SpatialDomains::Geometry2DSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
Nektar::SpatialDomains::TriGeomSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:740
Nektar::SpatialDomains::Composite::m_geomVec
std::vector< std::shared_ptr< Geometry > > m_geomVec
Definition: MeshGraph.h:133
Generated on Mon Aug 23 2021 06:57:44 for Nektar++ by doxygen 1.9.1