ProcessPerAlign.cpp 10.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
////////////////////////////////////////////////////////////////////////////////
//
//  File: ProcessPerAlign.cpp
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: Reorder composites to align periodic boundaries.
//
////////////////////////////////////////////////////////////////////////////////

36
#include <LibUtilities/BasicUtils/SharedArray.hpp>
37
38
39
40
41

#include <LocalRegions/SegExp.h>
#include <LocalRegions/QuadExp.h>
#include <LocalRegions/TriExp.h>
#include <LocalRegions/NodalTriExp.h>
42

43
#include <NekMeshUtils/MeshElements/Element.h>
44

45
46
#include "ProcessPerAlign.h"

Julian Marcon's avatar
Julian Marcon committed
47
48
#include <boost/algorithm/string.hpp>

49
using namespace std;
50
using namespace Nektar::NekMeshUtils;
51

52
53
namespace Nektar
{
54
55
namespace Utilities
{
Dave Moxey's avatar
Dave Moxey committed
56

57
58
ModuleKey ProcessPerAlign::className =
    GetModuleFactory().RegisterCreatorFunction(
59
        ModuleKey(eProcessModule, "peralign"), ProcessPerAlign::create);
Dave Moxey's avatar
Dave Moxey committed
60

61
62
63
/**
 * @class ProcessPerAlign
 */
64

65
66
67
68
69
/**
 * @brief Default constructor.
 */
ProcessPerAlign::ProcessPerAlign(MeshSharedPtr m) : ProcessModule(m)
{
70
71
72
73
74
    m_config["surf1"] =
        ConfigOption(false, "-1", "Tag identifying first surface.");
    m_config["surf2"] =
        ConfigOption(false, "-1", "Tag identifying first surface.");
    m_config["dir"] = ConfigOption(
Julian Marcon's avatar
Julian Marcon committed
75
76
        false, "", "Direction in which to align (either x, y, or z; "
                   "or vector with components separated by a comma)");
77
78
    m_config["orient"] =
        ConfigOption(true, "0", "Attempt to reorient tets and prisms");
79
}
80

81
82
83
84
85
86
/**
 * @brief Destructor.
 */
ProcessPerAlign::~ProcessPerAlign()
{
}
87

88
89
void ProcessPerAlign::Process()
{
90
91
92
93
    int surf1   = m_config["surf1"].as<int>();
    int surf2   = m_config["surf2"].as<int>();
    string dir  = m_config["dir"].as<string>();
    bool orient = m_config["orient"].as<bool>();
94

95
96
97
98
99
100
    if (surf1 == -1)
    {
        cerr << "WARNING: surf1 must be set to a positive integer. "
             << "Skipping periodic alignment." << endl;
        return;
    }
101

102
103
104
105
106
107
    if (surf2 == -1)
    {
        cerr << "WARNING: surf2 must be set to a positive integer. "
             << "Skipping periodic alignment." << endl;
        return;
    }
108

Julian Marcon's avatar
Julian Marcon committed
109
110
    vector<string> tmp1;
    boost::split(tmp1, dir, boost::is_any_of(","));
111

112
    NekDouble vec[3];
Julian Marcon's avatar
Julian Marcon committed
113
114
115

    if (tmp1.size() == 1)
    {
116
        if (!dir.size() && m_mesh->m_spaceDim == 2)
Julian Marcon's avatar
Julian Marcon committed
117
        {
118
119
120
            Array<OneD, NekDouble> T =
                m_mesh->m_cad->GetPeriodicTranslationVector(surf1, surf2);
            NekDouble mag = sqrt(T[0] * T[0] + T[1] * T[1]);
Julian Marcon's avatar
Julian Marcon committed
121

122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
            vec[0] = T[0] / mag;
            vec[1] = T[1] / mag;
            vec[2] = T[2] / mag;
        }
        else
        {
            if (dir != "x" && dir != "y" && dir != "z")
            {
                cerr << "WARNING: dir must be set to either x, y or z. "
                    << "Skipping periodic alignment." << endl;
                return;
            }

            vec[0] = dir == "x" ? 1.0 : 0.0;
            vec[1] = dir == "y" ? 1.0 : 0.0;
            vec[2] = dir == "z" ? 1.0 : 0.0;
        }
Julian Marcon's avatar
Julian Marcon committed
139
140
141
142
143
144
145
146
147
148
149
    }
    else if (tmp1.size() == 3)
    {
        vec[0] = boost::lexical_cast<NekDouble>(tmp1[0]);
        vec[1] = boost::lexical_cast<NekDouble>(tmp1[1]);
        vec[2] = boost::lexical_cast<NekDouble>(tmp1[2]);
    }
    else
    {
        ASSERTL0(false,"expected three components or letter for direction");
    }
150

151
152
    CompositeMap::iterator it1 = m_mesh->m_composite.find(surf1);
    CompositeMap::iterator it2 = m_mesh->m_composite.find(surf2);
Dave Moxey's avatar
Dave Moxey committed
153

154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
    if (it1 == m_mesh->m_composite.end())
    {
        cerr << "WARNING: Couldn't find surface " << surf1
             << ". Skipping periodic alignment." << endl;
        return;
    }

    if (it2 == m_mesh->m_composite.end())
    {
        cerr << "WARNING: Couldn't find surface " << surf2 << ", "
             << "skipping periodic alignment." << endl;
        return;
    }

    CompositeSharedPtr c1 = it1->second;
    CompositeSharedPtr c2 = it2->second;

    if (c1->m_items.size() != c2->m_items.size())
    {
        cerr << "WARNING: Surfaces " << surf1 << " and " << surf2
             << " have different numbers of elements. Skipping periodic"
             << " alignment." << endl;
        return;
    }

    c1->m_reorder = false;
    c2->m_reorder = false;

    map<int, pair<FaceSharedPtr, vector<int> > > perFaces;

    // Loop over elements, calculate centroids of elements in c2.
    map<int, Node> centroidMap;
    map<int, Node>::iterator it;
    for (int i = 0; i < c2->m_items.size(); ++i)
    {
        Node centroid;
        for (int j = 0; j < c2->m_items[i]->GetVertexCount(); ++j)
        {
            centroid += *(c2->m_items[i]->GetVertex(j));
        }
        centroid /= (NekDouble)c2->m_items[i]->GetVertexCount();
        centroidMap[i] = centroid;
    }
197

198
199
200
    boost::unordered_set<int> elmtDone;
    map<int, int> elmtPairs;
    map<int, int> vertCheck;
Dave Moxey's avatar
Dave Moxey committed
201

202
203
204
205
206
207
208
209
    for (int i = 0; i < c1->m_items.size(); ++i)
    {
        Node centroid;
        for (int j = 0; j < c1->m_items[i]->GetVertexCount(); ++j)
        {
            centroid += *(c1->m_items[i]->GetVertex(j));
        }
        centroid /= (NekDouble)c1->m_items[i]->GetVertexCount();
210

211
212
213
        for (it = centroidMap.begin(); it != centroidMap.end(); ++it)
        {
            if (elmtDone.count(it->first) > 0)
214
            {
215
                continue;
216
217
            }

218
            Node dx = it->second - centroid;
219
220
221
            if (fabs(fabs(dx.m_x * vec[0] + dx.m_y * vec[1] + dx.m_z * vec[2]) /
                         sqrt(dx.abs2()) -
                     1.0) < 1e-8)
222
            {
223
224
225
226
                // Found match
                int id1, id2;

                if (c1->m_items[i]->GetConf().m_e == LibUtilities::eSegment)
227
                {
228
                    id1 = c1->m_items[i]->GetEdgeLink()->m_id;
229
                    id2 = c2->m_items[it->first]->GetEdgeLink()->m_id;
230
                }
231
                else
232
                {
233
                    id1 = c1->m_items[i]->GetFaceLink()->m_id;
234
235
                    id2 = c2->m_items[it->first]->GetFaceLink()->m_id;
                }
236

237
238
                elmtDone.insert(it->first);
                elmtPairs[i] = it->first;
239

240
241
242
                // Identify periodic vertices
                int nVerts = c1->m_items[i]->GetVertexCount();
                vector<int> perVerts(nVerts, 0), perVertsInv(nVerts, 0);
243

244
245
246
247
                if (orient)
                {
                    for (int k = 0; k < nVerts; ++k)
                    {
248
249
                        NodeSharedPtr n1 =
                            c1->m_items[i]->GetFaceLink()->m_vertexList[k];
250
                        int l;
251

252
                        for (l = 0; l < nVerts; ++l)
253
                        {
254
255
256
                            NodeSharedPtr n2 = c2->m_items[it->first]
                                                   ->GetFaceLink()
                                                   ->m_vertexList[l];
257
258

                            Node dn = *n2 - *n1;
259
260
261
262
                            if (fabs(fabs(dn.m_x * vec[0] + dn.m_y * vec[1] +
                                          dn.m_z * vec[2]) /
                                         sqrt(dn.abs2()) -
                                     1.0) < 1e-8)
263
                            {
264
                                perVerts[k]    = l;
265
                                perVertsInv[l] = k;
266

267
268
269
                                int id1 = n1->m_id;
                                int id2 = n2->m_id;
                                if (vertCheck.count(id1) == 0)
270
                                {
271
                                    vertCheck[id1] = id2;
272
                                }
273
274
275
276
277
278
279
                                else
                                {
                                    ASSERTL0(vertCheck[id1] == id2,
                                             "Periodic vertex already "
                                             "identified!");
                                }
                                break;
280
                            }
281
                        }
282
283
284
                        ASSERTL1(l < nVerts,
                                 "Could not identify periodic vertices.");
                    }
285

286
287
288
                    int tot1 = 0, tot2 = 0;
                    for (int k = 0; k < nVerts; ++k)
                    {
289
                        tot1 += perVerts[k];
290
                        tot2 += perVertsInv[k];
291
                    }
292
293
                    ASSERTL0(tot1 == nVerts * (nVerts - 1) / 2 &&
                                 tot2 == nVerts * (nVerts - 1) / 2,
294
                             "Error identifying periodic vertices");
295
296
                }

297
                if (c2->m_items[i]->GetConf().m_e != LibUtilities::eSegment)
298
                {
299
300
                    perFaces[id1] = make_pair(
                        c2->m_items[it->first]->GetFaceLink(), perVerts);
301
302
                    perFaces[id2] =
                        make_pair(c1->m_items[i]->GetFaceLink(), perVertsInv);
303
                }
304
                break;
305
            }
306
        }
307

308
309
310
311
312
313
314
315
        if (it == centroidMap.end())
        {
            cerr << "WARNING: Could not find matching edge for surface "
                 << "element " << c1->m_items[i]->GetId() << ". "
                 << "Skipping periodic alignment." << endl;
            return;
        }
    }
316

317
318
    // Reorder vectors.
    vector<ElementSharedPtr> tmp = c2->m_items;
319

320
    map<int, int>::iterator mIt;
321

322
323
324
    for (int i = 0; i < tmp.size(); ++i)
    {
        c2->m_items[i] = tmp[elmtPairs[i]];
325
    }
326
327
328
329
330
331
332

    if (orient)
    {
        ReorderPrisms(perFaces);
    }
}
}
333
}