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