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