Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ProcessDetectSurf.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessDetectSurf.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 "ProcessDetectSurf.h"
38 
41 
42 #include <vector>
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace Utilities
48  {
49  ModuleKey ProcessDetectSurf::className =
51  ModuleKey(eProcessModule, "detect"), ProcessDetectSurf::create,
52  "Process elements to detect a surface.");
53 
54  ProcessDetectSurf::ProcessDetectSurf(MeshSharedPtr m) : ProcessModule(m)
55  {
56  m_config["vol"] = ConfigOption(false, "-1",
57  "Tag identifying surface to process.");
58  }
59 
61  {
62  }
63 
64  struct EdgeInfo {
65  EdgeInfo() : count(0) {}
66  int count;
68  unsigned int group;
69  };
70 
72  {
73  if (m_mesh->m_expDim > 2)
74  {
75  cerr << "Surface detection only implemented for 2D meshes" << endl;
76  return;
77  }
78 
79  int i, j;
80  string surf = m_config["vol"].as<string>();
81 
82  // Obtain vector of surface IDs from string.
83  vector<unsigned int> surfs;
84  if (surf != "-1")
85  {
86  ParseUtils::GenerateSeqVector(surf.c_str(), surfs);
87  sort(surfs.begin(), surfs.end());
88  }
89 
90  // If we're running in verbose mode print out a list of surfaces.
91  if (m_mesh->m_verbose)
92  {
93  cout << "ProcessDetectSurf: detecting surfaces";
94  if (surfs.size() > 0)
95  {
96  cout << " for surface" << (surfs.size() == 1 ? "" : "s")
97  << " " << surf << endl;
98  }
99  }
100 
101  vector<ElementSharedPtr> &el = m_mesh->m_element[m_mesh->m_expDim];
102  map<int, EdgeInfo> edgeCount;
103  set<int> doneIds;
104  map<int, int> idMap;
105 
106  // Iterate over list of surface elements.
107  for (i = 0; i < el.size(); ++i)
108  {
109  // Work out whether this lies on our surface of interest.
110  if (surfs.size() > 0)
111  {
112  vector<int> inter, tags = el[i]->GetTagList();
113 
114  sort(tags.begin(), tags.end());
115  set_intersection(surfs.begin(), surfs.end(),
116  tags .begin(), tags .end(),
117  back_inserter(inter));
118 
119  // It doesn't continue to next element.
120  if (inter.size() != 1)
121  {
122  continue;
123  }
124  }
125 
126  // List all edges.
127  ElementSharedPtr elmt = el[i];
128  for (j = 0; j < elmt->GetEdgeCount(); ++j)
129  {
130  EdgeSharedPtr e = elmt->GetEdge(j);
131  int eId = e->m_id;
132  edgeCount[eId].count++;
133  edgeCount[eId].edge = e;
134  }
135 
136  doneIds.insert(elmt->GetId());
137  ASSERTL0(idMap.count(elmt->GetId()) == 0, "Shouldn't happen");
138  idMap[elmt->GetId()] = i;
139  }
140 
142  unsigned int maxId = 0;
143 
144  for (cIt = m_mesh->m_composite.begin(); cIt != m_mesh->m_composite.end(); ++cIt)
145  {
146  maxId = (std::max)(cIt->first, maxId);
147  }
148 
149  ++maxId;
150 
152 
153  while (doneIds.size() > 0)
154  {
155  ElementSharedPtr start
156  = m_mesh->m_element[m_mesh->m_expDim][idMap[*(doneIds.begin())]];
157 
158  vector<ElementSharedPtr> block;
159  FindContiguousSurface(start, doneIds, block);
160  ASSERTL0(block.size() > 0, "Contiguous block not found");
161 
162  // Loop over all edges in block.
163  for (i = 0; i < block.size(); ++i)
164  {
165  // Find edge info.
166  ElementSharedPtr elmt = block[i];
167 
168  for (j = 0; j < elmt->GetEdgeCount(); ++j)
169  {
170  eIt = edgeCount.find(elmt->GetEdge(j)->m_id);
171  ASSERTL0(eIt != edgeCount.end(), "Couldn't find edge");
172  eIt->second.group = maxId;
173  }
174  }
175 
176  ++maxId;
177  }
178 
179  for (eIt = edgeCount.begin(); eIt != edgeCount.end(); ++eIt)
180  {
181  if (eIt->second.count > 1)
182  {
183  continue;
184  }
185 
186  unsigned int compId = eIt->second.group;
187  CompositeMap::iterator cIt = m_mesh->m_composite.find(compId);
188 
189  if (cIt == m_mesh->m_composite.end())
190  {
191  CompositeSharedPtr comp(new Composite());
192  comp->m_id = compId;
193  comp->m_tag = "E";
194  cIt = m_mesh->m_composite.insert(std::make_pair(compId, comp)).first;
195  }
196 
197  vector<int> tags(1);
198  tags[0] = compId;
199  vector<NodeSharedPtr> nodeList(2);
200  nodeList[0] = eIt->second.edge->m_n1;
201  nodeList[1] = eIt->second.edge->m_n2;
202 
203  ElmtConfig conf(LibUtilities::eSegment, 1, false, false);
205  CreateInstance(LibUtilities::eSegment,conf,nodeList,tags);
206  elmt->SetEdgeLink(eIt->second.edge);
207 
208  cIt->second->m_items.push_back(elmt);
209  }
210  }
211 
213  ElementSharedPtr start,
214  set<int> &doneIds,
215  vector<ElementSharedPtr> &block)
216  {
217  block.push_back(start);
218  doneIds.erase(start->GetId());
219 
220  vector<EdgeSharedPtr> edges = start->GetEdgeList();
221 
222  for (int i = 0; i < edges.size(); ++i)
223  {
224  for (int j = 0; j < edges[i]->m_elLink.size(); ++j)
225  {
226  ElementSharedPtr elmt = edges[i]->m_elLink[j].first;
227  if (elmt == start)
228  {
229  continue;
230  }
231 
232  if (doneIds.count(elmt->GetId()) == 0)
233  {
234  continue;
235  }
236 
237  FindContiguousSurface(elmt, doneIds, block);
238  }
239  }
240  }
241  }
242 }