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