Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ExpList3DHomogeneous2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList3DHomogeneous2D.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 2 directions 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  ExpList3DHomogeneous2D::ExpList3DHomogeneous2D():
49  {
51  }
52 
55  const LibUtilities::BasisKey &HomoBasis_y,
56  const LibUtilities::BasisKey &HomoBasis_z,
57  const NekDouble lhom_y,
58  const NekDouble lhom_z,
59  const bool useFFT,
60  const bool dealiasing,
61  const Collections::ImplementationType ImpType):
62  ExpListHomogeneous2D(pSession,HomoBasis_y,HomoBasis_z,lhom_y,lhom_z,useFFT,dealiasing)
63  {
65  }
66 
67  // Constructor for ExpList3DHomogeneous2D to act as a Explist1D field
70  const LibUtilities::BasisKey &HomoBasis_y,
71  const LibUtilities::BasisKey &HomoBasis_z,
72  const NekDouble lhom_y,
73  const NekDouble lhom_z,
74  const bool useFFT,
75  const bool dealiasing,
77  const Collections::ImplementationType ImpType):
78  ExpListHomogeneous2D(pSession,HomoBasis_y,HomoBasis_z,lhom_y,lhom_z,useFFT,dealiasing)
79  {
81 
82  int n,j,nel;
83  bool False = false;
84  ExpList1DSharedPtr line_zero;
85 
86  //
87  m_lines[0] = line_zero = MemoryManager<ExpList1D>::
88  AllocateSharedPtr(m_session,graph1D, False,ImpType);
89 
92  nel = m_lines[0]->GetExpSize();
93 
94  for(j = 0; j < nel; ++j)
95  {
96  (*m_exp).push_back(m_lines[0]->GetExp(j));
97  }
98 
99  int ny = m_homogeneousBasis_y->GetNumPoints();
100  int nz = m_homogeneousBasis_z->GetNumPoints();
101 
102  for(n = 1; n < (ny*nz); ++n)
103  {
105  for(j = 0; j < nel; ++j)
106  {
107  (*m_exp).push_back((*m_exp)[j]);
108  }
109  }
110 
111  // Setup Default optimisation information.
112  nel = GetExpSize();
115 
116  SetCoeffPhys();
117  }
118 
119 
120  /**
121  * @param In ExpList3DHomogeneous2D object to copy.
122  */
124  const ExpList3DHomogeneous2D &In,
125  const bool DeclareLinesSetCoeffPhys):
127  {
129 
130  if(DeclareLinesSetCoeffPhys)
131  {
132  bool False = false;
133  ExpList1DSharedPtr zero_line = boost::dynamic_pointer_cast<ExpList1D> (In.m_lines[0]);
134 
135  for(int n = 0; n < m_lines.num_elements(); ++n)
136  {
138  }
139 
140  SetCoeffPhys();
141  }
142  }
143 
144  /**
145  *
146  */
148  const ExpList3DHomogeneous2D &In,
149  const std::vector<unsigned int> &eIDs,
150  const bool DeclareLinesSetCoeffPhys,
151  const Collections::ImplementationType ImpType):
152  ExpListHomogeneous2D(In, eIDs)
153  {
155 
156  if(DeclareLinesSetCoeffPhys)
157  {
158  bool False = false;
159  std::vector<unsigned int> eIDsLine;
160  int nel = eIDs.size()/m_lines.num_elements();
161 
162  for (int i = 0; i < nel; ++i)
163  {
164  eIDsLine.push_back(eIDs[i]);
165  }
166 
167  ExpList1DSharedPtr zero_line_old =
168  boost::dynamic_pointer_cast<ExpList1D> (In.m_lines[0]);
169 
170  ExpList1DSharedPtr zero_line =
171  MemoryManager<ExpList1D>::AllocateSharedPtr(*(zero_line_old), eIDsLine, ImpType);
172 
173  for(int n = 0; n < m_lines.num_elements(); ++n)
174  {
176  }
177 
178  SetCoeffPhys();
179  }
180  }
181 
182  /**
183  * Destructor
184  */
186  {
187  }
188 
190  {
191  int i,n,cnt;
192  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
193  int npoints_per_line = m_lines[0]->GetTotPoints();
194 
195  int nyzlines = m_lines.num_elements();
196 
197  // Set total coefficients and points
198  m_ncoeffs = ncoeffs_per_line*nyzlines;
199  m_npoints = npoints_per_line*nyzlines;
200 
203 
204  int nel = m_lines[0]->GetExpSize();
205  m_coeff_offset = Array<OneD,int>(nel*nyzlines);
206  m_phys_offset = Array<OneD,int>(nel*nyzlines);
207  Array<OneD, NekDouble> tmparray;
208 
209  for(cnt = n = 0; n < nyzlines; ++n)
210  {
211  m_lines[n]->SetCoeffsArray(tmparray= m_coeffs + ncoeffs_per_line*n);
212  m_lines[n]->SetPhysArray(tmparray = m_phys + npoints_per_line*n);
213 
214  for(i = 0; i < nel; ++i)
215  {
216  m_coeff_offset[cnt] = m_lines[n]->GetCoeff_Offset(i) + n*ncoeffs_per_line;
217  m_phys_offset[cnt++] = m_lines[n]->GetPhys_Offset(i) + n*npoints_per_line;
218  }
219  }
220  }
221 
226  {
227  int n,m,j;
228  Array<OneD, NekDouble> tmp_xc;
229  int nylines = m_homogeneousBasis_y->GetNumPoints();
230  int nzlines = m_homogeneousBasis_z->GetNumPoints();
231 
232  int npoints = GetTotPoints(eid);
233 
234  // Fill x-y-z-direction
237 
238  Array<OneD, NekDouble> x(npoints);
239  Array<OneD, NekDouble> y(nylines);
240  Array<OneD, NekDouble> z(nzlines);
241 
242  Vmath::Smul(nylines,m_lhom_y/2.0,pts_y,1,y,1);
243  Vmath::Sadd(nylines,m_lhom_y/2.0,y,1,y,1);
244 
245  Vmath::Smul(nzlines,m_lhom_z/2.0,pts_z,1,z,1);
246  Vmath::Sadd(nzlines,m_lhom_z/2.0,z,1,z,1);
247 
248  (*m_exp)[eid]->GetCoords(x);
249 
250 
251  for(m = 0; m < nzlines; ++m)
252  {
253  for(j = 0; j < nylines; ++j)
254  {
255  for(n = 0; n < npoints; ++n)
256  {
257  Vmath::Fill(1,x[n],tmp_xc = xc0 + n +(j*npoints) + (m*npoints*nylines), 1);
258  Vmath::Fill(1,y[j],tmp_xc = xc1 + n +(j*npoints) + (m*npoints*nylines), 1);
259  Vmath::Fill(1,z[m],tmp_xc = xc2 + n +(j*npoints) + (m*npoints*nylines), 1);
260  }
261  }
262  }
263  }
264 
265  /**
266  * The operation calls the 2D plane coordinates through the
267  * function ExpList#GetCoords and then evaluated the third
268  * coordinate using the member \a m_lhom
269  *
270  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
271  * will be stored in this array.
272  *
273  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
274  * will be stored in this array.
275  *
276  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
277  * will be stored in this array. This
278  * coordinate is evaluated using the
279  * predefined value \a m_lhom
280  */
284  {
285  int n,m,j;
286  Array<OneD, NekDouble> tmp_xc;
287  int npoints = m_lines[0]->GetTotPoints();
288 
289  int nylines = m_homogeneousBasis_y->GetNumPoints();
290  int nzlines = m_homogeneousBasis_z->GetNumPoints();
291 
292  // Fill z-direction
295 
296  Array<OneD, NekDouble> x(npoints);
297  Array<OneD, NekDouble> y(nylines);
298  Array<OneD, NekDouble> z(nzlines);
299 
300  m_lines[0]->GetCoords(x);
301 
302  Vmath::Smul(nylines,m_lhom_y/2.0,pts_y,1,y,1);
303  Vmath::Sadd(nylines,m_lhom_y/2.0,y,1,y,1);
304 
305  Vmath::Smul(nzlines,m_lhom_z/2.0,pts_z,1,z,1);
306  Vmath::Sadd(nzlines,m_lhom_z/2.0,z,1,z,1);
307 
308  for(m = 0; m < nzlines; ++m)
309  {
310  for(j = 0; j < nylines; ++j)
311  {
312  for(n = 0; n < npoints; ++n)
313  {
314  Vmath::Fill(1,x[n],tmp_xc = xc0 + n +(j*npoints) + (m*npoints*nylines), 1);
315  Vmath::Fill(1,y[j],tmp_xc = xc1 + n +(j*npoints) + (m*npoints*nylines), 1);
316  Vmath::Fill(1,z[m],tmp_xc = xc2 + n +(j*npoints) + (m*npoints*nylines), 1);
317  }
318  }
319  }
320  }
321 
322 
323  /**
324  * Write Tecplot Files Zone
325  * @param outfile Output file name.
326  * @param expansion Expansion that is considered
327  */
328  void ExpList3DHomogeneous2D::v_WriteTecplotZone(std::ostream &outfile, int expansion)
329  {
330  int i,j;
331 
332  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
333  int nquad1 = m_homogeneousBasis_y->GetNumPoints();
334  int nquad2 = m_homogeneousBasis_z->GetNumPoints();
335 
336  Array<OneD,NekDouble> coords[3];
337 
338  coords[0] = Array<OneD,NekDouble>(3*nquad0*nquad1*nquad2);
339  coords[1] = coords[0] + nquad0*nquad1*nquad2;
340  coords[2] = coords[1] + nquad0*nquad1*nquad2;
341 
342  GetCoords(expansion,coords[0],coords[1],coords[2]);
343 
344  outfile << "Zone, I=" << nquad0 << ", J=" << nquad1 <<",K="
345  << nquad2 << ", F=Block" << std::endl;
346 
347  for(j = 0; j < 3; ++j)
348  {
349  for(i = 0; i < nquad0*nquad1*nquad2; ++i)
350  {
351  outfile << coords[j][i] << " ";
352  }
353  outfile << std::endl;
354  }
355  }
356 
357 
358  void ExpList3DHomogeneous2D::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int)
359  {
360  int i,j,k;
361  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
362  int nquad1 = m_homogeneousBasis_y->GetNumPoints();
363  int nquad2 = m_homogeneousBasis_z->GetNumPoints();
364  int ntot = nquad0*nquad1*nquad2;
365  int ntotminus = (nquad0-1)*(nquad1-1)*(nquad2-1);
366 
367  Array<OneD,NekDouble> coords[3];
368  coords[0] = Array<OneD,NekDouble>(ntot);
369  coords[1] = Array<OneD,NekDouble>(ntot);
370  coords[2] = Array<OneD,NekDouble>(ntot);
371  GetCoords(expansion,coords[0],coords[1],coords[2]);
372 
373  outfile << " <Piece NumberOfPoints=\""
374  << ntot << "\" NumberOfCells=\""
375  << ntotminus << "\">" << endl;
376  outfile << " <Points>" << endl;
377  outfile << " <DataArray type=\"Float64\" "
378  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
379  outfile << " ";
380  for (i = 0; i < ntot; ++i)
381  {
382  for (j = 0; j < 3; ++j)
383  {
384  outfile << coords[j][i] << " ";
385  }
386  outfile << endl;
387  }
388  outfile << endl;
389  outfile << " </DataArray>" << endl;
390  outfile << " </Points>" << endl;
391  outfile << " <Cells>" << endl;
392  outfile << " <DataArray type=\"Int32\" "
393  << "Name=\"connectivity\" format=\"ascii\">" << endl;
394  for (i = 0; i < nquad0-1; ++i)
395  {
396  for (j = 0; j < nquad1-1; ++j)
397  {
398  for (k = 0; k < nquad2-1; ++k)
399  {
400  outfile << k*nquad0*nquad1 + j*nquad0 + i << " "
401  << k*nquad0*nquad1 + j*nquad0 + i + 1 << " "
402  << k*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
403  << k*nquad0*nquad1 + (j+1)*nquad0 + i << " "
404  << (k+1)*nquad0*nquad1 + j*nquad0 + i << " "
405  << (k+1)*nquad0*nquad1 + j*nquad0 + i + 1 << " "
406  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
407  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i << " " << endl;
408  }
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*8+8 << " ";
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 << "12 ";
426  }
427  outfile << endl;
428  outfile << " </DataArray>" << endl;
429  outfile << " </Cells>" << endl;
430  outfile << " <PointData>" << endl;
431  }
432 
433 
435  const Array<OneD, const NekDouble> &inarray,
436  const Array<OneD, const NekDouble> &soln)
437  {
438  int cnt = 0;
439  NekDouble errL2,err = 0.0;
442 
443  int nylines = m_homogeneousBasis_y->GetNumPoints();
444  int nzlines = m_homogeneousBasis_z->GetNumPoints();
445 
446  for(int m = 0; m < nzlines; ++m)
447  {
448  for(int n = 0; n < nylines; ++n)
449  {
450  errL2 = m_lines[n+(m*nylines)]->L2(inarray + cnt, soln + cnt);
451  cnt += m_lines[n+(m*nylines)]->GetTotPoints();
452  err += errL2*errL2*w_y[n]*m_lhom_y*0.5*w_z[m]*m_lhom_z*0.5;
453  }
454  }
455 
456  return sqrt(err);
457  }
458  } //end of namespace
459 } //end of namespace
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
Definition: ExpList.h:1052
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:2067
STL namespace.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1015
NekDouble m_lhom_z
Width of homogeneous direction z.
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:2046
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, ExpListSharedPtr > m_lines
Vector of ExpList, will be filled with ExpList1D.
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1047
LibUtilities::BasisSharedPtr m_homogeneousBasis_z
Base expansion in z direction.
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:1535
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
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_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1050
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2)
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:976
NekDouble m_lhom_y
Width of homogeneous direction y.
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:969
double NekDouble
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
LibUtilities::BasisSharedPtr m_homogeneousBasis_y
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
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 . ...
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:277
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error with respect to soln of the global spectral/hp element approximat...
Definition: ExpList.h:528
Abstraction of a one-dimensional multi-elemental expansion which is merely a collection of local expa...
Describes the specification for a Basis.
Definition: Basis.h:50