Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ExpList2DHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList2DHomogeneous1D.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: An ExpList2D which is homogeneous in 1 direction and so
33 // uses much of the functionality from a ExpList1D and its daughters
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <MultiRegions/ExpList1D.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  namespace MultiRegions
45  {
46  // Forward declaration for typedefs
47  ExpList2DHomogeneous1D::ExpList2DHomogeneous1D():
49  {
50  }
51 
52  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
55  const LibUtilities::BasisKey &HomoBasis,
56  const NekDouble lhom,
57  const bool useFFT,
58  const bool dealiasing,
59  const Array<OneD, ExpListSharedPtr> &planes)
60  : ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
61  {
62  int i, n, cnt, nel;
63 
64  ASSERTL1(m_planes.num_elements() == planes.num_elements(),
65  "Size of basis number of points and number"
66  "of planes are not the same");
67 
68  // Set up expansion list with elements from all planes.
71  planes.num_elements() * planes[0]->GetExpSize());
72 
73  for(cnt = n = 0; n < planes.num_elements(); ++n)
74  {
75  m_planes[n] = planes[n];
76  for (i = 0; i < planes[n]->GetExpSize(); ++i)
77  {
78  (*m_exp)[cnt++] = planes[n]->GetExp(i);
79  }
80  }
81 
82  // Setup Default optimisation information.
83  nel = GetExpSize();
86 
87  SetCoeffPhys();
88  }
89 
90  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
93  const LibUtilities::BasisKey &HomoBasis,
94  const NekDouble lhom,
95  const bool useFFT,
96  const bool dealiasing,
97  const SpatialDomains::MeshGraphSharedPtr &graph1D):
98  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
99  {
100  int n, j, nel;
101  ExpList1DSharedPtr plane_zero;
102 
103  // note that nzplanes can be larger than nzmodes
104  m_planes[0] = plane_zero = MemoryManager<ExpList1D>::
105  AllocateSharedPtr(m_session, graph1D, false);
106 
108 
109  nel = m_planes[0]->GetExpSize();
110 
111  for (j = 0; j < nel; ++j)
112  {
113  (*m_exp).push_back(m_planes[0]->GetExp(j));
114  }
115 
116  for (n = 1; n < m_planes.num_elements(); ++n)
117  {
119  AllocateSharedPtr(*plane_zero, false);
120  for(j = 0; j < nel; ++j)
121  {
122  (*m_exp).push_back((*m_exp)[j]);
123  }
124  }
125 
126  // Setup Default optimisation information.
127  nel = GetExpSize();
130 
131  SetCoeffPhys();
132  }
133 
134 
135  /**
136  * @param In ExpList2DHomogeneous1D object to copy.
137  */
139  const ExpList2DHomogeneous1D &In):
141  {
142  ExpList1DSharedPtr zero_plane =
143  boost::dynamic_pointer_cast<ExpList1D> (In.m_planes[0]);
144 
145  for (int n = 0; n < m_planes.num_elements(); ++n)
146  {
148  AllocateSharedPtr(*zero_plane, false);
149  }
150 
151  SetCoeffPhys();
152  }
153 
154  /**
155  * Destructor
156  */
158  {
159  }
160 
162  {
163  int i, n, cnt;
164  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
165  int npoints_per_plane = m_planes[0]->GetTotPoints();
166 
167  int nzplanes = m_planes.num_elements();
168 
169  // Set total coefficients and points
170  m_ncoeffs = ncoeffs_per_plane*nzplanes;
171  m_npoints = npoints_per_plane*nzplanes;
172 
175 
176  int nel = m_planes[0]->GetExpSize();
177  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
178  m_phys_offset = Array<OneD,int>(nel*nzplanes);
179  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
180  Array<OneD, NekDouble> tmparray;
181 
182  for (cnt = n = 0; n < nzplanes; ++n)
183  {
184  m_planes[n]->SetCoeffsArray(
185  tmparray= m_coeffs + ncoeffs_per_plane*n);
186  m_planes[n]->SetPhysArray(
187  tmparray = m_phys + npoints_per_plane*n);
188 
189  for(i = 0; i < nel; ++i)
190  {
191  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i)
192  + n*ncoeffs_per_plane;
193  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i)
194  + n*npoints_per_plane;
195  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i)
196  + n*nel;
197  }
198  }
199  }
200 
202  const int eid,
206  {
207  int n, coordim;
208  Array<OneD, NekDouble> tmp_xc,xhom;
209  int nyplanes = m_planes.num_elements();
210  int npoints = GetTotPoints(eid);
211 
212  switch(coordim = GetCoordim(0))
213  {
214  case 1:
215  {
216  (*m_exp)[eid]->GetCoords(xc0);
217  xhom = xc1;
218  }
219  break;
220  case 2:
221  ASSERTL0(xc1.num_elements() != 0,
222  "output coord_1 is not defined");
223  {
224  (*m_exp)[eid]->GetCoords(xc0,xc1);
225  xhom = xc2;
226  }
227  break;
228  default:
229  ASSERTL0(false, "Cannot have coordim being three dimensions"
230  "in a homogeneous field");
231  break;
232  }
233 
234  // Fill homogeneous-direction
236  Array<OneD, NekDouble> z(nyplanes);
237 
238  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
239 
240  for(n = 0; n < m_planes.num_elements(); n++)
241  {
242  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
243  }
244 
245  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
246  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
247 
248  for (n = 0; n < nyplanes; ++n)
249  {
250  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
251  if(n)
252  {
253  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
254  if(coordim == 2)
255  {
256  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
257  }
258  }
259  }
260  }
261 
262  /**
263  * The operation calls the 2D plane coordinates through the
264  * function ExpList#GetCoords and then evaluated the third
265  * coordinate using the member \a m_lhom
266  *
267  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
268  * will be stored in this array.
269  *
270  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
271  * will be stored in this array. This
272  * coordinate might be evaluated using the
273  * predefined value \a m_lhom
274  *
275  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
276  * will be stored in this array. This
277  * coordinate is evaluated using the
278  * predefined value \a m_lhom
279  */
284  {
285  int n,coordim;
286  Array<OneD, NekDouble> tmp_xc, xhom;
287  int nyplanes = m_planes.num_elements();
288  int npoints = m_planes[0]->GetTotPoints();
289 
290  m_planes[0]->GetCoords(xc0,xc1);
291 
292  if((coordim = GetCoordim(0)) == 1)
293  {
294  xhom = xc1;
295  }
296  else
297  {
298  xhom = xc2;
299  }
300 
301  // Fill z-direction
303  Array<OneD, NekDouble> z(nyplanes);
304  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
305 
306  for(n = 0; n < m_planes.num_elements(); n++)
307  {
308  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
309  }
310 
311  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
312  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
313 
314  for(n = 0; n < nyplanes; ++n)
315  {
316  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
317  if(n)
318  {
319  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
320  if(coordim == 2)
321  {
322  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
323  }
324  }
325  }
326  }
327 
328 
329  /**
330  * Write Tecplot Files Zone
331  * @param outfile Output file name.
332  * @param expansion Expansion that is considered
333  */
335  std::ostream &outfile,
336  int expansion)
337  {
338  int i, j;
339 
340  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
341  int nquad1 = m_planes.num_elements();
342 
343  Array<OneD,NekDouble> coords[3];
344 
345  coords[0] = Array<OneD,NekDouble>(3*nquad0*nquad1);
346  coords[1] = coords[0] + nquad0*nquad1;
347  coords[2] = coords[1] + nquad0*nquad1;
348 
349  GetCoords(expansion,coords[0],coords[1],coords[2]);
350 
351  outfile << "Zone, I=" << nquad0 << ", J=" << nquad1
352  << ", F=Block" << std::endl;
353 
354  for(j = 0; j < GetCoordim(0)+1; ++j)
355  {
356  for(i = 0; i < nquad0*nquad1; ++i)
357  {
358  outfile << coords[j][i] << " ";
359  }
360  outfile << std::endl;
361  }
362  }
363 
364 
366  std::ostream &outfile,
367  int expansion,
368  int istrip)
369  {
370  // If there is only one plane (e.g. HalfMode), we write a 2D plane.
371  if (m_planes.num_elements() == 1)
372  {
373  m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
374  return;
375  }
376 
377  // If we are using Fourier points, output extra plane to fill domain
378  int outputExtraPlane = 0;
379  if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
380  && m_homogeneousBasis->GetPointsType() ==
382  {
383  outputExtraPlane = 1;
384  }
385  int i, j;
386  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
387  int nquad1 = m_planes.num_elements() + outputExtraPlane;
388  int ntot = nquad0*nquad1;
389  int ntotminus = (nquad0-1)*(nquad1-1);
390 
391  Array<OneD,NekDouble> coords[3];
392  coords[0] = Array<OneD,NekDouble>(ntot);
393  coords[1] = Array<OneD,NekDouble>(ntot);
394  coords[2] = Array<OneD,NekDouble>(ntot);
395  GetCoords(expansion,coords[0],coords[1],coords[2]);
396 
397  if (outputExtraPlane)
398  {
399  // Copy coords[0] and coords[1] to extra plane
401  Vmath::Vcopy (nquad0, coords[0], 1,
402  tmp = coords[0] + (nquad1-1)*nquad0, 1);
403  Vmath::Vcopy (nquad0, coords[1], 1,
404  tmp = coords[1] + (nquad1-1)*nquad0, 1);
405  // Fill coords[2] for extra plane
406  NekDouble z = coords[2][nquad0*m_planes.num_elements()-1] +
407  (coords[2][nquad0] - coords[2][0]);
408  Vmath::Fill(nquad0, z, tmp = coords[2] + (nquad1-1)*nquad0, 1);
409  }
410 
411  outfile << " <Piece NumberOfPoints=\""
412  << ntot << "\" NumberOfCells=\""
413  << ntotminus << "\">" << endl;
414  outfile << " <Points>" << endl;
415  outfile << " <DataArray type=\"Float64\" "
416  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
417  outfile << " ";
418  for (i = 0; i < ntot; ++i)
419  {
420  for (j = 0; j < 3; ++j)
421  {
422  outfile << coords[j][i] << " ";
423  }
424  outfile << endl;
425  }
426  outfile << endl;
427  outfile << " </DataArray>" << endl;
428  outfile << " </Points>" << endl;
429  outfile << " <Cells>" << endl;
430  outfile << " <DataArray type=\"Int32\" "
431  << "Name=\"connectivity\" format=\"ascii\">" << endl;
432  for (i = 0; i < nquad0-1; ++i)
433  {
434  for (j = 0; j < nquad1-1; ++j)
435  {
436  outfile << j*nquad0 + i << " "
437  << j*nquad0 + i + 1 << " "
438  << (j+1)*nquad0 + i + 1 << " "
439  << (j+1)*nquad0 + i << endl;
440  }
441  }
442  outfile << endl;
443  outfile << " </DataArray>" << endl;
444  outfile << " <DataArray type=\"Int32\" "
445  << "Name=\"offsets\" format=\"ascii\">" << endl;
446  for (i = 0; i < ntotminus; ++i)
447  {
448  outfile << i*4+4 << " ";
449  }
450  outfile << endl;
451  outfile << " </DataArray>" << endl;
452  outfile << " <DataArray type=\"UInt8\" "
453  << "Name=\"types\" format=\"ascii\">" << endl;
454  for (i = 0; i < ntotminus; ++i)
455  {
456  outfile << "9 ";
457  }
458  outfile << endl;
459  outfile << " </DataArray>" << endl;
460  outfile << " </Cells>" << endl;
461  outfile << " <PointData>" << endl;
462  }
463 
465  Array<OneD, Array<OneD, NekDouble> > &normals)
466  {
467  int nPlanes = m_planes.num_elements();
468  int nPtsPlane = m_planes[0]->GetNpoints();
469  int nDim = GetCoordim(0) + 1;
470 
471  ASSERTL1(normals.num_elements() >= nDim,
472  "Output vector does not have sufficient dimensions to"
473  "match coordim");
474  ASSERTL1(normals[0].num_elements() >= nPtsPlane,
475  "Output vector does not have sufficient dimensions to"
476  "match coordim");
477 
478  // Calculate normals from plane 0.
479  m_planes[0]->GetNormals(normals);
480 
481  // Copy remaining planes.
482  for (int i = 0; i < nDim; ++i)
483  {
484  for (int n = 1; n < nPlanes; ++n)
485  {
486  Vmath::Vcopy(nPtsPlane, &normals[i][0], 1,
487  &normals[i][n*nPtsPlane], 1);
488  }
489  }
490  }
491  } //end of namespace
492 } //end of namespace
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
#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:1060
LibUtilities::TranspositionSharedPtr m_transposition
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2075
STL namespace.
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
NekDouble m_lhom
Width of homogeneous direction.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1015
Fourier Expansion .
Definition: BasisType.h:52
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:2054
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1047
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:66
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1036
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1543
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1050
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:976
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:1058
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:969
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2)
Array< OneD, ExpListSharedPtr > m_planes
double NekDouble
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:315
void SetCoeffPhys(void)
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
This class is the abstraction of a one-dimensional multi-elemental expansions which is merely a colle...
Definition: ExpList1D.h:61
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:1906
void GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
This function calculates the coordinates of all the elemental quadrature points . ...
#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
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Describes the specification for a Basis.
Definition: Basis.h:50