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