Nektar++
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 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  // Forward declaration for typedefs
47  {
48  }
49 
50  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
53  const LibUtilities::BasisKey &HomoBasis,
54  const NekDouble lhom,
55  const bool useFFT,
56  const bool dealiasing,
57  const Array<OneD, ExpListSharedPtr> &planes)
58  : ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
59  {
60  int i, n, cnt, nel;
61 
62  ASSERTL1(m_planes.num_elements() == planes.num_elements(),
63  "Size of basis number of points and number"
64  "of planes are not the same");
65 
66  // Set up expansion list with elements from all planes.
69  planes.num_elements() * planes[0]->GetExpSize());
70 
71  for(cnt = n = 0; n < planes.num_elements(); ++n)
72  {
73  m_planes[n] = planes[n];
74  for (i = 0; i < planes[n]->GetExpSize(); ++i)
75  {
76  (*m_exp)[cnt++] = planes[n]->GetExp(i);
77  }
78  }
79 
80  // Setup Default optimisation information.
81  nel = GetExpSize();
84 
85  SetCoeffPhys();
86  }
87 
88  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
91  const LibUtilities::BasisKey &HomoBasis,
92  const NekDouble lhom,
93  const bool useFFT,
94  const bool dealiasing,
95  const SpatialDomains::MeshGraphSharedPtr &graph1D):
96  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
97  {
98  int n, j, nel;
99  ExpList1DSharedPtr plane_zero;
100 
101  // note that nzplanes can be larger than nzmodes
102  m_planes[0] = plane_zero = MemoryManager<ExpList1D>::
103  AllocateSharedPtr(m_session, graph1D, false);
104 
106 
107  nel = m_planes[0]->GetExpSize();
108 
109  for (j = 0; j < nel; ++j)
110  {
111  (*m_exp).push_back(m_planes[0]->GetExp(j));
112  }
113 
114  for (n = 1; n < m_planes.num_elements(); ++n)
115  {
117  AllocateSharedPtr(*plane_zero, false);
118  for(j = 0; j < nel; ++j)
119  {
120  (*m_exp).push_back((*m_exp)[j]);
121  }
122  }
123 
124  // Setup Default optimisation information.
125  nel = GetExpSize();
128 
129  SetCoeffPhys();
130  }
131 
132 
133  /**
134  * @param In ExpList2DHomogeneous1D object to copy.
135  */
137  const ExpList2DHomogeneous1D &In):
139  {
140  ExpList1DSharedPtr zero_plane =
141  boost::dynamic_pointer_cast<ExpList1D> (In.m_planes[0]);
142 
143  for (int n = 0; n < m_planes.num_elements(); ++n)
144  {
146  AllocateSharedPtr(*zero_plane, false);
147  }
148 
149  SetCoeffPhys();
150  }
151 
152  /**
153  * Destructor
154  */
156  {
157  }
158 
160  {
161  int i, n, cnt;
162  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
163  int npoints_per_plane = m_planes[0]->GetTotPoints();
164 
165  int nzplanes = m_planes.num_elements();
166 
167  // Set total coefficients and points
168  m_ncoeffs = ncoeffs_per_plane*nzplanes;
169  m_npoints = npoints_per_plane*nzplanes;
170 
173 
174  int nel = m_planes[0]->GetExpSize();
175  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
176  m_phys_offset = Array<OneD,int>(nel*nzplanes);
177  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
178  Array<OneD, NekDouble> tmparray;
179 
180  for (cnt = n = 0; n < nzplanes; ++n)
181  {
182  m_planes[n]->SetCoeffsArray(
183  tmparray= m_coeffs + ncoeffs_per_plane*n);
184  m_planes[n]->SetPhysArray(
185  tmparray = m_phys + npoints_per_plane*n);
186 
187  for(i = 0; i < nel; ++i)
188  {
189  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i)
190  + n*ncoeffs_per_plane;
191  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i)
192  + n*npoints_per_plane;
193  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i)
194  + n*nel;
195  }
196  }
197  }
198 
200  const int eid,
204  {
205  int n, coordim;
206  Array<OneD, NekDouble> tmp_xc,xhom;
207  int nyplanes = m_planes.num_elements();
208  int npoints = GetTotPoints(eid);
209 
210  switch(coordim = GetCoordim(0))
211  {
212  case 1:
213  {
214  (*m_exp)[eid]->GetCoords(xc0);
215  xhom = xc1;
216  }
217  break;
218  case 2:
219  ASSERTL0(xc1.num_elements() != 0,
220  "output coord_1 is not defined");
221  {
222  (*m_exp)[eid]->GetCoords(xc0,xc1);
223  xhom = xc2;
224  }
225  break;
226  default:
227  ASSERTL0(false, "Cannot have coordim being three dimensions"
228  "in a homogeneous field");
229  break;
230  }
231 
232  // Fill homogeneous-direction
234  Array<OneD, NekDouble> z(nyplanes);
235 
236  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
237 
238  for(n = 0; n < m_planes.num_elements(); n++)
239  {
240  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
241  }
242 
243  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
244  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
245 
246  for (n = 0; n < nyplanes; ++n)
247  {
248  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
249  if(n)
250  {
251  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
252  if(coordim == 2)
253  {
254  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
255  }
256  }
257  }
258  }
259 
260  /**
261  * The operation calls the 2D plane coordinates through the
262  * function ExpList#GetCoords and then evaluated the third
263  * coordinate using the member \a m_lhom
264  *
265  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
266  * will be stored in this array.
267  *
268  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
269  * will be stored in this array. This
270  * coordinate might be evaluated using the
271  * predefined value \a m_lhom
272  *
273  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
274  * will be stored in this array. This
275  * coordinate is evaluated using the
276  * predefined value \a m_lhom
277  */
282  {
283  int n,coordim;
284  Array<OneD, NekDouble> tmp_xc, xhom;
285  int nyplanes = m_planes.num_elements();
286  int npoints = m_planes[0]->GetTotPoints();
287 
288  m_planes[0]->GetCoords(xc0,xc1);
289 
290  if((coordim = GetCoordim(0)) == 1)
291  {
292  xhom = xc1;
293  }
294  else
295  {
296  xhom = xc2;
297  }
298 
299  // Fill z-direction
301  Array<OneD, NekDouble> z(nyplanes);
302  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
303 
304  for(n = 0; n < m_planes.num_elements(); n++)
305  {
306  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
307  }
308 
309  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
310  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
311 
312  for(n = 0; n < nyplanes; ++n)
313  {
314  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
315  if(n)
316  {
317  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
318  if(coordim == 2)
319  {
320  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
321  }
322  }
323  }
324  }
325 
326 
327  /**
328  * Write Tecplot Files Zone
329  * @param outfile Output file name.
330  * @param expansion Expansion that is considered
331  */
333  std::ostream &outfile,
334  int expansion)
335  {
336  int i, j;
337 
338  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
339  int nquad1 = m_planes.num_elements();
340 
341  Array<OneD,NekDouble> coords[3];
342 
343  coords[0] = Array<OneD,NekDouble>(3*nquad0*nquad1);
344  coords[1] = coords[0] + nquad0*nquad1;
345  coords[2] = coords[1] + nquad0*nquad1;
346 
347  GetCoords(expansion,coords[0],coords[1],coords[2]);
348 
349  outfile << "Zone, I=" << nquad0 << ", J=" << nquad1
350  << ", F=Block" << std::endl;
351 
352  for(j = 0; j < GetCoordim(0)+1; ++j)
353  {
354  for(i = 0; i < nquad0*nquad1; ++i)
355  {
356  outfile << coords[j][i] << " ";
357  }
358  outfile << std::endl;
359  }
360  }
361 
362 
364  std::ostream &outfile,
365  int expansion,
366  int istrip)
367  {
368  int i, j;
369  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
370  int nquad1 = m_planes.num_elements();
371  int ntot = nquad0*nquad1;
372  int ntotminus = (nquad0-1)*(nquad1-1);
373 
374  Array<OneD,NekDouble> coords[3];
375  coords[0] = Array<OneD,NekDouble>(ntot);
376  coords[1] = Array<OneD,NekDouble>(ntot);
377  coords[2] = Array<OneD,NekDouble>(ntot);
378  GetCoords(expansion,coords[0],coords[1],coords[2]);
379 
380  outfile << " <Piece NumberOfPoints=\""
381  << ntot << "\" NumberOfCells=\""
382  << ntotminus << "\">" << endl;
383  outfile << " <Points>" << endl;
384  outfile << " <DataArray type=\"Float64\" "
385  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
386  outfile << " ";
387  for (i = 0; i < ntot; ++i)
388  {
389  for (j = 0; j < 3; ++j)
390  {
391  outfile << coords[j][i] << " ";
392  }
393  outfile << endl;
394  }
395  outfile << endl;
396  outfile << " </DataArray>" << endl;
397  outfile << " </Points>" << endl;
398  outfile << " <Cells>" << endl;
399  outfile << " <DataArray type=\"Int32\" "
400  << "Name=\"connectivity\" format=\"ascii\">" << endl;
401  for (i = 0; i < nquad0-1; ++i)
402  {
403  for (j = 0; j < nquad1-1; ++j)
404  {
405  outfile << j*nquad0 + i << " "
406  << j*nquad0 + i + 1 << " "
407  << (j+1)*nquad0 + i + 1 << " "
408  << (j+1)*nquad0 + i << endl;
409  }
410  }
411  outfile << endl;
412  outfile << " </DataArray>" << endl;
413  outfile << " <DataArray type=\"Int32\" "
414  << "Name=\"offsets\" format=\"ascii\">" << endl;
415  for (i = 0; i < ntotminus; ++i)
416  {
417  outfile << i*4+4 << " ";
418  }
419  outfile << endl;
420  outfile << " </DataArray>" << endl;
421  outfile << " <DataArray type=\"UInt8\" "
422  << "Name=\"types\" format=\"ascii\">" << endl;
423  for (i = 0; i < ntotminus; ++i)
424  {
425  outfile << "9 ";
426  }
427  outfile << endl;
428  outfile << " </DataArray>" << endl;
429  outfile << " </Cells>" << endl;
430  outfile << " <PointData>" << endl;
431  }
432 
434  Array<OneD, Array<OneD, NekDouble> > &normals)
435  {
436  int nPlanes = m_planes.num_elements();
437  int nPtsPlane = m_planes[0]->GetNpoints();
438  int nDim = GetCoordim(0) + 1;
439 
440  ASSERTL1(normals.num_elements() >= nDim,
441  "Output vector does not have sufficient dimensions to"
442  "match coordim");
443  ASSERTL1(normals[0].num_elements() >= nPtsPlane,
444  "Output vector does not have sufficient dimensions to"
445  "match coordim");
446 
447  // Calculate normals from plane 0.
448  m_planes[0]->GetNormals(normals);
449 
450  // Copy remaining planes.
451  for (int i = 0; i < nDim; ++i)
452  {
453  for (int n = 1; n < nPlanes; ++n)
454  {
455  Vmath::Vcopy(nPtsPlane, &normals[i][0], 1,
456  &normals[i][n*nPtsPlane], 1);
457  }
458  }
459  }
460  } //end of namespace
461 } //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:161
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
Definition: ExpList.h:971
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:1858
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:926
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:909
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1837
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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:958
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
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:199
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:947
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1401
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:961
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:887
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
Definition: ExpList.h:969
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:880
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:301
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:1735
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:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
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:1038
Describes the specification for a Basis.
Definition: Basis.h:50