49 map<LibUtilities::ShapeType, DerivUtilSharedPtr> ProcessVarOpti::BuildDerivUtil(
53 map<LibUtilities::ShapeType, DerivUtilSharedPtr> ret;
58 typedef std::pair<LibUtilities::PointsType, LibUtilities::PointsType>
61 map<LibUtilities::ShapeType, PTypes> typeMap;
64 if (m_mesh->m_nummode + o <= 11)
79 for (it = typeMap.begin(); it != typeMap.end(); it++)
81 PTypes pType = it->second;
87 const int pDim = pkey1.GetPointsDim();
88 const int order = m_mesh->m_nummode - 1;
110 der->ptsStd = u1[0].num_elements();
111 der->pts = u2[0].num_elements();
142 "Unknown element type for derivative utility setup");
148 VandermondeI.Invert();
150 for (
int i = 0; i < pDim; ++i)
154 der->VdmD[i] = interp * der->VdmDStd[i];
162 ret[it->first] = der;
176 value_vector(val_vec) {}
180 return value_vector[i1] > value_vector[i2];
186 vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
187 vector<ElUtilSharedPtr> elLock)
192 for (
int i = 0; i < elLock.size(); i++)
194 vector<NodeSharedPtr> nodes;
195 elLock[i]->GetEl()->GetCurvedNodes(nodes);
196 for (
int j = 0; j < nodes.size(); j++)
198 ignoredNodes.insert(nodes[j]);
205 switch (m_mesh->m_spaceDim)
210 for(it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
212 if((*it)->m_elLink.size() == 2)
217 boundaryNodes.insert((*it)->m_n1);
218 boundaryNodes.insert((*it)->m_n2);
219 for(
int i = 0; i < (*it)->m_edgeNodes.size(); i++)
221 boundaryNodes.insert((*it)->m_edgeNodes[i]);
231 for (it = m_mesh->m_faceSet.begin();
232 it != m_mesh->m_faceSet.end(); it++)
234 if ((*it)->m_elLink.size() == 2)
239 vector<NodeSharedPtr> vs = (*it)->m_vertexList;
240 for (
int j = 0; j < vs.size(); j++)
242 boundaryNodes.insert(vs[j]);
245 vector<EdgeSharedPtr> es = (*it)->m_edgeList;
246 for (
int j = 0; j < es.size(); j++)
248 for (
int k = 0; k < es[j]->m_edgeNodes.size(); k++)
250 boundaryNodes.insert(es[j]->m_edgeNodes[k]);
254 for (
int i = 0; i < (*it)->m_faceNodes.size(); i++)
256 boundaryNodes.insert((*it)->m_faceNodes[i]);
264 for (nit = m_mesh->m_vertexSet.begin(); nit != m_mesh->m_vertexSet.end(); ++nit)
266 if((*nit)->GetNumCadCurve() > 1)
268 boundaryNodes.insert((*nit));
280 vector<NodeSharedPtr> remainEdgeVertex;
281 vector<NodeSharedPtr> remainFace;
282 vector<NodeSharedPtr> remainVolume;
287 for (nit = m_mesh->m_vertexSet.begin(); nit != m_mesh->m_vertexSet.end();
292 if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
294 remainEdgeVertex.push_back(*nit);
295 if ((*nit)->GetNumCadCurve() == 1)
299 else if ((*nit)->GetNumCADSurf() == 1)
305 m_res->nDoF += m_mesh->m_spaceDim;
312 for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
314 vector<NodeSharedPtr> n = (*eit)->m_edgeNodes;
315 for (
int j = 0; j < n.size(); j++)
319 if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
321 remainEdgeVertex.push_back(n[j]);
322 if (n[j]->GetNumCadCurve() == 1)
326 else if (n[j]->GetNumCADSurf() == 1)
332 m_res->nDoF += m_mesh->m_spaceDim;
340 for (fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
342 for (
int j = 0; j < (*fit)->m_faceNodes.size(); j++)
346 if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
348 remainFace.push_back((*fit)->m_faceNodes[j]);
349 if ((*fit)->m_faceNodes[j]->GetNumCADSurf() == 1)
355 m_res->nDoF += m_mesh->m_spaceDim;
362 for (
int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
364 vector<NodeSharedPtr> ns =
365 m_mesh->m_element[m_mesh->m_expDim][i]->GetVolumeNodes();
366 for (
int j = 0; j < ns.size(); j++)
370 if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
372 remainVolume.push_back(ns[j]);
373 m_res->nDoF += m_mesh->m_spaceDim;
379 m_res->n = remainEdgeVertex.size()
380 + remainFace.size() + remainVolume.size();
383 vector<vector<NodeSharedPtr> > ret;
384 vector<vector<NodeSharedPtr> > retPart;
389 vector<int> num_el(remainEdgeVertex.size());
390 for (
int i = 0; i < remainEdgeVertex.size(); i++)
394 vector<ElUtilSharedPtr> &elUtils = it->second;
395 num_el[i] = elUtils.size();
398 vector<int> permNode(remainEdgeVertex.size());
399 for (
int i = 0; i < remainEdgeVertex.size(); ++i)
403 std::sort(permNode.begin(), permNode.end(),
NodeComparator(num_el));
405 vector<NodeSharedPtr> remainEdgeVertexSort(remainEdgeVertex.size());
406 for (
int i = 0; i < remainEdgeVertex.size(); ++i)
409 remainEdgeVertexSort[i] = remainEdgeVertex[j];
412 retPart = CreateColoursets(remainEdgeVertexSort);
413 if(m_mesh->m_verbose)
415 cout <<
"Number of Edge/Vertex Coloursets: " << retPart.size() << endl;
417 for (
int i = 0; i < retPart.size(); i++)
419 ret.push_back(retPart[i]);
424 retPart = CreateColoursets(remainFace);
425 if(m_mesh->m_verbose)
427 cout <<
"Number of Face Coloursets: " << retPart.size() << endl;
429 for (
int i = 0; i < retPart.size(); i++)
431 ret.push_back(retPart[i]);
436 retPart = CreateColoursets(remainVolume);
437 if(m_mesh->m_verbose)
439 cout <<
"Number of Volume Coloursets: " << retPart.size() << endl;
441 for (
int i = 0; i < retPart.size(); i++)
443 ret.push_back(retPart[i]);
447 if(m_mesh->m_verbose)
455 vector<vector<NodeSharedPtr> > ProcessVarOpti::CreateColoursets(
456 vector<NodeSharedPtr> remain)
458 vector<vector<NodeSharedPtr> > retPart;
461 while (remain.size() > 0)
463 vector<NodeSharedPtr> layer;
466 for (
int i = 0; i < remain.size(); i++)
469 ASSERTL0(it != m_nodeElMap.end(),
"could not find node");
471 vector<ElUtilSharedPtr> &elUtils = it->second;
473 bool islocked =
false;
474 for (
int j = 0; j < elUtils.size(); j++)
477 if (sit != locked.end())
485 layer.push_back(remain[i]);
486 completed.insert(remain[i]->m_id);
487 for (
int j = 0; j < elUtils.size(); j++)
489 locked.insert(elUtils[j]->GetId());
495 vector<NodeSharedPtr> tmp = remain;
497 for (
int i = 0; i < tmp.size(); i++)
500 if (sit == completed.end())
502 remain.push_back(tmp[i]);
507 retPart.push_back(layer);
510 if(m_mesh->m_verbose)
520 void ProcessVarOpti::GetElementMap(
521 int o, map<LibUtilities::ShapeType, DerivUtilSharedPtr> derMap)
523 for (
int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
526 vector<NodeSharedPtr> ns;
527 el->GetCurvedNodes(ns);
529 el, derMap[el->GetShapeType()], m_res, m_mesh->m_nummode, o));
530 m_dataSet.push_back(d);
533 for (
int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
536 vector<NodeSharedPtr> ns;
537 el->GetCurvedNodes(ns);
539 for (
int j = 0; j < ns.size(); j++)
541 m_nodeElMap[ns[j]->m_id].push_back(m_dataSet[i]);
544 ASSERTL0(derMap[el->GetShapeType()]->ptsStd == ns.size(),
545 "mismatch node count");
549 vector<ElUtilSharedPtr> ProcessVarOpti::GetLockedElements(
NekDouble thres)
551 vector<ElUtilSharedPtr> elBelowThres;
552 for (
int i = 0; i < m_dataSet.size(); ++i)
554 if (m_dataSet[i]->GetScaledJac() < thres)
556 elBelowThres.push_back(m_dataSet[i]);
560 boost::unordered_set<int> inmesh;
562 vector<ElUtilSharedPtr> totest;
564 for (
int i = 0; i < elBelowThres.size(); i++)
566 t = inmesh.insert(elBelowThres[i]->GetId());
568 vector<FaceSharedPtr> f = elBelowThres[i]->GetEl()->GetFaceList();
569 for (
int j = 0; j < f.size(); j++)
571 for (
int k = 0; k < f[j]->m_elLink.size(); k++)
573 if (f[j]->m_elLink[k].first->GetId() ==
574 elBelowThres[i]->GetId())
577 t = inmesh.insert(f[j]->m_elLink[k].first->GetId());
581 m_dataSet[f[j]->m_elLink[k].first->GetId()]);
587 for (
int i = 0; i < 6; i++)
589 vector<ElUtilSharedPtr> tmp = totest;
591 for (
int j = 0; j < tmp.size(); j++)
593 vector<FaceSharedPtr> f = tmp[j]->GetEl()->GetFaceList();
594 for (
int k = 0; k < f.size(); k++)
596 for (
int l = 0; l < f[k]->m_elLink.size(); l++)
598 if (f[k]->m_elLink[l].first->GetId() == tmp[j]->GetId())
601 t = inmesh.insert(f[k]->m_elLink[l].first->GetId());
605 m_dataSet[f[k]->m_elLink[l].first->GetId()]);
613 vector<ElUtilSharedPtr> ret;
614 for (
int i = 0; i < m_dataSet.size(); ++i)
617 inmesh.find(m_dataSet[i]->GetId());
618 if (s == inmesh.end())
620 ret.push_back(m_dataSet[i]);
627 void ProcessVarOpti::RemoveLinearCurvature()
629 for(
int i = 0; i < m_dataSet.size(); i++)
631 if(m_dataSet[i]->GetScaledJac() > 0.999)
634 vector<NodeSharedPtr> ns;
635 el->SetVolumeNodes(ns);
639 map<int, vector<FaceSharedPtr> > edgeToFace;
642 for(fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
645 for(
int i = 0; i < (*fit)->m_elLink.size(); i++)
647 int id = (*fit)->m_elLink[i].first->GetId();
648 if(m_dataSet[
id]->GetScaledJac() <= 0.999)
656 (*fit)->m_faceNodes.clear();
659 vector<EdgeSharedPtr> es = (*fit)->m_edgeList;
660 for(
int i = 0; i < es.size(); i++)
662 edgeToFace[es[i]->m_id].push_back(*fit);
667 for(eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
669 map<int, vector<FaceSharedPtr> >
::iterator it;
670 it = edgeToFace.find((*eit)->m_id);
671 ASSERTL0(it != edgeToFace.end(),
"not found");
673 for(
int i = 0; i < it->second.size(); i++)
675 if(it->second[i]->m_faceNodes.size() > 0)
683 (*eit)->m_edgeNodes.clear();
#define ASSERTL0(condition, msg)
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
A class to assist in the construction of nodal simplex and hybrid elements in two and three dimension...
Specialisation of the NodalUtil class to support nodal triangular elements.
boost::unordered_set< NodeSharedPtr, NodeHash > NodeSet
3D Nodal Symmetric positive internal tet (Whitherden, Vincent)
SharedMatrix GetVandermondeForDeriv(int dir)
Return the Vandermonde matrix of the derivative of the basis functions for the nodal distribution...
SharedMatrix GetInterpolationMatrix(Array< OneD, Array< OneD, NekDouble > > &xi)
Construct the interpolation matrix used to evaluate the basis at the points xi inside the element...
const vector< int > & value_vector
Specialisation of the NodalUtil class to support nodal hex elements.
3D Nodal Electrostatic Points on a Tetrahedron
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
Specialisation of the NodalUtil class to support nodal prismatic elements.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
SharedMatrix GetVandermonde()
Return the Vandermonde matrix for the nodal distribution.
NodeComparator(const vector< int > &val_vec)
boost::shared_ptr< DerivUtil > DerivUtilSharedPtr
Specialisation of the NodalUtil class to support nodal quad elements.
boost::shared_ptr< Element > ElementSharedPtr
bool operator()(int i1, int i2)
boost::shared_ptr< ElUtil > ElUtilSharedPtr
2D Nodal Electrostatic Points on a Triangle
3D electrostatically spaced points on a Prism
Specialisation of the NodalUtil class to support nodal tetrahedral elements.
2D Nodal Symmetric positive internal triangle (Whitherden, Vincent)