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