Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Functions
FldAddScalGrad.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField2D.h>
Include dependency graph for FldAddScalGrad.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 14 of file FldAddScalGrad.cpp.

References ASSERTL0, Nektar::SpatialDomains::eDeformed, Nektar::StdRegions::StdExpansion::GetTotPoints(), Nektar::LibUtilities::Import(), sign, Vmath::Smul(), Vmath::Svtsvtp(), Vmath::Svtvp(), Vmath::Vvtvp(), Vmath::Vvtvvtp(), and Nektar::LibUtilities::Write().

15 {
16  int i,j;
17  int surfID;
18 
19  if(argc != 5)
20  {
21  fprintf(stderr,"Usage: FldAddScalGrad meshfile infld outfld BoundaryID\n");
22  exit(1);
23  }
24 
25  surfID = boost::lexical_cast<int>(argv[argc - 1]);
26 
27  argv[argc -1] = argv[argc - 2];
28 
30  = LibUtilities::SessionReader::CreateInstance(argc, argv);
31 
32  //----------------------------------------------
33  // Read in mesh from input file
34  string meshfile(argv[argc-4]);
35  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
36  //----------------------------------------------
37 
38  //----------------------------------------------
39  // Import field file.
40  string fieldfile(argv[argc-3]);
41  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
42  vector<vector<NekDouble> > fielddata;
43  LibUtilities::Import(fieldfile,fielddef,fielddata);
44  //----------------------------------------------
45 
46  //----------------------------------------------
47  // Define Expansion
48  int expdim = graphShPt->GetMeshDimension();
49  int nfields = 1;
50  int addfields = 7;
51  Array<OneD, MultiRegions::ExpListSharedPtr> exp(nfields + addfields);
52  MultiRegions::AssemblyMapCGSharedPtr m_locToGlobalMap;
53 
54  switch(expdim)
55  {
56  case 1:
57  {
58  ASSERTL0(false,"Expansion dimension not recognised");
59  }
60  break;
61  case 2:
62  {
63  ASSERTL0(false,"Expansion dimension not recognised");
64  }
65  break;
66  case 3:
67  {
70  ::AllocateSharedPtr(vSession, graphShPt,
71  vSession->GetVariable(0));
72 
73  m_locToGlobalMap = originalfield->GetLocalToGlobalMap();
74 
75  exp[0] = originalfield;
76  for (i=0; i<addfields; i++)
77  {
79  ::AllocateSharedPtr(*originalfield, graphShPt,
80  vSession->GetVariable(0));
81  }
82  }
83  break;
84  default:
85  ASSERTL0(false,"Expansion dimension not recognised");
86  break;
87  }
88  //----------------------------------------------
89 
90  //----------------------------------------------
91  // Copy data from field file
92  for(j = 0; j < nfields+addfields; ++j)
93  {
94  for(int i = 0; i < fielddata.size(); ++i)
95  {
96  exp[j]->ExtractDataToCoeffs(fielddef [i],
97  fielddata[i],
98  fielddef [i]->m_fields[0],
99  exp[j]->UpdateCoeffs());
100  }
101 
102  exp[j]->BwdTrans(exp[j]->GetCoeffs(),exp[j]->UpdatePhys());
103  }
104 
105 
106  //----------------------------------------------
107 
108  //----------------------------------------------
109  int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
110  nt = exp[0]->GetNpoints();
111  Array<OneD, Array<OneD, NekDouble> > grad(expdim);
112  Array<OneD, Array<OneD, NekDouble> > fgrad(expdim);
113  Array<OneD, Array<OneD, NekDouble> > values(addfields);
114 
115  // Set up mapping from Boundary condition to element details.
118  Array<OneD, int> BoundarytoElmtID;
119  Array<OneD, int> BoundarytoTraceID;
122  Array<OneD, NekDouble> outvalues;
123 
124  exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
125 
126  //get boundary expansions for each field
127  for (i = 0; i<addfields; i++)
128  {
129  BndExp[i] = exp[i]->GetBndCondExpansions();
130  }
131 
132  // loop over the types of boundary conditions
133  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
134  {
135  // identify boundary which the user wanted
136  if(n == surfID)
137  {
138  for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
139  {
140  // find element and face of this expansion.
141  elmtid = BoundarytoElmtID[cnt];
142  elmt = exp[0]->GetExp(elmtid);
143  nq = elmt->GetTotPoints();
144  offset = exp[0]->GetPhys_Offset(elmtid);
145 
146  // Initialise local arrays for the velocity gradients
147  // size of total number of quadrature points for each element (hence local).
148  for(j = 0; j < expdim; ++j)
149  {
150  grad[j] = Array<OneD, NekDouble>(nq);
151  }
152 
153  if(expdim == 2)
154  {
155  }
156  else
157  {
158  for (j = 0; j< addfields; j++)
159  {
160  values[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][n]->GetCoeff_Offset(i);
161  }
162 
163  // Get face 2D expansion from element expansion
164  bc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (BndExp[0][n]->GetExp(i));
165 
166  // Number of face quadrature points
167  nfq = bc->GetTotPoints();
168 
169  //identify boundary of element
170  boundary = BoundarytoTraceID[cnt];
171 
172  //Extract scalar field
173  U = exp[0]->GetPhys() + offset;
174 
175  //Compute gradients
176  elmt->PhysDeriv(U,grad[0],grad[1],grad[2]);
177 
178  if(i ==0)
179  {
180  for (j = 0; j< nq; j++)
181  {
182  cout << "element grad: " << grad[0][j] << endl;
183  }
184  }
185 
186  for(j = 0; j < expdim; ++j)
187  {
188  fgrad[j] = Array<OneD, NekDouble>(nfq);
189  }
190 
191 
192  // Get gradient at the quadrature points of the face
193  for(j = 0; j < expdim; ++j)
194  {
195  elmt->GetFacePhysVals(boundary,bc,grad[j],fgrad[j]);
196  bc->FwdTrans(fgrad[j],values[j]);
197  }
198 
199  if(i ==0)
200  {
201  for (j = 0; j< nfq; j++)
202  {
203  cout << "face grad: " << fgrad[0][j] << endl;
204  }
205  }
206 
207  const SpatialDomains::GeomFactorsSharedPtr m_metricinfo=bc->GetMetricInfo();
208 
210  = elmt->GetFaceNormal(boundary);
211 
212  Array<OneD, NekDouble> gradnorm(nfq);
213 
214  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
215  {
216  Vmath::Vvtvvtp(nfq,normals[0],1,fgrad[0],1,
217  normals[1],1,fgrad[1],1,gradnorm,1);
218  Vmath::Vvtvp (nfq,normals[2],1,fgrad[2],1,gradnorm,1,gradnorm,1);
219  }
220  else
221  {
222  Vmath::Svtsvtp(nfq,normals[0][0],fgrad[0],1,
223  normals[1][0],fgrad[1],1,gradnorm,1);
224  Vmath::Svtvp(nfq,normals[2][0],fgrad[2],1,gradnorm,1,gradnorm,1);
225  }
226 
227  for(j = 0; j<expdim; j++)
228  {
229  bc->FwdTrans(normals[j],values[j+expdim]);
230  }
231 
232  //gradient (grad(u) n)
233  Vmath::Smul(nfq,-1.0,gradnorm,1,gradnorm,1);
234  bc->FwdTrans(gradnorm,values[expdim*2]);
235  }
236  }
237 
238  }
239  else
240  {
241  cnt += BndExp[0][n]->GetExpSize();
242  }
243  }
244 
245  for(j = 0; j < addfields; ++j)
246  {
247  int ncoeffs = exp[0]->GetNcoeffs();
248  Array<OneD, NekDouble> output(ncoeffs);
249 
250  output=exp[j+1]->UpdateCoeffs();
251 
252  int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
253  Array<OneD, NekDouble> outarray(nGlobal,0.0);
254 
255  int bndcnt=0;
256 
257  const Array<OneD,const int>& map = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsMap();
258  NekDouble sign;
259 
260  for(int i = 0; i < BndExp[j].num_elements(); ++i)
261  {
262  if(i==surfID)
263  {
264  const Array<OneD,const NekDouble>& coeffs = BndExp[j][i]->GetCoeffs();
265  for(int k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
266  {
267  sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
268  outarray[map[bndcnt++]] = sign * coeffs[k];
269  }
270  }
271  else
272  {
273  bndcnt += BndExp[j][i]->GetNcoeffs();
274  }
275  }
276  m_locToGlobalMap->GlobalToLocal(outarray,output);
277  }
278 
279 
280  //-----------------------------------------------
281  // Write solution to file with additional computed fields
282  string out(argv[argc-2]);
283  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
284  = exp[0]->GetFieldDefinitions();
285  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
286 
287  vector<string > outname;
288 
289  outname.push_back("du/dx");
290  outname.push_back("du/dy");
291  outname.push_back("du/dz");
292  outname.push_back("nx");
293  outname.push_back("ny");
294  outname.push_back("nz");
295  outname.push_back("gradient");
296 
297 
298  for(j = 0; j < nfields+addfields; ++j)
299  {
300  for(i = 0; i < FieldDef.size(); ++i)
301  {
302  if (j >= nfields)
303  {
304  FieldDef[i]->m_fields.push_back(outname[j-nfields]);
305  }
306  else
307  {
308  FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
309  }
310  exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
311  }
312  }
313 
314  LibUtilities::Write(out, FieldDef, FieldData);
315  //-----------------------------------------------
316 
317  return 0;
318 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:279
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:485
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
double NekDouble
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
Definition: FieldIO.cpp:235
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:537
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:591
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:208
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:52
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Geometry is curved or has non-constant factors.