Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 
38 #include "ProcessExtractSurf.h"
39 
40 using namespace std;
41 using namespace Nektar::NekMeshUtils;
42 
43 namespace Nektar
44 {
45 namespace Utilities
46 {
47 
48 ModuleKey ProcessExtractSurf::className =
50  ModuleKey(eProcessModule, "extract"),
51  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(
57  false, "NotSet", "Tag identifying surface/composite to process.");
58  m_config["detectbnd"] =
59  ConfigOption(false, "-1", "Tag to detect on boundary composites");
60 }
61 
63 {
64 }
65 
67 {
68  int i, j;
69  string surf = m_config["surf"].as<string>();
70  bool detectbnd = m_config["detectbnd"].beenSet;
71 
72  // Obtain vector of surface IDs from string.
73  vector<unsigned int> surfs;
74  ASSERTL0(ParseUtils::GenerateSeqVector(surf.c_str(), surfs),
75  "Failed to interp surf string. Have you specified this string?");
76  sort(surfs.begin(), surfs.end());
77 
78  // If we're running in verbose mode print out a list of surfaces.
79  if (m_mesh->m_verbose)
80  {
81  cout << "ProcessExtractSurf: extracting surface"
82  << (surfs.size() > 1 ? "s" : "") << " " << surf << endl;
83  }
84 
85  // Make a copy of all existing elements of one dimension lower.
86  vector<ElementSharedPtr> el = m_mesh->m_element[m_mesh->m_expDim - 1];
87 
88  // Clear all elements.
89  m_mesh->m_element[m_mesh->m_expDim].clear();
90  m_mesh->m_element[m_mesh->m_expDim - 1].clear();
91 
92  // Clear existing vertices, edges and faces.
93  m_mesh->m_vertexSet.clear();
94  m_mesh->m_edgeSet.clear();
95  m_mesh->m_faceSet.clear();
96 
97  // Clear all edge -> element links.
98  for (i = 0; i < el.size(); ++i)
99  {
100  vector<EdgeSharedPtr> edges = el[i]->GetEdgeList();
101  for (j = 0; j < edges.size(); ++j)
102  {
103  edges[j]->m_elLink.clear();
104  }
105 
106  FaceSharedPtr f = el[i]->GetFaceLink();
107  if (f)
108  {
109  for (j = 0; j < f->m_edgeList.size(); ++j)
110  {
111  f->m_edgeList[j]->m_elLink.clear();
112  }
113  }
114  }
115 
116  // keptIds stores IDs of elements we processed earlier.
117  boost::unordered_set<int> keptIds;
118 
119  EdgeSet bndEdgeSet;
120 
121  // Iterate over list of surface elements.
122  for (i = 0; i < el.size(); ++i)
123  {
124  // Work out whether this lies on our surface of interest.
125  vector<int> inter, tags = el[i]->GetTagList();
126 
127  sort(tags.begin(), tags.end());
128  set_intersection(surfs.begin(),
129  surfs.end(),
130  tags.begin(),
131  tags.end(),
132  back_inserter(inter));
133 
134  // It doesn't continue to next element.
135  if (inter.size() != 1)
136  {
137  continue;
138  }
139 
140  // Get list of element vertices and edges.
141  ElementSharedPtr elmt = el[i];
142  vector<NodeSharedPtr> verts = elmt->GetVertexList();
143  vector<EdgeSharedPtr> edges = elmt->GetEdgeList();
144 
145  // Insert surface vertices.
146  for (j = 0; j < verts.size(); ++j)
147  {
148  m_mesh->m_vertexSet.insert(verts[j]);
149  }
150 
151  // Problem: edges and element IDs aren't enumerated with
152  // geometry IDs by some input modules/the Module ProcessEdges
153  // function. Get around this by replacing everything in the
154  // edge/face with information from edge/face link.
155  EdgeSharedPtr e = elmt->GetEdgeLink();
156  FaceSharedPtr f = elmt->GetFaceLink();
157  if (e)
158  {
159  elmt->SetId(e->m_id);
160  }
161  else if (f)
162  {
163  for (j = 0; j < f->m_vertexList.size(); j++)
164  {
165  elmt->SetVertex(j, f->m_vertexList[j]);
166  }
167 
168  EdgeSet::iterator edit;
169  for (j = 0; j < edges.size(); ++j)
170  {
171  m_mesh->m_edgeSet.insert(f->m_edgeList[j]);
172  elmt->SetEdge(j, f->m_edgeList[j]);
173  f->m_edgeList[j]->m_elLink.push_back(std::make_pair(elmt, j));
174 
175  // generate a list of edges on boundary of surfaces being
176  // extracted
177  edit = bndEdgeSet.find(f->m_edgeList[j]);
178  if (edit != bndEdgeSet.end())
179  {
180  // remove since visited more than once
181  bndEdgeSet.erase(edit);
182  }
183  else
184  {
185  bndEdgeSet.insert(f->m_edgeList[j]);
186  }
187  }
188  elmt->SetVolumeNodes(f->m_faceNodes);
189  elmt->SetId(f->m_id);
190  elmt->SetCurveType(f->m_curveType);
191  }
192  else
193  {
194  for (j = 0; j < edges.size(); ++j)
195  {
196  m_mesh->m_edgeSet.insert(edges[j]);
197  }
198  }
199 
200  // Nullify edge/face links to get correct tag
201  elmt->SetFaceLink(FaceSharedPtr());
202  elmt->SetEdgeLink(EdgeSharedPtr());
203  keptIds.insert(elmt->GetId());
204 
205  // Push element back into the list.
206  m_mesh->m_element[m_mesh->m_expDim - 1].push_back(elmt);
207  }
208 
209  // Decrement the expansion dimension to get manifold embedding.
210  m_mesh->m_expDim--;
211 
212  // Now process composites. This is necessary because 2D surfaces may
213  // contain both quadrilaterals and triangles and so need to be split
214  // up.
215  CompositeMap tmp = m_mesh->m_composite;
217 
218  m_mesh->m_composite.clear();
219  int maxId = -1;
220 
221  // Loop over composites for first time to determine any composites
222  // which don't have elements of the correct dimension.
223  for (it = tmp.begin(); it != tmp.end(); ++it)
224  {
225  if (it->second->m_items[0]->GetDim() != m_mesh->m_expDim)
226  {
227  continue;
228  }
229 
230  vector<ElementSharedPtr> el = it->second->m_items;
231  it->second->m_items.clear();
232 
233  for (i = 0; i < el.size(); ++i)
234  {
235  if (keptIds.count(el[i]->GetId()) > 0)
236  {
237  it->second->m_items.push_back(el[i]);
238  }
239  }
240 
241  if (it->second->m_items.size() == 0)
242  {
243  continue;
244  }
245 
246  m_mesh->m_composite.insert(*it);
247 
248  // Figure out the maximum ID so if we need to create new
249  // composites we can give them a unique ID.
250  maxId = (std::max)(maxId, (int)it->second->m_id) + 1;
251  }
252 
253  tmp = m_mesh->m_composite;
254  m_mesh->m_composite.clear();
256 
257  // Now do another loop over the composites to remove composites
258  // which don't contain any elements in the new mesh.
259  for (it = tmp.begin(); it != tmp.end(); ++it)
260  {
261  CompositeSharedPtr c = it->second;
262  vector<ElementSharedPtr> el = c->m_items;
263 
264  // Remove all but the first element from this composite.
265  string initialTag = el[0]->GetTag();
266  c->m_items.resize(1);
267  c->m_tag = initialTag;
268 
269  // newComps stores the new composites. The key is the composite
270  // type (e.g. Q for quad) and value is the composite.
271  map<string, CompositeSharedPtr> newComps;
272  newComps[initialTag] = c;
273 
274  // Loop over remaining elements in composite and figure out
275  // whether it needs to be split up.
276  for (i = 1; i < el.size(); ++i)
277  {
278  // See if tag exists. If it does, we append this to the
279  // composite, otherwise we create a new composite and store
280  // it in newComps.
281  string tag = el[i]->GetTag();
282  it2 = newComps.find(tag);
283  if (it2 == newComps.end())
284  {
285  CompositeSharedPtr newComp(new Composite());
286  newComp->m_id = maxId++;
287  newComp->m_tag = tag;
288  newComp->m_items.push_back(el[i]);
289  newComps[tag] = newComp;
290  }
291  else
292  {
293  it2->second->m_items.push_back(el[i]);
294  }
295  }
296 
297  // Print out mapping information if we remapped composite IDs.
298  if (m_mesh->m_verbose && newComps.size() > 1)
299  {
300  cout << "- Mapping composite " << it->first << " ->";
301  }
302 
303  // Insert new composites.
304  for (i = 0, it2 = newComps.begin(); it2 != newComps.end(); ++it2, ++i)
305  {
306  if (m_mesh->m_verbose && newComps.size() > 1)
307  {
308  cout << (i > 0 ? ", " : " ") << it2->second->m_id << "("
309  << it2->second->m_tag << ")";
310  }
311  m_mesh->m_composite[it2->second->m_id] = it2->second;
312  }
313 
314  if (m_mesh->m_verbose && newComps.size() > 1)
315  {
316  cout << endl;
317  }
318  }
319 
320  // Detect composites for boundaries. This is done by looping over all
321  // elements identifiying if they are not part of required surfaces and if
322  // not setting up a list of boundary edges (identified by only being visited
323  // once). This list is then compared against an earlier identification of
324  // boundary edges on the required surfaces and if the two overlap add a
325  // segment element and put segment element in composite as well
326  if (detectbnd)
327  {
328  if (m_mesh->m_expDim != 2)
329  {
330  cerr << "Surface boundary detection only implemented for 2D meshes"
331  << endl;
332  return;
333  }
334 
335  map<int, EdgeSet> surfBndEdgeSet;
336  EdgeSet::iterator edit;
337  map<int, string> surfLabels;
338 
339  // Iterate over list of surface elements.
340  for (i = 0; i < el.size(); ++i)
341  {
342  // Work out whether this lies on our surface of interest.
343  vector<int> inter, tags = el[i]->GetTagList();
344 
345  ASSERTL0(tags.size() == 1, "Not sure what mutliple tags implies");
346 
347  sort(tags.begin(), tags.end());
348  set_intersection(surfs.begin(),
349  surfs.end(),
350  tags.begin(),
351  tags.end(),
352  back_inserter(inter));
353 
354  // It does so continue to next element.
355  if (inter.size() == 1)
356  {
357  continue;
358  }
359 
360  int surf = tags[0];
361 
362  // gather surface labels if they exist.
363  if (m_mesh->m_faceLabels.count(surf))
364  {
365  surfLabels[surf] = m_mesh->m_faceLabels[surf];
366  }
367 
368  // Get list of element vertices and edges.
369  ElementSharedPtr elmt = el[i];
370  vector<EdgeSharedPtr> edges = elmt->GetEdgeList();
371 
372  FaceSharedPtr f = elmt->GetFaceLink();
373  if (f)
374  {
375  for (j = 0; j < edges.size(); ++j)
376  {
377  // generate a list of edges on boundary of surfaces being
378  // extracted
379  if (surfBndEdgeSet.count(surf))
380  {
381  edit = surfBndEdgeSet[surf].find(f->m_edgeList[j]);
382  if (edit != surfBndEdgeSet[surf].end())
383  {
384  // remove since visited more than once
385  surfBndEdgeSet[surf].erase(edit);
386  }
387  else
388  {
389  surfBndEdgeSet[surf].insert(f->m_edgeList[j]);
390  }
391  }
392  else
393  {
394  EdgeSet newEdgeSet;
395  surfBndEdgeSet[surf] = newEdgeSet;
396  surfBndEdgeSet[surf].insert(f->m_edgeList[j]);
397  }
398  }
399  }
400  }
401 
402  m_mesh->m_faceLabels.clear();
403 
404  // iteratve over surfBndEdgeSet and see if they are in BndEdgeSet
406  for (esetit = surfBndEdgeSet.begin(); esetit != surfBndEdgeSet.end();
407  ++esetit)
408  {
409  CompositeSharedPtr newComp(new Composite());
410  newComp->m_id = maxId;
411  newComp->m_tag = "E";
412  // set up labels if they exist
413  if (surfLabels.count(esetit->first))
414  {
415  newComp->m_label = surfLabels[esetit->first];
416  }
417 
418  for (edit = esetit->second.begin(); edit != esetit->second.end();
419  ++edit)
420  {
421  EdgeSet::iterator locit;
422  if ((locit = bndEdgeSet.find(*edit)) != bndEdgeSet.end())
423  {
424  // make 1D segment element
426 
427  vector<int> tags;
428  tags.push_back(maxId);
429 
430  // make unique node list
431  vector<NodeSharedPtr> nodeList;
432  nodeList.push_back((*locit)->m_n1);
433  nodeList.push_back((*locit)->m_n2);
434 
435  ElmtConfig conf(elType, 1, true, true);
437  elType, conf, nodeList, tags);
438  E->SetId((*locit)->m_id);
439  m_mesh->m_element[E->GetDim()].push_back(E);
440  newComp->m_items.push_back(E);
441  }
442  }
443 
444  if (newComp->m_items.size())
445  {
446  m_mesh->m_composite[maxId++] = newComp;
447  }
448  }
449  }
450 }
451 }
452 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
virtual void Process()
Write mesh to output file.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
std::map< unsigned int, CompositeSharedPtr > CompositeMap
Container of composites; key is the composite id, value is the composite.
Definition: Composite.h:124
STL namespace.
pair< ModuleType, string > ModuleKey
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
boost::shared_ptr< Composite > CompositeSharedPtr
Shared pointer to a composite.
Definition: Composite.h:121
Represents a command-line configuration option.
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:135
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Abstract base class for processing modules.
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:161
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:148
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215