Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessCyl.cpp
Go to the documentation of this file.
1 #include <string>
2 using namespace std;
3 
4 #include "MeshElements.h"
5 #include "ProcessCyl.h"
6 
7 #include <LocalRegions/SegExp.h>
8 #include <LocalRegions/QuadExp.h>
9 #include <LocalRegions/TriExp.h>
12 
13 namespace Nektar
14 {
15  namespace Utilities
16  {
17  ModuleKey ProcessCyl::className =
19  ModuleKey(eProcessModule, "cyl"),
20  ProcessCyl::create);
21 
22  /**
23  * @brief Default constructor.
24  */
25  ProcessCyl::ProcessCyl(MeshSharedPtr m) : ProcessModule(m)
26  {
27  m_config["surf"] = ConfigOption(false, "-1",
28  "Tag identifying surface to process.");
29  m_config["r"] = ConfigOption(false, "0.0",
30  "Radius of cylinder.");
31  m_config["N"] = ConfigOption(false, "7",
32  "Number of points along edge.");
33  }
34 
35 
36  /**
37  * @brief Destructor.
38  */
40  {
41  }
42 
43 
45  {
46  set<pair<int,int> > tmp;
47  int surfTag = m_config["surf"].as<int>();
48  int prismedge[2][3] = {{0,5,4}, {2,6,7}};
49  Node zax(0,0,0,1);
50  int dim = m_mesh->m_expDim;
51 
52  for (int i = 0; i < m_mesh->m_element[dim].size(); ++i)
53  {
54  ElementSharedPtr el = m_mesh->m_element[dim][i];
55  int nSurf = dim == 3 ? el->GetFaceCount() : el->GetEdgeCount();
56 
57  for (int j = 0; j < nSurf; ++j)
58  {
59  int bl = el->GetBoundaryLink(j);
60  if (bl == -1)
61  {
62  continue;
63  }
64 
65  ElementSharedPtr bEl = m_mesh->m_element[dim - 1][bl];
66  vector<int> tags = bEl->GetTagList();
67 
68  if (find(tags.begin(), tags.end(), surfTag) ==
69  tags.end())
70  {
71  continue;
72  }
73 
74  ASSERTL0(j == 1 || j == 3, "rofl");
75 
76  // Check all edge interior points.
77  for (int k = 0; k < 3; ++k)
78  {
79  EdgeSharedPtr edge = el->GetEdge(prismedge[(j-1)/2][k]);
80  GenerateEdgeNodes(edge);
81  }
82  }
83  }
84  }
85 
86 
88  {
89  NodeSharedPtr n1 = edge->m_n1;
90  NodeSharedPtr n2 = edge->m_n2;
91 
92  int nq = m_config["N"].as<int>();
93  double r = m_config["r"].as<double>();
94  double t1 = atan2(n1->m_y, n1->m_x);
95  double t2 = atan2(n2->m_y, n2->m_x);
96  double dt;
97  double dz;
98 
99  if (t1 < -M_PI/2.0 && t2 > 0.0)
100  {
101  t1 += 2*M_PI;
102  }
103  if (t2 < -M_PI/2.0 && t1 > 0.0)
104  {
105  t2 += 2*M_PI;
106  }
107 
108  dt = (t2-t1) / (nq-1);
109  dz = (n2->m_z - n1->m_z) / (nq-1);
110 
111  edge->m_edgeNodes.resize(nq-2);
112  Node dx = (*n2-*n1) * (1.0/(nq-1.0));
113  for (int i = 1; i < nq-1; ++i)
114  {
115  edge->m_edgeNodes[i-1] = NodeSharedPtr(
116  new Node(0, r*cos(t1 + i*dt),
117  r*sin(t1 + i*dt),
118  n1->m_z + i*dz));
119  }
120  edge->m_curveType = LibUtilities::ePolyEvenlySpaced;
121  }
122  }
123 }