Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
OptimiseFunctions.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OptimseFunctions.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: functions and gradients for optimisation
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
38 
39 using namespace std;
40 namespace Nektar
41 {
42 namespace NekMeshUtils
43 {
44 
46 {
47  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
48 }
50 {
52  ret[0] = a[0] - b[0];
53  ret[1] = a[1] - b[1];
54  ret[2] = a[2] - b[2];
55  return ret;
56 }
58 {
60  ret[0] = a[0] * t;
61  ret[1] = a[1] * t;
62  ret[2] = a[2] * t;
63  return ret;
64 }
66 {
68  ret[0] = a[0] + b[0];
69  ret[1] = a[1] + b[1];
70  ret[2] = a[2] + b[2];
71  return ret;
72 }
73 
74 Array<OneD, NekDouble> OptiEdge::Getxi()
75 {
77 
78  switch (o->GetType())
79  {
80  case curve:
81  xi = Array<OneD, NekDouble>(all.num_elements() - 2);
82  for (int i = 1; i < all.num_elements() - 1; i++)
83  {
84  xi[i - 1] = all[i];
85  }
86  break;
87 
88  case surf:
89  xi = Array<OneD, NekDouble>(all.num_elements() - 4);
90  for (int i = 2; i < all.num_elements() - 2; i++)
91  {
92  xi[i - 2] = all[i];
93  }
94  break;
95 
96  case vert:
97  ASSERTL0(false, "Should not be able to pass vert");
98  }
99  return xi;
100 }
101 
102 Array<OneD, NekDouble> OptiEdge::Getli()
103 {
106  switch (o->GetType())
107  {
108  case curve:
109  li = Array<OneD, NekDouble>(all.num_elements() - 2);
110  bnds = boost::dynamic_pointer_cast<CADCurve>(o)->Bounds();
111  for (int i = 1; i < all.num_elements() - 1; i++)
112  {
113  li[i - 1] = bnds[0];
114  }
115  break;
116 
117  case surf:
118  li = Array<OneD, NekDouble>(all.num_elements() - 4);
119  bnds = boost::dynamic_pointer_cast<CADSurf>(o)->GetBounds();
120  for (int i = 2; i < all.num_elements() - 2; i++)
121  {
122  if (i % 2 == 0)
123  {
124  li[i - 2] = bnds[0];
125  }
126  else
127  {
128  li[i - 2] = bnds[2];
129  }
130  }
131  break;
132 
133  case vert:
134  ASSERTL0(false, "Should not be able to pass vert");
135  }
136  return li;
137 }
138 
139 Array<OneD, NekDouble> OptiEdge::Getui()
140 {
143  switch (o->GetType())
144  {
145  case curve:
146  ui = Array<OneD, NekDouble>(all.num_elements() - 2);
147  bnds = boost::dynamic_pointer_cast<CADCurve>(o)->Bounds();
148  for (int i = 1; i < all.num_elements() - 1; i++)
149  {
150  ui[i - 1] = bnds[1];
151  }
152  break;
153 
154  case surf:
155  ui = Array<OneD, NekDouble>(all.num_elements() - 4);
156  bnds = boost::dynamic_pointer_cast<CADSurf>(o)->GetBounds();
157  for (int i = 2; i < all.num_elements() - 2; i++)
158  {
159  if (i % 2 == 0)
160  {
161  ui[i - 2] = bnds[1];
162  }
163  else
164  {
165  ui[i - 2] = bnds[3];
166  }
167  }
168  break;
169 
170  case vert:
171  ASSERTL0(false, "Should not be able to pass vert");
172  }
173  return ui;
174 }
175 
177 {
178  Array<OneD, NekDouble> val(all.num_elements());
179 
180  if (o->GetType() == curve)
181  {
182  val[0] = all[0];
183  for (int i = 0; i < xitst.num_elements(); i++)
184  {
185  val[i + 1] = xitst[i];
186  }
187  val[all.num_elements() - 1] = all[all.num_elements() - 1];
188  }
189  else if (o->GetType() == surf)
190  {
191  val[0] = all[0];
192  val[1] = all[1];
193  for (int i = 0; i < xitst.num_elements(); i++)
194  {
195  val[i + 2] = xitst[i];
196  }
197  val[all.num_elements() - 2] = all[all.num_elements() - 2];
198  val[all.num_elements() - 1] = all[all.num_elements() - 1];
199  }
200 
201  NekDouble ret = 0.0;
202  if (o->GetType() == curve)
203  {
204  CADCurveSharedPtr c = boost::dynamic_pointer_cast<CADCurve>(o);
205 
206  for (int i = 0; i < all.num_elements() - 1; i++)
207  {
208  Array<OneD, NekDouble> dis = Take(c->P(val[i + 1]), c->P(val[i]));
209  NekDouble norm =
210  dis[0] * dis[0] + dis[1] * dis[1] + dis[2] * dis[2];
211  ret += norm / (z[i + 1] - z[i]);
212  }
213  }
214  else if (o->GetType() == surf)
215  {
216  CADSurfSharedPtr s = boost::dynamic_pointer_cast<CADSurf>(o);
217  // need to organise the val array
218  Array<OneD, Array<OneD, NekDouble> > uv(val.num_elements() / 2);
219  for (int i = 0; i < val.num_elements() / 2; i++)
220  {
221  uv[i] = Array<OneD, NekDouble>(2);
222  uv[i][0] = val[i * 2 + 0];
223  uv[i][1] = val[i * 2 + 1];
224  }
225  for (int i = 0; i < uv.num_elements() - 1; i++)
226  {
227  Array<OneD, NekDouble> dis = Take(s->P(uv[i + 1]), s->P(uv[i]));
228  NekDouble norm =
229  dis[0] * dis[0] + dis[1] * dis[1] + dis[2] * dis[2];
230  ret += norm / (z[i + 1] - z[i]);
231  }
232  }
233  return ret;
234 }
235 
236 DNekMat OptiEdge::dF(Array<OneD, NekDouble> xitst)
237 {
238  Array<OneD, NekDouble> val(all.num_elements());
239 
240  if (o->GetType() == curve)
241  {
242  val[0] = all[0];
243  for (int i = 0; i < xitst.num_elements(); i++)
244  {
245  val[i + 1] = xitst[i];
246  }
247  val[all.num_elements() - 1] = all[all.num_elements() - 1];
248  }
249  else if (o->GetType() == surf)
250  {
251  val[0] = all[0];
252  val[1] = all[1];
253  for (int i = 0; i < xitst.num_elements(); i++)
254  {
255  val[i + 2] = xitst[i];
256  }
257  val[all.num_elements() - 2] = all[all.num_elements() - 2];
258  val[all.num_elements() - 1] = all[all.num_elements() - 1];
259  }
260 
261  DNekMat ret;
262 
263  if (o->GetType() == curve)
264  {
265  CADCurveSharedPtr c = boost::dynamic_pointer_cast<CADCurve>(o);
266  vector<Array<OneD, NekDouble> > r;
267  vector<Array<OneD, NekDouble> > dr;
268 
269  for (int i = 0; i < all.num_elements(); i++)
270  {
271  Array<OneD, NekDouble> ri(3), dri(3);
272  Array<OneD, NekDouble> d2 = c->D2(val[i]);
273  for (int j = 0; j < 3; j++)
274  {
275  ri[j] = d2[j];
276  dri[j] = d2[j + 3];
277  }
278  r.push_back(ri);
279  dr.push_back(dri);
280  }
281 
282  DNekMat J(all.num_elements() - 2, 1, 0.0);
283  for (int i = 0; i < all.num_elements() - 2; i++)
284  {
285  J(i, 0) =
286  2.0 / (z[i + 1] - z[i]) * Dot(dr[i + 1], Take(r[i + 1], r[i])) -
287  2.0 / (z[i + 2] - z[i + 1]) *
288  Dot(dr[i + 1], Take(r[i + 2], r[i + 1]));
289  }
290 
291  ret = J;
292  }
293  else if (o->GetType() == surf)
294  {
295  CADSurfSharedPtr s = boost::dynamic_pointer_cast<CADSurf>(o);
296  // need to organise the all array
297  Array<OneD, Array<OneD, NekDouble> > uv(val.num_elements() / 2);
298  for (int i = 0; i < val.num_elements() / 2; i++)
299  {
300  uv[i] = Array<OneD, NekDouble>(2);
301  uv[i][0] = val[i * 2 + 0];
302  uv[i][1] = val[i * 2 + 1];
303  }
304 
305  vector<Array<OneD, NekDouble> > r;
306  vector<Array<OneD, NekDouble> > dru, drv;
307  for (int i = 0; i < uv.num_elements(); i++)
308  {
309  Array<OneD, NekDouble> ri(3), drui(3), drvi(3);
310  Array<OneD, NekDouble> d2 = s->D1(uv[i]);
311  for (int j = 0; j < 3; j++)
312  {
313  ri[j] = d2[j];
314  drui[j] = d2[j + 3];
315  drvi[j] = d2[j + 6];
316  }
317  r.push_back(ri);
318  dru.push_back(drui);
319  drv.push_back(drvi);
320  }
321 
322  DNekMat J(2 * (uv.num_elements() - 2), 1, 0.0);
323  for (int i = 0; i < uv.num_elements() - 2; i++)
324  {
325  J(2 * i + 0, 0) = 2.0 / (z[i + 1] - z[i]) *
326  Dot(dru[i + 1], Take(r[i + 1], r[i])) +
327  2.0 / (z[i + 2] - z[i + 1]) *
328  Dot(dru[i + 1], Take(r[i + 1], r[i + 2]));
329 
330  J(2 * i + 1, 0) = 2.0 / (z[i + 1] - z[i]) *
331  Dot(drv[i + 1], Take(r[i + 1], r[i])) +
332  2.0 / (z[i + 2] - z[i + 1]) *
333  Dot(drv[i + 1], Take(r[i + 1], r[i + 2]));
334  }
335 
336  ret = J;
337  }
338 
339  return ret;
340 }
341 
342 void OptiEdge::Update(Array<OneD, NekDouble> xinew)
343 {
344  if (o->GetType() == curve)
345  {
346  for (int i = 0; i < xinew.num_elements(); i++)
347  {
348  all[i + 1] = xinew[i];
349  }
350  }
351  else if (o->GetType() == surf)
352  {
353  for (int i = 0; i < xinew.num_elements(); i++)
354  {
355  all[i + 2] = xinew[i];
356  }
357  }
358 }
359 
360 Array<OneD, NekDouble> OptiFace::Getxi()
361 {
362  Array<OneD, NekDouble> ret(ni * 2);
363  for (int i = np - ni; i < np; i++)
364  {
365  ret[2 * (i - np + ni) + 0] = uv[i][0];
366  ret[2 * (i - np + ni) + 1] = uv[i][1];
367  }
368  return ret;
369 }
370 
371 Array<OneD, NekDouble> OptiFace::Getli()
372 {
373  Array<OneD, NekDouble> li(ni * 2);
374  Array<OneD, NekDouble> bnds = s->GetBounds();
375  for (int i = np - ni; i < np; i++)
376  {
377  li[2 * (i - np + ni) + 0] = bnds[0];
378  li[2 * (i - np + ni) + 1] = bnds[2];
379  }
380  return li;
381 }
382 
383 Array<OneD, NekDouble> OptiFace::Getui()
384 {
385  Array<OneD, NekDouble> ui(ni * 2);
386  Array<OneD, NekDouble> bnds = s->GetBounds();
387  for (int i = np - ni; i < np; i++)
388  {
389  ui[2 * (i - np + ni) + 0] = bnds[1];
390  ui[2 * (i - np + ni) + 1] = bnds[3];
391  }
392  return ui;
393 }
394 
396 {
398  for (int i = np - ni; i < np; i++)
399  {
400  val[i][0] = xitst[(i - np + ni) * 2 + 0];
401  val[i][1] = xitst[(i - np + ni) * 2 + 1];
402  }
403 
404  NekDouble ret = 0.0;
405  set<pair<int, int> >::iterator it;
406  for (it = spring.begin(); it != spring.end(); it++)
407  {
409  Take(s->P(uv[(*it).first]), s->P(uv[(*it).second]));
410  NekDouble norm = dis[0] * dis[0] + dis[1] * dis[1] + dis[2] * dis[2];
411  ret += norm / z[(*it)];
412  }
413 
414  return ret;
415 }
416 
417 DNekMat OptiFace::dF(Array<OneD, NekDouble> xitst)
418 {
420  for (int i = np - ni; i < np; i++)
421  {
422  val[i][0] = xitst[(i - np + ni) * 2 + 0];
423  val[i][1] = xitst[(i - np + ni) * 2 + 1];
424  }
425 
426  DNekMat ret(ni * 2, 1, 0.0);
427 
428  vector<Array<OneD, NekDouble> > r, dru, drv;
429 
430  for (int i = 0; i < val.num_elements(); i++)
431  {
432  Array<OneD, NekDouble> ri(3), du(3), dv(3);
433  Array<OneD, NekDouble> d1 = s->D1(val[i]);
434  for (int i = 0; i < 3; i++)
435  {
436  ri[i] = d1[i];
437  du[i] = d1[i + 3];
438  dv[i] = d1[i + 6];
439  }
440  r.push_back(ri);
441  dru.push_back(du);
442  drv.push_back(dv);
443  }
444 
445  for (int i = 0; i < ni * 2; i++) // for each of the varibles
446  {
447  int var = floor(i / 2) + np - ni;
448  int tp = i % 2; // 0 is a u 1 is a v
449 
450  set<pair<int, int> >::iterator it;
451  for (it = spring.begin(); it != spring.end();
452  it++) // for each of the springs
453  {
454  Array<OneD, NekDouble> dr1, dr2;
455 
456  if ((*it).first == var)
457  {
458  if (tp == 0)
459  {
460  dr1 = dru[(*it).first];
461  }
462  else
463  {
464  dr1 = drv[(*it).first];
465  }
466  }
467  else
468  {
469  dr1 = Array<OneD, NekDouble>(3, 0.0);
470  }
471 
472  if ((*it).second == var)
473  {
474  if (tp == 0)
475  {
476  dr2 = dru[(*it).second];
477  }
478  else
479  {
480  dr2 = drv[(*it).second];
481  }
482  }
483  else
484  {
485  dr2 = Array<OneD, NekDouble>(3, 0.0);
486  }
487 
488  ret(i, 0) +=
489  2.0 / z[(*it)] *
490  Dot(Take(r[(*it).first], r[(*it).second]), Take(dr1, dr2));
491  }
492  }
493  return ret;
494 }
495 
496 void OptiFace::Update(Array<OneD, NekDouble> xinew)
497 {
498  for (int i = np - ni; i < np; i++)
499  {
500  uv[i][0] = xinew[2 * (i - np + ni) + 0];
501  uv[i][1] = xinew[2 * (i - np + ni) + 1];
502  }
503 }
504 }
505 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
class for CAD curves.
Definition: CADCurve.h:61
Array< OneD, NekDouble > Take(Array< OneD, NekDouble > a, Array< OneD, NekDouble > b)
STL namespace.
Array< OneD, NekDouble > Times(NekDouble t, Array< OneD, NekDouble > a)
NekDouble Dot(Array< OneD, NekDouble > a, Array< OneD, NekDouble > b)
class for handleing a cad surface
Definition: CADSurf.h:72
double NekDouble
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:217
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Array< OneD, NekDouble > Add(Array< OneD, NekDouble > a, Array< OneD, NekDouble > b)
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:168