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