Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessSurfDistance.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessSurfDistance.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 height of elements connected to a surface.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <string>
38 using namespace std;
39 
40 #include "ProcessSurfDistance.h"
41 
42 namespace Nektar
43 {
44 namespace FieldUtils
45 {
46 
47 ModuleKey ProcessSurfDistance::className =
49  ModuleKey(eProcessModule, "surfdistance"),
50  ProcessSurfDistance::create,
51  "Computes height of element connected to a surface.");
52 
53 ProcessSurfDistance::ProcessSurfDistance(FieldSharedPtr f) : ProcessModule(f)
54 {
55  m_config["bnd"] =
56  ConfigOption(false, "-1", "Boundary region to calculate height");
57  f->m_writeBndFld = true;
58  f->m_declareExpansionAsContField = true;
59  m_f->m_fldToBnd = false;
60 }
61 
63 {
64 }
65 
66 void ProcessSurfDistance::Process(po::variables_map &vm)
67 {
68  if (m_f->m_verbose)
69  {
70  if (m_f->m_comm->TreatAsRankZero())
71  {
72  cout << "ProcessSurfDistance: Calculating distance to surface..."
73  << endl;
74  }
75  }
76 
77  int i, j, k, cnt;
78  int surf = m_config["bnd"].as<int>();
79  int expdim = m_f->m_graph->GetMeshDimension();
80 
81  ASSERTL0(surf >= 0, "Invalid surface " + boost::lexical_cast<string>(surf));
82 
83  // Add this boundary region to the list that we will output.
84  m_f->m_bndRegionsToWrite.push_back(surf);
85 
86  // Remove existing fields.
87  m_f->m_exp.resize(1);
88 
89  // Grab boundary expansions.
91  m_f->m_exp[0]->GetBndCondExpansions();
92 
93  // Get map that takes us from boundary element to element.
94  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
95  m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID, BoundarytoTraceID);
96 
97  if (m_f->m_fielddef.size() == 0)
98  {
99  m_f->m_fielddef = m_f->m_exp[0]->GetFieldDefinitions();
100  m_f->m_fielddef[0]->m_fields.push_back("dist");
101  }
102  else
103  {
104  // Override field variable
105  m_f->m_fielddef[0]->m_fields[0] = "dist";
106  }
107 
108  ASSERTL0(!(m_f->m_fielddef[0]->m_numHomogeneousDir),
109  "Homogeneous expansions not supported");
110 
111  for (i = cnt = 0; i < BndExp.num_elements(); ++i)
112  {
113  if (i != surf)
114  {
115  cnt += BndExp[i]->GetExpSize();
116  continue;
117  }
118 
119  for (j = 0; j < BndExp[i]->GetExpSize(); ++j, ++cnt)
120  {
121  int elmtNum = BoundarytoElmtID[cnt];
122  int facetNum = BoundarytoTraceID[cnt];
123  int oppositeNum;
124 
125  // Get boundary and element expansions.
126  LocalRegions::ExpansionSharedPtr bndElmt = BndExp[i]->GetExp(j);
128  m_f->m_exp[0]->GetExp(elmtNum);
129 
130  // Determine which face is opposite to the surface
131  switch (elmt->DetShapeType())
132  {
134  {
135  oppositeNum = (facetNum + 2) % 4;
136  }
137  break;
138 
140  {
141  switch (facetNum)
142  {
143  case 1:
144  oppositeNum = 3;
145  break;
146  case 3:
147  oppositeNum = 1;
148  break;
149  default:
150  ASSERTL0(false, "Surface must be on a triangular "
151  "face of the prism.");
152  }
153  }
154  break;
155 
157  {
158  switch (facetNum)
159  {
160  case 0:
161  oppositeNum = 5;
162  break;
163  case 1:
164  oppositeNum = 3;
165  break;
166  case 2:
167  oppositeNum = 4;
168  break;
169  case 3:
170  oppositeNum = 1;
171  break;
172  case 4:
173  oppositeNum = 2;
174  break;
175  case 5:
176  oppositeNum = 0;
177  break;
178  default:
179  ASSERTL0(false, "Face out of bound.");
180  }
181  }
182  break;
183 
184  default:
185  ASSERTL0(false, "Element not supported");
186  }
187 
188  int nq = elmt->GetTotPoints();
189  int nqBnd = bndElmt->GetTotPoints();
190 
192  x[0] = Array<OneD, NekDouble>(nq);
193  x[1] = Array<OneD, NekDouble>(nq);
194  x[2] = Array<OneD, NekDouble>(nq);
195  elmt->GetCoords(x[0], x[1], x[2]);
196 
197  Array<OneD, NekDouble> face1(nqBnd), face2(nqBnd);
199  BndExp[i]->UpdatePhys() + BndExp[i]->GetPhys_Offset(j);
200 
201  // Zero existing value.
202  Vmath::Zero(nqBnd, dist, 1);
203 
204  // Calculate distance between two faces of the element
205  for (k = 0; k < expdim; ++k)
206  {
207  switch (expdim)
208  {
209  case 2:
210  {
211  elmt->GetEdgePhysVals(facetNum, bndElmt, x[k], face1);
212  elmt->GetEdgePhysVals(oppositeNum, bndElmt, x[k],
213  face2);
214  // Consider edge orientation
215  if (elmt->GetEorient(facetNum) ==
216  elmt->GetEorient(oppositeNum))
217  {
218  Vmath::Reverse(nqBnd, face2, 1, face2, 1);
219  }
220  }
221  break;
222  case 3:
223  {
224  // Use orientation from the surface for both faces
225  StdRegions::Orientation orientation =
226  elmt->GetForient(facetNum);
227  elmt->GetFacePhysVals(facetNum, bndElmt, x[k], face1,
228  orientation);
229  elmt->GetFacePhysVals(oppositeNum, bndElmt, x[k], face2,
230  orientation);
231  }
232  break;
233  default:
234  ASSERTL0(false, "Expansion not supported");
235  }
236  Vmath::Vsub(nqBnd, face1, 1, face2, 1, face1, 1);
237  Vmath::Vvtvp(nqBnd, face1, 1, face1, 1, dist, 1, dist, 1);
238  }
239  Vmath::Vsqrt(nqBnd, dist, 1, dist, 1);
240  }
241 
242  BndExp[i]->FwdTrans(BndExp[i]->GetPhys(), BndExp[i]->UpdateCoeffs());
243  }
244 }
245 }
246 }
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual void Process(po::variables_map &vm)
Write mesh to output file.
Represents a command-line configuration option.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:408
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
STL namespace.
pair< ModuleType, string > ModuleKey
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1085
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:767
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:343
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Abstract base class for processing modules.
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215