Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessLinear.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessLinear.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: linearises mesh.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 #include "ProcessLinear.h"
38 
39 using namespace std;
40 
41 namespace Nektar
42 {
43 namespace Utilities
44 {
45 
46 using namespace Nektar::NekMeshUtils;
47 
48 ModuleKey ProcessLinear::className = GetModuleFactory().RegisterCreatorFunction(
49  ModuleKey(eProcessModule, "linearise"),
50  ProcessLinear::create,
51  "Linearises mesh.");
52 
53 ProcessLinear::ProcessLinear(MeshSharedPtr m) : ProcessModule(m)
54 {
55  m_config["all"] =
56  ConfigOption(true, "0", "remove curve nodes for all elements.");
57  m_config["invalid"] =
58  ConfigOption(false, "0", "remove curve nodes if element is invalid.");
59  m_config["prismonly"] =
60  ConfigOption(true, "0", "only acts on prims");
61  m_config["extract"] =
62  ConfigOption(false, "", "dump a mesh of the extracted elements");
63 }
64 
66 {
67 }
68 
70 {
71  if (m_mesh->m_verbose)
72  {
73  cout << "ProcessLinear: Linearising mesh... " << endl;
74  }
75 
76  bool all = m_config["all"].as<bool>();
77  bool invalid = m_config["invalid"].beenSet;
78  NekDouble thr = m_config["invalid"].as<NekDouble>();
79 
80  ASSERTL0(all || invalid,
81  "Must specify an option: all (to remove all curvature) or invalid "
82  "(to remove curvature that makes elements invalid)");
83 
84  if (all)
85  {
87  for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end();
88  eit++)
89  {
90  (*eit)->m_edgeNodes.clear();
91  }
92 
94  for (fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end();
95  fit++)
96  {
97  (*fit)->m_faceNodes.clear();
98  }
99 
100  for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
101  {
102  vector<NodeSharedPtr> empty;
103  m_mesh->m_element[m_mesh->m_expDim][i]->SetVolumeNodes(empty);
104  }
105 
106  if (m_mesh->m_verbose)
107  {
108  cerr << "Removed all element curvature" << endl;
109  }
110  }
111  else if (invalid)
112  {
113  map<int,vector<FaceSharedPtr> > eidToFace;
114  map<int,vector<ElementSharedPtr> > eidToElm;
115 
116  vector<ElementSharedPtr> els = m_mesh->m_element[m_mesh->m_expDim];
117  vector<ElementSharedPtr> el = els;
118 
119  for(int i = 0; i < el.size(); i++)
120  {
121  vector<EdgeSharedPtr> e = el[i]->GetEdgeList();
122  for(int j = 0; j < e.size(); j++)
123  {
124  eidToElm[e[j]->m_id].push_back(el[i]);
125  }
126  }
127 
128  if(m_mesh->m_expDim > 2)
129  {
131  for(it = m_mesh->m_faceSet.begin();
132  it != m_mesh->m_faceSet.end(); it++)
133  {
134  vector<EdgeSharedPtr> es = (*it)->m_edgeList;
135  for(int i = 0; i < es.size(); i++)
136  {
137  eidToFace[es[i]->m_id].push_back((*it));
138  }
139  }
140  }
141 
142  set<int> neigh;
143  vector<NodeSharedPtr> zeroNodes;
144  boost::unordered_set<int> clearedEdges, clearedFaces, clearedElmts;
145 
146  vector<ElementSharedPtr> dumpEls;
147 
148  // Iterate over list of elements of expansion dimension.
149  while(el.size() > 0)
150  {
151  for (int i = 0; i < el.size(); ++i)
152  {
153  if(m_config["prismonly"].beenSet)
154  {
155  if(el[i]->GetConf().m_e != LibUtilities::ePrism)
156  {
157  continue;
158  }
159  }
160 
161  if (Invalid(el[i],thr))//(!gfac->IsValid())
162  {
163  dumpEls.push_back(el[i]);
164  clearedElmts.insert(el[i]->GetId());;
165  el[i]->SetVolumeNodes(zeroNodes);
166 
167  vector<FaceSharedPtr> f = el[i]->GetFaceList();
168  for (int j = 0; j < f.size(); j++)
169  {
170  f[j]->m_faceNodes.clear();
171  clearedFaces.insert(f[j]->m_id);
172  }
173  vector<EdgeSharedPtr> e = el[i]->GetEdgeList();
174  for(int j = 0; j < e.size(); j++)
175  {
176  e[j]->m_edgeNodes.clear();
177  clearedEdges.insert(e[j]->m_id);
178  }
179 
180  if(m_mesh->m_expDim > 2)
181  {
182  for(int j = 0; j < e.size(); j++)
183  {
184  map<int,vector<FaceSharedPtr> >::iterator it =
185  eidToFace.find(e[j]->m_id);
186  for(int k = 0; k < it->second.size(); k++)
187  {
188  clearedEdges.insert(it->second[k]->m_id);
189  it->second[k]->m_faceNodes.clear();
190  }
191  }
192  }
193 
194  for(int j = 0; j < e.size(); j++)
195  {
196  map<int,vector<ElementSharedPtr> >::iterator it =
197  eidToElm.find(e[j]->m_id);
198  for(int k = 0; k < it->second.size(); k++)
199  {
200  neigh.insert(it->second[k]->GetId());
201  }
202  }
203  }
204  }
205 
206  el.clear();
207  set<int>::iterator it1;
209  for(int i = 0; i < els.size(); i++)
210  {
211  it1 = neigh.find(els[i]->GetId());
212  it2 = clearedElmts.find(els[i]->GetId());
213  if(it1 != neigh.end() && it2 == clearedElmts.end())
214  {
215  el.push_back(els[i]);
216  }
217  }
218  neigh.clear();
219  }
220 
221  if (m_mesh->m_verbose)
222  {
223  cerr << "Removed curvature from " << clearedElmts.size()
224  << " elements (" << clearedEdges.size() << " edges, "
225  << clearedFaces.size() << " faces)" << endl;
226  }
227 
228  if(m_config["extract"].beenSet)
229  {
230  MeshSharedPtr dmp = boost::shared_ptr<Mesh>(new Mesh());
231  dmp->m_expDim = 3;
232  dmp->m_spaceDim = 3;
233  dmp->m_nummode = 2;
234 
235  dmp->m_element[3] = dumpEls;
236 
238  ModuleKey(eOutputModule, "xml"), dmp);
239  mod->RegisterConfig("outfile", m_config["extract"].as<string>().c_str());
240  mod->ProcessVertices();
241  mod->ProcessEdges();
242  mod->ProcessFaces();
243  mod->ProcessElements();
244  mod->ProcessComposites();
245  mod->Process();
246  }
247  }
248 }
249 
251 {
252  // Create elemental geometry.
254  el->GetGeom(m_mesh->m_spaceDim);
255 
256  // Generate geometric factors.
258  geom->GetGeomFactors();
259 
260  if(!gfac->IsValid())
261  {
262  return true;
263  }
264 
265  vector<NodeSharedPtr> ns = el->GetVertexList();
266  ElmtConfig c = el->GetConf();
267  c.m_order = 1;
268  c.m_faceNodes = false;
269  c.m_volumeNodes = false;
270  c.m_reorient = false;
271 
273  c.m_e, c, ns, el->GetTagList());
274  SpatialDomains::GeometrySharedPtr geomL = elL->GetGeom(m_mesh->m_spaceDim);
275  SpatialDomains::GeomFactorsSharedPtr gfacL = geomL->GetGeomFactors();
276 
277  LibUtilities::PointsKeyVector p = geom->GetPointsKeys();
278  SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
279  SpatialDomains::DerivStorage derivL = gfacL->GetDeriv(p);
280  const int pts = deriv[0][0].num_elements();
281  Array<OneD,NekDouble> jc(pts);
282  Array<OneD,NekDouble> jcL(pts);
283  for (int k = 0; k < pts; ++k)
284  {
285  DNekMat jac(m_mesh->m_expDim, m_mesh->m_expDim, 0.0, eFULL);
286  DNekMat jacL(m_mesh->m_expDim, m_mesh->m_expDim, 0.0, eFULL);
287 
288  for (int l = 0; l < m_mesh->m_expDim; ++l)
289  {
290  for (int j = 0; j < m_mesh->m_expDim; ++j)
291  {
292  jac(j,l) = deriv[l][j][k];
293  jacL(j,l) = derivL[l][j][k];
294  }
295  }
296 
297  if(m_mesh->m_expDim == 2)
298  {
299  jc[k] = jac(0,0) * jac(1,1) - jac(0,1)*jac(1,0);
300  jcL[k] = jacL(0,0) * jacL(1,1) - jacL(0,1)*jacL(1,0);
301  }
302  else if(m_mesh->m_expDim == 3)
303  {
304  jc[k] = jac(0,0) * (jac(1,1)*jac(2,2) - jac(2,1)*jac(1,2)) -
305  jac(0,1) * (jac(1,0)*jac(2,2) - jac(2,0)*jac(1,2)) +
306  jac(0,2) * (jac(1,0)*jac(2,1) - jac(2,0)*jac(1,1));
307  jcL[k] = jacL(0,0) * (jacL(1,1)*jacL(2,2) - jacL(2,1)*jacL(1,2)) -
308  jacL(0,1) * (jacL(1,0)*jacL(2,2) - jacL(2,0)*jacL(1,2)) +
309  jacL(0,2) * (jacL(1,0)*jacL(2,1) - jacL(2,0)*jacL(1,1));
310  }
311  }
312 
313  Array<OneD, NekDouble> j(pts);
314  Vmath::Vdiv(jc.num_elements(),jc,1,jcL,1,j,1);
315 
316  NekDouble scaledJac = Vmath::Vmin(j.num_elements(),j,1) /
317  Vmath::Vmax(j.num_elements(),j,1);
318 
319  //cout << scaledJac << endl;
320 
321  if(scaledJac < thr)
322  {
323  return true;
324  }
325 
326  return false;
327 }
328 }
329 }
bool m_faceNodes
Denotes whether the element contains face nodes. For 2D elements, if this is true then the element co...
Definition: ElementConfig.h:80
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Basic information about an element.
Definition: ElementConfig.h:50
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
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
virtual void Process()
Write mesh to output file.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
STL namespace.
pair< ModuleType, string > ModuleKey
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:241
ElementFactory & GetElementFactory()
Definition: Element.cpp:47
unsigned int m_order
Order of the element.
Definition: ElementConfig.h:87
bool Invalid(NekMeshUtils::ElementSharedPtr el, NekDouble thr)
bool m_volumeNodes
Denotes whether the element contains volume (i.e. interior) nodes. These are not supported by either ...
Definition: ElementConfig.h:85
boost::shared_ptr< Module > ModuleSharedPtr
Represents a command-line configuration option.
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
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< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
bool m_reorient
Denotes whether the element needs to be re-orientated for a spectral element framework.
Definition: ElementConfig.h:90
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
std::pair< ModuleType, std::string > ModuleKey
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
LibUtilities::ShapeType m_e
Element type (e.g. triangle, quad, etc).
Definition: ElementConfig.h:77
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215