Nektar++
ExpList1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList1D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Expansion list 1D definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 #include <MultiRegions/ExpList1D.h>
39 #include <LocalRegions/SegExp.h>
40 #include <LocalRegions/QuadExp.h>
43 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
45 
46 namespace Nektar
47 {
48  namespace MultiRegions
49  {
50  /**
51  * @class ExpList1D
52  * This multi-elemental expansion, which does not exhibit any coupling
53  * between the expansion on the separate elements, can be formulated
54  * as, \f[u^{\delta}(x_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
55  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(x_i).\f]
56  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
57  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
58  *
59  * Instances of this class may be optionally constructed which use
60  * generalised segment expansions (LocalRegions#GenSegExp), rather than
61  * the standard segment expansions (LocalRegions#SegExp).
62  * LocalRegions#GenSegExp provides additional spatial
63  * data including segment normals and is enabled using the \a
64  * UseGenSegExp flag.
65  *
66  * This class inherits all its variables and member functions from the
67  * base class MultiRegions#ExpList.
68  */
69 
70  /**
71  * Assumes use of standard segment expansions only. All data storage
72  * areas are initialised to empty arrays by the default ExpList
73  * constructor.
74  */
76  ExpList()
77  {
78  SetExpType(e1D);
79  }
80 
81 
82  /**
83  * Creates an identical copy of another ExpList1D object.
84  */
85  ExpList1D::ExpList1D(const ExpList1D &In, const bool DeclareCoeffPhysArrays):
86  ExpList(In,DeclareCoeffPhysArrays)
87  {
88  SetExpType(e1D);
89  }
90 
91 
92  /**
93  * After initialising the data inherited through MultiRegions#ExpList,
94  * populate the expansion list from the segments defined in the supplied
95  * SpatialDomains#MeshGraph1D. All expansions in the graph are defined
96  * using the same LibUtilities#BasisKey which overrides that specified
97  * in \a graph1D.
98  *
99  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
100  * for details.
101  * @deprecated This constructor is no longer required since
102  * the basis key is now available from the graph.
103  * @param Ba BasisKey describing quadrature points and
104  * number of modes.
105  * @param graph1D Domain and expansion definitions.
106  */
108  const LibUtilities::BasisKey &Ba,
109  const SpatialDomains::MeshGraphSharedPtr &graph1D):
110  ExpList(pSession,graph1D)
111  {
112  SetExpType(e1D);
113 
114  int id=0;
117 
118  const SpatialDomains::ExpansionMap &expansions
119  = graph1D->GetExpansions();
120 
121  // For each element in the mesh, create a segment expansion using
122  // the supplied BasisKey and segment geometry.
123  SpatialDomains::ExpansionMap::const_iterator expIt;
124  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
125  {
126  if ((SegmentGeom = boost::dynamic_pointer_cast<
128  expIt->second->m_geomShPtr)))
129  {
131  ::AllocateSharedPtr(Ba,SegmentGeom);
132  seg->SetElmtId(id++);
133  (*m_exp).push_back(seg);
134  }
135  else
136  {
137  ASSERTL0(false,"dynamic cast to a SegGeom failed");
138  }
139  }
140 
141  // Setup Default optimisation information.
142  int nel = GetExpSize();
145 
146  // Allocate storage for data and populate element offset lists.
148 
151 
154  }
155 
156 
157  /**
158  * Given a mesh \a graph1D, containing information about the domain and
159  * the spectral/hp element expansion, this constructor fills the list
160  * of local expansions \texttt{m_exp} with the proper expansions,
161  * calculates the total number of quadrature points \f$x_i\f$ and local
162  * expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory for
163  * the arrays #m_coeffs and #m_phys.
164  *
165  * For each element its corresponding LibUtilities#BasisKey is
166  * retrieved and this is used to construct either a standard segment
167  * (LocalRegions#SegExp) or a generalised segment
168  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
169  * Finally, ExpList#SetCoeffPhys is called to initialise the data
170  * storage areas and set up the offset arrays.
171  *
172  * @param graph1D A mesh, containing information about the
173  * domain and the spectral/hp element expansion.
174  * @param UseGenSegExp If true, create general segment expansions
175  * instead of just normal segment expansions.
176  */
178  const SpatialDomains::MeshGraphSharedPtr &graph1D,
179  const bool DeclareCoeffPhysArrays):
180  ExpList(pSession,graph1D)
181  {
182  SetExpType(e1D);
183 
184  int id=0;
187 
188  // Retrieve the list of expansions
189  const SpatialDomains::ExpansionMap &expansions
190  = graph1D->GetExpansions();
191 
192  // Process each expansion in the graph
193  SpatialDomains::ExpansionMap::const_iterator expIt;
194  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
195  {
196  // Retrieve basis key from expansion
197  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
198 
199  if ((SegmentGeom = boost::dynamic_pointer_cast<
201  expIt->second->m_geomShPtr)))
202  {
204  ::AllocateSharedPtr(bkey, SegmentGeom);
205 
206  // Assign next ID
207  seg->SetElmtId(id++);
208 
209  // Add the expansion
210  (*m_exp).push_back(seg);
211  }
212  else
213  {
214  ASSERTL0(false,"dynamic cast to a SegGeom failed");
215  }
216  }
217 
218  // Setup Default optimisation information.
219  int nel = GetExpSize();
222 
223  // set up offset arrays.
225 
226  if(DeclareCoeffPhysArrays)
227  {
228  // Set up m_coeffs, m_phys.
231  }
232 
235  }
236 
237 
238 
239  /**
240  * Given a mesh \a graph1D, and the spectral/hp element
241  * expansion as well as and separate information about a \a
242  * domain, this constructor fills the list of local
243  * expansions \texttt{m_exp} with the proper expansions,
244  * calculates the total number of quadrature points \f$x_i\f$
245  * and local expansion coefficients \f$\hat{u}^e_n\f$ and
246  * allocates memory for the arrays #m_coeffs and #m_phys.
247  *
248  * For each element its corresponding LibUtilities#BasisKey is
249  * retrieved and this is used to construct either a standard segment
250  * (LocalRegions#SegExp) or a generalised segment
251  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
252  * Finally, ExpList#SetCoeffPhys is called to initialise the data
253  * storage areas and set up the offset arrays.
254  *
255  * @param graph1D A mesh, containing information about the
256  * domain and the spectral/hp element expansion.
257  * @param UseGenSegExp If true, create general segment expansions
258  * instead of just normal segment expansions.
259  */
261  const SpatialDomains::MeshGraphSharedPtr &graph1D,
262  const SpatialDomains::CompositeMap &domain,
263  const bool DeclareCoeffPhysArrays,
264  const std::string var,
265  bool SetToOneSpaceDimension):
266  ExpList(pSession,graph1D)
267  {
268  int j,id=0;
272  SpatialDomains::CompositeMap::const_iterator compIt;
273 
274  // Retrieve the list of expansions
275  const SpatialDomains::ExpansionMap &expansions
276  = graph1D->GetExpansions(var);
277 
278  // Process each composite region in domain
279  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
280  {
281  comp = compIt->second;
282 
283  // Process each expansion in the graph
284  for(j = 0; j < compIt->second->size(); ++j)
285  {
286  SpatialDomains::ExpansionMap::const_iterator expIt;
287 
288  if((SegmentGeom = boost::dynamic_pointer_cast<
290  (*compIt->second)[j])))
291  {
292  // Retrieve basis key from expansion and define expansion
293  if((expIt = expansions.find(SegmentGeom->GetGlobalID())) != expansions.end())
294  {
295  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
296 
297  if(SetToOneSpaceDimension)
298  {
299  SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
300  SegmentGeom->GenerateOneSpaceDimGeom();
301 
303  ::AllocateSharedPtr(bkey, OneDSegmentGeom);
304  }
305  else
306  {
308  ::AllocateSharedPtr(bkey, SegmentGeom);
309  }
310  }
311  else
312  {
313  ASSERTL0(false,"Failed to find basis key");
314  }
315 
316  }
317  else
318  {
319  ASSERTL0(false,"Failed to dynamic cast geometry to SegGeom");
320  }
321 
322 
323  // Assign next ID
324  seg->SetElmtId(id++);
325 
326  // Add the expansion
327  (*m_exp).push_back(seg);
328  }
329  }
330 
331  // Setup Default optimisation information.
332  int nel = GetExpSize();
335 
336  // set up offset arrays.
338 
339  if(DeclareCoeffPhysArrays)
340  {
341  // Set up m_coeffs, m_phys.
344  }
345 
348  }
349 
350 
351  /**
352  * Fills the list of local expansions with the segments from the 2D
353  * mesh specified by \a domain. This CompositeMap contains a list of
354  * Composites which define the Neumann boundary.
355  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
356  * for details.
357  * @param domain A domain, comprising of one or more composite
358  * regions,
359  * @param graph2D A mesh, containing information about the
360  * domain and the spectral/hp element expansion.
361  * @param UseGenSegExp If true, create general segment expansions
362  * instead of just normal segment expansions.
363  */
365  const SpatialDomains::CompositeMap &domain,
366  const SpatialDomains::MeshGraphSharedPtr &graph2D,
367  const bool DeclareCoeffPhysArrays,
368  const std::string variable):
369  ExpList(pSession,graph2D)
370  {
371  SetExpType(e1D);
372 
373  m_graph = graph2D;
374 
375  int j, id=0;
377  SpatialDomains::CompositeMap::const_iterator compIt;
380 
381  // Process each composite region.
382  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
383  {
384  comp = compIt->second;
385  // Process each expansion in the region.
386  for(j = 0; j < compIt->second->size(); ++j)
387  {
388  if((SegmentGeom = boost::dynamic_pointer_cast<
390  (*compIt->second)[j])))
391  {
392  // Retrieve the basis key from the expansion.
394  = boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(graph2D)->GetEdgeBasisKey(SegmentGeom, variable);
395 
397  ::AllocateSharedPtr(bkey, SegmentGeom);
398 
399  // Add the segment to the expansion list.
400  seg->SetElmtId(id++);
401  (*m_exp).push_back(seg);
402  }
403  else
404  {
405  ASSERTL0(false,"dynamic cast to a SegGeom failed");
406  }
407  }
408 
409  }
410 
411  // Setup Default optimisation information.
412  int nel = GetExpSize();
415 
416  // Allocate storage for data and populate element offset lists.
418 
419  // Set up m_coeffs, m_phys.
420  if(DeclareCoeffPhysArrays)
421  {
424  }
425 
427  }
428 
429  /**
430  * Store expansions for the trace space expansions used in
431  * DisContField2D.
432  *
433  * @param bndConstraint Array of ExpList1D objects each containing a
434  * 1D spectral/hp element expansion on a single
435  * boundary region.
436  * @param bndCond Array of BoundaryCondition objects which contain
437  * information about the boundary conditions on the
438  * different boundary regions.
439  * @param locexp Complete domain expansion list.
440  * @param graph2D 2D mesh corresponding to the expansion list.
441  * @param periodicEdges List of periodic edges.
442  * @param UseGenSegExp If true, create general segment expansions
443  * instead of just normal segment expansions.
444  */
446  const LibUtilities::SessionReaderSharedPtr &pSession,
447  const Array<OneD,const ExpListSharedPtr> &bndConstraint,
449  const LocalRegions::ExpansionVector &locexp,
450  const SpatialDomains::MeshGraphSharedPtr &graph2D,
451  const PeriodicMap &periodicEdges,
452  const bool DeclareCoeffPhysArrays,
453  const std::string variable):
454  ExpList(pSession,graph2D)
455  {
456  int i, j, id, elmtid = 0;
457  set<int> edgesDone;
458 
465 
466  SetExpType(e1D);
467 
468  map<int,int> EdgeDone;
469  map<int,int> NormalSet;
470 
472 
473  // First loop over boundary conditions to renumber
474  // Dirichlet boundaries
475  for(i = 0; i < bndCond.num_elements(); ++i)
476  {
477  if(bndCond[i]->GetBoundaryConditionType()
479  {
480  for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
481  {
482  LibUtilities::BasisKey bkey = bndConstraint[i]
483  ->GetExp(j)->GetBasis(0)->GetBasisKey();
484  exp1D = bndConstraint[i]->GetExp(j)->
485  as<LocalRegions::Expansion1D>();
486  segGeom = exp1D->GetGeom1D();
487 
489  ::AllocateSharedPtr(bkey, segGeom);
490  edgesDone.insert(segGeom->GetEid());
491 
492  seg->SetElmtId(elmtid++);
493  (*m_exp).push_back(seg);
494  }
495  }
496  }
497 
499  LibUtilities::BasisKey> > edgeOrders;
502 
503  for(i = 0; i < locexp.size(); ++i)
504  {
505  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
506 
507  for(j = 0; j < locexp[i]->GetNedges(); ++j)
508  {
509  segGeom = exp2D->GetGeom2D()->GetEdge(j);
510  id = segGeom->GetEid();
511  // Ignore Dirichlet edges
512  if (edgesDone.count(id) != 0)
513  {
514  continue;
515  }
516 
517  it = edgeOrders.find(id);
518 
519  if (it == edgeOrders.end())
520  {
521  edgeOrders.insert(std::make_pair(id, std::make_pair(
522  segGeom, locexp[i]->DetEdgeBasisKey(j))));
523  }
524  else // variable modes/points
525  {
527  = locexp[i]->DetEdgeBasisKey(j);
528  LibUtilities::BasisKey existing
529  = it->second.second;
530 
531  int np1 = edge .GetNumPoints();
532  int np2 = existing.GetNumPoints();
533  int nm1 = edge .GetNumModes ();
534  int nm2 = existing.GetNumModes ();
535 
536  if (np2 >= np1 && nm2 >= nm1)
537  {
538  continue;
539  }
540  else if (np2 < np1 && nm2 < nm1)
541  {
542  it->second.second = edge;
543  }
544  else
545  {
546  ASSERTL0(false,
547  "inappropriate number of points/modes (max "
548  "num of points is not set with max order)");
549  }
550  }
551  }
552  }
553 
555  pSession->GetComm()->GetRowComm();
556  int nproc = vComm->GetSize(); // number of processors
557  int edgepr = vComm->GetRank(); // ID processor
558 
559  if (nproc > 1)
560  {
561  int eCnt = 0;
562 
563  // Count the number of edges on each partition
564  for(i = 0; i < locexp.size(); ++i)
565  {
566  eCnt += locexp[i]->GetNedges();
567  }
568 
569  // Set up the offset and the array that will contain the list of
570  // edge IDs, then reduce this across processors.
571  Array<OneD, int> edgesCnt(nproc, 0);
572  edgesCnt[edgepr] = eCnt;
573  vComm->AllReduce(edgesCnt, LibUtilities::ReduceSum);
574 
575  // Set up offset array.
576  int totEdgeCnt = Vmath::Vsum(nproc, edgesCnt, 1);
577  Array<OneD, int> eTotOffsets(nproc,0);
578  for (i = 1; i < nproc; ++i)
579  {
580  eTotOffsets[i] = eTotOffsets[i-1] + edgesCnt[i-1];
581  }
582 
583  // Local list of the edges per element
584  Array<OneD, int> EdgesTotID(totEdgeCnt, 0);
585  Array<OneD, int> EdgesTotNm(totEdgeCnt, 0);
586  Array<OneD, int> EdgesTotPnts(totEdgeCnt, 0);
587 
588  int cntr = eTotOffsets[edgepr];
589 
590  for(i = 0; i < locexp.size(); ++i)
591  {
592  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
593 
594  int nedges = locexp[i]->GetNedges();
595 
596  for(j = 0; j < nedges; ++j, ++cntr)
597  {
598  LibUtilities::BasisKey bkeyEdge =
599  locexp[i]->DetEdgeBasisKey(j);
600  EdgesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
601  EdgesTotNm [cntr] = bkeyEdge.GetNumModes();
602  EdgesTotPnts[cntr] = bkeyEdge.GetNumPoints();
603  }
604  }
605 
606  vComm->AllReduce(EdgesTotID, LibUtilities::ReduceSum);
607  vComm->AllReduce(EdgesTotNm, LibUtilities::ReduceSum);
608  vComm->AllReduce(EdgesTotPnts, LibUtilities::ReduceSum);
609 
610  for (i = 0; i < totEdgeCnt; ++i)
611  {
612  it = edgeOrders.find(EdgesTotID[i]);
613 
614  if (it == edgeOrders.end())
615  {
616  continue;
617  }
618 
619  LibUtilities::BasisKey existing
620  = it->second.second;
622  existing.GetBasisType(), EdgesTotNm[i],
623  LibUtilities::PointsKey(EdgesTotPnts[i],
624  existing.GetPointsType()));
625 
626 
627  int np1 = edge .GetNumPoints();
628  int np2 = existing.GetNumPoints();
629  int nm1 = edge .GetNumModes ();
630  int nm2 = existing.GetNumModes ();
631 
632  if (np2 >= np1 && nm2 >= nm1)
633  {
634  continue;
635  }
636  else if (np2 < np1 && nm2 < nm1)
637  {
638  it->second.second = edge;
639  }
640  else
641  {
642  ASSERTL0(false,
643  "inappropriate number of points/modes (max "
644  "num of points is not set with max order)");
645  }
646  }
647  }
648 
649  for (it = edgeOrders.begin(); it != edgeOrders.end(); ++it)
650  {
652  ::AllocateSharedPtr(it->second.second, it->second.first);
653  seg->SetElmtId(elmtid++);
654  (*m_exp).push_back(seg);
655  }
656 
657  // Setup Default optimisation information.
658  int nel = GetExpSize();
661 
662  // Set up offset information and array sizes
664 
665  // Set up m_coeffs, m_phys.
666  if(DeclareCoeffPhysArrays)
667  {
670  }
671 
673  }
674 
675  /**
676  * Each expansion (local element) is processed in turn to
677  * determine the number of coefficients and physical data
678  * points it contributes to the domain. Three arrays,
679  * #m_coeff_offset, #m_phys_offset and #m_offset_elmt_id, are
680  * also initialised and updated to store the data offsets of
681  * each element in the #m_coeffs and #m_phys arrays, and the
682  * element id that each consecutive block is associated
683  * respectively.
684  */
686  {
687  int i;
688 
689  // Set up offset information and array sizes
693 
694  m_ncoeffs = m_npoints = 0;
695 
696  for(i = 0; i < m_exp->size(); ++i)
697  {
699  m_phys_offset [i] = m_npoints;
700  m_offset_elmt_id[i] = i;
701  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
702  m_npoints += (*m_exp)[i]->GetTotPoints();
703  }
704  }
705 
706  /**
707  *
708  */
710  {
711  }
712 
713  /**
714  * To perform post-processing on the entire domain use \a elmtId = 0.
715  * @param kernel The post-processing kernel.
716  * @param inarray The set of evaluation points.
717  * @param outarray Contains the resulting post-processed
718  * solution for element \a elmId.
719  * @param h The mesh spacing.
720  * @param elmId Optionally specifies which element to perform
721  * the post-processing on (0=whole domain).
722  */
724  Array<OneD,NekDouble> &inarray,
725  Array<OneD,NekDouble> &outarray,
726  NekDouble h,
727  int elmId)
728 
729  {
730  int i,j,r;
731 
732  // get the local element expansion of the elmId element
734 
735  // Get the quadrature points and weights required for integration
736  int quad_npoints = elmExp->GetTotPoints();
737  LibUtilities::PointsKey quadPointsKey(quad_npoints,
738  elmExp->GetPointsType(0));
739  Array<OneD,NekDouble> quad_points
740  = LibUtilities::PointsManager()[quadPointsKey]->GetZ();
741  Array<OneD,NekDouble> quad_weights
742  = LibUtilities::PointsManager()[quadPointsKey]->GetW();
743 
744  // Declare variable for the local kernel breaks
745  int kernel_width = kernel->GetKernelWidth();
746  Array<OneD,NekDouble> local_kernel_breaks(kernel_width+1);
747 
748  // Declare variable for the transformed quadrature points
749  Array<OneD,NekDouble> mapped_quad_points(quad_npoints);
750 
751  // For each evaluation point
752  for(i = 0; i < inarray.num_elements(); i++)
753  {
754  // Move the center of the kernel to the current point
755  kernel->MoveKernelCenter(inarray[i],local_kernel_breaks);
756 
757  // Find the mesh breaks under the kernel support
758  Array<OneD,NekDouble> mesh_breaks;
759  kernel->FindMeshUnderKernel(local_kernel_breaks,h,mesh_breaks);
760 
761  // Sort the total breaks for integration purposes
762  int total_nbreaks = local_kernel_breaks.num_elements() +
763  mesh_breaks.num_elements();
764  // number of the total breaks
765  Array<OneD,NekDouble> total_breaks(total_nbreaks);
766  kernel->Sort(local_kernel_breaks,mesh_breaks,total_breaks);
767 
768  // Integrate the product of kernel and function over the total
769  // breaks
770  NekDouble integral_value = 0.0;
771  for(j = 0; j < total_breaks.num_elements()-1; j++)
772  {
773  NekDouble a = total_breaks[j];
774  NekDouble b = total_breaks[j+1];
775 
776  // Map the quadrature points to the appropriate interval
777  for(r = 0; r < quad_points.num_elements(); r++)
778  {
779  mapped_quad_points[r]
780  = (quad_points[r] + 1.0) * 0.5 * (b - a) + a;
781  }
782 
783  // Evaluate the function at the transformed quadrature
784  // points
785  Array<OneD,NekDouble> u_value(quad_npoints);
786  Array<OneD,NekDouble> coeffs = GetCoeffs();
787 
788  PeriodicEval(coeffs,mapped_quad_points,h,
789  elmExp->GetBasisNumModes(0),u_value);
790 
791  // Evaluate the kernel at the transformed quadrature points
792  Array<OneD,NekDouble> k_value(quad_npoints);
793  kernel->EvaluateKernel(mapped_quad_points,h,k_value);
794 
795  // Integrate
796  for(r = 0; r < quad_npoints; r++)
797  {
798  integral_value += (b - a) * 0.5 * k_value[r]
799  * u_value[r] * quad_weights[r];
800  }
801  }
802  outarray[i] = integral_value/h;
803  }
804  }
805 
806 
807  /**
808  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of an expansion,
809  * periodically evaluate the spectral/hp expansion
810  * \f$u^{\delta}(\boldsymbol{x})\f$ at arbitrary points.
811  * @param inarray1 An array of size \f$N_{\mathrm{eof}}\f$
812  * containing the local coefficients
813  * \f$\hat{u}_n^e\f$.
814  * @param inarray2 Contains the set of evaluation points.
815  * @param h The mesh spacing.
816  * @param nmodes The number of polynomial modes for each element
817  * (we consider that each element has the same
818  * number of polynomial modes).
819  * @param outarray Contains the resulting values at the
820  * evaluation points
821  */
823  Array<OneD,NekDouble> &inarray2,
824  NekDouble h, int nmodes,
825  Array<OneD,NekDouble> &outarray)
826  {
827  int i,j,r;
828 
829  // Get the number of elements in the domain
830  int num_elm = GetExpSize();
831 
832  // initializing the outarray
833  for(i = 0; i < outarray.num_elements(); i++)
834  {
835  outarray[i] = 0.0;
836  }
837 
838  // Make a copy for further modification
839  int x_size = inarray2.num_elements();
840  Array<OneD,NekDouble> x_values_cp(x_size);
841 
842  // Determining the element to which the x belongs
843  Array<OneD,int> x_elm(x_size);
844  for(i = 0; i < x_size; i++ )
845  {
846  x_elm[i] = (int)floor(inarray2[i]/h);
847  }
848 
849  // Clamp indices periodically
850  for(i = 0; i < x_size; i++)
851  {
852  while(x_elm[i] < 0)
853  {
854  x_elm[i] += num_elm;
855  }
856  while(x_elm[i] >= num_elm)
857  {
858  x_elm[i] -= num_elm ;
859  }
860  }
861 
862  // Map the values of x to [-1 1] on its interval
863  for(i = 0; i < x_size; i++)
864  {
865  x_values_cp[i] = (inarray2[i]/h - floor(inarray2[i]/h))*2 - 1.0;
866  }
867 
868  // Evaluate the jocobi polynomials
869  // (Evaluating the base at some points other than the quadrature
870  // points). Should it be added to the base class????
871  Array<TwoD,NekDouble> jacobi_poly(nmodes,x_size);
872  for(i = 0; i < nmodes; i++)
873  {
874  Polylib::jacobfd(x_size,x_values_cp.get(),
875  jacobi_poly.get()+i*x_size,NULL,i,0.0,0.0);
876  }
877 
878  // Evaluate the function values
879  for(r = 0; r < nmodes; r++)
880  {
881  for(j = 0; j < x_size; j++)
882  {
883  int index = ((x_elm[j])*nmodes)+r;
884  outarray[j] += inarray1[index]*jacobi_poly[r][j];
885  }
886  }
887 
888  }
889 
890 
891  /**
892  * Sets up the normals on all edges of expansions in the domain.
893  * @param locexp Complete list of domain expansions.
894  */
895 // void ExpList1D::SetUpPhysNormals(
896 // const StdRegions::StdExpansionVector &locexp)
897 // {
898 // map<int, int> EdgeGID;
899 // int i,cnt,n,id;
900 //
901 // // setup map of all global ids along boundary
902 // for(cnt = i = 0; i < (*m_exp).size(); ++i)
903 // {
904 // id = (*m_exp)[i]->GetGeom1D()->GetEid();
905 // EdgeGID[id] = cnt++;
906 // }
907 //
908 // // Loop over elements and find edges that match;
909 // for(cnt = n = 0; n < locexp.size(); ++n)
910 // {
911 // for(i = 0; i < locexp[n]->GetNedges(); ++i)
912 // {
913 // id = locexp[n]->GetGeom2D()->GetEid(i);
914 //
915 // if(EdgeGID.count(id) > 0)
916 // {
917 // (*m_exp)[EdgeGID.find(id)->second]
918 // ->SetUpPhysNormals(locexp[n],i);
919 // }
920 // }
921 // }
922 // }
923 
924  //croth
926  {
927  int i, j;
928  for (i = 0; i < m_exp->size(); ++i)
929  {
930  for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
931  {
932  (*m_exp)[i]->ComputeVertexNormal(j);
933  }
934  }
935  }
936 
937  /**
938  * Upwind the left and right states given by the Arrays Fwd and Bwd
939  * using the vector quantity Vec and ouput the upwinded value in the
940  * array upwind.
941  *
942  * @param Vec Velocity field.
943  * @param Fwd Left state.
944  * @param Bwd Right state.
945  * @param Upwind Output vector.
946  */
948  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
949  const Array<OneD, const NekDouble> &Fwd,
950  const Array<OneD, const NekDouble> &Bwd,
951  Array<OneD, NekDouble> &Upwind)
952  {
953  int i,j,k,e_npoints,offset;
954  Array<OneD,NekDouble> normals;
955  NekDouble Vn;
956 
957  // Assume whole array is of same coordimate dimension
958  int coordim = GetCoordim(0);
959 
960  ASSERTL1(Vec.num_elements() >= coordim,
961  "Input vector does not have sufficient dimensions to "
962  "match coordim");
963 
964  // Process each expansion
965  for(i = 0; i < m_exp->size(); ++i)
966  {
967  // Get the number of points in the expansion and the normals.
968  e_npoints = (*m_exp)[i]->GetNumPoints(0);
969  normals = (*m_exp)[i]->GetPhysNormals();
970 
971  // Get the physical data offset of the expansion in m_phys.
972  offset = m_phys_offset[i];
973 
974  // Compute each data point.
975  for(j = 0; j < e_npoints; ++j)
976  {
977  // Calculate normal velocity.
978  Vn = 0.0;
979  for(k = 0; k < coordim; ++k)
980  {
981  Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
982  }
983 
984  // Upwind based on direction of normal velocity.
985  if(Vn > 0.0)
986  {
987  Upwind[offset + j] = Fwd[offset + j];
988  }
989  else
990  {
991  Upwind[offset + j] = Bwd[offset + j];
992  }
993  }
994  }
995  }
996 
997  /**
998  * One-dimensional upwind.
999  * \see ExpList1D::Upwind(
1000  * const Array<OneD, const Array<OneD, NekDouble> >,
1001  * const Array<OneD, const NekDouble>,
1002  * const Array<OneD, const NekDouble>,
1003  * Array<OneD, NekDouble>, int)
1004  *
1005  * @param Vn Velocity field.
1006  * @param Fwd Left state.
1007  * @param Bwd Right state.
1008  * @param Upwind Output vector.
1009  */
1011  const Array<OneD, const NekDouble> &Vn,
1012  const Array<OneD, const NekDouble> &Fwd,
1013  const Array<OneD, const NekDouble> &Bwd,
1014  Array<OneD, NekDouble> &Upwind)
1015  {
1016  int i,j,e_npoints,offset;
1017  Array<OneD,NekDouble> normals;
1018 
1019  // Process each expansion.
1020  for(i = 0; i < m_exp->size(); ++i)
1021  {
1022  // Get the number of points and the data offset.
1023  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1024  offset = m_phys_offset[i];
1025 
1026  // Process each point in the expansion.
1027  for(j = 0; j < e_npoints; ++j)
1028  {
1029  // Upwind based on one-dimensional velocity.
1030  if(Vn[offset + j] > 0.0)
1031  {
1032  Upwind[offset + j] = Fwd[offset + j];
1033  }
1034  else
1035  {
1036  Upwind[offset + j] = Bwd[offset + j];
1037  }
1038  }
1039  }
1040  }
1041 
1042 
1043  /**
1044  * For each local element, copy the normals stored in the element list
1045  * into the array \a normals.
1046  * @param normals Multidimensional array in which to copy normals
1047  * to. Must have dimension equal to or larger than
1048  * the spatial dimension of the elements.
1049  */
1051  Array<OneD, Array<OneD, NekDouble> > &normals)
1052  {
1053  int i,j,k,e_npoints,offset;
1055  Array<OneD,Array<OneD,NekDouble> > locnormals;
1056  Array<OneD,Array<OneD,NekDouble> > locnormals2;
1058  // Assume whole array is of same coordinate dimension
1059  int coordim = GetCoordim(0);
1060 
1061  ASSERTL1(normals.num_elements() >= coordim,
1062  "Output vector does not have sufficient dimensions to "
1063  "match coordim");
1064 
1065  for (i = 0; i < m_exp->size(); ++i)
1066  {
1068 
1070  loc_exp->GetLeftAdjacentElementExp();
1071 
1072  int edgeNumber = loc_exp->GetLeftAdjacentElementEdge();
1073 
1074  // Get the number of points and normals for this expansion.
1075  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1076 
1077  locnormals = loc_elmt->GetEdgeNormal(edgeNumber);
1078  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
1079  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
1080 
1081  if (e_nmodes != loc_nmodes)
1082  {
1083  if (loc_exp->GetRightAdjacentElementEdge() >= 0)
1084  {
1086  loc_exp->GetRightAdjacentElementExp();
1087 
1088  int EdgeNumber = loc_exp->GetRightAdjacentElementEdge();
1089  // Serial case: right element is connected so we can
1090  // just grab that normal.
1091  locnormals = loc_elmt->GetEdgeNormal(EdgeNumber);
1092 
1093  offset = m_phys_offset[i];
1094 
1095  // Process each point in the expansion.
1096  for (j = 0; j < e_npoints; ++j)
1097  {
1098  // Process each spatial dimension and copy the values
1099  // into the output array.
1100  for (k = 0; k < coordim; ++k)
1101  {
1102  normals[k][offset + j] = -locnormals[k][j];
1103  }
1104  }
1105  }
1106  else
1107  {
1108  // Parallel case: need to interpolate normal.
1109  Array<OneD, Array<OneD, NekDouble> > normal(coordim);
1110 
1111  for (int p = 0; p < coordim; ++p)
1112  {
1113  normal[p] = Array<OneD, NekDouble>(e_npoints,0.0);
1114  LibUtilities::PointsKey to_key =
1115  loc_exp->GetBasis(0)->GetPointsKey();
1116  LibUtilities::PointsKey from_key =
1117  loc_elmt->GetBasis(0)->GetPointsKey();
1118  LibUtilities::Interp1D(from_key,
1119  locnormals[p],
1120  to_key,
1121  normal[p]);
1122  }
1123 
1124  offset = m_phys_offset[i];
1125 
1126  // Process each point in the expansion.
1127  for (j = 0; j < e_npoints; ++j)
1128  {
1129  // Process each spatial dimension and copy the values
1130  // into the output array.
1131  for (k = 0; k < coordim; ++k)
1132  {
1133  normals[k][offset + j] = normal[k][j];
1134  }
1135  }
1136  }
1137  }
1138  else
1139  {
1140  // Get the physical data offset for this expansion.
1141  offset = m_phys_offset[i];
1142 
1143  // Process each point in the expansion.
1144  for (j = 0; j < e_npoints; ++j)
1145  {
1146  // Process each spatial dimension and copy the values
1147  // into the output array.
1148  for (k = 0; k < coordim; ++k)
1149  {
1150  normals[k][offset + j] = locnormals[k][j];
1151  }
1152  }
1153  }
1154  }
1155  }
1156 
1157  /**
1158  *
1159  */
1161  {
1162 // Array<OneD, int> NumShape(1,0);
1163 // NumShape[0] = GetExpSize();
1164 //
1165 // int one = 1;
1166 // m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
1167 // ::AllocateSharedPtr(m_session,one,NumShape);
1168  }
1169 
1170 
1171  /**
1172  * Writes out the header for a <PIECE> VTK XML segment describing the
1173  * geometric information which comprises this element. This includes
1174  * vertex coordinates for each quadrature point, vertex connectivity
1175  * information, cell types and cell offset data.
1176  *
1177  * @param outfile Output stream to write data to.
1178  */
1179  void ExpList1D::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int)
1180  {
1181  int i,j;
1182  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1183  int ntot = nquad0;
1184  int ntotminus = (nquad0-1);
1185 
1186  Array<OneD,NekDouble> coords[3];
1187  coords[0] = Array<OneD,NekDouble>(ntot);
1188  coords[1] = Array<OneD,NekDouble>(ntot);
1189  coords[2] = Array<OneD,NekDouble>(ntot);
1190  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1191 
1192  outfile << " <Piece NumberOfPoints=\""
1193  << ntot << "\" NumberOfCells=\""
1194  << ntotminus << "\">" << endl;
1195  outfile << " <Points>" << endl;
1196  outfile << " <DataArray type=\"Float32\" "
1197  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1198  outfile << " ";
1199  for (i = 0; i < ntot; ++i)
1200  {
1201  for (j = 0; j < 3; ++j)
1202  {
1203  outfile << setprecision(8) << scientific
1204  << (float)coords[j][i] << " ";
1205  }
1206  outfile << endl;
1207  }
1208  outfile << endl;
1209  outfile << " </DataArray>" << endl;
1210  outfile << " </Points>" << endl;
1211  outfile << " <Cells>" << endl;
1212  outfile << " <DataArray type=\"Int32\" "
1213  << "Name=\"connectivity\" format=\"ascii\">" << endl;
1214  for (i = 0; i < nquad0-1; ++i)
1215  {
1216  outfile << i << " " << i+1 << endl;
1217  }
1218  outfile << endl;
1219  outfile << " </DataArray>" << endl;
1220  outfile << " <DataArray type=\"Int32\" "
1221  << "Name=\"offsets\" format=\"ascii\">" << endl;
1222  for (i = 0; i < ntotminus; ++i)
1223  {
1224  outfile << i*2+2 << " ";
1225  }
1226  outfile << endl;
1227  outfile << " </DataArray>" << endl;
1228  outfile << " <DataArray type=\"UInt8\" "
1229  << "Name=\"types\" format=\"ascii\">" << endl;
1230  for (i = 0; i < ntotminus; ++i)
1231  {
1232  outfile << "3 ";
1233  }
1234  outfile << endl;
1235  outfile << " </DataArray>" << endl;
1236  outfile << " </Cells>" << endl;
1237  outfile << " <PointData>" << endl;
1238  }
1239 
1240  } //end of namespace
1241 } //end of namespace
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:1775
ExpList1D()
The default constructor.
Definition: ExpList1D.cpp:75
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::map< int, vector< PeriodicEntity > > PeriodicMap
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
Definition: ExpList.h:971
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:1858
PointsType GetPointsType() const
Return type of quadrature.
Definition: Basis.h:151
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:926
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:909
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1837
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual void v_SetUpPhysNormals()
Set up the normals on each expansion.
Definition: ExpList1D.cpp:925
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion1D.h:116
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)
Upwind the Fwd and Bwd states based on the velocity field given by Vec.
Definition: ExpList1D.cpp:947
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:958
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:255
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:231
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:947
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:961
void PeriodicEval(Array< OneD, NekDouble > &inarray1, Array< OneD, NekDouble > &inarray2, NekDouble h, int nmodes, Array< OneD, NekDouble > &outarray)
Evaluates the global spectral/hp expansion at some arbitray set of points.
Definition: ExpList1D.cpp:822
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:887
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
Definition: ExpList.h:969
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
boost::shared_ptr< Kernel > KernelSharedPtr
Definition: kernel.h:215
PointsManagerT & PointsManager(void)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:883
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:111
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:112
void SetCoeffPhysOffsets(void)
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
Definition: ExpList1D.cpp:685
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Populate normals with the normals of all expansions.
Definition: ExpList1D.cpp:1050
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:54
This class is the abstraction of a one-dimensional multi-elemental expansions which is merely a colle...
Definition: ExpList1D.h:61
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry1D.h:48
void PostProcess(LibUtilities::KernelSharedPtr kernel, Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, NekDouble h, int elmId=0)
Performs the post-processing on a specified element.
Definition: ExpList1D.cpp:723
virtual void v_ReadGlobalOptimizationParameters()
Definition: ExpList1D.cpp:1160
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
Definition: ExpList.cpp:2723
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:1735
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
const StdRegions::StdExpansionVector &locexp);
Definition: ExpList1D.cpp:1179
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:211
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:714
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
int GetNedges() const
This function returns the number of edges of the expansion domain.
Definition: StdExpansion.h:272
virtual ~ExpList1D()
Destructor.
Definition: ExpList1D.cpp:709
Describes the specification for a Basis.
Definition: Basis.h:50
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1946
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49