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