51 #define TOL_BLEND 1.0e-8
94 ConfigOption(
false,
"5",
"Number of points to add to face edges.");
96 ConfigOption(
false,
"-1",
"Tag identifying surface to process.");
98 true,
"-1",
"Curve both triangular faces of prism on boundary.");
100 false,
"NoFile",
"Use alternative file for Spherigon definition");
102 false,
"1.0",
"Apply scaling factor to coordinates in file ");
106 "Add randowm noise to normals of amplitude AMP "
107 "in specified region. input string is "
108 "Amp,xmin,xmax,ymin,ymax,zmin,zmax");
146 return sqrt(tmp1 * tmp1 + tmp2 * tmp2 + tmp3 * tmp3);
176 vector<double> &blend)
178 vector<double> tmp(r.size());
179 double totBlend = 0.0;
183 for (i = 0; i < nV; ++i)
186 tmp[i] = (Q[i] - P).abs2();
189 for (i = 0; i < nV; ++i)
191 int ip = (i + 1) % nV, im = (i - 1 + nV) % nV;
196 r[i] * r[i] * (r[im] * r[im] * tmp[im] / (tmp[im] + tmp[i]) +
197 r[ip] * r[ip] * tmp[ip] / (tmp[ip] + tmp[i]));
198 totBlend += blend[i];
202 for (i = 0; i < nV; ++i)
204 blend[i] /= totBlend;
234 for (
int i = 0; i < el.size(); ++i)
242 "Spherigon expansions must be lines, triangles or "
246 int nV = e->GetVertexCount();
247 vector<NodeSharedPtr> node(nV);
249 for (
int j = 0; j < nV; ++j)
251 node[j] = e->GetVertex(j);
256 if (mesh->m_spaceDim == 3)
259 Node v1 = *(node[1]) - *(node[0]);
260 Node v2 = *(node[2]) - *(node[0]);
266 Node dx = *(node[1]) - *(node[0]);
267 dx /= sqrt(dx.
abs2());
275 for (
int j = 0; j < nV; ++j)
277 nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
278 if (nIt == mesh->m_vertexNormals.end())
280 mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
290 for (nIt = mesh->m_vertexNormals.begin();
291 nIt != mesh->m_vertexNormals.end();
294 Node &n = mesh->m_vertexNormals[nIt->first];
305 "Spherigon implementation only valid in 2D/3D.");
308 boost::unordered_set<int> visitedEdges;
311 vector<ElementSharedPtr> el;
315 cout <<
"ProcessSpherigon: Smoothing mesh..." << endl;
332 string surfTag =
m_config[
"surf"].as<
string>();
333 bool prismTag =
m_config[
"BothTriFacesOnPrism"].beenSet;
337 vector<unsigned int> surfs;
339 sort(surfs.begin(), surfs.end());
341 m_mesh->m_spherigonSurfs.clear();
342 for (
int i = 0; i <
m_mesh->m_element[
m_mesh->m_expDim].size(); ++i)
345 int nSurf =
m_mesh->m_expDim == 3 ? el->GetFaceCount()
346 : el->GetEdgeCount();
348 for (
int j = 0; j < nSurf; ++j)
350 int bl = el->GetBoundaryLink(j);
358 vector<int> tags = bEl->GetTagList();
361 sort(tags.begin(), tags.end());
362 set_intersection(surfs.begin(),
366 back_inserter(inter));
368 if (inter.size() == 1)
370 m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
374 if (nSurf == 5 && prismTag)
378 int triFace = j == 1 ? 3 : 1;
379 m_mesh->m_spherigonSurfs.insert(
380 make_pair(i, triFace));
387 if (
m_mesh->m_spherigonSurfs.size() == 0)
389 cerr <<
"WARNING: Spherigon surfaces have not been defined "
390 <<
"-- ignoring smoothing." << endl;
394 if (
m_mesh->m_expDim == 3)
396 for (it =
m_mesh->m_spherigonSurfs.begin();
397 it !=
m_mesh->m_spherigonSurfs.end();
401 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetFace(
403 vector<NodeSharedPtr> nodes = f->m_vertexList;
413 for (
int i = 0; i < f->m_vertexList.size(); ++i)
415 elmt->SetVertex(i, f->m_vertexList[i]);
417 for (
int i = 0; i < f->m_edgeList.size(); ++i)
419 elmt->SetEdge(i, f->m_edgeList[i]);
427 for (it =
m_mesh->m_spherigonSurfs.begin();
428 it !=
m_mesh->m_spherigonSurfs.end();
432 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetEdge(
434 vector<NodeSharedPtr> nodes;
438 nodes.push_back(edge->m_n1);
439 nodes.push_back(edge->m_n2);
446 elmt->SetVertex(0, nodes[0]);
447 elmt->SetVertex(1, nodes[1]);
450 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetEdge(
458 ASSERTL0(
false,
"Spherigon expansions must be 2/3 dimensional");
463 bool normalsGenerated =
false;
466 std::string normalfile =
m_config[
"usenormalfile"].as<
string>();
467 if (normalfile.compare(
"NoFile") != 0)
473 cout <<
"Inputing normal file: " << normalfile
474 <<
" with scaling of " << scale << endl;
480 inply.open(normalfile.c_str());
481 ASSERTL0(inply,
string(
"Could not open input ply file: ") + normalfile);
485 plyfile = boost::shared_ptr<InputPly>(
new InputPly(m));
486 plyfile->ReadPly(inply, scale);
487 plyfile->ProcessVertices();
495 Node minx(0, 0.0, 0.0, 0.0), tmp, tmpsav;
498 map<int, NodeSharedPtr> surfverts;
501 for (
int i = 0; i < el.size(); ++i)
504 int nV = e->GetVertexCount();
505 for (
int j = 0; j < nV; ++j)
507 int id = e->GetVertex(j)->m_id;
508 surfverts[id] = e->GetVertex(j);
520 cout <<
"\t Processing surface normals " << endl;
523 map<int, int> locnorm;
524 for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt, ++cnt)
528 for (j = 0, it = plymesh->m_vertexSet.begin();
529 it != plymesh->m_vertexSet.end();
532 tmp = *(vIt->second) - *(*it);
538 cntmin = (*it)->m_id;
542 locnorm[cntmin] = vIt->first;
544 ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),
545 "cntmin is out of range");
546 m_mesh->m_vertexNormals[vIt->first] =
547 plymesh->m_vertexNormals[cntmin];
551 cout <<
"\t end of processing surface normals " << endl;
553 normalsGenerated =
true;
555 else if (
m_mesh->m_vertexNormals.size() == 0)
558 normalsGenerated =
true;
562 std::string normalnoise =
m_config[
"normalnoise"].as<
string>();
563 if (normalnoise.compare(
"NotSpecified") != 0)
565 vector<NekDouble> values;
568 "Failed to interpret normal noise string");
570 int nvalues = values.size() / 2;
575 cout <<
"\t adding noise to normals of amplitude " << amp
577 for (
int i = 0; i < nvalues; ++i)
579 cout << values[2 * i + 1] <<
"," << values[2 * i + 2] <<
" ";
585 map<int, NodeSharedPtr> surfverts;
588 for (
int i = 0; i < el.size(); ++i)
591 int nV = e->GetVertexCount();
592 for (
int j = 0; j < nV; ++j)
594 int id = e->GetVertex(j)->m_id;
595 surfverts[id] = e->GetVertex(j);
599 for (vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
601 bool AddNoise =
false;
603 for (
int i = 0; i < nvalues; ++i)
610 if (((vIt->second)->m_x > values[2 * i + 1]) &&
611 ((vIt->second)->m_x < values[2 * i + 2]))
619 if (((vIt->second)->m_x > values[2 * i + 1]) &&
620 ((vIt->second)->m_x < values[2 * i + 2]) &&
621 ((vIt->second)->m_y > values[2 * i + 3]) &&
622 ((vIt->second)->m_y < values[2 * i + 4]))
630 if (((vIt->second)->m_x > values[2 * i + 1]) &&
631 ((vIt->second)->m_x < values[2 * i + 2]) &&
632 ((vIt->second)->m_y > values[2 * i + 3]) &&
633 ((vIt->second)->m_y < values[2 * i + 4]) &&
634 ((vIt->second)->m_z > values[2 * i + 5]) &&
635 ((vIt->second)->m_z < values[2 * i + 6]))
647 Node rvec(0, rand(), rand(), rand());
648 rvec *= values[0] / sqrt(rvec.
abs2());
650 Node normal =
m_mesh->m_vertexNormals[vIt->first];
653 normal /= sqrt(normal.abs2());
655 m_mesh->m_vertexNormals[vIt->first] = normal;
663 int nquad =
m_mesh->m_spaceDim == 3 ? nq * nq : nq;
672 ASSERTL0(nq > 2,
"Number of points must be greater than 2.");
687 stdtri->GetNodalPoints(xnodal, ynodal);
689 int edgeMap[3][4][2] = {
690 {{0, 1}, {-1, -1}, {-1, -1}, {-1, -1}},
691 {{0, 1}, {nq - 1, nq}, {nq * (nq - 1), -nq}, {-1, -1}},
692 {{0, 1}, {nq - 1, nq}, {nq * nq - 1, -1}, {nq * (nq - 1), -nq}}
695 int vertMap[3][4][2] = {
696 {{0, 1}, {0, 0}, {0, 0}, {0, 0}},
697 {{0, 1}, {1, 2}, {2, 3}, {0, 0}},
698 {{0, 1}, {1, 2}, {2, 3}, {3, 0}},
701 for (
int i = 0; i < el.size(); ++i)
717 e->GetGeom(
m_mesh->m_spaceDim));
721 seg->GetCoords(x, y, z);
734 tri->GetCoords(xc, yc, zc);
735 nquad = nq * (nq + 1) / 2;
737 for (
int j = 0; j < nquad; ++j)
739 coord[0] = xnodal[j];
740 coord[1] = ynodal[j];
741 x[j] = stdtri->PhysEvaluate(coord, xc);
744 for (
int j = 0; j < nquad; ++j)
746 coord[0] = xnodal[j];
747 coord[1] = ynodal[j];
748 y[j] = stdtri->PhysEvaluate(coord, yc);
751 for (
int j = 0; j < nquad; ++j)
753 coord[0] = xnodal[j];
754 coord[1] = ynodal[j];
755 z[j] = stdtri->PhysEvaluate(coord, zc);
766 quad->GetCoords(x, y, z);
771 ASSERTL0(
false,
"Unknown expansion type.");
775 if (
m_mesh->m_spaceDim == 2)
781 int nV = e->GetVertexCount();
783 for (
int j = 0; j < nV; ++j)
785 v.push_back(*(e->GetVertex(j)));
787 "Normal has not been defined");
788 vN.push_back(
m_mesh->m_vertexNormals[v[j].m_id]);
791 vector<Node> tmp(nV);
792 vector<double> r(nV);
796 vector<double> blend(nV);
797 vector<Node> out(nquad);
800 double segLength = sqrt((v[0] - v[1]).abs2());
803 for (
int j = 0; j < nquad; ++j)
805 Node P(0, x[j], y[j], z[j]);
810 if (
m_mesh->m_spaceDim == 2)
815 r[0] = sqrt((P - v[0]).abs2()) / segLength;
816 r[0] = max(min(1.0, r[0]), 0.0);
820 N = vN[0] * r[0] + vN[1] * r[1];
822 else if (
m_mesh->m_spaceDim == 3)
824 for (
int k = 0; k < nV; ++k)
832 for (
int k = 0; k < nV; ++k)
835 for (
int l = 0; l < nV - 2; ++l)
838 tmp[(k + l + 2) % nV]);
844 for (
int k = 0; k < nV; ++k)
854 for (
int k = 0; k < nV; ++k)
859 K[k] = P + N * ((v[k] - P).dot(N));
860 tmp1 = (v[k] - K[k]).dot(vN[k]) / (1.0 + N.
dot(vN[k]));
861 Q[k] = K[k] + N * tmp1;
862 Qp[k] = v[k] - N * ((v[k] - P).dot(N));
871 for (
int k = 0; k < nV; ++k)
873 P += Q[k] * blend[k];
876 if ((boost::math::isnan)(P.
m_x) || (boost::math::isnan)(P.
m_y) ||
877 (boost::math::isnan)(P.
m_z))
880 "spherigon point is a nan. Check to see if "
881 "ply file is correct if using input normal file");
891 int offset = (int)e->GetConf().m_e - 2;
893 for (
int edge = 0; edge < e->GetEdgeCount(); ++edge)
895 eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
896 if (eIt == visitedEdges.end())
899 !(v[vertMap[offset][edge][0]] == *(e->GetEdge(edge)->m_n1));
902 e->GetEdge(edge)->m_edgeNodes.clear();
906 for (
int j = 1; j < nq - 1; ++j)
908 int v = edgeMap[offset][edge][0] +
909 j * edgeMap[offset][edge][1];
910 e->GetEdge(edge)->m_edgeNodes.push_back(
916 for (
int j = 0; j < nq - 2; ++j)
918 int v = 3 + edge * (nq - 2) + j;
919 e->GetEdge(edge)->m_edgeNodes.push_back(
926 reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
927 e->GetEdge(edge)->m_edgeNodes.end());
930 e->GetEdge(edge)->m_curveType =
933 visitedEdges.insert(e->GetEdge(edge)->m_id);
938 if (
m_mesh->m_spaceDim == 3)
940 vector<NodeSharedPtr> volNodes;
944 volNodes.resize((nq - 2) * (nq - 2));
945 for (
int j = 1; j < nq - 1; ++j)
947 for (
int k = 1; k < nq - 1; ++k)
950 volNodes[(j - 1) * (nq - 2) + (k - 1)] =
957 for (
int j = 3 + 3 * (nq - 2); j < nquad; ++j)
963 e->SetVolumeNodes(volNodes);
968 if (
m_mesh->m_expDim == 3)
972 for (it =
m_mesh->m_spherigonSurfs.begin();
973 it !=
m_mesh->m_spherigonSurfs.end();
977 m_mesh->m_element[
m_mesh->m_expDim][it->first]->GetFace(
980 f->m_faceNodes = el[elmt]->GetVolumeNodes();
981 f->m_curveType = f->m_vertexList.size() == 3
987 if (normalsGenerated)
989 m_mesh->m_vertexNormals.clear();
NEKMESHUTILS_EXPORT NekDouble dot(const Node &pSrc) const
#define ASSERTL0(condition, msg)
Basic information about an element.
pair< ModuleType, string > ModuleKey
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.
virtual ~ProcessSpherigon()
Destructor.
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
map< string, ConfigOption > m_config
List of configuration values.
Principle Modified Functions .
MeshSharedPtr m_mesh
Mesh object.
NekDouble m_y
Y-coordinate.
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
boost::shared_ptr< Node > NodeSharedPtr
Principle Orthogonal Functions .
Defines a specification for a set of points.
boost::shared_ptr< QuadExp > QuadExpSharedPtr
double CrossProdMag(Node &a, Node &b)
Calculate the magnitude of the cross product .
void UnitCrossProd(Node &a, Node &b, Node &c)
boost::shared_ptr< InputPly > InputPlySharedPtr
Represents a command-line configuration option.
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.
void SuperBlend(std::vector< double > &r, std::vector< Node > &Q, Node &P, std::vector< double > &blend)
Calculate the blending function for spherigon implementation.
void GenerateNormals(std::vector< ElementSharedPtr > &el, MeshSharedPtr &mesh)
Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hy...
boost::shared_ptr< TriGeom > TriGeomSharedPtr
boost::shared_ptr< Element > ElementSharedPtr
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
NekDouble m_z
Z-coordinate.
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.
double Blend(double r)
Calculate the blending function for spherigon implementation.
1D Gauss-Lobatto-Legendre quadrature points
2D Nodal Electrostatic Points on a Triangle
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
boost::shared_ptr< NodalTriExp > NodalTriExpSharedPtr