Nektar++
ProcessScalGrad.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessScalGrad.cpp
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: Computes scalar gradient field.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessScalGrad.h"
41 
44 #include <MultiRegions/ExpList.h>
45 
46 namespace Nektar
47 {
48 namespace Utilities
49 {
50 
51 ModuleKey ProcessScalGrad::className =
53  ModuleKey(eProcessModule, "scalargrad"),
54  ProcessScalGrad::create, "Computes scalar gradient field.");
55 
56 ProcessScalGrad::ProcessScalGrad(FieldSharedPtr f) : ProcessModule(f)
57 {
58  m_config["bnd"] = ConfigOption(false,"All","Boundary to be extracted");
59  f->m_writeBndFld = true;
60  f->m_declareExpansionAsContField = true;
61  m_f->m_fldToBnd = false;
62 }
63 
65 {
66 }
67 
68 void ProcessScalGrad::Process(po::variables_map &vm)
69 {
70  int i, j, k;
71  if (m_f->m_verbose)
72  {
73  cout << "ProcessScalGrad: Calculating scalar gradient..." << endl;
74  }
75 
76  // Set up Field options to output boundary fld
77  string bvalues = m_config["bnd"].as<string>();
78 
79  if(bvalues.compare("All") == 0)
80  {
82  BndExp = m_f->m_exp[0]->GetBndCondExpansions();
83 
84  for(i = 0; i < BndExp.num_elements(); ++i)
85  {
86  m_f->m_bndRegionsToWrite.push_back(i);
87  }
88  }
89  else
90  {
92  m_f->m_bndRegionsToWrite),"Failed to interpret range string");
93  }
94 
95  int spacedim = m_f->m_graph->GetSpaceDimension();
96  if ((m_f->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
97  (m_f->m_fielddef[0]->m_numHomogeneousDir) == 2)
98  {
99  spacedim = 3;
100  }
101 
102  int nfields = m_f->m_fielddef[0]->m_fields.size();
103  //ASSERTL0(nfields == 1,"Implicit assumption that input is in ADR format of (u)");
104 
105  if (spacedim == 1)
106  {
107  ASSERTL0(false, "Error: scalar gradient for a 1D problem cannot "
108  "be computed");
109  }
110 
111 
112  int ngrad = spacedim;
113  int n, cnt, elmtid, nq, offset, boundary, nfq;
114  int npoints = m_f->m_exp[0]->GetNpoints();
115  string var;
116  Array<OneD, NekDouble> scalar;
117  Array<OneD, Array<OneD, NekDouble> > grad(ngrad), fgrad(ngrad), outfield(nfields);
118 
121  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
123 
124  m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID, BoundarytoTraceID);
125 
126  for (i = 0; i < nfields; i++)
127  {
128  var = m_f->m_fielddef[0]->m_fields[i];
129  stringstream filename;
130  filename << var << "_scalar_gradient";
131  filename >> var;
132  m_f->m_fielddef[0]->m_fields[i] = var;
133 
134  BndExp[i] = m_f->m_exp[i]->GetBndCondExpansions();
135  outfield[i] = Array<OneD, NekDouble>(npoints);
136  }
137 
138  // loop over the types of boundary conditions
139  for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
140  {
141  bool doneBnd = false;
142  // identify if boundary has been defined
143  for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
144  {
145  if(n == m_f->m_bndRegionsToWrite[b])
146  {
147  doneBnd = true;
148  for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
149  {
150  // find element and face of this expansion.
151  elmtid = BoundarytoElmtID[cnt];
152  elmt = m_f->m_exp[0]->GetExp(elmtid);
153  nq = elmt->GetTotPoints();
154  offset = m_f->m_exp[0]->GetPhys_Offset(elmtid);
155 
156  // Initialise local arrays for the velocity gradients, and stress components
157  // size of total number of quadrature points for each element (hence local).
158  for(j = 0; j < ngrad; ++j)
159  {
160  grad[j] = Array<OneD, NekDouble>(nq);
161  }
162 
163  if(spacedim == 2)
164  {
165  //Not implemented in 2D.
166  }
167  else
168  {
169  for(j = 0; j < nfields; j++)
170  {
171  outfield[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][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  nfq = bc->GetTotPoints();
177 
178  //identify boundary of element looking at.
179  boundary = BoundarytoTraceID[cnt];
180 
181  //Get face normals
182  const SpatialDomains::GeomFactorsSharedPtr m_metricinfo = bc->GetMetricInfo();
183 
185  = elmt->GetFaceNormal(boundary);
186 
187  // initialise arrays
188  for(j = 0; j < ngrad; ++j)
189  {
190  fgrad[j] = Array<OneD, NekDouble>(nfq);
191  }
192  Array<OneD, NekDouble> gradnorm(nfq);
193 
194  for(k = 0; k < nfields; k++)
195  {
196  Vmath::Zero(nfq, gradnorm, 1);
197 
198  scalar = m_f->m_exp[k]->GetPhys() + offset;
199  elmt->PhysDeriv(scalar, grad[0],grad[1],grad[2]);
200 
201  for(j = 0; j < ngrad; ++j)
202  {
203  elmt->GetFacePhysVals(boundary,bc,grad[j],fgrad[j]);
204  }
205 
206  //surface curved
207  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
208  {
209  for (j=0; j<ngrad; j++)
210  {
211  Vmath::Vvtvp(nfq, normals[j], 1, fgrad[j], 1, gradnorm, 1, gradnorm, 1);
212  }
213  }
214  else
215  {
216  for (j=0; j<ngrad; j++)
217  {
218  Vmath::Svtvp(nfq, normals[j][0], fgrad[j], 1, gradnorm, 1, gradnorm, 1);
219  }
220  }
221  bc->FwdTrans(gradnorm, outfield[k]);
222  }
223 
224  }
225  }
226  }
227  }
228  if(doneBnd == false)
229  {
230  cnt += BndExp[0][n]->GetExpSize();
231  }
232  }
233 
234  for(j = 0; j < nfields; ++j)
235  {
236  for(int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
237  {
238  m_f->m_exp[j]->UpdateBndCondExpansion(m_f->m_bndRegionsToWrite[b]) = BndExp[j][m_f->m_bndRegionsToWrite[b]];
239  }
240 
241  }
242 }
243 
244 }
245 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:96
virtual void Process()=0
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
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
FieldSharedPtr m_f
Field object.
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
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:677
Represents a command-line configuration option.
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Geometry is curved or has non-constant factors.
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215