Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BGFS-B.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 
38 
39 #include <limits>
40 #include <set>
41 
42 using namespace std;
43 namespace Nektar
44 {
45 namespace NekMeshUtils
46 {
47 
48 // this function will perform an update on the solution vector contained within
49 // opti
51 {
52 
53  Array<OneD, NekDouble> xi = opti->Getxi();
54  Array<OneD, NekDouble> ui = opti->Getui();
55  Array<OneD, NekDouble> li = opti->Getli();
56 
57  set<int> Fset;
58  Array<OneD, NekDouble> ti(xi.num_elements());
59  for (int i = 0; i < ti.num_elements(); i++)
60  {
61  if (J(i, 0) < 0)
62  {
63  ti[i] = (xi[i] - ui[i]) / J(i, 0);
64  }
65  else if (J(i, 0) > 0)
66  {
67  ti[i] = (xi[i] - li[i]) / J(i, 0);
68  }
69  else
70  {
71  ti[i] = numeric_limits<double>::max();
72  }
73  if (ti[i] > 0)
74  {
75  Fset.insert(i);
76  }
77  }
78 
79  // intitialise d
80  DNekMat d(xi.num_elements(), 1);
81  for (int i = 0; i < xi.num_elements(); i++)
82  {
83  if (fabs(ti[i]) < 1E-10)
84  {
85  d(i, 0) = 0;
86  }
87  else
88  {
89  d(i, 0) = -1.0 * J(i, 0);
90  }
91  }
92 
93  Array<OneD, NekDouble> xci(xi.num_elements());
94 
95  bool hitbounded = false;
96 
97  for (int i = 0; i < xci.num_elements(); i++)
98  {
99  if (xi[i] + d(i, 0) < li[i])
100  {
101  xci[i] = li[i];
102  Fset.erase(i);
103  cout << "hit bounded" << endl;
104  hitbounded = true;
105  continue;
106  }
107  else
108  {
109  xci[i] = xi[i] + d(i, 0);
110  }
111 
112  if (xi[i] + d(i, 0) > ui[i])
113  {
114  xci[i] = ui[i];
115  Fset.erase(i);
116  cout << "hit bounded" << endl;
117  hitbounded = true;
118  continue;
119  }
120  else
121  {
122  xci[i] = xi[i] + d(i, 0);
123  }
124  }
125 
126  DNekMat Z(xci.num_elements(), xci.num_elements(), 0.0);
127 
129  for (int i = 0; i < xci.num_elements(); i++)
130  {
131  it = Fset.find(i);
132  if (it != Fset.end())
133  {
134  Z(i, i) = 1.0;
135  }
136  }
137 
138  DNekMat dx(xci.num_elements(), 1, 0.0);
139  for (int i = 0; i < xci.num_elements(); i++)
140  {
141  dx(i, 0) = xci[i] - xi[i];
142  }
143 
144  DNekMat rg = Z * (J + B * dx);
145 
146  DNekMat du = -1.0 * H * rg;
147 
148  NekDouble alpha = 1.0;
149  for (it = Fset.begin(); it != Fset.end(); it++)
150  {
151  int i = (*it);
152  if (li[i] - xci[i] > alpha * du(i, 0))
153  {
154  alpha = min(alpha, (li[i] - xci[i]) / du(i, 0));
155  }
156  else if (ui[i] - xci[i] < alpha * du(i, 0))
157  {
158  alpha = min(alpha, (ui[i] - xci[i]) / du(i, 0));
159  }
160  }
161 
162  DNekMat grad = alpha * du;
163 
164  Array<OneD, NekDouble> dk(xci.num_elements()), xibar(xci.num_elements());
165  for (int i = 0; i < xci.num_elements(); i++)
166  {
167  set<int>::iterator f = Fset.find(i);
168  if (f != Fset.end())
169  {
170  xibar[i] = xci[i] + grad(i, 0);
171  }
172  else
173  {
174  xibar[i] = xci[i];
175  }
176  }
177 
178  /*if(hitbounded)
179  {
180  cout << endl << endl << Z << endl << endl;
181  cout << rg << endl << endl;
182  cout << du << endl << endl;
183  cout << alpha << endl << endl;
184  cout << grad << endl << endl;
185 
186  for(int i = 0; i < xi.num_elements(); i++)
187  {
188  cout << xci[i] << " " << xibar[i] << endl;
189  }
190 
191  exit(-1);
192  }*/
193 
194  Vmath::Vsub(xci.num_elements(), &xibar[0], 1, &xi[0], 1, &dk[0], 1);
195 
196  NekDouble c = 0.0;
197  NekDouble r = 0.0;
198  NekDouble l = 0.0;
199  for (int i = 0; i < dk.num_elements(); i++)
200  {
201  c += 1E-4 * J(i, 0) * dk[i];
202  r += J(i, 0) * dk[i];
203  }
204 
205  /*cout << endl << J << endl << endl;
206 
207  for(int i = 0; i < dk.num_elements(); i++)
208  {
209  cout << dk[i] << endl;
210  }
211  cout << endl;*/
212 
213  // this section needs a case evaluation for edges on faces
214  NekDouble lam = 2.0;
215  int iterct = 0;
216  NekDouble fo = opti->F(xi);
217  NekDouble fn;
218  Array<OneD, NekDouble> tst(xi.num_elements());
219  do
220  {
221  if (iterct > 20)
222  {
223  // cout << "failed line search" << endl;
224  return false;
225  }
226  iterct++;
227 
228  lam *= 0.5;
229 
230  for (int i = 0; i < xi.num_elements(); i++)
231  {
232  tst[i] = xi[i] + lam * dk[i];
233  }
234 
235  fn = opti->F(tst);
236 
237  DNekMat jn = opti->dF(tst);
238 
239  l = 0.0;
240  for (int i = 0; i < dk.num_elements(); i++)
241  {
242  l += jn(i, 0) * dk[i];
243  }
244 
245  } while (fn > fo + c || fabs(l) > 0.9 * fabs(r));
246  // wolfe conditions
247 
248  // tst at this point is the new all vector
249  // now need to update hessians
250  DNekMat Jn = opti->dF(tst);
251  DNekMat y = Jn - J;
252  DNekMat yT = Jn - J;
253  yT.Transpose();
254  DNekMat s(dk.num_elements(), 1, 0.0);
255  for (int i = 0; i < dk.num_elements(); i++)
256  {
257  s(i, 0) = lam * dk[i];
258  }
259  DNekMat sT = s;
260  sT.Transpose();
261 
262  DNekMat d1 = yT * s;
263  DNekMat d2 = sT * B * s;
264  DNekMat d3 = sT * y;
265  DNekMat n1 = yT * H * y;
266 
267  NekDouble ynorm = 0.0;
268  for (int i = 0; i < dk.num_elements(); i++)
269  {
270  ynorm += y(i, 0) * y(i, 0);
271  }
272 
273  if (d3(0, 0) > 2.2E-16 * ynorm)
274  {
275  B = B + y * yT * (1.0 / d1(0, 0)) - B * s * sT * B * (1.0 / d2(0, 0));
276  H = H + (d3(0, 0) + n1(0, 0)) / d3(0, 0) / d3(0, 0) * s * sT -
277  1.0 / d3(0, 0) * (H * y * sT + s * yT * H);
278  }
279 
280  J = Jn;
281  opti->Update(tst);
282 
283  return true;
284 }
285 }
286 }
STL namespace.
double NekDouble
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< OptiObj > OptiObjSharedPtr
Definition: OptimiseObj.h:90
bool BGFSUpdate(OptiObjSharedPtr opti, DNekMat &J, DNekMat &B, DNekMat &H)
Definition: BGFS-B.cpp:50