Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpList2DHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList2DHomogeneous1D.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 ExpList2D which is homogeneous in 1 direction and so
33 // uses much of the functionality from a ExpList1D and its daughters
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #include <MultiRegions/ExpList2DHomogeneous1D.h>
38 #include <MultiRegions/ExpList1D.h>
39 
40 namespace Nektar
41 {
42  namespace MultiRegions
43  {
44  // Forward declaration for typedefs
45  ExpList2DHomogeneous1D::ExpList2DHomogeneous1D():
46  ExpListHomogeneous1D()
47  {
48  }
49 
50  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
51  ExpList2DHomogeneous1D::ExpList2DHomogeneous1D(
52  const LibUtilities::SessionReaderSharedPtr &pSession,
53  const LibUtilities::BasisKey &HomoBasis,
54  const NekDouble lhom,
55  const bool useFFT,
56  const bool dealiasing,
57  const Array<OneD, ExpListSharedPtr> &planes)
58  : ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
59  {
60  int i, n, cnt, nel;
61 
62  ASSERTL1(m_planes.num_elements() == planes.num_elements(),
63  "Size of basis number of points and number"
64  "of planes are not the same");
65 
66  // Set up expansion list with elements from all planes.
67  m_exp = MemoryManager<LocalRegions::ExpansionVector>
68  ::AllocateSharedPtr(
69  planes.num_elements() * planes[0]->GetExpSize());
70 
71  for(cnt = n = 0; n < planes.num_elements(); ++n)
72  {
73  m_planes[n] = planes[n];
74  for (i = 0; i < planes[n]->GetExpSize(); ++i)
75  {
76  (*m_exp)[cnt++] = planes[n]->GetExp(i);
77  }
78  }
79 
80  // Setup Default optimisation information.
81  nel = GetExpSize();
82  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
83  ::AllocateSharedPtr(nel);
84 
85  SetCoeffPhys();
86  }
87 
88  // Constructor for ExpList2DHomogeneous1D to act as a Explist2D field
89  ExpList2DHomogeneous1D::ExpList2DHomogeneous1D(
90  const LibUtilities::SessionReaderSharedPtr &pSession,
91  const LibUtilities::BasisKey &HomoBasis,
92  const NekDouble lhom,
93  const bool useFFT,
94  const bool dealiasing,
95  const SpatialDomains::MeshGraphSharedPtr &graph1D):
96  ExpListHomogeneous1D(pSession,HomoBasis,lhom,useFFT,dealiasing)
97  {
98  int n, j, nel;
99  ExpList1DSharedPtr plane_zero;
100 
101  // note that nzplanes can be larger than nzmodes
102  m_planes[0] = plane_zero = MemoryManager<ExpList1D>::
103  AllocateSharedPtr(m_session, graph1D, false);
104 
105  m_exp = MemoryManager<LocalRegions::ExpansionVector>::AllocateSharedPtr();
106 
107  nel = m_planes[0]->GetExpSize();
108 
109  for (j = 0; j < nel; ++j)
110  {
111  (*m_exp).push_back(m_planes[0]->GetExp(j));
112  }
113 
114  for (n = 1; n < m_planes.num_elements(); ++n)
115  {
116  m_planes[n] = MemoryManager<ExpList1D>::
117  AllocateSharedPtr(*plane_zero, false);
118  for(j = 0; j < nel; ++j)
119  {
120  (*m_exp).push_back((*m_exp)[j]);
121  }
122  }
123 
124  // Setup Default optimisation information.
125  nel = GetExpSize();
126  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
127  ::AllocateSharedPtr(nel);
128 
129  SetCoeffPhys();
130  }
131 
132 
133  /**
134  * @param In ExpList2DHomogeneous1D object to copy.
135  */
136  ExpList2DHomogeneous1D::ExpList2DHomogeneous1D(
137  const ExpList2DHomogeneous1D &In):
138  ExpListHomogeneous1D(In)
139  {
140  ExpList1DSharedPtr zero_plane =
141  boost::dynamic_pointer_cast<ExpList1D> (In.m_planes[0]);
142 
143  for (int n = 0; n < m_planes.num_elements(); ++n)
144  {
145  m_planes[n] = MemoryManager<ExpList1D>::
146  AllocateSharedPtr(*zero_plane, false);
147  }
148 
149  SetCoeffPhys();
150  }
151 
152  /**
153  * Destructor
154  */
155  ExpList2DHomogeneous1D::~ExpList2DHomogeneous1D()
156  {
157  }
158 
159  void ExpList2DHomogeneous1D::SetCoeffPhys(void)
160  {
161  int i, n, cnt;
162  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
163  int npoints_per_plane = m_planes[0]->GetTotPoints();
164 
165  int nzplanes = m_planes.num_elements();
166 
167  // Set total coefficients and points
168  m_ncoeffs = ncoeffs_per_plane*nzplanes;
169  m_npoints = npoints_per_plane*nzplanes;
170 
171  m_coeffs = Array<OneD, NekDouble> (m_ncoeffs);
172  m_phys = Array<OneD, NekDouble> (m_npoints);
173 
174  int nel = m_planes[0]->GetExpSize();
175  m_coeff_offset = Array<OneD,int>(nel*nzplanes);
176  m_phys_offset = Array<OneD,int>(nel*nzplanes);
177  m_offset_elmt_id = Array<OneD,int>(nel*nzplanes);
178  Array<OneD, NekDouble> tmparray;
179 
180  for (cnt = n = 0; n < nzplanes; ++n)
181  {
182  m_planes[n]->SetCoeffsArray(
183  tmparray= m_coeffs + ncoeffs_per_plane*n);
184  m_planes[n]->SetPhysArray(
185  tmparray = m_phys + npoints_per_plane*n);
186 
187  for(i = 0; i < nel; ++i)
188  {
189  m_coeff_offset[cnt] = m_planes[n]->GetCoeff_Offset(i)
190  + n*ncoeffs_per_plane;
191  m_phys_offset[cnt] = m_planes[n]->GetPhys_Offset(i)
192  + n*npoints_per_plane;
193  m_offset_elmt_id[cnt++] = m_planes[n]->GetOffset_Elmt_Id(i)
194  + n*nel;
195  }
196  }
197  }
198 
199  void ExpList2DHomogeneous1D::GetCoords(
200  const int eid,
201  Array<OneD, NekDouble> &xc0,
202  Array<OneD, NekDouble> &xc1,
203  Array<OneD, NekDouble> &xc2)
204  {
205  int n, coordim;
206  Array<OneD, NekDouble> tmp_xc,xhom;
207  int nyplanes = m_planes.num_elements();
208  int npoints = GetTotPoints(eid);
209 
210  switch(coordim = GetCoordim(0))
211  {
212  case 1:
213  {
214  (*m_exp)[eid]->GetCoords(xc0);
215  xhom = xc1;
216  }
217  break;
218  case 2:
219  ASSERTL0(xc1.num_elements() != 0,
220  "output coord_1 is not defined");
221  {
222  (*m_exp)[eid]->GetCoords(xc0,xc1);
223  xhom = xc2;
224  }
225  break;
226  default:
227  ASSERTL0(false, "Cannot have coordim being three dimensions"
228  "in a homogeneous field");
229  break;
230  }
231 
232  // Fill homogeneous-direction
233  Array<OneD, const NekDouble> pts = m_homogeneousBasis->GetZ();
234  Array<OneD, NekDouble> z(nyplanes);
235 
236  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
237 
238  for(n = 0; n < m_planes.num_elements(); n++)
239  {
240  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
241  }
242 
243  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
244  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
245 
246  for (n = 0; n < nyplanes; ++n)
247  {
248  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
249  if(n)
250  {
251  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
252  if(coordim == 2)
253  {
254  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
255  }
256  }
257  }
258  }
259 
260  /**
261  * The operation calls the 2D plane coordinates through the
262  * function ExpList#GetCoords and then evaluated the third
263  * coordinate using the member \a m_lhom
264  *
265  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
266  * will be stored in this array.
267  *
268  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
269  * will be stored in this array. This
270  * coordinate might be evaluated using the
271  * predefined value \a m_lhom
272  *
273  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
274  * will be stored in this array. This
275  * coordinate is evaluated using the
276  * predefined value \a m_lhom
277  */
278  void ExpList2DHomogeneous1D::v_GetCoords(
279  Array<OneD, NekDouble> &xc0,
280  Array<OneD, NekDouble> &xc1,
281  Array<OneD, NekDouble> &xc2)
282  {
283  int n,coordim;
284  Array<OneD, NekDouble> tmp_xc, xhom;
285  int nyplanes = m_planes.num_elements();
286  int npoints = m_planes[0]->GetTotPoints();
287 
288  m_planes[0]->GetCoords(xc0,xc1);
289 
290  if((coordim = GetCoordim(0)) == 1)
291  {
292  xhom = xc1;
293  }
294  else
295  {
296  xhom = xc2;
297  }
298 
299  // Fill z-direction
300  Array<OneD, const NekDouble> pts = m_homogeneousBasis->GetZ();
301  Array<OneD, NekDouble> z(nyplanes);
302  Array<OneD, NekDouble> local_pts(m_planes.num_elements());
303 
304  for(n = 0; n < m_planes.num_elements(); n++)
305  {
306  local_pts[n] = pts[m_transposition->GetPlaneID(n)];
307  }
308 
309  Vmath::Smul(nyplanes,m_lhom/2.0,local_pts,1,z,1);
310  Vmath::Sadd(nyplanes,m_lhom/2.0,z,1,z,1);
311 
312  for(n = 0; n < nyplanes; ++n)
313  {
314  Vmath::Fill(npoints,z[n],tmp_xc = xhom + npoints*n,1);
315  if(n)
316  {
317  Vmath::Vcopy(npoints,xc0,1,tmp_xc = xc0+npoints*n,1);
318  if(coordim == 2)
319  {
320  Vmath::Vcopy(npoints,xc1,1,tmp_xc = xc1+npoints*n,1);
321  }
322  }
323  }
324  }
325 
326 
327  /**
328  * Write Tecplot Files Zone
329  * @param outfile Output file name.
330  * @param expansion Expansion that is considered
331  */
332  void ExpList2DHomogeneous1D::v_WriteTecplotZone(
333  std::ostream &outfile,
334  int expansion)
335  {
336  int i, j;
337 
338  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
339  int nquad1 = m_planes.num_elements();
340 
341  Array<OneD,NekDouble> coords[3];
342 
343  coords[0] = Array<OneD,NekDouble>(3*nquad0*nquad1);
344  coords[1] = coords[0] + nquad0*nquad1;
345  coords[2] = coords[1] + nquad0*nquad1;
346 
347  GetCoords(expansion,coords[0],coords[1],coords[2]);
348 
349  outfile << "Zone, I=" << nquad0 << ", J=" << nquad1
350  << ", F=Block" << std::endl;
351 
352  for(j = 0; j < GetCoordim(0)+1; ++j)
353  {
354  for(i = 0; i < nquad0*nquad1; ++i)
355  {
356  outfile << coords[j][i] << " ";
357  }
358  outfile << std::endl;
359  }
360  }
361 
362 
363  void ExpList2DHomogeneous1D::v_WriteVtkPieceHeader(
364  std::ostream &outfile,
365  int expansion)
366  {
367  int i, j;
368  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
369  int nquad1 = m_planes.num_elements();
370  int ntot = nquad0*nquad1;
371  int ntotminus = (nquad0-1)*(nquad1-1);
372 
373  Array<OneD,NekDouble> coords[3];
374  coords[0] = Array<OneD,NekDouble>(ntot);
375  coords[1] = Array<OneD,NekDouble>(ntot);
376  coords[2] = Array<OneD,NekDouble>(ntot);
377  GetCoords(expansion,coords[0],coords[1],coords[2]);
378 
379  outfile << " <Piece NumberOfPoints=\""
380  << ntot << "\" NumberOfCells=\""
381  << ntotminus << "\">" << endl;
382  outfile << " <Points>" << endl;
383  outfile << " <DataArray type=\"Float32\" "
384  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
385  outfile << " ";
386  for (i = 0; i < ntot; ++i)
387  {
388  for (j = 0; j < 3; ++j)
389  {
390  outfile << coords[j][i] << " ";
391  }
392  outfile << endl;
393  }
394  outfile << endl;
395  outfile << " </DataArray>" << endl;
396  outfile << " </Points>" << endl;
397  outfile << " <Cells>" << endl;
398  outfile << " <DataArray type=\"Int32\" "
399  << "Name=\"connectivity\" format=\"ascii\">" << endl;
400  for (i = 0; i < nquad0-1; ++i)
401  {
402  for (j = 0; j < nquad1-1; ++j)
403  {
404  outfile << j*nquad0 + i << " "
405  << j*nquad0 + i + 1 << " "
406  << (j+1)*nquad0 + i + 1 << " "
407  << (j+1)*nquad0 + i << endl;
408  }
409  }
410  outfile << endl;
411  outfile << " </DataArray>" << endl;
412  outfile << " <DataArray type=\"Int32\" "
413  << "Name=\"offsets\" format=\"ascii\">" << endl;
414  for (i = 0; i < ntotminus; ++i)
415  {
416  outfile << i*4+4 << " ";
417  }
418  outfile << endl;
419  outfile << " </DataArray>" << endl;
420  outfile << " <DataArray type=\"UInt8\" "
421  << "Name=\"types\" format=\"ascii\">" << endl;
422  for (i = 0; i < ntotminus; ++i)
423  {
424  outfile << "9 ";
425  }
426  outfile << endl;
427  outfile << " </DataArray>" << endl;
428  outfile << " </Cells>" << endl;
429  outfile << " <PointData>" << endl;
430  }
431 
432  void ExpList2DHomogeneous1D::v_GetNormals(
433  Array<OneD, Array<OneD, NekDouble> > &normals)
434  {
435  int nPlanes = m_planes.num_elements();
436  int nPtsPlane = m_planes[0]->GetNpoints();
437  int nDim = GetCoordim(0) + 1;
438 
439  ASSERTL1(normals.num_elements() >= nDim,
440  "Output vector does not have sufficient dimensions to"
441  "match coordim");
442  ASSERTL1(normals[0].num_elements() >= nPtsPlane,
443  "Output vector does not have sufficient dimensions to"
444  "match coordim");
445 
446  // Calculate normals from plane 0.
447  m_planes[0]->GetNormals(normals);
448 
449  // Copy remaining planes.
450  for (int i = 0; i < nDim; ++i)
451  {
452  for (int n = 1; n < nPlanes; ++n)
453  {
454  Vmath::Vcopy(nPtsPlane, &normals[i][0], 1,
455  &normals[i][n*nPtsPlane], 1);
456  }
457  }
458  }
459  } //end of namespace
460 } //end of namespace
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2