Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessExtractSurf.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessExtractSurf.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: Extract one or more surfaces from mesh.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "MeshElements.h"
37 #include "ProcessExtractSurf.h"
38 
41 
42 #include <vector>
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace Utilities
48  {
49  ModuleKey ProcessExtractSurf::className =
51  ModuleKey(eProcessModule, "extract"), ProcessExtractSurf::create,
52  "Process elements to extract a specified surface(s) or composites(s).");
53 
54  ProcessExtractSurf::ProcessExtractSurf(MeshSharedPtr m) : ProcessModule(m)
55  {
56  m_config["surf"] = ConfigOption(false, "-1",
57  "Tag identifying surface/composite to process.");
58  }
59 
61  {
62  }
63 
65  {
66  int i, j;
67  string surf = m_config["surf"].as<string>();
68 
69  // Obtain vector of surface IDs from string.
70  vector<unsigned int> surfs;
71  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
72  sort(surfs.begin(), surfs.end());
73 
74  // If we're running in verbose mode print out a list of surfaces.
75  if (m_mesh->m_verbose)
76  {
77  cout << "ProcessExtractSurf: extracting surface"
78  << (surfs.size() > 1 ? "s" : "") << " " << surf << endl;
79  }
80 
81  // Make a copy of all existing elements of one dimension lower.
82  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim-1];
83 
84  // Clear all elements.
85  m_mesh->m_element[m_mesh->m_expDim] .clear();
86  m_mesh->m_element[m_mesh->m_expDim-1].clear();
87 
88  // Clear existing vertices, edges and faces.
89  m_mesh->m_vertexSet.clear();
90  m_mesh->m_edgeSet.clear();
91  m_mesh->m_faceSet.clear();
92 
93  // Clear all edge -> element links.
94  for (i = 0; i < el.size(); ++i)
95  {
96  vector<EdgeSharedPtr> edges = el[i]->GetEdgeList();
97  for (j = 0; j < edges.size(); ++j)
98  {
99  edges[j]->m_elLink.clear();
100  }
101 
102  FaceSharedPtr f = el[i]->GetFaceLink();
103  if (f)
104  {
105  for (j = 0; j < f->m_edgeList.size(); ++j)
106  {
107  f->m_edgeList[j]->m_elLink.clear();
108  }
109  }
110  }
111 
112  // keptIds stores IDs of elements we processed earlier.
113  boost::unordered_set<int> keptIds;
114 
115  // Iterate over list of surface elements.
116  for (i = 0; i < el.size(); ++i)
117  {
118  // Work out whether this lies on our surface of interest.
119  vector<int> inter, tags = el[i]->GetTagList();
120 
121  sort(tags.begin(), tags.end());
122  set_intersection(surfs.begin(), surfs.end(),
123  tags .begin(), tags .end(),
124  back_inserter(inter));
125 
126  // It doesn't continue to next element.
127  if (inter.size() != 1)
128  {
129  continue;
130  }
131 
132  // Get list of element vertices and edges.
133  ElementSharedPtr elmt = el[i];
134  vector<NodeSharedPtr> verts = elmt->GetVertexList();
135  vector<EdgeSharedPtr> edges = elmt->GetEdgeList();
136 
137  // Insert surface vertices.
138  for (j = 0; j < verts.size(); ++j)
139  {
140  m_mesh->m_vertexSet.insert(verts[j]);
141  }
142 
143  // Problem: edges and element IDs aren't enumerated with
144  // geometry IDs by some input modules/the Module ProcessEdges
145  // function. Get around this by replacing everything in the
146  // edge/face with information from edge/face link.
147  EdgeSharedPtr e = elmt->GetEdgeLink();
148  FaceSharedPtr f = elmt->GetFaceLink();
149  if (e)
150  {
151  elmt->SetId(e->m_id);
152  }
153  else if (f)
154  {
155  for (j = 0; j < edges.size(); ++j)
156  {
157  m_mesh->m_edgeSet.insert(f->m_edgeList[j]);
158  elmt->SetEdge(j, f->m_edgeList[j]);
159  f->m_edgeList[j]->m_elLink.push_back(
160  std::make_pair(elmt, j));
161  }
162  elmt->SetId(f->m_id);
163  }
164  else
165  {
166  for (j = 0; j < edges.size(); ++j)
167  {
168  m_mesh->m_edgeSet.insert(edges[j]);
169  }
170  }
171 
172  // Nullify edge/face links to get correct tag
173  elmt->SetFaceLink(FaceSharedPtr());
174  elmt->SetEdgeLink(EdgeSharedPtr());
175  keptIds.insert(elmt->GetId());
176 
177  // Push element back into the list.
178  m_mesh->m_element[m_mesh->m_expDim-1].push_back(elmt);
179  }
180 
181  // Decrement the expansion dimension to get manifold embedding.
182  m_mesh->m_expDim--;
183 
184  // Now process composites. This is necessary because 2D surfaces may
185  // contain both quadrilaterals and triangles and so need to be split
186  // up.
187  CompositeMap tmp = m_mesh->m_composite;
189 
190  m_mesh->m_composite.clear();
191  int maxId = -1;
192 
193  // Loop over composites for first time to determine any composites
194  // which don't have elements of the correct dimension.
195  for (it = tmp.begin(); it != tmp.end(); ++it)
196  {
197  if (it->second->m_items[0]->GetDim() != m_mesh->m_expDim)
198  {
199  continue;
200  }
201 
202  vector<ElementSharedPtr> el = it->second->m_items;
203  it->second->m_items.clear();
204 
205  for (i = 0; i < el.size(); ++i)
206  {
207  if (keptIds.count(el[i]->GetId()) > 0)
208  {
209  it->second->m_items.push_back(el[i]);
210  }
211  }
212 
213  if (it->second->m_items.size() == 0)
214  {
215  continue;
216  }
217 
218  m_mesh->m_composite.insert(*it);
219 
220  // Figure out the maximum ID so if we need to create new
221  // composites we can give them a unique ID.
222  maxId = (std::max)(maxId, (int)it->second->m_id) + 1;
223  }
224 
225  tmp = m_mesh->m_composite;
226  m_mesh->m_composite.clear();
228 
229  // Now do another loop over the composites to remove composites
230  // which don't contain any elements in the new mesh.
231  for (it = tmp.begin(); it != tmp.end(); ++it)
232  {
233  CompositeSharedPtr c = it->second;
234  vector<ElementSharedPtr> el = c->m_items;
235 
236  // Remove all but the first element from this composite.
237  string initialTag = el[0]->GetTag();
238  c->m_items.resize(1);
239  c->m_tag = initialTag;
240 
241  // newComps stores the new composites. The key is the composite
242  // type (e.g. Q for quad) and value is the composite.
243  map<string, CompositeSharedPtr> newComps;
244  newComps[initialTag] = c;
245 
246  // Loop over remaining elements in composite and figure out
247  // whether it needs to be split up.
248  for (i = 1; i < el.size(); ++i)
249  {
250  // See if tag exists. If it does, we append this to the
251  // composite, otherwise we create a new composite and store
252  // it in newComps.
253  string tag = el[i]->GetTag();
254  it2 = newComps.find(tag);
255  if (it2 == newComps.end())
256  {
257  CompositeSharedPtr newComp(new Composite());
258  newComp->m_id = maxId++;
259  newComp->m_tag = tag;
260  newComp->m_items.push_back(el[i]);
261  newComps[tag] = newComp;
262  }
263  else
264  {
265  it2->second->m_items.push_back(el[i]);
266  }
267  }
268 
269  // Print out mapping information if we remapped composite IDs.
270  if (m_mesh->m_verbose && newComps.size() > 1)
271  {
272  cout << "- Mapping composite " << it->first << " ->";
273  }
274 
275  // Insert new composites.
276  for (i = 0, it2 = newComps.begin(); it2 != newComps.end();
277  ++it2, ++i)
278  {
279  if (m_mesh->m_verbose && newComps.size() > 1)
280  {
281  cout << (i > 0 ? ", " : " ") << it2->second->m_id << "("
282  << it2->second->m_tag << ")";
283  }
284  m_mesh->m_composite[it2->second->m_id] = it2->second;
285  }
286 
287  if (m_mesh->m_verbose && newComps.size() > 1)
288  {
289  cout << endl;
290  }
291  }
292  }
293  }
294 }