Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HOSurfaceMesh.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: SurfaceMeshing.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: surfacemeshing object methods.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <list>
37 #include <algorithm>
38 
44 
46 #include <LocalRegions/MatrixKey.h>
48 
49 using namespace std;
50 namespace Nektar
51 {
52 namespace NekMeshUtils
53 {
54 
55 ModuleKey HOSurfaceMesh::className = GetModuleFactory().RegisterCreatorFunction(
56  ModuleKey(eProcessModule, "hosurface"),
57  HOSurfaceMesh::create,
58  "Generates a high-order surface mesh based on CAD");
59 
60 HOSurfaceMesh::HOSurfaceMesh(MeshSharedPtr m) : ProcessModule(m)
61 {
62  m_config["opti"] =
63  ConfigOption(true, "0", "Perform edge node optimisation.");
64 }
65 
67 {
68 }
69 
71 {
72  if (m_mesh->m_verbose)
73  cout << endl << "\tHigh-Order Surface meshing" << endl;
74 
75  LibUtilities::PointsKey ekey(m_mesh->m_nummode,
78 
79  LibUtilities::PointsManager()[ekey]->GetPoints(gll);
80 
81  LibUtilities::PointsKey pkey(m_mesh->m_nummode,
83 
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;
283  for (int k = 0; k < nq; k++)
284  {
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  {
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  {
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  {
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 
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  {
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  {
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 
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  {
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 }
506 }
507 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
boost::shared_ptr< OptiEdge > OptiEdgeSharedPtr
base class for CAD curves.
Definition: CADCurve.h:52
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.
int GetId()
Return ID of the vertex.
Definition: CADObject.h:87
Represents a face comprised of three or more edges.
Definition: Face.h:61
STL namespace.
pair< ModuleType, string > ModuleKey
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
Represents a command-line configuration option.
base class for a cad surface
Definition: CADSurf.h:55
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Definition: Points.h:58
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:135
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
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.
bool BGFSUpdate(OptiObjSharedPtr opti, DNekMat &J, DNekMat &B, DNekMat &H)
Definition: BGFS-B.cpp:50
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:161
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:148
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
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215