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
66 {
72 };
73 
74 struct FieldPts
75 {
76  FieldPts(void): m_ptsDim(0),
77  m_ptype(ePtsFile) {}
78 
79  int m_ptsDim;
80  int m_nFields;
81  Array<OneD, Array<OneD, NekDouble> > m_pts;
83  vector<int> m_npts;
84  vector<std::string> m_fields;
85  vector<Array<OneD, int> > m_ptsConn;
86 
87  // Interpolate field_id (which is the id after the coordinates)
88  void Interp1DPts(const NekDouble coord,
89  Array<OneD, NekDouble > &intfields)
90  {
91  // currently assume first field is coordinate
92  WARNINGL1(m_ptsDim == 1,
93  "Assumed only one coordinate given taking first coordinate "
94  "for interpolation");
95  int npts = m_pts[0].num_elements();
96  int i;
97 
98  for(i = 0; i < npts-1; ++i)
99  {
100  if((m_pts[0][i] <= coord) && (coord <= m_pts[0][i+1]))
101  {
102  NekDouble pdiff = m_pts[0][i+1]-m_pts[0][i];
103 
104  if(npts <= 2)
105  {
106  // linear interpolation
107  for(int j = 0; j < m_nFields; ++j)
108  {
109  intfields[j] = m_pts[m_ptsDim+j][i]
110  * (m_pts[0][i+1] - coord) / pdiff
111  + m_pts[m_ptsDim+j][i+1]
112  * (coord - m_pts[0][i]) / pdiff;
113  }
114  }
115  else // quadratic interpolation
116  {
117  if(i < npts-2)
118  { // forwards stencil
119  NekDouble pdiff2 = m_pts[0][i+2] - m_pts[0][i+1];
120 
121  NekDouble h1 = (m_pts[0][i+1]-coord)
122  * (m_pts[0][i+2] - coord)
123  / (pdiff * (pdiff+pdiff2));
124  NekDouble h2 = (coord-m_pts[0][i])
125  * (m_pts[0][i+2] - coord)
126  / (pdiff * pdiff2);
127  NekDouble h3 = (coord-m_pts[0][i])
128  * (coord - m_pts[0][i+1])
129  / ((pdiff + pdiff2) * pdiff2);
130  for(int j = 0; j < m_nFields; ++j)
131  {
132  intfields[j] = m_pts[m_ptsDim+j][i] * h1
133  + m_pts[m_ptsDim+j][i+1] * h2
134  + m_pts[m_ptsDim+j][i+2] * h3;
135  }
136  }
137  else
138  { // backwards stencil
139  NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i-1];
140 
141  NekDouble h1 = (m_pts[0][i+1]-coord)
142  * (coord - m_pts[0][i-1])
143  / (pdiff * pdiff2);
144  NekDouble h2 = (coord - m_pts[0][i])
145  * (coord - m_pts[0][i-1])
146  / (pdiff * (pdiff + pdiff2));
147  NekDouble h3 = (m_pts[0][i]-coord)
148  * (m_pts[0][i+1] - coord)
149  / ((pdiff + pdiff2) * pdiff);
150  for(int j = 0; j < m_nFields; ++j)
151  {
152  intfields[j] = m_pts[m_ptsDim+j][i] * h1
153  + m_pts[m_ptsDim+j][i+1] * h2
154  + m_pts[m_ptsDim+j][i-1] * h3;
155  }
156  }
157  }
158  break;
159  }
160  }
161 
162  ASSERTL0(i != npts-1, "Failed to find coordinate " +
163  boost::lexical_cast<string>(coord) +
164  " within provided input points");
165  };
166 
167 };
168 
169 typedef boost::shared_ptr<FieldPts> FieldPtsSharedPtr;
171 
172 struct Field {
173  Field() : m_verbose(false),
174  m_declareExpansionAsContField(false),
175  m_declareExpansionAsDisContField(false),
176  m_writeBndFld(false),
177  m_fieldPts(NullFieldPts){}
178 
179  ~Field()
180  {
181  if (m_session)
182  {
183  m_session->Finalise();
184  }
185  }
186  bool m_verbose;
187  vector<LibUtilities::FieldDefinitionsSharedPtr> m_fielddef;
188  vector<vector<double> > m_data;
189  vector<MultiRegions::ExpListSharedPtr> m_exp;
190 
193 
198  map<string, vector<string> > m_inputfiles;
199 
201  vector<unsigned int> m_bndRegionsToWrite;
203 
206 
207 
208  MultiRegions::ExpListSharedPtr SetUpFirstExpList(int NumHomogeneousDir,
209  bool fldfilegiven = false)
210  {
211 
213 
214  // Set up expansion list
215  int expdim = m_graph->GetMeshDimension();
216 
217  bool useFFT = false;
218  bool dealiasing = false;
219 
220  switch (expdim)
221  {
222  case 1:
223  {
224  ASSERTL0(NumHomogeneousDir <= 2,
225  "Quasi-3D approach is only set up for 1 or 2 "
226  "homogeneous directions");
227 
228  if (NumHomogeneousDir == 1)
229  {
231 
232  // Define Homogeneous expansion
233  int nplanes;
234  NekDouble ly;
236 
237  if(fldfilegiven)
238  {
239  nplanes = m_fielddef[0]->m_numModes[1];
240  ly = m_fielddef[0]->m_homogeneousLengths[0];
241  btype = m_fielddef[0]->m_basis[1];
242  }
243  else
244  {
245  m_session->LoadParameter("HomModesZ", nplanes);
246  m_session->LoadParameter("LY",ly);
247  btype = LibUtilities::eFourier;
248  }
249 
250  // Choose points to be at evenly spaced points at
251  // nplanes points
253  Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
254 
255  const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
256 
257 
258 
259  if(m_declareExpansionAsContField||
260  m_declareExpansionAsDisContField)
261  {
262  ASSERTL0(false,"ContField2DHomogeneous1D or "
263  "DisContField2DHomogenenous1D has "
264  "not been implemented");
265  }
266 
267  Exp2DH1 = MemoryManager<MultiRegions::
269  AllocateSharedPtr(m_session, Bkey, ly,
270  useFFT, dealiasing,
271  m_graph);
272  exp = Exp2DH1;
273  }
274  else if (NumHomogeneousDir == 2)
275  {
277 
278  int nylines,nzlines;
279  NekDouble ly,lz;
280  LibUtilities::BasisType btype1,btype2;
281 
282  if(fldfilegiven)
283  {
284  nylines = m_fielddef[0]->m_numModes[1];
285  nzlines = m_fielddef[0]->m_numModes[2];
286  ly = m_fielddef[0]->m_homogeneousLengths[0];
287  lz = m_fielddef[0]->m_homogeneousLengths[1];
288  btype1 = m_fielddef[0]->m_basis[1];
289  btype2 = m_fielddef[0]->m_basis[2];
290  }
291  else
292  {
293  m_session->LoadParameter("HomModesY", nylines);
294  m_session->LoadParameter("HomModesZ", nzlines);
295  m_session->LoadParameter("LY",ly);
296  m_session->LoadParameter("LZ",lz);
297  btype1 = LibUtilities::eFourier;
298  btype2 = LibUtilities::eFourier;
299  }
300 
301  // Choose points to be at evenly spaced points at
302  // nplanes points
304  PkeyY(nylines, LibUtilities::ePolyEvenlySpaced);
305  const LibUtilities::BasisKey BkeyY(btype1, nylines, PkeyY);
306 
308  PkeyZ(nzlines, LibUtilities::ePolyEvenlySpaced);
309  const LibUtilities::BasisKey BkeyZ(btype2, nzlines, PkeyZ);
310 
311  if(m_declareExpansionAsContField)
312  {
313  Exp3DH2 = MemoryManager<MultiRegions::
315  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
316  ly, lz, useFFT, dealiasing,
317  m_graph,
318  m_session->GetVariable(0));
319  }
320  else if(m_declareExpansionAsDisContField)
321  {
322  Exp3DH2 = MemoryManager<MultiRegions::
324  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
325  ly, lz, useFFT, dealiasing,
326  m_graph,
327  m_session->GetVariable(0));
328  }
329  else
330  {
331  Exp3DH2 = MemoryManager<MultiRegions::
333  AllocateSharedPtr(m_session, BkeyY, BkeyZ,
334  ly, lz, useFFT, dealiasing,
335  m_graph);
336  }
337 
338  exp = Exp3DH2;
339  }
340  else
341  {
343 
344  if(m_declareExpansionAsContField)
345  {
347  ::AllocateSharedPtr(m_session, m_graph,
348  m_session->GetVariable(0));
349  }
350  else if(m_declareExpansionAsDisContField)
351  {
353  ::AllocateSharedPtr(m_session, m_graph,
354  m_session->GetVariable(0));
355  }
356  else
357  {
359  ::AllocateSharedPtr(m_session, m_graph);
360  }
361 
362  exp = Exp1D;
363  }
364  }
365  break;
366  case 2:
367  {
368  ASSERTL0(NumHomogeneousDir <= 1,
369  "NumHomogeneousDir is only set up for 1");
370 
371  if (NumHomogeneousDir == 1)
372  {
374 
375  // Define Homogeneous expansion
376  int nplanes;
377  NekDouble lz;
379 
380  if(fldfilegiven)
381  {
382  nplanes = m_fielddef[0]->m_numModes[2];
383  lz = m_fielddef[0]->m_homogeneousLengths[0];
384  btype = m_fielddef[0]->m_basis[2];
385  }
386  else
387  {
388  m_session->LoadParameter("HomModesZ", nplanes);
389  m_session->LoadParameter("LZ",lz);
390  btype = LibUtilities::eFourier;
391  }
392 
393  // Choose points to be at evenly spaced points at
394  // nplanes points
396  Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
397 
398  const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
399 
400  if(m_declareExpansionAsContField)
401  {
402  Exp3DH1 = MemoryManager<MultiRegions::
404  AllocateSharedPtr(m_session, Bkey, lz, useFFT,
405  dealiasing, m_graph,
406  m_session->GetVariable(0));
407  }
408  else if (m_declareExpansionAsDisContField)
409  {
410  Exp3DH1 = MemoryManager<MultiRegions::
412  AllocateSharedPtr(m_session,
413  Bkey, lz, useFFT,
414  dealiasing, m_graph,
415  m_session->GetVariable(0));
416  }
417  else
418  {
419  Exp3DH1 = MemoryManager<MultiRegions::
421  AllocateSharedPtr(m_session, Bkey, lz, useFFT,
422  dealiasing, m_graph);
423  }
424  exp = Exp3DH1;
425  }
426  else
427  {
429 
430  if(m_declareExpansionAsContField)
431  {
433  ::AllocateSharedPtr(m_session,m_graph,
434  m_session->GetVariable(0));
435  }
436  else if(m_declareExpansionAsDisContField)
437  {
439  ::AllocateSharedPtr(m_session,m_graph,
440  m_session->GetVariable(0));
441  }
442  else
443  {
445  ::AllocateSharedPtr(m_session,m_graph);
446  }
447 
448  exp = Exp2D;
449  }
450  }
451  break;
452  case 3:
453  {
455 
456  if(m_declareExpansionAsContField)
457  {
459  ::AllocateSharedPtr(m_session,m_graph,
460  m_session->GetVariable(0));
461  }
462  else if(m_declareExpansionAsDisContField)
463  {
465  ::AllocateSharedPtr(m_session,m_graph,
466  m_session->GetVariable(0));
467  }
468  else
469  {
471  ::AllocateSharedPtr(m_session, m_graph);
472  }
473 
474  exp = Exp3D;
475  }
476  break;
477  default:
478  ASSERTL0(false, "Expansion dimension not recognised");
479  break;
480  }
481 
482  return exp;
483  };
484 
485  MultiRegions::ExpListSharedPtr AppendExpList(int NumHomogeneousDir,
486  string var = "DefaultVar",
487  bool NewField = false)
488  {
490  switch (m_graph->GetMeshDimension())
491  {
492  case 1:
493  {
494  if (NumHomogeneousDir == 1)
495  {
496  ASSERTL0(m_declareExpansionAsContField ||
497  m_declareExpansionAsDisContField,
498  "ContField2DHomogeneous1D or "
499  "DisContField2DHomogenenous1D has not been "
500  "implemented");
501 
503  boost::dynamic_pointer_cast<MultiRegions::
504  ExpList2DHomogeneous1D>(m_exp[0]);
505 
508  AllocateSharedPtr(*tmp2);
509 
510  }
511  else if (NumHomogeneousDir == 2)
512  {
513  if(m_declareExpansionAsContField)
514  {
516  boost::dynamic_pointer_cast<MultiRegions::
517  ContField3DHomogeneous2D>(m_exp[0]);
518 
521  AllocateSharedPtr(*tmp2);
522  }
523  else if(m_declareExpansionAsDisContField)
524  {
526  boost::dynamic_pointer_cast<MultiRegions::
528 
531  AllocateSharedPtr(*tmp2);
532  }
533  else
534  {
536  boost::dynamic_pointer_cast<MultiRegions::
537  ExpList3DHomogeneous2D>(m_exp[0]);
538 
541  AllocateSharedPtr(*tmp2);
542  }
543 
544 
545  }
546  else
547  {
548  if(m_declareExpansionAsContField)
549  {
551  boost::dynamic_pointer_cast<MultiRegions::
552  ContField1D>(m_exp[0]);
553 
555  AllocateSharedPtr(m_session,m_graph,var);
556  }
557  else if(m_declareExpansionAsDisContField)
558  {
560  boost::dynamic_pointer_cast<MultiRegions::
561  DisContField1D>(m_exp[0]);
562 
564  AllocateSharedPtr(m_session,m_graph,var);
565  }
566  else
567  {
569  boost::dynamic_pointer_cast<MultiRegions::
570  ExpList1D>(m_exp[0]);
571 
573  AllocateSharedPtr(*tmp2);
574  }
575 
576  }
577  }
578  break;
579  case 2:
580  {
581  if (NumHomogeneousDir == 1)
582  {
583  if(m_declareExpansionAsContField)
584  {
585  if(NewField)
586  {
587  bool useFFT = false;
588  bool dealiasing = false;
589 
591  ContField3DHomogeneous1D>::AllocateSharedPtr(
592  m_session,
593  m_exp[0]->GetHomogeneousBasis()
594  ->GetBasisKey(),
595  m_exp[0]->GetHomoLen(),
596  useFFT, dealiasing, m_graph, var);
597  }
598  else
599  {
601  boost::dynamic_pointer_cast<MultiRegions::
602  ContField3DHomogeneous1D>(m_exp[0]);
603 
604  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
607  AllocateSharedPtr(*tmp2);
608  }
609  }
610  else if(m_declareExpansionAsDisContField)
611  {
612  if(NewField)
613  {
614  bool useFFT = false;
615  bool dealiasing = false;
616 
618  DisContField3DHomogeneous1D>::AllocateSharedPtr(
619  m_session,
620  m_exp[0]->GetHomogeneousBasis()
621  ->GetBasisKey(),
622  m_exp[0]->GetHomoLen(),
623  useFFT, dealiasing, m_graph,var);
624  }
625  else
626  {
628  boost::dynamic_pointer_cast<MultiRegions::
630  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
631 
634  AllocateSharedPtr(*tmp2);
635  }
636  }
637  else
638  {
639  if(NewField)
640  {
641  bool useFFT = false;
642  bool dealiasing = false;
643 
645  ExpList3DHomogeneous1D>::AllocateSharedPtr(
646  m_session,
647  m_exp[0]->GetHomogeneousBasis()
648  ->GetBasisKey(),
649  m_exp[0]->GetHomoLen(),
650  useFFT, dealiasing, m_graph);
651  }
652  else
653  {
655  boost::dynamic_pointer_cast<MultiRegions::
656  ExpList3DHomogeneous1D>(m_exp[0]);
657  ASSERTL0(tmp2,"Failed to type cast m_exp[0]");
658 
661  AllocateSharedPtr(*tmp2);
662  }
663  }
664 
665  }
666  else
667  {
668  if(m_declareExpansionAsContField)
669  {
670  if(NewField)
671  {
673  AllocateSharedPtr(m_session,m_graph,var);
674  }
675  else // call copy constructor
676  {
677 
679  boost::dynamic_pointer_cast<MultiRegions::
680  ContField2D>(m_exp[0]);
681 
683  AllocateSharedPtr(*tmp2,m_graph,var);
684  }
685  }
686  else if(m_declareExpansionAsDisContField)
687  {
688  if(NewField)
689  {
691  AllocateSharedPtr(m_session,m_graph,var);
692  }
693  else // call copy constructor
694  {
696  boost::dynamic_pointer_cast<MultiRegions::
697  DisContField2D>(m_exp[0]);
698 
700  AllocateSharedPtr(*tmp2,m_graph,var);
701  }
702  }
703  else
704  {
706  boost::dynamic_pointer_cast<MultiRegions::
707  ExpList2D>(m_exp[0]);
708 
710  AllocateSharedPtr(*tmp2);
711  }
712  }
713  }
714  break;
715  case 3:
716  {
717  if(m_declareExpansionAsContField)
718  {
719  if(NewField)
720  {
722  AllocateSharedPtr(m_session,m_graph,var);
723  }
724  else
725  {
727  boost::dynamic_pointer_cast<MultiRegions::
728  ContField3D>(m_exp[0]);
729 
731  AllocateSharedPtr(*tmp2,m_graph,var);
732  }
733  }
734  else if(m_declareExpansionAsDisContField)
735  {
736  if(NewField)
737  {
739  AllocateSharedPtr(m_session,m_graph,var);
740  }
741  else
742  {
744  boost::dynamic_pointer_cast<MultiRegions::
745  DisContField3D>(m_exp[0]);
746 
748  AllocateSharedPtr(*tmp2,m_graph,var);
749  }
750  }
751  else
752  {
754  boost::dynamic_pointer_cast<MultiRegions::
755  ExpList3D>(m_exp[0]);
756 
758  AllocateSharedPtr(*tmp2);
759  }
760  }
761  break;
762  default:
763  ASSERTL0(false, "Expansion dimension not recognised");
764  break;
765  }
766 
767  return tmp;
768  }
769 
770 };
771 
772 typedef boost::shared_ptr<Field> FieldSharedPtr;
773 
774 }
775 }
776