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