Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Field.hpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Field.hpp
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: Field converter module base classes.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <boost/shared_ptr.hpp>
37 
42 
43 #include <MultiRegions/ExpList.h>
50 
56 
57 
58 using namespace std;
59 
60 namespace Nektar
61 {
62 namespace Utilities
63 {
64 
65 enum PtsType{
70 };
71 
72 struct FieldPts
73 {
74  FieldPts(void): m_ptsDim(0),
75  m_ptype(ePtsFile) {}
76 
77  int m_ptsDim;
78  int m_nFields;
79  Array<OneD, Array<OneD, NekDouble> > m_pts;
81  vector<int> m_npts;
82  vector<std::string> m_fields;
83 
84 
85  // Interpolate field_id (which is the id after the coordinates)
86  void Interp1DPts(
87  const NekDouble coord,
88  Array<OneD, NekDouble > &intfields)
89  {
90  // currently assume first field is coordinate
91  WARNINGL1(m_ptsDim == 1,
92  "Assumed only one coordinate given taking first coordinate "
93  "for interpolation");
94  int npts = m_pts[0].num_elements();
95  int i;
96 
97  for(i = 0; i < npts-1; ++i)
98  {
99  if((m_pts[0][i] <= coord) && (coord <= m_pts[0][i+1]))
100  {
101  NekDouble pdiff = m_pts[0][i+1]-m_pts[0][i];
102 
103  if(npts <= 2)
104  {
105  // linear interpolation
106  for(int j = 0; j < m_nFields; ++j)
107  {
108  intfields[j] = m_pts[m_ptsDim+j][i]
109  * (m_pts[0][i+1] - coord) / pdiff
110  + m_pts[m_ptsDim+j][i+1]
111  * (coord - m_pts[0][i]) / pdiff;
112  }
113  }
114  else // quadratic interpolation
115  {
116  if(i < npts-2)
117  { // forwards stencil
118  NekDouble pdiff2 = m_pts[0][i+2] - m_pts[0][i+1];
119 
120  NekDouble h1 = (m_pts[0][i+1]-coord)
121  * (m_pts[0][i+2] - coord)
122  / (pdiff * (pdiff+pdiff2));
123  NekDouble h2 = (coord-m_pts[0][i])
124  * (m_pts[0][i+2] - coord)
125  / (pdiff * pdiff2);
126  NekDouble h3 = (coord-m_pts[0][i])
127  * (coord - m_pts[0][i+1])
128  / ((pdiff + pdiff2) * pdiff2);
129  for(int j = 0; j < m_nFields; ++j)
130  {
131  intfields[j] = m_pts[m_ptsDim+j][i] * h1
132  + m_pts[m_ptsDim+j][i+1] * h2
133  + m_pts[m_ptsDim+j][i+2] * h3;
134  }
135  }
136  else
137  { // backwards stencil
138  NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i-1];
139 
140  NekDouble h1 = (m_pts[0][i+1]-coord)
141  * (coord - m_pts[0][i-1])
142  / (pdiff * pdiff2);
143  NekDouble h2 = (coord - m_pts[0][i])
144  * (coord - m_pts[0][i-1])
145  / (pdiff * (pdiff + pdiff2));
146  NekDouble h3 = (m_pts[0][i]-coord)
147  * (m_pts[0][i+1] - coord)
148  / ((pdiff + pdiff2) * pdiff);
149  for(int j = 0; j < m_nFields; ++j)
150  {
151  intfields[j] = m_pts[m_ptsDim+j][i] * h1
152  + m_pts[m_ptsDim+j][i+1] * h2
153  + m_pts[m_ptsDim+j][i-1] * h3;
154  }
155  }
156  }
157  break;
158  }
159  }
160  ASSERTL0(i != npts-1, "Failed to find coordinate " +
161  boost::lexical_cast<string>(coord) +
162  " within provided input points");
163  };
164 
165 };
166 
167 typedef boost::shared_ptr<FieldPts> FieldPtsSharedPtr;
169 
170 struct Field {
171  Field() : m_verbose(false),
172  m_declareExpansionAsContField(false),
173  m_declareExpansionAsDisContField(false),
174  m_writeBndFld(false),
175  m_fieldPts(NullFieldPts){}
176 
177  ~Field()
178  {
179  if (m_session)
180  {
181  m_session->Finalise();
182  }
183  }
184  bool m_verbose;
185  vector<LibUtilities::FieldDefinitionsSharedPtr> m_fielddef;
186  vector<vector<double> > m_data;
187  vector<MultiRegions::ExpListSharedPtr> m_exp;
188 
191 
196  map<string, vector<string> > m_inputfiles;
197 
199  vector<unsigned int> m_bndRegionsToWrite;
201 
203 
204 
205  MultiRegions::ExpListSharedPtr SetUpFirstExpList(int NumHomogeneousDir,
206  bool fldfilegiven = false)
207  {
208 
210 
211  // Set up expansion list
212  int expdim = m_graph->GetMeshDimension();
213 
214  bool useFFT = false;
215  bool dealiasing = false;
216 
217  switch (expdim)
218  {
219  case 1:
220  {
221  ASSERTL0(NumHomogeneousDir <= 2,
222  "Quasi-3D approach is only set up for 1 or 2 "
223  "homogeneous directions");
224 
225  if (NumHomogeneousDir == 1)
226  {
228 
229  // Define Homogeneous expansion
230  int nplanes;
231  NekDouble ly;
233 
234  if(fldfilegiven)
235  {
236  nplanes = m_fielddef[0]->m_numModes[1];
237  ly = m_fielddef[0]->m_homogeneousLengths[0];
238  btype = m_fielddef[0]->m_basis[1];
239  }
240  else
241  {
242  m_session->LoadParameter("HomModesZ", nplanes);
243  m_session->LoadParameter("LY",ly);
244  btype = LibUtilities::eFourier;
245  }
246 
247  // Choose points to be at evenly spaced points at
248  // nplanes points
250  Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
251 
252  const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
253 
254 
255 
256  if(m_declareExpansionAsContField||
257  m_declareExpansionAsDisContField)
258  {
259  ASSERTL0(false,"ContField2DHomogeneous1D or "
260  "DisContField2DHomogenenous1D has "
261  "not been implemented");
262  }
263 
264  Exp2DH1 = MemoryManager<MultiRegions::
266  AllocateSharedPtr(m_session, Bkey, ly,
267  useFFT, dealiasing,
268  m_graph);
269  exp = Exp2DH1;
270  }
271  else if (NumHomogeneousDir == 2)
272  {
274 
275  int nylines,nzlines;
276  NekDouble ly,lz;
277  LibUtilities::BasisType btype1,btype2;
278 
279  if(fldfilegiven)
280  {
281  nylines = m_fielddef[0]->m_numModes[1];
282  nzlines = m_fielddef[0]->m_numModes[2];
283  ly = m_fielddef[0]->m_homogeneousLengths[0];
284  lz = m_fielddef[0]->m_homogeneousLengths[1];
285  btype1 = m_fielddef[0]->m_basis[1];
286  btype2 = m_fielddef[0]->m_basis[2];
287  }
288  else
289  {
290  m_session->LoadParameter("HomModesY", nylines);
291  m_session->LoadParameter("HomModesZ", nzlines);
292  m_session->LoadParameter("LY",ly);
293  m_session->LoadParameter("LZ",lz);
294  btype1 = LibUtilities::eFourier;
295  btype2 = LibUtilities::eFourier;
296  }
297 
298  // Choose points to be at evenly spaced points at
299  // nplanes points
301  PkeyY(nylines, LibUtilities::ePolyEvenlySpaced);
302  const LibUtilities::BasisKey BkeyY(btype1, nylines, PkeyY);
303 
305  PkeyZ(nzlines, LibUtilities::ePolyEvenlySpaced);
306  const LibUtilities::BasisKey BkeyZ(btype2, nzlines, PkeyZ);
307 
308  if(m_declareExpansionAsContField)
309  {
310  Exp3DH2 = MemoryManager<MultiRegions::
312  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
313  ly, lz, useFFT, dealiasing,
314  m_graph,
315  m_session->GetVariable(0));
316  }
317  else if(m_declareExpansionAsDisContField)
318  {
319  Exp3DH2 = MemoryManager<MultiRegions::
321  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
322  ly, lz, useFFT, dealiasing,
323  m_graph,
324  m_session->GetVariable(0));
325  }
326  else
327  {
328  Exp3DH2 = MemoryManager<MultiRegions::
330  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
331  ly, lz, useFFT, dealiasing,
332  m_graph);
333  }
334 
335  exp = Exp3DH2;
336  }
337  else
338  {
340 
341  if(m_declareExpansionAsContField)
342  {
344  ::AllocateSharedPtr(m_session, m_graph,
345  m_session->GetVariable(0));
346  }
347  else if(m_declareExpansionAsDisContField)
348  {
350  ::AllocateSharedPtr(m_session, m_graph,
351  m_session->GetVariable(0));
352  }
353  else
354  {
356  ::AllocateSharedPtr(m_session, m_graph);
357  }
358 
359  exp = Exp1D;
360  }
361  }
362  break;
363  case 2:
364  {
365  ASSERTL0(NumHomogeneousDir <= 1,
366  "NumHomogeneousDir is only set up for 1");
367 
368  if (NumHomogeneousDir == 1)
369  {
371 
372  // Define Homogeneous expansion
373  int nplanes;
374  NekDouble lz;
376 
377  if(fldfilegiven)
378  {
379  nplanes = m_fielddef[0]->m_numModes[2];
380  lz = m_fielddef[0]->m_homogeneousLengths[0];
381  btype = m_fielddef[0]->m_basis[2];
382  }
383  else
384  {
385  m_session->LoadParameter("HomModesZ", nplanes);
386  m_session->LoadParameter("LZ",lz);
387  btype = LibUtilities::eFourier;
388  }
389 
390  // Choose points to be at evenly spaced points at
391  // nplanes points
393  Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
394 
395  const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
396 
397  if(m_declareExpansionAsContField)
398  {
399  Exp3DH1 = MemoryManager<MultiRegions::
401  AllocateSharedPtr(m_session, Bkey, lz, useFFT,
402  dealiasing, m_graph,
403  m_session->GetVariable(0));
404  }
405  else if (m_declareExpansionAsDisContField)
406  {
407  Exp3DH1 = MemoryManager<MultiRegions::
409  AllocateSharedPtr(m_session,
410  Bkey, lz, useFFT,
411  dealiasing, m_graph,
412  m_session->GetVariable(0));
413  }
414  else
415  {
416  Exp3DH1 = MemoryManager<MultiRegions::
418  AllocateSharedPtr(m_session, Bkey, lz, useFFT,
419  dealiasing, m_graph);
420  }
421  exp = Exp3DH1;
422  }
423  else
424  {
426 
427  if(m_declareExpansionAsContField)
428  {
430  ::AllocateSharedPtr(m_session,m_graph,
431  m_session->GetVariable(0));
432  }
433  else if(m_declareExpansionAsDisContField)
434  {
436  ::AllocateSharedPtr(m_session,m_graph,
437  m_session->GetVariable(0));
438  }
439  else
440  {
442  ::AllocateSharedPtr(m_session,m_graph);
443  }
444 
445  exp = Exp2D;
446  }
447  }
448  break;
449  case 3:
450  {
452 
453  if(m_declareExpansionAsContField)
454  {
456  ::AllocateSharedPtr(m_session,m_graph,
457  m_session->GetVariable(0));
458  }
459  else if(m_declareExpansionAsDisContField)
460  {
462  ::AllocateSharedPtr(m_session,m_graph,
463  m_session->GetVariable(0));
464  }
465  else
466  {
468  ::AllocateSharedPtr(m_session, m_graph);
469  }
470 
471  exp = Exp3D;
472  }
473  break;
474  default:
475  ASSERTL0(false, "Expansion dimension not recognised");
476  break;
477  }
478 
479  return exp;
480  };
481 
482  MultiRegions::ExpListSharedPtr AppendExpList(int NumHomogeneousDir,
483  string var = "DefaultVar",
484  bool NewField = false)
485  {
487  switch (m_graph->GetMeshDimension())
488  {
489  case 1:
490  {
491  if (NumHomogeneousDir == 1)
492  {
493  ASSERTL0(m_declareExpansionAsContField ||
494  m_declareExpansionAsDisContField,
495  "ContField2DHomogeneous1D or "
496  "DisContField2DHomogenenous1D has not been "
497  "implemented");
498 
500  boost::dynamic_pointer_cast<MultiRegions::
501  ExpList2DHomogeneous1D>(m_exp[0]);
502 
505  AllocateSharedPtr(*tmp2);
506 
507  }
508  else if (NumHomogeneousDir == 2)
509  {
510  if(m_declareExpansionAsContField)
511  {
513  boost::dynamic_pointer_cast<MultiRegions::
514  ContField3DHomogeneous2D>(m_exp[0]);
515 
518  AllocateSharedPtr(*tmp2);
519  }
520  else if(m_declareExpansionAsDisContField)
521  {
523  boost::dynamic_pointer_cast<MultiRegions::
525 
528  AllocateSharedPtr(*tmp2);
529  }
530  else
531  {
533  boost::dynamic_pointer_cast<MultiRegions::
534  ExpList3DHomogeneous2D>(m_exp[0]);
535 
538  AllocateSharedPtr(*tmp2);
539  }
540 
541 
542  }
543  else
544  {
545  if(m_declareExpansionAsContField)
546  {
548  boost::dynamic_pointer_cast<MultiRegions::
549  ContField1D>(m_exp[0]);
550 
552  AllocateSharedPtr(m_session,m_graph,var);
553  }
554  else if(m_declareExpansionAsDisContField)
555  {
557  boost::dynamic_pointer_cast<MultiRegions::
558  DisContField1D>(m_exp[0]);
559 
561  AllocateSharedPtr(m_session,m_graph,var);
562  }
563  else
564  {
566  boost::dynamic_pointer_cast<MultiRegions::
567  ExpList1D>(m_exp[0]);
568 
570  AllocateSharedPtr(*tmp2);
571  }
572 
573  }
574  }
575  break;
576  case 2:
577  {
578  if (NumHomogeneousDir == 1)
579  {
580  if(m_declareExpansionAsContField)
581  {
582  if(NewField)
583  {
584  bool useFFT = false;
585  bool dealiasing = false;
586 
588  ContField3DHomogeneous1D>::AllocateSharedPtr(
589  m_session,
590  m_exp[0]->GetHomogeneousBasis()
591  ->GetBasisKey(),
592  m_exp[0]->GetHomoLen(),
593  useFFT, dealiasing, m_graph, var);
594  }
595  else
596  {
598  boost::dynamic_pointer_cast<MultiRegions::
599  ContField3DHomogeneous1D>(m_exp[0]);
600 
601  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
604  AllocateSharedPtr(*tmp2);
605  }
606  }
607  else if(m_declareExpansionAsDisContField)
608  {
609  if(NewField)
610  {
611  bool useFFT = false;
612  bool dealiasing = false;
613 
615  DisContField3DHomogeneous1D>::AllocateSharedPtr(
616  m_session,
617  m_exp[0]->GetHomogeneousBasis()
618  ->GetBasisKey(),
619  m_exp[0]->GetHomoLen(),
620  useFFT, dealiasing, m_graph,var);
621  }
622  else
623  {
625  boost::dynamic_pointer_cast<MultiRegions::
627  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
628 
631  AllocateSharedPtr(*tmp2);
632  }
633  }
634  else
635  {
636  if(NewField)
637  {
638  bool useFFT = false;
639  bool dealiasing = false;
640 
642  ExpList3DHomogeneous1D>::AllocateSharedPtr(
643  m_session,
644  m_exp[0]->GetHomogeneousBasis()
645  ->GetBasisKey(),
646  m_exp[0]->GetHomoLen(),
647  useFFT, dealiasing, m_graph);
648  }
649  else
650  {
652  boost::dynamic_pointer_cast<MultiRegions::
653  ExpList3DHomogeneous1D>(m_exp[0]);
654  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
655 
658  AllocateSharedPtr(*tmp2);
659  }
660  }
661 
662  }
663  else
664  {
665  if(m_declareExpansionAsContField)
666  {
667  if(NewField)
668  {
670  AllocateSharedPtr(m_session,m_graph,var);
671  }
672  else // call copy constructor
673  {
674 
676  boost::dynamic_pointer_cast<MultiRegions::
677  ContField2D>(m_exp[0]);
678 
680  AllocateSharedPtr(*tmp2,m_graph,var);
681  }
682  }
683  else if(m_declareExpansionAsDisContField)
684  {
685  if(NewField)
686  {
688  AllocateSharedPtr(m_session,m_graph,var);
689  }
690  else // call copy constructor
691  {
693  boost::dynamic_pointer_cast<MultiRegions::
694  DisContField2D>(m_exp[0]);
695 
697  AllocateSharedPtr(*tmp2,m_graph,var);
698  }
699  }
700  else
701  {
703  boost::dynamic_pointer_cast<MultiRegions::
704  ExpList2D>(m_exp[0]);
705 
707  AllocateSharedPtr(*tmp2);
708  }
709  }
710  }
711  break;
712  case 3:
713  {
714  if(m_declareExpansionAsContField)
715  {
716  if(NewField)
717  {
719  AllocateSharedPtr(m_session,m_graph,var);
720  }
721  else
722  {
724  boost::dynamic_pointer_cast<MultiRegions::
725  ContField3D>(m_exp[0]);
726 
728  AllocateSharedPtr(*tmp2,m_graph,var);
729  }
730  }
731  else if(m_declareExpansionAsDisContField)
732  {
733  if(NewField)
734  {
736  AllocateSharedPtr(m_session,m_graph,var);
737  }
738  else
739  {
741  boost::dynamic_pointer_cast<MultiRegions::
742  DisContField3D>(m_exp[0]);
743 
745  AllocateSharedPtr(*tmp2,m_graph,var);
746  }
747  }
748  else
749  {
751  boost::dynamic_pointer_cast<MultiRegions::
752  ExpList3D>(m_exp[0]);
753 
755  AllocateSharedPtr(*tmp2);
756  }
757  }
758  break;
759  default:
760  ASSERTL0(false, "Expansion dimension not recognised");
761  break;
762  }
763 
764  return tmp;
765  }
766 
767 };
768 
769 typedef boost::shared_ptr<Field> FieldSharedPtr;
770 
771 }
772 }
773