44 #include <boost/thread/mutex.hpp>
65 "Optimise mesh locations.");
71 ConfigOption(
true,
"",
"Optimise for linear elasticity");
83 ConfigOption(
false,
"500",
"Maximum number of iterations");
87 ConfigOption(
false,
"0.0",
"create regions based on target");
89 ConfigOption(
false,
"",
"writes residual values to file");
107 cout <<
"ProcessVarOpti: Optimising... " << endl;
110 if (
m_config[
"linearelastic"].beenSet)
114 else if (
m_config[
"winslow"].beenSet)
122 else if (
m_config[
"hyperelastic"].beenSet)
128 ASSERTL0(
false,
"not opti type set");
131 const int maxIter =
m_config[
"maxiter"].as<
int>();
147 for (eit =
m_mesh->m_edgeSet.begin(); eit !=
m_mesh->m_edgeSet.end();
150 if ((*eit)->m_edgeNodes.size() > 0)
152 m_mesh->m_nummode = (*eit)->m_edgeNodes.size() + 2;
158 ASSERTL0(fd,
"failed to find order of mesh");
162 int intOrder =
m_config[
"overint"].as<
int>();
163 intOrder =
m_mesh->m_nummode + intOrder <= 11 ?
164 intOrder : 11 -
m_mesh->m_nummode;
168 cout <<
"Identified mesh order as: " <<
m_mesh->m_nummode - 1 << endl;
173 ASSERTL0(
false,
"cannot deal with manifolds");
189 map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
195 m_res->worstJac = numeric_limits<double>::max();
196 for (
int i = 0; i <
m_dataSet.size(); i++)
202 vector<ElUtilSharedPtr> elLock;
211 vector<vector<NodeOptiSharedPtr> > optiNodes;
215 for (
int i = 0; i < freenodes.size(); i++)
217 vector<NodeOptiSharedPtr> ns;
218 for (
int j = 0; j < freenodes[i].size(); j++)
223 int optiKind =
m_mesh->m_spaceDim;
225 if (freenodes[i][j]->GetNumCadCurve() == 1)
229 else if (freenodes[i][j]->GetNumCADSurf() == 1)
235 optiKind += 10 *
m_mesh->m_expDim;
239 ASSERTL0(c == check.end(),
"duplicate node");
240 check.insert(freenodes[i][j]->m_id);
243 optiKind, freenodes[i][j], it->second,
m_res, derivUtils,
246 optiNodes.push_back(ns);
249 int nset = optiNodes.size();
251 int mn = numeric_limits<int>::max();
253 for (
int i = 0; i < nset; i++)
255 p += optiNodes[i].size();
256 mn = min(mn,
int(optiNodes[i].size()));
257 mx = max(mx,
int(optiNodes[i].size()));
263 string name =
m_config[
"histfile"].as<
string>() +
"_start.txt";
264 histFile.open(name.c_str());
266 for (
int i = 0; i <
m_dataSet.size(); i++)
268 histFile <<
m_dataSet[i]->GetScaledJac() << endl;
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;
288 int nThreads =
m_config[
"numthreads"].as<
int>();
302 resFile.open(
m_config[
"resfile"].as<string>().c_str());
305 for (
int i = 0; i < optiNodes.size(); i++)
307 vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
308 for (
int j = 0; j < optiNodes[i].size(); j++)
310 optiNodes[i][j]->CalcMinJac();
314 while (
m_res->val > restol)
319 m_res->nReset[0] = 0;
320 m_res->nReset[1] = 0;
321 m_res->nReset[2] = 0;
323 for (
int i = 0; i < optiNodes.size(); i++)
325 vector<Thread::ThreadJob *> jobs(optiNodes[i].size());
326 for (
int j = 0; j < optiNodes[i].size(); j++)
328 jobs[j] = optiNodes[i][j]->GetJob();
331 tm->SetNumWorkers(0);
333 tm->SetNumWorkers(nThreads);
338 m_res->worstJac = numeric_limits<double>::max();
340 vector<Thread::ThreadJob *> elJobs(
m_dataSet.size());
341 for (
int i = 0; i <
m_dataSet.size(); i++)
346 tm->SetNumWorkers(0);
347 tm->QueueJobs(elJobs);
348 tm->SetNumWorkers(nThreads);
353 resFile <<
m_res->val <<
" " <<
m_res->worstJac <<
" "
354 <<
m_res->func << endl;
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
377 string name =
m_config[
"histfile"].as<
string>() +
"_end.txt";
378 histFile.open(name.c_str());
380 for (
int i = 0; i <
m_dataSet.size(); i++)
382 histFile <<
m_dataSet[i]->GetScaledJac() << endl;
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;
426 pow(
m_xi[0][i], modes.first) * pow(
m_xi[1][i], modes.second);
453 vector<NodeSharedPtr> nodes;
454 elmt->GetCurvedNodes(nodes);
460 const int nPoints = 200;
461 const int overInt = 40;
465 const NekDouble dx = length / (nPoints - 1);
467 cout <<
"# overint = " << overInt << endl;
468 cout <<
"# Columns: x, y, over-integration orders (0 -> " << overInt - 1
470 <<
" min(scaledJac)" << endl;
473 for (
int k = 0; k < nPoints; ++k)
475 node->m_y = originY + k * dx;
476 for (
int j = 0; j < nPoints; ++j)
478 node->m_x = originX + j * dx;
479 cout << node->m_x <<
" " << node->m_y <<
" ";
483 for (
int i = 0; i < overInt; ++i)
490 map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
496 for (
int j = 0; j <
m_dataSet.size(); j++)
505 m_mesh->m_spaceDim * 11, node,
512 nodeOpti->CalcMinJac();
513 cout << nodeOpti->GetFunctional<2>(minJacNew) <<
" ";
520 cout << minJacNew << endl;
#define ASSERTL0(condition, msg)
virtual ~NodalUtilTriMonomial()
NodeOptiFactory & GetNodeOptiFactory()
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.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Specialisation of the NodalUtil class to support nodal triangular elements.
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
pair< ModuleType, string > ModuleKey
MeshSharedPtr m_mesh
Mesh object.
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)
std::vector< std::pair< int, int > > m_ordering
Mapping from the indexing of the basis to a continuous ordering.
int m_degree
Degree of the nodal element.
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< NodeOpti > NodeOptiSharedPtr
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...
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
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.
void SetThreadingType(const std::string &p_type)
Sets what ThreadManagers will be created in CreateInstance.
int m_numPoints
Total number of nodal points.
NodalUtilTriangle(int degree, Array< OneD, NekDouble > r, Array< OneD, NekDouble > s)
Construct the nodal utility class for a triangle.
boost::shared_ptr< Element > ElementSharedPtr
std::pair< ModuleType, std::string > ModuleKey
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
virtual ~ProcessVarOpti()
1D Gauss-Lobatto-Legendre quadrature points
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
boost::shared_ptr< ThreadManager > ThreadManagerSharedPtr