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