Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  for (int i = 0; i < xci.num_elements(); i++)
96  {
97  if (xi[i] + d(i, 0) < li[i])
98  {
99  xci[i] = li[i];
100  Fset.erase(i);
101  continue;
102  }
103  else
104  {
105  xci[i] = xi[i] + d(i, 0);
106  }
107 
108  if (xi[i] + d(i, 0) > ui[i])
109  {
110  xci[i] = ui[i];
111  Fset.erase(i);
112  continue;
113  }
114  else
115  {
116  xci[i] = xi[i] + d(i, 0);
117  }
118  }
119 
120  DNekMat Z(xci.num_elements(), xci.num_elements(), 0.0);
121 
123  for (int i = 0; i < xci.num_elements(); i++)
124  {
125  it = Fset.find(i);
126  if (it != Fset.end())
127  {
128  Z(i, i) = 1.0;
129  }
130  }
131 
132  DNekMat dx(xci.num_elements(), 1, 0.0);
133  for (int i = 0; i < xci.num_elements(); i++)
134  {
135  dx(i, 0) = xci[i] - xi[i];
136  }
137 
138  DNekMat rg = Z * (J + B * dx);
139 
140  DNekMat du = -1.0 * H * rg;
141 
142  NekDouble alpha = 1.0;
143  for (it = Fset.begin(); it != Fset.end(); it++)
144  {
145  int i = (*it);
146  if (li[i] - xci[i] > alpha * du(i, 0))
147  {
148  alpha = min(alpha, (li[i] - xci[i]) / du(i, 0));
149  }
150  else if (ui[i] - xci[i] < alpha * du(i, 0))
151  {
152  alpha = min(alpha, (ui[i] - xci[i]) / du(i, 0));
153  }
154  }
155 
156  DNekMat grad = alpha * du;
157 
158  Array<OneD, NekDouble> dk(xci.num_elements()), xibar(xci.num_elements());
159  for (int i = 0; i < xci.num_elements(); i++)
160  {
161  set<int>::iterator f = Fset.find(i);
162  if (f != Fset.end())
163  {
164  xibar[i] = xci[i] + grad(i, 0);
165  }
166  else
167  {
168  xibar[i] = xci[i];
169  }
170  }
171 
172  Vmath::Vsub(xci.num_elements(), &xibar[0], 1, &xi[0], 1, &dk[0], 1);
173 
174  NekDouble c = 0.0;
175  NekDouble r = 0.0;
176  NekDouble l = 0.0;
177  for (int i = 0; i < dk.num_elements(); i++)
178  {
179  c += 1E-4 * J(i, 0) * dk[i];
180  r += J(i, 0) * dk[i];
181  }
182 
183  /*cout << endl << J << endl << endl;
184 
185  for(int i = 0; i < dk.num_elements(); i++)
186  {
187  cout << dk[i] << endl;
188  }
189  cout << endl;*/
190 
191  // this section needs a case evaluation for edges on faces
192  NekDouble lam = 2.0;
193  int iterct = 0;
194  NekDouble fo = opti->F(xi);
195  NekDouble fn;
196  Array<OneD, NekDouble> tst(xi.num_elements());
197  do
198  {
199  if (iterct > 100)
200  {
201  // cout << "failed line search" << endl;
202  return false;
203  }
204  iterct++;
205 
206  lam *= 0.5;
207 
208  for (int i = 0; i < xi.num_elements(); i++)
209  {
210  tst[i] = xi[i] + lam * dk[i];
211  }
212 
213  fn = opti->F(tst);
214 
215  DNekMat jn = opti->dF(tst);
216 
217  l = 0.0;
218  for (int i = 0; i < dk.num_elements(); i++)
219  {
220  l += jn(i, 0) * dk[i];
221  }
222 
223  } while (fn > fo + c || fabs(l) > 1.0 * fabs(r));
224  // wolfe conditions
225 
226  // tst at this point is the new all vector
227  // now need to update hessians
228  DNekMat Jn = opti->dF(tst);
229  DNekMat y = Jn - J;
230  DNekMat yT = Jn - J;
231  yT.Transpose();
232  DNekMat s(dk.num_elements(), 1, 0.0);
233  for (int i = 0; i < dk.num_elements(); i++)
234  {
235  s(i, 0) = lam * dk[i];
236  }
237  DNekMat sT = s;
238  sT.Transpose();
239 
240  DNekMat d1 = yT * s;
241  DNekMat d2 = sT * B * s;
242  DNekMat d3 = sT * y;
243  DNekMat n1 = yT * H * y;
244 
245  NekDouble ynorm = 0.0;
246  for (int i = 0; i < dk.num_elements(); i++)
247  {
248  ynorm += y(i, 0) * y(i, 0);
249  }
250 
251  if (d3(0, 0) > 2.2E-16 * ynorm)
252  {
253  B = B + y * yT * (1.0 / d1(0, 0)) - B * s * sT * B * (1.0 / d2(0, 0));
254  H = H + (d3(0, 0) + n1(0, 0)) / d3(0, 0) / d3(0, 0) * s * sT -
255  1.0 / d3(0, 0) * (H * y * sT + s * yT * H);
256  }
257 
258  J = Jn;
259  opti->Update(tst);
260 
261  return true;
262 }
263 }
264 }
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:343
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