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