Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  */
95  const std::vector<unsigned int> &eIDs,
96  const bool DeclareCoeffPhysArrays):
97  ExpList(In, eIDs, DeclareCoeffPhysArrays)
98  {
99  SetExpType(e1D);
100 
101  // Setup Default optimisation information.
102  int nel = GetExpSize();
105 
106  // Allocate storage for data and populate element offset lists.
108 
111  }
112 
113 
114  /**
115  * After initialising the data inherited through MultiRegions#ExpList,
116  * populate the expansion list from the segments defined in the supplied
117  * SpatialDomains#MeshGraph1D. All expansions in the graph are defined
118  * using the same LibUtilities#BasisKey which overrides that specified
119  * in \a graph1D.
120  *
121  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
122  * for details.
123  * @deprecated This constructor is no longer required since
124  * the basis key is now available from the graph.
125  * @param Ba BasisKey describing quadrature points and
126  * number of modes.
127  * @param graph1D Domain and expansion definitions.
128  */
130  const LibUtilities::BasisKey &Ba,
131  const SpatialDomains::MeshGraphSharedPtr &graph1D):
132  ExpList(pSession,graph1D)
133  {
134  SetExpType(e1D);
135 
136  int id=0;
139 
140  const SpatialDomains::ExpansionMap &expansions
141  = graph1D->GetExpansions();
142 
143  // For each element in the mesh, create a segment expansion using
144  // the supplied BasisKey and segment geometry.
145  SpatialDomains::ExpansionMap::const_iterator expIt;
146  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
147  {
148  if ((SegmentGeom = boost::dynamic_pointer_cast<
150  expIt->second->m_geomShPtr)))
151  {
153  ::AllocateSharedPtr(Ba,SegmentGeom);
154  seg->SetElmtId(id++);
155  (*m_exp).push_back(seg);
156  }
157  else
158  {
159  ASSERTL0(false,"dynamic cast to a SegGeom failed");
160  }
161  }
162 
163  // Setup Default optimisation information.
164  int nel = GetExpSize();
167 
168  // Allocate storage for data and populate element offset lists.
170 
173 
176  }
177 
178 
179  /**
180  * Given a mesh \a graph1D, containing information about the domain and
181  * the spectral/hp element expansion, this constructor fills the list
182  * of local expansions \texttt{m_exp} with the proper expansions,
183  * calculates the total number of quadrature points \f$x_i\f$ and local
184  * expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory for
185  * the arrays #m_coeffs and #m_phys.
186  *
187  * For each element its corresponding LibUtilities#BasisKey is
188  * retrieved and this is used to construct either a standard segment
189  * (LocalRegions#SegExp) or a generalised segment
190  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
191  * Finally, ExpList#SetCoeffPhys is called to initialise the data
192  * storage areas and set up the offset arrays.
193  *
194  * @param graph1D A mesh, containing information about the
195  * domain and the spectral/hp element expansion.
196  * @param UseGenSegExp If true, create general segment expansions
197  * instead of just normal segment expansions.
198  */
200  const SpatialDomains::MeshGraphSharedPtr &graph1D,
201  const bool DeclareCoeffPhysArrays):
202  ExpList(pSession,graph1D)
203  {
204  SetExpType(e1D);
205 
206  int id=0;
209 
210  // Retrieve the list of expansions
211  const SpatialDomains::ExpansionMap &expansions
212  = graph1D->GetExpansions();
213 
214  // Process each expansion in the graph
215  SpatialDomains::ExpansionMap::const_iterator expIt;
216  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
217  {
218  // Retrieve basis key from expansion
219  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
220 
221  if ((SegmentGeom = boost::dynamic_pointer_cast<
223  expIt->second->m_geomShPtr)))
224  {
226  ::AllocateSharedPtr(bkey, SegmentGeom);
227 
228  // Assign next ID
229  seg->SetElmtId(id++);
230 
231  // Add the expansion
232  (*m_exp).push_back(seg);
233  }
234  else
235  {
236  ASSERTL0(false,"dynamic cast to a SegGeom failed");
237  }
238  }
239 
240  // Setup Default optimisation information.
241  int nel = GetExpSize();
244 
245  // set up offset arrays.
247 
248  if(DeclareCoeffPhysArrays)
249  {
250  // Set up m_coeffs, m_phys.
253  }
254 
257  }
258 
259 
260 
261  /**
262  * Given a mesh \a graph1D, and the spectral/hp element
263  * expansion as well as and separate information about a \a
264  * domain, this constructor fills the list of local
265  * expansions \texttt{m_exp} with the proper expansions,
266  * calculates the total number of quadrature points \f$x_i\f$
267  * and local expansion coefficients \f$\hat{u}^e_n\f$ and
268  * allocates memory for the arrays #m_coeffs and #m_phys.
269  *
270  * For each element its corresponding LibUtilities#BasisKey is
271  * retrieved and this is used to construct either a standard segment
272  * (LocalRegions#SegExp) or a generalised segment
273  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
274  * Finally, ExpList#SetCoeffPhys is called to initialise the data
275  * storage areas and set up the offset arrays.
276  *
277  * @param graph1D A mesh, containing information about the
278  * domain and the spectral/hp element expansion.
279  * @param UseGenSegExp If true, create general segment expansions
280  * instead of just normal segment expansions.
281  */
283  const SpatialDomains::MeshGraphSharedPtr &graph1D,
284  const SpatialDomains::CompositeMap &domain,
285  const bool DeclareCoeffPhysArrays,
286  const std::string var,
287  bool SetToOneSpaceDimension):
288  ExpList(pSession,graph1D)
289  {
290  int j,id=0;
294  SpatialDomains::CompositeMap::const_iterator compIt;
295 
296  // Retrieve the list of expansions
297  const SpatialDomains::ExpansionMap &expansions
298  = graph1D->GetExpansions(var);
299 
300  // Process each composite region in domain
301  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
302  {
303  comp = compIt->second;
304 
305  // Process each expansion in the graph
306  for(j = 0; j < compIt->second->size(); ++j)
307  {
308  SpatialDomains::ExpansionMap::const_iterator expIt;
309 
310  if((SegmentGeom = boost::dynamic_pointer_cast<
312  (*compIt->second)[j])))
313  {
314  // Retrieve basis key from expansion and define expansion
315  if((expIt = expansions.find(SegmentGeom->GetGlobalID())) != expansions.end())
316  {
317  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
318 
319  if(SetToOneSpaceDimension)
320  {
321  SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
322  SegmentGeom->GenerateOneSpaceDimGeom();
323 
325  ::AllocateSharedPtr(bkey, OneDSegmentGeom);
326  }
327  else
328  {
330  ::AllocateSharedPtr(bkey, SegmentGeom);
331  }
332  }
333  else
334  {
335  ASSERTL0(false,"Failed to find basis key");
336  }
337 
338  }
339  else
340  {
341  ASSERTL0(false,"Failed to dynamic cast geometry to SegGeom");
342  }
343 
344 
345  // Assign next ID
346  seg->SetElmtId(id++);
347 
348  // Add the expansion
349  (*m_exp).push_back(seg);
350  }
351  }
352 
353  // Setup Default optimisation information.
354  int nel = GetExpSize();
357 
358  // set up offset arrays.
360 
361  if(DeclareCoeffPhysArrays)
362  {
363  // Set up m_coeffs, m_phys.
366  }
367 
370  }
371 
372 
373  /**
374  * Fills the list of local expansions with the segments from the 2D
375  * mesh specified by \a domain. This CompositeMap contains a list of
376  * Composites which define the Neumann boundary.
377  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
378  * for details.
379  * @param domain A domain, comprising of one or more composite
380  * regions,
381  * @param graph2D A mesh, containing information about the
382  * domain and the spectral/hp element expansion.
383  * @param UseGenSegExp If true, create general segment expansions
384  * instead of just normal segment expansions.
385  */
387  const SpatialDomains::CompositeMap &domain,
388  const SpatialDomains::MeshGraphSharedPtr &graph2D,
389  const bool DeclareCoeffPhysArrays,
390  const std::string variable):
391  ExpList(pSession,graph2D)
392  {
393  SetExpType(e1D);
394 
395  m_graph = graph2D;
396 
397  int j, id=0;
399  SpatialDomains::CompositeMap::const_iterator compIt;
402 
403  // Process each composite region.
404  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
405  {
406  comp = compIt->second;
407  // Process each expansion in the region.
408  for(j = 0; j < compIt->second->size(); ++j)
409  {
410  if((SegmentGeom = boost::dynamic_pointer_cast<
412  (*compIt->second)[j])))
413  {
414  // Retrieve the basis key from the expansion.
416  = boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(graph2D)->GetEdgeBasisKey(SegmentGeom, variable);
417 
419  ::AllocateSharedPtr(bkey, SegmentGeom);
420 
421  // Add the segment to the expansion list.
422  seg->SetElmtId(id++);
423  (*m_exp).push_back(seg);
424  }
425  else
426  {
427  ASSERTL0(false,"dynamic cast to a SegGeom failed");
428  }
429  }
430 
431  }
432 
433  // Setup Default optimisation information.
434  int nel = GetExpSize();
437 
438  // Allocate storage for data and populate element offset lists.
440 
441  // Set up m_coeffs, m_phys.
442  if(DeclareCoeffPhysArrays)
443  {
446  }
447 
449  }
450 
451  /**
452  * Store expansions for the trace space expansions used in
453  * DisContField2D.
454  *
455  * @param bndConstraint Array of ExpList1D objects each containing a
456  * 1D spectral/hp element expansion on a single
457  * boundary region.
458  * @param bndCond Array of BoundaryCondition objects which contain
459  * information about the boundary conditions on the
460  * different boundary regions.
461  * @param locexp Complete domain expansion list.
462  * @param graph2D 2D mesh corresponding to the expansion list.
463  * @param periodicEdges List of periodic edges.
464  * @param UseGenSegExp If true, create general segment expansions
465  * instead of just normal segment expansions.
466  */
468  const LibUtilities::SessionReaderSharedPtr &pSession,
469  const Array<OneD,const ExpListSharedPtr> &bndConstraint,
471  const LocalRegions::ExpansionVector &locexp,
472  const SpatialDomains::MeshGraphSharedPtr &graph2D,
473  const PeriodicMap &periodicEdges,
474  const bool DeclareCoeffPhysArrays,
475  const std::string variable):
476  ExpList(pSession,graph2D)
477  {
478  int i, j, id, elmtid = 0;
479  set<int> edgesDone;
480 
487 
488  SetExpType(e1D);
489 
490  map<int,int> EdgeDone;
491  map<int,int> NormalSet;
492 
494 
495  // First loop over boundary conditions to renumber
496  // Dirichlet boundaries
497  for(i = 0; i < bndCond.num_elements(); ++i)
498  {
499  if(bndCond[i]->GetBoundaryConditionType()
501  {
502  for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
503  {
504  LibUtilities::BasisKey bkey = bndConstraint[i]
505  ->GetExp(j)->GetBasis(0)->GetBasisKey();
506  exp1D = bndConstraint[i]->GetExp(j)->
507  as<LocalRegions::Expansion1D>();
508  segGeom = exp1D->GetGeom1D();
509 
511  ::AllocateSharedPtr(bkey, segGeom);
512  edgesDone.insert(segGeom->GetEid());
513 
514  seg->SetElmtId(elmtid++);
515  (*m_exp).push_back(seg);
516  }
517  }
518  }
519 
521  LibUtilities::BasisKey> > edgeOrders;
524 
525  for(i = 0; i < locexp.size(); ++i)
526  {
527  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
528 
529  for(j = 0; j < locexp[i]->GetNedges(); ++j)
530  {
531  segGeom = exp2D->GetGeom2D()->GetEdge(j);
532  id = segGeom->GetEid();
533  // Ignore Dirichlet edges
534  if (edgesDone.count(id) != 0)
535  {
536  continue;
537  }
538 
539  it = edgeOrders.find(id);
540 
541  if (it == edgeOrders.end())
542  {
543  edgeOrders.insert(std::make_pair(id, std::make_pair(
544  segGeom, locexp[i]->DetEdgeBasisKey(j))));
545  }
546  else // variable modes/points
547  {
549  = locexp[i]->DetEdgeBasisKey(j);
550  LibUtilities::BasisKey existing
551  = it->second.second;
552 
553  int np1 = edge .GetNumPoints();
554  int np2 = existing.GetNumPoints();
555  int nm1 = edge .GetNumModes ();
556  int nm2 = existing.GetNumModes ();
557 
558  if (np2 >= np1 && nm2 >= nm1)
559  {
560  continue;
561  }
562  else if (np2 < np1 && nm2 < nm1)
563  {
564  it->second.second = edge;
565  }
566  else
567  {
568  ASSERTL0(false,
569  "inappropriate number of points/modes (max "
570  "num of points is not set with max order)");
571  }
572  }
573  }
574  }
575 
577  pSession->GetComm()->GetRowComm();
578  int nproc = vComm->GetSize(); // number of processors
579  int edgepr = vComm->GetRank(); // ID processor
580 
581  if (nproc > 1)
582  {
583  int eCnt = 0;
584 
585  // Count the number of edges on each partition
586  for(i = 0; i < locexp.size(); ++i)
587  {
588  eCnt += locexp[i]->GetNedges();
589  }
590 
591  // Set up the offset and the array that will contain the list of
592  // edge IDs, then reduce this across processors.
593  Array<OneD, int> edgesCnt(nproc, 0);
594  edgesCnt[edgepr] = eCnt;
595  vComm->AllReduce(edgesCnt, LibUtilities::ReduceSum);
596 
597  // Set up offset array.
598  int totEdgeCnt = Vmath::Vsum(nproc, edgesCnt, 1);
599  Array<OneD, int> eTotOffsets(nproc,0);
600  for (i = 1; i < nproc; ++i)
601  {
602  eTotOffsets[i] = eTotOffsets[i-1] + edgesCnt[i-1];
603  }
604 
605  // Local list of the edges per element
606  Array<OneD, int> EdgesTotID(totEdgeCnt, 0);
607  Array<OneD, int> EdgesTotNm(totEdgeCnt, 0);
608  Array<OneD, int> EdgesTotPnts(totEdgeCnt, 0);
609 
610  int cntr = eTotOffsets[edgepr];
611 
612  for(i = 0; i < locexp.size(); ++i)
613  {
614  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
615 
616  int nedges = locexp[i]->GetNedges();
617 
618  for(j = 0; j < nedges; ++j, ++cntr)
619  {
620  LibUtilities::BasisKey bkeyEdge =
621  locexp[i]->DetEdgeBasisKey(j);
622  EdgesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
623  EdgesTotNm [cntr] = bkeyEdge.GetNumModes();
624  EdgesTotPnts[cntr] = bkeyEdge.GetNumPoints();
625  }
626  }
627 
628  vComm->AllReduce(EdgesTotID, LibUtilities::ReduceSum);
629  vComm->AllReduce(EdgesTotNm, LibUtilities::ReduceSum);
630  vComm->AllReduce(EdgesTotPnts, LibUtilities::ReduceSum);
631 
632  for (i = 0; i < totEdgeCnt; ++i)
633  {
634  it = edgeOrders.find(EdgesTotID[i]);
635 
636  if (it == edgeOrders.end())
637  {
638  continue;
639  }
640 
641  LibUtilities::BasisKey existing
642  = it->second.second;
644  existing.GetBasisType(), EdgesTotNm[i],
645  LibUtilities::PointsKey(EdgesTotPnts[i],
646  existing.GetPointsType()));
647 
648 
649  int np1 = edge .GetNumPoints();
650  int np2 = existing.GetNumPoints();
651  int nm1 = edge .GetNumModes ();
652  int nm2 = existing.GetNumModes ();
653 
654  if (np2 >= np1 && nm2 >= nm1)
655  {
656  continue;
657  }
658  else if (np2 < np1 && nm2 < nm1)
659  {
660  it->second.second = edge;
661  }
662  else
663  {
664  ASSERTL0(false,
665  "inappropriate number of points/modes (max "
666  "num of points is not set with max order)");
667  }
668  }
669  }
670 
671  for (it = edgeOrders.begin(); it != edgeOrders.end(); ++it)
672  {
674  ::AllocateSharedPtr(it->second.second, it->second.first);
675  seg->SetElmtId(elmtid++);
676  (*m_exp).push_back(seg);
677  }
678 
679  // Setup Default optimisation information.
680  int nel = GetExpSize();
683 
684  // Set up offset information and array sizes
686 
687  // Set up m_coeffs, m_phys.
688  if(DeclareCoeffPhysArrays)
689  {
692  }
693 
695  }
696 
697  /**
698  * Each expansion (local element) is processed in turn to
699  * determine the number of coefficients and physical data
700  * points it contributes to the domain. Three arrays,
701  * #m_coeff_offset, #m_phys_offset and #m_offset_elmt_id, are
702  * also initialised and updated to store the data offsets of
703  * each element in the #m_coeffs and #m_phys arrays, and the
704  * element id that each consecutive block is associated
705  * respectively.
706  */
708  {
709  int i;
710 
711  // Set up offset information and array sizes
715 
716  m_ncoeffs = m_npoints = 0;
717 
718  for(i = 0; i < m_exp->size(); ++i)
719  {
721  m_phys_offset [i] = m_npoints;
722  m_offset_elmt_id[i] = i;
723  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
724  m_npoints += (*m_exp)[i]->GetTotPoints();
725  }
726  }
727 
728  /**
729  *
730  */
732  {
733  }
734 
735  /**
736  * To perform post-processing on the entire domain use \a elmtId = 0.
737  * @param kernel The post-processing kernel.
738  * @param inarray The set of evaluation points.
739  * @param outarray Contains the resulting post-processed
740  * solution for element \a elmId.
741  * @param h The mesh spacing.
742  * @param elmId Optionally specifies which element to perform
743  * the post-processing on (0=whole domain).
744  */
746  Array<OneD,NekDouble> &inarray,
747  Array<OneD,NekDouble> &outarray,
748  NekDouble h,
749  int elmId)
750 
751  {
752  int i,j,r;
753 
754  // get the local element expansion of the elmId element
756 
757  // Get the quadrature points and weights required for integration
758  int quad_npoints = elmExp->GetTotPoints();
759  LibUtilities::PointsKey quadPointsKey(quad_npoints,
760  elmExp->GetPointsType(0));
761  Array<OneD,NekDouble> quad_points
762  = LibUtilities::PointsManager()[quadPointsKey]->GetZ();
763  Array<OneD,NekDouble> quad_weights
764  = LibUtilities::PointsManager()[quadPointsKey]->GetW();
765 
766  // Declare variable for the local kernel breaks
767  int kernel_width = kernel->GetKernelWidth();
768  Array<OneD,NekDouble> local_kernel_breaks(kernel_width+1);
769 
770  // Declare variable for the transformed quadrature points
771  Array<OneD,NekDouble> mapped_quad_points(quad_npoints);
772 
773  // For each evaluation point
774  for(i = 0; i < inarray.num_elements(); i++)
775  {
776  // Move the center of the kernel to the current point
777  kernel->MoveKernelCenter(inarray[i],local_kernel_breaks);
778 
779  // Find the mesh breaks under the kernel support
780  Array<OneD,NekDouble> mesh_breaks;
781  kernel->FindMeshUnderKernel(local_kernel_breaks,h,mesh_breaks);
782 
783  // Sort the total breaks for integration purposes
784  int total_nbreaks = local_kernel_breaks.num_elements() +
785  mesh_breaks.num_elements();
786  // number of the total breaks
787  Array<OneD,NekDouble> total_breaks(total_nbreaks);
788  kernel->Sort(local_kernel_breaks,mesh_breaks,total_breaks);
789 
790  // Integrate the product of kernel and function over the total
791  // breaks
792  NekDouble integral_value = 0.0;
793  for(j = 0; j < total_breaks.num_elements()-1; j++)
794  {
795  NekDouble a = total_breaks[j];
796  NekDouble b = total_breaks[j+1];
797 
798  // Map the quadrature points to the appropriate interval
799  for(r = 0; r < quad_points.num_elements(); r++)
800  {
801  mapped_quad_points[r]
802  = (quad_points[r] + 1.0) * 0.5 * (b - a) + a;
803  }
804 
805  // Evaluate the function at the transformed quadrature
806  // points
807  Array<OneD,NekDouble> u_value(quad_npoints);
808  Array<OneD,NekDouble> coeffs = GetCoeffs();
809 
810  PeriodicEval(coeffs,mapped_quad_points,h,
811  elmExp->GetBasisNumModes(0),u_value);
812 
813  // Evaluate the kernel at the transformed quadrature points
814  Array<OneD,NekDouble> k_value(quad_npoints);
815  kernel->EvaluateKernel(mapped_quad_points,h,k_value);
816 
817  // Integrate
818  for(r = 0; r < quad_npoints; r++)
819  {
820  integral_value += (b - a) * 0.5 * k_value[r]
821  * u_value[r] * quad_weights[r];
822  }
823  }
824  outarray[i] = integral_value/h;
825  }
826  }
827 
828 
829  /**
830  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of an expansion,
831  * periodically evaluate the spectral/hp expansion
832  * \f$u^{\delta}(\boldsymbol{x})\f$ at arbitrary points.
833  * @param inarray1 An array of size \f$N_{\mathrm{eof}}\f$
834  * containing the local coefficients
835  * \f$\hat{u}_n^e\f$.
836  * @param inarray2 Contains the set of evaluation points.
837  * @param h The mesh spacing.
838  * @param nmodes The number of polynomial modes for each element
839  * (we consider that each element has the same
840  * number of polynomial modes).
841  * @param outarray Contains the resulting values at the
842  * evaluation points
843  */
845  Array<OneD,NekDouble> &inarray2,
846  NekDouble h, int nmodes,
847  Array<OneD,NekDouble> &outarray)
848  {
849  int i,j,r;
850 
851  // Get the number of elements in the domain
852  int num_elm = GetExpSize();
853 
854  // initializing the outarray
855  for(i = 0; i < outarray.num_elements(); i++)
856  {
857  outarray[i] = 0.0;
858  }
859 
860  // Make a copy for further modification
861  int x_size = inarray2.num_elements();
862  Array<OneD,NekDouble> x_values_cp(x_size);
863 
864  // Determining the element to which the x belongs
865  Array<OneD,int> x_elm(x_size);
866  for(i = 0; i < x_size; i++ )
867  {
868  x_elm[i] = (int)floor(inarray2[i]/h);
869  }
870 
871  // Clamp indices periodically
872  for(i = 0; i < x_size; i++)
873  {
874  while(x_elm[i] < 0)
875  {
876  x_elm[i] += num_elm;
877  }
878  while(x_elm[i] >= num_elm)
879  {
880  x_elm[i] -= num_elm ;
881  }
882  }
883 
884  // Map the values of x to [-1 1] on its interval
885  for(i = 0; i < x_size; i++)
886  {
887  x_values_cp[i] = (inarray2[i]/h - floor(inarray2[i]/h))*2 - 1.0;
888  }
889 
890  // Evaluate the jocobi polynomials
891  // (Evaluating the base at some points other than the quadrature
892  // points). Should it be added to the base class????
893  Array<TwoD,NekDouble> jacobi_poly(nmodes,x_size);
894  for(i = 0; i < nmodes; i++)
895  {
896  Polylib::jacobfd(x_size,x_values_cp.get(),
897  jacobi_poly.get()+i*x_size,NULL,i,0.0,0.0);
898  }
899 
900  // Evaluate the function values
901  for(r = 0; r < nmodes; r++)
902  {
903  for(j = 0; j < x_size; j++)
904  {
905  int index = ((x_elm[j])*nmodes)+r;
906  outarray[j] += inarray1[index]*jacobi_poly[r][j];
907  }
908  }
909 
910  }
911 
912 
913  /**
914  * Sets up the normals on all edges of expansions in the domain.
915  * @param locexp Complete list of domain expansions.
916  */
917 // void ExpList1D::SetUpPhysNormals(
918 // const StdRegions::StdExpansionVector &locexp)
919 // {
920 // map<int, int> EdgeGID;
921 // int i,cnt,n,id;
922 //
923 // // setup map of all global ids along boundary
924 // for(cnt = i = 0; i < (*m_exp).size(); ++i)
925 // {
926 // id = (*m_exp)[i]->GetGeom1D()->GetEid();
927 // EdgeGID[id] = cnt++;
928 // }
929 //
930 // // Loop over elements and find edges that match;
931 // for(cnt = n = 0; n < locexp.size(); ++n)
932 // {
933 // for(i = 0; i < locexp[n]->GetNedges(); ++i)
934 // {
935 // id = locexp[n]->GetGeom2D()->GetEid(i);
936 //
937 // if(EdgeGID.count(id) > 0)
938 // {
939 // (*m_exp)[EdgeGID.find(id)->second]
940 // ->SetUpPhysNormals(locexp[n],i);
941 // }
942 // }
943 // }
944 // }
945 
946  //croth
948  {
949  int i, j;
950  for (i = 0; i < m_exp->size(); ++i)
951  {
952  for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
953  {
954  (*m_exp)[i]->ComputeVertexNormal(j);
955  }
956  }
957  }
958 
959  /**
960  * Upwind the left and right states given by the Arrays Fwd and Bwd
961  * using the vector quantity Vec and ouput the upwinded value in the
962  * array upwind.
963  *
964  * @param Vec Velocity field.
965  * @param Fwd Left state.
966  * @param Bwd Right state.
967  * @param Upwind Output vector.
968  */
970  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
971  const Array<OneD, const NekDouble> &Fwd,
972  const Array<OneD, const NekDouble> &Bwd,
973  Array<OneD, NekDouble> &Upwind)
974  {
975  int i,j,k,e_npoints,offset;
976  Array<OneD,NekDouble> normals;
977  NekDouble Vn;
978 
979  // Assume whole array is of same coordimate dimension
980  int coordim = GetCoordim(0);
981 
982  ASSERTL1(Vec.num_elements() >= coordim,
983  "Input vector does not have sufficient dimensions to "
984  "match coordim");
985 
986  // Process each expansion
987  for(i = 0; i < m_exp->size(); ++i)
988  {
989  // Get the number of points in the expansion and the normals.
990  e_npoints = (*m_exp)[i]->GetNumPoints(0);
991  normals = (*m_exp)[i]->GetPhysNormals();
992 
993  // Get the physical data offset of the expansion in m_phys.
994  offset = m_phys_offset[i];
995 
996  // Compute each data point.
997  for(j = 0; j < e_npoints; ++j)
998  {
999  // Calculate normal velocity.
1000  Vn = 0.0;
1001  for(k = 0; k < coordim; ++k)
1002  {
1003  Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
1004  }
1005 
1006  // Upwind based on direction of normal velocity.
1007  if(Vn > 0.0)
1008  {
1009  Upwind[offset + j] = Fwd[offset + j];
1010  }
1011  else
1012  {
1013  Upwind[offset + j] = Bwd[offset + j];
1014  }
1015  }
1016  }
1017  }
1018 
1019  /**
1020  * One-dimensional upwind.
1021  * \see ExpList1D::Upwind(
1022  * const Array<OneD, const Array<OneD, NekDouble> >,
1023  * const Array<OneD, const NekDouble>,
1024  * const Array<OneD, const NekDouble>,
1025  * Array<OneD, NekDouble>, int)
1026  *
1027  * @param Vn Velocity field.
1028  * @param Fwd Left state.
1029  * @param Bwd Right state.
1030  * @param Upwind Output vector.
1031  */
1033  const Array<OneD, const NekDouble> &Vn,
1034  const Array<OneD, const NekDouble> &Fwd,
1035  const Array<OneD, const NekDouble> &Bwd,
1036  Array<OneD, NekDouble> &Upwind)
1037  {
1038  int i,j,e_npoints,offset;
1039  Array<OneD,NekDouble> normals;
1040 
1041  // Process each expansion.
1042  for(i = 0; i < m_exp->size(); ++i)
1043  {
1044  // Get the number of points and the data offset.
1045  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1046  offset = m_phys_offset[i];
1047 
1048  // Process each point in the expansion.
1049  for(j = 0; j < e_npoints; ++j)
1050  {
1051  // Upwind based on one-dimensional velocity.
1052  if(Vn[offset + j] > 0.0)
1053  {
1054  Upwind[offset + j] = Fwd[offset + j];
1055  }
1056  else
1057  {
1058  Upwind[offset + j] = Bwd[offset + j];
1059  }
1060  }
1061  }
1062  }
1063 
1064 
1065  /**
1066  * For each local element, copy the normals stored in the element list
1067  * into the array \a normals.
1068  * @param normals Multidimensional array in which to copy normals
1069  * to. Must have dimension equal to or larger than
1070  * the spatial dimension of the elements.
1071  */
1073  Array<OneD, Array<OneD, NekDouble> > &normals)
1074  {
1075  int i,j,k,e_npoints,offset;
1077  Array<OneD,Array<OneD,NekDouble> > locnormals;
1078  Array<OneD,Array<OneD,NekDouble> > locnormals2;
1080  // Assume whole array is of same coordinate dimension
1081  int coordim = GetCoordim(0);
1082 
1083  ASSERTL1(normals.num_elements() >= coordim,
1084  "Output vector does not have sufficient dimensions to "
1085  "match coordim");
1086 
1087  for (i = 0; i < m_exp->size(); ++i)
1088  {
1090 
1092  loc_exp->GetLeftAdjacentElementExp();
1093 
1094  int edgeNumber = loc_exp->GetLeftAdjacentElementEdge();
1095 
1096  // Get the number of points and normals for this expansion.
1097  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1098 
1099  locnormals = loc_elmt->GetEdgeNormal(edgeNumber);
1100  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
1101  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
1102 
1103  if (e_nmodes != loc_nmodes)
1104  {
1105  if (loc_exp->GetRightAdjacentElementEdge() >= 0)
1106  {
1108  loc_exp->GetRightAdjacentElementExp();
1109 
1110  int EdgeNumber = loc_exp->GetRightAdjacentElementEdge();
1111  // Serial case: right element is connected so we can
1112  // just grab that normal.
1113  locnormals = loc_elmt->GetEdgeNormal(EdgeNumber);
1114 
1115  offset = m_phys_offset[i];
1116 
1117  // Process each point in the expansion.
1118  for (j = 0; j < e_npoints; ++j)
1119  {
1120  // Process each spatial dimension and copy the values
1121  // into the output array.
1122  for (k = 0; k < coordim; ++k)
1123  {
1124  normals[k][offset + j] = -locnormals[k][j];
1125  }
1126  }
1127  }
1128  else
1129  {
1130  // Parallel case: need to interpolate normal.
1131  Array<OneD, Array<OneD, NekDouble> > normal(coordim);
1132 
1133  for (int p = 0; p < coordim; ++p)
1134  {
1135  normal[p] = Array<OneD, NekDouble>(e_npoints,0.0);
1136  LibUtilities::PointsKey to_key =
1137  loc_exp->GetBasis(0)->GetPointsKey();
1138  LibUtilities::PointsKey from_key =
1139  loc_elmt->GetBasis(0)->GetPointsKey();
1140  LibUtilities::Interp1D(from_key,
1141  locnormals[p],
1142  to_key,
1143  normal[p]);
1144  }
1145 
1146  offset = m_phys_offset[i];
1147 
1148  // Process each point in the expansion.
1149  for (j = 0; j < e_npoints; ++j)
1150  {
1151  // Process each spatial dimension and copy the values
1152  // into the output array.
1153  for (k = 0; k < coordim; ++k)
1154  {
1155  normals[k][offset + j] = normal[k][j];
1156  }
1157  }
1158  }
1159  }
1160  else
1161  {
1162  // Get the physical data offset for this expansion.
1163  offset = m_phys_offset[i];
1164 
1165  // Process each point in the expansion.
1166  for (j = 0; j < e_npoints; ++j)
1167  {
1168  // Process each spatial dimension and copy the values
1169  // into the output array.
1170  for (k = 0; k < coordim; ++k)
1171  {
1172  normals[k][offset + j] = locnormals[k][j];
1173  }
1174  }
1175  }
1176  }
1177  }
1178 
1179  /**
1180  *
1181  */
1183  {
1184 // Array<OneD, int> NumShape(1,0);
1185 // NumShape[0] = GetExpSize();
1186 //
1187 // int one = 1;
1188 // m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
1189 // ::AllocateSharedPtr(m_session,one,NumShape);
1190  }
1191 
1192 
1193  /**
1194  * Writes out the header for a <PIECE> VTK XML segment describing the
1195  * geometric information which comprises this element. This includes
1196  * vertex coordinates for each quadrature point, vertex connectivity
1197  * information, cell types and cell offset data.
1198  *
1199  * @param outfile Output stream to write data to.
1200  */
1201  void ExpList1D::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int)
1202  {
1203  int i,j;
1204  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1205  int ntot = nquad0;
1206  int ntotminus = (nquad0-1);
1207 
1208  Array<OneD,NekDouble> coords[3];
1209  coords[0] = Array<OneD,NekDouble>(ntot, 0.0);
1210  coords[1] = Array<OneD,NekDouble>(ntot, 0.0);
1211  coords[2] = Array<OneD,NekDouble>(ntot, 0.0);
1212  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1213 
1214  outfile << " <Piece NumberOfPoints=\""
1215  << ntot << "\" NumberOfCells=\""
1216  << ntotminus << "\">" << endl;
1217  outfile << " <Points>" << endl;
1218  outfile << " <DataArray type=\"Float32\" "
1219  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1220  outfile << " ";
1221  for (i = 0; i < ntot; ++i)
1222  {
1223  for (j = 0; j < 3; ++j)
1224  {
1225  outfile << setprecision(8) << scientific
1226  << (float)coords[j][i] << " ";
1227  }
1228  outfile << endl;
1229  }
1230  outfile << endl;
1231  outfile << " </DataArray>" << endl;
1232  outfile << " </Points>" << endl;
1233  outfile << " <Cells>" << endl;
1234  outfile << " <DataArray type=\"Int32\" "
1235  << "Name=\"connectivity\" format=\"ascii\">" << endl;
1236  for (i = 0; i < nquad0-1; ++i)
1237  {
1238  outfile << i << " " << i+1 << endl;
1239  }
1240  outfile << endl;
1241  outfile << " </DataArray>" << endl;
1242  outfile << " <DataArray type=\"Int32\" "
1243  << "Name=\"offsets\" format=\"ascii\">" << endl;
1244  for (i = 0; i < ntotminus; ++i)
1245  {
1246  outfile << i*2+2 << " ";
1247  }
1248  outfile << endl;
1249  outfile << " </DataArray>" << endl;
1250  outfile << " <DataArray type=\"UInt8\" "
1251  << "Name=\"types\" format=\"ascii\">" << endl;
1252  for (i = 0; i < ntotminus; ++i)
1253  {
1254  outfile << "3 ";
1255  }
1256  outfile << endl;
1257  outfile << " </DataArray>" << endl;
1258  outfile << " </Cells>" << endl;
1259  outfile << " <PointData>" << endl;
1260  }
1261 
1262  } //end of namespace
1263 } //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:1834
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:1001
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:1917
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:956
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:939
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1896
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual void v_SetUpPhysNormals()
Set up the normals on each expansion.
Definition: ExpList1D.cpp:947
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion1D.h:123
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:969
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:988
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
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:269
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:977
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:991
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:844
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:917
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
Definition: ExpList.h:999
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:913
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:115
void SetCoeffPhysOffsets(void)
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
Definition: ExpList1D.cpp:707
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:1072
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:745
virtual void v_ReadGlobalOptimizationParameters()
Definition: ExpList1D.cpp:1182
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:2950
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:1794
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:1201
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:251
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:723
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:442
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:731
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:174
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49