Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessVarOpti.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessJac.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 
38 
39 #include "ElUtil.h"
40 #include "NodeOpti.h"
41 #include "ProcessVarOpti.h"
43 
44 #include <boost/thread/mutex.hpp>
45 
46 #include <StdRegions/StdPrismExp.h>
47 #include <StdRegions/StdQuadExp.h>
48 #include <StdRegions/StdTetExp.h>
49 #include <StdRegions/StdTriExp.h>
50 
53 
54 using namespace std;
55 using namespace Nektar::NekMeshUtils;
56 
57 namespace Nektar
58 {
59 namespace Utilities
60 {
61 
62 ModuleKey ProcessVarOpti::className =
64  ModuleKey(eProcessModule, "varopti"), ProcessVarOpti::create,
65  "Optimise mesh locations.");
66 
67 ProcessVarOpti::ProcessVarOpti(MeshSharedPtr m) : ProcessModule(m)
68 {
69  // clang-format off
70  m_config["linearelastic"] =
71  ConfigOption(true, "", "Optimise for linear elasticity");
72  m_config["winslow"] =
73  ConfigOption(true, "", "Optimise for winslow");
74  m_config["roca"] =
75  ConfigOption(true, "", "Optimise for roca method");
76  m_config["hyperelastic"] =
77  ConfigOption(true, "", "Optimise for hyper elasticity");
78  m_config["numthreads"] =
79  ConfigOption(false, "1", "Number of threads");
80  m_config["restol"] =
81  ConfigOption(false, "1e-6", "Tolerance criterion");
82  m_config["maxiter"] =
83  ConfigOption(false, "500", "Maximum number of iterations");
84  m_config["nq"] =
85  ConfigOption(false, "-1", "Order of mesh");
86  m_config["region"] =
87  ConfigOption(false, "0.0", "create regions based on target");
88  m_config["resfile"] =
89  ConfigOption(false, "", "writes residual values to file");
90  m_config["histfile"] =
91  ConfigOption(false, "", "histogram of scaled jac");
92  m_config["overint"] =
93  ConfigOption(false, "6", "over integration order");
94  m_config["analytics"] =
95  ConfigOption(false, "", "basic analytics module");
96  // clang-format on
97 }
98 
100 {
101 }
102 
104 {
105  if (m_mesh->m_verbose)
106  {
107  cout << "ProcessVarOpti: Optimising... " << endl;
108  }
109 
110  if (m_config["linearelastic"].beenSet)
111  {
112  m_opti = eLinEl;
113  }
114  else if (m_config["winslow"].beenSet)
115  {
116  m_opti = eWins;
117  }
118  else if (m_config["roca"].beenSet)
119  {
120  m_opti = eRoca;
121  }
122  else if (m_config["hyperelastic"].beenSet)
123  {
124  m_opti = eHypEl;
125  }
126  else
127  {
128  ASSERTL0(false, "not opti type set");
129  }
130 
131  const int maxIter = m_config["maxiter"].as<int>();
132  const NekDouble restol = m_config["restol"].as<NekDouble>();
133 
134  // m_mesh->m_nummode = m_config["nq"].as<int>();
135 
136  EdgeSet::iterator eit;
137  bool fd = false;
138 
139  if (m_config["nq"].beenSet)
140  {
141  m_mesh->m_nummode = m_config["nq"].as<int>();
142  fd = true;
143  }
144 
145  if (!fd)
146  {
147  for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end();
148  eit++)
149  {
150  if ((*eit)->m_edgeNodes.size() > 0)
151  {
152  m_mesh->m_nummode = (*eit)->m_edgeNodes.size() + 2;
153  fd = true;
154  break;
155  }
156  }
157  }
158  ASSERTL0(fd, "failed to find order of mesh");
159 
160  // Safety feature: limit over-integration order for high-order triangles
161  // over order 5.
162  int intOrder = m_config["overint"].as<int>();
163  intOrder = m_mesh->m_nummode + intOrder <= 11 ?
164  intOrder : 11 - m_mesh->m_nummode;
165 
166  if (m_mesh->m_verbose)
167  {
168  cout << "Identified mesh order as: " << m_mesh->m_nummode - 1 << endl;
169  }
170 
171  if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
172  {
173  ASSERTL0(false, "cannot deal with manifolds");
174  }
175 
176  m_res = boost::shared_ptr<Residual>(new Residual);
177  m_res->val = 1.0;
178 
179 
180  m_mesh->MakeOrder(m_mesh->m_nummode - 1,
182 
183  if (m_config["analytics"].beenSet)
184  {
185  Analytics();
186  return;
187  }
188 
189  map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
190  BuildDerivUtil(intOrder);
191 
192  GetElementMap(intOrder, derivUtils);
193 
194  m_res->startInv = 0;
195  m_res->worstJac = numeric_limits<double>::max();
196  for (int i = 0; i < m_dataSet.size(); i++)
197  {
198  m_dataSet[i]->Evaluate();
199  m_dataSet[i]->InitialMinJac();
200  }
201 
202  vector<ElUtilSharedPtr> elLock;
203 
204  if (m_config["region"].beenSet)
205  {
206  elLock = GetLockedElements(m_config["region"].as<NekDouble>());
207  }
208 
209 
210  vector<vector<NodeSharedPtr> > freenodes = GetColouredNodes(elLock);
211  vector<vector<NodeOptiSharedPtr> > optiNodes;
212 
213  // turn the free nodes into optimisable objects with all required data
214  set<int> check;
215  for (int i = 0; i < freenodes.size(); i++)
216  {
217  vector<NodeOptiSharedPtr> ns;
218  for (int j = 0; j < freenodes[i].size(); j++)
219  {
220  NodeElMap::iterator it = m_nodeElMap.find(freenodes[i][j]->m_id);
221  ASSERTL0(it != m_nodeElMap.end(), "could not find");
222 
223  int optiKind = m_mesh->m_spaceDim;
224 
225  if (freenodes[i][j]->GetNumCadCurve() == 1)
226  {
227  optiKind += 10;
228  }
229  else if (freenodes[i][j]->GetNumCADSurf() == 1)
230  {
231  optiKind += 20;
232  }
233  else
234  {
235  optiKind += 10 * m_mesh->m_expDim;
236  }
237 
238  set<int>::iterator c = check.find(freenodes[i][j]->m_id);
239  ASSERTL0(c == check.end(), "duplicate node");
240  check.insert(freenodes[i][j]->m_id);
241 
242  ns.push_back(GetNodeOptiFactory().CreateInstance(
243  optiKind, freenodes[i][j], it->second, m_res, derivUtils,
244  m_opti));
245  }
246  optiNodes.push_back(ns);
247  }
248 
249  int nset = optiNodes.size();
250  int p = 0;
251  int mn = numeric_limits<int>::max();
252  int mx = 0;
253  for (int i = 0; i < nset; i++)
254  {
255  p += optiNodes[i].size();
256  mn = min(mn, int(optiNodes[i].size()));
257  mx = max(mx, int(optiNodes[i].size()));
258  }
259 
260  if (m_config["histfile"].beenSet)
261  {
262  ofstream histFile;
263  string name = m_config["histfile"].as<string>() + "_start.txt";
264  histFile.open(name.c_str());
265 
266  for (int i = 0; i < m_dataSet.size(); i++)
267  {
268  histFile << m_dataSet[i]->GetScaledJac() << endl;
269  }
270  histFile.close();
271  }
272 
273  if(m_mesh->m_verbose)
274  {
275  cout << scientific << endl;
276  cout << "N elements:\t\t" << m_mesh->m_element[m_mesh->m_expDim].size() - elLock.size() << endl
277  << "N elements invalid:\t" << m_res->startInv << endl
278  << "Worst jacobian:\t\t" << m_res->worstJac << endl
279  << "N free nodes:\t\t" << m_res->n << endl
280  << "N Dof:\t\t\t" << m_res->nDoF << endl
281  << "N color sets:\t\t" << nset << endl
282  << "Avg set colors:\t\t" << p/nset << endl
283  << "Min set:\t\t" << mn << endl
284  << "Max set:\t\t" << mx << endl
285  << "Residual tolerance:\t" << restol << endl;
286  }
287 
288  int nThreads = m_config["numthreads"].as<int>();
289 
290  int ctr = 0;
292  tms.SetThreadingType("ThreadManagerBoost");
295 
296  Timer t;
297  t.Start();
298 
299  ofstream resFile;
300  if (m_config["resfile"].beenSet)
301  {
302  resFile.open(m_config["resfile"].as<string>().c_str());
303  }
304 
305  for (int i = 0; i < optiNodes.size(); i++)
306  {
307  vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
308  for (int j = 0; j < optiNodes[i].size(); j++)
309  {
310  optiNodes[i][j]->CalcMinJac();
311  }
312  }
313 
314  while (m_res->val > restol)
315  {
316  ctr++;
317  m_res->val = 0.0;
318  m_res->func = 0.0;
319  m_res->nReset[0] = 0;
320  m_res->nReset[1] = 0;
321  m_res->nReset[2] = 0;
322  m_res->alphaI = 0;
323  for (int i = 0; i < optiNodes.size(); i++)
324  {
325  vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
326  for (int j = 0; j < optiNodes[i].size(); j++)
327  {
328  jobs[j] = optiNodes[i][j]->GetJob();
329  }
330 
331  tm->SetNumWorkers(0);
332  tm->QueueJobs(jobs);
333  tm->SetNumWorkers(nThreads);
334  tm->Wait();
335  }
336 
337  m_res->startInv = 0;
338  m_res->worstJac = numeric_limits<double>::max();
339 
340  vector<Thread::ThreadJob *> elJobs(m_dataSet.size());
341  for (int i = 0; i < m_dataSet.size(); i++)
342  {
343  elJobs[i] = m_dataSet[i]->GetJob();
344  }
345 
346  tm->SetNumWorkers(0);
347  tm->QueueJobs(elJobs);
348  tm->SetNumWorkers(nThreads);
349  tm->Wait();
350 
351  if (m_config["resfile"].beenSet)
352  {
353  resFile << m_res->val << " " << m_res->worstJac << " "
354  << m_res->func << endl;
355  }
356 
357  if(m_mesh->m_verbose)
358  {
359  cout << ctr
360  << "\tResidual: " << m_res->val
361  << "\tMin Jac: " << m_res->worstJac
362  << "\tInvalid: " << m_res->startInv
363  << "\tReset nodes: " << m_res->nReset[0] << "/" << m_res->nReset[1]
364  << "/" << m_res->nReset[2] << "\tFunctional: " << m_res->func
365  << endl;
366  }
367 
368  if(ctr >= maxIter)
369  {
370  break;
371  }
372  }
373 
374  if (m_config["histfile"].beenSet)
375  {
376  ofstream histFile;
377  string name = m_config["histfile"].as<string>() + "_end.txt";
378  histFile.open(name.c_str());
379 
380  for (int i = 0; i < m_dataSet.size(); i++)
381  {
382  histFile << m_dataSet[i]->GetScaledJac() << endl;
383  }
384  histFile.close();
385  }
386  if (m_config["resfile"].beenSet)
387  {
388  resFile.close();
389  }
390 
391  t.Stop();
392 
393  //RemoveLinearCurvature();
394 
395  if(m_mesh->m_verbose)
396  {
397  cout << "Time to compute: " << t.TimePerTest(1) << endl;
398  cout << "Invalid at end: " << m_res->startInv << endl;
399  cout << "Worst at end: " << m_res->worstJac << endl;
400  }
401 }
402 
404 {
405 public:
408  : NodalUtilTriangle(degree, r, s)
409  {
410  }
411 
413  {
414  }
415 
416 protected:
417  virtual NekVector<NekDouble> v_OrthoBasis(const int mode)
418  {
419  // Monomial basis.
420  std::pair<int, int> modes = m_ordering[mode];
422 
423  for (int i = 0; i < m_numPoints; ++i)
424  {
425  ret(i) =
426  pow(m_xi[0][i], modes.first) * pow(m_xi[1][i], modes.second);
427  }
428 
429  return ret;
430  }
431 
433  const int mode)
434  {
435  ASSERTL0(false, "not supported");
436  return NekVector<NekDouble>();
437  }
438 
439  virtual boost::shared_ptr<NodalUtil> v_CreateUtil(
441  {
443  m_degree, xi[0], xi[1]);
444  }
445 };
446 
448 {
449  // Grab the first element from the list
450  ElementSharedPtr elmt = m_mesh->m_element[m_mesh->m_expDim][0];
451 
452  // Get curved nodes
453  vector<NodeSharedPtr> nodes;
454  elmt->GetCurvedNodes(nodes);
455 
456  // We're going to investigate only the first node (corner node)
457  NodeSharedPtr node = nodes[4];
458 
459  // Loop over overintegration orders
460  const int nPoints = 200;
461  const int overInt = 40;
462  const NekDouble originX = -1.0;
463  const NekDouble originY = -1.0;
464  const NekDouble length = 2.0;
465  const NekDouble dx = length / (nPoints - 1);
466 
467  cout << "# overint = " << overInt << endl;
468  cout << "# Columns: x, y, over-integration orders (0 -> " << overInt - 1
469  << "), "
470  << " min(scaledJac)" << endl;
471 
472  // Loop over square defined by (originX, originY), length
473  for (int k = 0; k < nPoints; ++k)
474  {
475  node->m_y = originY + k * dx;
476  for (int j = 0; j < nPoints; ++j)
477  {
478  node->m_x = originX + j * dx;
479  cout << node->m_x << " " << node->m_y << " ";
480 
481  NekDouble minJacNew;
482 
483  for (int i = 0; i < overInt; ++i)
484  {
485  // Clear any existing node to element mapping.
486  m_dataSet.clear();
487  m_nodeElMap.clear();
488 
489  // Build deriv utils and element map.
490  map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
491  BuildDerivUtil(i);
492 
493  // Reconstruct element map
494  GetElementMap(i, derivUtils);
495 
496  for (int j = 0; j < m_dataSet.size(); j++)
497  {
498  m_dataSet[j]->Evaluate();
499  m_dataSet[j]->InitialMinJac();
500  }
501 
502  // Create NodeOpti object.
503  NodeOptiSharedPtr nodeOpti =
505  m_mesh->m_spaceDim * 11, node,
506  m_nodeElMap.find(node->m_id)->second, m_res, derivUtils,
507  m_opti);
508 
509  minJacNew = 0.0;
510 
511  // Evaluate functional.
512  nodeOpti->CalcMinJac();
513  cout << nodeOpti->GetFunctional<2>(minJacNew) << " ";
514  // NekDouble eigen;
515  // nodeOpti->GetFunctional<2>(minJacNew);
516  // nodeOpti->MinEigen<2>(eigen);
517  // cout << eigen << " ";
518  }
519 
520  cout << minJacNew << endl;
521  }
522  }
523 }
524 }
525 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
NodeOptiFactory & GetNodeOptiFactory()
Definition: NodeOpti.cpp:52
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Specialisation of the NodalUtil class to support nodal triangular elements.
Definition: NodalUtil.h:170
virtual NekVector< NekDouble > v_OrthoBasis(const int mode)
Return the value of the modal functions for the triangular element at the nodal points m_xi for a giv...
std::vector< ElUtilSharedPtr > GetLockedElements(NekDouble thres)
void GetElementMap(int o, std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > derMap)
std::map< LibUtilities::ShapeType, DerivUtilSharedPtr > BuildDerivUtil(int o)
std::vector< ElUtilSharedPtr > m_dataSet
STL namespace.
pair< ModuleType, string > ModuleKey
virtual boost::shared_ptr< NodalUtil > v_CreateUtil(Array< OneD, Array< OneD, NekDouble > > &xi)
Construct a NodalUtil object of the appropriate element type for a given set of points.
std::vector< std::vector< NodeSharedPtr > > GetColouredNodes(std::vector< ElUtilSharedPtr > elLock)
NodalUtilTriMonomial(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s)
void Stop()
Definition: Timer.cpp:62
std::vector< std::pair< int, int > > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
Definition: NodalUtil.h:184
int m_degree
Degree of the nodal element.
Definition: NodalUtil.h:108
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< NodeOpti > NodeOptiSharedPtr
Definition: NodeOpti.h:135
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t, unsigned int nThr)
Creates an instance of a ThreadManager (which one is determined by a previous call to SetThreadingTyp...
Definition: Thread.cpp:200
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Abstract base class for processing modules.
virtual NekVector< NekDouble > v_OrthoBasisDeriv(const int dir, const int mode)
Return the value of the derivative of the modal functions for the triangular element at the nodal poi...
Array< OneD, Array< OneD, NekDouble > > m_xi
Coordinates of the nodal points defining the basis.
Definition: NodalUtil.h:112
void SetThreadingType(const std::string &p_type)
Sets what ThreadManagers will be created in CreateInstance.
Definition: Thread.cpp:162
int m_numPoints
Total number of nodal points.
Definition: NodalUtil.h:110
NodalUtilTriangle(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s)
Construct the nodal utility class for a triangle.
Definition: NodalUtil.cpp:239
boost::shared_ptr< Element > ElementSharedPtr
Definition: Edge.h:49
void Start()
Definition: Timer.cpp:51
std::pair< ModuleType, std::string > ModuleKey
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:108
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
boost::shared_ptr< ThreadManager > ThreadManagerSharedPtr
Definition: Thread.h:74