Nektar++
ExpList.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Expansion list definition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 
37 #include <boost/core/ignore_unused.hpp>
38 
40 #include <MultiRegions/ExpList.h>
42 
45 #include <LocalRegions/HexExp.h>
46 #include <LocalRegions/MatrixKey.h> // for MatrixKey
48 #include <LocalRegions/PointExp.h>
49 #include <LocalRegions/PrismExp.h>
50 #include <LocalRegions/PyrExp.h>
51 #include <LocalRegions/QuadExp.h>
52 #include <LocalRegions/SegExp.h>
53 #include <LocalRegions/TetExp.h>
54 #include <LocalRegions/TriExp.h>
55 
56 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h> // for AssemblyMapCG, etc
57 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h> // for AssemblyMapCG, etc
58 #include <MultiRegions/GlobalLinSysKey.h> // for GlobalLinSysKey
59 #include <MultiRegions/GlobalMatrix.h> // for GlobalMatrix, etc
60 #include <MultiRegions/GlobalMatrixKey.h> // for GlobalMatrixKey
61 
63 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
69 
71 #include <Collections/Operator.h>
72 
75 
76 using namespace std;
77 
78 namespace Nektar
79 {
80 namespace MultiRegions
81 {
82 /**
83  * @class ExpList
84  * All multi-elemental expansions \f$u^{\delta}(\boldsymbol{x})\f$ can
85  * be considered as the assembly of the various elemental contributions.
86  * On a discrete level, this yields,
87  * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
88  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
89  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
90  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
91  * As it is the lowest level class, it contains the definition of the
92  * common data and common routines to all multi-elemental expansions.
93  *
94  * The class stores a vector of expansions, \a m_exp, (each derived from
95  * StdRegions#StdExpansion) which define the constituent components of
96  * the domain. The coefficients from these expansions are concatenated
97  * in \a m_coeffs, while the expansion evaluated at the quadrature
98  * points is stored in \a m_phys.
99  */
100 
101 /**
102  * Creates an empty expansion list.
103  */
104 ExpList::ExpList(const ExpansionType type)
105  : m_expType(type), m_ncoeffs(0), m_npoints(0), m_physState(false),
106  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
107  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
108  m_WaveSpace(false)
109 {
110 }
111 
112 /*----------------------------------------------------------------*/
113 /* Copy Construtor */
114 /*-----------------------------------------------------------------*/
115 
116 /**
117  * Copies an existing expansion list.
118  * @param in Source expansion list.
119  */
120 ExpList::ExpList(const ExpList &in, const bool DeclareCoeffPhysArrays)
121  : std::enable_shared_from_this<ExpList>(in), m_expType(in.m_expType),
122  m_comm(in.m_comm), m_session(in.m_session), m_graph(in.m_graph),
123  m_ncoeffs(in.m_ncoeffs), m_npoints(in.m_npoints), m_physState(false),
124  m_exp(in.m_exp), m_collections(in.m_collections),
125  m_collectionsDoInit(in.m_collectionsDoInit),
126  m_coll_coeff_offset(in.m_coll_coeff_offset),
127  m_coll_phys_offset(in.m_coll_phys_offset),
128  m_coeff_offset(in.m_coeff_offset), m_phys_offset(in.m_phys_offset),
129  m_blockMat(in.m_blockMat), m_WaveSpace(false)
130 {
131 
132  // Set up m_coeffs, m_phys and offset arrays.
133  // use this to keep memory declaration in one place
134  SetupCoeffPhys(DeclareCoeffPhysArrays, false);
135 }
136 
137 /**
138  * Copies the eIds elements from an existing expansion list.
139  * @param in Source expansion list.
140  * @param in elements that will be in the new exp list.
141  */
142 ExpList::ExpList(const ExpList &in, const std::vector<unsigned int> &eIDs,
143  const bool DeclareCoeffPhysArrays,
144  const Collections::ImplementationType ImpType)
145  : m_expType(in.m_expType), m_comm(in.m_comm), m_session(in.m_session),
146  m_graph(in.m_graph), m_physState(false),
147  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
148  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
149  m_WaveSpace(false)
150 {
151  for (int i = 0; i < eIDs.size(); ++i)
152  {
153  (*m_exp).push_back((*(in.m_exp))[eIDs[i]]);
154  }
155 
156  // Set up m_coeffs, m_phys and offset arrays.
157  SetupCoeffPhys(DeclareCoeffPhysArrays);
158 
159  // set up collections
160  CreateCollections(ImpType);
161 }
162 
163 /**
164  * Given a meshgraph \a graph, containing information about
165  * the domain and the spectral/hp element expansion, this
166  * constructor fills the list of local expansions
167  * \texttt{m_exp} with the proper expansions, calculates the
168  * total number of quadrature points \f$x_i\f$ and local
169  * expansion coefficients \f$\hat{u}^e_n\f$ and allocates
170  * memory for the arrays #m_coeffs and #m_phys.
171  *
172  * @param pSession A session within information about expansion
173  *
174  * @param graph A meshgraph, containing information about the
175  * domain and the spectral/hp element expansion.
176  *
177  * @param DeclareCoeffPhysArrays Declare the coefficient and
178  * phys space arrays
179  *
180  * @param ImpType Detail about the implementation type to use
181  * in operators
182  */
185  const bool DeclareCoeffPhysArrays, const std::string &var,
186  const Collections::ImplementationType ImpType)
187  : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
188  m_physState(false),
189  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
190  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
191  m_WaveSpace(false)
192 {
193  // Retrieve the list of expansions
194  const SpatialDomains::ExpansionInfoMap &expansions =
195  graph->GetExpansionInfo(var);
196 
197  // Initialise Expansionn Vector
198  InitialiseExpVector(expansions);
199 
200  // Setup phys coeff space
201  SetupCoeffPhys(DeclareCoeffPhysArrays);
202 
203  // Initialise collection
204  CreateCollections(ImpType);
205 }
206 
207 /**
208  * Given an expansion vector \a expansions, containing
209  * information about the domain and the spectral/hp element
210  * expansion, this constructor fills the list of local
211  * expansions \texttt{m_exp} with the proper expansions,
212  * calculates the total number of quadrature points
213  * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
214  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays
215  * #m_coeffs and #m_phys.
216  *
217  * @param pSession A session within information about expansion
218  * @param expansions A vector containing information about the
219  * domain and the spectral/hp element
220  * expansion.
221  * @param DeclareCoeffPhysArrays Declare the coefficient and
222  * phys space arrays
223  * @param ImpType Detail about the implementation type to use
224  * in operators
225  */
227  const SpatialDomains::ExpansionInfoMap &expansions,
228  const bool DeclareCoeffPhysArrays,
229  const Collections::ImplementationType ImpType)
230  : m_comm(pSession->GetComm()), m_session(pSession), m_physState(false),
231  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
232  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
233  m_WaveSpace(false)
234 {
235 
236  // Initialise expansion vector
237  InitialiseExpVector(expansions);
238 
239  // Set up m_coeffs, m_phys and offset arrays.
240  SetupCoeffPhys(DeclareCoeffPhysArrays);
241 
242  // Setup Collection
243  CreateCollections(ImpType);
244 }
245 
246 //----------------------------------------------------------------------
247 // 0D Expansion Constructors
248 //----------------------------------------------------------------------
250  : m_expType(e0D), m_ncoeffs(1), m_npoints(1), m_physState(false),
251  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
252  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
253  m_WaveSpace(false)
254 {
257  (*m_exp).push_back(Point);
258 
259  SetupCoeffPhys();
260 }
261 
262 /**
263  * Store expansions for the trace space expansions used in
264  * DisContField
265  *
266  * @param pSession A session within information about expansion
267  * @param bndConstraint Array of ExpList1D objects each containing a
268  * 1D spectral/hp element expansion on a single
269  * boundary region.
270  * @param bndCond Array of BoundaryCondition objects which contain
271  * information about the boundary conditions on the
272  * different boundary regions.
273  * @param locexp Complete domain expansion list.
274  * @param graph mesh corresponding to the expansion list.
275  * @param DeclareCoeffPhysArrays Declare the coefficient and
276  * phys space arrays
277  * @param variable The variable name associated with the expansion
278  * @param ImpType Detail about the implementation type to use
279  * in operators
280  */
282  const LibUtilities::SessionReaderSharedPtr &pSession,
283  const Array<OneD, const ExpListSharedPtr> &bndConstraint,
285  const LocalRegions::ExpansionVector &locexp,
287  const LibUtilities::CommSharedPtr &comm, const bool DeclareCoeffPhysArrays,
288  const std::string variable, const Collections::ImplementationType ImpType)
289  : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
290  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
291  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
292  m_WaveSpace(false)
293 {
294  boost::ignore_unused(variable, ImpType);
295  int i, j, id, elmtid = 0;
296  set<int> tracesDone;
297 
304 
310 
311  // First loop over boundary conditions to reorder
312  // Dirichlet boundaries
313  for (i = 0; i < bndCond.size(); ++i)
314  {
315  if (bndCond[i]->GetBoundaryConditionType() ==
317  {
318  for (j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
319  {
320  if ((exp0D =
321  std::dynamic_pointer_cast<LocalRegions::Expansion0D>(
322  bndConstraint[i]->GetExp(j))))
323  {
324  m_expType = e0D;
325 
326  PointGeom = exp0D->GetGeom()->GetVertex(0);
327  exp = MemoryManager<
328  LocalRegions::PointExp>::AllocateSharedPtr(PointGeom);
329  tracesDone.insert(PointGeom->GetVid());
330  }
331  else if ((exp1D = std::dynamic_pointer_cast<
333  bndConstraint[i]->GetExp(j))))
334  {
335  m_expType = e1D;
336 
338  exp1D->GetBasis(0)->GetBasisKey();
339  segGeom = exp1D->GetGeom1D();
340  exp =
342  bkey, segGeom);
343  tracesDone.insert(segGeom->GetGlobalID());
344  }
345  else if ((exp2D = std::dynamic_pointer_cast<
346  LocalRegions::Expansion2D>(
347  bndConstraint[i]->GetExp(j))))
348  {
349  m_expType = e2D;
350 
351  LibUtilities::BasisKey bkey0 =
352  exp2D->GetBasis(0)->GetBasisKey();
353  LibUtilities::BasisKey bkey1 =
354  exp2D->GetBasis(1)->GetBasisKey();
355  FaceGeom = exp2D->GetGeom2D();
356 
357  // if face is a quad
358  if ((QuadGeom = std::dynamic_pointer_cast<
359  SpatialDomains::QuadGeom>(FaceGeom)))
360  {
361  exp = MemoryManager<
362  LocalRegions::QuadExp>::AllocateSharedPtr(bkey0,
363  bkey1,
364  QuadGeom);
365  tracesDone.insert(QuadGeom->GetGlobalID());
366  }
367  // if face is a triangle
368  else if ((TriGeom = std::dynamic_pointer_cast<
369  SpatialDomains::TriGeom>(FaceGeom)))
370  {
371  exp = MemoryManager<
372  LocalRegions::TriExp>::AllocateSharedPtr(bkey0,
373  bkey1,
374  TriGeom);
375  tracesDone.insert(TriGeom->GetGlobalID());
376  }
377  else
378  {
380  "dynamic cast to a "
381  "proper face geometry failed");
382  }
383  }
384  // Assign next id
385  exp->SetElmtId(elmtid++);
386 
387  // Add the expansion
388  (*m_exp).push_back(exp);
389  }
390  }
391  }
392 
393  map<int, pair<SpatialDomains::Geometry1DSharedPtr, LibUtilities::BasisKey>>
394  edgeOrders;
395 
397  pair<LibUtilities::BasisKey, LibUtilities::BasisKey>>>
398  faceOrders;
399 
400  for (i = 0; i < locexp.size(); ++i)
401  {
402  if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
403  locexp[i])))
404  {
405  m_expType = e0D;
406 
407  for (j = 0; j < 2; ++j)
408  {
409  PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
410  id = PointGeom->GetVid();
411 
412  // Ignore Dirichlet edges
413  if (tracesDone.count(id) != 0)
414  {
415  continue;
416  }
417 
419  PointGeom);
420  tracesDone.insert(id);
421  exp->SetElmtId(elmtid++);
422  (*m_exp).push_back(exp);
423  }
424  }
425  else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
426  locexp[i])))
427  {
428  m_expType = e1D;
429  for (j = 0; j < locexp[i]->GetNtraces(); ++j)
430  {
431  segGeom = exp2D->GetGeom2D()->GetEdge(j);
432  id = segGeom->GetGlobalID();
433  // Ignore Dirichlet edges
434  if (tracesDone.count(id) != 0)
435  {
436  continue;
437  }
438 
439  auto it = edgeOrders.find(id);
440 
441  if (it == edgeOrders.end())
442  {
443  edgeOrders.insert(std::make_pair(
444  id, std::make_pair(segGeom,
445  locexp[i]->GetTraceBasisKey(j))));
446  }
447  else // variable modes/points
448  {
449  LibUtilities::BasisKey edge =
450  locexp[i]->GetTraceBasisKey(j);
451  LibUtilities::BasisKey existing = it->second.second;
452 
453  int np1 = edge.GetNumPoints();
454  int np2 = existing.GetNumPoints();
455  int nm1 = edge.GetNumModes();
456  int nm2 = existing.GetNumModes();
457 
458  if (np2 >= np1 && nm2 >= nm1)
459  {
460  continue;
461  }
462  else if (np2 < np1 && nm2 < nm1)
463  {
464  it->second.second = edge;
465  }
466  else
467  {
469  "inappropriate number of points/modes (max"
470  "num of points is not set with max order)");
471  }
472  }
473  }
474  }
475  else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
476  locexp[i])))
477  {
478  m_expType = e2D;
479  for (j = 0; j < exp3D->GetNtraces(); ++j)
480  {
481  FaceGeom = exp3D->GetGeom3D()->GetFace(j);
482  id = FaceGeom->GetGlobalID();
483 
484  if (tracesDone.count(id) != 0)
485  {
486  continue;
487  }
488  auto it = faceOrders.find(id);
489 
490  if (it == faceOrders.end())
491  {
492  LibUtilities::BasisKey face_dir0 =
493  locexp[i]->GetTraceBasisKey(j, 0);
494  LibUtilities::BasisKey face_dir1 =
495  locexp[i]->GetTraceBasisKey(j, 1);
496 
497  faceOrders.insert(std::make_pair(
498  id,
499  std::make_pair(FaceGeom,
500  std::make_pair(face_dir0, face_dir1))));
501  }
502  else // variable modes/points
503  {
504  LibUtilities::BasisKey face0 =
505  locexp[i]->GetTraceBasisKey(j, 0);
506  LibUtilities::BasisKey face1 =
507  locexp[i]->GetTraceBasisKey(j, 1);
508  LibUtilities::BasisKey existing0 = it->second.second.first;
509  LibUtilities::BasisKey existing1 = it->second.second.second;
510 
511  int np11 = face0.GetNumPoints();
512  int np12 = face1.GetNumPoints();
513  int np21 = existing0.GetNumPoints();
514  int np22 = existing1.GetNumPoints();
515  int nm11 = face0.GetNumModes();
516  int nm12 = face1.GetNumModes();
517  int nm21 = existing0.GetNumModes();
518  int nm22 = existing1.GetNumModes();
519 
520  if ((np22 >= np12 || np21 >= np11) &&
521  (nm22 >= nm12 || nm21 >= nm11))
522  {
523  continue;
524  }
525  else if ((np22 < np12 || np21 < np11) &&
526  (nm22 < nm12 || nm21 < nm11))
527  {
528  it->second.second.first = face0;
529  it->second.second.second = face1;
530  }
531  else
532  {
534  "inappropriate number of points/modes (max"
535  "num of points is not set with max order)");
536  }
537  }
538  }
539  }
540  }
541 
542  int nproc = m_comm->GetRowComm()->GetSize(); // number of processors
543  int tracepr = m_comm->GetRowComm()->GetRank(); // ID processor
544 
545  if (nproc > 1)
546  {
547  int tCnt = 0;
548 
549  // Count the number of traces on each partition
550  for (i = 0; i < locexp.size(); ++i)
551  {
552  tCnt += locexp[i]->GetNtraces();
553  }
554 
555  // Set up the offset and the array that will contain the list of
556  // edge IDs, then reduce this across processors.
557  Array<OneD, int> tracesCnt(nproc, 0);
558  tracesCnt[tracepr] = tCnt;
559  m_comm->GetRowComm()->AllReduce(tracesCnt, LibUtilities::ReduceSum);
560 
561  // Set up offset array.
562  int totTraceCnt = Vmath::Vsum(nproc, tracesCnt, 1);
563  Array<OneD, int> tTotOffsets(nproc, 0);
564 
565  for (i = 1; i < nproc; ++i)
566  {
567  tTotOffsets[i] = tTotOffsets[i - 1] + tracesCnt[i - 1];
568  }
569 
570  // Local list of the edges per element
571  Array<OneD, int> TracesTotID(totTraceCnt, 0);
572  Array<OneD, int> TracesTotNm0(totTraceCnt, 0);
573  Array<OneD, int> TracesTotNm1(totTraceCnt, 0);
574  Array<OneD, int> TracesTotPnts0(totTraceCnt, 0);
575  Array<OneD, int> TracesTotPnts1(totTraceCnt, 0);
576 
577  int cntr = tTotOffsets[tracepr];
578 
579  for (i = 0; i < locexp.size(); ++i)
580  {
581  if ((exp2D = locexp[i]->as<LocalRegions::Expansion2D>()))
582  {
583 
584  int nedges = locexp[i]->GetNtraces();
585 
586  for (j = 0; j < nedges; ++j, ++cntr)
587  {
588  LibUtilities::BasisKey bkeyEdge =
589  locexp[i]->GetTraceBasisKey(j);
590  TracesTotID[cntr] = exp2D->GetGeom2D()->GetEid(j);
591  TracesTotNm0[cntr] = bkeyEdge.GetNumModes();
592  TracesTotPnts0[cntr] = bkeyEdge.GetNumPoints();
593  }
594  }
595  else if ((exp3D = locexp[i]->as<LocalRegions::Expansion3D>()))
596  {
597  int nfaces = locexp[i]->GetNtraces();
598 
599  for (j = 0; j < nfaces; ++j, ++cntr)
600  {
601  LibUtilities::BasisKey face_dir0 =
602  locexp[i]->GetTraceBasisKey(j, 0);
603  LibUtilities::BasisKey face_dir1 =
604  locexp[i]->GetTraceBasisKey(j, 1);
605 
606  TracesTotID[cntr] = exp3D->GetGeom3D()->GetFid(j);
607  TracesTotNm0[cntr] = face_dir0.GetNumModes();
608  TracesTotNm1[cntr] = face_dir1.GetNumModes();
609  TracesTotPnts0[cntr] = face_dir0.GetNumPoints();
610  TracesTotPnts1[cntr] = face_dir1.GetNumPoints();
611  }
612  }
613  }
614 
615  m_comm->GetRowComm()->AllReduce(TracesTotID, LibUtilities::ReduceSum);
616  m_comm->GetRowComm()->AllReduce(TracesTotNm0, LibUtilities::ReduceSum);
617  m_comm->GetRowComm()->AllReduce(TracesTotPnts0,
619  if (m_expType == e2D)
620  {
621  m_comm->AllReduce(TracesTotNm1, LibUtilities::ReduceSum);
622  m_comm->AllReduce(TracesTotPnts1, LibUtilities::ReduceSum);
623  }
624 
625  if (edgeOrders.size())
626  {
627  for (i = 0; i < totTraceCnt; ++i)
628  {
629  auto it = edgeOrders.find(TracesTotID[i]);
630 
631  if (it == edgeOrders.end())
632  {
633  continue;
634  }
635 
636  LibUtilities::BasisKey existing = it->second.second;
637  LibUtilities::BasisKey edge(
638  existing.GetBasisType(), TracesTotNm0[i],
639  LibUtilities::PointsKey(TracesTotPnts0[i],
640  existing.GetPointsType()));
641 
642  int np1 = edge.GetNumPoints();
643  int np2 = existing.GetNumPoints();
644  int nm1 = edge.GetNumModes();
645  int nm2 = existing.GetNumModes();
646 
647  if (np2 >= np1 && nm2 >= nm1)
648  {
649  continue;
650  }
651  else if (np2 < np1 && nm2 < nm1)
652  {
653  it->second.second = edge;
654  }
655  else
656  {
658  "inappropriate number of points/modes (max "
659  "num of points is not set with max order)");
660  }
661  }
662  }
663  else if (faceOrders.size())
664  {
665  for (i = 0; i < totTraceCnt; ++i)
666  {
667  auto it = faceOrders.find(TracesTotID[i]);
668 
669  if (it == faceOrders.end())
670  {
671  continue;
672  }
673 
674  LibUtilities::BasisKey existing0 = it->second.second.first;
675  LibUtilities::BasisKey existing1 = it->second.second.second;
676  LibUtilities::BasisKey face0(
677  existing0.GetBasisType(), TracesTotNm0[i],
678  LibUtilities::PointsKey(TracesTotPnts0[i],
679  existing0.GetPointsType()));
680  LibUtilities::BasisKey face1(
681  existing1.GetBasisType(), TracesTotNm1[i],
682  LibUtilities::PointsKey(TracesTotPnts1[i],
683  existing1.GetPointsType()));
684 
685  int np11 = face0.GetNumPoints();
686  int np12 = face1.GetNumPoints();
687  int np21 = existing0.GetNumPoints();
688  int np22 = existing1.GetNumPoints();
689  int nm11 = face0.GetNumModes();
690  int nm12 = face1.GetNumModes();
691  int nm21 = existing0.GetNumModes();
692  int nm22 = existing1.GetNumModes();
693 
694  if ((np22 >= np12 || np21 >= np11) &&
695  (nm22 >= nm12 || nm21 >= nm11))
696  {
697  continue;
698  }
699  else if ((np22 < np12 || np21 < np11) &&
700  (nm22 < nm12 || nm21 < nm11))
701  {
702  it->second.second.first = face0;
703  it->second.second.second = face1;
704  }
705  else
706  {
708  "inappropriate number of points/modes (max "
709  "num of points is not set with max order)");
710  }
711  }
712  }
713  }
714 
715  if (edgeOrders.size())
716  {
717  for (auto &it : edgeOrders)
718  {
720  it.second.second, it.second.first);
721  exp->SetElmtId(elmtid++);
722  (*m_exp).push_back(exp);
723  }
724  }
725  else
726  {
727  for (auto &it : faceOrders)
728  {
729  FaceGeom = it.second.first;
730 
731  if ((QuadGeom = std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
732  FaceGeom)))
733  {
735  it.second.second.first, it.second.second.second, QuadGeom);
736  }
737  else if ((TriGeom =
738  std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
739  FaceGeom)))
740  {
742  it.second.second.first, it.second.second.second, TriGeom);
743  }
744  exp->SetElmtId(elmtid++);
745  (*m_exp).push_back(exp);
746  }
747  }
748 
749  // Set up m_coeffs, m_phys and offset arrays.
750  SetupCoeffPhys(DeclareCoeffPhysArrays);
751 
752  // Set up collections
753  if (m_expType != e0D)
754  {
755  CreateCollections(ImpType);
756  }
757 }
758 
759 /**
760  * Set expansions for localtrace space expansions used in
761  * DisContField as part of Gradient Jump Penalisation
762  *
763  * @param pSession A session within information about expansion
764  * @param locexp Complete domain expansion list.
765  * @param graph mesh corresponding to the expansion list.
766  * @param DeclareCoeffPhysArrays Declare the coefficient and
767  * phys space arrays
768  * @param variable The variable name associated with the expansion
769  * @param ImpType Detail about the implementation type to use
770  * in operators
771  */
773  const LocalRegions::ExpansionVector &locexp,
775  const bool DeclareCoeffPhysArrays, const std::string variable,
776  const Collections::ImplementationType ImpType)
777  : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
778  m_physState(false),
779  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
780  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
781  m_WaveSpace(false)
782 {
783  boost::ignore_unused(variable, ImpType);
784  int i, j, elmtid = 0;
785 
792 
798 
799  for (i = 0; i < locexp.size(); ++i)
800  {
801  if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
802  locexp[i])))
803  {
804  m_expType = e0D;
805 
806  for (j = 0; j < 2; ++j)
807  {
808  PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
809 
811  PointGeom);
812  exp->SetElmtId(elmtid++);
813  (*m_exp).push_back(exp);
814  }
815  }
816  else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
817  locexp[i])))
818  {
819  m_expType = e1D;
820  LibUtilities::BasisKey edgeKey0 =
821  locexp[i]->GetBasis(0)->GetBasisKey();
822 
823  for (j = 0; j < locexp[i]->GetNtraces(); ++j)
824  {
825  segGeom = exp2D->GetGeom2D()->GetEdge(j);
826 
827  int dir = exp2D->GetGeom2D()->GetDir(j);
828 
829  if (locexp[i]->GetNtraces() == 3)
830  {
831  LibUtilities::BasisKey edgeKey =
832  locexp[i]->GetBasis(dir)->GetBasisKey();
833 
834  LibUtilities::BasisKey nEdgeKey(edgeKey0.GetBasisType(),
835  edgeKey.GetNumModes(),
836  edgeKey.GetPointsKey());
837 
838  exp =
840  nEdgeKey, segGeom);
841  }
842  else
843  {
844  exp =
846  locexp[i]->GetBasis(dir)->GetBasisKey(), segGeom);
847  }
848 
849  exp->SetElmtId(elmtid++);
850  (*m_exp).push_back(exp);
851  }
852  }
853  else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
854  locexp[i])))
855  {
856  m_expType = e2D;
857 
858  LibUtilities::BasisKey face0_dir0 =
859  locexp[i]->GetBasis(0)->GetBasisKey();
860  LibUtilities::BasisKey face0_dir1 =
861  locexp[i]->GetBasis(1)->GetBasisKey();
862 
863  for (j = 0; j < exp3D->GetNtraces(); ++j)
864  {
865  FaceGeom = exp3D->GetGeom3D()->GetFace(j);
866 
867  int dir0 = exp3D->GetGeom3D()->GetDir(j, 0);
868  int dir1 = exp3D->GetGeom3D()->GetDir(j, 1);
869 
870  LibUtilities::BasisKey face_dir0 =
871  locexp[i]->GetBasis(dir0)->GetBasisKey();
872  LibUtilities::BasisKey face_dir1 =
873  locexp[i]->GetBasis(dir1)->GetBasisKey();
874 
875  if ((QuadGeom =
876  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
877  FaceGeom)))
878  {
879  exp =
881  face_dir0, face_dir1, QuadGeom);
882  }
883  else if ((TriGeom = std::dynamic_pointer_cast<
884  SpatialDomains::TriGeom>(FaceGeom)))
885  {
886 
887  LibUtilities::BasisKey nface_dir0(face0_dir0.GetBasisType(),
888  face_dir0.GetNumModes(),
889  face_dir0.GetPointsKey());
890  LibUtilities::BasisKey nface_dir1(face0_dir1.GetBasisType(),
891  face_dir1.GetNumModes(),
892  face_dir1.GetPointsKey());
893  exp =
895  nface_dir0, nface_dir1, TriGeom);
896  }
897  exp->SetElmtId(elmtid++);
898  (*m_exp).push_back(exp);
899  }
900  }
901  }
902 
903  // Set up m_coeffs, m_phys and offset arrays.
904  SetupCoeffPhys(DeclareCoeffPhysArrays);
905 
906  // Set up collections
907  if (m_expType != e0D)
908  {
909  CreateCollections(ImpType);
910  }
911 }
912 
913 /**
914  * Fills the list of local expansions with the trace from the mesh
915  * specified by \a domain. This CompositeMap contains a list of
916  * Composites which define the boundary. It is also used to set up
917  * expansion domains in the 1D Pulse Wave solver.
918  *
919  * @param pSession A session within information about expansion
920  * @param domain A domain, comprising of one or more composite
921  * regions,
922  * @param graph A mesh, containing information about the
923  * domain and the spectral/hp element expansion.
924  * @param DeclareCoeffPhysArrays Declare the coefficient and
925  * phys space arrays. Default is true.
926  * @param variable The variable name associated with the expansion
927  * @param SetToOneSpaceDimension Reduce to one space dimension expansion
928  * @param comm An optional communicator that can be used with the
929  * boundary expansion in case of more global
930  * parallel operations. Default to a Null Communicator
931  * @param ImpType Detail about the implementation type to use
932  * in operators. Default is eNoImpType.
933  *
934  */
936  const SpatialDomains::CompositeMap &domain,
938  const bool DeclareCoeffPhysArrays, const std::string variable,
939  bool SetToOneSpaceDimension,
940  const LibUtilities::CommSharedPtr comm,
941  const Collections::ImplementationType ImpType)
942  : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
943  m_exp(MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr()),
944  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
945  m_WaveSpace(false)
946 {
947  int j, elmtid = 0;
952 
954 
956 
957  int meshdim = graph->GetMeshDimension();
958 
959  // Retrieve the list of expansions (needed of meshdim == 1
960  const SpatialDomains::ExpansionInfoMap &expansions =
961  graph->GetExpansionInfo(variable);
962 
963  // Retrieve the list of expansions
964  // Process each composite region.
965  for (auto &compIt : domain)
966  {
967  // Process each expansion in the region.
968  for (j = 0; j < compIt.second->m_geomVec.size(); ++j)
969  {
970  if ((PtGeom = std::dynamic_pointer_cast<SpatialDomains::PointGeom>(
971  compIt.second->m_geomVec[j])))
972  {
973  m_expType = e0D;
974 
976  PtGeom);
977  }
978  else if ((SegGeom =
979  std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
980  compIt.second->m_geomVec[j])))
981  {
982  m_expType = e1D;
983 
984  // Retrieve the basis key from the expansion.
986 
987  if (meshdim == 1)
988  {
989  auto expIt = expansions.find(SegGeom->GetGlobalID());
990  ASSERTL0(expIt != expansions.end(),
991  "Failed to find basis key");
992  bkey = expIt->second->m_basisKeyVector[0];
993  }
994  else
995  {
996  bkey = graph->GetEdgeBasisKey(SegGeom, variable);
997  }
998 
999  if (SetToOneSpaceDimension)
1000  {
1001  SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
1002  SegGeom->GenerateOneSpaceDimGeom();
1003 
1004  exp =
1006  bkey, OneDSegmentGeom);
1007  }
1008  else
1009  {
1010 
1011  exp =
1013  bkey, SegGeom);
1014  }
1015  }
1016  else if ((TriGeom =
1017  std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
1018  compIt.second->m_geomVec[j])))
1019  {
1020  m_expType = e2D;
1021 
1022  LibUtilities::BasisKey TriBa =
1023  graph->GetFaceBasisKey(TriGeom, 0, variable);
1024  LibUtilities::BasisKey TriBb =
1025  graph->GetFaceBasisKey(TriGeom, 1, variable);
1026 
1027  if (graph->GetExpansionInfo()
1028  .begin()
1029  ->second->m_basisKeyVector[0]
1030  .GetBasisType() == LibUtilities::eGLL_Lagrange)
1031  {
1032  NEKERROR(ErrorUtil::efatal, "This method needs sorting");
1034 
1035  exp = MemoryManager<
1036  LocalRegions::NodalTriExp>::AllocateSharedPtr(TriBa,
1037  TriBb,
1038  TriNb,
1039  TriGeom);
1040  }
1041  else
1042  {
1043  exp =
1045  TriBa, TriBb, TriGeom);
1046  }
1047  }
1048  else if ((QuadGeom =
1049  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1050  compIt.second->m_geomVec[j])))
1051  {
1052  m_expType = e2D;
1053 
1054  LibUtilities::BasisKey QuadBa =
1055  graph->GetFaceBasisKey(QuadGeom, 0, variable);
1056  LibUtilities::BasisKey QuadBb =
1057  graph->GetFaceBasisKey(QuadGeom, 1, variable);
1058 
1060  QuadBa, QuadBb, QuadGeom);
1061  }
1062  else
1063  {
1065  "dynamic cast to a Geom (possibly 3D) failed");
1066  }
1067 
1068  exp->SetElmtId(elmtid++);
1069  (*m_exp).push_back(exp);
1070  }
1071  }
1072 
1073  // Set up m_coeffs, m_phys and offset arrays.
1074  SetupCoeffPhys(DeclareCoeffPhysArrays);
1075 
1076  if (m_expType != e0D)
1077  {
1078  CreateCollections(ImpType);
1079  }
1080 }
1081 
1082 /**
1083  * Each expansion (local element) is processed in turn to
1084  * determine the number of coefficients and physical data
1085  * points it contributes to the domain. Twoe arrays,
1086  * #m_coeff_offset are #m_phys_offset are also initialised and
1087  * updated to store the data offsets of each element in the
1088  * #m_coeffs and #m_phys arrays, and the element id that each
1089  * consecutive block is associated respectively.
1090  * Finally we initialise #m_coeffs and #m_phys
1091  */
1092 void ExpList::SetupCoeffPhys(bool DeclareCoeffPhysArrays, bool SetupOffsets)
1093 {
1094  if (SetupOffsets)
1095  {
1096  int i;
1097 
1098  // Set up offset information and array sizes
1101 
1102  m_ncoeffs = m_npoints = 0;
1103 
1104  for (i = 0; i < m_exp->size(); ++i)
1105  {
1107  m_phys_offset[i] = m_npoints;
1108  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
1109  m_npoints += (*m_exp)[i]->GetTotPoints();
1110  }
1111  }
1112 
1113  if (DeclareCoeffPhysArrays)
1114  {
1117  }
1118 
1120 
1121  for (int i = 0; i < m_exp->size(); ++i)
1122  {
1123  int coeffs_offset = m_coeff_offset[i];
1124 
1125  int loccoeffs = (*m_exp)[i]->GetNcoeffs();
1126 
1127  for (int j = 0; j < loccoeffs; ++j)
1128  {
1129  m_coeffsToElmt[coeffs_offset + j].first = i;
1130  m_coeffsToElmt[coeffs_offset + j].second = j;
1131  }
1132  }
1133 }
1134 
1136  const SpatialDomains::ExpansionInfoMap &expmap)
1137 {
1138 
1140  SpatialDomains::TriGeomSharedPtr TriangleGeom;
1141  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
1146 
1147  int id = 0;
1149 
1150  m_expType = eNoType;
1151  // Process each expansion in the graph
1152  for (auto &expIt : expmap)
1153  {
1154  const SpatialDomains::ExpansionInfoShPtr expInfo = expIt.second;
1155 
1156  switch (expInfo->m_basisKeyVector.size())
1157  {
1158  case 1: // Segment Expansions
1159  {
1161  "Cannot mix expansion dimensions in one vector");
1162  m_expType = e1D;
1163 
1164  if ((SegmentGeom =
1165  std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1166  expInfo->m_geomShPtr)))
1167  {
1168  // Retrieve basis key from expansion
1169  LibUtilities::BasisKey bkey = expInfo->m_basisKeyVector[0];
1170 
1171  exp =
1173  bkey, SegmentGeom);
1174  }
1175  else
1176  {
1178  "dynamic cast to a 1D Geom failed");
1179  }
1180  }
1181  break;
1182  case 2:
1183  {
1185  "Cannot mix expansion dimensions in one vector");
1186  m_expType = e2D;
1187 
1188  LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1189  LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1190 
1191  if ((TriangleGeom =
1192  std::dynamic_pointer_cast<SpatialDomains ::TriGeom>(
1193  expInfo->m_geomShPtr)))
1194  {
1195 
1196  // This is not elegantly implemented needs re-thinking.
1198  {
1200  Ba.GetNumModes(),
1201  Ba.GetPointsKey());
1202 
1203  LibUtilities::PointsType TriNb =
1206  AllocateSharedPtr(newBa, Bb, TriNb, TriangleGeom);
1207  }
1208  else
1209  {
1211  AllocateSharedPtr(Ba, Bb, TriangleGeom);
1212  }
1213  }
1214  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
1215  SpatialDomains::QuadGeom>(expInfo->m_geomShPtr)))
1216  {
1217  exp =
1219  Ba, Bb, QuadrilateralGeom);
1220  }
1221  else
1222  {
1224  "dynamic cast to a 2D Geom failed");
1225  }
1226  }
1227  break;
1228  case 3:
1229  {
1231  "Cannot mix expansion dimensions in one vector");
1232  m_expType = e3D;
1233 
1234  LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1235  LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1236  LibUtilities::BasisKey Bc = expInfo->m_basisKeyVector[2];
1237 
1238  if ((TetGeom =
1239  std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
1240  expInfo->m_geomShPtr)))
1241  {
1242 
1245  {
1247  "LocalRegions::NodalTetExp is not implemented "
1248  "yet");
1249  }
1250  else
1251  {
1252  exp = MemoryManager<
1253  LocalRegions::TetExp>::AllocateSharedPtr(Ba, Bb, Bc,
1254  TetGeom);
1255  }
1256  }
1257  else if ((PrismGeom = std::dynamic_pointer_cast<
1258  SpatialDomains ::PrismGeom>(
1259  expInfo->m_geomShPtr)))
1260  {
1261  exp = MemoryManager<
1262  LocalRegions::PrismExp>::AllocateSharedPtr(Ba, Bb, Bc,
1263  PrismGeom);
1264  }
1265  else if ((PyrGeom = std::dynamic_pointer_cast<
1266  SpatialDomains::PyrGeom>(expInfo->m_geomShPtr)))
1267  {
1268 
1269  exp =
1271  Ba, Bb, Bc, PyrGeom);
1272  }
1273  else if ((HexGeom = std::dynamic_pointer_cast<
1274  SpatialDomains::HexGeom>(expInfo->m_geomShPtr)))
1275  {
1276  exp =
1278  Ba, Bb, Bc, HexGeom);
1279  }
1280  else
1281  {
1283  "dynamic cast to a Geom failed");
1284  }
1285  }
1286  break;
1287  default:
1289  "Dimension of basis key is greater than 3");
1290  }
1291 
1292  // Assign next id
1293  exp->SetElmtId(id++);
1294 
1295  // Add the expansion
1296  (*m_exp).push_back(exp);
1297  }
1298 }
1299 
1300 /**
1301  *
1302  */
1304 {
1305  return m_expType;
1306 }
1307 
1309 {
1310 }
1311 
1312 /**
1313  * Retrieves the block matrix specified by \a bkey, and computes
1314  * \f$ y=Mx \f$.
1315  * @param gkey GlobalMatrixKey specifying the block matrix to
1316  * use in the matrix-vector multiply.
1317  * @param inarray Input vector \f$ x \f$.
1318  * @param outarray Output vector \f$ y \f$.
1319  */
1321  const Array<OneD, const NekDouble> &inarray,
1322  Array<OneD, NekDouble> &outarray)
1323 {
1324  // Retrieve the block matrix using the given key.
1325  const DNekScalBlkMatSharedPtr &blockmat = GetBlockMatrix(gkey);
1326  int nrows = blockmat->GetRows();
1327  int ncols = blockmat->GetColumns();
1328 
1329  // Create NekVectors from the given data arrays
1330  NekVector<NekDouble> in(ncols, inarray, eWrapper);
1331  NekVector<NekDouble> out(nrows, outarray, eWrapper);
1332 
1333  // Perform matrix-vector multiply.
1334  out = (*blockmat) * in;
1335 }
1336 
1337 /**
1338  * multiply the metric jacobi and quadrature weights
1339  */
1341  const Array<OneD, const NekDouble> &inarray,
1342  Array<OneD, NekDouble> &outarray)
1343 {
1344  Array<OneD, NekDouble> e_outarray;
1345 
1346  for (int i = 0; i < (*m_exp).size(); ++i)
1347  {
1348  (*m_exp)[i]->MultiplyByQuadratureMetric(inarray + m_phys_offset[i],
1349  e_outarray = outarray +
1350  m_phys_offset[i]);
1351  }
1352 }
1353 
1354 /**
1355  * Divided by the metric jacobi and quadrature weights
1356  */
1358  const Array<OneD, const NekDouble> &inarray,
1359  Array<OneD, NekDouble> &outarray)
1360 {
1361  Array<OneD, NekDouble> e_outarray;
1362 
1363  for (int i = 0; i < (*m_exp).size(); ++i)
1364  {
1365  (*m_exp)[i]->DivideByQuadratureMetric(inarray + m_phys_offset[i],
1366  e_outarray =
1367  outarray + m_phys_offset[i]);
1368  }
1369 }
1370 
1371 /**
1372  * The operation is evaluated locally for every element by the function
1373  * StdRegions#Expansion#IProductWRTDerivBase.
1374  *
1375  * @param dir {0,1} is the direction in which the
1376  * derivative of the basis should be taken
1377  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1378  * containing the values of the function
1379  * \f$f(\boldsymbol{x})\f$ at the quadrature
1380  * points \f$\boldsymbol{x}_i\f$.
1381  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1382  * used to store the result.
1383  */
1385  const Array<OneD, const NekDouble> &inarray,
1386  Array<OneD, NekDouble> &outarray)
1387 {
1388  int i;
1389 
1390  Array<OneD, NekDouble> e_outarray;
1391 
1392  for (i = 0; i < (*m_exp).size(); ++i)
1393  {
1394  (*m_exp)[i]->IProductWRTDerivBase(dir, inarray + m_phys_offset[i],
1395  e_outarray =
1396  outarray + m_coeff_offset[i]);
1397  }
1398 }
1399 
1400 /**
1401  * @brief Directional derivative along a given direction
1402  *
1403  */
1405  const Array<OneD, const NekDouble> &direction,
1406  const Array<OneD, const NekDouble> &inarray,
1407  Array<OneD, NekDouble> &outarray)
1408 {
1409  int npts_e;
1410  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1411  int nq = direction.size() / coordim;
1412 
1413  Array<OneD, NekDouble> e_outarray;
1414  Array<OneD, NekDouble> e_MFdiv;
1415 
1416  Array<OneD, NekDouble> locdir;
1417 
1418  for (int i = 0; i < (*m_exp).size(); ++i)
1419  {
1420  npts_e = (*m_exp)[i]->GetTotPoints();
1421  locdir = Array<OneD, NekDouble>(npts_e * coordim);
1422 
1423  for (int k = 0; k < coordim; ++k)
1424  {
1425  Vmath::Vcopy(npts_e, &direction[k * nq + m_phys_offset[i]], 1,
1426  &locdir[k * npts_e], 1);
1427  }
1428 
1429  (*m_exp)[i]->IProductWRTDirectionalDerivBase(
1430  locdir, inarray + m_phys_offset[i],
1431  e_outarray = outarray + m_coeff_offset[i]);
1432  }
1433 }
1434 
1435 /**
1436  * The operation is evaluated locally for every element by the function
1437  * StdRegions#StdExpansion#IProductWRTDerivBase.
1438  *
1439  * @param inarray An array of arrays of size \f$Q_{\mathrm{tot}}\f$
1440  * containing the values of the function
1441  * \f$f(\boldsymbol{x})\f$ at the quadrature
1442  * points \f$\boldsymbol{x}_i\f$ in dir directions.
1443  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1444  * used to store the result.
1445  */
1447  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1448  Array<OneD, NekDouble> &outarray)
1449 {
1450  Array<OneD, NekDouble> tmp0, tmp1, tmp2;
1451  // assume coord dimension defines the size of Deriv Base
1452  int dim = GetCoordim(0);
1453 
1454  ASSERTL1(inarray.size() >= dim, "inarray is not of sufficient dimension");
1455 
1456  // initialise if required
1458  {
1459  for (int i = 0; i < m_collections.size(); ++i)
1460  {
1462  }
1464  }
1465 
1466  LibUtilities::Timer timer;
1467 
1468  LIKWID_MARKER_START("IProductWRTDerivBase_coll");
1469  timer.Start();
1470 
1471  switch (dim)
1472  {
1473  case 1:
1474  for (int i = 0; i < m_collections.size(); ++i)
1475  {
1476  m_collections[i].ApplyOperator(
1478  inarray[0] + m_coll_phys_offset[i],
1479  tmp0 = outarray + m_coll_coeff_offset[i]);
1480  }
1481  break;
1482  case 2:
1483  for (int i = 0; i < m_collections.size(); ++i)
1484  {
1485  m_collections[i].ApplyOperator(
1487  inarray[0] + m_coll_phys_offset[i],
1488  tmp0 = inarray[1] + m_coll_phys_offset[i],
1489  tmp1 = outarray + m_coll_coeff_offset[i]);
1490  }
1491  break;
1492  case 3:
1493  for (int i = 0; i < m_collections.size(); ++i)
1494  {
1495  m_collections[i].ApplyOperator(
1497  inarray[0] + m_coll_phys_offset[i],
1498  tmp0 = inarray[1] + m_coll_phys_offset[i],
1499  tmp1 = inarray[2] + m_coll_phys_offset[i],
1500  tmp2 = outarray + m_coll_coeff_offset[i]);
1501  }
1502  break;
1503  default:
1504  NEKERROR(ErrorUtil::efatal, "Dimension of inarray not correct");
1505  break;
1506  }
1507 
1508  timer.Stop();
1509  LIKWID_MARKER_STOP("IProductWRTDerivBase_coll");
1510 
1511  // Elapsed time
1512  timer.AccumulateRegion("Collections:IProductWRTDerivBase", 10);
1513 }
1514 /**
1515  * Given a function \f$f(\boldsymbol{x})\f$ evaluated at
1516  * the quadrature points, this function calculates the
1517  * derivatives \f$\frac{d}{dx_1}\f$, \f$\frac{d}{dx_2}\f$
1518  * and \f$\frac{d}{dx_3}\f$ of the function
1519  * \f$f(\boldsymbol{x})\f$ at the same quadrature
1520  * points. The local distribution of the quadrature points
1521  * allows an elemental evaluation of the derivative. This
1522  * is done by a call to the function
1523  * StdRegions#StdExpansion#PhysDeriv.
1524  *
1525  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1526  * containing the values of the function
1527  * \f$f(\boldsymbol{x})\f$ at the quadrature
1528  * points \f$\boldsymbol{x}_i\f$.
1529  * @param out_d0 The discrete evaluation of the
1530  * derivative\f$\frac{d}{dx_1}\f$ will
1531  * be stored in this array of size
1532  * \f$Q_{\mathrm{tot}}\f$.
1533  * @param out_d1 The discrete evaluation of the
1534  * derivative\f$\frac{d}{dx_2}\f$ will be
1535  * stored in this array of size
1536  * \f$Q_{\mathrm{tot}}\f$. Note that if no
1537  * memory is allocated for \a out_d1, the
1538  * derivative \f$\frac{d}{dx_2}\f$ will not be
1539  * calculated.
1540  * @param out_d2 The discrete evaluation of the
1541  * derivative\f$\frac{d}{dx_3}\f$ will be
1542  * stored in this array of size
1543  * \f$Q_{\mathrm{tot}}\f$. Note that if no
1544  * memory is allocated for \a out_d2, the
1545  * derivative \f$\frac{d}{dx_3}\f$ will not be
1546  * calculated.
1547  */
1549  Array<OneD, NekDouble> &out_d0,
1550  Array<OneD, NekDouble> &out_d1,
1551  Array<OneD, NekDouble> &out_d2)
1552 {
1553  Array<OneD, NekDouble> e_out_d0;
1554  Array<OneD, NekDouble> e_out_d1;
1555  Array<OneD, NekDouble> e_out_d2;
1556  int offset;
1557 
1558  // initialise if required
1560  {
1561  for (int i = 0; i < m_collections.size(); ++i)
1562  {
1563  m_collections[i].Initialise(Collections::ePhysDeriv);
1564  }
1566  }
1567 
1568  LibUtilities::Timer timer;
1569  timer.Start();
1570  for (int i = 0; i < m_collections.size(); ++i)
1571  {
1572  offset = m_coll_phys_offset[i];
1573  e_out_d0 = out_d0 + offset;
1574  e_out_d1 = out_d1 + offset;
1575  e_out_d2 = out_d2 + offset;
1576 
1577  m_collections[i].ApplyOperator(Collections::ePhysDeriv,
1578  inarray + offset, e_out_d0, e_out_d1,
1579  e_out_d2);
1580  }
1581  timer.Stop();
1582  // Elapsed time
1583  timer.AccumulateRegion("Collections:PhysDeriv", 10);
1584 }
1585 
1586 void ExpList::v_PhysDeriv(const int dir,
1587  const Array<OneD, const NekDouble> &inarray,
1588  Array<OneD, NekDouble> &out_d)
1589 {
1590  Direction edir = DirCartesianMap[dir];
1591  v_PhysDeriv(edir, inarray, out_d);
1592 }
1593 
1595  const Array<OneD, const NekDouble> &inarray,
1596  Array<OneD, NekDouble> &out_d)
1597 {
1598  int i;
1599  if (edir == MultiRegions::eS)
1600  {
1601  Array<OneD, NekDouble> e_out_ds;
1602  for (i = 0; i < (*m_exp).size(); ++i)
1603  {
1604  e_out_ds = out_d + m_phys_offset[i];
1605  (*m_exp)[i]->PhysDeriv_s(inarray + m_phys_offset[i], e_out_ds);
1606  }
1607  }
1608  else if (edir == MultiRegions::eN)
1609  {
1610  Array<OneD, NekDouble> e_out_dn;
1611  for (i = 0; i < (*m_exp).size(); i++)
1612  {
1613  e_out_dn = out_d + m_phys_offset[i];
1614  (*m_exp)[i]->PhysDeriv_n(inarray + m_phys_offset[i], e_out_dn);
1615  }
1616  }
1617  else
1618  {
1619 
1620  // initialise if required
1622  {
1623  for (int i = 0; i < m_collections.size(); ++i)
1624  {
1625  m_collections[i].Initialise(Collections::ePhysDeriv);
1626  }
1628  }
1629 
1630  // convert enum into int
1631  int intdir = (int)edir;
1632  Array<OneD, NekDouble> e_out_d;
1633  int offset;
1634  for (int i = 0; i < m_collections.size(); ++i)
1635  {
1636  offset = m_coll_phys_offset[i];
1637  e_out_d = out_d + offset;
1638 
1639  m_collections[i].ApplyOperator(Collections::ePhysDeriv, intdir,
1640  inarray + offset, e_out_d);
1641  }
1642  }
1643 }
1644 
1647 {
1648  int nq = GetTotPoints();
1649  Array<OneD, NekDouble> Vx(nq);
1650  Array<OneD, NekDouble> Uy(nq);
1652 
1653  bool halfMode = false;
1654  if (GetExpType() == e3DH1D)
1655  {
1656  m_session->MatchSolverInfo("ModeType", "HalfMode", halfMode, false);
1657  }
1658 
1659  switch (m_expType)
1660  {
1661  case e2D:
1662  {
1663  PhysDeriv(xDir, Vel[yDir], Vx);
1664  PhysDeriv(yDir, Vel[xDir], Uy);
1665 
1666  Vmath::Vsub(nq, Vx, 1, Uy, 1, Dummy, 1);
1667 
1668  PhysDeriv(Dummy, Q[1], Q[0]);
1669 
1670  Vmath::Smul(nq, -1.0, Q[1], 1, Q[1], 1);
1671  }
1672  break;
1673 
1674  case e3D:
1675  case e3DH1D:
1676  case e3DH2D:
1677  {
1678  Array<OneD, NekDouble> Vz(nq);
1679  Array<OneD, NekDouble> Uz(nq);
1680  Array<OneD, NekDouble> Wx(nq);
1681  Array<OneD, NekDouble> Wy(nq);
1682 
1683  PhysDeriv(Vel[xDir], Dummy, Uy, Uz);
1684  PhysDeriv(Vel[yDir], Vx, Dummy, Vz);
1685  PhysDeriv(Vel[zDir], Wx, Wy, Dummy);
1686 
1687  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
1688  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
1689  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
1690 
1691  PhysDeriv(Q[0], Dummy, Uy, Uz);
1692  PhysDeriv(Q[1], Vx, Dummy, Vz);
1693  PhysDeriv(Q[2], Wx, Wy, Dummy);
1694 
1695  // For halfmode, need to change the sign of z derivatives
1696  if (halfMode)
1697  {
1698  Vmath::Neg(nq, Uz, 1);
1699  Vmath::Neg(nq, Vz, 1);
1700  }
1701 
1702  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
1703  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
1704  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
1705  }
1706  break;
1707  default:
1708  ASSERTL0(0, "Dimension not supported");
1709  break;
1710  }
1711 }
1712 
1714  const Array<OneD, const NekDouble> &direction,
1715  const Array<OneD, const NekDouble> &inarray,
1716  Array<OneD, NekDouble> &outarray)
1717 {
1718  int npts_e;
1719  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1720  int nq = direction.size() / coordim;
1721 
1722  Array<OneD, NekDouble> e_outarray;
1723  Array<OneD, NekDouble> e_MFdiv;
1724  Array<OneD, NekDouble> locdir;
1725 
1726  for (int i = 0; i < (*m_exp).size(); ++i)
1727  {
1728  npts_e = (*m_exp)[i]->GetTotPoints();
1729  locdir = Array<OneD, NekDouble>(npts_e * coordim);
1730 
1731  for (int k = 0; k < coordim; ++k)
1732  {
1733  Vmath::Vcopy(npts_e, &direction[k * nq + m_phys_offset[i]], 1,
1734  &locdir[k * npts_e], 1);
1735  }
1736 
1737  (*m_exp)[i]->PhysDirectionalDeriv(inarray + m_phys_offset[i], locdir,
1738  e_outarray =
1739  outarray + m_phys_offset[i]);
1740  }
1741 }
1742 
1744  const NekDouble alpha, const NekDouble exponent,
1745  const NekDouble cutoff)
1746 {
1747  Array<OneD, NekDouble> e_array;
1748 
1749  for (int i = 0; i < (*m_exp).size(); ++i)
1750  {
1751  (*m_exp)[i]->ExponentialFilter(e_array = array + m_phys_offset[i],
1752  alpha, exponent, cutoff);
1753  }
1754 }
1755 
1756 /**
1757  * The coefficients of the function to be acted upon
1758  * should be contained in the \param inarray. The
1759  * resulting coefficients are stored in \param outarray
1760  *
1761  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
1762  * containing the inner product.
1763  */
1765  Array<OneD, NekDouble> &outarray)
1766 {
1768  const DNekScalBlkMatSharedPtr &InvMass = GetBlockMatrix(mkey);
1769 
1770  // Inverse mass matrix
1771  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
1772  if (inarray.get() == outarray.get())
1773  {
1774  NekVector<NekDouble> in(m_ncoeffs, inarray); // copy data
1775  out = (*InvMass) * in;
1776  }
1777  else
1778  {
1779  NekVector<NekDouble> in(m_ncoeffs, inarray, eWrapper);
1780  out = (*InvMass) * in;
1781  }
1782 }
1783 
1784 /**
1785  * Given a function \f$u(\boldsymbol{x})\f$ defined at the
1786  * quadrature points, this function determines the
1787  * transformed elemental coefficients \f$\hat{u}_n^e\f$
1788  * employing a discrete elemental Galerkin projection from
1789  * physical space to coefficient space. For each element,
1790  * the operation is evaluated locally by the function
1791  * StdRegions#StdExpansion#IproductWRTBase followed by a
1792  * call to #MultiRegions#MultiplyByElmtInvMass.
1793  *
1794  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1795  * containing the values of the function
1796  * \f$f(\boldsymbol{x})\f$ at the quadrature
1797  * points \f$\boldsymbol{x}_i\f$.
1798  * @param outarray The resulting coefficients
1799  * \f$\hat{u}_n^e\f$ will be stored in this
1800  * array of size \f$N_{\mathrm{eof}}\f$.
1801  */
1803  Array<OneD, NekDouble> &outarray)
1804 {
1806 
1807  IProductWRTBase(inarray, f);
1808  MultiplyByElmtInvMass(f, outarray);
1809 }
1810 
1812  const Array<OneD, const NekDouble> &inarray,
1813  Array<OneD, NekDouble> &outarray)
1814 {
1815  int i;
1816 
1817  Array<OneD, NekDouble> e_outarray;
1818 
1819  for (i = 0; i < (*m_exp).size(); ++i)
1820  {
1821  (*m_exp)[i]->FwdTransBndConstrained(inarray + m_phys_offset[i],
1822  e_outarray =
1823  outarray + m_coeff_offset[i]);
1824  }
1825 }
1826 
1827 /**
1828  * This function smooth a field after some calculaitons which have
1829  * been done elementally.
1830  *
1831  * @param field An array containing the field in physical space
1832  *
1833  */
1835 {
1836  boost::ignore_unused(field);
1837  // Do nothing unless the method is implemented in the appropriate
1838  // class, i.e. ContField1D,ContField2D, etc.
1839 
1840  // So far it has been implemented just for ContField2D and
1841  // ContField3DHomogeneous1D
1842 
1843  // Block in case users try the smoothing with a modal expansion.
1844  // Maybe a different techique for the smoothing require
1845  // implementation for modal basis.
1846 
1847  ASSERTL0((*m_exp)[0]->GetBasisType(0) == LibUtilities::eGLL_Lagrange ||
1848  (*m_exp)[0]->GetBasisType(0) == LibUtilities::eGauss_Lagrange,
1849  "Smoothing is currently not allowed unless you are using "
1850  "a nodal base for efficiency reasons. The implemented "
1851  "smoothing technique requires the mass matrix inversion "
1852  "which is trivial just for GLL_LAGRANGE_SEM and "
1853  "GAUSS_LAGRANGE_SEMexpansions.");
1854 }
1855 
1856 /**
1857  * This function assembles the block diagonal matrix
1858  * \f$\underline{\boldsymbol{M}}^e\f$, which is the
1859  * concatenation of the local matrices
1860  * \f$\boldsymbol{M}^e\f$ of the type \a mtype, that is
1861  *
1862  * \f[
1863  * \underline{\boldsymbol{M}}^e = \left[
1864  * \begin{array}{cccc}
1865  * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
1866  * 0 & \boldsymbol{M}^2 & 0 & 0 \\
1867  * 0 & 0 & \ddots & 0 \\
1868  * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array}\right].\f]
1869  *
1870  * @param mtype the type of matrix to be assembled
1871  * @param scalar an optional parameter
1872  * @param constant an optional parameter
1873  */
1875  const GlobalMatrixKey &gkey)
1876 {
1877  int i, cnt1;
1878  int n_exp = 0;
1879  DNekScalMatSharedPtr loc_mat;
1880  DNekScalBlkMatSharedPtr BlkMatrix;
1881  map<int, int> elmt_id;
1883 
1885  {
1886  for (i = 0; i < (*m_exp).size(); ++i)
1887  {
1888  if ((*m_exp)[i]->DetShapeType() == ShapeType)
1889  {
1890  elmt_id[n_exp++] = i;
1891  }
1892  }
1893  }
1894  else
1895  {
1896  n_exp = (*m_exp).size();
1897  for (i = 0; i < n_exp; ++i)
1898  {
1899  elmt_id[i] = i;
1900  }
1901  }
1902 
1903  Array<OneD, unsigned int> nrows(n_exp);
1904  Array<OneD, unsigned int> ncols(n_exp);
1905 
1906  switch (gkey.GetMatrixType())
1907  {
1908  case StdRegions::eBwdTrans:
1909  {
1910  // set up an array of integers for block matrix construction
1911  for (i = 0; i < n_exp; ++i)
1912  {
1913  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1914  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1915  }
1916  }
1917  break;
1919  {
1920  // set up an array of integers for block matrix construction
1921  for (i = 0; i < n_exp; ++i)
1922  {
1923  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1924  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1925  }
1926  }
1927  break;
1928  case StdRegions::eMass:
1929  case StdRegions::eInvMass:
1933  {
1934  // set up an array of integers for block matrix construction
1935  for (i = 0; i < n_exp; ++i)
1936  {
1937  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1938  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1939  }
1940  }
1941  break;
1942 
1944  {
1945  // set up an array of integers for block matrix construction
1946  for (i = 0; i < n_exp; ++i)
1947  {
1948  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1949  ncols[i] =
1950  (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1951  }
1952  }
1953  break;
1955  {
1956  // set up an array of integers for block matrix construction
1957  for (i = 0; i < n_exp; ++i)
1958  {
1959  nrows[i] =
1960  (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1961  ncols[i] =
1962  (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1963  }
1964  }
1965  break;
1966  default:
1967  {
1969  "Global Matrix creation not defined for this "
1970  "type of matrix");
1971  }
1972  }
1973 
1974  MatrixStorage blkmatStorage = eDIAGONAL;
1975  BlkMatrix = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nrows, ncols,
1976  blkmatStorage);
1977 
1978  int nvarcoeffs = gkey.GetNVarCoeffs();
1979  int eid;
1980  Array<OneD, NekDouble> varcoeffs_wk;
1981 
1982  for (i = cnt1 = 0; i < n_exp; ++i)
1983  {
1984  // need to be initialised with zero size for non
1985  // variable coefficient case
1986  StdRegions::VarCoeffMap varcoeffs;
1987 
1988  eid = elmt_id[i];
1989  if (nvarcoeffs > 0)
1990  {
1991  for (auto &x : gkey.GetVarCoeffs())
1992  {
1993  varcoeffs[x.first] = x.second + m_phys_offset[eid];
1994  }
1995  }
1996 
1997  LocalRegions::MatrixKey matkey(
1998  gkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(), *(*m_exp)[eid],
1999  gkey.GetConstFactors(), varcoeffs);
2000 
2001  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>(
2002  (*m_exp)[elmt_id.find(i)->second])
2003  ->GetLocMatrix(matkey);
2004  BlkMatrix->SetBlock(i, i, loc_mat);
2005  }
2006 
2007  return BlkMatrix;
2008 }
2009 
2011  const GlobalMatrixKey &gkey)
2012 {
2013  auto matrixIter = m_blockMat->find(gkey);
2014 
2015  if (matrixIter == m_blockMat->end())
2016  {
2017  return ((*m_blockMat)[gkey] = GenBlockMatrix(gkey));
2018  }
2019  else
2020  {
2021  return matrixIter->second;
2022  }
2023 }
2024 
2025 // Routines for continous matrix solution
2026 /**
2027  * This operation is equivalent to the evaluation of
2028  * \f$\underline{\boldsymbol{M}}^e\boldsymbol{\hat{u}}_l\f$, that is,
2029  * \f[ \left[
2030  * \begin{array}{cccc}
2031  * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
2032  * 0 & \boldsymbol{M}^2 & 0 & 0 \\
2033  * 0 & 0 & \ddots & 0 \\
2034  * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array} \right]
2035  *\left [ \begin{array}{c}
2036  * \boldsymbol{\hat{u}}^{1} \\
2037  * \boldsymbol{\hat{u}}^{2} \\
2038  * \vdots \\
2039  * \boldsymbol{\hat{u}}^{{{N_{\mathrm{el}}}}} \end{array} \right ]\f]
2040  * where \f$\boldsymbol{M}^e\f$ are the local matrices of type
2041  * specified by the key \a mkey. The decoupling of the local matrices
2042  * allows for a local evaluation of the operation. However, rather than
2043  * a local matrix-vector multiplication, the local operations are
2044  * evaluated as implemented in the function
2045  * StdRegions#StdExpansion#GeneralMatrixOp.
2046  *
2047  * @param mkey This key uniquely defines the type matrix
2048  * required for the operation.
2049  * @param inarray The vector \f$\boldsymbol{\hat{u}}_l\f$ of
2050  * size \f$N_{\mathrm{eof}}\f$.
2051  * @param outarray The resulting vector of size
2052  * \f$N_{\mathrm{eof}}\f$.
2053  */
2055  const Array<OneD, const NekDouble> &inarray,
2056  Array<OneD, NekDouble> &outarray)
2057 {
2058  int nvarcoeffs = gkey.GetNVarCoeffs();
2059 
2060  if ((nvarcoeffs == 0) && (gkey.GetMatrixType() == StdRegions::eHelmholtz))
2061  {
2062  // initialise if required
2063  if (m_collections.size() &&
2065  {
2066  for (int i = 0; i < m_collections.size(); ++i)
2067  {
2068  m_collections[i].Initialise(Collections::eHelmholtz,
2069  gkey.GetConstFactors());
2070  }
2072  }
2073  else
2074  {
2075  for (int i = 0; i < m_collections.size(); ++i)
2076  {
2077  m_collections[i].CheckFactors(Collections::eHelmholtz,
2078  gkey.GetConstFactors(),
2079  m_coll_phys_offset[i]);
2080  }
2081  }
2082 
2084  for (int i = 0; i < m_collections.size(); ++i)
2085  {
2086  m_collections[i].ApplyOperator(
2088  tmp = outarray + m_coll_coeff_offset[i]);
2089  }
2090  }
2091  else
2092  {
2093  Array<OneD, NekDouble> tmp_outarray;
2094  for (int i = 0; i < (*m_exp).size(); ++i)
2095  {
2096  // need to be initialised with zero size for non
2097  // variable coefficient case
2098  StdRegions::VarCoeffMap varcoeffs;
2099 
2100  if (nvarcoeffs > 0)
2101  {
2102  for (auto &x : gkey.GetVarCoeffs())
2103  {
2104  varcoeffs[x.first] = x.second + m_phys_offset[i];
2105  }
2106  }
2107 
2109  gkey.GetMatrixType(), (*m_exp)[i]->DetShapeType(),
2110  *((*m_exp)[i]), gkey.GetConstFactors(), varcoeffs);
2111 
2112  (*m_exp)[i]->GeneralMatrixOp(
2113  inarray + m_coeff_offset[i],
2114  tmp_outarray = outarray + m_coeff_offset[i], mkey);
2115  }
2116  }
2117 }
2118 
2119 /**
2120  * Retrieves local matrices from each expansion in the expansion list
2121  * and combines them together to generate a global matrix system.
2122  * @param mkey Matrix key for the matrix to be generated.
2123  * @param locToGloMap Local to global mapping.
2124  * @returns Shared pointer to the generated global matrix.
2125  */
2127  const GlobalMatrixKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2128 {
2129  int i, j, n, gid1, gid2, cntdim1, cntdim2;
2130  NekDouble sign1, sign2;
2131  DNekScalMatSharedPtr loc_mat;
2132 
2133  unsigned int glob_rows = 0;
2134  unsigned int glob_cols = 0;
2135  unsigned int loc_rows = 0;
2136  unsigned int loc_cols = 0;
2137 
2138  bool assembleFirstDim = false;
2139  bool assembleSecondDim = false;
2140 
2141  switch (mkey.GetMatrixType())
2142  {
2143  case StdRegions::eBwdTrans:
2144  {
2145  glob_rows = m_npoints;
2146  glob_cols = locToGloMap->GetNumGlobalCoeffs();
2147 
2148  assembleFirstDim = false;
2149  assembleSecondDim = true;
2150  }
2151  break;
2153  {
2154  glob_rows = locToGloMap->GetNumGlobalCoeffs();
2155  glob_cols = m_npoints;
2156 
2157  assembleFirstDim = true;
2158  assembleSecondDim = false;
2159  }
2160  break;
2161  case StdRegions::eMass:
2165  {
2166  glob_rows = locToGloMap->GetNumGlobalCoeffs();
2167  glob_cols = locToGloMap->GetNumGlobalCoeffs();
2168 
2169  assembleFirstDim = true;
2170  assembleSecondDim = true;
2171  }
2172  break;
2173  default:
2174  {
2176  "Global Matrix creation not defined for this "
2177  "type of matrix");
2178  }
2179  }
2180 
2181  COOMatType spcoomat;
2182  CoordType coord;
2183 
2184  int nvarcoeffs = mkey.GetNVarCoeffs();
2185  int eid;
2186 
2187  // fill global matrix
2188  for (n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
2189  {
2190  // need to be initialised with zero size for non
2191  // variable coefficient case
2192  StdRegions::VarCoeffMap varcoeffs;
2193 
2194  eid = n;
2195  if (nvarcoeffs > 0)
2196  {
2197  for (auto &x : mkey.GetVarCoeffs())
2198  {
2199  varcoeffs[x.first] = x.second + m_phys_offset[eid];
2200  }
2201  }
2202 
2203  LocalRegions::MatrixKey matkey(
2204  mkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(),
2205  *((*m_exp)[eid]), mkey.GetConstFactors(), varcoeffs);
2206 
2207  loc_mat =
2208  std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])
2209  ->GetLocMatrix(matkey);
2210 
2211  loc_rows = loc_mat->GetRows();
2212  loc_cols = loc_mat->GetColumns();
2213 
2214  for (i = 0; i < loc_rows; ++i)
2215  {
2216  if (assembleFirstDim)
2217  {
2218  gid1 = locToGloMap->GetLocalToGlobalMap(cntdim1 + i);
2219  sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
2220  }
2221  else
2222  {
2223  gid1 = cntdim1 + i;
2224  sign1 = 1.0;
2225  }
2226 
2227  for (j = 0; j < loc_cols; ++j)
2228  {
2229  if (assembleSecondDim)
2230  {
2231  gid2 = locToGloMap->GetLocalToGlobalMap(cntdim2 + j);
2232  sign2 = locToGloMap->GetLocalToGlobalSign(cntdim2 + j);
2233  }
2234  else
2235  {
2236  gid2 = cntdim2 + j;
2237  sign2 = 1.0;
2238  }
2239 
2240  // sparse matrix fill
2241  coord = make_pair(gid1, gid2);
2242  if (spcoomat.count(coord) == 0)
2243  {
2244  spcoomat[coord] = sign1 * sign2 * (*loc_mat)(i, j);
2245  }
2246  else
2247  {
2248  spcoomat[coord] += sign1 * sign2 * (*loc_mat)(i, j);
2249  }
2250  }
2251  }
2252  cntdim1 += loc_rows;
2253  cntdim2 += loc_cols;
2254  }
2255 
2257  glob_cols, spcoomat);
2258 }
2259 
2261  const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2262 {
2263  int i, j, n, gid1, gid2, loc_lda, eid;
2264  NekDouble sign1, sign2, value;
2265  DNekScalMatSharedPtr loc_mat;
2266 
2267  int totDofs = locToGloMap->GetNumGlobalCoeffs();
2268  int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
2269 
2270  unsigned int rows = totDofs - NumDirBCs;
2271  unsigned int cols = totDofs - NumDirBCs;
2272  NekDouble zero = 0.0;
2273 
2274  DNekMatSharedPtr Gmat;
2275  int bwidth = locToGloMap->GetFullSystemBandWidth();
2276 
2277  int nvarcoeffs = mkey.GetNVarCoeffs();
2278  MatrixStorage matStorage;
2279 
2280  map<int, RobinBCInfoSharedPtr> RobinBCInfo = GetRobinBCInfo();
2281 
2282  switch (mkey.GetMatrixType())
2283  {
2284  // case for all symmetric matices
2287  if ((2 * (bwidth + 1)) < rows)
2288  {
2291  rows, cols, zero, matStorage, bwidth, bwidth);
2292  }
2293  else
2294  {
2295  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
2297  rows, cols, zero, matStorage);
2298  }
2299 
2300  break;
2301  default: // Assume general matrix - currently only set up
2302  // for full invert
2303  {
2304  matStorage = eFULL;
2305  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows, cols, zero,
2306  matStorage);
2307  }
2308  }
2309 
2310  // fill global symmetric matrix
2311  for (n = 0; n < (*m_exp).size(); ++n)
2312  {
2313  // need to be initialised with zero size for non
2314  // variable coefficient case
2315  StdRegions::VarCoeffMap varcoeffs;
2316 
2317  eid = n;
2318  if (nvarcoeffs > 0)
2319  {
2320  for (auto &x : mkey.GetVarCoeffs())
2321  {
2322  varcoeffs[x.first] = x.second + m_phys_offset[eid];
2323  }
2324  }
2325 
2326  LocalRegions::MatrixKey matkey(
2327  mkey.GetMatrixType(), (*m_exp)[eid]->DetShapeType(),
2328  *((*m_exp)[eid]), mkey.GetConstFactors(), varcoeffs);
2329 
2330  loc_mat =
2331  std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])
2332  ->GetLocMatrix(matkey);
2333 
2334  if (RobinBCInfo.count(n) != 0) // add robin mass matrix
2335  {
2337 
2338  // declare local matrix from scaled matrix.
2339  int rows = loc_mat->GetRows();
2340  int cols = loc_mat->GetColumns();
2341  const NekDouble *dat = loc_mat->GetRawPtr();
2342  DNekMatSharedPtr new_mat =
2344  Blas::Dscal(rows * cols, loc_mat->Scale(), new_mat->GetRawPtr(), 1);
2345 
2346  // add local matrix contribution
2347  for (rBC = RobinBCInfo.find(n)->second; rBC; rBC = rBC->next)
2348  {
2349  (*m_exp)[n]->AddRobinMassMatrix(
2350  rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
2351  }
2352 
2353  NekDouble one = 1.0;
2354  // redeclare loc_mat to point to new_mat plus the scalar.
2355  loc_mat =
2357  }
2358 
2359  loc_lda = loc_mat->GetColumns();
2360 
2361  for (i = 0; i < loc_lda; ++i)
2362  {
2363  gid1 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + i) -
2364  NumDirBCs;
2365  sign1 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + i);
2366  if (gid1 >= 0)
2367  {
2368  for (j = 0; j < loc_lda; ++j)
2369  {
2370  gid2 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] +
2371  j) -
2372  NumDirBCs;
2373  sign2 = locToGloMap->GetLocalToGlobalSign(
2374  m_coeff_offset[n] + j);
2375  if (gid2 >= 0)
2376  {
2377  // When global matrix is symmetric,
2378  // only add the value for the upper
2379  // triangular part in order to avoid
2380  // entries to be entered twice
2381  if ((matStorage == eFULL) || (gid2 >= gid1))
2382  {
2383  value = Gmat->GetValue(gid1, gid2) +
2384  sign1 * sign2 * (*loc_mat)(i, j);
2385  Gmat->SetValue(gid1, gid2, value);
2386  }
2387  }
2388  }
2389  }
2390  }
2391  }
2392 
2393  return Gmat;
2394 }
2395 
2396 /**
2397  * Consider a linear system
2398  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{f}\f$ to be solved. Dependent
2399  * on the solution method, this function constructs
2400  * - <b>The full linear system</b><BR>
2401  * A call to the function #GenGlobalLinSysFullDirect
2402  * - <b>The statically condensed linear system</b><BR>
2403  * A call to the function #GenGlobalLinSysStaticCond
2404  *
2405  * @param mkey A key which uniquely defines the global
2406  * matrix to be constructed.
2407  * @param locToGloMap Contains the mapping array and required
2408  * information for the transformation from
2409  * local to global degrees of freedom.
2410  * @return (A shared pointer to) the global linear system in
2411  * required format.
2412  */
2414  const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2415 {
2416  GlobalLinSysSharedPtr returnlinsys;
2417  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
2418 
2420 
2421  if (vType >= eSIZE_GlobalSysSolnType)
2422  {
2423  NEKERROR(ErrorUtil::efatal, "Matrix solution type not defined");
2424  }
2425  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
2426 
2427  return GetGlobalLinSysFactory().CreateInstance(vSolnType, mkey, vExpList,
2428  locToGloMap);
2429 }
2430 
2432  const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
2433 {
2434  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
2435  const map<int, RobinBCInfoSharedPtr> vRobinBCInfo = GetRobinBCInfo();
2436 
2438 
2439  if (vType >= eSIZE_GlobalSysSolnType)
2440  {
2441  NEKERROR(ErrorUtil::efatal, "Matrix solution type not defined");
2442  }
2443  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
2444 
2445  return GetGlobalLinSysFactory().CreateInstance(vSolnType, mkey, vExpList,
2446  locToGloMap);
2447 }
2448 
2449 /**
2450  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of
2451  * an expansion, this function evaluates the spectral/hp
2452  * expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
2453  * quadrature points \f$\boldsymbol{x}_i\f$. The operation
2454  * is evaluated locally by the elemental function
2455  * StdRegions#StdExpansion#BwdTrans.
2456  *
2457  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
2458  * containing the local coefficients
2459  * \f$\hat{u}_n^e\f$.
2460  * @param outarray The resulting physical values at the
2461  * quadrature points
2462  * \f$u^{\delta}(\boldsymbol{x}_i)\f$
2463  * will be stored in this array of size
2464  * \f$Q_{\mathrm{tot}}\f$.
2465  */
2466 
2468  Array<OneD, NekDouble> &outarray)
2469 {
2470  LibUtilities::Timer timer;
2471 
2472  if (m_expType == e0D)
2473  {
2474  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
2475  }
2476  else
2477  {
2478  // initialise if required
2480  {
2481  for (int i = 0; i < m_collections.size(); ++i)
2482  {
2483  m_collections[i].Initialise(Collections::eBwdTrans);
2484  }
2486  }
2487 
2488  LIKWID_MARKER_START("v_BwdTrans");
2489  timer.Start();
2490 
2492  for (int i = 0; i < m_collections.size(); ++i)
2493  {
2494  m_collections[i].ApplyOperator(
2496  tmp = outarray + m_coll_phys_offset[i]);
2497  }
2498 
2499  timer.Stop();
2500  LIKWID_MARKER_STOP("v_BwdTrans");
2501  }
2502 
2503  // Elapsed time
2504  timer.AccumulateRegion("Collections:BwdTrans", 10);
2505 }
2506 
2508  const Array<OneD, const NekDouble> &gloCoord)
2509 {
2510  return GetExp(GetExpIndex(gloCoord));
2511 }
2512 
2513 /**
2514  * @todo need a smarter search here that first just looks at bounding
2515  * vertices - suggest first seeing if point is within 10% of
2516  * region defined by vertices. The do point search.
2517  */
2519  NekDouble tol, bool returnNearestElmt, int cachedId,
2520  NekDouble maxDistance)
2521 {
2522  Array<OneD, NekDouble> Lcoords(gloCoord.size());
2523 
2524  return GetExpIndex(gloCoord, Lcoords, tol, returnNearestElmt, cachedId,
2525  maxDistance);
2526 }
2527 
2529  Array<OneD, NekDouble> &locCoords, NekDouble tol,
2530  bool returnNearestElmt, int cachedId,
2531  NekDouble maxDistance)
2532 {
2533  if (GetNumElmts() == 0)
2534  {
2535  return -1;
2536  }
2537 
2538  if (m_elmtToExpId.size() == 0)
2539  {
2540  // Loop in reverse order so that in case where using a
2541  // Homogeneous expansion it sets geometry ids to first part of
2542  // m_exp list. Otherwise will set to second (complex) expansion
2543  for (int i = (*m_exp).size() - 1; i >= 0; --i)
2544  {
2545  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2546  }
2547  }
2548 
2549  NekDouble nearpt = 1e6;
2550  NekDouble nearpt_min = 1e6;
2551  int min_id = -1;
2552  Array<OneD, NekDouble> savLocCoords(locCoords.size());
2553 
2554  if (cachedId >= 0 && cachedId < (*m_exp).size())
2555  {
2556  nearpt = 1e12;
2557  if ((*m_exp)[cachedId]->GetGeom()->MinMaxCheck(gloCoords) &&
2558  (*m_exp)[cachedId]->GetGeom()->ContainsPoint(gloCoords, locCoords,
2559  tol, nearpt))
2560  {
2561  return cachedId;
2562  }
2563  else if (returnNearestElmt && (nearpt < nearpt_min))
2564  {
2565  // If it does not lie within, keep track of which element
2566  // is nearest.
2567  min_id = cachedId;
2568  nearpt_min = nearpt;
2569  Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
2570  }
2571  }
2572 
2573  NekDouble x = (gloCoords.size() > 0 ? gloCoords[0] : 0.0);
2574  NekDouble y = (gloCoords.size() > 1 ? gloCoords[1] : 0.0);
2575  NekDouble z = (gloCoords.size() > 2 ? gloCoords[2] : 0.0);
2578  GetExp(0)->GetCoordim(), -1, x, y, z);
2579 
2580  // Get the list of elements whose bounding box contains the desired
2581  // point.
2582  std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
2583 
2584  // Check each element in turn to see if point lies within it.
2585  for (int i = 0; i < elmts.size(); ++i)
2586  {
2587  int id = m_elmtToExpId[elmts[i]];
2588  if (id == cachedId)
2589  {
2590  continue;
2591  }
2592  if ((*m_exp)[id]->GetGeom()->ContainsPoint(gloCoords, locCoords, tol,
2593  nearpt))
2594  {
2595  return id;
2596  }
2597  else if (returnNearestElmt && (nearpt < nearpt_min))
2598  {
2599  // If it does not lie within, keep track of which element
2600  // is nearest.
2601  min_id = id;
2602  nearpt_min = nearpt;
2603  Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
2604  }
2605  }
2606 
2607  // If the calling function is with just the nearest element, return
2608  // that. Otherwise return -1 to indicate no matching elemenet found.
2609  if (returnNearestElmt && nearpt_min <= maxDistance)
2610  {
2611 
2612  std::string msg =
2613  "Failed to find point within element to "
2614  "tolerance of " +
2615  boost::lexical_cast<std::string>(tol) + " using local point (" +
2616  boost::lexical_cast<std::string>(locCoords[0]) + "," +
2617  boost::lexical_cast<std::string>(locCoords[1]) + "," +
2618  boost::lexical_cast<std::string>(locCoords[1]) +
2619  ") in element: " + boost::lexical_cast<std::string>(min_id);
2620  WARNINGL1(false, msg.c_str());
2621 
2622  Vmath::Vcopy(locCoords.size(), savLocCoords, 1, locCoords, 1);
2623  return min_id;
2624  }
2625  else
2626  {
2627  return -1;
2628  }
2629 }
2630 
2631 /**
2632  * Given some coordinates, output the expansion field value at that
2633  * point
2634  */
2636  const Array<OneD, const NekDouble> &phys)
2637 {
2638  int dim = GetCoordim(0);
2639  ASSERTL0(dim == coords.size(), "Invalid coordinate dimension.");
2640 
2641  // Grab the element index corresponding to coords.
2642  Array<OneD, NekDouble> xi(dim);
2643  int elmtIdx = GetExpIndex(coords, xi);
2644  ASSERTL0(elmtIdx > 0, "Unable to find element containing point.");
2645 
2646  // Grab that element's physical storage.
2647  Array<OneD, NekDouble> elmtPhys = phys + m_phys_offset[elmtIdx];
2648 
2649  // Evaluate the element at the appropriate point.
2650  return (*m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
2651 }
2652 
2653 /**
2654  * Configures geometric info, such as tangent direction, on each
2655  * expansion.
2656  * @param graph2D Mesh
2657  */
2659 {
2660 }
2661 
2662 /**
2663  * @brief Reset geometry information, metrics, matrix managers and
2664  * geometry information.
2665  *
2666  * This routine clears all matrix managers and resets all geometry
2667  * information, which allows the geometry information to be dynamically
2668  * updated as the solver is run.
2669  */
2671 {
2672  // Reset matrix managers.
2674  LocalRegions::MatrixKey::opLess>::ClearManager();
2676  LocalRegions::MatrixKey::opLess>::ClearManager();
2677 
2678  // Loop over all elements and reset geometry information.
2679  for (int i = 0; i < m_exp->size(); ++i)
2680  {
2681  (*m_exp)[i]->GetGeom()->Reset(m_graph->GetCurvedEdges(),
2682  m_graph->GetCurvedFaces());
2683  }
2684 
2685  // Loop over all elements and rebuild geometric factors.
2686  for (int i = 0; i < m_exp->size(); ++i)
2687  {
2688  (*m_exp)[i]->Reset();
2689  }
2690 }
2691 
2692 /**
2693  * Write Tecplot Files Header
2694  * @param outfile Output file name.
2695  * @param var variables names
2696  */
2697 void ExpList::v_WriteTecplotHeader(std::ostream &outfile, std::string var)
2698 {
2699  if (GetNumElmts() == 0)
2700  {
2701  return;
2702  }
2703 
2704  int coordim = GetExp(0)->GetCoordim();
2705  char vars[3] = {'x', 'y', 'z'};
2706 
2707  if (m_expType == e3DH1D)
2708  {
2709  coordim += 1;
2710  }
2711  else if (m_expType == e3DH2D)
2712  {
2713  coordim += 2;
2714  }
2715 
2716  outfile << "Variables = x";
2717  for (int i = 1; i < coordim; ++i)
2718  {
2719  outfile << ", " << vars[i];
2720  }
2721 
2722  if (var.size() > 0)
2723  {
2724  outfile << ", " << var;
2725  }
2726 
2727  outfile << std::endl << std::endl;
2728 }
2729 
2730 /**
2731  * Write Tecplot Files Zone
2732  * @param outfile Output file name.
2733  * @param expansion Expansion that is considered
2734  */
2735 void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
2736 {
2737  int i, j;
2738  int coordim = GetCoordim(0);
2739  int nPoints = GetTotPoints();
2740  int nBases = (*m_exp)[0]->GetNumBases();
2741  int numBlocks = 0;
2742 
2744 
2745  if (expansion == -1)
2746  {
2747  nPoints = GetTotPoints();
2748 
2749  coords[0] = Array<OneD, NekDouble>(nPoints);
2750  coords[1] = Array<OneD, NekDouble>(nPoints);
2751  coords[2] = Array<OneD, NekDouble>(nPoints);
2752 
2753  GetCoords(coords[0], coords[1], coords[2]);
2754 
2755  for (i = 0; i < m_exp->size(); ++i)
2756  {
2757  int numInt = 1;
2758 
2759  for (j = 0; j < nBases; ++j)
2760  {
2761  numInt *= (*m_exp)[i]->GetNumPoints(j) - 1;
2762  }
2763 
2764  numBlocks += numInt;
2765  }
2766  }
2767  else
2768  {
2769  nPoints = (*m_exp)[expansion]->GetTotPoints();
2770 
2771  coords[0] = Array<OneD, NekDouble>(nPoints);
2772  coords[1] = Array<OneD, NekDouble>(nPoints);
2773  coords[2] = Array<OneD, NekDouble>(nPoints);
2774 
2775  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
2776 
2777  numBlocks = 1;
2778  for (j = 0; j < nBases; ++j)
2779  {
2780  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j) - 1;
2781  }
2782  }
2783 
2784  if (m_expType == e3DH1D)
2785  {
2786  nBases += 1;
2787  coordim += 1;
2788  int nPlanes = GetZIDs().size();
2789  NekDouble tmp = numBlocks * (nPlanes - 1.0) / nPlanes;
2790  numBlocks = (int)tmp;
2791  }
2792  else if (m_expType == e3DH2D)
2793  {
2794  nBases += 2;
2795  coordim += 1;
2796  }
2797 
2798  outfile << "Zone, N=" << nPoints << ", E=" << numBlocks << ", F=FEBlock";
2799 
2800  switch (nBases)
2801  {
2802  case 2:
2803  outfile << ", ET=QUADRILATERAL" << std::endl;
2804  break;
2805  case 3:
2806  outfile << ", ET=BRICK" << std::endl;
2807  break;
2808  default:
2809  NEKERROR(ErrorUtil::efatal, "Not set up for this type of output");
2810  break;
2811  }
2812 
2813  // Write out coordinates
2814  for (j = 0; j < coordim; ++j)
2815  {
2816  for (i = 0; i < nPoints; ++i)
2817  {
2818  outfile << coords[j][i] << " ";
2819  if (i % 1000 == 0 && i)
2820  {
2821  outfile << std::endl;
2822  }
2823  }
2824  outfile << std::endl;
2825  }
2826 }
2827 
2828 void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
2829 {
2830  int i, j, k, l;
2831  int nbase = (*m_exp)[0]->GetNumBases();
2832  int cnt = 0;
2833 
2834  std::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
2835 
2836  if (expansion != -1)
2837  {
2838  exp = std::shared_ptr<LocalRegions::ExpansionVector>(
2840  (*exp)[0] = (*m_exp)[expansion];
2841  }
2842 
2843  if (nbase == 2)
2844  {
2845  for (i = 0; i < (*exp).size(); ++i)
2846  {
2847  const int np0 = (*exp)[i]->GetNumPoints(0);
2848  const int np1 = (*exp)[i]->GetNumPoints(1);
2849 
2850  for (j = 1; j < np1; ++j)
2851  {
2852  for (k = 1; k < np0; ++k)
2853  {
2854  outfile << cnt + (j - 1) * np0 + k << " ";
2855  outfile << cnt + (j - 1) * np0 + k + 1 << " ";
2856  outfile << cnt + j * np0 + k + 1 << " ";
2857  outfile << cnt + j * np0 + k << endl;
2858  }
2859  }
2860 
2861  cnt += np0 * np1;
2862  }
2863  }
2864  else if (nbase == 3)
2865  {
2866  for (i = 0; i < (*exp).size(); ++i)
2867  {
2868  const int np0 = (*exp)[i]->GetNumPoints(0);
2869  const int np1 = (*exp)[i]->GetNumPoints(1);
2870  const int np2 = (*exp)[i]->GetNumPoints(2);
2871  const int np01 = np0 * np1;
2872 
2873  for (j = 1; j < np2; ++j)
2874  {
2875  for (k = 1; k < np1; ++k)
2876  {
2877  for (l = 1; l < np0; ++l)
2878  {
2879  outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l
2880  << " ";
2881  outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l + 1
2882  << " ";
2883  outfile << cnt + (j - 1) * np01 + k * np0 + l + 1
2884  << " ";
2885  outfile << cnt + (j - 1) * np01 + k * np0 + l << " ";
2886  outfile << cnt + j * np01 + (k - 1) * np0 + l << " ";
2887  outfile << cnt + j * np01 + (k - 1) * np0 + l + 1
2888  << " ";
2889  outfile << cnt + j * np01 + k * np0 + l + 1 << " ";
2890  outfile << cnt + j * np01 + k * np0 + l << endl;
2891  }
2892  }
2893  }
2894  cnt += np0 * np1 * np2;
2895  }
2896  }
2897  else
2898  {
2899  NEKERROR(ErrorUtil::efatal, "Not set up for this dimension");
2900  }
2901 }
2902 
2903 /**
2904  * Write Tecplot Files Field
2905  * @param outfile Output file name.
2906  * @param expansion Expansion that is considered
2907  */
2908 void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
2909 {
2910  if (expansion == -1)
2911  {
2912  int totpoints = GetTotPoints();
2913  if (m_physState == false)
2914  {
2916  }
2917 
2918  for (int i = 0; i < totpoints; ++i)
2919  {
2920  outfile << m_phys[i] << " ";
2921  if (i % 1000 == 0 && i)
2922  {
2923  outfile << std::endl;
2924  }
2925  }
2926  outfile << std::endl;
2927  }
2928  else
2929  {
2930  int nPoints = (*m_exp)[expansion]->GetTotPoints();
2931 
2932  for (int i = 0; i < nPoints; ++i)
2933  {
2934  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
2935  }
2936 
2937  outfile << std::endl;
2938  }
2939 }
2940 
2941 void ExpList::WriteVtkHeader(std::ostream &outfile)
2942 {
2943  outfile << "<?xml version=\"1.0\"?>" << endl;
2944  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
2945  << "byte_order=\"LittleEndian\">" << endl;
2946  outfile << " <UnstructuredGrid>" << endl;
2947 }
2948 
2949 void ExpList::WriteVtkFooter(std::ostream &outfile)
2950 {
2951  outfile << " </UnstructuredGrid>" << endl;
2952  outfile << "</VTKFile>" << endl;
2953 }
2954 
2955 void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion,
2956  int istrip)
2957 {
2958  boost::ignore_unused(istrip);
2959  int i, j, k;
2960  int nbase = (*m_exp)[expansion]->GetNumBases();
2961  int ntot = (*m_exp)[expansion]->GetTotPoints();
2962  int nquad[3];
2963 
2964  int ntotminus = 1;
2965  for (i = 0; i < nbase; ++i)
2966  {
2967  nquad[i] = (*m_exp)[expansion]->GetNumPoints(i);
2968  ntotminus *= (nquad[i] - 1);
2969  }
2970 
2971  Array<OneD, NekDouble> coords[3];
2972  coords[0] = Array<OneD, NekDouble>(ntot, 0.0);
2973  coords[1] = Array<OneD, NekDouble>(ntot, 0.0);
2974  coords[2] = Array<OneD, NekDouble>(ntot, 0.0);
2975  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
2976 
2977  outfile << " <Piece NumberOfPoints=\"" << ntot << "\" NumberOfCells=\""
2978  << ntotminus << "\">" << endl;
2979  outfile << " <Points>" << endl;
2980  outfile << " <DataArray type=\"Float64\" "
2981  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
2982  outfile << " ";
2983  for (i = 0; i < ntot; ++i)
2984  {
2985  for (j = 0; j < 3; ++j)
2986  {
2987  outfile << setprecision(8) << scientific << (float)coords[j][i]
2988  << " ";
2989  }
2990  outfile << endl;
2991  }
2992  outfile << endl;
2993  outfile << " </DataArray>" << endl;
2994  outfile << " </Points>" << endl;
2995  outfile << " <Cells>" << endl;
2996  outfile << " <DataArray type=\"Int32\" "
2997  << "Name=\"connectivity\" format=\"ascii\">" << endl;
2998 
2999  int ns = 0; // pow(2,dim) for later usage
3000  string ostr;
3001  switch (m_expType)
3002  {
3003  case e1D:
3004  {
3005  ns = 2;
3006  ostr = "3 ";
3007  for (i = 0; i < nquad[0] - 1; ++i)
3008  {
3009  outfile << i << " " << i + 1 << endl;
3010  }
3011  }
3012  break;
3013  case e2D:
3014  {
3015  ns = 4;
3016  ostr = "9 ";
3017  for (i = 0; i < nquad[0] - 1; ++i)
3018  {
3019  for (j = 0; j < nquad[1] - 1; ++j)
3020  {
3021  outfile << j * nquad[0] + i << " " << j * nquad[0] + i + 1
3022  << " " << (j + 1) * nquad[0] + i + 1 << " "
3023  << (j + 1) * nquad[0] + i << endl;
3024  }
3025  }
3026  }
3027  break;
3028  case e3D:
3029  {
3030  ns = 8;
3031  ostr = "12 ";
3032  for (i = 0; i < nquad[0] - 1; ++i)
3033  {
3034  for (j = 0; j < nquad[1] - 1; ++j)
3035  {
3036  for (k = 0; k < nquad[2] - 1; ++k)
3037  {
3038  outfile
3039  << k * nquad[0] * nquad[1] + j * nquad[0] + i << " "
3040  << k * nquad[0] * nquad[1] + j * nquad[0] + i + 1
3041  << " "
3042  << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] +
3043  i + 1
3044  << " "
3045  << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] + i
3046  << " "
3047  << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] + i
3048  << " "
3049  << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] +
3050  i + 1
3051  << " "
3052  << (k + 1) * nquad[0] * nquad[1] +
3053  (j + 1) * nquad[0] + i + 1
3054  << " "
3055  << (k + 1) * nquad[0] * nquad[1] +
3056  (j + 1) * nquad[0] + i
3057  << " " << endl;
3058  }
3059  }
3060  }
3061  }
3062  break;
3063  default:
3064  break;
3065  }
3066 
3067  outfile << endl;
3068  outfile << " </DataArray>" << endl;
3069  outfile << " <DataArray type=\"Int32\" "
3070  << "Name=\"offsets\" format=\"ascii\">" << endl;
3071  for (i = 0; i < ntotminus; ++i)
3072  {
3073  outfile << i * ns + ns << " ";
3074  }
3075  outfile << endl;
3076  outfile << " </DataArray>" << endl;
3077  outfile << " <DataArray type=\"UInt8\" "
3078  << "Name=\"types\" format=\"ascii\">" << endl;
3079  for (i = 0; i < ntotminus; ++i)
3080  {
3081  outfile << ostr;
3082  }
3083  outfile << endl;
3084  outfile << " </DataArray>" << endl;
3085  outfile << " </Cells>" << endl;
3086  outfile << " <PointData>" << endl;
3087 }
3088 
3089 void ExpList::WriteVtkPieceFooter(std::ostream &outfile, int expansion)
3090 {
3091  boost::ignore_unused(expansion);
3092  outfile << " </PointData>" << endl;
3093  outfile << " </Piece>" << endl;
3094 }
3095 
3096 void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
3097  std::string var)
3098 {
3099  int i;
3100  int nq = (*m_exp)[expansion]->GetTotPoints();
3101 
3102  // printing the fields of that zone
3103  outfile << " <DataArray type=\"Float64\" Name=\"" << var << "\">"
3104  << endl;
3105  outfile << " ";
3106  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
3107  for (i = 0; i < nq; ++i)
3108  {
3109  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i])
3110  << " ";
3111  }
3112  outfile << endl;
3113  outfile << " </DataArray>" << endl;
3114 }
3115 
3116 /**
3117  * Given a spectral/hp approximation
3118  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
3119  * (which should be contained in #m_phys), this function calculates the
3120  * \f$L_\infty\f$ error of this approximation with respect to an exact
3121  * solution. The local distribution of the quadrature points allows an
3122  * elemental evaluation of this operation through the functions
3123  * StdRegions#StdExpansion#Linf.
3124  *
3125  * The exact solution, also evaluated at the quadrature
3126  * points, should be contained in the variable #m_phys of
3127  * the ExpList object \a Sol.
3128  *
3129  * @param soln A 1D array, containing the discrete
3130  * evaluation of the exact solution at the
3131  * quadrature points in its array #m_phys.
3132  * @return The \f$L_\infty\f$ error of the approximation.
3133  */
3135  const Array<OneD, const NekDouble> &soln)
3136 {
3137  NekDouble err = 0.0;
3138 
3139  if (soln == NullNekDouble1DArray)
3140  {
3141  err = Vmath::Vmax(m_npoints, inarray, 1);
3142  }
3143  else
3144  {
3145  for (int i = 0; i < m_npoints; ++i)
3146  {
3147  err = max(err, abs(inarray[i] - soln[i]));
3148  }
3149  }
3150 
3151  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
3152 
3153  return err;
3154 }
3155 
3156 /**
3157  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
3158  * evaluated at the quadrature points (which should be contained in
3159  * #m_phys), this function calculates the \f$L_2\f$ error of this
3160  * approximation with respect to an exact solution. The local
3161  * distribution of the quadrature points allows an elemental evaluation
3162  * of this operation through the functions StdRegions#StdExpansion#L2.
3163  *
3164  * The exact solution, also evaluated at the quadrature points, should
3165  * be contained in the variable #m_phys of the ExpList object \a Sol.
3166  *
3167  * @param Sol An ExpList, containing the discrete
3168  * evaluation of the exact solution at the
3169  * quadrature points in its array #m_phys.
3170  * @return The \f$L_2\f$ error of the approximation.
3171  */
3173  const Array<OneD, const NekDouble> &soln)
3174 {
3175  NekDouble err = 0.0, errl2;
3176  int i;
3177 
3178  if (soln == NullNekDouble1DArray)
3179  {
3180  for (i = 0; i < (*m_exp).size(); ++i)
3181  {
3182  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
3183  err += errl2 * errl2;
3184  }
3185  }
3186  else
3187  {
3188  for (i = 0; i < (*m_exp).size(); ++i)
3189  {
3190  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
3191  soln + m_phys_offset[i]);
3192  err += errl2 * errl2;
3193  }
3194  }
3195 
3196  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
3197 
3198  return sqrt(err);
3199 }
3200 
3201 /**
3202  * The integration is evaluated locally, that is
3203  * \f[\int
3204  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
3205  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
3206  * where the integration over the separate elements is done by the
3207  * function StdRegions#StdExpansion#Integral, which discretely
3208  * evaluates the integral using Gaussian quadrature.
3209  *
3210  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
3211  * containing the values of the function
3212  * \f$f(\boldsymbol{x})\f$ at the quadrature
3213  * points \f$\boldsymbol{x}_i\f$.
3214  * @return The value of the discretely evaluated integral
3215  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
3216  */
3218 {
3219  NekDouble sum = 0.0;
3220  int i = 0;
3221 
3222  for (i = 0; i < (*m_exp).size(); ++i)
3223  {
3224  sum += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
3225  }
3226  m_comm->GetRowComm()->AllReduce(sum, LibUtilities::ReduceSum);
3227 
3228  return sum;
3229 }
3230 
3232  const Array<OneD, Array<OneD, NekDouble>> &inarray)
3233 {
3234  NekDouble flux = 0.0;
3235  int i = 0;
3236  int j;
3237 
3238  for (i = 0; i < (*m_exp).size(); ++i)
3239  {
3240  Array<OneD, Array<OneD, NekDouble>> tmp(inarray.size());
3241  for (j = 0; j < inarray.size(); ++j)
3242  {
3243  tmp[j] = Array<OneD, NekDouble>(inarray[j] + m_phys_offset[i]);
3244  }
3245  flux += (*m_exp)[i]->VectorFlux(tmp);
3246  }
3247 
3248  return flux;
3249 }
3250 
3252 {
3254  "This method is not defined or valid for this class type");
3255  Array<OneD, NekDouble> NoEnergy(1, 0.0);
3256  return NoEnergy;
3257 }
3258 
3260 {
3262  "This method is not defined or valid for this class type");
3264  return trans;
3265 }
3266 
3268 {
3270  "This method is not defined or valid for this class type");
3271  NekDouble len = 0.0;
3272  return len;
3273 }
3274 
3276 {
3277  boost::ignore_unused(lhom);
3279  "This method is not defined or valid for this class type");
3280 }
3281 
3283 {
3285  "This method is not defined or valid for this class type");
3286  Array<OneD, unsigned int> NoModes(1);
3287  return NoModes;
3288 }
3289 
3291 {
3293  "This method is not defined or valid for this class type");
3294  Array<OneD, unsigned int> NoModes(1);
3295  return NoModes;
3296 }
3297 
3299 {
3301  "This method is not defined or valid for this class type");
3302 }
3303 
3304 void ExpList::ExtractFileBCs(const std::string &fileName,
3306  const std::string &varName,
3307  const std::shared_ptr<ExpList> locExpList)
3308 {
3309  string varString = fileName.substr(0, fileName.find_last_of("."));
3310  int j, k, len = varString.length();
3311  varString = varString.substr(len - 1, len);
3312 
3313  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
3314  std::vector<std::vector<NekDouble>> FieldData;
3315 
3316  std::string ft = LibUtilities::FieldIO::GetFileType(fileName, comm);
3319  ft, comm, m_session->GetSharedFilesystem());
3320 
3321  f->Import(fileName, FieldDef, FieldData);
3322 
3323  bool found = false;
3324  for (j = 0; j < FieldDef.size(); ++j)
3325  {
3326  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
3327  {
3328  if (FieldDef[j]->m_fields[k] == varName)
3329  {
3330  // Copy FieldData into locExpList
3331  locExpList->ExtractDataToCoeffs(FieldDef[j], FieldData[j],
3332  FieldDef[j]->m_fields[k],
3333  locExpList->UpdateCoeffs());
3334  found = true;
3335  }
3336  }
3337  }
3338 
3339  ASSERTL0(found, "Could not find variable '" + varName +
3340  "' in file boundary condition " + fileName);
3341  locExpList->BwdTrans(locExpList->GetCoeffs(), locExpList->UpdatePhys());
3342 }
3343 
3344 /**
3345  * Given a spectral/hp approximation
3346  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
3347  * (which should be contained in #m_phys), this function calculates the
3348  * \f$H^1_2\f$ error of this approximation with respect to an exact
3349  * solution. The local distribution of the quadrature points allows an
3350  * elemental evaluation of this operation through the functions
3351  * StdRegions#StdExpansion#H1.
3352  *
3353  * The exact solution, also evaluated at the quadrature points, should
3354  * be contained in the variable #m_phys of the ExpList object \a Sol.
3355  *
3356  * @param soln An 1D array, containing the discrete evaluation
3357  * of the exact solution at the quadrature points.
3358  *
3359  * @return The \f$H^1_2\f$ error of the approximation.
3360  */
3362  const Array<OneD, const NekDouble> &soln)
3363 {
3364  NekDouble err = 0.0, errh1;
3365  int i;
3366 
3367  for (i = 0; i < (*m_exp).size(); ++i)
3368  {
3369  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
3370  soln + m_phys_offset[i]);
3371  err += errh1 * errh1;
3372  }
3373 
3374  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
3375 
3376  return sqrt(err);
3377 }
3378 
3380  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
3381  int NumHomoDir, Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
3382  std::vector<NekDouble> &HomoLen, bool homoStrips,
3383  std::vector<unsigned int> &HomoSIDs, std::vector<unsigned int> &HomoZIDs,
3384  std::vector<unsigned int> &HomoYIDs)
3385 {
3386  int startenum = (int)LibUtilities::eSegment;
3387  int endenum = (int)LibUtilities::eHexahedron;
3388  int s = 0;
3390 
3391  ASSERTL1(NumHomoDir == HomoBasis.size(),
3392  "Homogeneous basis is not the same length as NumHomoDir");
3393  ASSERTL1(NumHomoDir == HomoLen.size(),
3394  "Homogeneous length vector is not the same length as NumHomDir");
3395 
3396  // count number of shapes
3397  switch ((*m_exp)[0]->GetShapeDimension())
3398  {
3399  case 1:
3400  startenum = (int)LibUtilities::eSegment;
3401  endenum = (int)LibUtilities::eSegment;
3402  break;
3403  case 2:
3404  startenum = (int)LibUtilities::eTriangle;
3405  endenum = (int)LibUtilities::eQuadrilateral;
3406  break;
3407  case 3:
3408  startenum = (int)LibUtilities::eTetrahedron;
3409  endenum = (int)LibUtilities::eHexahedron;
3410  break;
3411  }
3412 
3413  for (s = startenum; s <= endenum; ++s)
3414  {
3415  std::vector<unsigned int> elementIDs;
3416  std::vector<LibUtilities::BasisType> basis;
3417  std::vector<unsigned int> numModes;
3418  std::vector<std::string> fields;
3419 
3420  bool first = true;
3421  bool UniOrder = true;
3422  int n;
3423 
3424  shape = (LibUtilities::ShapeType)s;
3425 
3426  for (int i = 0; i < (*m_exp).size(); ++i)
3427  {
3428  if ((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
3429  {
3430  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
3431  if (first)
3432  {
3433  for (int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3434  {
3435  basis.push_back(
3436  (*m_exp)[i]->GetBasis(j)->GetBasisType());
3437  numModes.push_back(
3438  (*m_exp)[i]->GetBasis(j)->GetNumModes());
3439  }
3440 
3441  // add homogeneous direction details if defined
3442  for (n = 0; n < NumHomoDir; ++n)
3443  {
3444  basis.push_back(HomoBasis[n]->GetBasisType());
3445  numModes.push_back(HomoBasis[n]->GetNumModes());
3446  }
3447 
3448  first = false;
3449  }
3450  else
3451  {
3452  ASSERTL0(
3453  (*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],
3454  "Routine is not set up for multiple bases definitions");
3455 
3456  for (int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3457  {
3458  numModes.push_back(
3459  (*m_exp)[i]->GetBasis(j)->GetNumModes());
3460  if (numModes[j] !=
3461  (*m_exp)[i]->GetBasis(j)->GetNumModes())
3462  {
3463  UniOrder = false;
3464  }
3465  }
3466  // add homogeneous direction details if defined
3467  for (n = 0; n < NumHomoDir; ++n)
3468  {
3469  numModes.push_back(HomoBasis[n]->GetNumModes());
3470  }
3471  }
3472  }
3473  }
3474 
3475  if (elementIDs.size() > 0)
3476  {
3479  AllocateSharedPtr(shape, elementIDs, basis, UniOrder,
3480  numModes, fields, NumHomoDir, HomoLen,
3481  homoStrips, HomoSIDs, HomoZIDs, HomoYIDs);
3482  fielddef.push_back(fdef);
3483  }
3484  }
3485 }
3486 
3487 //
3488 // Virtual functions
3489 //
3490 std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::
3492 {
3493  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
3494  v_GetFieldDefinitions(returnval);
3495  return returnval;
3496 }
3497 
3499  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
3500 {
3501  GeneralGetFieldDefinitions(fielddef);
3502 }
3503 
3504 // Append the element data listed in elements
3505 // fielddef->m_ElementIDs onto fielddata
3508  std::vector<NekDouble> &fielddata)
3509 {
3510  v_AppendFieldData(fielddef, fielddata, m_coeffs);
3511 }
3512 
3515  std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
3516 {
3517  int i;
3518  // Determine mapping from element ids to location in
3519  // expansion list
3520  // Determine mapping from element ids to location in
3521  // expansion list
3522  map<int, int> ElmtID_to_ExpID;
3523  for (i = 0; i < (*m_exp).size(); ++i)
3524  {
3525  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3526  }
3527 
3528  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
3529  {
3530  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
3531  int datalen = (*m_exp)[eid]->GetNcoeffs();
3532  fielddata.insert(fielddata.end(), &coeffs[m_coeff_offset[eid]],
3533  &coeffs[m_coeff_offset[eid]] + datalen);
3534  }
3535 }
3536 
3537 /// Extract the data in fielddata into the coeffs
3540  std::vector<NekDouble> &fielddata, std::string &field,
3541  Array<OneD, NekDouble> &coeffs)
3542 {
3543  v_ExtractDataToCoeffs(fielddef, fielddata, field, coeffs);
3544 }
3545 
3547  const std::shared_ptr<ExpList> &fromExpList,
3548  const Array<OneD, const NekDouble> &fromCoeffs,
3549  Array<OneD, NekDouble> &toCoeffs)
3550 {
3551  v_ExtractCoeffsToCoeffs(fromExpList, fromCoeffs, toCoeffs);
3552 }
3553 
3554 /**
3555  * @brief Extract data from raw field data into expansion list.
3556  *
3557  * @param fielddef Field definitions.
3558  * @param fielddata Data for associated field.
3559  * @param field Field variable name.
3560  * @param coeffs Resulting coefficient array.
3561  */
3564  std::vector<NekDouble> &fielddata, std::string &field,
3565  Array<OneD, NekDouble> &coeffs)
3566 {
3567  int i, expId;
3568  int offset = 0;
3569  int modes_offset = 0;
3570  int datalen = fielddata.size() / fielddef->m_fields.size();
3571 
3572  // Find data location according to field definition
3573  for (i = 0; i < fielddef->m_fields.size(); ++i)
3574  {
3575  if (fielddef->m_fields[i] == field)
3576  {
3577  break;
3578  }
3579  offset += datalen;
3580  }
3581 
3582  ASSERTL0(i != fielddef->m_fields.size(),
3583  "Field (" + field + ") not found in file.");
3584 
3585  if (m_elmtToExpId.size() == 0)
3586  {
3587  // Loop in reverse order so that in case where using a
3588  // Homogeneous expansion it sets geometry ids to first part of
3589  // m_exp list. Otherwise will set to second (complex) expansion
3590  for (i = (*m_exp).size() - 1; i >= 0; --i)
3591  {
3592  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3593  }
3594  }
3595 
3596  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
3597  {
3598  // Reset modes_offset in the case where all expansions of
3599  // the same order.
3600  if (fielddef->m_uniOrder == true)
3601  {
3602  modes_offset = 0;
3603  }
3604 
3606  fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
3607 
3608  const int elmtId = fielddef->m_elementIDs[i];
3609  auto eIt = m_elmtToExpId.find(elmtId);
3610 
3611  if (eIt == m_elmtToExpId.end())
3612  {
3613  offset += datalen;
3614  modes_offset += (*m_exp)[0]->GetNumBases();
3615  continue;
3616  }
3617 
3618  expId = eIt->second;
3619 
3620  bool sameBasis = true;
3621  for (int j = 0; j < fielddef->m_basis.size(); ++j)
3622  {
3623  if (fielddef->m_basis[j] != (*m_exp)[expId]->GetBasisType(j))
3624  {
3625  sameBasis = false;
3626  break;
3627  }
3628  }
3629 
3630  if (datalen == (*m_exp)[expId]->GetNcoeffs() && sameBasis)
3631  {
3632  Vmath::Vcopy(datalen, &fielddata[offset], 1,
3633  &coeffs[m_coeff_offset[expId]], 1);
3634  }
3635  else
3636  {
3637  (*m_exp)[expId]->ExtractDataToCoeffs(
3638  &fielddata[offset], fielddef->m_numModes, modes_offset,
3639  &coeffs[m_coeff_offset[expId]], fielddef->m_basis);
3640  }
3641 
3642  offset += datalen;
3643  modes_offset += (*m_exp)[0]->GetNumBases();
3644  }
3645 
3646  return;
3647 }
3648 
3650  const std::shared_ptr<ExpList> &fromExpList,
3651  const Array<OneD, const NekDouble> &fromCoeffs,
3652  Array<OneD, NekDouble> &toCoeffs)
3653 {
3654  int i;
3655  int offset = 0;
3656 
3657  for (i = 0; i < (*m_exp).size(); ++i)
3658  {
3659  std::vector<unsigned int> nummodes;
3660  vector<LibUtilities::BasisType> basisTypes;
3661  for (int j = 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
3662  {
3663  nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
3664  basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
3665  }
3666 
3667  (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes, 0,
3668  &toCoeffs[m_coeff_offset[i]],
3669  basisTypes);
3670 
3671  offset += fromExpList->GetExp(i)->GetNcoeffs();
3672  }
3673 }
3674 
3675 /**
3676  * Get the weight value on boundaries
3677  */
3679  Array<OneD, NekDouble> &weightJump)
3680 {
3681  size_t nTracePts = weightAver.size();
3682  // average for interior traces
3683  for (int i = 0; i < nTracePts; ++i)
3684  {
3685  weightAver[i] = 0.5;
3686  weightJump[i] = 1.0;
3687  }
3688  FillBwdWithBwdWeight(weightAver, weightJump);
3689 }
3690 
3692  const Array<OneD, const NekDouble> &CircCentre,
3693  Array<OneD, Array<OneD, NekDouble>> &outarray)
3694 {
3695  int npts;
3696 
3697  int MFdim = 3;
3698  int nq = outarray[0].size() / MFdim;
3699 
3700  // Assume whole array is of same coordinate dimension
3701  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
3702 
3703  Array<OneD, Array<OneD, NekDouble>> MFloc(MFdim * coordim);
3704  // Process each expansion.
3705  for (int i = 0; i < m_exp->size(); ++i)
3706  {
3707  npts = (*m_exp)[i]->GetTotPoints();
3708 
3709  for (int j = 0; j < MFdim * coordim; ++j)
3710  {
3711  MFloc[j] = Array<OneD, NekDouble>(npts, 0.0);
3712  }
3713 
3714  // MF from LOCALREGIONS
3715  (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
3716  (*m_exp)[i]->GetPointsKeys(), MMFdir, CircCentre, MFloc);
3717 
3718  // Get the physical data offset for this expansion.
3719  for (int j = 0; j < MFdim; ++j)
3720  {
3721  for (int k = 0; k < coordim; ++k)
3722  {
3723  Vmath::Vcopy(npts, &MFloc[j * coordim + k][0], 1,
3724  &outarray[j][k * nq + m_phys_offset[i]], 1);
3725  }
3726  }
3727  }
3728 }
3729 
3730 /**
3731  * @brief Generate vector v such that v[i] = scalar1 if i is in the
3732  * element < ElementID. Otherwise, v[i] = scalar2.
3733  *
3734  */
3735 void ExpList::GenerateElementVector(const int ElementID,
3736  const NekDouble scalar1,
3737  const NekDouble scalar2,
3738  Array<OneD, NekDouble> &outarray)
3739 {
3740  int npoints_e;
3741  NekDouble coeff;
3742 
3743  Array<OneD, NekDouble> outarray_e;
3744 
3745  for (int i = 0; i < (*m_exp).size(); ++i)
3746  {
3747  npoints_e = (*m_exp)[i]->GetTotPoints();
3748 
3749  if (i <= ElementID)
3750  {
3751  coeff = scalar1;
3752  }
3753  else
3754  {
3755  coeff = scalar2;
3756  }
3757 
3758  outarray_e = Array<OneD, NekDouble>(npoints_e, coeff);
3759  Vmath::Vcopy(npoints_e, &outarray_e[0], 1, &outarray[m_phys_offset[i]],
3760  1);
3761  }
3762 }
3763 
3766 {
3768  "This method is not defined or valid for this class type");
3770  return result;
3771 }
3772 
3773 std::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
3774 {
3775  boost::ignore_unused(i);
3777  "This method is not defined or valid for this class type");
3778  static std::shared_ptr<ExpList> result;
3779  return result;
3780 }
3781 
3782 /**
3783  * Upwind the left and right states given by the Arrays Fwd and Bwd
3784  * using the vector quantity Vec and ouput the upwinded value in the
3785  * array upwind.
3786  *
3787  * @param Vec Velocity field.
3788  * @param Fwd Left state.
3789  * @param Bwd Right state.
3790  * @param Upwind Output vector.
3791  */
3793  const Array<OneD, const NekDouble> &Fwd,
3794  const Array<OneD, const NekDouble> &Bwd,
3795  Array<OneD, NekDouble> &Upwind)
3796 {
3797  switch (m_expType)
3798  {
3799  case e1D:
3800  {
3801  int i, j, k, e_npoints, offset;
3802  Array<OneD, NekDouble> normals;
3803  NekDouble Vn;
3804 
3805  // Assume whole array is of same coordimate dimension
3806  int coordim = GetCoordim(0);
3807 
3808  ASSERTL1(Vec.size() >= coordim,
3809  "Input vector does not have sufficient dimensions to "
3810  "match coordim");
3811 
3812  // Process each expansion
3813  for (i = 0; i < m_exp->size(); ++i)
3814  {
3815  // Get the number of points in the expansion and the normals.
3816  e_npoints = (*m_exp)[i]->GetNumPoints(0);
3817  normals = (*m_exp)[i]->GetPhysNormals();
3818 
3819  // Get the physical data offset of the expansion in m_phys.
3820  offset = m_phys_offset[i];
3821 
3822  // Compute each data point.
3823  for (j = 0; j < e_npoints; ++j)
3824  {
3825  // Calculate normal velocity.
3826  Vn = 0.0;
3827  for (k = 0; k < coordim; ++k)
3828  {
3829  Vn += Vec[k][offset + j] * normals[k * e_npoints + j];
3830  }
3831 
3832  // Upwind based on direction of normal velocity.
3833  if (Vn > 0.0)
3834  {
3835  Upwind[offset + j] = Fwd[offset + j];
3836  }
3837  else
3838  {
3839  Upwind[offset + j] = Bwd[offset + j];
3840  }
3841  }
3842  }
3843  }
3844  break;
3845  default:
3847  "This method is not defined or valid for this class type");
3848  break;
3849  }
3850 }
3851 
3852 /**
3853  * One-dimensional upwind.
3854  * \see ExpList::Upwind(
3855  * const Array<OneD, const Array<OneD, NekDouble> >,
3856  * const Array<OneD, const NekDouble>,
3857  * const Array<OneD, const NekDouble>,
3858  * Array<OneD, NekDouble>, int)
3859  *
3860  * @param Vn Velocity field.
3861  * @param Fwd Left state.
3862  * @param Bwd Right state.
3863  * @param Upwind Output vector.
3864  */
3866  const Array<OneD, const NekDouble> &Fwd,
3867  const Array<OneD, const NekDouble> &Bwd,
3868  Array<OneD, NekDouble> &Upwind)
3869 {
3870  ASSERTL1(Vn.size() >= m_npoints, "Vn is not of sufficient length");
3871  ASSERTL1(Fwd.size() >= m_npoints, "Fwd is not of sufficient length");
3872  ASSERTL1(Bwd.size() >= m_npoints, "Bwd is not of sufficient length");
3873  ASSERTL1(Upwind.size() >= m_npoints, "Upwind is not of sufficient length");
3874 
3875  // Process each point in the expansion.
3876  for (int j = 0; j < m_npoints; ++j)
3877  {
3878  // Upwind based on one-dimensional velocity.
3879  if (Vn[j] > 0.0)
3880  {
3881  Upwind[j] = Fwd[j];
3882  }
3883  else
3884  {
3885  Upwind[j] = Bwd[j];
3886  }
3887  }
3888 }
3889 
3890 std::shared_ptr<ExpList> &ExpList::v_GetTrace()
3891 {
3893  "This method is not defined or valid for this class type");
3894  static std::shared_ptr<ExpList> returnVal;
3895  return returnVal;
3896 }
3897 
3898 std::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
3899 {
3901  "This method is not defined or valid for this class type");
3902  static std::shared_ptr<AssemblyMapDG> result;
3903  return result;
3904 }
3905 
3907 {
3908  return GetTraceMap()->GetBndCondIDToGlobalTraceID();
3909 }
3910 
3912 {
3914  "This method is not defined or valid for this class type");
3915  static std::vector<bool> result;
3916  return result;
3917 }
3918 
3919 /**
3920  * @brief Helper function to re-align face to a given orientation.
3921  */
3922 void AlignFace(const StdRegions::Orientation orient, const int nquad1,
3923  const int nquad2, const Array<OneD, const NekDouble> &in,
3925 {
3926  // Copy transpose.
3927  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
3931  {
3932  for (int i = 0; i < nquad2; ++i)
3933  {
3934  for (int j = 0; j < nquad1; ++j)
3935  {
3936  out[i * nquad1 + j] = in[j * nquad2 + i];
3937  }
3938  }
3939  }
3940  else
3941  {
3942  for (int i = 0; i < nquad2; ++i)
3943  {
3944  for (int j = 0; j < nquad1; ++j)
3945  {
3946  out[i * nquad1 + j] = in[i * nquad1 + j];
3947  }
3948  }
3949  }
3950 
3951  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
3955  {
3956  // Reverse x direction
3957  for (int i = 0; i < nquad2; ++i)
3958  {
3959  for (int j = 0; j < nquad1 / 2; ++j)
3960  {
3961  swap(out[i * nquad1 + j], out[i * nquad1 + nquad1 - j - 1]);
3962  }
3963  }
3964  }
3965 
3966  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
3970  {
3971  // Reverse y direction
3972  for (int j = 0; j < nquad1; ++j)
3973  {
3974  for (int i = 0; i < nquad2 / 2; ++i)
3975  {
3976  swap(out[i * nquad1 + j], out[(nquad2 - i - 1) * nquad1 + j]);
3977  }
3978  }
3979  }
3980 }
3981 
3982 /**
3983  * For each local element, copy the normals stored in the element list
3984  * into the array \a normals.
3985  * @param normals Multidimensional array in which to copy normals
3986  * to. Must have dimension equal to or larger than
3987  * the spatial dimension of the elements.
3988  */
3990 {
3991  int i, j, k, e_npoints, offset;
3993 
3994  // Assume whole array is of same coordinate dimension
3995  int coordim = GetCoordim(0);
3996 
3997  ASSERTL1(normals.size() >= coordim,
3998  "Output vector does not have sufficient dimensions to "
3999  "match coordim");
4000 
4001  switch (m_expType)
4002  {
4003  case e0D:
4004  {
4005  // Process each expansion.
4006  for (i = 0; i < m_exp->size(); ++i)
4007  {
4008  LocalRegions::ExpansionSharedPtr loc_exp = (*m_exp)[i];
4009 
4011  loc_exp->GetLeftAdjacentElementExp();
4012 
4013  // Get the number of points and normals for this expansion.
4014  e_npoints = 1;
4015  locnormals = loc_elmt->GetTraceNormal(
4016  loc_exp->GetLeftAdjacentElementTrace());
4017 
4018  // Get the physical data offset for this expansion.
4019  offset = m_phys_offset[i];
4020 
4021  // Process each point in the expansion.
4022  for (j = 0; j < e_npoints; ++j)
4023  {
4024  // Process each spatial dimension and copy the
4025  // values into the output array.
4026  for (k = 0; k < coordim; ++k)
4027  {
4028  normals[k][offset] = locnormals[k][0];
4029  }
4030  }
4031  }
4032  }
4033  break;
4034  case e1D:
4035  {
4039 
4040  for (i = 0; i < m_exp->size(); ++i)
4041  {
4042  LocalRegions::ExpansionSharedPtr loc_exp = (*m_exp)[i];
4043 
4045  loc_exp->GetLeftAdjacentElementExp();
4046 
4047  int edgeNumber = loc_exp->GetLeftAdjacentElementTrace();
4048 
4049  // Get the number of points and normals for this expansion.
4050  e_npoints = (*m_exp)[i]->GetNumPoints(0);
4051 
4052  locnormals = loc_elmt->GetTraceNormal(edgeNumber);
4053  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
4054  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
4055 
4056  if (e_nmodes != loc_nmodes)
4057  {
4058  if (loc_exp->GetRightAdjacentElementTrace() >= 0)
4059  {
4061  loc_exp->GetRightAdjacentElementExp();
4062 
4063  int EdgeNumber =
4064  loc_exp->GetRightAdjacentElementTrace();
4065 
4066  // Serial case: right element is connected so we can
4067  // just grab that normal.
4068  locnormals = loc_elmt->GetTraceNormal(EdgeNumber);
4069 
4070  offset = m_phys_offset[i];
4071 
4072  // Process each point in the expansion.
4073  for (j = 0; j < e_npoints; ++j)
4074  {
4075  // Process each spatial dimension and
4076  // copy the values into the output
4077  // array.
4078  for (k = 0; k < coordim; ++k)
4079  {
4080  normals[k][offset + j] = -locnormals[k][j];
4081  }
4082  }
4083  }
4084  else
4085  {
4086  // Parallel case: need to interpolate normal.
4087  Array<OneD, Array<OneD, NekDouble>> normal(coordim);
4088 
4089  for (int p = 0; p < coordim; ++p)
4090  {
4091  normal[p] = Array<OneD, NekDouble>(e_npoints, 0.0);
4092  LibUtilities::PointsKey to_key =
4093  loc_exp->GetBasis(0)->GetPointsKey();
4094  LibUtilities::PointsKey from_key =
4095  loc_elmt->GetBasis(0)->GetPointsKey();
4096  LibUtilities::Interp1D(from_key, locnormals[p],
4097  to_key, normal[p]);
4098  }
4099 
4100  offset = m_phys_offset[i];
4101 
4102  // Process each point in the expansion.
4103  for (j = 0; j < e_npoints; ++j)
4104  {
4105  // Process each spatial dimension and copy the
4106  // values into the output array.
4107  for (k = 0; k < coordim; ++k)
4108  {
4109  normals[k][offset + j] = normal[k][j];
4110  }
4111  }
4112  }
4113  }
4114  else
4115  {
4116  // Get the physical data offset for this expansion.
4117  offset = m_phys_offset[i];
4118 
4119  // Process each point in the expansion.
4120  for (j = 0; j < e_npoints; ++j)
4121  {
4122  // Process each spatial dimension and copy the values
4123  // into the output array.
4124  for (k = 0; k < coordim; ++k)
4125  {
4126  normals[k][offset + j] = locnormals[k][j];
4127  }
4128  }
4129  }
4130  }
4131  }
4132  break;
4133  case e2D:
4134  {
4136 
4137  // Process each expansion.
4138  for (i = 0; i < m_exp->size(); ++i)
4139  {
4140  LocalRegions::ExpansionSharedPtr traceExp = (*m_exp)[i];
4142  traceExp->GetLeftAdjacentElementExp();
4143 
4144  // Get the number of points and normals for this expansion.
4145  int faceNum = traceExp->GetLeftAdjacentElementTrace();
4146  int offset = m_phys_offset[i];
4147 
4148  const Array<OneD, const Array<OneD, NekDouble>> &locNormals =
4149  exp3D->GetTraceNormal(faceNum);
4150 
4151  // Project normals from 3D element onto the same
4152  // orientation as the trace expansion.
4153  StdRegions::Orientation orient = exp3D->GetTraceOrient(faceNum);
4154 
4155  int fromid0, fromid1;
4156 
4158  {
4159  fromid0 = 0;
4160  fromid1 = 1;
4161  }
4162  else
4163  {
4164  fromid0 = 1;
4165  fromid1 = 0;
4166  }
4167 
4168  LibUtilities::BasisKey faceBasis0 =
4169  exp3D->GetTraceBasisKey(faceNum, fromid0);
4170  LibUtilities::BasisKey faceBasis1 =
4171  exp3D->GetTraceBasisKey(faceNum, fromid1);
4172  LibUtilities::BasisKey traceBasis0 =
4173  traceExp->GetBasis(0)->GetBasisKey();
4174  LibUtilities::BasisKey traceBasis1 =
4175  traceExp->GetBasis(1)->GetBasisKey();
4176 
4177  const int faceNq0 = faceBasis0.GetNumPoints();
4178  const int faceNq1 = faceBasis1.GetNumPoints();
4179 
4180  Array<OneD, int> faceids;
4181  exp3D->ReOrientTracePhysMap(orient, faceids, faceNq0, faceNq1);
4182 
4183  Array<OneD, NekDouble> traceNormals(faceNq0 * faceNq1);
4184 
4185  for (j = 0; j < coordim; ++j)
4186  {
4187  Vmath::Scatr(faceNq0 * faceNq1, locNormals[j], faceids,
4188  traceNormals);
4189 
4191  faceBasis0.GetPointsKey(), faceBasis1.GetPointsKey(),
4192  traceNormals, traceBasis0.GetPointsKey(),
4193  traceBasis1.GetPointsKey(), tmp = normals[j] + offset);
4194  }
4195  }
4196  }
4197  break;
4198  default:
4199  {
4201  "This method is not defined or valid for this class type");
4202  }
4203  }
4204 }
4205 
4207  Array<OneD, NekDouble> &lengthsBwd)
4208 {
4209  int e_npoints;
4210 
4211  Array<OneD, NekDouble> locLeng;
4214  Array<OneD, int> LRbndnumbs(2);
4216  lengLR[0] = lengthsFwd;
4217  lengLR[1] = lengthsBwd;
4221  int e_npoints0 = -1;
4222  if (m_expType == e1D)
4223  {
4224  for (int i = 0; i < m_exp->size(); ++i)
4225  {
4226  loc_exp = (*m_exp)[i];
4227  int offset = m_phys_offset[i];
4228 
4229  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
4230  e_npoints = (*m_exp)[i]->GetNumPoints(0);
4231  if (e_npoints0 < e_npoints)
4232  {
4233  for (int nlr = 0; nlr < 2; nlr++)
4234  {
4235  lengintp[nlr] = Array<OneD, NekDouble>(e_npoints, 0.0);
4236  }
4237  e_npoints0 = e_npoints;
4238  }
4239 
4240  LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4241  LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4242 
4243  LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4244  LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4245  for (int nlr = 0; nlr < 2; ++nlr)
4246  {
4247  Vmath::Zero(e_npoints0, lengintp[nlr], 1);
4248  lengAdd[nlr] = lengintp[nlr];
4249  int bndNumber = LRbndnumbs[nlr];
4250  loc_elmt = LRelmts[nlr];
4251  if (bndNumber >= 0)
4252  {
4253  locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
4254  lengAdd[nlr] = locLeng;
4255 
4256  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
4257  if (e_nmodes != loc_nmodes)
4258  {
4259  // Parallel case: need to interpolate.
4260  LibUtilities::PointsKey to_key =
4261  loc_exp->GetBasis(0)->GetPointsKey();
4262  LibUtilities::PointsKey from_key =
4263  loc_elmt->GetBasis(0)->GetPointsKey();
4264  LibUtilities::Interp1D(from_key, locLeng, to_key,
4265  lengintp[nlr]);
4266  lengAdd[nlr] = lengintp[nlr];
4267  }
4268  }
4269 
4270  for (int j = 0; j < e_npoints; ++j)
4271  {
4272  lengLR[nlr][offset + j] = lengAdd[nlr][j];
4273  }
4274  }
4275  }
4276  }
4277  else if (m_expType == e2D)
4278  {
4279  for (int i = 0; i < m_exp->size(); ++i)
4280  {
4281  loc_exp = (*m_exp)[i];
4282  int offset = m_phys_offset[i];
4283 
4284  LibUtilities::BasisKey traceBasis0 =
4285  loc_exp->GetBasis(0)->GetBasisKey();
4286  LibUtilities::BasisKey traceBasis1 =
4287  loc_exp->GetBasis(1)->GetBasisKey();
4288  const int TraceNq0 = traceBasis0.GetNumPoints();
4289  const int TraceNq1 = traceBasis1.GetNumPoints();
4290  e_npoints = TraceNq0 * TraceNq1;
4291  if (e_npoints0 < e_npoints)
4292  {
4293  for (int nlr = 0; nlr < 2; nlr++)
4294  {
4295  lengintp[nlr] = Array<OneD, NekDouble>(e_npoints, 0.0);
4296  }
4297  e_npoints0 = e_npoints;
4298  }
4299 
4300  LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4301  LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4302 
4303  LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4304  LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4305  for (int nlr = 0; nlr < 2; ++nlr)
4306  {
4307  Vmath::Zero(e_npoints0, lengintp[nlr], 1);
4308  int bndNumber = LRbndnumbs[nlr];
4309  loc_elmt = LRelmts[nlr];
4310  if (bndNumber >= 0)
4311  {
4312  locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
4313  // Project normals from 3D element onto the
4314  // same orientation as the trace expansion.
4315  StdRegions::Orientation orient =
4316  loc_elmt->GetTraceOrient(bndNumber);
4317 
4318  int fromid0, fromid1;
4320  {
4321  fromid0 = 0;
4322  fromid1 = 1;
4323  }
4324  else
4325  {
4326  fromid0 = 1;
4327  fromid1 = 0;
4328  }
4329 
4330  LibUtilities::BasisKey faceBasis0 =
4331  loc_elmt->GetTraceBasisKey(bndNumber, fromid0);
4332  LibUtilities::BasisKey faceBasis1 =
4333  loc_elmt->GetTraceBasisKey(bndNumber, fromid1);
4334  const int faceNq0 = faceBasis0.GetNumPoints();
4335  const int faceNq1 = faceBasis1.GetNumPoints();
4336  Array<OneD, NekDouble> alignedLeng(faceNq0 * faceNq1);
4337 
4338  AlignFace(orient, faceNq0, faceNq1, locLeng, alignedLeng);
4340  faceBasis0.GetPointsKey(), faceBasis1.GetPointsKey(),
4341  alignedLeng, traceBasis0.GetPointsKey(),
4342  traceBasis1.GetPointsKey(), lengintp[nlr]);
4343  }
4344 
4345  for (int j = 0; j < e_npoints; ++j)
4346  {
4347  lengLR[nlr][offset + j] = lengintp[nlr][j];
4348  }
4349  }
4350  }
4351  }
4352 }
4353 
4355  const Array<OneD, const NekDouble> &Fy,
4356  Array<OneD, NekDouble> &outarray)
4357 {
4358  boost::ignore_unused(Fx, Fy, outarray);
4360  "This method is not defined or valid for this class type");
4361 }
4362 
4364  Array<OneD, NekDouble> &outarray)
4365 {
4366  boost::ignore_unused(Fn, outarray);
4368  "This method is not defined or valid for this class type");
4369 }
4370 
4372  const Array<OneD, const NekDouble> &Bwd,
4373  Array<OneD, NekDouble> &outarray)
4374 {
4375  boost::ignore_unused(Fwd, Bwd, outarray);
4377  "This method is not defined or valid for this class type");
4378 }
4379 
4382 {
4383  boost::ignore_unused(Fwd, Bwd);
4385  "This method is not defined or valid for this class type");
4386 }
4387 
4390  Array<OneD, NekDouble> &Bwd, bool FillBnd,
4391  bool PutFwdInBwdOnBCs, bool DoExchange)
4392 {
4393  boost::ignore_unused(field, Fwd, Bwd, FillBnd, PutFwdInBwdOnBCs,
4394  DoExchange);
4396  "This method is not defined or valid for this class type");
4397 }
4398 
4401  bool PutFwdInBwdOnBCs)
4402 {
4403  boost::ignore_unused(Fwd, Bwd, PutFwdInBwdOnBCs);
4405  "This method is not defined or valid for this class type");
4406 }
4407 
4409  const Array<OneD, const NekDouble> &Bwd,
4410  Array<OneD, NekDouble> &field)
4411 {
4412  boost::ignore_unused(field, Fwd, Bwd);
4414  "v_AddTraceQuadPhysToField is not defined for this class type");
4415 }
4416 
4418  const Array<OneD, const NekDouble> &Fwd,
4420 {
4421  boost::ignore_unused(field, Fwd, Bwd);
4423  "v_AddTraceQuadPhysToOffDiag is not defined for this class");
4424 }
4425 
4427  const Array<OneD, const NekDouble> &Bwd,
4428  Array<OneD, NekDouble> &locTraceFwd,
4429  Array<OneD, NekDouble> &locTraceBwd)
4430 {
4431  boost::ignore_unused(Fwd, Bwd, locTraceFwd, locTraceBwd);
4433  "v_GetLocTraceFromTracePts is not defined for this class");
4434 }
4435 
4437 {
4439  "v_GetBndCondBwdWeight is not defined for this class type");
4440  static Array<OneD, NekDouble> tmp;
4441  return tmp;
4442 }
4443 
4444 void ExpList::v_SetBndCondBwdWeight(const int index, const NekDouble value)
4445 {
4446  boost::ignore_unused(index, value);
4448  "v_setBndCondBwdWeight is not defined for this class type");
4449 }
4450 
4451 const vector<bool> &ExpList::v_GetLeftAdjacentFaces(void) const
4452 {
4454  "This method is not defined or valid for this class type");
4455  static vector<bool> tmp;
4456  return tmp;
4457 }
4458 
4460 {
4461  boost::ignore_unused(outarray);
4463  "This method is not defined or valid for this class type");
4464 }
4465 
4467  Array<OneD, NekDouble> &outarray)
4468 {
4469  boost::ignore_unused(inarray, outarray);
4471  "This method is not defined or valid for this class type");
4472 }
4473 
4475  const Array<OneD, const NekDouble> &inarray,
4476  Array<OneD, NekDouble> &outarray)
4477 {
4478  boost::ignore_unused(inarray, outarray);
4480  "This method is not defined or valid for this class type");
4481 }
4482 
4484  Array<OneD, NekDouble> &outarray,
4485  const StdRegions::ConstFactorMap &factors,
4486  const StdRegions::VarCoeffMap &varcoeff,
4487  const MultiRegions::VarFactorsMap &varfactors,
4488  const Array<OneD, const NekDouble> &dirForcing,
4489  const bool PhysSpaceForcing)
4490 {
4491  boost::ignore_unused(inarray, outarray, factors, varcoeff, varfactors,
4492  dirForcing, PhysSpaceForcing);
4493  NEKERROR(ErrorUtil::efatal, "HelmSolve not implemented.");
4494 }
4495 
4497  const Array<OneD, Array<OneD, NekDouble>> &velocity,
4498  const Array<OneD, const NekDouble> &inarray,
4499  Array<OneD, NekDouble> &outarray, const NekDouble lambda,
4500  const Array<OneD, const NekDouble> &dirForcing)
4501 {
4502  boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4504  "This method is not defined or valid for this class type");
4505 }
4506 
4508  const Array<OneD, Array<OneD, NekDouble>> &velocity,
4509  const Array<OneD, const NekDouble> &inarray,
4510  Array<OneD, NekDouble> &outarray, const NekDouble lambda,
4511  const Array<OneD, const NekDouble> &dirForcing)
4512 {
4513  boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4515  "This method is not defined or valid for this class type");
4516 }
4517 
4519  Array<OneD, NekDouble> &outarray,
4520  bool Shuff, bool UnShuff)
4521 {
4522  boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4524  "This method is not defined or valid for this class type");
4525 }
4526 
4528  Array<OneD, NekDouble> &outarray,
4529  bool Shuff, bool UnShuff)
4530 {
4531  boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4533  "This method is not defined or valid for this class type");
4534 }
4535 
4537  const Array<OneD, NekDouble> &inarray2,
4538  Array<OneD, NekDouble> &outarray)
4539 {
4540  boost::ignore_unused(inarray1, inarray2, outarray);
4542  "This method is not defined or valid for this class type");
4543 }
4544 
4546  const Array<OneD, Array<OneD, NekDouble>> &inarray1,
4547  const Array<OneD, Array<OneD, NekDouble>> &inarray2,
4548  Array<OneD, Array<OneD, NekDouble>> &outarray)
4549 {
4550  boost::ignore_unused(inarray1, inarray2, outarray);
4552  "This method is not defined or valid for this class type");
4553 }
4554 
4556  const Array<OneD, NekDouble> &TotField, int BndID)
4557 {
4558  boost::ignore_unused(BndVals, TotField, BndID);
4560  "This method is not defined or valid for this class type");
4561 }
4562 
4565  Array<OneD, NekDouble> &outarray,
4566  int BndID)
4567 {
4568  boost::ignore_unused(V1, V2, outarray, BndID);
4570  "This method is not defined or valid for this class type");
4571 }
4572 
4575 {
4577  switch (GetCoordim(0))
4578  {
4579  case 1:
4580  {
4581  for (int i = 0; i < GetExpSize(); ++i)
4582  {
4583  (*m_exp)[i]->NormVectorIProductWRTBase(
4584  V[0] + GetPhys_Offset(i),
4585  tmp = outarray + GetCoeff_Offset(i));
4586  }
4587  }
4588  break;
4589  case 2:
4590  {
4591  for (int i = 0; i < GetExpSize(); ++i)
4592  {
4593  (*m_exp)[i]->NormVectorIProductWRTBase(
4594  V[0] + GetPhys_Offset(i), V[1] + GetPhys_Offset(i),
4595  tmp = outarray + GetCoeff_Offset(i));
4596  }
4597  }
4598  break;
4599  case 3:
4600  {
4601  for (int i = 0; i < GetExpSize(); ++i)
4602  {
4603  (*m_exp)[i]->NormVectorIProductWRTBase(
4604  V[0] + GetPhys_Offset(i), V[1] + GetPhys_Offset(i),
4605  V[2] + GetPhys_Offset(i),
4606  tmp = outarray + GetCoeff_Offset(i));
4607  }
4608  }
4609  break;
4610  default:
4611  NEKERROR(ErrorUtil::efatal, "Dimension not supported");
4612  break;
4613  }
4614 }
4615 
4617 {
4618  boost::ignore_unused(outarray);
4620  "This method is not defined or valid for this class type");
4621 }
4622 
4623 /**
4624  */
4626 {
4628  "This method is not defined or valid for this class type");
4629 }
4630 
4631 /**
4632  */
4634 {
4635  boost::ignore_unused(nreg);
4637  "This method is not defined or valid for this class type");
4638 }
4639 
4640 void ExpList::v_LocalToGlobal(bool useComm)
4641 {
4642  boost::ignore_unused(useComm);
4644  "This method is not defined or valid for this class type");
4645 }
4646 
4648  Array<OneD, NekDouble> &outarray, bool useComm)
4649 {
4650  boost::ignore_unused(inarray, outarray, useComm);
4652  "This method is not defined or valid for this class type");
4653 }
4654 
4656 {
4658  "This method is not defined or valid for this class type");
4659 }
4660 
4662  Array<OneD, NekDouble> &outarray)
4663 {
4664  boost::ignore_unused(inarray, outarray);
4666  "This method is not defined or valid for this class type");
4667 }
4668 
4670  Array<OneD, NekDouble> &outarray)
4671 {
4672  v_FwdTransLocalElmt(inarray, outarray);
4673 }
4674 
4675 /**
4676  * The operation is evaluated locally for every element by the function
4677  * StdRegions#StdExpansion#IProductWRTBase.
4678  *
4679  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
4680  * containing the values of the function
4681  * \f$f(\boldsymbol{x})\f$ at the quadrature
4682  * points \f$\boldsymbol{x}_i\f$.
4683  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
4684  * used to store the result.
4685  */
4687  Array<OneD, NekDouble> &outarray)
4688 {
4689  // initialise if required
4691  {
4692  for (int i = 0; i < m_collections.size(); ++i)
4693  {
4695  }
4697  }
4698 
4700  for (int i = 0; i < m_collections.size(); ++i)
4701  {
4702 
4704  inarray + m_coll_phys_offset[i],
4705  tmp = outarray + m_coll_coeff_offset[i]);
4706  }
4707 }
4708 
4709 /**
4710  * The operation is evaluated locally by the elemental
4711  * function StdRegions#StdExpansion#GetCoords.
4712  *
4713  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
4714  * will be stored in this array.
4715  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
4716  * will be stored in this array.
4717  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
4718  * will be stored in this array.
4719  */
4721  Array<OneD, NekDouble> &coord_1,
4722  Array<OneD, NekDouble> &coord_2)
4723 {
4724  if (GetNumElmts() == 0)
4725  {
4726  return;
4727  }
4728 
4729  int i;
4730  Array<OneD, NekDouble> e_coord_0;
4731  Array<OneD, NekDouble> e_coord_1;
4732  Array<OneD, NekDouble> e_coord_2;
4733 
4734  switch (GetExp(0)->GetCoordim())
4735  {
4736  case 1:
4737  for (i = 0; i < (*m_exp).size(); ++i)
4738  {
4739  e_coord_0 = coord_0 + m_phys_offset[i];
4740  (*m_exp)[i]->GetCoords(e_coord_0);
4741  }
4742  break;
4743  case 2:
4744  ASSERTL0(coord_1.size() != 0, "output coord_1 is not defined");
4745 
4746  for (i = 0; i < (*m_exp).size(); ++i)
4747  {
4748  e_coord_0 = coord_0 + m_phys_offset[i];
4749  e_coord_1 = coord_1 + m_phys_offset[i];
4750  (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1);
4751  }
4752  break;
4753  case 3:
4754  ASSERTL0(coord_1.size() != 0, "output coord_1 is not defined");
4755  ASSERTL0(coord_2.size() != 0, "output coord_2 is not defined");
4756 
4757  for (i = 0; i < (*m_exp).size(); ++i)
4758  {
4759  e_coord_0 = coord_0 + m_phys_offset[i];
4760  e_coord_1 = coord_1 + m_phys_offset[i];
4761  e_coord_2 = coord_2 + m_phys_offset[i];
4762  (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1, e_coord_2);
4763  }
4764  break;
4765  }
4766 }
4767 
4771 {
4772  (*m_exp)[eid]->GetCoords(xc0, xc1, xc2);
4773 }
4774 
4775 /**
4776  * @brief: Set up a normal along the trace elements between
4777  * two elements at elemental level
4778  *
4779  */
4781 {
4782  for (int i = 0; i < m_exp->size(); ++i)
4783  {
4784  for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
4785  {
4786  (*m_exp)[i]->ComputeTraceNormal(j);
4787  }
4788  }
4789 }
4790 
4791 /**
4792  */
4793 void ExpList::v_GetBndElmtExpansion(int i, std::shared_ptr<ExpList> &result,
4794  const bool DeclareCoeffPhysArrays)
4795 {
4796  boost::ignore_unused(i, result, DeclareCoeffPhysArrays);
4798  "This method is not defined or valid for this class type");
4799 }
4800 
4801 /**
4802  */
4804  const Array<OneD, NekDouble> &element,
4805  Array<OneD, NekDouble> &boundary)
4806 {
4807  int n, cnt;
4808  Array<OneD, NekDouble> tmp1, tmp2;
4810 
4811  Array<OneD, int> ElmtID, EdgeID;
4812  GetBoundaryToElmtMap(ElmtID, EdgeID);
4813 
4814  // Initialise result
4815  boundary =
4817 
4818  // Skip other boundary regions
4819  for (cnt = n = 0; n < i; ++n)
4820  {
4821  cnt += GetBndCondExpansions()[n]->GetExpSize();
4822  }
4823 
4824  int offsetBnd;
4825  int offsetElmt = 0;
4826  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4827  {
4828  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4829 
4830  elmt = GetExp(ElmtID[cnt + n]);
4831  elmt->GetTracePhysVals(
4832  EdgeID[cnt + n], GetBndCondExpansions()[i]->GetExp(n),
4833  tmp1 = element + offsetElmt, tmp2 = boundary + offsetBnd);
4834 
4835  offsetElmt += elmt->GetTotPoints();
4836  }
4837 }
4838 
4839 /**
4840  */
4842  const Array<OneD, const NekDouble> &phys,
4843  Array<OneD, NekDouble> &bndElmt)
4844 {
4845  int n, cnt, nq;
4846 
4847  Array<OneD, int> ElmtID, EdgeID;
4848  GetBoundaryToElmtMap(ElmtID, EdgeID);
4849 
4850  // Skip other boundary regions
4851  for (cnt = n = 0; n < i; ++n)
4852  {
4853  cnt += GetBndCondExpansions()[n]->GetExpSize();
4854  }
4855 
4856  // Count number of points
4857  int npoints = 0;
4858  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4859  {
4860  npoints += GetExp(ElmtID[cnt + n])->GetTotPoints();
4861  }
4862 
4863  // Initialise result
4864  bndElmt = Array<OneD, NekDouble>(npoints, 0.0);
4865 
4866  // Extract data
4867  int offsetPhys;
4868  int offsetElmt = 0;
4869  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4870  {
4871  nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
4872  offsetPhys = GetPhys_Offset(ElmtID[cnt + n]);
4873  Vmath::Vcopy(nq, &phys[offsetPhys], 1, &bndElmt[offsetElmt], 1);
4874  offsetElmt += nq;
4875  }
4876 }
4877 
4878 /**
4879  */
4881  const Array<OneD, const NekDouble> &phys,
4883 {
4884  int n, cnt;
4887 
4888  Array<OneD, int> ElmtID, EdgeID;
4889  GetBoundaryToElmtMap(ElmtID, EdgeID);
4890 
4891  // Initialise result
4892  bnd =
4894 
4895  // Skip other boundary regions
4896  for (cnt = n = 0; n < i; ++n)
4897  {
4898  cnt += GetBndCondExpansions()[n]->GetExpSize();
4899  }
4900 
4901  int offsetBnd;
4902  int offsetPhys;
4903  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4904  {
4905  offsetPhys = GetPhys_Offset(ElmtID[cnt + n]);
4906  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4907 
4908  elmt = GetExp(ElmtID[cnt + n]);
4909  elmt->GetTracePhysVals(EdgeID[cnt + n],
4910  GetBndCondExpansions()[i]->GetExp(n),
4911  phys + offsetPhys, tmp1 = bnd + offsetBnd);
4912  }
4913 }
4914 
4915 /**
4916  */
4918  Array<OneD, Array<OneD, NekDouble>> &normals)
4919 {
4920  int j, n, cnt, nq;
4921  int coordim = GetCoordim(0);
4924 
4925  Array<OneD, int> ElmtID, EdgeID;
4926  GetBoundaryToElmtMap(ElmtID, EdgeID);
4927 
4928  // Initialise result
4929  normals = Array<OneD, Array<OneD, NekDouble>>(coordim);
4930  for (j = 0; j < coordim; ++j)
4931  {
4932  normals[j] = Array<OneD, NekDouble>(
4933  GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
4934  }
4935 
4936  // Skip other boundary regions
4937  for (cnt = n = 0; n < i; ++n)
4938  {
4939  cnt += GetBndCondExpansions()[n]->GetExpSize();
4940  }
4941 
4942  int offset;
4943  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4944  {
4945  offset = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4946  nq = GetBndCondExpansions()[i]->GetExp(n)->GetTotPoints();
4947 
4948  elmt = GetExp(ElmtID[cnt + n]);
4949  const Array<OneD, const Array<OneD, NekDouble>> normalsElmt =
4950  elmt->GetTraceNormal(EdgeID[cnt + n]);
4951  // Copy to result
4952  for (j = 0; j < coordim; ++j)
4953  {
4954  Vmath::Vcopy(nq, normalsElmt[j], 1, tmp = normals[j] + offset, 1);
4955  }
4956  }
4957 }
4958 
4959 /**
4960  */
4962  Array<OneD, int> &EdgeID)
4963 {
4964  boost::ignore_unused(ElmtID, EdgeID);
4966  "This method is not defined or valid for this class type");
4967 }
4968 
4970  Array<OneD, NekDouble> &weightjmp)
4971 {
4972  boost::ignore_unused(weightave, weightjmp);
4973  NEKERROR(ErrorUtil::efatal, "v_FillBwdWithBwdWeight not defined");
4974 }
4975 
4978 {
4979  boost::ignore_unused(Fwd, Bwd);
4980  NEKERROR(ErrorUtil::efatal, "v_PeriodicBwdCopy not defined");
4981 }
4982 
4983 /**
4984  */
4987 {
4989  "This method is not defined or valid for this class type");
4991  return result;
4992 }
4993 
4994 /**
4995  */
4998 {
5000  "This method is not defined or valid for this class type");
5002  return result;
5003 }
5004 
5005 /**
5006  */
5008  const std::string varName,
5009  const NekDouble x2_in,
5010  const NekDouble x3_in)
5011 {
5012  boost::ignore_unused(time, varName, x2_in, x3_in);
5014  "This method is not defined or valid for this class type");
5015 }
5016 
5017 /**
5018  */
5019 map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
5020 {
5022  "This method is not defined or valid for this class type");
5023  static map<int, RobinBCInfoSharedPtr> result;
5024  return result;
5025 }
5026 
5027 /**
5028  */
5030  PeriodicMap &periodicEdges,
5031  PeriodicMap &periodicFaces)
5032 {
5033  boost::ignore_unused(periodicVerts, periodicEdges, periodicFaces);
5035  "This method is not defined or valid for this class type");
5036 }
5037 
5040  unsigned int regionId, const std::string &variable)
5041 {
5042  auto collectionIter = collection.find(regionId);
5043  ASSERTL1(collectionIter != collection.end(),
5044  "Unable to locate collection " +
5045  boost::lexical_cast<string>(regionId));
5046 
5048  (*collectionIter).second;
5049  auto conditionMapIter = bndCondMap->find(variable);
5050  ASSERTL1(conditionMapIter != bndCondMap->end(),
5051  "Unable to locate condition map.");
5052 
5053  const SpatialDomains::BoundaryConditionShPtr boundaryCondition =
5054  (*conditionMapIter).second;
5055 
5056  return boundaryCondition;
5057 }
5058 
5060 {
5061  boost::ignore_unused(n);
5063  "This method is not defined or valid for this class type");
5064  return NullExpListSharedPtr;
5065 }
5066 
5067 /**
5068  * @brief Construct collections of elements containing a single element
5069  * type and polynomial order from the list of expansions.
5070  */
5072 {
5074  vector<std::pair<LocalRegions::ExpansionSharedPtr, int>>>
5075  collections;
5076 
5077  // Set up initialisation flags
5079  std::vector<bool>(Collections::SIZE_OperatorType, true);
5080 
5081  // Figure out optimisation parameters if provided in
5082  // session file or default given
5084  m_session, (*m_exp)[0]->GetShapeDimension(), ImpType);
5085  // ImpType = colOpt.GetDefaultImplementationType();
5086 
5087  // turn on autotuning if explicitly specified in xml file
5088  // or command line option is set but only do optimisation
5089  // for volumetric elements (not boundary condition)
5090  bool autotuning = colOpt.IsUsingAutotuning();
5091  if ((autotuning == false) && (ImpType == Collections::eNoImpType))
5092  {
5093  // turn on autotuning if writeoptfile specified
5094  // if m_graph available
5095  if (m_session->GetUpdateOptFile() && m_graph)
5096  {
5097  // only turn on autotuning for volumetric elements
5098  // where Mesh Dimension is equal to the Shape
5099  // Dimension of element.
5100  if (m_graph->GetMeshDimension() == (*m_exp)[0]->GetShapeDimension())
5101  {
5102  autotuning = true;
5103  }
5104  }
5105  }
5106  bool verbose = (m_session->DefinesCmdLineArgument("verbose")) &&
5107  (m_session->GetComm()->GetRank() == 0);
5108  int collmax =
5109  (colOpt.GetMaxCollectionSize() > 0 ? colOpt.GetMaxCollectionSize()
5110  : 2 * m_exp->size());
5111 
5112  // clear vectors in case previously called
5113  m_collections.clear();
5114  m_coll_coeff_offset.clear();
5115  m_coll_phys_offset.clear();
5116 
5117  // Loop over expansions, and create collections for each element type
5118  for (int i = 0; i < m_exp->size(); ++i)
5119  {
5120  collections[(*m_exp)[i]->DetShapeType()].push_back(
5121  std::pair<LocalRegions::ExpansionSharedPtr, int>((*m_exp)[i], i));
5122  }
5123 
5124  for (auto &it : collections)
5125  {
5126  LocalRegions::ExpansionSharedPtr exp = it.second[0].first;
5127 
5128  Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(exp);
5129 
5130  vector<StdRegions::StdExpansionSharedPtr> collExp;
5131 
5132  int prevCoeffOffset = m_coeff_offset[it.second[0].second];
5133  int prevPhysOffset = m_phys_offset[it.second[0].second];
5134  int collcnt;
5135 
5136  m_coll_coeff_offset.push_back(prevCoeffOffset);
5137  m_coll_phys_offset.push_back(prevPhysOffset);
5138 
5139  int collsize = 0;
5140 
5141  if (it.second.size() == 1) // single element case
5142  {
5143  collExp.push_back(it.second[0].first);
5144 
5145  // if no Imp Type provided and No
5146  // setting in xml file. reset
5147  // impTypes using timings
5148  if (autotuning)
5149  {
5150  if (collExp.size() > collsize)
5151  {
5152  impTypes =
5153  colOpt.SetWithTimings(collExp, impTypes, verbose);
5154  collsize = collExp.size();
5155  }
5156  }
5157 
5158  Collections::Collection tmp(collExp, impTypes);
5159  m_collections.push_back(tmp);
5160  }
5161  else
5162  {
5163  // set up first geometry
5164  collExp.push_back(it.second[0].first);
5165  int prevnCoeff = it.second[0].first->GetNcoeffs();
5166  int prevnPhys = it.second[0].first->GetTotPoints();
5167  bool prevDeformed =
5168  it.second[0].first->GetMetricInfo()->GetGtype() ==
5170  collcnt = 1;
5171 
5172  for (int i = 1; i < it.second.size(); ++i)
5173  {
5174  int nCoeffs = it.second[i].first->GetNcoeffs();
5175  int nPhys = it.second[i].first->GetTotPoints();
5176  bool Deformed =
5177  it.second[i].first->GetMetricInfo()->GetGtype() ==
5179  int coeffOffset = m_coeff_offset[it.second[i].second];
5180  int physOffset = m_phys_offset[it.second[i].second];
5181 
5182  // check to see if next elmt is different or
5183  // collmax reached and if so end collection
5184  // and start new one
5185  if (prevCoeffOffset + nCoeffs != coeffOffset ||
5186  prevnCoeff != nCoeffs ||
5187  prevPhysOffset + nPhys != physOffset ||
5188  prevDeformed != Deformed || prevnPhys != nPhys ||
5189  collcnt >= collmax)
5190  {
5191 
5192  // if no Imp Type provided and No
5193  // settign in xml file. reset
5194  // impTypes using timings
5195  if (autotuning)
5196  {
5197  if (collExp.size() > collsize)
5198  {
5199  impTypes = colOpt.SetWithTimings(collExp, impTypes,
5200  verbose);
5201  collsize = collExp.size();
5202  }
5203  }
5204 
5205  Collections::Collection tmp(collExp, impTypes);
5206  m_collections.push_back(tmp);
5207 
5208  // start new geom list
5209  collExp.clear();
5210 
5211  m_coll_coeff_offset.push_back(coeffOffset);
5212  m_coll_phys_offset.push_back(physOffset);
5213  collExp.push_back(it.second[i].first);
5214  collcnt = 1;
5215  }
5216  else // add to list of collections
5217  {
5218  collExp.push_back(it.second[i].first);
5219  collcnt++;
5220  }
5221 
5222  // if end of list finish up collection
5223  if (i == it.second.size() - 1)
5224  {
5225  // if no Imp Type provided and No
5226  // settign in xml file.
5227  if (autotuning)
5228  {
5229  if (collExp.size() > collsize)
5230  {
5231  impTypes = colOpt.SetWithTimings(collExp, impTypes,
5232  verbose);
5233  collsize = collExp.size();
5234  }
5235  }
5236 
5237  Collections::Collection tmp(collExp, impTypes);
5238  m_collections.push_back(tmp);
5239 
5240  collExp.clear();
5241  collcnt = 0;
5242  }
5243 
5244  prevCoeffOffset = coeffOffset;
5245  prevPhysOffset = physOffset;
5246  prevDeformed = Deformed;
5247  prevnCoeff = nCoeffs;
5248  prevnPhys = nPhys;
5249  }
5250  }
5251  }
5252 
5253  // update optimisation file
5254  if ((m_session->GetUpdateOptFile()) && (ImpType == Collections::eNoImpType))
5255  {
5256  colOpt.UpdateOptFile(m_session->GetSessionName(), m_comm);
5257  // turn off writeoptfile option so only first
5258  // instance is timed
5259  m_session->SetUpdateOptFile(false);
5260  }
5261 }
5262 
5264 {
5266 }
5267 
5269  const Array<OneD, NekDouble> &inarray,
5270  Array<OneD, NekDouble> &outarray)
5271 {
5272  int cnt, cnt1;
5273 
5274  cnt = cnt1 = 0;
5275 
5276  switch (m_expType)
5277  {
5278  case e2D:
5279  {
5280  for (int i = 0; i < GetExpSize(); ++i)
5281  {
5282  // get new points key
5283  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5284  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5285  int npt0 = (int)pt0 * scale;
5286  int npt1 = (int)pt1 * scale;
5287 
5288  LibUtilities::PointsKey newPointsKey0(
5289  npt0, (*m_exp)[i]->GetPointsType(0));
5290  LibUtilities::PointsKey newPointsKey1(
5291  npt1, (*m_exp)[i]->GetPointsType(1));
5292 
5293  // Interpolate points;
5294  LibUtilities::Interp2D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5295  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5296  &inarray[cnt], newPointsKey0,
5297  newPointsKey1, &outarray[cnt1]);
5298 
5299  cnt += pt0 * pt1;
5300  cnt1 += npt0 * npt1;
5301  }
5302  }
5303  break;
5304  case e3D:
5305  {
5306  for (int i = 0; i < GetExpSize(); ++i)
5307  {
5308  // get new points key
5309  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5310  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5311  int pt2 = (*m_exp)[i]->GetNumPoints(2);
5312  int npt0 = (int)pt0 * scale;
5313  int npt1 = (int)pt1 * scale;
5314  int npt2 = (int)pt2 * scale;
5315 
5316  LibUtilities::PointsKey newPointsKey0(
5317  npt0, (*m_exp)[i]->GetPointsType(0));
5318  LibUtilities::PointsKey newPointsKey1(
5319  npt1, (*m_exp)[i]->GetPointsType(1));
5320  LibUtilities::PointsKey newPointsKey2(
5321  npt2, (*m_exp)[i]->GetPointsType(2));
5322 
5323  // Interpolate points;
5324  LibUtilities::Interp3D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5325  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5326  (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
5327  &inarray[cnt], newPointsKey0,
5328  newPointsKey1, newPointsKey2,
5329  &outarray[cnt1]);
5330 
5331  cnt += pt0 * pt1 * pt2;
5332  cnt1 += npt0 * npt1 * npt2;
5333  }
5334  }
5335  break;
5336  default:
5337  {
5338  NEKERROR(ErrorUtil::efatal, "This expansion is not set");
5339  }
5340  break;
5341  }
5342 }
5343 
5345  const Array<OneD, const NekDouble> &FwdFlux,
5346  const Array<OneD, const NekDouble> &BwdFlux,
5347  Array<OneD, NekDouble> &outarray)
5348 {
5349  boost::ignore_unused(FwdFlux, BwdFlux, outarray);
5350  NEKERROR(ErrorUtil::efatal, "AddTraceIntegralToOffDiag not defined");
5351 }
5352 
5354  const Array<OneD, const Array<OneD, NekDouble>> &inarray, const int nDirctn,
5355  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5356 {
5357  int nTotElmt = (*m_exp).size();
5358  int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5359  int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5360 
5361  Array<OneD, NekDouble> tmpCoef(nElmtCoef, 0.0);
5362  Array<OneD, NekDouble> tmpPhys(nElmtPnt, 0.0);
5363 
5364  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
5365  {
5366  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5367  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5368 
5369  if (tmpPhys.size() != nElmtPnt || tmpCoef.size() != nElmtCoef)
5370  {
5371  tmpPhys = Array<OneD, NekDouble>(nElmtPnt, 0.0);
5372  tmpCoef = Array<OneD, NekDouble>(nElmtCoef, 0.0);
5373  }
5374 
5375  for (int ncl = 0; ncl < nElmtPnt; ncl++)
5376  {
5377  tmpPhys[ncl] = inarray[nelmt][ncl];
5378 
5379  (*m_exp)[nelmt]->IProductWRTDerivBase(nDirctn, tmpPhys, tmpCoef);
5380 
5381  for (int nrw = 0; nrw < nElmtCoef; nrw++)
5382  {
5383  (*mtxPerVar[nelmt])(nrw, ncl) = tmpCoef[nrw];
5384  }
5385  // to maintain all the other columes are zero.
5386  tmpPhys[ncl] = 0.0;
5387  }
5388  }
5389 }
5390 
5392  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5393 {
5394  int nTotElmt = (*m_exp).size();
5395 
5396  int nspacedim = m_graph->GetSpaceDimension();
5397  Array<OneD, Array<OneD, NekDouble>> projectedpnts(nspacedim);
5398  Array<OneD, Array<OneD, NekDouble>> tmppnts(nspacedim);
5399  Array<OneD, DNekMatSharedPtr> ArrayStdMat(nspacedim);
5400  Array<OneD, Array<OneD, NekDouble>> ArrayStdMat_data(nspacedim);
5401 
5402  Array<OneD, NekDouble> clmnArray, clmnStdMatArray;
5403 
5405  int nElmtPntPrevious = 0;
5406  int nElmtCoefPrevious = 0;
5407 
5408  int nElmtPnt, nElmtCoef;
5409  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
5410  {
5411  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5412  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5413  LibUtilities::ShapeType ElmtTypeNow = (*m_exp)[nelmt]->DetShapeType();
5414 
5415  if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
5416  (ElmtTypeNow != ElmtTypePrevious))
5417  {
5418  if (nElmtPntPrevious != nElmtPnt)
5419  {
5420  for (int ndir = 0; ndir < nspacedim; ndir++)
5421  {
5422  projectedpnts[ndir] = Array<OneD, NekDouble>(nElmtPnt, 0.0);
5423  tmppnts[ndir] = Array<OneD, NekDouble>(nElmtPnt, 0.0);
5424  }
5425  }
5427  stdExp = (*m_exp)[nelmt]->GetStdExp();
5429  stdExp->DetShapeType(), *stdExp);
5430 
5431  ArrayStdMat[0] = stdExp->GetStdMatrix(matkey);
5432  ArrayStdMat_data[0] = ArrayStdMat[0]->GetPtr();
5433 
5434  if (nspacedim > 1)
5435  {
5436  StdRegions::StdMatrixKey matkey(
5437  StdRegions::eDerivBase1, stdExp->DetShapeType(), *stdExp);
5438 
5439  ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
5440  ArrayStdMat_data[1] = ArrayStdMat[1]->GetPtr();
5441 
5442  if (nspacedim > 2)
5443  {
5445  stdExp->DetShapeType(),
5446  *stdExp);
5447 
5448  ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
5449  ArrayStdMat_data[2] = ArrayStdMat[2]->GetPtr();
5450  }
5451  }
5452 
5453  ElmtTypePrevious = ElmtTypeNow;
5454  nElmtPntPrevious = nElmtPnt;
5455  nElmtCoefPrevious = nElmtCoef;
5456  }
5457  else
5458  {
5459  for (int ndir = 0; ndir < nspacedim; ndir++)
5460  {
5461  Vmath::Zero(nElmtPnt, projectedpnts[ndir], 1);
5462  }
5463  }
5464 
5465  for (int ndir = 0; ndir < nspacedim; ndir++)
5466  {
5467  (*m_exp)[nelmt]->AlignVectorToCollapsedDir(
5468  ndir, inarray[ndir][nelmt], tmppnts);
5469  for (int n = 0; n < nspacedim; n++)
5470  {
5471  Vmath::Vadd(nElmtPnt, tmppnts[n], 1, projectedpnts[n], 1,
5472  projectedpnts[n], 1);
5473  }
5474  }
5475 
5476  for (int ndir = 0; ndir < nspacedim; ndir++)
5477  {
5478  // weight with metric
5479  (*m_exp)[nelmt]->MultiplyByQuadratureMetric(projectedpnts[ndir],
5480  projectedpnts[ndir]);
5481  Array<OneD, NekDouble> MatDataArray = mtxPerVar[nelmt]->GetPtr();
5482 
5483  for (int np = 0; np < nElmtPnt; np++)
5484  {
5485  NekDouble factor = projectedpnts[ndir][np];
5486  clmnArray = MatDataArray + np * nElmtCoef;
5487  clmnStdMatArray = ArrayStdMat_data[ndir] + np * nElmtCoef;
5488  Vmath::Svtvp(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray,
5489  1, clmnArray, 1);
5490  }
5491  }
5492  }
5493 }
5494 
5495 // TODO: Reduce cost by getting Bwd Matrix directly
5497  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
5498  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5499 {
5501  int nElmtPntPrevious = 0;
5502  int nElmtCoefPrevious = 0;
5503  int nTotElmt = (*m_exp).size();
5504  int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5505  int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5506 
5507  Array<OneD, NekDouble> tmpPhys;
5508  Array<OneD, NekDouble> clmnArray, clmnStdMatArray;
5509  Array<OneD, NekDouble> stdMat_data;
5510 
5511  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
5512  {
5513  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5514  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5515  LibUtilities::ShapeType ElmtTypeNow = (*m_exp)[nelmt]->DetShapeType();
5516 
5517  if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
5518  (ElmtTypeNow != ElmtTypePrevious))
5519  {
5521  stdExp = (*m_exp)[nelmt]->GetStdExp();
5523  stdExp->DetShapeType(), *stdExp);
5524 
5525  DNekMatSharedPtr BwdMat = stdExp->GetStdMatrix(matkey);
5526  stdMat_data = BwdMat->GetPtr();
5527 
5528  if (nElmtPntPrevious != nElmtPnt)
5529  {
5530  tmpPhys = Array<OneD, NekDouble>(nElmtPnt, 0.0);
5531  }
5532 
5533  ElmtTypePrevious = ElmtTypeNow;
5534  nElmtPntPrevious = nElmtPnt;
5535  nElmtCoefPrevious = nElmtCoef;
5536  }
5537 
5538  (*m_exp)[nelmt]->MultiplyByQuadratureMetric(
5539  inarray[nelmt],
5540  tmpPhys); // weight with metric
5541 
5542  Array<OneD, NekDouble> MatDataArray = mtxPerVar[nelmt]->GetPtr();
5543 
5544  for (int np = 0; np < nElmtPnt; np++)
5545  {
5546  NekDouble factor = tmpPhys[np];
5547  clmnArray = MatDataArray + np * nElmtCoef;
5548  clmnStdMatArray = stdMat_data + np * nElmtCoef;
5549  Vmath::Smul(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray, 1);
5550  }
5551  }
5552 }
5553 
5554 /**
5555  * @brief inverse process of v_GetFwdBwdTracePhys. Given Trace integration of
5556  * Fwd and Bwd Jacobian, with dimension NtotalTrace*TraceCoef*TracePhys. return
5557  * Elemental Jacobian matrix with dimension
5558  * NtotalElement*ElementCoef*ElementPhys.
5559  *
5560  *
5561  * @param field is a NekDouble array which contains the 2D data
5562  * from which we wish to extract the backward and
5563  * forward orientated trace/edge arrays.
5564  * @param Fwd The resulting forwards space.
5565  * @param Bwd The resulting backwards space.
5566  */
5571 {
5573  std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
5574  tracelist->GetExp();
5575  int ntotTrace = (*traceExp).size();
5576  int nTracePnt, nTraceCoef;
5577 
5578  std::shared_ptr<LocalRegions::ExpansionVector> fieldExp = GetExp();
5579  int nElmtCoef;
5580 
5581  const MultiRegions::LocTraceToTraceMapSharedPtr locTraceToTraceMap =
5583  const Array<OneD, const Array<OneD, int>> LRAdjExpid =
5584  locTraceToTraceMap->GetLeftRightAdjacentExpId();
5585  const Array<OneD, const Array<OneD, bool>> LRAdjflag =
5586  locTraceToTraceMap->GetLeftRightAdjacentExpFlag();
5587 
5589  locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffMap();
5591  locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffSign();
5592  DNekMatSharedPtr ElmtMat;
5593  Array<OneD, NekDouble> ElmtMat_data;
5594  // int nclAdjExp;
5595  int nrwAdjExp;
5596  int MatIndex, nPnts;
5597  NekDouble sign = 1.0;
5598 
5599  int nTracePntsTtl = tracelist->GetTotPoints();
5600  int nlocTracePts = locTraceToTraceMap->GetNLocTracePts();
5601  int nlocTracePtsFwd = locTraceToTraceMap->GetNFwdLocTracePts();
5602  int nlocTracePtsBwd = nlocTracePts - nlocTracePtsFwd;
5603 
5604  Array<OneD, int> nlocTracePtsLR(2);
5605  nlocTracePtsLR[0] = nlocTracePtsFwd;
5606  nlocTracePtsLR[1] = nlocTracePtsBwd;
5607 
5608  size_t nFwdBwdNonZero = 0;
5609  Array<OneD, int> tmpIndex{2, -1};
5610  for (int i = 0; i < 2; ++i)
5611  {
5612  if (nlocTracePtsLR[i] > 0)
5613  {
5614  tmpIndex[nFwdBwdNonZero] = i;
5615  nFwdBwdNonZero++;
5616  }
5617  }
5618 
5619  Array<OneD, int> nlocTracePtsNonZeroIndex{nFwdBwdNonZero};
5620  for (int i = 0; i < nFwdBwdNonZero; ++i)
5621  {
5622  nlocTracePtsNonZeroIndex[i] = tmpIndex[i];
5623  }
5624 
5625  Array<OneD, NekDouble> TraceFwdPhy(nTracePntsTtl);
5626  Array<OneD, NekDouble> TraceBwdPhy(nTracePntsTtl);
5627  Array<OneD, Array<OneD, NekDouble>> tmplocTrace(2);
5628  for (int k = 0; k < 2; ++k)
5629  {
5630  tmplocTrace[k] = NullNekDouble1DArray;
5631  }
5632 
5633  for (int k = 0; k < nFwdBwdNonZero; ++k)
5634  {
5635  size_t i = nlocTracePtsNonZeroIndex[k];
5636  tmplocTrace[i] = Array<OneD, NekDouble>(nlocTracePtsLR[i]);
5637  }
5638 
5639  int nNumbElmt = fieldMat.size();
5640  Array<OneD, Array<OneD, NekDouble>> ElmtMatDataArray(nNumbElmt);
5641  Array<OneD, int> ElmtCoefArray(nNumbElmt);
5642  for (int i = 0; i < nNumbElmt; i++)
5643  {
5644  ElmtMatDataArray[i] = fieldMat[i]->GetPtr();
5645  ElmtCoefArray[i] = GetNcoeffs(i);
5646  }
5647 
5648  int nTraceCoefMax = 0;
5649  int nTraceCoefMin = std::numeric_limits<int>::max();
5650  Array<OneD, int> TraceCoefArray(ntotTrace);
5651  Array<OneD, int> TracePntArray(ntotTrace);
5652  Array<OneD, int> TraceOffArray(ntotTrace);
5653  Array<OneD, Array<OneD, NekDouble>> FwdMatData(ntotTrace);
5654  Array<OneD, Array<OneD, NekDouble>> BwdMatData(ntotTrace);
5655  for (int nt = 0; nt < ntotTrace; nt++)
5656  {
5657  nTraceCoef = (*traceExp)[nt]->GetNcoeffs();
5658  nTracePnt = tracelist->GetTotPoints(nt);
5659  int noffset = tracelist->GetPhys_Offset(nt);
5660  TraceCoefArray[nt] = nTraceCoef;
5661  TracePntArray[nt] = nTracePnt;
5662  TraceOffArray[nt] = noffset;
5663  FwdMatData[nt] = FwdMat[nt]->GetPtr();
5664  BwdMatData[nt] = BwdMat[nt]->GetPtr();
5665  if (nTraceCoef > nTraceCoefMax)
5666  {
5667  nTraceCoefMax = nTraceCoef;
5668  }
5669  if (nTraceCoef < nTraceCoefMin)
5670  {
5671  nTraceCoefMin = nTraceCoef;
5672  }
5673  }
5674  WARNINGL1(nTraceCoefMax == nTraceCoefMin,
5675  "nTraceCoefMax!=nTraceCoefMin: Effeciency may be low ");
5676 
5677  int traceID, nfieldPnts, ElmtId, noffset;
5678  const Array<OneD, const Array<OneD, int>> LocTracephysToTraceIDMap =
5679  locTraceToTraceMap->GetLocTracephysToTraceIDMap();
5680  const Array<OneD, const int> fieldToLocTraceMap =
5681  locTraceToTraceMap->GetLocTraceToFieldMap();
5682  Array<OneD, Array<OneD, int>> fieldToLocTraceMapLR(2);
5683  noffset = 0;
5684  for (int k = 0; k < nFwdBwdNonZero; ++k)
5685  {
5686  size_t i = nlocTracePtsNonZeroIndex[k];
5687  fieldToLocTraceMapLR[i] = Array<OneD, int>(nlocTracePtsLR[i]);
5688  Vmath::Vcopy(nlocTracePtsLR[i], &fieldToLocTraceMap[0] + noffset, 1,
5689  &fieldToLocTraceMapLR[i][0], 1);
5690  noffset += nlocTracePtsLR[i];
5691  }
5692 
5693  Array<OneD, Array<OneD, int>> MatIndexArray(2);
5694  for (int k = 0; k < nFwdBwdNonZero; ++k)
5695  {
5696  size_t nlr = nlocTracePtsNonZeroIndex[k];
5697  MatIndexArray[nlr] = Array<OneD, int>(nlocTracePtsLR[nlr]);
5698  for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5699  {
5700  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5701  nTraceCoef = TraceCoefArray[traceID];
5702  ElmtId = LRAdjExpid[nlr][traceID];
5703  noffset = GetPhys_Offset(ElmtId);
5704  nElmtCoef = ElmtCoefArray[ElmtId];
5705  nfieldPnts = fieldToLocTraceMapLR[nlr][nloc];
5706  nPnts = nfieldPnts - noffset;
5707 
5708  MatIndexArray[nlr][nloc] = nPnts * nElmtCoef;
5709  }
5710  }
5711 
5712  for (int nc = 0; nc < nTraceCoefMin; nc++)
5713  {
5714  for (int nt = 0; nt < ntotTrace; nt++)
5715  {
5716  nTraceCoef = TraceCoefArray[nt];
5717  nTracePnt = TracePntArray[nt];
5718  noffset = TraceOffArray[nt];
5719  Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
5720  &TraceFwdPhy[noffset], 1);
5721  Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
5722  &TraceBwdPhy[noffset], 1);
5723  }
5724 
5725  for (int k = 0; k < nFwdBwdNonZero; ++k)
5726  {
5727  size_t i = nlocTracePtsNonZeroIndex[k];
5728  Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
5729  }
5730 
5731  GetLocTraceFromTracePts(TraceFwdPhy, TraceBwdPhy, tmplocTrace[0],
5732  tmplocTrace[1]);
5733 
5734  for (int k = 0; k < nFwdBwdNonZero; ++k)
5735  {
5736  size_t nlr = nlocTracePtsNonZeroIndex[k];
5737  for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5738  {
5739  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5740  nTraceCoef = TraceCoefArray[traceID];
5741  ElmtId = LRAdjExpid[nlr][traceID];
5742  nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5743  sign = elmtLRSign[nlr][traceID][nc];
5744  MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5745 
5746  ElmtMatDataArray[ElmtId][MatIndex] -=
5747  sign * tmplocTrace[nlr][nloc];
5748  }
5749  }
5750  }
5751 
5752  for (int nc = nTraceCoefMin; nc < nTraceCoefMax; nc++)
5753  {
5754  for (int nt = 0; nt < ntotTrace; nt++)
5755  {
5756  nTraceCoef = TraceCoefArray[nt];
5757  nTracePnt = TracePntArray[nt];
5758  noffset = TraceOffArray[nt];
5759  if (nc < nTraceCoef)
5760  {
5761  Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
5762  &TraceFwdPhy[noffset], 1);
5763  Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
5764  &TraceBwdPhy[noffset], 1);
5765  }
5766  else
5767  {
5768  Vmath::Zero(nTracePnt, &TraceFwdPhy[noffset], 1);
5769  Vmath::Zero(nTracePnt, &TraceBwdPhy[noffset], 1);
5770  }
5771  }
5772 
5773  for (int k = 0; k < nFwdBwdNonZero; ++k)
5774  {
5775  size_t i = nlocTracePtsNonZeroIndex[k];
5776  Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
5777  }
5778  GetLocTraceFromTracePts(TraceFwdPhy, TraceBwdPhy, tmplocTrace[0],
5779  tmplocTrace[1]);
5780 
5781  for (int k = 0; k < nFwdBwdNonZero; ++k)
5782  {
5783  size_t nlr = nlocTracePtsNonZeroIndex[k];
5784  for (int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5785  {
5786  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5787  nTraceCoef = TraceCoefArray[traceID];
5788  if (nc < nTraceCoef)
5789  {
5790  ElmtId = LRAdjExpid[nlr][traceID];
5791  nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5792  sign = -elmtLRSign[nlr][traceID][nc];
5793  MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5794 
5795  ElmtMatDataArray[ElmtId][MatIndex] +=
5796  sign * tmplocTrace[nlr][nloc];
5797  }
5798  }
5799  }
5800  }
5801 }
5802 
5804  const int dir, const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
5805  Array<OneD, DNekMatSharedPtr> ElmtJacCoef)
5806 {
5807  int nelmt;
5808  int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
5809 
5810  nelmtcoef = GetNcoeffs(0);
5811  nelmtpnts = GetTotPoints(0);
5812 
5813  Array<OneD, NekDouble> innarray(nelmtpnts, 0.0);
5814  Array<OneD, NekDouble> outarray(nelmtcoef, 0.0);
5815 
5816  Array<OneD, NekDouble> MatQ_data;
5817  Array<OneD, NekDouble> MatC_data;
5818 
5819  DNekMatSharedPtr tmpMatQ, tmpMatC;
5820 
5821  nelmtcoef0 = nelmtcoef;
5822  nelmtpnts0 = nelmtpnts;
5823 
5824  for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5825  {
5826  nelmtcoef = GetNcoeffs(nelmt);
5827  nelmtpnts = GetTotPoints(nelmt);
5828 
5829  tmpMatQ = ElmtJacQuad[nelmt];
5830  tmpMatC = ElmtJacCoef[nelmt];
5831 
5832  MatQ_data = tmpMatQ->GetPtr();
5833  MatC_data = tmpMatC->GetPtr();
5834 
5835  if (nelmtcoef != nelmtcoef0)
5836  {
5837  outarray = Array<OneD, NekDouble>(nelmtcoef, 0.0);
5838  nelmtcoef0 = nelmtcoef;
5839  }
5840 
5841  if (nelmtpnts != nelmtpnts0)
5842  {
5843  innarray = Array<OneD, NekDouble>(nelmtpnts, 0.0);
5844  nelmtpnts0 = nelmtpnts;
5845  }
5846 
5847  for (int np = 0; np < nelmtcoef; np++)
5848  {
5849  Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
5850  1);
5851  (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
5852  (*m_exp)[nelmt]->IProductWRTDerivBase(dir, innarray, outarray);
5853 
5854  Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
5855  nelmtcoef, &MatC_data[0] + np, nelmtcoef);
5856  }
5857  }
5858 }
5859 
5861  const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
5862  Array<OneD, DNekMatSharedPtr> ElmtJacCoef)
5863 {
5864  int nelmt;
5865  int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
5866 
5867  nelmtcoef = GetNcoeffs(0);
5868  nelmtpnts = GetTotPoints(0);
5869 
5870  Array<OneD, NekDouble> innarray(nelmtpnts, 0.0);
5871  Array<OneD, NekDouble> outarray(nelmtcoef, 0.0);
5872 
5873  Array<OneD, NekDouble> MatQ_data;
5874  Array<OneD, NekDouble> MatC_data;
5875 
5876  DNekMatSharedPtr tmpMatQ, tmpMatC;
5877 
5878  nelmtcoef0 = nelmtcoef;
5879  nelmtpnts0 = nelmtpnts;
5880 
5881  for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5882  {
5883  nelmtcoef = GetNcoeffs(nelmt);
5884  nelmtpnts = GetTotPoints(nelmt);
5885 
5886  tmpMatQ = ElmtJacQuad[nelmt];
5887  tmpMatC = ElmtJacCoef[nelmt];
5888 
5889  MatQ_data = tmpMatQ->GetPtr();
5890  MatC_data = tmpMatC->GetPtr();
5891 
5892  if (nelmtcoef != nelmtcoef0)
5893  {
5894  outarray = Array<OneD, NekDouble>(nelmtcoef, 0.0);
5895  nelmtcoef0 = nelmtcoef;
5896  }
5897 
5898  if (nelmtpnts != nelmtpnts0)
5899  {
5900  innarray = Array<OneD, NekDouble>(nelmtpnts, 0.0);
5901  nelmtpnts0 = nelmtpnts;
5902  }
5903 
5904  for (int np = 0; np < nelmtcoef; np++)
5905  {
5906  Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
5907  1);
5908  (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
5909  (*m_exp)[nelmt]->IProductWRTBase(innarray, outarray);
5910 
5911  Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
5912  nelmtcoef, &MatC_data[0] + np, nelmtcoef);
5913  }
5914  }
5915 }
5916 
5918  const NekDouble scale, const Array<OneD, NekDouble> &inarray,
5919  Array<OneD, NekDouble> &outarray)
5920 {
5921  int cnt, cnt1;
5922 
5923  cnt = cnt1 = 0;
5924 
5925  switch (m_expType)
5926  {
5927  case e2D:
5928  {
5929  for (int i = 0; i < GetExpSize(); ++i)
5930  {
5931  // get new points key
5932  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5933  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5934  int npt0 = (int)pt0 * scale;
5935  int npt1 = (int)pt1 * scale;
5936 
5937  LibUtilities::PointsKey newPointsKey0(
5938  npt0, (*m_exp)[i]->GetPointsType(0));
5939  LibUtilities::PointsKey newPointsKey1(
5940  npt1, (*m_exp)[i]->GetPointsType(1));
5941 
5942  // Project points;
5944  newPointsKey0, newPointsKey1, &inarray[cnt],
5945  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5946  (*m_exp)[i]->GetBasis(1)->GetPointsKey(), &outarray[cnt1]);
5947 
5948  cnt += npt0 * npt1;
5949  cnt1 += pt0 * pt1;
5950  }
5951  }
5952  break;
5953  case e3D:
5954  {
5955  for (int i = 0; i < GetExpSize(); ++i)
5956  {
5957  // get new points key
5958  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5959  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5960  int pt2 = (*m_exp)[i]->GetNumPoints(2);
5961  int npt0 = (int)pt0 * scale;
5962  int npt1 = (int)pt1 * scale;
5963  int npt2 = (int)pt2 * scale;
5964 
5965  LibUtilities::PointsKey newPointsKey0(
5966  npt0, (*m_exp)[i]->GetPointsType(0));
5967  LibUtilities::PointsKey newPointsKey1(
5968  npt1, (*m_exp)[i]->GetPointsType(1));
5969  LibUtilities::PointsKey newPointsKey2(
5970  npt2, (*m_exp)[i]->GetPointsType(2));
5971 
5972  // Project points;
5974  newPointsKey0, newPointsKey1, newPointsKey2, &inarray[cnt],
5975  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5976  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5977  (*m_exp)[i]->GetBasis(2)->GetPointsKey(), &outarray[cnt1]);
5978 
5979  cnt += npt0 * npt1 * npt2;
5980  cnt1 += pt0 * pt1 * pt2;
5981  }
5982  }
5983  break;
5984  default:
5985  {
5986  NEKERROR(ErrorUtil::efatal, "not setup for this expansion");
5987  }
5988  break;
5989  }
5990 }
5991 
5993 {
5994  NEKERROR(ErrorUtil::efatal, "v_GetLocTraceToTraceMap not coded");
5996 }
5997 
5998 } // namespace MultiRegions
5999 } // namespace Nektar
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:250
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define LIKWID_MARKER_START(regionTag)
Definition: Likwid.hpp:46
#define LIKWID_MARKER_STOP(regionTag)
Definition: Likwid.hpp:47
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
COLLECTIONS_EXPORT OperatorImpMap SetWithTimings(std::vector< StdRegions::StdExpansionSharedPtr > pGeom, OperatorImpMap &impTypes, bool verbose=true)
COLLECTIONS_EXPORT void UpdateOptFile(std::string sessName, LibUtilities::CommSharedPtr &comm)
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Describes the specification for a Basis.
Definition: Basis.h:50
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:130
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:141
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:147
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:83
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
Definition: FieldIO.cpp:97
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Defines a specification for a set of points.
Definition: Points.h:59
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:73
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2670
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1158
void AddTraceJacToElmtJac(const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian,...
Definition: ExpList.cpp:5567
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1787
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions()
Definition: ExpList.h:2270
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2319
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1641
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Definition: ExpList.cpp:3649
virtual void v_Upwind(const Array< OneD, const Array< OneD, NekDouble >> &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.cpp:3792
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2010
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4371
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:104
void ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract the data in fielddata into the coeffs.
Definition: ExpList.cpp:3538
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:5071
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.cpp:4961
void ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: ExpList.cpp:1743
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:4527
std::shared_ptr< GlobalMatrix > GenGlobalMatrix(const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Generates a global matrix from the given key and map.
Definition: ExpList.cpp:2126
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:4518
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4686
void AddRightIPTPhysDerivBase(const int dir, const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:5803
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4459
int GetExpIndex(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
This function returns the index of the local elemental expansion containing the arbitrary point given...
Definition: ExpList.cpp:2518
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2828
BlockMatrixMapShPtr m_blockMat
Definition: ExpList.h:1218
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
Definition: ExpList.cpp:3290
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
multiply the metric jacobi and quadrature weights
Definition: ExpList.cpp:1340
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2450
void WriteVtkPieceFooter(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:3089
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
Definition: ExpList.cpp:4880
void GetBwdWeight(Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
Get the weight value for boundary conditions for boundary average and jump calculations.
Definition: ExpList.cpp:3678
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5344
virtual std::vector< bool > & v_GetLeftAdjacentTraces(void)
Definition: ExpList.cpp:3911
void InitialiseExpVector(const SpatialDomains::ExpansionInfoMap &expmap)
Define a list of elements using the geometry and basis key information in expmap;.
Definition: ExpList.cpp:1135
void SetupCoeffPhys(bool DeclareCoeffPhysArrays=true, bool SetupOffsets=true)
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
Definition: ExpList.cpp:1092
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:5038
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble >> &inarray1, const Array< OneD, Array< OneD, NekDouble >> &inarray2, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: ExpList.cpp:4545
virtual void v_CurlCurl(Array< OneD, Array< OneD, NekDouble >> &Vel, Array< OneD, Array< OneD, NekDouble >> &Q)
Definition: ExpList.cpp:1645
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1713
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
Definition: ExpList.cpp:4555
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.cpp:4426
virtual void v_WriteTecplotHeader(std::ostream &outfile, std::string var="")
Definition: ExpList.cpp:2697
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
Definition: ExpList.cpp:5059
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n.
Definition: ExpList.h:2232
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2204
virtual void v_AddTraceQuadPhysToOffDiag(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:4417
virtual void v_ClearGlobalLinSysManager(void)
Definition: ExpList.cpp:3298
void GeneralGetFieldDefinitions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, int NumHomoDir=0, Array< OneD, LibUtilities::BasisSharedPtr > &HomoBasis=LibUtilities::NullBasisSharedPtr1DArray, std::vector< NekDouble > &HomoLen=LibUtilities::NullNekDoubleVector, bool homoStrips=false, std::vector< unsigned int > &HomoSIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoZIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoYIDs=LibUtilities::NullUnsignedIntVector)
Definition: ExpList.cpp:3379
virtual void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
Definition: ExpList.cpp:5029
void IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to the derivative (in directio...
Definition: ExpList.cpp:1384
virtual void v_LocalToGlobal(bool UseComm)
Definition: ExpList.cpp:4640
int GetNumElmts(void)
This function returns the number of elements in the expansion which may be different for a homogeoeno...
Definition: ExpList.h:694
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions()
Definition: ExpList.cpp:4986
virtual ~ExpList()
The default destructor.
Definition: ExpList.cpp:1308
virtual void v_FillBndCondFromField()
Definition: ExpList.cpp:4625
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1802
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4616
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
Definition: ExpList.cpp:5007
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
Definition: ExpList.cpp:4451
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble >> &normals)
Populate normals with the normals of all expansions.
Definition: ExpList.cpp:3989
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1210
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble >> &inarray)
Definition: ExpList.cpp:3231
std::shared_ptr< AssemblyMapDG > & GetTraceMap(void)
Definition: ExpList.h:2324
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap()
Definition: ExpList.cpp:3898
Array< OneD, std::pair< int, int > > m_coeffsToElmt
m_coeffs to elemental value map
Definition: ExpList.h:1216
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:4380
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
Definition: ExpList.cpp:2635
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Definition: ExpList.cpp:3217
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo(void)
Definition: ExpList.cpp:5019
virtual const Array< OneD, const int > & v_GetTraceBndMap()
Definition: ExpList.cpp:3906
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
Definition: ExpList.cpp:3251
virtual NekDouble v_GetHomoLen(void)
Definition: ExpList.cpp:3267
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:3304
const LocTraceToTraceMapSharedPtr & GetLocTraceToTraceMap() const
Definition: ExpList.h:2283
virtual void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:4507
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
Definition: ExpList.cpp:3282
void ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Extract the data from fromField using fromExpList the coeffs using the basic ExpList Elemental expans...
Definition: ExpList.cpp:3546
void GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
This function calculates the coordinates of all the elemental quadrature points .
Definition: ExpList.h:1907
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:2431
virtual void v_SetBndCondBwdWeight(const int index, const NekDouble value)
Definition: ExpList.cpp:4444
void GetMatIpwrtDeriveBase(const Array< OneD, const Array< OneD, NekDouble >> &inarray, const int nDirctn, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:5353
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:1834
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
Definition: ExpList.cpp:3506
std::vector< bool > m_collectionsDoInit
Vector of bools to act as an initialise on first call flag.
Definition: ExpList.h:1201
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1184
void AddRightIPTBaseMatrix(const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:5860
void GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.h:2388
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
Definition: ExpList.cpp:3096
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Definition: ExpList.cpp:2955
int GetShapeDimension()
This function returns the dimension of the shape of the element eid.
Definition: ExpList.h:2064
NekDouble H1(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:3361
std::vector< int > m_coll_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1207
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1126
std::shared_ptr< GlobalLinSys > GenGlobalLinSys(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
This operation constructs the global linear system of type mkey.
Definition: ExpList.cpp:2413
void WriteVtkHeader(std::ostream &outfile)
Definition: ExpList.cpp:2941
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1196
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
Definition: ExpList.cpp:4720
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1136
void GenerateElementVector(const int ElementID, const NekDouble scalar1, const NekDouble scalar2, Array< OneD, NekDouble > &outarray)
Generate vector v such that v[i] = scalar1 if i is in the element < ElementID. Otherwise,...
Definition: ExpList.cpp:3735
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Definition: ExpList.cpp:2260
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
Definition: ExpList.cpp:4563
virtual void v_ExtractElmtToBndPhys(const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
Definition: ExpList.cpp:4803
virtual void v_ExtractPhysToBndElmt(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
Definition: ExpList.cpp:4841
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4354
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Definition: ExpList.cpp:4969
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5917
void GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the result of the multiplication of a matrix of type specified by mkey with ...
Definition: ExpList.cpp:2054
virtual void v_SetUpPhysNormals()
: Set up a normal along the trace elements between two elements at elemental level
Definition: ExpList.cpp:4780
virtual void v_GlobalToLocal(void)
Definition: ExpList.cpp:4655
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
Definition: ExpList.cpp:3259
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2223
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Definition: ExpList.cpp:4483
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1132
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4474
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:4408
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight()
Definition: ExpList.cpp:4436
virtual void v_GetMovingFrames(const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: ExpList.cpp:3691
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Definition: ExpList.cpp:3172
void WriteVtkFooter(std::ostream &outfile)
Definition: ExpList.cpp:2949
void MultiplyByBlockMatrix(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1320
virtual std::shared_ptr< ExpList > & v_GetTrace()
Definition: ExpList.cpp:3890
void GetElmtNormalLength(Array< OneD, NekDouble > &lengthsFwd, Array< OneD, NekDouble > &lengthsBwd)
Get the length of elements in boundary normal direction.
Definition: ExpList.cpp:4206
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
Definition: ExpList.cpp:3491
std::vector< int > m_coll_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1204
void IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directional derivative along a given direction.
Definition: ExpList.cpp:1404
void FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill Bwd with boundary conditions.
Definition: ExpList.h:2396
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2467
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
Definition: ExpList.cpp:4997
virtual void v_LinearAdvectionDiffusionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:4496
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2908
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: ExpList.h:1034
Collections::CollectionVector m_collections
Definition: ExpList.h:1198
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5268
virtual void v_GetBoundaryNormals(int i, Array< OneD, Array< OneD, NekDouble >> &normals)
Definition: ExpList.cpp:4917
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2658
Array< OneD, const unsigned int > GetZIDs(void)
This function returns a vector containing the wave numbers in z-direction associated with the 3D homo...
Definition: ExpList.h:631
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: ExpList.cpp:1548
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1811
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1227
void GetDiagMatIpwrtBase(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:5496
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2735
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:4976
std::map< int, RobinBCInfoSharedPtr > GetRobinBCInfo()
Definition: ExpList.h:963
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
Definition: ExpList.cpp:4399
void DivideByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Divided by the metric jacobi and quadrature weights.
Definition: ExpList.cpp:1357
virtual const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions(void)
Definition: ExpList.cpp:3765
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1213
virtual void v_SetHomoLen(const NekDouble lhom)
Definition: ExpList.cpp:3275
void Upwind(const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.h:2311
const DNekScalBlkMatSharedPtr GenBlockMatrix(const GlobalMatrixKey &gkey)
This function assembles the block diagonal matrix of local matrices of the type mtype.
Definition: ExpList.cpp:1874
void PhysDeriv(Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
Definition: ExpList.h:1954
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally mulplies the coefficient space of Sin my the elemental inverse of the mass ...
Definition: ExpList.cpp:1764
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
Definition: ExpList.h:1849
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Definition: ExpList.cpp:4793
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4669
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
Definition: ExpList.cpp:3773
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n.
Definition: ExpList.h:2240
virtual const std::shared_ptr< LocTraceToTraceMap > & v_GetLocTraceToTraceMap(void) const
Definition: ExpList.cpp:5992
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
Definition: ExpList.cpp:3562
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1119
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1175
ExpansionType GetExpType(void)
Returns the type of the expansion.
Definition: ExpList.cpp:1303
NekDouble Linf(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:3134
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4536
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1688
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:2054
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
LibUtilities::ShapeType GetShapeType() const
Return the expansion type associated with key.
const StdRegions::VarCoeffMap & GetVarCoeffs() const
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:112
void PhysGalerkinProject3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:52
std::shared_ptr< SessionReader > SessionReaderSharedPtr
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
Definition: Interp.cpp:167
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:177
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:106
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset=0)
Definition: ShapeType.hpp:308
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:72
std::shared_ptr< Transposition > TranspositionSharedPtr
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:83
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:59
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
std::shared_ptr< PointExp > PointExpSharedPtr
Definition: PointExp.h:135
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:48
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:50
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:89
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: ExpList.cpp:3922
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1633
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:51
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
const char *const GlobalSysSolnTypeMap[]
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
Definition: GlobalMatrix.h:86
GlobalLinSysFactory & GetGlobalLinSysFactory()
static LocTraceToTraceMapSharedPtr NullLocTraceToTraceMapSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:51
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< GlobalMatrixKey, DNekScalBlkMatSharedPtr > BlockMatrixMap
A map between global matrix keys and their associated block matrices.
Definition: ExpList.h:93
static const NekDouble kNekZeroTol
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
GeomMMF
Principle direction for MMF.
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:212
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:84
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:138
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:219
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:77
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:84
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:218
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:86
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:63
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:143
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:68
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::pair< IndexType, IndexType > CoordType
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
std::map< CoordType, NekDouble > COOMatType
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:822
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:945
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:71
std::shared_ptr< RobinBCInfo > next