Nektar++
ExpList3DHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList3DHomogeneous1D.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 ExpList which is homogeneous in 1 direction and so
33 // uses much of the functionality from a ExpList2D and its daughters
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <MultiRegions/ExpList2D.h>
39 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  // Forward declaration for typedefs
47  {
49  }
50 
53  const LibUtilities::BasisKey &HomoBasis,
54  const NekDouble lhom,
55  const bool useFFT,
56  const bool dealiasing):
57  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
58  {
60  }
61 
62  // Constructor for ExpList3DHomogeneous1D to act as a Explist2D field
65  const LibUtilities::BasisKey &HomoBasis,
66  const NekDouble lhom,
67  const bool useFFT,
68  const bool dealiasing,
70  const std::string &var):
71  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
72  {
74 
75  GenExpList3DHomogeneous1D(graph2D->GetExpansions(var));
76  }
77 
78  // Constructor for ExpList3DHomogeneous1D to act as a Explist2D field
81  const LibUtilities::BasisKey &HomoBasis,
82  const NekDouble lhom,
83  const bool useFFT,
84  const bool dealiasing,
85  const SpatialDomains::ExpansionMap &expansions):
86  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
87  {
89 
90  GenExpList3DHomogeneous1D(expansions);
91  }
92 
94  {
95  int n,j,nel;
96  bool False = false;
97  ExpList2DSharedPtr plane_zero;
98 
99  // note that nzplanes can be larger than nzmodes
100  m_planes[0] = plane_zero = MemoryManager<ExpList2D>::AllocateSharedPtr(m_session, expansions, False);
101 
103  nel = m_planes[0]->GetExpSize();
104 
105  for(j = 0; j < nel; ++j)
106  {
107  (*m_exp).push_back(m_planes[0]->GetExp(j));
108  }
109 
110  for(n = 1; n < m_homogeneousBasis->GetNumPoints(); ++n)
111  {
113  for(j = 0; j < nel; ++j)
114  {
115  (*m_exp).push_back((*m_exp)[j]);
116  }
117  }
118 
119  // Setup Default optimisation information.
120  nel = GetExpSize();
123 
124  SetCoeffPhys();
125  }
126 
127 
128  /**
129  * @param In ExpList3DHomogeneous1D object to copy.
130  */
131  ExpList3DHomogeneous1D::ExpList3DHomogeneous1D(const ExpList3DHomogeneous1D &In, bool DeclarePlanesSetCoeffPhys):
133  {
135 
136  if(DeclarePlanesSetCoeffPhys)
137  {
138  bool False = false;
139  ExpList2DSharedPtr zero_plane = boost::dynamic_pointer_cast<ExpList2D> (In.m_planes[0]);
140 
141  for(int n = 0; n < m_planes.num_elements(); ++n)
142  {
144  }
145 
146  SetCoeffPhys();
147  }
148  }
149 
150  /**
151  * Destructor
152  */
154  {
155  }
156 
158  {
159  int i,n,cnt;
160  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
161  int npoints_per_plane = m_planes[0]->GetTotPoints();
162 
163  int nzplanes = m_planes.num_elements();
164 
165  // Set total coefficients and points
166  m_ncoeffs = ncoeffs_per_plane*nzplanes;
167  m_npoints = npoints_per_plane*nzplanes;
168 
171 
172  int nel = m_planes[0]->GetExpSize();
173  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
174  m_phys_offset = Array<OneD,int>(nel*nzplanes);
175  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
176  Array<OneD, NekDouble> tmparray;
177 
178  for(cnt = n = 0; n < nzplanes; ++n)
179  {
180  m_planes[n]->SetCoeffsArray(tmparray= m_coeffs + ncoeffs_per_plane*n);
181  m_planes[n]->SetPhysArray(tmparray = m_phys + npoints_per_plane*n);
182 
183  for(i = 0; i < nel; ++i)
184  {
185  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i) + n*ncoeffs_per_plane;
186  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i) + n*npoints_per_plane;
187  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i) + n*nel;
188  }
189  }
190  }
191 
196  {
197  int n;
198  Array<OneD, NekDouble> tmp_xc;
199  int nzplanes = m_planes.num_elements();
200  int npoints = GetTotPoints(eid);
201 
202  (*m_exp)[eid]->GetCoords(xc0,xc1);
203 
204  // Fill z-direction
206  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
207 
208  for(n = 0; n < m_planes.num_elements(); n++)
209  {
210  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
211  }
212 
213  Array<OneD, NekDouble> z(nzplanes);
214 
215  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
216  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
217 
218  for(n = 0; n < nzplanes; ++n)
219  {
220  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
221  if(n)
222  {
223  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
224  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
225  }
226  }
227  }
228 
229  /**
230  * The operation calls the 2D plane coordinates through the
231  * function ExpList#GetCoords and then evaluated the third
232  * coordinate using the member \a m_lhom
233  *
234  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
235  * will be stored in this array.
236  *
237  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
238  * will be stored in this array.
239  *
240  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
241  * will be stored in this array. This
242  * coordinate is evaluated using the
243  * predefined value \a m_lhom
244  */
248  {
249  int n;
250  Array<OneD, NekDouble> tmp_xc;
251  int nzplanes = m_planes.num_elements();
252  int npoints = m_planes[0]->GetTotPoints();
253 
254  m_planes[0]->GetCoords(xc0,xc1);
255 
256  // Fill z-direction
258 
259  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
260 
261  for(n = 0; n < m_planes.num_elements(); n++)
262  {
263  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
264  }
265 
266  Array<OneD, NekDouble> z(nzplanes);
267 
268  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
269  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
270 
271  for(n = 0; n < nzplanes; ++n)
272  {
273  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
274  if(n)
275  {
276  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
277  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
278  }
279  }
280  }
281 
282  void ExpList3DHomogeneous1D::v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
283  {
284  ASSERTL0(expansion == -1, "Multi-zone output not supported for homogeneous expansions.");
285 
286  const int nPtsPlane = m_planes[0]->GetNpoints();
287  const int nElmt = m_planes[0]->GetExpSize();
288  const int nPlanes = m_planes.num_elements();
289 
290  int cnt = 0;
291  int cnt2 = 0;
292  for (int i = 0; i < nElmt; ++i)
293  {
294  const int np0 = (*m_exp)[i]->GetNumPoints(0);
295  const int np1 = (*m_exp)[i]->GetNumPoints(1);
296 
297  for (int n = 1; n < nPlanes; ++n)
298  {
299  const int o1 = (n-1) * nPtsPlane;
300  const int o2 = n * nPtsPlane;
301  for (int j = 1; j < np1; ++j)
302  {
303  for(int k = 1; k < np0; ++k)
304  {
305  outfile << cnt + (j-1)*np0 + (k-1) + o1 + 1 << " ";
306  outfile << cnt + (j-1)*np0 + (k-1) + o2 + 1 << " ";
307  outfile << cnt + (j-1)*np0 + k + o2 + 1 << " ";
308  outfile << cnt + (j-1)*np0 + k + o1 + 1 << " ";
309  outfile << cnt + j *np0 + (k-1) + o1 + 1 << " ";
310  outfile << cnt + j *np0 + (k-1) + o2 + 1 << " ";
311  outfile << cnt + j *np0 + k + o2 + 1 << " ";
312  outfile << cnt + j *np0 + k + o1 + 1 << endl;
313  cnt2++;
314  }
315  }
316  }
317 
318  cnt += np0*np1;
319  }
320  }
321 
322 
324  std::ostream &outfile,
325  int expansion,
326  int istrip)
327  {
328  int i,j,k;
329  int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
330  int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
331  int nq2 = m_planes.num_elements();
332  int ntot = nq0*nq1*nq2;
333  int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
334 
335  Array<OneD,NekDouble> coords[3];
336  coords[0] = Array<OneD,NekDouble>(ntot);
337  coords[1] = Array<OneD,NekDouble>(ntot);
338  coords[2] = Array<OneD,NekDouble>(ntot);
339  GetCoords(expansion,coords[0],coords[1],coords[2]);
340 
341  NekDouble DistStrip;
342  m_session->LoadParameter("DistStrip", DistStrip, 0);
343  // Reset the z-coords for homostrips
344  for(int i = 0; i < ntot; i++)
345  {
346  coords[2][i] += istrip*DistStrip;
347  }
348 
349  outfile << " <Piece NumberOfPoints=\""
350  << ntot << "\" NumberOfCells=\""
351  << ntotminus << "\">" << endl;
352  outfile << " <Points>" << endl;
353  outfile << " <DataArray type=\"Float32\" "
354  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
355  outfile << " ";
356  for (i = 0; i < ntot; ++i)
357  {
358  for (j = 0; j < 3; ++j)
359  {
360  outfile << coords[j][i] << " ";
361  }
362  outfile << endl;
363  }
364  outfile << endl;
365  outfile << " </DataArray>" << endl;
366  outfile << " </Points>" << endl;
367  outfile << " <Cells>" << endl;
368  outfile << " <DataArray type=\"Int32\" "
369  << "Name=\"connectivity\" format=\"ascii\">" << endl;
370  for (i = 0; i < nq0-1; ++i)
371  {
372  for (j = 0; j < nq1-1; ++j)
373  {
374  for (k = 0; k < nq2-1; ++k)
375  {
376  outfile << k*nq0*nq1 + j*nq0 + i << " "
377  << k*nq0*nq1 + j*nq0 + i + 1 << " "
378  << k*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
379  << k*nq0*nq1 + (j+1)*nq0 + i << " "
380  << (k+1)*nq0*nq1 + j*nq0 + i << " "
381  << (k+1)*nq0*nq1 + j*nq0 + i + 1 << " "
382  << (k+1)*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
383  << (k+1)*nq0*nq1 + (j+1)*nq0 + i << endl;
384  }
385  }
386  }
387  outfile << endl;
388  outfile << " </DataArray>" << endl;
389  outfile << " <DataArray type=\"Int32\" "
390  << "Name=\"offsets\" format=\"ascii\">" << endl;
391  for (i = 0; i < ntotminus; ++i)
392  {
393  outfile << i*8+8 << " ";
394  }
395  outfile << endl;
396  outfile << " </DataArray>" << endl;
397  outfile << " <DataArray type=\"UInt8\" "
398  << "Name=\"types\" format=\"ascii\">" << endl;
399  for (i = 0; i < ntotminus; ++i)
400  {
401  outfile << "12 ";
402  }
403  outfile << endl;
404  outfile << " </DataArray>" << endl;
405  outfile << " </Cells>" << endl;
406  outfile << " <PointData>" << endl;
407  }
408 
410  const Array<OneD, const NekDouble> &inarray,
411  const Array<OneD, const NekDouble> &soln)
412  {
413  int cnt = 0;
414  NekDouble errL2, err = 0.0;
416  Array<OneD, NekDouble> local_w(m_planes.num_elements());
417 
418  for(int n = 0; n < m_planes.num_elements(); n++)
419  {
420  local_w[n] = w[m_transposition->GetPlaneID(n)];
421  }
422 
423  if (soln == NullNekDouble1DArray)
424  {
425  for(int n = 0; n < m_planes.num_elements(); ++n)
426  {
427  errL2 = m_planes[n]->L2(inarray + cnt);
428  cnt += m_planes[n]->GetTotPoints();
429  err += errL2*errL2*local_w[n]*m_lhom*0.5;
430  }
431  }
432  else
433  {
434  for(int n = 0; n < m_planes.num_elements(); ++n)
435  {
436  errL2 = m_planes[n]->L2(inarray + cnt, soln + cnt);
437  cnt += m_planes[n]->GetTotPoints();
438  err += errL2*errL2*local_w[n]*m_lhom*0.5;
439  }
440  }
441  m_comm->GetColumnComm()->AllReduce(err, LibUtilities::ReduceSum);
442 
443  return sqrt(err);
444  }
445 
447  {
448  Array<OneD, NekDouble> energy(m_planes.num_elements()/2);
449  NekDouble area = 0.0;
450 
451  // Calculate total area of elements.
452  for (int n = 0; n < m_planes[0]->GetExpSize(); ++n)
453  {
454  Array<OneD, NekDouble> inarray(m_planes[0]->GetExp(n)->GetTotPoints(), 1.0);
455  area += m_planes[0]->GetExp(n)->Integral(inarray);
456  }
457 
458  m_comm->GetRowComm()->AllReduce(area, LibUtilities::ReduceSum);
459 
460  // Calculate L2 norm of real/imaginary planes.
461  for (int n = 0; n < m_planes.num_elements(); n += 2)
462  {
463  NekDouble err;
464 
465  energy[n/2] = 0;
466 
467  for(int i = 0; i < m_planes[n]->GetExpSize(); ++i)
468  {
469  StdRegions::StdExpansionSharedPtr exp = m_planes[n]->GetExp(i);
470  Array<OneD, NekDouble> phys(exp->GetTotPoints());
471  exp->BwdTrans(m_planes[n]->GetCoeffs()+m_planes[n]->GetCoeff_Offset(i),
472  phys);
473  err = exp->L2(phys);
474  energy[n/2] += err*err;
475 
476  exp = m_planes[n+1]->GetExp(i);
477  exp->BwdTrans(m_planes[n+1]->GetCoeffs()+m_planes[n+1]->GetCoeff_Offset(i),
478  phys);
479  err = exp->L2(phys);
480  energy[n/2] += err*err;
481  }
482 
483  m_comm->GetRowComm()->AllReduce(energy[n/2], LibUtilities::ReduceSum);
484  energy[n/2] /= 2.0*area;
485  }
486 
487  return energy;
488  }
489  } //end of namespace
490 } //end of namespace
491 
492 
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
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:1775
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:1867
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
static Array< OneD, NekDouble > NullNekDouble1DArray
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
void GenExpList3DHomogeneous1D(const SpatialDomains::ExpansionMap &expansions)
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
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2)
void SetCoeffPhys(void)
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
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
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
Array< OneD, ExpListSharedPtr > m_planes
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
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
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:877
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 . ...
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
Definition: ExpList2D.h:60
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:211
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
Describes the specification for a Basis.
Definition: Basis.h:50
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171