Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 <string>
37 #include <iostream>
38 using namespace std;
39 
40 #include "ProcessSurfDistance.h"
41 
42 namespace Nektar
43 {
44 namespace Utilities
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)
54  : ProcessModule(f)
55 {
56  m_config["bnd"] = 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  int i, j, k, cnt;
69  int surf = m_config["bnd"].as<int>();
70  int expdim = m_f->m_graph->GetMeshDimension();
71 
72  ASSERTL0(surf >= 0, "Invalid surface "+boost::lexical_cast<string>(surf));
73 
74  // Add this boundary region to the list that we will output.
75  m_f->m_bndRegionsToWrite.push_back(surf);
76 
77  // Remove existing fields.
78  m_f->m_exp.resize(1);
79 
80  // Grab boundary expansions.
82  m_f->m_exp[0]->GetBndCondExpansions();
83 
84  // Get map that takes us from boundary element to element.
85  Array<OneD, int> BoundarytoElmtID, BoundarytoTraceID;
86  m_f->m_exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID, BoundarytoTraceID);
87 
88  if (m_f->m_fielddef.size() == 0)
89  {
90  m_f->m_fielddef = m_f->m_exp[0]->GetFieldDefinitions();
91  m_f->m_fielddef[0]->m_fields.push_back("dist");
92  }
93  else
94  {
95  // Override field variable
96  m_f->m_fielddef[0]->m_fields[0] = "dist";
97  }
98 
99  ASSERTL0(!(m_f->m_fielddef[0]->m_numHomogeneousDir),
100  "Homogeneous expansions not supported");
101 
102  for (i = cnt = 0; i < BndExp.num_elements(); ++i)
103  {
104  if (i != surf)
105  {
106  cnt += BndExp[i]->GetExpSize();
107  continue;
108  }
109 
110  for (j = 0; j < BndExp[i]->GetExpSize(); ++j, ++cnt)
111  {
112  int elmtNum = BoundarytoElmtID [cnt];
113  int facetNum = BoundarytoTraceID[cnt];
114  int oppositeNum;
115 
116  // Get boundary and element expansions.
117  LocalRegions::ExpansionSharedPtr bndElmt = BndExp[i]->GetExp(j);
119  m_f->m_exp[0]->GetExp(elmtNum);
120 
121  // Determine which face is opposite to the surface
122  switch(elmt->DetShapeType())
123  {
125  {
126  oppositeNum = (facetNum+2)%4;
127  }
128  break;
129 
131  {
132  switch(facetNum)
133  {
134  case 1:
135  oppositeNum = 3;
136  break;
137  case 3:
138  oppositeNum = 1;
139  break;
140  default:
141  ASSERTL0(false,
142  "Surface must be on a triangular face of the prism.");
143  }
144  }
145  break;
146 
148  {
149  switch(facetNum)
150  {
151  case 0:
152  oppositeNum = 5;
153  break;
154  case 1:
155  oppositeNum = 3;
156  break;
157  case 2:
158  oppositeNum = 4;
159  break;
160  case 3:
161  oppositeNum = 1;
162  break;
163  case 4:
164  oppositeNum = 2;
165  break;
166  case 5:
167  oppositeNum = 0;
168  break;
169  default:
170  ASSERTL0(false, "Face out of bound.");
171  }
172  }
173  break;
174 
175  default:
176  ASSERTL0(false, "Element not supported");
177  }
178 
179  int nq = elmt ->GetTotPoints();
180  int nqBnd = bndElmt->GetTotPoints();
181 
183  x[0] = Array<OneD, NekDouble>(nq);
184  x[1] = Array<OneD, NekDouble>(nq);
185  x[2] = Array<OneD, NekDouble>(nq);
186  elmt->GetCoords(x[0], x[1], x[2]);
187 
188  Array<OneD, NekDouble> face1(nqBnd), face2(nqBnd);
190  BndExp[i]->UpdatePhys() + BndExp[i]->GetPhys_Offset(j);
191 
192  // Zero existing value.
193  Vmath::Zero(nqBnd, dist, 1);
194 
195  // Calculate distance between two faces of the element
196  for (k = 0; k < expdim; ++k)
197  {
198  switch(expdim)
199  {
200  case 2:
201  {
202  elmt->GetEdgePhysVals(facetNum, bndElmt, x[k], face1);
203  elmt->GetEdgePhysVals(oppositeNum, bndElmt, x[k], face2);
204  // Consider edge orientation
205  if (elmt->GetEorient(facetNum) ==
206  elmt->GetEorient(oppositeNum))
207  {
208  Vmath::Reverse(nqBnd, face2, 1, face2, 1);
209  }
210  }
211  break;
212  case 3:
213  {
214  // Use orientation from the surface for both faces
215  StdRegions::Orientation orientation =
216  elmt->GetForient(facetNum);
217  elmt->GetFacePhysVals(facetNum, bndElmt,
218  x[k], face1, orientation);
219  elmt->GetFacePhysVals(oppositeNum, bndElmt,
220  x[k], face2, orientation);
221  }
222  break;
223  default:
224  ASSERTL0(false, "Expansion not supported");
225  }
226  Vmath::Vsub (nqBnd, face1, 1, face2, 1, face1, 1);
227  Vmath::Vvtvp(nqBnd, face1, 1, face1, 1, dist, 1, dist, 1);
228  }
229  Vmath::Vsqrt(nqBnd, dist, 1, dist, 1);
230  }
231 
232  BndExp[i]->FwdTrans(BndExp[i]->GetPhys(), BndExp[i]->UpdateCoeffs());
233  }
234 }
235 
236 }
237 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
pair< ModuleType, string > ModuleKey
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
virtual void Process()=0
map< string, ConfigOption > m_config
List of configuration values.
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
STL namespace.
FieldSharedPtr m_f
Field object.
void Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1071
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:695
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:329
Represents a command-line configuration option.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
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