Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
54  const LibUtilities::BasisKey &HomoBasis_y,
55  const LibUtilities::BasisKey &HomoBasis_z,
56  const NekDouble lhom_y,
57  const NekDouble lhom_z,
58  const bool useFFT,
59  const bool dealiasing):
60  ExpListHomogeneous2D(pSession,HomoBasis_y,HomoBasis_z,lhom_y,lhom_z,useFFT,dealiasing)
61  {
63  }
64 
65  // Constructor for ExpList3DHomogeneous2D to act as a Explist1D field
67  const LibUtilities::BasisKey &HomoBasis_y,
68  const LibUtilities::BasisKey &HomoBasis_z,
69  const NekDouble lhom_y,
70  const NekDouble lhom_z,
71  const bool useFFT,
72  const bool dealiasing,
73  const SpatialDomains::MeshGraphSharedPtr &graph1D):
74  ExpListHomogeneous2D(pSession,HomoBasis_y,HomoBasis_z,lhom_y,lhom_z,useFFT,dealiasing)
75  {
77 
78  int n,j,nel;
79  bool False = false;
80  ExpList1DSharedPtr line_zero;
81 
82  //
84  False);
85 
87  nel = m_lines[0]->GetExpSize();
88 
89  for(j = 0; j < nel; ++j)
90  {
91  (*m_exp).push_back(m_lines[0]->GetExp(j));
92  }
93 
94  int ny = m_homogeneousBasis_y->GetNumPoints();
95  int nz = m_homogeneousBasis_z->GetNumPoints();
96 
97  for(n = 1; n < (ny*nz); ++n)
98  {
100  for(j = 0; j < nel; ++j)
101  {
102  (*m_exp).push_back((*m_exp)[j]);
103  }
104  }
105 
106  // Setup Default optimisation information.
107  nel = GetExpSize();
110 
111  SetCoeffPhys();
112  }
113 
114 
115  /**
116  * @param In ExpList3DHomogeneous2D object to copy.
117  */
118  ExpList3DHomogeneous2D::ExpList3DHomogeneous2D(const ExpList3DHomogeneous2D &In, const bool DeclareLinesSetCoeffPhys):
120  {
122 
123  if(DeclareLinesSetCoeffPhys)
124  {
125  bool False = false;
126  ExpList1DSharedPtr zero_line = boost::dynamic_pointer_cast<ExpList1D> (In.m_lines[0]);
127 
128  for(int n = 0; n < m_lines.num_elements(); ++n)
129  {
131  }
132 
133  SetCoeffPhys();
134  }
135  }
136 
137  /**
138  *
139  */
141  const std::vector<unsigned int> &eIDs,
142  const bool DeclareLinesSetCoeffPhys):
144  {
146 
147  if(DeclareLinesSetCoeffPhys)
148  {
149  bool False = false;
150  std::vector<unsigned int> eIDsLine;
151  int nel = eIDs.size()/m_lines.num_elements();
152 
153  for (int i = 0; i < nel; ++i)
154  {
155  eIDsLine.push_back(eIDs[i]);
156  }
157 
158  ExpList1DSharedPtr zero_line_old =
159  boost::dynamic_pointer_cast<ExpList1D> (In.m_lines[0]);
160 
161  ExpList1DSharedPtr zero_line =
162  MemoryManager<ExpList1D>::AllocateSharedPtr(*(zero_line_old), eIDsLine);
163 
164  for(int n = 0; n < m_lines.num_elements(); ++n)
165  {
167  }
168 
169  SetCoeffPhys();
170  }
171  }
172 
173  /**
174  * Destructor
175  */
177  {
178  }
179 
181  {
182  int i,n,cnt;
183  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
184  int npoints_per_line = m_lines[0]->GetTotPoints();
185 
186  int nyzlines = m_lines.num_elements();
187 
188  // Set total coefficients and points
189  m_ncoeffs = ncoeffs_per_line*nyzlines;
190  m_npoints = npoints_per_line*nyzlines;
191 
194 
195  int nel = m_lines[0]->GetExpSize();
196  m_coeff_offset = Array<OneD,int>(nel*nyzlines);
197  m_phys_offset = Array<OneD,int>(nel*nyzlines);
198  m_offset_elmt_id = Array<OneD,int>(nel*nyzlines);
199  Array<OneD, NekDouble> tmparray;
200 
201  for(cnt = n = 0; n < nyzlines; ++n)
202  {
203  m_lines[n]->SetCoeffsArray(tmparray= m_coeffs + ncoeffs_per_line*n);
204  m_lines[n]->SetPhysArray(tmparray = m_phys + npoints_per_line*n);
205 
206  for(i = 0; i < nel; ++i)
207  {
208  m_coeff_offset[cnt] = m_lines[n]->GetCoeff_Offset(i) + n*ncoeffs_per_line;
209  m_phys_offset[cnt] = m_lines[n]->GetPhys_Offset(i) + n*npoints_per_line;
210  m_offset_elmt_id[cnt++] = m_lines[n]->GetOffset_Elmt_Id(i) + n*nel;
211  }
212  }
213  }
214 
219  {
220  int n,m,j;
221  Array<OneD, NekDouble> tmp_xc;
222  int nylines = m_homogeneousBasis_y->GetNumPoints();
223  int nzlines = m_homogeneousBasis_z->GetNumPoints();
224 
225  int npoints = GetTotPoints(eid);
226 
227  // Fill x-y-z-direction
230 
231  Array<OneD, NekDouble> x(npoints);
232  Array<OneD, NekDouble> y(nylines);
233  Array<OneD, NekDouble> z(nzlines);
234 
235  Vmath::Smul(nylines,m_lhom_y/2.0,pts_y,1,y,1);
236  Vmath::Sadd(nylines,m_lhom_y/2.0,y,1,y,1);
237 
238  Vmath::Smul(nzlines,m_lhom_z/2.0,pts_z,1,z,1);
239  Vmath::Sadd(nzlines,m_lhom_z/2.0,z,1,z,1);
240 
241  (*m_exp)[eid]->GetCoords(x);
242 
243 
244  for(m = 0; m < nzlines; ++m)
245  {
246  for(j = 0; j < nylines; ++j)
247  {
248  for(n = 0; n < npoints; ++n)
249  {
250  Vmath::Fill(1,x[n],tmp_xc = xc0 + n +(j*npoints) + (m*npoints*nylines), 1);
251  Vmath::Fill(1,y[j],tmp_xc = xc1 + n +(j*npoints) + (m*npoints*nylines), 1);
252  Vmath::Fill(1,z[m],tmp_xc = xc2 + n +(j*npoints) + (m*npoints*nylines), 1);
253  }
254  }
255  }
256  }
257 
258  /**
259  * The operation calls the 2D plane coordinates through the
260  * function ExpList#GetCoords and then evaluated the third
261  * coordinate using the member \a m_lhom
262  *
263  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
264  * will be stored in this array.
265  *
266  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
267  * will be stored in this array.
268  *
269  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
270  * will be stored in this array. This
271  * coordinate is evaluated using the
272  * predefined value \a m_lhom
273  */
277  {
278  int n,m,j;
279  Array<OneD, NekDouble> tmp_xc;
280  int npoints = m_lines[0]->GetTotPoints();
281 
282  int nylines = m_homogeneousBasis_y->GetNumPoints();
283  int nzlines = m_homogeneousBasis_z->GetNumPoints();
284 
285  // Fill z-direction
288 
289  Array<OneD, NekDouble> x(npoints);
290  Array<OneD, NekDouble> y(nylines);
291  Array<OneD, NekDouble> z(nzlines);
292 
293  m_lines[0]->GetCoords(x);
294 
295  Vmath::Smul(nylines,m_lhom_y/2.0,pts_y,1,y,1);
296  Vmath::Sadd(nylines,m_lhom_y/2.0,y,1,y,1);
297 
298  Vmath::Smul(nzlines,m_lhom_z/2.0,pts_z,1,z,1);
299  Vmath::Sadd(nzlines,m_lhom_z/2.0,z,1,z,1);
300 
301  for(m = 0; m < nzlines; ++m)
302  {
303  for(j = 0; j < nylines; ++j)
304  {
305  for(n = 0; n < npoints; ++n)
306  {
307  Vmath::Fill(1,x[n],tmp_xc = xc0 + n +(j*npoints) + (m*npoints*nylines), 1);
308  Vmath::Fill(1,y[j],tmp_xc = xc1 + n +(j*npoints) + (m*npoints*nylines), 1);
309  Vmath::Fill(1,z[m],tmp_xc = xc2 + n +(j*npoints) + (m*npoints*nylines), 1);
310  }
311  }
312  }
313  }
314 
315 
316  /**
317  * Write Tecplot Files Zone
318  * @param outfile Output file name.
319  * @param expansion Expansion that is considered
320  */
321  void ExpList3DHomogeneous2D::v_WriteTecplotZone(std::ostream &outfile, int expansion)
322  {
323  int i,j;
324 
325  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
326  int nquad1 = m_homogeneousBasis_y->GetNumPoints();
327  int nquad2 = m_homogeneousBasis_z->GetNumPoints();
328 
329  Array<OneD,NekDouble> coords[3];
330 
331  coords[0] = Array<OneD,NekDouble>(3*nquad0*nquad1*nquad2);
332  coords[1] = coords[0] + nquad0*nquad1*nquad2;
333  coords[2] = coords[1] + nquad0*nquad1*nquad2;
334 
335  GetCoords(expansion,coords[0],coords[1],coords[2]);
336 
337  outfile << "Zone, I=" << nquad0 << ", J=" << nquad1 <<",K="
338  << nquad2 << ", F=Block" << std::endl;
339 
340  for(j = 0; j < 3; ++j)
341  {
342  for(i = 0; i < nquad0*nquad1*nquad2; ++i)
343  {
344  outfile << coords[j][i] << " ";
345  }
346  outfile << std::endl;
347  }
348  }
349 
350 
351  void ExpList3DHomogeneous2D::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int)
352  {
353  int i,j,k;
354  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
355  int nquad1 = m_homogeneousBasis_y->GetNumPoints();
356  int nquad2 = m_homogeneousBasis_z->GetNumPoints();
357  int ntot = nquad0*nquad1*nquad2;
358  int ntotminus = (nquad0-1)*(nquad1-1)*(nquad2-1);
359 
360  Array<OneD,NekDouble> coords[3];
361  coords[0] = Array<OneD,NekDouble>(ntot);
362  coords[1] = Array<OneD,NekDouble>(ntot);
363  coords[2] = Array<OneD,NekDouble>(ntot);
364  GetCoords(expansion,coords[0],coords[1],coords[2]);
365 
366  outfile << " <Piece NumberOfPoints=\""
367  << ntot << "\" NumberOfCells=\""
368  << ntotminus << "\">" << endl;
369  outfile << " <Points>" << endl;
370  outfile << " <DataArray type=\"Float64\" "
371  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
372  outfile << " ";
373  for (i = 0; i < ntot; ++i)
374  {
375  for (j = 0; j < 3; ++j)
376  {
377  outfile << coords[j][i] << " ";
378  }
379  outfile << endl;
380  }
381  outfile << endl;
382  outfile << " </DataArray>" << endl;
383  outfile << " </Points>" << endl;
384  outfile << " <Cells>" << endl;
385  outfile << " <DataArray type=\"Int32\" "
386  << "Name=\"connectivity\" format=\"ascii\">" << endl;
387  for (i = 0; i < nquad0-1; ++i)
388  {
389  for (j = 0; j < nquad1-1; ++j)
390  {
391  for (k = 0; k < nquad2-1; ++k)
392  {
393  outfile << k*nquad0*nquad1 + j*nquad0 + i << " "
394  << k*nquad0*nquad1 + j*nquad0 + i + 1 << " "
395  << k*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
396  << k*nquad0*nquad1 + (j+1)*nquad0 + i << " "
397  << (k+1)*nquad0*nquad1 + j*nquad0 + i << " "
398  << (k+1)*nquad0*nquad1 + j*nquad0 + i + 1 << " "
399  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i + 1 << " "
400  << (k+1)*nquad0*nquad1 + (j+1)*nquad0 + i << " " << endl;
401  }
402  }
403  }
404  outfile << endl;
405  outfile << " </DataArray>" << endl;
406  outfile << " <DataArray type=\"Int32\" "
407  << "Name=\"offsets\" format=\"ascii\">" << endl;
408  for (i = 0; i < ntotminus; ++i)
409  {
410  outfile << i*8+8 << " ";
411  }
412  outfile << endl;
413  outfile << " </DataArray>" << endl;
414  outfile << " <DataArray type=\"UInt8\" "
415  << "Name=\"types\" format=\"ascii\">" << endl;
416  for (i = 0; i < ntotminus; ++i)
417  {
418  outfile << "12 ";
419  }
420  outfile << endl;
421  outfile << " </DataArray>" << endl;
422  outfile << " </Cells>" << endl;
423  outfile << " <PointData>" << endl;
424  }
425 
426 
428  const Array<OneD, const NekDouble> &inarray,
429  const Array<OneD, const NekDouble> &soln)
430  {
431  int cnt = 0;
432  NekDouble errL2,err = 0.0;
435 
436  int nylines = m_homogeneousBasis_y->GetNumPoints();
437  int nzlines = m_homogeneousBasis_z->GetNumPoints();
438 
439  for(int m = 0; m < nzlines; ++m)
440  {
441  for(int n = 0; n < nylines; ++n)
442  {
443  errL2 = m_lines[n+(m*nylines)]->L2(inarray + cnt, soln + cnt);
444  cnt += m_lines[n+(m*nylines)]->GetTotPoints();
445  err += errL2*errL2*w_y[n]*m_lhom_y*0.5*w_z[m]*m_lhom_z*0.5;
446  }
447  }
448 
449  return sqrt(err);
450  }
451  } //end of namespace
452 } //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:1001
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.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:956
NekDouble m_lhom_z
Width of homogeneous direction z.
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
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:988
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: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
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:991
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: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
NekDouble m_lhom_y
Width of homogeneous direction y.
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:910
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:301
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:253
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:477
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