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