36 #include <boost/geometry.hpp>
37 #include <boost/geometry/geometries/point.hpp>
38 #include <boost/geometry/geometries/box.hpp>
39 #include <boost/geometry/index/rtree.hpp>
54 namespace bg = boost::geometry;
55 namespace bgi = boost::geometry::index;
60 #define TOL_BLEND 1.0e-8
102 ConfigOption(
false,
"5",
"Number of points to add to face edges.");
104 ConfigOption(
false,
"-1",
"Tag identifying surface to process.");
106 true,
"-1",
"Curve both triangular faces of prism on boundary.");
108 false,
"NoFile",
"Use alternative file for Spherigon definition");
110 false,
"1.0",
"Apply scaling factor to coordinates in file ");
114 "Add randowm noise to normals of amplitude AMP "
115 "in specified region. input string is "
116 "Amp,xmin,xmax,ymin,ymax,zmin,zmax");
154 return sqrt(tmp1 * tmp1 + tmp2 * tmp2 + tmp3 * tmp3);
184 vector<double> &blend)
186 vector<double> tmp(r.size());
187 double totBlend = 0.0;
191 for (i = 0; i < nV; ++i)
194 tmp[i] = (Q[i] -
P).abs2();
197 for (i = 0; i < nV; ++i)
199 int ip = (i + 1) % nV, im = (i - 1 + nV) % nV;
204 r[i] * r[i] * (r[im] * r[im] * tmp[im] / (tmp[im] + tmp[i]) +
205 r[ip] * r[ip] * tmp[ip] / (tmp[ip] + tmp[i]));
206 totBlend += blend[i];
210 for (i = 0; i < nV; ++i)
212 blend[i] /= totBlend;
225 map<int,NodeSharedPtr> &surfverts)
233 typedef bg::model::point<NekDouble, 3, bg::cs::cartesian>
Point;
234 typedef pair<Point, unsigned int> PointI;
238 map<int,int> TreeidtoPlyid;
241 vector<PointI> dataPts;
242 for (j = 0, it = plymesh->m_vertexSet.begin();
243 it != plymesh->m_vertexSet.end();
246 dataPts.push_back(make_pair(Point( (*it)->m_x,
249 TreeidtoPlyid[j] = (*it)->m_id;
253 bgi::rtree<PointI, bgi::rstar<16> > rtree;
254 rtree.insert(dataPts.begin(), dataPts.end());
257 for (cnt = 0, vIt = surfverts.begin(); vIt != surfverts.end();
263 "Nearest ply verts",prog);
270 Point queryPt(vIt->second->m_x, vIt->second->m_y, vIt->second->m_z);
272 vector<PointI> result;
273 rtree.query(bgi::nearest(queryPt, n_neighbs), std::back_inserter(result));
275 ASSERTL1(bg::distance(result[0].first,queryPt) < bg::distance(result[1].first,queryPt),
276 "Assumption that dist values are ordered from smallest to largest is not correct");
278 cntmin = TreeidtoPlyid[result[0].second];
280 ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
281 "cntmin is out of range");
283 m_mesh->m_vertexNormals[vIt->first] =
284 plymesh->m_vertexNormals[cntmin];
306 for (
int i = 0; i < el.size(); ++i)
314 "Spherigon expansions must be lines, triangles or "
318 int nV = e->GetVertexCount();
319 vector<NodeSharedPtr> node(nV);
321 for (
int j = 0; j < nV; ++j)
323 node[j] = e->GetVertex(j);
328 if (mesh->m_spaceDim == 3)
331 Node v1 = *(node[1]) - *(node[0]);
332 Node v2 = *(node[2]) - *(node[0]);
338 Node dx = *(node[1]) - *(node[0]);
348 for (
int j = 0; j < nV; ++j)
350 nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
351 if (nIt == mesh->m_vertexNormals.end())
353 mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
363 for (nIt = mesh->m_vertexNormals.begin();
364 nIt != mesh->m_vertexNormals.end();
367 Node &n = mesh->m_vertexNormals[nIt->first];
379 "Spherigon implementation only valid in 2D/3D.");
382 boost::unordered_set<int> visitedEdges;
385 vector<ElementSharedPtr> el;
389 cout <<
"ProcessSpherigon: Smoothing mesh..." << endl;
406 string surfTag =
m_config[
"surf"].as<
string>();
407 bool prismTag =
m_config[
"BothTriFacesOnPrism"].beenSet;
411 vector<unsigned int> surfs;
413 sort(surfs.begin(), surfs.end());
415 m_mesh->m_spherigonSurfs.clear();
416 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
419 int nSurf =
m_mesh->m_expDim == 3 ? el->GetFaceCount()
420 : el->GetEdgeCount();
422 for (
int j = 0; j < nSurf; ++j)
424 int bl = el->GetBoundaryLink(j);
432 vector<int> tags = bEl->GetTagList();
435 sort(tags.begin(), tags.end());
436 set_intersection(surfs.begin(),
440 back_inserter(inter));
442 if (inter.size() == 1)
444 m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
448 if (nSurf == 5 && prismTag)
452 int triFace = j == 1 ? 3 : 1;
453 m_mesh->m_spherigonSurfs.insert(
454 make_pair(i, triFace));
461 if (
m_mesh->m_spherigonSurfs.size() == 0)
463 cerr <<
"WARNING: Spherigon surfaces have not been defined "
464 <<
"-- ignoring smoothing." << endl;
468 if (
m_mesh->m_expDim == 3)
470 for (it =
m_mesh->m_spherigonSurfs.begin();
471 it !=
m_mesh->m_spherigonSurfs.end();
475 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetFace(
477 vector<NodeSharedPtr> nodes = f->m_vertexList;
487 for (
int i = 0; i < f->m_vertexList.size(); ++i)
489 elmt->SetVertex(i, f->m_vertexList[i]);
491 for (
int i = 0; i < f->m_edgeList.size(); ++i)
493 elmt->SetEdge(i, f->m_edgeList[i]);
501 for (it =
m_mesh->m_spherigonSurfs.begin();
502 it !=
m_mesh->m_spherigonSurfs.end();
506 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetEdge(
508 vector<NodeSharedPtr> nodes;
512 nodes.push_back(edge->m_n1);
513 nodes.push_back(edge->m_n2);
520 elmt->SetVertex(0, nodes[0]);
521 elmt->SetVertex(1, nodes[1]);
524 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetEdge(
532 ASSERTL0(
false,
"Spherigon expansions must be 2/3 dimensional");
537 bool normalsGenerated =
false;
540 std::string normalfile =
m_config[
"usenormalfile"].as<
string>();
541 if (normalfile.compare(
"NoFile") != 0)
547 cout <<
"Inputing normal file: " << normalfile
548 <<
" with scaling of " << scale << endl;
552 io::filtering_istream inply;
555 inplyTmp.open(normalfile.c_str());
557 string(
"Could not open input ply file: ") + normalfile);
559 inply.push(inplyTmp);
562 plyfile = boost::shared_ptr<InputPly>(
new InputPly(m));
563 plyfile->ReadPly(inply, scale);
564 plyfile->ProcessVertices();
569 cout <<
"\t Generating ply normals" << endl;
575 Node minx(0, 0.0, 0.0, 0.0), tmp, tmpsav;
578 map<int, NodeSharedPtr> surfverts;
581 for (
int i = 0; i < el.size(); ++i)
584 int nV = e->GetVertexCount();
585 for (
int j = 0; j < nV; ++j)
587 int id = e->GetVertex(j)->m_id;
588 surfverts[id] = e->GetVertex(j);
595 cout <<
"\t Processing surface normals " << endl;
602 normalsGenerated =
true;
604 else if (
m_mesh->m_vertexNormals.size() == 0)
607 normalsGenerated =
true;
611 std::string normalnoise =
m_config[
"normalnoise"].as<
string>();
612 if (normalnoise.compare(
"NotSpecified") != 0)
614 vector<NekDouble> values;
617 "Failed to interpret normal noise string");
619 int nvalues = values.size() / 2;
624 cout <<
"\t adding noise to normals of amplitude " << amp
626 for (
int i = 0; i < nvalues; ++i)
628 cout << values[2 * i + 1] <<
"," << values[2 * i + 2] <<
" ";
634 map<int, NodeSharedPtr> surfverts;
637 for (
int i = 0; i < el.size(); ++i)
640 int nV = e->GetVertexCount();
641 for (
int j = 0; j < nV; ++j)
643 int id = e->GetVertex(j)->m_id;
644 surfverts[id] = e->GetVertex(j);
648 for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
650 bool AddNoise =
false;
652 for (
int i = 0; i < nvalues; ++i)
659 if (((vIt->second)->m_x > values[2 * i + 1]) &&
660 ((vIt->second)->m_x < values[2 * i + 2]))
668 if (((vIt->second)->m_x > values[2 * i + 1]) &&
669 ((vIt->second)->m_x < values[2 * i + 2]) &&
670 ((vIt->second)->m_y > values[2 * i + 3]) &&
671 ((vIt->second)->m_y < values[2 * i + 4]))
679 if (((vIt->second)->m_x > values[2 * i + 1]) &&
680 ((vIt->second)->m_x < values[2 * i + 2]) &&
681 ((vIt->second)->m_y > values[2 * i + 3]) &&
682 ((vIt->second)->m_y < values[2 * i + 4]) &&
683 ((vIt->second)->m_z > values[2 * i + 5]) &&
684 ((vIt->second)->m_z < values[2 * i + 6]))
696 Node rvec(0, rand(), rand(), rand());
697 rvec *= values[0] / sqrt(rvec.
abs2());
699 Node normal =
m_mesh->m_vertexNormals[vIt->first];
702 normal /= sqrt(normal.abs2());
704 m_mesh->m_vertexNormals[vIt->first] = normal;
712 int nquad =
m_mesh->m_spaceDim == 3 ? nq * nq : nq;
721 ASSERTL0(nq > 2,
"Number of points must be greater than 2.");
736 stdtri->GetNodalPoints(xnodal, ynodal);
738 int edgeMap[3][4][2] = {
739 {{0, 1}, {-1, -1}, {-1, -1}, {-1, -1}},
740 {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}},
741 {{0, 1}, {nq - 1, nq}, {nq * nq - 1, -1}, {nq * (nq - 1), -nq}}
744 int vertMap[3][4][2] = {
745 {{0, 1}, {0, 0}, {0, 0}, {0, 0}},
746 {{0, 1}, {1, 2}, {2, 3}, {0, 0}},
747 {{0, 1}, {1, 2}, {2, 3}, {3, 0}},
750 for (
int i = 0; i < el.size(); ++i)
766 e->GetGeom(
m_mesh->m_spaceDim));
770 seg->GetCoords(x, y, z);
783 tri->GetCoords(xc, yc, zc);
784 nquad = nq * (nq + 1) / 2;
786 for (
int j = 0; j < nquad; ++j)
788 coord[0] = xnodal[j];
789 coord[1] = ynodal[j];
790 x[j] = stdtri->PhysEvaluate(coord, xc);
793 for (
int j = 0; j < nquad; ++j)
795 coord[0] = xnodal[j];
796 coord[1] = ynodal[j];
797 y[j] = stdtri->PhysEvaluate(coord, yc);
800 for (
int j = 0; j < nquad; ++j)
802 coord[0] = xnodal[j];
803 coord[1] = ynodal[j];
804 z[j] = stdtri->PhysEvaluate(coord, zc);
815 quad->GetCoords(x, y, z);
820 ASSERTL0(
false,
"Unknown expansion type.");
824 if (
m_mesh->m_spaceDim == 2)
830 int nV = e->GetVertexCount();
832 for (
int j = 0; j < nV; ++j)
834 v.push_back(*(e->GetVertex(j)));
836 "Normal has not been defined");
837 vN.push_back(
m_mesh->m_vertexNormals[v[j].m_id]);
840 vector<Node> tmp(nV);
841 vector<NekDouble> r(nV);
845 vector<NekDouble> blend(nV);
846 vector<Node> out(nquad);
849 NekDouble segLength = sqrt((v[0] - v[1]).abs2());
852 for (
int j = 0; j < nquad; ++j)
854 Node P(0, x[j], y[j], z[j]);
859 if (
m_mesh->m_spaceDim == 2)
864 r[0] = sqrt((P - v[0]).abs2()) / segLength;
865 r[0] = max(min(1.0, r[0]), 0.0);
869 N = vN[0] * r[0] + vN[1] * r[1];
871 else if (
m_mesh->m_spaceDim == 3)
873 for (
int k = 0; k < nV; ++k)
881 for (
int k = 0; k < nV; ++k)
884 for (
int l = 0; l < nV - 2; ++l)
887 tmp[(k + l + 2) % nV]);
893 for (
int k = 0; k < nV; ++k)
903 for (
int k = 0; k < nV; ++k)
908 K[k] = P + N * ((v[k] -
P).dot(N));
909 tmp1 = (v[k] - K[k]).dot(vN[k]) / (1.0 + N.
dot(vN[k]));
910 Q[k] = K[k] + N * tmp1;
911 Qp[k] = v[k] - N * ((v[k] -
P).dot(N));
920 for (
int k = 0; k < nV; ++k)
922 P += Q[k] * blend[k];
925 if ((boost::math::isnan)(P.
m_x) || (boost::math::isnan)(P.
m_y) ||
926 (boost::math::isnan)(P.
m_z))
929 "spherigon point is a nan. Check to see if "
930 "ply file is correct if using input normal file");
940 int offset = (int)e->GetConf().m_e - 2;
942 for (
int edge = 0; edge < e->GetEdgeCount(); ++edge)
944 eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
945 if (eIt == visitedEdges.end())
948 !(v[vertMap[offset][edge][0]] == *(e->GetEdge(edge)->m_n1));
951 e->GetEdge(edge)->m_edgeNodes.clear();
955 for (
int j = 1; j < nq - 1; ++j)
957 int v = edgeMap[offset][edge][0] +
958 j * edgeMap[offset][edge][1];
959 e->GetEdge(edge)->m_edgeNodes.push_back(
965 for (
int j = 0; j < nq - 2; ++j)
967 int v = 3 + edge * (nq - 2) + j;
968 e->GetEdge(edge)->m_edgeNodes.push_back(
975 reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
976 e->GetEdge(edge)->m_edgeNodes.end());
979 e->GetEdge(edge)->m_curveType =
982 visitedEdges.insert(e->GetEdge(edge)->m_id);
987 if (
m_mesh->m_spaceDim == 3)
989 vector<NodeSharedPtr> volNodes;
993 volNodes.resize((nq - 2) * (nq - 2));
994 for (
int j = 1; j < nq - 1; ++j)
996 for (
int k = 1; k < nq - 1; ++k)
999 volNodes[(j - 1) * (nq - 2) + (k - 1)] =
1006 for (
int j = 3 + 3 * (nq - 2); j < nquad; ++j)
1012 e->SetVolumeNodes(volNodes);
1017 if (
m_mesh->m_expDim == 3)
1021 for (it =
m_mesh->m_spherigonSurfs.begin();
1022 it !=
m_mesh->m_spherigonSurfs.end();
1026 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetFace(
1029 f->m_faceNodes = el[elmt]->GetVolumeNodes();
1030 f->m_curveType = f->m_vertexList.size() == 3
1036 if (normalsGenerated)
1038 m_mesh->m_vertexNormals.clear();
void UnitCrossProd(NekMeshUtils::Node &a, NekMeshUtils::Node &b, NekMeshUtils::Node &c)
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
NekDouble CrossProdMag(NekMeshUtils::Node &a, NekMeshUtils::Node &b)
Calculate the magnitude of the cross product .
#define ASSERTL0(condition, msg)
Basic information about an element.
void FindNormalFromPlyFile(NekMeshUtils::MeshSharedPtr &plymesh, std::map< int, NekMeshUtils::NodeSharedPtr > &surfverts)
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.
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual ~ProcessSpherigon()
Destructor.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
void SuperBlend(std::vector< NekDouble > &r, std::vector< NekMeshUtils::Node > &Q, NekMeshUtils::Node &P, std::vector< NekDouble > &blend)
Calculate the blending function for spherigon implementation.
void GenerateNormals(std::vector< NekMeshUtils::ElementSharedPtr > &el, NekMeshUtils::MeshSharedPtr &mesh)
Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hy...
Principle Modified Functions .
pair< ModuleType, string > ModuleKey
NekDouble m_y
Y-coordinate.
MeshSharedPtr m_mesh
Mesh object.
ElementFactory & GetElementFactory()
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
boost::shared_ptr< StdNodalTriExp > StdNodalTriExpSharedPtr
Represents a point in the domain.
virtual void Process()
Write mesh to output file.
Gauss Radau pinned at x=-1, .
Principle Orthogonal Functions .
NEKMESHUTILS_EXPORT NekDouble abs2() const
boost::shared_ptr< SegExp > SegExpSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Principle Orthogonal Functions .
Defines a specification for a set of points.
boost::shared_ptr< QuadExp > QuadExpSharedPtr
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< InputPly > InputPlySharedPtr
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
NekDouble m_x
X-coordinate.
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Abstract base class for processing modules.
NekDouble Blend(NekDouble r)
Calculate the blending function for spherigon implementation.
boost::shared_ptr< TriGeom > TriGeomSharedPtr
boost::shared_ptr< Element > ElementSharedPtr
boost::shared_ptr< Face > FaceSharedPtr
NekDouble m_z
Z-coordinate.
std::pair< ModuleType, std::string > ModuleKey
void Zero(int n, T *x, const int incx)
Zero vector.
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Describes the specification for a Basis.
1D Gauss-Lobatto-Legendre quadrature points
2D Nodal Electrostatic Points on a Triangle
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
boost::shared_ptr< NodalTriExp > NodalTriExpSharedPtr