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 using namespace std;
41 
42 namespace Nektar
43 {
44  namespace MultiRegions
45  {
46  // Forward declaration for typedefs
47  ExpList3DHomogeneous1D::ExpList3DHomogeneous1D():
49  {
51  }
52 
55  const LibUtilities::BasisKey &HomoBasis,
56  const NekDouble lhom,
57  const bool useFFT,
58  const bool dealiasing):
59  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
60  {
62  }
63 
64  // Constructor for ExpList3DHomogeneous1D to act as a Explist2D field
67  const LibUtilities::BasisKey &HomoBasis,
68  const NekDouble lhom,
69  const bool useFFT,
70  const bool dealiasing,
72  const std::string &var):
73  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
74  {
76 
77  GenExpList3DHomogeneous1D(graph2D->GetExpansions(var));
78  }
79 
80  // Constructor for ExpList3DHomogeneous1D to act as a Explist2D field
83  const LibUtilities::BasisKey &HomoBasis,
84  const NekDouble lhom,
85  const bool useFFT,
86  const bool dealiasing,
87  const SpatialDomains::ExpansionMap &expansions):
88  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
89  {
91 
92  GenExpList3DHomogeneous1D(expansions);
93  }
94 
96  {
97  int n,j,nel;
98  bool False = false;
99  ExpList2DSharedPtr plane_zero;
100 
101  // note that nzplanes can be larger than nzmodes
102  m_planes[0] = plane_zero = MemoryManager<ExpList2D>::AllocateSharedPtr(m_session, expansions, False);
103 
105  nel = m_planes[0]->GetExpSize();
106 
107  for(j = 0; j < nel; ++j)
108  {
109  (*m_exp).push_back(m_planes[0]->GetExp(j));
110  }
111 
112  for(n = 1; n < m_homogeneousBasis->GetNumPoints(); ++n)
113  {
115  for(j = 0; j < nel; ++j)
116  {
117  (*m_exp).push_back((*m_exp)[j]);
118  }
119  }
120 
121  // Setup Default optimisation information.
122  nel = GetExpSize();
125 
126  SetCoeffPhys();
127  }
128 
129 
130  /**
131  * @param In ExpList3DHomogeneous1D object to copy.
132  */
133  ExpList3DHomogeneous1D::ExpList3DHomogeneous1D(const ExpList3DHomogeneous1D &In, bool DeclarePlanesSetCoeffPhys):
135  {
137 
138  if(DeclarePlanesSetCoeffPhys)
139  {
140  bool False = false;
141  ExpList2DSharedPtr zero_plane = boost::dynamic_pointer_cast<ExpList2D> (In.m_planes[0]);
142 
143  for(int n = 0; n < m_planes.num_elements(); ++n)
144  {
146  }
147 
148  SetCoeffPhys();
149  }
150  }
151 
152  /**
153  * @param In ExpList3DHomogeneous1D object to copy.
154  * @param eIDs Id of elements that should be copied.
155  */
157  const std::vector<unsigned int> &eIDs,
158  bool DeclarePlanesSetCoeffPhys):
159  ExpListHomogeneous1D(In, eIDs)
160  {
162 
163  if(DeclarePlanesSetCoeffPhys)
164  {
165  bool False = false;
166  std::vector<unsigned int> eIDsPlane;
167  int nel = eIDs.size()/m_planes.num_elements();
168 
169  for (int i = 0; i < nel; ++i)
170  {
171  eIDsPlane.push_back(eIDs[i]);
172  }
173 
174  ExpList2DSharedPtr zero_plane_old =
175  boost::dynamic_pointer_cast<ExpList2D> (In.m_planes[0]);
176 
177  ExpList2DSharedPtr zero_plane =
178  MemoryManager<ExpList2D>::AllocateSharedPtr(*(zero_plane_old), eIDsPlane);
179 
180  for(int n = 0; n < m_planes.num_elements(); ++n)
181  {
183  }
184 
185  SetCoeffPhys();
186  }
187  }
188 
189  /**
190  * Destructor
191  */
193  {
194  }
195 
197  {
198  int i,n,cnt;
199  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
200  int npoints_per_plane = m_planes[0]->GetTotPoints();
201 
202  int nzplanes = m_planes.num_elements();
203 
204  // Set total coefficients and points
205  m_ncoeffs = ncoeffs_per_plane*nzplanes;
206  m_npoints = npoints_per_plane*nzplanes;
207 
210 
211  int nel = m_planes[0]->GetExpSize();
212  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
213  m_phys_offset = Array<OneD,int>(nel*nzplanes);
214  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
215  Array<OneD, NekDouble> tmparray;
216 
217  for(cnt = n = 0; n < nzplanes; ++n)
218  {
219  m_planes[n]->SetCoeffsArray(tmparray= m_coeffs + ncoeffs_per_plane*n);
220  m_planes[n]->SetPhysArray(tmparray = m_phys + npoints_per_plane*n);
221 
222  for(i = 0; i < nel; ++i)
223  {
224  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i) + n*ncoeffs_per_plane;
225  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i) + n*npoints_per_plane;
226  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i) + n*nel;
227  }
228  }
229  }
230 
235  {
236  int n;
237  Array<OneD, NekDouble> tmp_xc;
238  int nzplanes = m_planes.num_elements();
239  int npoints = GetTotPoints(eid);
240 
241  (*m_exp)[eid]->GetCoords(xc0,xc1);
242 
243  // Fill z-direction
245  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
246 
247  for(n = 0; n < m_planes.num_elements(); n++)
248  {
249  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
250  }
251 
252  Array<OneD, NekDouble> z(nzplanes);
253 
254  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
255  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
256 
257  for(n = 0; n < nzplanes; ++n)
258  {
259  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
260  if(n)
261  {
262  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
263  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
264  }
265  }
266  }
267 
268  /**
269  * The operation calls the 2D plane coordinates through the
270  * function ExpList#GetCoords and then evaluated the third
271  * coordinate using the member \a m_lhom
272  *
273  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
274  * will be stored in this array.
275  *
276  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
277  * will be stored in this array.
278  *
279  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
280  * will be stored in this array. This
281  * coordinate is evaluated using the
282  * predefined value \a m_lhom
283  */
287  {
288  int n;
289  Array<OneD, NekDouble> tmp_xc;
290  int nzplanes = m_planes.num_elements();
291  int npoints = m_planes[0]->GetTotPoints();
292 
293  m_planes[0]->GetCoords(xc0,xc1);
294 
295  // Fill z-direction
297 
298  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
299 
300  for(n = 0; n < m_planes.num_elements(); n++)
301  {
302  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
303  }
304 
305  Array<OneD, NekDouble> z(nzplanes);
306 
307  Vmath::Smul(nzplanes,m_lhom/2.0,local_pts,1,z,1);
308  Vmath::Sadd(nzplanes,m_lhom/2.0,z,1,z,1);
309 
310  for(n = 0; n < nzplanes; ++n)
311  {
312  Vmath::Fill(npoints,z[n],tmp_xc = xc2 + npoints*n,1);
313  if(n)
314  {
315  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
316  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
317  }
318  }
319  }
320 
321  void ExpList3DHomogeneous1D::v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
322  {
323  ASSERTL0(expansion == -1, "Multi-zone output not supported for homogeneous expansions.");
324 
325  const int nPtsPlane = m_planes[0]->GetNpoints();
326  const int nElmt = m_planes[0]->GetExpSize();
327  const int nPlanes = m_planes.num_elements();
328 
329  int cnt = 0;
330  int cnt2 = 0;
331  for (int i = 0; i < nElmt; ++i)
332  {
333  const int np0 = (*m_exp)[i]->GetNumPoints(0);
334  const int np1 = (*m_exp)[i]->GetNumPoints(1);
335 
336  for (int n = 1; n < nPlanes; ++n)
337  {
338  const int o1 = (n-1) * nPtsPlane;
339  const int o2 = n * nPtsPlane;
340  for (int j = 1; j < np1; ++j)
341  {
342  for(int k = 1; k < np0; ++k)
343  {
344  outfile << cnt + (j-1)*np0 + (k-1) + o1 + 1 << " ";
345  outfile << cnt + (j-1)*np0 + (k-1) + o2 + 1 << " ";
346  outfile << cnt + (j-1)*np0 + k + o2 + 1 << " ";
347  outfile << cnt + (j-1)*np0 + k + o1 + 1 << " ";
348  outfile << cnt + j *np0 + (k-1) + o1 + 1 << " ";
349  outfile << cnt + j *np0 + (k-1) + o2 + 1 << " ";
350  outfile << cnt + j *np0 + k + o2 + 1 << " ";
351  outfile << cnt + j *np0 + k + o1 + 1 << endl;
352  cnt2++;
353  }
354  }
355  }
356 
357  cnt += np0*np1;
358  }
359  }
360 
361 
363  std::ostream &outfile,
364  int expansion,
365  int istrip)
366  {
367  // If there is only one plane (e.g. HalfMode), we write a 2D plane.
368  if (m_planes.num_elements() == 1)
369  {
370  m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
371  return;
372  }
373 
374  int i,j,k;
375  int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
376  int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
377  int nq2 = m_planes.num_elements();
378  int ntot = nq0*nq1*nq2;
379  int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
380 
381  Array<OneD,NekDouble> coords[3];
382  coords[0] = Array<OneD,NekDouble>(ntot);
383  coords[1] = Array<OneD,NekDouble>(ntot);
384  coords[2] = Array<OneD,NekDouble>(ntot);
385  GetCoords(expansion,coords[0],coords[1],coords[2]);
386 
387  NekDouble DistStrip;
388  m_session->LoadParameter("DistStrip", DistStrip, 0);
389  // Reset the z-coords for homostrips
390  for(int i = 0; i < ntot; i++)
391  {
392  coords[2][i] += istrip*DistStrip;
393  }
394 
395  outfile << " <Piece NumberOfPoints=\""
396  << ntot << "\" NumberOfCells=\""
397  << ntotminus << "\">" << endl;
398  outfile << " <Points>" << endl;
399  outfile << " <DataArray type=\"Float32\" "
400  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
401  outfile << " ";
402  for (i = 0; i < ntot; ++i)
403  {
404  for (j = 0; j < 3; ++j)
405  {
406  outfile << coords[j][i] << " ";
407  }
408  outfile << endl;
409  }
410  outfile << endl;
411  outfile << " </DataArray>" << endl;
412  outfile << " </Points>" << endl;
413  outfile << " <Cells>" << endl;
414  outfile << " <DataArray type=\"Int32\" "
415  << "Name=\"connectivity\" format=\"ascii\">" << endl;
416  for (i = 0; i < nq0-1; ++i)
417  {
418  for (j = 0; j < nq1-1; ++j)
419  {
420  for (k = 0; k < nq2-1; ++k)
421  {
422  outfile << k*nq0*nq1 + j*nq0 + i << " "
423  << k*nq0*nq1 + j*nq0 + i + 1 << " "
424  << k*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
425  << k*nq0*nq1 + (j+1)*nq0 + i << " "
426  << (k+1)*nq0*nq1 + j*nq0 + i << " "
427  << (k+1)*nq0*nq1 + j*nq0 + i + 1 << " "
428  << (k+1)*nq0*nq1 + (j+1)*nq0 + i + 1 << " "
429  << (k+1)*nq0*nq1 + (j+1)*nq0 + i << endl;
430  }
431  }
432  }
433  outfile << endl;
434  outfile << " </DataArray>" << endl;
435  outfile << " <DataArray type=\"Int32\" "
436  << "Name=\"offsets\" format=\"ascii\">" << endl;
437  for (i = 0; i < ntotminus; ++i)
438  {
439  outfile << i*8+8 << " ";
440  }
441  outfile << endl;
442  outfile << " </DataArray>" << endl;
443  outfile << " <DataArray type=\"UInt8\" "
444  << "Name=\"types\" format=\"ascii\">" << endl;
445  for (i = 0; i < ntotminus; ++i)
446  {
447  outfile << "12 ";
448  }
449  outfile << endl;
450  outfile << " </DataArray>" << endl;
451  outfile << " </Cells>" << endl;
452  outfile << " <PointData>" << endl;
453  }
454 
456  const Array<OneD, const NekDouble> &inarray,
457  const Array<OneD, const NekDouble> &soln)
458  {
459  int cnt = 0;
460  NekDouble errL2, err = 0.0;
462  Array<OneD, NekDouble> local_w(m_planes.num_elements());
463 
464  for(int n = 0; n < m_planes.num_elements(); n++)
465  {
466  local_w[n] = w[m_transposition->GetPlaneID(n)];
467  }
468 
469  if (soln == NullNekDouble1DArray)
470  {
471  for(int n = 0; n < m_planes.num_elements(); ++n)
472  {
473  errL2 = m_planes[n]->L2(inarray + cnt);
474  cnt += m_planes[n]->GetTotPoints();
475  err += errL2*errL2*local_w[n]*m_lhom*0.5;
476  }
477  }
478  else
479  {
480  for(int n = 0; n < m_planes.num_elements(); ++n)
481  {
482  errL2 = m_planes[n]->L2(inarray + cnt, soln + cnt);
483  cnt += m_planes[n]->GetTotPoints();
484  err += errL2*errL2*local_w[n]*m_lhom*0.5;
485  }
486  }
487  m_comm->GetColumnComm()->AllReduce(err, LibUtilities::ReduceSum);
488 
489  return sqrt(err);
490  }
491 
493  {
494  Array<OneD, NekDouble> energy(m_planes.num_elements()/2);
495  NekDouble area = 0.0;
496 
497  // Calculate total area of elements.
498  for (int n = 0; n < m_planes[0]->GetExpSize(); ++n)
499  {
500  Array<OneD, NekDouble> inarray(m_planes[0]->GetExp(n)->GetTotPoints(), 1.0);
501  area += m_planes[0]->GetExp(n)->Integral(inarray);
502  }
503 
504  m_comm->GetRowComm()->AllReduce(area, LibUtilities::ReduceSum);
505 
506  // Calculate L2 norm of real/imaginary planes.
507  for (int n = 0; n < m_planes.num_elements(); n += 2)
508  {
509  NekDouble err;
510 
511  energy[n/2] = 0;
512 
513  for(int i = 0; i < m_planes[n]->GetExpSize(); ++i)
514  {
515  StdRegions::StdExpansionSharedPtr exp = m_planes[n]->GetExp(i);
516  Array<OneD, NekDouble> phys(exp->GetTotPoints());
517  exp->BwdTrans(m_planes[n]->GetCoeffs()+m_planes[n]->GetCoeff_Offset(i),
518  phys);
519  err = exp->L2(phys);
520  energy[n/2] += err*err;
521 
522  exp = m_planes[n+1]->GetExp(i);
523  exp->BwdTrans(m_planes[n+1]->GetCoeffs()+m_planes[n+1]->GetCoeff_Offset(i),
524  phys);
525  err = exp->L2(phys);
526  energy[n/2] += err*err;
527  }
528 
529  m_comm->GetRowComm()->AllReduce(energy[n/2], LibUtilities::ReduceSum);
530  energy[n/2] /= 2.0*area;
531  }
532 
533  return energy;
534  }
535  } //end of namespace
536 } //end of namespace
537 
538 
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
STL namespace.
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:253
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