Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FieldToField.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <iomanip>
4 #include <LibUtilities/BasicUtils/SessionReader.h>
5 #include <LibUtilities/BasicUtils/SharedArray.hpp>
6 
7 
8 #include <MultiRegions/ExpList.h>
9 #include <MultiRegions/ContField1D.h>
10 #include <MultiRegions/ContField2D.h>
11 #include <MultiRegions/ContField3D.h>
12 #include <MultiRegions/ContField3DHomogeneous1D.h>
13 
14 #include <tinyxml.h>
15 
16 #include <boost/math/special_functions/fpclassify.hpp>
17 
18 using namespace Nektar;
19 
20 int main(int argc, char *argv[])
21 {
22  void SetFields(
23  SpatialDomains::MeshGraphSharedPtr &mesh,
24  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef,
25  LibUtilities::SessionReaderSharedPtr &session,
26  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,
27  int nvariables,
28  const vector<std::string> &variables,
29  bool homogeneous);
30 
31  void InterpolateField(
32  Array<OneD, MultiRegions::ExpListSharedPtr> &field0,
33  Array<OneD, MultiRegions::ExpListSharedPtr> &field1,
34  Array<OneD, NekDouble> x1,
35  Array<OneD, NekDouble> y1,
36  Array<OneD, NekDouble> z1 = NullNekDouble1DArray,
37  NekDouble clamp_low = -10000000,
38  NekDouble clamp_up = 10000000);
39 
40  void InterpolateFieldHomo(
41  MultiRegions::ExpListSharedPtr field0,
42  Array<OneD, NekDouble> x1,
43  Array<OneD, NekDouble> y1,
44  MultiRegions::ExpListSharedPtr field1);
45 
46  bool Checkbndmeshes2D(
47  Array<OneD, NekDouble> x0,
48  Array<OneD, NekDouble> y0,
49  Array<OneD, NekDouble> x1,
50  Array<OneD, NekDouble> y1);
51 
52  bool Checkbndmeshes3D(
53  Array<OneD, NekDouble> x0,
54  Array<OneD, NekDouble> y0,
55  Array<OneD, NekDouble> z0,
56  Array<OneD, NekDouble> x1,
57  Array<OneD, NekDouble> y1,
58  Array<OneD, NekDouble> z1);
59 
60  void Writefield(
61  LibUtilities::SessionReaderSharedPtr vSession,
62  const vector<std::string> &variables,
63  string fieldfile,
64  SpatialDomains::MeshGraphSharedPtr &graph,
65  Array<OneD, MultiRegions::ExpListSharedPtr> &outfield);
66 
67  if (argc != 5)
68  {
69  fprintf(stderr, "Usage: ./FieldToField meshfile0 fieldfile0 "
70  "meshfile1 fieldfile1\n");
71  exit(1);
72  }
73 
74  //----------------------------------------------
75  string meshfile0(argv[argc-4]);
76  string fieldfile0(argv[argc-3]);
77 
78  std::vector<std::string> filenames0;
79  filenames0.push_back(meshfile0);
80  filenames0.push_back(meshfile0);
81  LibUtilities::SessionReaderSharedPtr vSession
82  = LibUtilities::SessionReader::CreateInstance(argc, argv, filenames0);
83 
84  // Read in mesh from input file0
85  SpatialDomains::MeshGraphSharedPtr graphShPt
86  = SpatialDomains::MeshGraph::Read(meshfile0);
87  int expdim = graphShPt->GetMeshDimension();
88 
89  //----------------------------------------------
90  // Import fieldfile0.
91  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
92  vector<vector<NekDouble> > fielddata;
93  LibUtilities::Import(fieldfile0, fielddef, fielddata);
94  //----------------------------------------------
95 
96  //read info from fldfile
97  // const std::vector<std::string> variables = fielddef[0]->m_fields;
98  const std::vector<std::string> variables = vSession->GetVariables();
99  bool homo = (fielddef[0]->m_numHomogeneousDir > 0)? true: false;
100 
101  // Define Expansion
102  int nfields;
103  nfields = variables.size();
104  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
105  fields= Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
106  SetFields(graphShPt,fielddef, vSession,fields,nfields,variables,homo);
107  int nq; //pointsper plane
108 
109  //-----------------------------------------------
110  // Copy data from file:fill fields with the fielddata
111  if (fielddef[0]->m_numHomogeneousDir == 1)
112  {
113  nq = fields[0]->GetPlane(1)->GetTotPoints();
114 
115  //THE IM PHYS VALUES ARE WRONG USING bwdTrans !!!
116  for (int j = 0; j < nfields; ++j)
117  {
118  for (int i = 0; i < fielddata.size(); ++i)
119  {
120  fields[j]->ExtractDataToCoeffs(fielddef[i],
121  fielddata[i],
122  fielddef[i]->m_fields[j],
123  fields[j]->UpdateCoeffs());
124  }
125 
126  //bwd plane 0
127  fields[j]->GetPlane(0)->BwdTrans_IterPerExp(
128  fields[j]->GetPlane(0)->GetCoeffs(),
129  fields[j]->GetPlane(0)->UpdatePhys());
130 
131  //bwd plane 1
132  fields[j]->GetPlane(1)->BwdTrans_IterPerExp(
133  fields[j]->GetPlane(1)->GetCoeffs(),
134  fields[j]->GetPlane(1)->UpdatePhys());
135  }
136  }
137  else
138  {
139  nq = fields[0]->GetTotPoints();
140  for (int j = 0; j < nfields; ++j)
141  {
142  for (int i = 0; i < fielddata.size(); ++i)
143  {
144  fields[j]->ExtractDataToCoeffs(fielddef[i],
145  fielddata[i],
146  fielddef[i]->m_fields[j],
147  fields[j]->UpdateCoeffs());
148  }
149  fields[j]->BwdTrans_IterPerExp(fields[j]->GetCoeffs(),
150  fields[j]->UpdatePhys());
151  }
152  }
153 
154 
155  // store mesh0 quadrature points
156  Array<OneD, NekDouble> x0(nq);
157  Array<OneD, NekDouble> y0(nq);
158  Array<OneD, NekDouble> z0(nq);
159 
160  if (fielddef[0]->m_numHomogeneousDir == 1)
161  {
162  fields[0]->GetPlane(1)->GetCoords(x0, y0, z0);
163  }
164  else
165  {
166  if (expdim == 2)
167  {
168  fields[0]->GetCoords(x0, y0);
169  }
170  else if (expdim == 3)
171  {
172  fields[0]->GetCoords(x0, y0, z0);
173  }
174  }
175 
176 
177  //----------------------------------------------
178 
179  //Read in the mesh1
180  string meshfile1(argv[argc-2]);
181  //Name of the output fieldfile
182  string fieldfile1(argv[argc-1]);
183  //-----------------------------------------------
184 
185  //define the output field:
186  Array<OneD, MultiRegions::ExpListSharedPtr> outfield;
187 
188  outfield = Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
189 
190  //set output fields over the graphShPt1 graph
191  SpatialDomains::MeshGraphSharedPtr graphShPt1;
192 
193  //remark: homo cases malloc() error or segmentation fault..
194  //remember: there is static cnt in stfields function to define
195  //the homo quantities only for the mesh0 case
196 
197  graphShPt1 = SpatialDomains::MeshGraph::Read(meshfile1);
198 
199  //define second session1..
200  std::vector<std::string> filenames1;
201  filenames1.push_back(meshfile1);
202  filenames1.push_back(meshfile1);
203 
204  argc = 2;
205  argv[1] = argv[3];
206  argv[2] = argv[4];
207 
208  LibUtilities::SessionReaderSharedPtr vSession1
209  = LibUtilities::SessionReader::CreateInstance(argc, argv,
210  filenames1,
211  vSession->GetComm());
212  //--------------------------------------------------
213 
214  //set explist 1
215  SetFields(graphShPt1, fielddef, vSession1, outfield, nfields,
216  variables, homo);
217  //------------------------------------------------
218  // store the new points:
219  int nq1;
220  if (fielddef[0]->m_numHomogeneousDir == 1)
221  {
222  nq1 = outfield[0]->GetPlane(1)->GetTotPoints();
223  }
224  else
225  {
226  nq1 = outfield[0]->GetTotPoints();
227  }
228  Array<OneD, NekDouble> x1(nq1);
229  Array<OneD, NekDouble> y1(nq1);
230  Array<OneD, NekDouble> z1(nq1);
231 
232  if (fielddef[0]->m_numHomogeneousDir == 1)
233  {
234  outfield[0]->GetPlane(1)->GetCoords(x1, y1, z1);
235  }
236  else
237  {
238  if (expdim == 2)
239  {
240  outfield[0]->GetCoords(x1, y1);
241  }
242  else if (expdim == 3)
243  {
244  outfield[0]->GetCoords(x1, y1, z1);
245  }
246  }
247 
248  //------------------------------------------------
249  //check 2Dmeshes compatibilities
250  // same max min values x,y...
251  //bool check = false;
252  bool check = true;
253  if (expdim == 2)
254  {
255  //check = Checkbndmeshes2D(x0, y0, x1, y1);
256  }
257  else if (expdim == 3)
258  {
259  check = Checkbndmeshes3D(x0, y0, z0, x1, y1, z1);
260  }
261  ASSERTL0(check, "meshes not compatible (different borders)");
262  //-----------------------------------------------
263 
264  if (fielddef[0]->m_numHomogeneousDir == 1)
265  {
266  for (int t = 0; t < nfields; t++)
267  {
268  InterpolateFieldHomo(fields[t], x1, y1, outfield[t]);
269  }
270  }
271  else
272  {
273  cout << "Interpolating [" << flush;
274 
275  if (expdim == 2)
276  {
277  InterpolateField(fields, outfield, x1, y1);
278  }
279  else if (expdim == 3)
280  {
281  NekDouble clamp_up, clamp_low;
282  vSession->LoadParameter("ClampToUpperValue", clamp_up, 10000000);
283  vSession->LoadParameter("ClampToLowerValue", clamp_low, -10000000);
284  InterpolateField(fields, outfield, x1, y1, z1, clamp_low, clamp_up);
285  }
286  cout << "]" << endl;
287  }
288 
289 
290  //------------------------------------------------
291  //write fieldfile
292  Writefield(vSession, variables, fieldfile1, graphShPt1, outfield);
293  //------------------------------------------------
294 
295 }
296 
297 
298 
299 // Define Expansion
300 void SetFields(
301  SpatialDomains::MeshGraphSharedPtr &mesh,
302  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef,
303  LibUtilities::SessionReaderSharedPtr &session,
304  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,
305  int nvariables,
306  const vector<std::string> &variables,
307  bool homogeneous)
308 {
309  // Session reader stuff has to be evaluated only from the
310  // first session which refers to mesh0
311  static int cnt = 0;
312 
313  // Setting parameteres for homogenous problems
314  NekDouble static LhomY;///< physical length in Y direction (if homogeneous)
315  NekDouble static LhomZ;///< physical length in Z direction (if homogeneous)
316 
317  bool static DeclareCoeffPhysArrays = true;
318 
319  int static npointsY;///< number of points in Y direction (if homogeneous)
320  int static npointsZ;///< number of points in Z direction (if homogeneous)
321 
322  bool static useFFT = false;
323  bool static deal = false;
324 
325  cnt++;
326  int i;
327  int expdim = mesh->GetMeshDimension();
328 
329  switch (expdim)
330  {
331  case 1:
332  {
333  if (fielddef[0]->m_numHomogeneousDir == 1)
334  {
335  const LibUtilities::PointsKey PkeyY(
336  npointsY,
337  LibUtilities::eFourierEvenlySpaced);
338  const LibUtilities::BasisKey BkeyY(
339  LibUtilities::eFourier,
340  npointsY, PkeyY);
341 
342  for (i = 0 ; i < nvariables; i++)
343  {
344  Exp[i] = MemoryManager<MultiRegions::
345  ContField3DHomogeneous1D>::AllocateSharedPtr(
346  session, BkeyY, LhomY,
347  useFFT, deal, mesh,
348  variables[i]);
349  }
350  }
351  else
352  {
353  for (i = 0 ; i < nvariables; i++)
354  {
355  Exp[i] = MemoryManager<MultiRegions::ContField1D>
356  ::AllocateSharedPtr(session, mesh, variables[i]);
357  }
358  }
359  break;
360  }
361  case 2:
362  {
363  if (fielddef[0]->m_numHomogeneousDir == 1)
364  {
365  const LibUtilities::PointsKey PkeyZ(
366  npointsZ,
367  LibUtilities::eFourierEvenlySpaced);
368  const LibUtilities::BasisKey BkeyZ(
369  LibUtilities::eFourier,
370  npointsZ, PkeyZ);
371 
372  MultiRegions::ContField3DHomogeneous1DSharedPtr firstfield =
373  MemoryManager<MultiRegions::ContField3DHomogeneous1D>::
374  AllocateSharedPtr(session, BkeyZ, LhomZ, useFFT,
375  deal, mesh, variables[0]);
376 
377  Exp[0] = firstfield;
378  for (i = 1 ; i < nvariables; i++)
379  {
380  Exp[i] = MemoryManager<MultiRegions::
381  ContField3DHomogeneous1D>::
382  AllocateSharedPtr(*firstfield);
383  }
384  }
385  else
386  {
387  i = 0;
388  MultiRegions::DisContField2DSharedPtr firstfield;
389  firstfield = MemoryManager<MultiRegions::DisContField2D>::
390  AllocateSharedPtr(session, mesh, variables[i],
391  DeclareCoeffPhysArrays);
392 
393  Exp[0] = firstfield;
394  for (i = 1 ; i < nvariables; i++)
395  {
396  Exp[i] = MemoryManager<MultiRegions::DisContField2D>
397  ::AllocateSharedPtr(*firstfield, mesh, variables[i],
398  DeclareCoeffPhysArrays);
399  }
400  }
401  break;
402  }
403  case 3:
404  {
405  if (fielddef[0]->m_numHomogeneousDir == 1)
406  {
407  ASSERTL0(false,
408  "3D fully periodic problems not implemented yet");
409  }
410  else
411  {
412  i = 0;
413  MultiRegions::ExpList3DSharedPtr firstfield =
414  MemoryManager<MultiRegions::ExpList3D>::
415  AllocateSharedPtr(session, mesh);
416 
417  Exp[0] = firstfield;
418  for (i = 1 ; i < nvariables; i++)
419  {
420  Exp[i] = MemoryManager<MultiRegions::ExpList3D>
421  ::AllocateSharedPtr(*firstfield);
422  }
423  }
424  break;
425  }
426  default:
427  ASSERTL0(false, "Expansion dimension not recognised");
428  break;
429  }
430 
431 }
432 
433 
434 
435 void InterpolateField(
436  Array<OneD, MultiRegions::ExpListSharedPtr> &field0,
437  Array<OneD, MultiRegions::ExpListSharedPtr> &field1,
438  Array<OneD, NekDouble> x,
439  Array<OneD, NekDouble> y,
440  Array<OneD, NekDouble> z,
441  NekDouble clamp_low,
442  NekDouble clamp_up)
443 {
444  int expdim = (z == NullNekDouble1DArray)? 2: 3;
445 
446  Array<OneD, NekDouble> coords(expdim), Lcoords(expdim);
447  int nq1 = field1[0]->GetTotPoints();
448  int elmtid, offset;
449  int r, f;
450  static int intpts = 0;
451 
452  ASSERTL0(field0.num_elements() == field1.num_elements(),
453  "Input field dimension must be same as output dimension");
454 
455  for (r = 0; r < nq1; r++)
456  {
457  coords[0] = x[r];
458  coords[1] = y[r];
459  if (expdim == 3)
460  {
461  coords[2] = z[r];
462  }
463 
464  // Obtain Element and LocalCoordinate to interpolate
465  elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3);
466 
467  offset = field0[0]->GetPhys_Offset(field0[0]->
468  GetOffset_Elmt_Id(elmtid));
469 
470  for (f = 0; f < field1.num_elements(); ++f)
471  {
472  NekDouble value;
473  value = field0[f]->GetExp(elmtid)->
474  StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
475 
476  if ((boost::math::isnan)(value))
477  {
478  ASSERTL0(false, "new value is not a number");
479  }
480  else
481  {
482  value = (value > clamp_up)? clamp_up :
483  ((value < clamp_low)? clamp_low :
484  value);
485 
486  field1[f]->UpdatePhys()[r] = value;
487  }
488  }
489 
490  if (intpts%1000 == 0)
491  {
492  cout <<"." << flush;
493  }
494  intpts ++;
495  }
496 }
497 
498 
499 
500 void InterpolateFieldHomo(
501  MultiRegions::ExpListSharedPtr field0,
502  Array<OneD, NekDouble> x1,
503  Array<OneD, NekDouble> y1,
504  MultiRegions::ExpListSharedPtr field1)
505 {
506  Array<OneD, NekDouble> coords(2);
507  int nq1 = field1->GetPlane(1)->GetTotPoints();
508  int elmtid, offset;
509 
510  //plane 0
511  for (int r = 0; r < nq1; r++)
512  {
513  coords[0] = x1[r];
514  coords[1] = y1[r];
515 
516  elmtid = field0->GetPlane(0)->GetExpIndex(coords, 1e-3);
517  offset = field0->GetPlane(0)->GetPhys_Offset(elmtid);
518  field1->GetPlane(0)->UpdatePhys()[r] = field0->GetPlane(0)->
519  GetExp(elmtid)->PhysEvaluate(
520  coords,
521  field0->GetPlane(0)->GetPhys() +offset);
522 
523  if ((boost::math::isnan)(field1->GetPlane(0)->UpdatePhys()[r]))
524  {
525  //ASSERTL0(abs(field1->UpdatePhys()[r])<10000000000,
526  // "interp failed");
527  }
528  }
529 
530 
531  // Plane1
532  for (int r = 0; r < nq1; r++)
533  {
534  coords[0] = x1[r];
535  coords[1] = y1[r];
536 
537  elmtid = field0->GetPlane(1)->GetExpIndex(coords, 1e-3);
538  offset = field0->GetPlane(1)->GetPhys_Offset(elmtid);
539  field1->GetPlane(1)->UpdatePhys()[r] = field0->GetPlane(1)->
540  GetExp(elmtid)->PhysEvaluate(
541  coords,
542  field0->GetPlane(1)->GetPhys() +offset);
543 
544  if((boost::math::isnan)(field1->GetPlane(1)->UpdatePhys()[r]))
545  {
546  //ASSERTL0(abs(field1->UpdatePhys()[r])<10000000000,
547  // "interp failed");
548  }
549 
550  }
551 }
552 
553 
554 
555 bool Checkbndmeshes2D(
556  Array<OneD, NekDouble> x0,
557  Array<OneD, NekDouble> y0,
558  Array<OneD, NekDouble> x1,
559  Array<OneD, NekDouble> y1)
560 {
561  NekDouble x0min, x0max, y0min, y0max;
562  NekDouble x1min, x1max, y1min, y1max;
563  NekDouble tol = 0.0000001;
564  NekDouble tol1 = 0.00001;
565  x0min = Vmath::Vmin(x0.num_elements(), x0, 1);
566  x0max = Vmath::Vmax(x0.num_elements(), x0, 1);
567  y0min = Vmath::Vmin(y0.num_elements(), y0, 1);
568  y0max = Vmath::Vmax(y0.num_elements(), y0, 1);
569 
570  x1min = Vmath::Vmin(x1.num_elements(), x1, 1);
571  x1max = Vmath::Vmax(x1.num_elements(), x1, 1);
572  y1min = Vmath::Vmin(y1.num_elements(), y1, 1);
573  y1max = Vmath::Vmax(y1.num_elements(), y1, 1);
574 
575  if (abs(x0min-x1min) < tol1 &&
576  abs(x0max-x1max) < tol1 &&
577  abs(y0min-y1min) < tol1 &&
578  abs(y0max-y1max) < tol1)
579  {
580  if (abs(x0min-x1min) > tol ||
581  abs(x0max-x1max) > tol ||
582  abs(y0min-y1min) > tol ||
583  abs(y0max-y1max) > tol)
584  {
585  cout<<"Warning: mesh boundary points differ more than 10^-7"<<endl;
586  }
587  return true;
588  }
589  else
590  {
591  return false;
592  }
593 }
594 
595 
596 
597 bool Checkbndmeshes3D(
598  Array<OneD, NekDouble> x0,
599  Array<OneD, NekDouble> y0,
600  Array<OneD, NekDouble> z0,
601  Array<OneD, NekDouble> x1,
602  Array<OneD, NekDouble> y1,
603  Array<OneD, NekDouble> z1)
604 {
605  NekDouble x0min, x0max, y0min, y0max, z0min, z0max;
606  NekDouble x1min, x1max, y1min, y1max, z1min, z1max;
607  NekDouble tol = 0.0000001;
608  NekDouble tol1 = 0.00001;
609  x0min = Vmath::Vmin(x0.num_elements(), x0, 1);
610  x0max = Vmath::Vmax(x0.num_elements(), x0, 1);
611  y0min = Vmath::Vmin(y0.num_elements(), y0, 1);
612  y0max = Vmath::Vmax(y0.num_elements(), y0, 1);
613  z0min = Vmath::Vmin(z0.num_elements(), z0, 1);
614  z0max = Vmath::Vmax(z0.num_elements(), z0, 1);
615 
616  x1min = Vmath::Vmin(x1.num_elements(), x1, 1);
617  x1max = Vmath::Vmax(x1.num_elements(), x1, 1);
618  y1min = Vmath::Vmin(y1.num_elements(), y1, 1);
619  y1max = Vmath::Vmax(y1.num_elements(), y1, 1);
620  z1min = Vmath::Vmin(z1.num_elements(), z1, 1);
621  z1max = Vmath::Vmax(z1.num_elements(), z1, 1);
622 
623  if (abs(x0min-x1min) < tol1 &&
624  abs(x0max-x1max) < tol1 &&
625  abs(y0min-y1min) < tol1 &&
626  abs(y0max-y1max) < tol1 &&
627  abs(z0min-z1min) < tol1 &&
628  abs(z0max-z1max) < tol1)
629  {
630  if (abs(x0min-x1min) > tol ||
631  abs(x0max-x1max) > tol ||
632  abs(y0min-y1min) > tol ||
633  abs(y0max-y1max) > tol ||
634  abs(z0min-z1min) > tol ||
635  abs(z0max-z1max) > tol)
636  {
637  cout<<"Warning: mesh boundary points differ more than 10^-7"<<endl;
638  }
639  return true;
640  }
641  else
642  {
643  return false;
644  }
645 }
646 
647 
648 
649 void Writefield(
650  LibUtilities::SessionReaderSharedPtr vSession,
651  const std::vector<std::string> &variables,
652  string fieldfile,
653  SpatialDomains::MeshGraphSharedPtr &graph,
654  Array<OneD, MultiRegions::ExpListSharedPtr> &outfield)
655 {
656  string var;
657  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
658  = outfield[0]->GetFieldDefinitions();
659  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
660  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(outfield.num_elements());
661 
662  for (int j = 0; j < fieldcoeffs.num_elements(); ++j)
663  {
664  if (FieldDef[0]->m_numHomogeneousDir == 1)
665  {
666  // plane 0
667  outfield[j]->GetPlane(0)->FwdTrans_IterPerExp(
668  outfield[j]->GetPlane(0)->GetPhys(),
669  outfield[j]->GetPlane(0)->UpdateCoeffs());
670 
671  // plane 1
672  outfield[j]->GetPlane(1)->FwdTrans_IterPerExp(
673  outfield[j]->GetPlane(1)->GetPhys(),
674  outfield[j]->GetPlane(1)->UpdateCoeffs());
675  }
676  else
677  {
678  outfield[j]->FwdTrans_IterPerExp(
679  outfield[j]->GetPhys(),
680  outfield[j]->UpdateCoeffs());
681  }
682 
683  fieldcoeffs[j] = outfield[j]->UpdateCoeffs();
684 
685  for (int i = 0; i < FieldDef.size(); i++)
686  {
687  FieldDef[i]->m_fields.push_back(variables[j]);
688  outfield[0]->AppendFieldData(FieldDef[i], FieldData[i],
689  fieldcoeffs[j]);
690  }
691  }
692  LibUtilities::Write(fieldfile, FieldDef, FieldData);
693 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2