Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
Nektar::NekMeshUtils::HOSurfaceMesh Class Reference

#include <HOSurfaceMesh.h>

Inheritance diagram for Nektar::NekMeshUtils::HOSurfaceMesh:
Inheritance graph
[legend]
Collaboration diagram for Nektar::NekMeshUtils::HOSurfaceMesh:
Collaboration graph
[legend]

Public Member Functions

 HOSurfaceMesh (MeshSharedPtr m)
 
virtual ~HOSurfaceMesh ()
 
virtual void Process ()
 
- Public Member Functions inherited from Nektar::NekMeshUtils::ProcessModule
NEKMESHUTILS_EXPORT ProcessModule (MeshSharedPtr p_m)
 
- Public Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT Module (MeshSharedPtr p_m)
 
NEKMESHUTILS_EXPORT void RegisterConfig (std::string key, std::string value)
 Register a configuration option with a module. More...
 
NEKMESHUTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
NEKMESHUTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
NEKMESHUTILS_EXPORT MeshSharedPtr GetMesh ()
 
virtual NEKMESHUTILS_EXPORT void ProcessVertices ()
 Extract element vertices. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessElements ()
 Generate element IDs. More...
 
virtual NEKMESHUTILS_EXPORT void ProcessComposites ()
 Generate composites. More...
 
virtual NEKMESHUTILS_EXPORT void ClearElementLinks ()
 

Static Public Member Functions

static boost::shared_ptr< Modulecreate (MeshSharedPtr m)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 

Additional Inherited Members

- Protected Member Functions inherited from Nektar::NekMeshUtils::Module
NEKMESHUTILS_EXPORT void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly. More...
 
NEKMESHUTILS_EXPORT void PrismLines (int prism, PerMap &perFaces, std::set< int > &prismsDone, std::vector< ElementSharedPtr > &line)
 
- Protected Attributes inherited from Nektar::NekMeshUtils::Module
MeshSharedPtr m_mesh
 Mesh object. More...
 
std::map< std::string,
ConfigOption
m_config
 List of configuration values. More...
 

Detailed Description

Definition at line 46 of file HOSurfaceMesh.h.

Constructor & Destructor Documentation

Nektar::NekMeshUtils::HOSurfaceMesh::HOSurfaceMesh ( MeshSharedPtr  m)

Definition at line 60 of file HOSurfaceMesh.cpp.

References Nektar::NekMeshUtils::Module::m_config.

60  : ProcessModule(m)
61 {
62  m_config["opti"] =
63  ConfigOption(true, "0", "Perform edge node optimisation.");
64 }
NEKMESHUTILS_EXPORT ProcessModule(MeshSharedPtr p_m)
std::map< std::string, ConfigOption > m_config
List of configuration values.
Nektar::NekMeshUtils::HOSurfaceMesh::~HOSurfaceMesh ( )
virtual

Definition at line 66 of file HOSurfaceMesh.cpp.

67 {
68 }

Member Function Documentation

static boost::shared_ptr<Module> Nektar::NekMeshUtils::HOSurfaceMesh::create ( MeshSharedPtr  m)
inlinestatic

Creates an instance of this class.

Definition at line 50 of file HOSurfaceMesh.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

51  {
53  }
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Nektar::NekMeshUtils::HOSurfaceMesh::Process ( )
virtual

Implements Nektar::NekMeshUtils::Module.

Definition at line 70 of file HOSurfaceMesh.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::NekMeshUtils::BGFSUpdate(), Nektar::NekMeshUtils::CADType::eCurve, Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eNodalTriElec, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::NekMeshUtils::CADObject::GetId(), Nektar::iterator, Nektar::NekMeshUtils::Module::m_config, Nektar::NekMeshUtils::Module::m_mesh, class_topology::Node, Nektar::LibUtilities::PointsManager(), and Nektar::LibUtilities::PrintProgressbar().

71 {
72  if (m_mesh->m_verbose)
73  cout << endl << "\tHigh-Order Surface meshing" << endl;
74 
75  LibUtilities::PointsKey ekey(m_mesh->m_nummode,
77  Array<OneD, NekDouble> gll;
78 
79  LibUtilities::PointsManager()[ekey]->GetPoints(gll);
80 
81  LibUtilities::PointsKey pkey(m_mesh->m_nummode,
83 
84  Array<OneD, NekDouble> u, v;
85 
86  int nq = m_mesh->m_nummode;
87 
88  int np = nq * (nq + 1) / 2;
89 
90  int ni = (nq-2)*(nq-3) / 2;
91 
92  int npq = nq * nq;
93 
94  int niq = npq - 4 - 4*(nq-2);
95 
96  LibUtilities::PointsManager()[pkey]->GetPoints(u, v);
97 
98  bool qOpti = m_config["opti"].beenSet;
99 
100  // loop over all the faces in the surface mesh, check all three edges for
101  // high order info, if nothing high-order the edge.
102 
103  EdgeSet surfaceEdges;
104  EdgeSet completedEdges;
105 
106  for(int i = 0; i < m_mesh->m_element[2].size(); i++)
107  {
108  vector<EdgeSharedPtr> es = m_mesh->m_element[2][i]->GetEdgeList();
109  for(int j = 0; j < es.size(); j++)
110  surfaceEdges.insert(es[j]);
111  }
112 
113  for (int i = 0; i < m_mesh->m_element[2].size(); i++)
114  {
115  if (m_mesh->m_verbose)
116  {
118  i, m_mesh->m_element[2].size(), "\t\tSurface elements");
119  }
120 
121  CADObjectSharedPtr o = m_mesh->m_element[2][i]->m_parentCAD;
122  CADSurfSharedPtr s = boost::dynamic_pointer_cast<CADSurf>(o);
123  int surf = s->GetId();
124 
125  FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
126 
127  if(!f)
128  {
129  f = boost::shared_ptr<Face>(new Face(m_mesh->m_element[2][i]->GetVertexList(),
130  vector<NodeSharedPtr>(),
131  m_mesh->m_element[2][i]->GetEdgeList(),
133  }
134 
135  f->m_parentCAD = s;
136 
137  vector<EdgeSharedPtr> edges = f->m_edgeList;
138  for (int j = 0; j < edges.size(); j++)
139  {
140  EdgeSharedPtr e = edges[j];
141  // test insert the edge into completedEdges
142  // if the edge already exists move on
143  // if not figure out its high-order information
144 
145  EdgeSet::iterator test = completedEdges.find(e);
146 
147  if (test != completedEdges.end())
148  {
149  continue;
150  }
151 
152  // the edges in the element are different to those in the face
153  // the cad information is stored in the element edges which are not
154  // in the m_mesh->m_edgeSet groups
155  // need to link them together and copy the cad information to be
156  // able to identify how to make it high-order
157  EdgeSet::iterator it = surfaceEdges.find(e);
158  ASSERTL0(it != surfaceEdges.end(),"could not find edge in surface");
159 
160  if((*it)->m_parentCAD)
161  {
162  e->m_parentCAD = (*it)->m_parentCAD;
163  }
164  else
165  {
166  e->m_parentCAD = s;
167  }
168 
169  vector<NodeSharedPtr> honodes(m_mesh->m_nummode - 2);
170 
171  if (e->m_parentCAD->GetType() == CADType::eCurve)
172  {
173  int cid = e->m_parentCAD->GetId();
174  CADCurveSharedPtr c = boost::dynamic_pointer_cast<CADCurve>(e->m_parentCAD);
175  NekDouble tb = e->m_n1->GetCADCurveInfo(cid);
176  NekDouble te = e->m_n2->GetCADCurveInfo(cid);
177 
178  // distrobute points along curve as inital guess
179  Array<OneD, NekDouble> ti(m_mesh->m_nummode);
180  for (int k = 0; k < m_mesh->m_nummode; k++)
181  {
182  ti[k] =
183  tb * (1.0 - gll[k]) / 2.0 + te * (1.0 + gll[k]) / 2.0;
184  }
185 
186  if(qOpti)
187  {
188  Array<OneD, NekDouble> xi(nq - 2);
189  for (int k = 1; k < nq - 1; k++)
190  {
191  xi[k - 1] = ti[k];
192  }
193 
194  OptiEdgeSharedPtr opti =
196 
197  DNekMat B(
198  nq - 2, nq - 2, 0.0); // approximate hessian (I to start)
199  for (int k = 0; k < nq - 2; k++)
200  {
201  B(k, k) = 1.0;
202  }
203  DNekMat H(nq - 2,
204  nq - 2,
205  0.0); // approximate inverse hessian (I to start)
206  for (int k = 0; k < nq - 2; k++)
207  {
208  H(k, k) = 1.0;
209  }
210 
211  DNekMat J = opti->dF(xi);
212 
213  Array<OneD, NekDouble> bnds = c->GetBounds();
214 
215  bool repeat = true;
216  int itct = 0;
217  while (repeat)
218  {
219  NekDouble Norm = 0;
220  for (int k = 0; k < nq - 2; k++)
221  {
222  Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
223  (bnds[1] - bnds[0]);
224  }
225  Norm = sqrt(Norm);
226 
227  if (Norm < 1E-8)
228  {
229  repeat = false;
230  break;
231  }
232  if (itct > 1000)
233  {
234  cout << "failed to optimise on curve" << endl;
235  for (int k = 0; k < nq; k++)
236  {
237  ti[k] = tb * (1.0 - gll[k]) / 2.0 +
238  te * (1.0 + gll[k]) / 2.0;
239  }
240  break;
241  }
242  itct++;
243 
244  if (!BGFSUpdate(opti, J, B, H))
245  {
246  if(m_mesh->m_verbose)
247  {
248  cout << "BFGS reported no update, curve on "
249  << c->GetId() << endl;
250  }
251  break;
252  }
253  }
254  // need to pull the solution out of opti
255  ti = opti->GetSolution();
256  }
257  vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s = c->GetAdjSurf();
258 
259  for (int k = 1; k < m_mesh->m_nummode - 1; k++)
260  {
261  Array<OneD, NekDouble> loc = c->P(ti[k]);
262  NodeSharedPtr nn = boost::shared_ptr<Node>(
263  new Node(0, loc[0], loc[1], loc[2]));
264 
265  nn->SetCADCurve(cid, c, ti[k]);
266  for(int m = 0; m < s.size(); m++)
267  {
268  Array<OneD, NekDouble> uv = s[m].first->locuv(loc);
269  nn->SetCADSurf(s[m].first->GetId(), s[m].first, uv);
270  }
271 
272  honodes[k - 1] = nn;
273  }
274  }
275  else
276  {
277  // edge is on surface and needs 2d optimisation
278  Array<OneD, NekDouble> uvb, uve;
279  uvb = e->m_n1->GetCADSurfInfo(surf);
280  uve = e->m_n2->GetCADSurfInfo(surf);
281  e->m_parentCAD = s;
282  Array<OneD, Array<OneD, NekDouble> > uvi(nq);
283  for (int k = 0; k < nq; k++)
284  {
285  Array<OneD, NekDouble> uv(2);
286  uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
287  uve[0] * (1.0 + gll[k]) / 2.0;
288  uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
289  uve[1] * (1.0 + gll[k]) / 2.0;
290  uvi[k] = uv;
291  }
292 
293  if(qOpti)
294  {
295  Array<OneD, NekDouble> bnds = s->GetBounds();
296  Array<OneD, NekDouble> all(2 * nq);
297  for (int k = 0; k < nq; k++)
298  {
299  all[k * 2 + 0] = uvi[k][0];
300  all[k * 2 + 1] = uvi[k][1];
301  }
302 
303  Array<OneD, NekDouble> xi(2 * (nq - 2));
304  for (int k = 1; k < nq - 1; k++)
305  {
306  xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
307  xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
308  }
309 
310  OptiEdgeSharedPtr opti =
312 
313  DNekMat B(2 * (nq - 2),
314  2 * (nq - 2),
315  0.0); // approximate hessian (I to start)
316  for (int k = 0; k < 2 * (nq - 2); k++)
317  {
318  B(k, k) = 1.0;
319  }
320  DNekMat H(2 * (nq - 2),
321  2 * (nq - 2),
322  0.0); // approximate inverse hessian (I to start)
323  for (int k = 0; k < 2 * (nq - 2); k++)
324  {
325  H(k, k) = 1.0;
326  }
327 
328  DNekMat J = opti->dF(xi);
329 
330  bool repeat = true;
331  int itct = 0;
332  while (repeat)
333  {
334  NekDouble Norm = 0;
335  for (int k = 0; k < 2 * (nq - 2); k++)
336  {
337  if (k % 2 == 0)
338  {
339  Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
340  (bnds[1] - bnds[0]);
341  }
342  else
343  {
344  Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
345  (bnds[3] - bnds[2]);
346  }
347  }
348  Norm = sqrt(Norm);
349 
350  if (Norm < 1E-8)
351  {
352  repeat = false;
353  break;
354  }
355 
356  if (itct > 1000)
357  {
358  cout << "failed to optimise on edge" << endl;
359  for (int k = 0; k < nq; k++)
360  {
361  Array<OneD, NekDouble> uv(2);
362  uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
363  uve[0] * (1.0 + gll[k]) / 2.0;
364  uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
365  uve[1] * (1.0 + gll[k]) / 2.0;
366  uvi[k] = uv;
367  }
368  break;
369  }
370  itct++;
371 
372  if (!BGFSUpdate(opti, J, B, H))
373  {
374  if(m_mesh->m_verbose)
375  {
376  cout << "BFGS reported no update, edge on " << surf
377  << endl;
378  }
379  // exit(-1);
380  break;
381  }
382  }
383 
384  all = opti->GetSolution();
385 
386  // need to put all backinto uv
387  for (int k = 0; k < nq; k++)
388  {
389  uvi[k][0] = all[k * 2 + 0];
390  uvi[k][1] = all[k * 2 + 1];
391  }
392  }
393 
394  for (int k = 1; k < nq - 1; k++)
395  {
396  Array<OneD, NekDouble> loc;
397  loc = s->P(uvi[k]);
398  NodeSharedPtr nn = boost::shared_ptr<Node>(
399  new Node(0, loc[0], loc[1], loc[2]));
400  nn->SetCADSurf(s->GetId(), s, uvi[k]);
401  honodes[k - 1] = nn;
402  }
403  }
404 
405  e->m_edgeNodes = honodes;
406  e->m_curveType = LibUtilities::eGaussLobattoLegendre;
407  completedEdges.insert(e);
408  }
409 
410  //just add the face interior nodes through interp and project
411  vector<NodeSharedPtr> vertices = f->m_vertexList;
412 
413  SpatialDomains::GeometrySharedPtr geom = f->GetGeom(3);
414  geom->FillGeom();
415  StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
416  Array<OneD, NekDouble> coeffs0 = geom->GetCoeffs(0);
417  Array<OneD, NekDouble> coeffs1 = geom->GetCoeffs(1);
418  Array<OneD, NekDouble> coeffs2 = geom->GetCoeffs(2);
419 
420  Array<OneD, NekDouble> xc(xmap->GetTotPoints());
421  Array<OneD, NekDouble> yc(xmap->GetTotPoints());
422  Array<OneD, NekDouble> zc(xmap->GetTotPoints());
423 
424  xmap->BwdTrans(coeffs0,xc);
425  xmap->BwdTrans(coeffs1,yc);
426  xmap->BwdTrans(coeffs2,zc);
427 
428  if(vertices.size() == 3)
429  {
430  //build an array of all uvs
431  vector<Array<OneD, NekDouble> > uvi;
432  for(int j = np-ni; j < np; j++)
433  {
434  Array<OneD, NekDouble> xp(2);
435  xp[0] = u[j];
436  xp[1] = v[j];
437 
438  Array<OneD, NekDouble> loc(3);
439  loc[0] = xmap->PhysEvaluate(xp, xc);
440  loc[1] = xmap->PhysEvaluate(xp, yc);
441  loc[2] = xmap->PhysEvaluate(xp, zc);
442 
443  Array<OneD, NekDouble> uv(2);
444  s->ProjectTo(loc,uv);
445  uvi.push_back(uv);
446 
447  }
448 
449  vector<NodeSharedPtr> honodes;
450  for(int j = 0; j < ni; j++)
451  {
452  Array<OneD, NekDouble> loc;
453  loc = s->P(uvi[j]);
454  NodeSharedPtr nn = boost::shared_ptr<Node>(new
455  Node(0,loc[0],loc[1],loc[2]));
456  nn->SetCADSurf(surf, s, uvi[j]);
457  honodes.push_back(nn);
458  }
459 
460  f->m_faceNodes = honodes;
461  f->m_curveType = LibUtilities::eNodalTriElec;
462  }
463  else if(vertices.size() == 4)
464  {
465  //build an array of all uvs
466  vector<Array<OneD, NekDouble> > uvi;
467  for(int i = 1; i < nq - 1; i++)
468  {
469  for(int j = 1; j < nq - 1; j++)
470  {
471  Array<OneD, NekDouble> xp(2);
472  xp[0] = gll[j];
473  xp[1] = gll[i];
474 
475  Array<OneD, NekDouble> loc(3);
476  loc[0] = xmap->PhysEvaluate(xp, xc);
477  loc[1] = xmap->PhysEvaluate(xp, yc);
478  loc[2] = xmap->PhysEvaluate(xp, zc);
479 
480  Array<OneD, NekDouble> uv(2);
481  s->ProjectTo(loc,uv);
482  uvi.push_back(uv);
483 
484  }
485  }
486 
487  vector<NodeSharedPtr> honodes;
488  for(int j = 0; j < niq; j++)
489  {
490  Array<OneD, NekDouble> loc;
491  loc = s->P(uvi[j]);
492  NodeSharedPtr nn = boost::shared_ptr<Node>(new
493  Node(0,loc[0],loc[1],loc[2]));
494  nn->SetCADSurf(surf, s, uvi[j]);
495  honodes.push_back(nn);
496  }
497 
498  f->m_faceNodes = honodes;
499  f->m_curveType = LibUtilities::eGaussLobattoLegendre;
500  }
501  }
502 
503  if (m_mesh->m_verbose)
504  cout << endl;
505 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
boost::shared_ptr< OptiEdge > OptiEdgeSharedPtr
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
Definition: Progressbar.hpp:69
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
Definition: NekTypeDefs.hpp:52
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< CADObject > CADObjectSharedPtr
Definition: CADObject.h:112
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:136
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool BGFSUpdate(OptiObjSharedPtr opti, DNekMat &J, DNekMat &B, DNekMat &H)
Definition: BGFS-B.cpp:50
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:162
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:153
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:70
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:194

Member Data Documentation

ModuleKey Nektar::NekMeshUtils::HOSurfaceMesh::className
static
Initial value:
ModuleKey(eProcessModule, "hosurface"),
"Generates a high-order surface mesh based on CAD")

Definition at line 54 of file HOSurfaceMesh.h.