Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 CADType::eCurve:
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 CADType::eSurf:
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 CADType::eVert:
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 CADType::eCurve:
109  li = Array<OneD, NekDouble>(all.num_elements() - 2);
110  bnds = boost::dynamic_pointer_cast<CADCurve>(o)->GetBounds();
111  for (int i = 1; i < all.num_elements() - 1; i++)
112  {
113  li[i - 1] = bnds[0];
114  }
115  break;
116 
117  case CADType::eSurf:
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 CADType::eVert:
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 CADType::eCurve:
146  ui = Array<OneD, NekDouble>(all.num_elements() - 2);
147  bnds = boost::dynamic_pointer_cast<CADCurve>(o)->GetBounds();
148  for (int i = 1; i < all.num_elements() - 1; i++)
149  {
150  ui[i - 1] = bnds[1];
151  }
152  break;
153 
154  case CADType::eSurf:
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 CADType::eVert:
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() == CADType::eCurve)
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() == CADType::eSurf)
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() == CADType::eCurve)
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() == CADType::eSurf)
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() == CADType::eCurve)
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() == CADType::eSurf)
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() == CADType::eCurve)
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() == CADType::eSurf)
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() == CADType::eCurve)
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() == CADType::eSurf)
352  {
353  for (int i = 0; i < xinew.num_elements(); i++)
354  {
355  all[i + 2] = xinew[i];
356  }
357  }
358 }
359 }
360 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
base class for CAD curves.
Definition: CADCurve.h:52
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)
base class for a cad surface
Definition: CADSurf.h:55
double NekDouble
boost::shared_ptr< CADSurf > CADSurfSharedPtr
Definition: CADSurf.h:172
Array< OneD, NekDouble > Add(Array< OneD, NekDouble > a, Array< OneD, NekDouble > b)
boost::shared_ptr< CADCurve > CADCurveSharedPtr
Definition: CADCurve.h:194