Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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  * @param In ExpList3DHomogeneous1D object to copy.
152  * @param eIDs Id of elements that should be copied.
153  */
155  const std::vector<unsigned int> &eIDs,
156  bool DeclarePlanesSetCoeffPhys):
157  ExpListHomogeneous1D(In, eIDs)
158  {
160 
161  if(DeclarePlanesSetCoeffPhys)
162  {
163  bool False = false;
164  std::vector<unsigned int> eIDsPlane;
165  int nel = eIDs.size()/m_planes.num_elements();
166 
167  for (int i = 0; i < nel; ++i)
168  {
169  eIDsPlane.push_back(eIDs[i]);
170  }
171 
172  ExpList2DSharedPtr zero_plane_old =
173  boost::dynamic_pointer_cast<ExpList2D> (In.m_planes[0]);
174 
175  ExpList2DSharedPtr zero_plane =
176  MemoryManager<ExpList2D>::AllocateSharedPtr(*(zero_plane_old), eIDsPlane);
177 
178  for(int n = 0; n < m_planes.num_elements(); ++n)
179  {
181  }
182 
183  SetCoeffPhys();
184  }
185  }
186 
187  /**
188  * Destructor
189  */
191  {
192  }
193 
195  {
196  int i,n,cnt;
197  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
198  int npoints_per_plane = m_planes[0]->GetTotPoints();
199 
200  int nzplanes = m_planes.num_elements();
201 
202  // Set total coefficients and points
203  m_ncoeffs = ncoeffs_per_plane*nzplanes;
204  m_npoints = npoints_per_plane*nzplanes;
205 
208 
209  int nel = m_planes[0]->GetExpSize();
210  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
211  m_phys_offset = Array<OneD,int>(nel*nzplanes);
212  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
213  Array<OneD, NekDouble> tmparray;
214 
215  for(cnt = n = 0; n < nzplanes; ++n)
216  {
217  m_planes[n]->SetCoeffsArray(tmparray= m_coeffs + ncoeffs_per_plane*n);
218  m_planes[n]->SetPhysArray(tmparray = m_phys + npoints_per_plane*n);
219 
220  for(i = 0; i < nel; ++i)
221  {
222  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i) + n*ncoeffs_per_plane;
223  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i) + n*npoints_per_plane;
224  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i) + n*nel;
225  }
226  }
227  }
228 
233  {
234  int n;
235  Array<OneD, NekDouble> tmp_xc;
236  int nzplanes = m_planes.num_elements();
237  int npoints = GetTotPoints(eid);
238 
239  (*m_exp)[eid]->GetCoords(xc0,xc1);
240 
241  // Fill z-direction
243  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
244 
245  for(n = 0; n < m_planes.num_elements(); n++)
246  {
247  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
248  }
249 
250  Array<OneD, NekDouble> z(nzplanes);
251 
252  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
253  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
254 
255  for(n = 0; n < nzplanes; ++n)
256  {
257  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
258  if(n)
259  {
260  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
261  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
262  }
263  }
264  }
265 
266  /**
267  * The operation calls the 2D plane coordinates through the
268  * function ExpList#GetCoords and then evaluated the third
269  * coordinate using the member \a m_lhom
270  *
271  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
272  * will be stored in this array.
273  *
274  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
275  * will be stored in this array.
276  *
277  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
278  * will be stored in this array. This
279  * coordinate is evaluated using the
280  * predefined value \a m_lhom
281  */
285  {
286  int n;
287  Array<OneD, NekDouble> tmp_xc;
288  int nzplanes = m_planes.num_elements();
289  int npoints = m_planes[0]->GetTotPoints();
290 
291  m_planes[0]->GetCoords(xc0,xc1);
292 
293  // Fill z-direction
295 
296  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
297 
298  for(n = 0; n < m_planes.num_elements(); n++)
299  {
300  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
301  }
302 
303  Array<OneD, NekDouble> z(nzplanes);
304 
305  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
306  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
307 
308  for(n = 0; n < nzplanes; ++n)
309  {
310  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
311  if(n)
312  {
313  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
314  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
315  }
316  }
317  }
318 
319  void ExpList3DHomogeneous1D::v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
320  {
321  ASSERTL0(expansion == -1, "Multi-zone output not supported for homogeneous expansions.");
322 
323  const int nPtsPlane = m_planes[0]->GetNpoints();
324  const int nElmt = m_planes[0]->GetExpSize();
325  const int nPlanes = m_planes.num_elements();
326 
327  int cnt = 0;
328  int cnt2 = 0;
329  for (int i = 0; i < nElmt; ++i)
330  {
331  const int np0 = (*m_exp)[i]->GetNumPoints(0);
332  const int np1 = (*m_exp)[i]->GetNumPoints(1);
333 
334  for (int n = 1; n < nPlanes; ++n)
335  {
336  const int o1 = (n-1) * nPtsPlane;
337  const int o2 = n * nPtsPlane;
338  for (int j = 1; j < np1; ++j)
339  {
340  for(int k = 1; k < np0; ++k)
341  {
342  outfile << cnt + (j-1)*np0 + (k-1) + o1 + 1 << " ";
343  outfile << cnt + (j-1)*np0 + (k-1) + o2 + 1 << " ";
344  outfile << cnt + (j-1)*np0 + k + o2 + 1 << " ";
345  outfile << cnt + (j-1)*np0 + k + o1 + 1 << " ";
346  outfile << cnt + j *np0 + (k-1) + o1 + 1 << " ";
347  outfile << cnt + j *np0 + (k-1) + o2 + 1 << " ";
348  outfile << cnt + j *np0 + k + o2 + 1 << " ";
349  outfile << cnt + j *np0 + k + o1 + 1 << endl;
350  cnt2++;
351  }
352  }
353  }
354 
355  cnt += np0*np1;
356  }
357  }
358 
359 
361  std::ostream &outfile,
362  int expansion,
363  int istrip)
364  {
365  // If there is only one plane (e.g. HalfMode), we write a 2D plane.
366  if (m_planes.num_elements() == 1)
367  {
368  m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
369  return;
370  }
371 
372  int i,j,k;
373  int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
374  int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
375  int nq2 = m_planes.num_elements();
376  int ntot = nq0*nq1*nq2;
377  int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
378 
379  Array<OneD,NekDouble> coords[3];
380  coords[0] = Array<OneD,NekDouble>(ntot);
381  coords[1] = Array<OneD,NekDouble>(ntot);
382  coords[2] = Array<OneD,NekDouble>(ntot);
383  GetCoords(expansion,coords[0],coords[1],coords[2]);
384 
385  NekDouble DistStrip;
386  m_session->LoadParameter("DistStrip", DistStrip, 0);
387  // Reset the z-coords for homostrips
388  for(int i = 0; i < ntot; i++)
389  {
390  coords[2][i] += istrip*DistStrip;
391  }
392 
393  outfile << " <Piece NumberOfPoints=\""
394  << ntot << "\" NumberOfCells=\""
395  << ntotminus << "\">" << endl;
396  outfile << " <Points>" << endl;
397  outfile << " <DataArray type=\"Float32\" "
398  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
399  outfile << " ";
400  for (i = 0; i < ntot; ++i)
401  {
402  for (j = 0; j < 3; ++j)
403  {
404  outfile << coords[j][i] << " ";
405  }
406  outfile << endl;
407  }
408  outfile << endl;
409  outfile << " </DataArray>" << endl;
410  outfile << " </Points>" << endl;
411  outfile << " <Cells>" << endl;
412  outfile << " <DataArray type=\"Int32\" "
413  << "Name=\"connectivity\" format=\"ascii\">" << endl;
414  for (i = 0; i < nq0-1; ++i)
415  {
416  for (j = 0; j < nq1-1; ++j)
417  {
418  for (k = 0; k < nq2-1; ++k)
419  {
420  outfile << k*nq0*nq1 + j*nq0 + i << " "
421  << k*nq0*nq1 + j*nq0 + i + 1 << " "
422  << k*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
423  << k*nq0*nq1 + (j+1)*nq0 + i << " "
424  << (k+1)*nq0*nq1 + j*nq0 + i << " "
425  << (k+1)*nq0*nq1 + j*nq0 + i + 1 << " "
426  << (k+1)*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
427  << (k+1)*nq0*nq1 + (j+1)*nq0 + i << endl;
428  }
429  }
430  }
431  outfile << endl;
432  outfile << " </DataArray>" << endl;
433  outfile << " <DataArray type=\"Int32\" "
434  << "Name=\"offsets\" format=\"ascii\">" << endl;
435  for (i = 0; i < ntotminus; ++i)
436  {
437  outfile << i*8+8 << " ";
438  }
439  outfile << endl;
440  outfile << " </DataArray>" << endl;
441  outfile << " <DataArray type=\"UInt8\" "
442  << "Name=\"types\" format=\"ascii\">" << endl;
443  for (i = 0; i < ntotminus; ++i)
444  {
445  outfile << "12 ";
446  }
447  outfile << endl;
448  outfile << " </DataArray>" << endl;
449  outfile << " </Cells>" << endl;
450  outfile << " <PointData>" << endl;
451  }
452 
454  const Array<OneD, const NekDouble> &inarray,
455  const Array<OneD, const NekDouble> &soln)
456  {
457  int cnt = 0;
458  NekDouble errL2, err = 0.0;
460  Array<OneD, NekDouble> local_w(m_planes.num_elements());
461 
462  for(int n = 0; n < m_planes.num_elements(); n++)
463  {
464  local_w[n] = w[m_transposition->GetPlaneID(n)];
465  }
466 
467  if (soln == NullNekDouble1DArray)
468  {
469  for(int n = 0; n < m_planes.num_elements(); ++n)
470  {
471  errL2 = m_planes[n]->L2(inarray + cnt);
472  cnt += m_planes[n]->GetTotPoints();
473  err += errL2*errL2*local_w[n]*m_lhom*0.5;
474  }
475  }
476  else
477  {
478  for(int n = 0; n < m_planes.num_elements(); ++n)
479  {
480  errL2 = m_planes[n]->L2(inarray + cnt, soln + cnt);
481  cnt += m_planes[n]->GetTotPoints();
482  err += errL2*errL2*local_w[n]*m_lhom*0.5;
483  }
484  }
485  m_comm->GetColumnComm()->AllReduce(err, LibUtilities::ReduceSum);
486 
487  return sqrt(err);
488  }
489 
491  {
492  Array<OneD, NekDouble> energy(m_planes.num_elements()/2);
493  NekDouble area = 0.0;
494 
495  // Calculate total area of elements.
496  for (int n = 0; n < m_planes[0]->GetExpSize(); ++n)
497  {
498  Array<OneD, NekDouble> inarray(m_planes[0]->GetExp(n)->GetTotPoints(), 1.0);
499  area += m_planes[0]->GetExp(n)->Integral(inarray);
500  }
501 
502  m_comm->GetRowComm()->AllReduce(area, LibUtilities::ReduceSum);
503 
504  // Calculate L2 norm of real/imaginary planes.
505  for (int n = 0; n < m_planes.num_elements(); n += 2)
506  {
507  NekDouble err;
508 
509  energy[n/2] = 0;
510 
511  for(int i = 0; i < m_planes[n]->GetExpSize(); ++i)
512  {
513  StdRegions::StdExpansionSharedPtr exp = m_planes[n]->GetExp(i);
514  Array<OneD, NekDouble> phys(exp->GetTotPoints());
515  exp->BwdTrans(m_planes[n]->GetCoeffs()+m_planes[n]->GetCoeff_Offset(i),
516  phys);
517  err = exp->L2(phys);
518  energy[n/2] += err*err;
519 
520  exp = m_planes[n+1]->GetExp(i);
521  exp->BwdTrans(m_planes[n+1]->GetCoeffs()+m_planes[n+1]->GetCoeff_Offset(i),
522  phys);
523  err = exp->L2(phys);
524  energy[n/2] += err*err;
525  }
526 
527  m_comm->GetRowComm()->AllReduce(energy[n/2], LibUtilities::ReduceSum);
528  energy[n/2] /= 2.0*area;
529  }
530 
531  return energy;
532  }
533  } //end of namespace
534 } //end of namespace
535 
536 
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:1834
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:1926
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:1001
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:1917
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:956
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:939
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1896
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual void v_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:988
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:977
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1453
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:991
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:917
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
Definition: ExpList.h:999
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:910
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:907
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:251
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
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:174