Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessPerAlign.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessPerAlign.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: Reorder composites to align periodic boundaries.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <string>
37 using namespace std;
38 
39 #include "MeshElements.h"
40 #include "ProcessPerAlign.h"
41 
42 #include <LocalRegions/SegExp.h>
43 #include <LocalRegions/QuadExp.h>
44 #include <LocalRegions/TriExp.h>
47 
48 namespace Nektar
49 {
50  namespace Utilities
51  {
52  ModuleKey ProcessPerAlign::className =
54  ModuleKey(eProcessModule, "peralign"),
55  ProcessPerAlign::create);
56 
57  /**
58  * @class ProcessPerAlign
59  */
60 
61  /**
62  * @brief Default constructor.
63  */
64  ProcessPerAlign::ProcessPerAlign(MeshSharedPtr m) : ProcessModule(m)
65  {
66  m_config["surf1"] = ConfigOption(false, "-1",
67  "Tag identifying first surface.");
68  m_config["surf2"] = ConfigOption(false, "-1",
69  "Tag identifying first surface.");
70  m_config["dir"] = ConfigOption(false, "",
71  "Direction in which to align (either x, y, or z)");
72  m_config["orient"] = ConfigOption(true, "0",
73  "Attempt to reorient tets and prisms");
74  }
75 
76  /**
77  * @brief Destructor.
78  */
80  {
81 
82  }
83 
85  {
86  int surf1 = m_config["surf1"]. as<int> ();
87  int surf2 = m_config["surf2"]. as<int> ();
88  string dir = m_config["dir"]. as<string>();
89  bool orient = m_config["orient"].as<bool> ();
90 
91  if (surf1 == -1)
92  {
93  cerr << "WARNING: surf1 must be set to a positive integer. "
94  << "Skipping periodic alignment." << endl;
95  return;
96  }
97 
98  if (surf2 == -1)
99  {
100  cerr << "WARNING: surf2 must be set to a positive integer. "
101  << "Skipping periodic alignment." << endl;
102  return;
103  }
104 
105  if (dir != "x" && dir != "y" && dir != "z")
106  {
107  cerr << "WARNING: dir must be set to either x, y or z. "
108  << "Skipping periodic alignment." << endl;
109  return;
110  }
111 
112  NekDouble vec[3];
113  vec[0] = dir == "x" ? 1.0 : 0.0;
114  vec[1] = dir == "y" ? 1.0 : 0.0;
115  vec[2] = dir == "z" ? 1.0 : 0.0;
116 
117  CompositeMap::iterator it1 = m_mesh->m_composite.find(surf1);
118  CompositeMap::iterator it2 = m_mesh->m_composite.find(surf2);
119 
120  if (it1 == m_mesh->m_composite.end())
121  {
122  cerr << "WARNING: Couldn't find surface " << surf1
123  << ". Skipping periodic alignment." << endl;
124  return;
125  }
126 
127  if (it2 == m_mesh->m_composite.end())
128  {
129  cerr << "WARNING: Couldn't find surface " << surf2 << ", "
130  << "skipping periodic alignment." << endl;
131  return;
132  }
133 
134  CompositeSharedPtr c1 = it1->second;
135  CompositeSharedPtr c2 = it2->second;
136 
137  if (c1->m_items.size() != c2->m_items.size())
138  {
139  cerr << "WARNING: Surfaces " << surf1 << " and " << surf2
140  << " have different numbers of elements. Skipping periodic"
141  << " alignment." << endl;
142  return;
143  }
144 
145  c1->m_reorder = false;
146  c2->m_reorder = false;
147 
148  map<int, pair<FaceSharedPtr, vector<int> > > perFaces;
149 
150  // Loop over elements, calculate centroids of elements in c2.
151  map<int, Node> centroidMap;
153  for (int i = 0; i < c2->m_items.size(); ++i)
154  {
155  Node centroid;
156  for (int j = 0; j < c2->m_items[i]->GetVertexCount(); ++j)
157  {
158  centroid += *(c2->m_items[i]->GetVertex(j));
159  }
160  centroid /= (NekDouble)c2->m_items[i]->GetVertexCount();
161  centroidMap[i] = centroid;
162  }
163 
164  boost::unordered_set<int> elmtDone;
165  map<int, int> elmtPairs;
166  map<int, int> vertCheck;
167 
168  for (int i = 0; i < c1->m_items.size(); ++i)
169  {
170  Node centroid;
171  for (int j = 0; j < c1->m_items[i]->GetVertexCount(); ++j)
172  {
173  centroid += *(c1->m_items[i]->GetVertex(j));
174  }
175  centroid /= (NekDouble)c1->m_items[i]->GetVertexCount();
176 
177  for (it = centroidMap.begin(); it != centroidMap.end(); ++it)
178  {
179  if (elmtDone.count(it->first) > 0)
180  {
181  continue;
182  }
183 
184  Node dx = it->second - centroid;
185  if (fabs(fabs(dx.m_x*vec[0] + dx.m_y*vec[1] + dx.m_z*vec[2])/
186  sqrt(dx.abs2()) - 1.0) < 1e-8)
187  {
188  // Found match
189  int id1, id2;
190 
191  if (c1->m_items[i]->GetConf().m_e == LibUtilities::eSegment)
192  {
193  id1 = c1->m_items[i] ->GetEdgeLink()->m_id;
194  id2 = c2->m_items[it->first]->GetEdgeLink()->m_id;
195  }
196  else
197  {
198  id1 = c1->m_items[i] ->GetFaceLink()->m_id;
199  id2 = c2->m_items[it->first]->GetFaceLink()->m_id;
200  }
201 
202  elmtDone.insert(it->first);
203  elmtPairs[i] = it->first;
204 
205  // Identify periodic vertices
206  int nVerts = c1->m_items[i]->GetVertexCount();
207  vector<int> perVerts(nVerts, 0), perVertsInv(nVerts, 0);
208 
209  if (orient)
210  {
211  for (int k = 0; k < nVerts; ++k)
212  {
213  NodeSharedPtr n1 = c1->m_items[i]->GetFaceLink()->m_vertexList[k];
214  int l;
215 
216  for (l = 0; l < nVerts; ++l)
217  {
218  NodeSharedPtr n2 =
219  c2->m_items[it->first]->GetFaceLink()->m_vertexList[l];
220 
221  Node dn = *n2 - *n1;
222  if (fabs(fabs(dn.m_x*vec[0] + dn.m_y*vec[1] +
223  dn.m_z*vec[2])/
224  sqrt(dn.abs2()) - 1.0) < 1e-8)
225  {
226  perVerts [k] = l;
227  perVertsInv[l] = k;
228 
229  int id1 = n1->m_id;
230  int id2 = n2->m_id;
231  if (vertCheck.count(id1) == 0)
232  {
233  vertCheck[id1] = id2;
234  }
235  else
236  {
237  ASSERTL0(vertCheck[id1] == id2,
238  "Periodic vertex already "
239  "identified!");
240  }
241  break;
242  }
243  }
244  ASSERTL1(l < nVerts,
245  "Could not identify periodic vertices.");
246  }
247 
248  int tot1 = 0, tot2 = 0;
249  for (int k = 0; k < nVerts; ++k)
250  {
251  tot1 += perVerts [k];
252  tot2 += perVertsInv[k];
253  }
254  ASSERTL0(tot1 == nVerts*(nVerts-1)/2 &&
255  tot2 == nVerts*(nVerts-1)/2,
256  "Error identifying periodic vertices");
257  }
258 
259  if (c2->m_items[i]->GetConf().m_e != LibUtilities::eSegment)
260  {
261  perFaces[id1] = make_pair(
262  c2->m_items[it->first]->GetFaceLink(), perVerts);
263  perFaces[id2] = make_pair(
264  c1->m_items[i] ->GetFaceLink(), perVertsInv);
265  }
266  break;
267  }
268  }
269 
270  if (it == centroidMap.end())
271  {
272  cerr << "WARNING: Could not find matching edge for surface "
273  << "element " << c1->m_items[i]->GetId() << ". "
274  << "Skipping periodic alignment." << endl;
275  return;
276  }
277  }
278 
279  // Reorder vectors.
280  vector<ElementSharedPtr> tmp = c2->m_items;
281 
283 
284  for (int i = 0; i < tmp.size(); ++i)
285  {
286  c2->m_items[i] = tmp[elmtPairs[i]];
287  }
288 
289  if (orient)
290  {
291  ReorderPrisms(perFaces);
292  }
293  }
294  }
295 }