Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
NodeOptiCAD.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: NodeOpti.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: Calculate jacobians of elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "NodeOptiCAD.h"
37 #include "Evaluator.hxx"
38 
39 using namespace std;
40 using namespace Nektar::NekMeshUtils;
41 
42 namespace Nektar
43 {
44 namespace Utilities
45 {
46 
48 
49 int NodeOpti1D3D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
50  13, NodeOpti1D3D::create, "1D3D");
51 
52 void NodeOpti1D3D::Optimise()
53 {
54  NekDouble minJacNew;
55  NekDouble currentW = GetFunctional<3>(minJacNew);
56  NekDouble newVal = currentW;
57 
58  if (m_grad[0] * m_grad[0] + m_grad[1] * m_grad[1] + m_grad[2] * m_grad[2] >
59  gradTol())
60  {
61  // modify the gradient to be on the cad system
62  ProcessGradient();
63 
64  // needs to optimise
65  NekDouble tc = m_node->GetCADCurveInfo(curve->GetId());
66  NekDouble xc = m_node->m_x;
67  NekDouble yc = m_node->m_y;
68  NekDouble zc = m_node->m_z;
69  NekDouble nt;
71 
73 
74  if (m_grad[1] < 1e-6)
75  {
76  m_grad[1] = 1e-6 - m_grad[1];
77  }
78 
79  sk[0] = m_grad[0] / m_grad[1] * -1.0;
80 
81  Array<OneD, NekDouble> bd = curve->GetBounds();
82 
83  bool found = false;
84 
85  NekDouble pg = m_grad[0] * sk[0];
86 
87  // normal gradient line Search
88  NekDouble alpha = 1.0;
89 
90  while (alpha > alphaTol())
91  {
92  // Update node
93  nt = tc + alpha * sk[0];
94  if (nt < bd[0] || nt > bd[1])
95  {
96  alpha /= 2.0;
97  continue;
98  }
99 
100  p = curve->P(nt);
101  m_node->m_x = p[0];
102  m_node->m_y = p[1];
103  m_node->m_z = p[2];
104 
105  newVal = GetFunctional<3>(minJacNew, false);
106 
107  // Wolfe conditions
108  if (newVal <= currentW + c1() * alpha * pg)
109  {
110  found = true;
111  break;
112  }
113 
114  alpha /= 2.0;
115  }
116 
117  if (!found)
118  {
119  // reset the node
120  nt = tc;
121  p = curve->P(nt);
122  m_node->m_x = p[0];
123  m_node->m_y = p[1];
124  m_node->m_z = p[2];
125 
126  mtx.lock();
127  m_res->nReset[0]++;
128  mtx.unlock();
129  }
130  else
131  {
132  m_minJac = minJacNew;
133  m_node->MoveCurve(p, curve->GetId(), nt);
134  }
135  mtx.lock();
136  m_res->val = max(sqrt((m_node->m_x - xc) * (m_node->m_x - xc) +
137  (m_node->m_y - yc) * (m_node->m_y - yc) +
138  (m_node->m_z - zc) * (m_node->m_z - zc)),
139  m_res->val);
140  mtx.unlock();
141  }
142 
143  mtx.lock();
144  m_res->func += newVal;
145  mtx.unlock();
146 }
147 
148 int NodeOpti2D3D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
149  23, NodeOpti2D3D::create, "2D3D");
150 
151 void NodeOpti2D3D::Optimise()
152 {
153  NekDouble minJacNew;
154  NekDouble currentW = GetFunctional<3>(minJacNew);
155  NekDouble newVal = currentW;
156 
157  if (m_grad[0] * m_grad[0] + m_grad[1] * m_grad[1] + m_grad[2] * m_grad[2] >
158  gradTol())
159  {
160  // modify the gradient to be on the cad system
161  ProcessGradient();
162 
163  // needs to optimise
164  Array<OneD, NekDouble> uvc = m_node->GetCADSurfInfo(surf->GetId());
165  NekDouble xc = m_node->m_x;
166  NekDouble yc = m_node->m_y;
167  NekDouble zc = m_node->m_z;
168  Array<OneD, NekDouble> uvt(2);
170  Array<OneD, NekDouble> bd = surf->GetBounds();
171 
173  NekDouble val;
174 
175  // Calculate minimum eigenvalue
176  MinEigen<2>(val);
177 
178  if (val < 1e-6)
179  {
180  // Add constant identity to Hessian matrix.
181  m_grad[2] += 1e-6 - val;
182  m_grad[4] += 1e-6 - val;
183  }
184 
185  sk[0] = -1.0 / (m_grad[2] * m_grad[4] - m_grad[3] * m_grad[3]) *
186  (m_grad[4] * m_grad[0] - m_grad[3] * m_grad[1]);
187  sk[1] = -1.0 / (m_grad[2] * m_grad[4] - m_grad[3] * m_grad[3]) *
188  (m_grad[2] * m_grad[1] - m_grad[3] * m_grad[0]);
189 
190  bool found = false;
191  NekDouble pg = (m_grad[0] * sk[0] + m_grad[1] * sk[1]);
192 
193  // normal gradient line Search
194  NekDouble alpha = 1.0;
195  while (alpha > alphaTol())
196  {
197  uvt[0] = uvc[0] + alpha * sk[0];
198  uvt[1] = uvc[1] + alpha * sk[1];
199 
200  if (uvt[0] < bd[0] || uvt[0] > bd[1] || uvt[1] < bd[2] ||
201  uvt[1] > bd[3])
202  {
203  alpha /= 2.0;
204  continue;
205  }
206 
207  p = surf->P(uvt);
208  m_node->m_x = p[0];
209  m_node->m_y = p[1];
210  m_node->m_z = p[2];
211 
212  newVal = GetFunctional<3>(minJacNew, false);
213 
214  // Wolfe conditions
215  if (newVal <= currentW + c1() * (alpha * pg))
216  {
217  found = true;
218  break;
219  }
220 
221  alpha /= 2.0;
222  }
223 
224  if (!found)
225  {
226  // reset the node
227  p = surf->P(uvc);
228  m_node->m_x = p[0];
229  m_node->m_y = p[1];
230  m_node->m_z = p[2];
231 
232  mtx.lock();
233  m_res->nReset[1]++;
234  mtx.unlock();
235  }
236  else
237  {
238  m_minJac = minJacNew;
239  m_node->Move(p, surf->GetId(), uvt);
240  }
241 
242  mtx.lock();
243  m_res->val = max(sqrt((m_node->m_x - xc) * (m_node->m_x - xc) +
244  (m_node->m_y - yc) * (m_node->m_y - yc) +
245  (m_node->m_z - zc) * (m_node->m_z - zc)),
246  m_res->val);
247  mtx.unlock();
248 
249  Array<OneD, NekDouble> uva = m_node->GetCADSurfInfo(surf->GetId());
250  if (uva[0] < bd[0] || uva[0] > bd[1] || uva[1] < bd[2] ||
251  uva[1] > bd[3])
252  {
253  ASSERTL0(false, "something when very wrong and the node finished "
254  "out of its parameter plane");
255  }
256  }
257 
258  mtx.lock();
259  m_res->func += newVal;
260  mtx.unlock();
261 }
262 
263 void NodeOpti1D3D::ProcessGradient()
264 {
265  NekDouble tc = m_node->GetCADCurveInfo(curve->GetId());
266  Array<OneD, NekDouble> grad = m_grad;
267  m_grad = Array<OneD, NekDouble>(2, 0.0);
268 
269  // Grab first and second order CAD derivatives
270  Array<OneD, NekDouble> d2 = curve->D2(tc);
271 
272  // Multiply gradient by derivative of CAD
273  m_grad[0] = grad[0] * d2[3] + grad[1] * d2[4] + grad[2] * d2[5];
274 
275  // Second order: product rule of above, so multiply gradient by second
276  // order CAD derivatives and Hessian by gradient of CAD
277  m_grad[1] = grad[0] * d2[6] + grad[1] * d2[7] + grad[2] * d2[8] +
278  d2[3] * (grad[3] * d2[3] + grad[4] * d2[4] + grad[5] * d2[5]) +
279  d2[4] * (grad[4] * d2[3] + grad[6] * d2[4] + grad[7] * d2[5]) +
280  d2[5] * (grad[5] * d2[3] + grad[7] * d2[4] + grad[8] * d2[5]);
281 }
282 
283 void NodeOpti2D3D::ProcessGradient()
284 {
285  Array<OneD, NekDouble> uvc = m_node->GetCADSurfInfo(surf->GetId());
286 
287  Array<OneD, NekDouble> grad = m_grad;
288  m_grad = Array<OneD, NekDouble>(5, 0.0);
289 
290  Array<OneD, NekDouble> d2 = surf->D2(uvc);
291  // r[0] x
292  // r[1] y
293  // r[2] z
294  // r[3] dx/du a
295  // r[4] dy/du b
296  // r[5] dz/du c
297  // r[6] dx/dv d
298  // r[7] dy/dv e
299  // r[8] dz/dv f
300  // r[9] d2x/du2
301  // r[10] d2y/du2
302  // r[11] d2z/du2
303  // r[12] d2x/dv2
304  // r[13] d2y/dv2
305  // r[14] d2z/dv2
306  // r[15] d2x/dudv
307  // r[16] d2y/dudv
308  // r[17] d2z/dudv
309  //
310  // grad[0] d/dx
311  // grad[1] d/dy
312  // grad[2] d/dz
313  // grad[3] d2/dx2
314  // grad[4] d2/dxdy
315  // grad[5] d2/dxdz
316  // grad[6] d2/dy2
317  // grad[7] d2/dydz
318  // grad[8] d2/dz2
319 
320  // Gradients
321  m_grad[0] = d2[3] * grad[0] + d2[4] * grad[1] + d2[5] * grad[2];
322  m_grad[1] = d2[6] * grad[0] + d2[7] * grad[1] + d2[8] * grad[2];
323 
324  m_grad[2] = grad[0] * d2[9] + grad[1] * d2[10] + grad[2] * d2[11] +
325  grad[3] * d2[3] * d2[3] + grad[6] * d2[4] * d2[4] +
326  grad[8] * d2[5] * d2[5] + 2.0 * grad[4] * d2[3] * d2[4] +
327  2.0 * grad[5] * d2[3] * d2[5] + 2.0 * grad[7] * d2[4] * d2[5];
328 
329  m_grad[4] = grad[0] * d2[12] + grad[1] * d2[13] + grad[2] * d2[14] +
330  grad[3] * d2[6] * d2[6] + grad[6] * d2[7] * d2[7] +
331  grad[8] * d2[8] * d2[8] + 2.0 * grad[4] * d2[6] * d2[7] +
332  2.0 * grad[5] * d2[6] * d2[8] + 2.0 * grad[7] * d2[7] * d2[8];
333 
334  m_grad[3] = grad[0] * d2[15] + grad[1] * d2[16] + grad[2] * d2[17] +
335  grad[3] * d2[3] * d2[6] + grad[6] * d2[4] * d2[7] +
336  grad[8] * d2[5] * d2[8] +
337  grad[4] * (d2[3] * d2[7] + d2[4] * d2[6]) +
338  grad[5] * (d2[3] * d2[8] + d2[5] * d2[6]) +
339  grad[7] * (d2[4] * d2[8] + d2[5] * d2[7]);
340 }
341 }
342 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
NodeOptiFactory & GetNodeOptiFactory()
Definition: NodeOpti.cpp:52
STL namespace.
boost::mutex mtx
Definition: NodeOptiCAD.cpp:47
double NekDouble
static boost::mutex mutex
Definition: Vmath.cpp:134
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215