Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | List of all members
Nektar::Utilities::ProcessSpherigon Class Reference

#include <ProcessSpherigon.h>

Inheritance diagram for Nektar::Utilities::ProcessSpherigon:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessSpherigon:
Collaboration graph
[legend]

Public Member Functions

 ProcessSpherigon (MeshSharedPtr m)
 Default constructor.
virtual ~ProcessSpherigon ()
 Destructor.
virtual void Process ()
 Write mesh to output file.
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule ()
 ProcessModule (FieldSharedPtr p_f)
 ProcessModule (MeshSharedPtr p_m)
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
virtual void Process (po::variables_map &vm)=0
void RegisterConfig (string key, string value)
 Register a configuration option with a module.
void PrintConfig ()
 Print out all configuration options for a module.
void SetDefaults ()
 Sets default configuration options for those which have not been set.
bool GetRequireEquiSpaced (void)
void SetRequireEquiSpaced (bool pVal)
void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 Module (MeshSharedPtr p_m)
void RegisterConfig (string key, string value)
void PrintConfig ()
void SetDefaults ()
MeshSharedPtr GetMesh ()
virtual void ProcessVertices ()
 Extract element vertices.

Static Public Member Functions

static boost::shared_ptr< Module > create (MeshSharedPtr m)
 Creates an instance of this class.

Static Public Attributes

static ModuleKey className

Protected Member Functions

void GenerateNormals (vector< ElementSharedPtr > &el, MeshSharedPtr &mesh)
 Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hybrid triangular/quadrilateral mesh in 3D.
double CrossProdMag (Node &a, Node &b)
 Calculate the magnitude of the cross product $ \vec{a}\times\vec{b} $.
void UnitCrossProd (Node &a, Node &b, Node &c)
double Blend (double r)
 Calculate the $ C^0 $ blending function for spherigon implementation.
void SuperBlend (vector< double > &r, vector< Node > &Q, Node &P, vector< double > &blend)
 Calculate the $ C^1 $ blending function for spherigon implementation.

Additional Inherited Members

- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object.
map< string, ConfigOption > m_config
 List of configuration values.
bool m_requireEquiSpaced
MeshSharedPtr m_mesh
 Mesh object.

Detailed Description

This class implements the spherigon surface smoothing technique which is documented in

"The SPHERIGON: A Simple Polygon Patch for Smoothing Quickly your Polygonal Meshes": P. Volino and N. Magnenat Thalmann, Computer Animation Proceedings (1998).

This implementation works in both a 2D manifold setting (for triangles and quadrilaterals embedded in 3-space) and in a full 3D enviroment (prisms, tetrahedra and hexahedra).

No additional information needs to be supplied directly to the module in order for it to work. However, 3D elements do rely on the mapping Mesh::spherigonSurfs to be populated by the relevant input modules.

Additionally, since the algorithm assumes normals are supplied which are perpendicular to the true surface defined at each vertex. If these are specified in Mesh::vertexNormals by the input module, better smoothing results can be obtained. Otherwise, normals are estimated by taking the average of all edge/face normals which connect to the vertex.

Definition at line 47 of file ProcessSpherigon.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessSpherigon::ProcessSpherigon ( MeshSharedPtr  m)

Default constructor.

Definition at line 90 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Module::m_config.

: ProcessModule(m)
{
m_config["N"] = ConfigOption(false, "5",
"Number of points to add to face edges.");
m_config["surf"] = ConfigOption(false, "-1",
"Tag identifying surface to process.");
m_config["BothTriFacesOnPrism"] = ConfigOption(true, "-1",
"Curve both triangular faces of prism on boundary.");
m_config["usenormalfile"] = ConfigOption(false,"NoFile",
"Use alternative file for Spherigon definition");
m_config["scalefile"] = ConfigOption(false,"1.0",
"Apply scaling factor to coordinates in file ");
m_config["normalnoise"] = ConfigOption(false,"NotSpecified",
"Add randowm noise to normals of amplitude AMP in specified region. input string is Amp,xmin,xmax,ymin,ymax,zmin,zmax");
}
Nektar::Utilities::ProcessSpherigon::~ProcessSpherigon ( )
virtual

Destructor.

Definition at line 109 of file ProcessSpherigon.cpp.

{
}

Member Function Documentation

double Nektar::Utilities::ProcessSpherigon::Blend ( double  r)
protected

Calculate the $ C^0 $ blending function for spherigon implementation.

See equation (5) of the paper.

Parameters
rBarycentric coordinate.

Definition at line 153 of file ProcessSpherigon.cpp.

{
return r*r;
}
static boost::shared_ptr<Module> Nektar::Utilities::ProcessSpherigon::create ( MeshSharedPtr  m)
inlinestatic

Creates an instance of this class.

Definition at line 51 of file ProcessSpherigon.h.

{
return MemoryManager<ProcessSpherigon>::AllocateSharedPtr(m);
}
double Nektar::Utilities::ProcessSpherigon::CrossProdMag ( Node &  a,
Node &  b 
)
protected

Calculate the magnitude of the cross product $ \vec{a}\times\vec{b} $.

Definition at line 137 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Node::m_x, Nektar::Utilities::Node::m_y, and Nektar::Utilities::Node::m_z.

Referenced by Process().

{
double tmp1 = a.m_y*b.m_z - a.m_z*b.m_y;
double tmp2 = a.m_z*b.m_x - a.m_x*b.m_z;
double tmp3 = a.m_x*b.m_y - a.m_y*b.m_x;
return sqrt(tmp1*tmp1 + tmp2*tmp2 + tmp3*tmp3);
}
void Nektar::Utilities::ProcessSpherigon::GenerateNormals ( std::vector< ElementSharedPtr > &  el,
MeshSharedPtr &  mesh 
)
protected

Generate a set of approximate vertex normals to a surface represented by line segments in 2D and a hybrid triangular/quadrilateral mesh in 3D.

This routine approximates the true vertex normals to a surface by averaging the normals of all edges/faces which connect to the vertex. It is better to use the exact surface normals which can be set in Mesh::vertexNormals, but where they are not supplied this routine calculates an approximation for the spherigon implementation.

Parameters
elVector of elements denoting the surface mesh.

Definition at line 223 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Node::abs2(), ASSERTL0, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTriangle, Nektar::iterator, Nektar::Utilities::Node::m_x, Nektar::Utilities::Node::m_y, Nektar::Utilities::Node::m_z, and UnitCrossProd().

Referenced by Process().

{
boost::unordered_map<int, Node>::iterator nIt;
for (int i = 0; i < el.size(); ++i)
{
ElementSharedPtr e = el[i];
// Ensure that element is a line, triangle or quad.
ASSERTL0(e->GetConf().m_e == LibUtilities::eSegment ||
e->GetConf().m_e == LibUtilities::eTriangle ||
e->GetConf().m_e == LibUtilities::eQuadrilateral,
"Spherigon expansions must be lines, triangles or "
"quadrilaterals.");
// Calculate normal for this element.
int nV = e->GetVertexCount();
vector<NodeSharedPtr> node(nV);
for (int j = 0; j < nV; ++j)
{
node[j] = e->GetVertex(j);
}
Node n;
if (mesh->m_spaceDim == 3)
{
// Create two tangent vectors and take unit cross product.
Node v1 = *(node[1]) - *(node[0]);
Node v2 = *(node[2]) - *(node[0]);
UnitCrossProd(v1, v2, n);
}
else
{
// Calculate gradient vector and invert.
Node dx = *(node[1]) - *(node[0]);
dx /= sqrt(dx.abs2());
n.m_x = -dx.m_y;
n.m_y = dx.m_x;
n.m_z = 0;
}
// Insert face normal into vertex normal list or add to existing
// value.
for (int j = 0; j < nV; ++j)
{
nIt = mesh->m_vertexNormals.find(e->GetVertex(j)->m_id);
if (nIt == mesh->m_vertexNormals.end())
{
mesh->m_vertexNormals[e->GetVertex(j)->m_id] = n;
}
else
{
nIt->second += n;
}
}
}
// Normalize resulting vectors.
for (nIt = mesh->m_vertexNormals.begin();
nIt != mesh->m_vertexNormals.end (); ++nIt)
{
Node &n = mesh->m_vertexNormals[nIt->first];
n /= sqrt(n.abs2());
}
}
void Nektar::Utilities::ProcessSpherigon::Process ( )
virtual

Write mesh to output file.

Perform the spherigon smoothing technique on the mesh.

Implements Nektar::Utilities::Module.

Definition at line 297 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Node::abs2(), ASSERTL0, ASSERTL1, CrossProdMag(), Nektar::Utilities::Node::dot(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eGaussRadauMAlpha1Beta0, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eNodalTriElec, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eSegment, Nektar::LibUtilities::eTriangle, GenerateNormals(), Nektar::ParseUtils::GenerateSeqVector(), Nektar::ParseUtils::GenerateUnOrderedVector(), GetEdge(), Nektar::Utilities::GetElementFactory(), Nektar::iterator, Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_mesh, Nektar::Utilities::Node::m_x, Nektar::Utilities::Node::m_y, Nektar::Utilities::Node::m_z, SuperBlend(), and Vmath::Zero().

{
ASSERTL0(m_mesh->m_spaceDim == 3 || m_mesh->m_spaceDim == 2,
"Spherigon implementation only valid in 2D/3D.");
boost::unordered_set<int>::iterator eIt;
boost::unordered_set<int> visitedEdges;
// First construct vector of elements to process.
vector<ElementSharedPtr> el;
if (m_mesh->m_verbose)
{
cout << "ProcessSpherigon: Smoothing mesh..." << endl;
}
if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
{
// Manifold case - copy expansion dimension.
el = m_mesh->m_element[m_mesh->m_expDim];
}
else if (m_mesh->m_spaceDim == m_mesh->m_expDim)
{
// Full 2D or 3D case - iterate over stored edges/faces and
// create segments/triangles representing those edges/faces.
set<pair<int,int> >::iterator it;
vector<int> t;
t.push_back(0);
// Construct list of spherigon edges/faces from a tag.
string surfTag = m_config["surf"].as<string>();
bool prismTag = m_config["BothTriFacesOnPrism"].beenSet;
if (surfTag != "")
{
vector<unsigned int> surfs;
ParseUtils::GenerateSeqVector(surfTag.c_str(), surfs);
sort(surfs.begin(), surfs.end());
m_mesh->m_spherigonSurfs.clear();
for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
{
ElementSharedPtr el = m_mesh->m_element[m_mesh->m_expDim][i];
int nSurf = m_mesh->m_expDim == 3 ? el->GetFaceCount() :
el->GetEdgeCount();
for (int j = 0; j < nSurf; ++j)
{
int bl = el->GetBoundaryLink(j);
if (bl == -1)
{
continue;
}
ElementSharedPtr bEl = m_mesh->m_element[m_mesh->m_expDim-1][bl];
vector<int> tags = bEl->GetTagList();
vector<int> inter;
sort(tags.begin(), tags.end());
set_intersection(surfs.begin(), surfs.end(),
tags .begin(), tags .end(),
back_inserter(inter));
if (inter.size() == 1)
{
m_mesh->m_spherigonSurfs.insert(make_pair(i, j));
// Curve other tri face on Prism. Note could be
// problem on pyramid when implemented.
if(nSurf == 5 && prismTag)
{
// add other end of prism on boundary for
// smoothing
int triFace = j == 1 ? 3 : 1;
m_mesh->m_spherigonSurfs.insert(
make_pair(i, triFace));
}
}
}
}
}
if (m_mesh->m_spherigonSurfs.size() == 0)
{
cerr << "WARNING: Spherigon surfaces have not been defined "
<< "-- ignoring smoothing." << endl;
return;
}
if (m_mesh->m_expDim == 3)
{
for (it = m_mesh->m_spherigonSurfs.begin();
it != m_mesh->m_spherigonSurfs.end (); ++it)
{
FaceSharedPtr f = m_mesh->m_element[m_mesh->m_expDim][it->first]->
GetFace(it->second);
vector<NodeSharedPtr> nodes = f->m_vertexList;
LibUtilities::ShapeType eType = (LibUtilities::ShapeType)(nodes.size());
ElmtConfig conf(eType, 1, false, false);
// Create 2D element.
ElementSharedPtr elmt = GetElementFactory().
CreateInstance(eType,conf,nodes,t);
// Copy vertices/edges from face.
for (int i = 0; i < f->m_vertexList.size(); ++i)
{
elmt->SetVertex(i, f->m_vertexList[i]);
}
for (int i = 0; i < f->m_edgeList.size(); ++i)
{
elmt->SetEdge(i, f->m_edgeList[i]);
}
el.push_back(elmt);
}
}
else
{
for (it = m_mesh->m_spherigonSurfs.begin();
it != m_mesh->m_spherigonSurfs.end (); ++it)
{
EdgeSharedPtr edge = m_mesh->m_element[m_mesh->m_expDim][it->first]->
GetEdge(it->second);
vector<NodeSharedPtr> nodes;
LibUtilities::ShapeType eType = LibUtilities::eSegment;
ElmtConfig conf(eType, 1, false, false);
nodes.push_back(edge->m_n1);
nodes.push_back(edge->m_n2);
// Create 2D element.
ElementSharedPtr elmt = GetElementFactory().
CreateInstance(eType,conf,nodes,t);
// Copy vertices/edges from original element.
elmt->SetVertex(0, nodes[0]);
elmt->SetVertex(1, nodes[1]);
elmt->SetEdge(0, m_mesh->m_element[m_mesh->m_expDim][it->first]->
GetEdge(it->second));
el.push_back(elmt);
}
}
}
else
{
ASSERTL0(false, "Spherigon expansions must be 2/3 dimensional");
}
// See if vertex normals have been generated. If they have not,
// approximate them by summing normals of surrounding elements.
bool normalsGenerated = false;
// Read Normal file if one exists
std::string normalfile = m_config["usenormalfile"].as<string>();
if(normalfile.compare("NoFile") != 0)
{
NekDouble scale = m_config["scalefile"].as<NekDouble>();
if (m_mesh->m_verbose)
{
cout << "Inputing normal file: " << normalfile << " with scaling of " << scale << endl;
}
ifstream inply;
InputPlySharedPtr plyfile;
inply.open(normalfile.c_str());
int j;
MeshSharedPtr m = boost::shared_ptr<Mesh>(new Mesh());
plyfile = boost::shared_ptr<InputPly>(new InputPly(m));
plyfile->ReadPly(inply,scale);
plyfile->ProcessVertices();
MeshSharedPtr plymesh = plyfile->GetMesh();
GenerateNormals(plymesh->m_element[plymesh->m_expDim],plymesh);
// finaly find nearest vertex and set normal to mesh surface file normal.
// probably should have a hex tree search ?
Array<OneD, NekDouble> len2(plymesh->m_vertexSet.size());
Node minx(0,0.0,0.0,0.0), tmp,tmpsav;
NodeSet::iterator it;
map<int,NodeSharedPtr>::iterator vIt;
map<int,NodeSharedPtr> surfverts;
// make a map of normal vertices to visit based on elements el
for (int i = 0; i < el.size(); ++i)
{
ElementSharedPtr e = el[i];
int nV = e->GetVertexCount();
for (int j = 0; j < nV; ++j)
{
int id = e->GetVertex(j)->m_id;
surfverts[id] = e->GetVertex(j);
}
}
//loop over all element in ply mesh and determine
//xmin,xmax,ymin,ymax as search criterion
NekDouble mindiff,diff;
int cntmin;
if (m_mesh->m_verbose)
{
cout << "\t Processing surface normals " << endl;
}
int cnt = 0;
map<int,int> locnorm;
for(vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt,++cnt)
{
mindiff = 1e12;
for(j = 0, it = plymesh->m_vertexSet.begin(); it != plymesh->m_vertexSet.end(); ++it, ++j)
{
tmp = *(vIt->second)- *(*it);
diff = tmp.abs2();
if(diff < mindiff)
{
mindiff = diff;
cntmin = (*it)->m_id;
tmpsav = tmp;
}
}
locnorm[cntmin] = vIt->first;
ASSERTL1(cntmin < plymesh->m_vertexNormals.size(),"cntmin is out of range");
m_mesh->m_vertexNormals[vIt->first] = plymesh->m_vertexNormals[cntmin];
}
if (m_mesh->m_verbose)
{
cout << "\t end of processing surface normals " << endl;
}
normalsGenerated = true;
}
else if (m_mesh->m_vertexNormals.size() == 0)
{
GenerateNormals(el,m_mesh);
normalsGenerated = true;
}
// See if we should add noise to normals
std::string normalnoise = m_config["normalnoise"].as<string>();
if(normalnoise.compare("NotSpecified") != 0)
{
vector<NekDouble> values;
ASSERTL0(ParseUtils::GenerateUnOrderedVector(normalnoise.c_str(),values),"Failed to interpret normal noise string");
int nvalues = values.size()/2;
NekDouble amp = values[0];
if (m_mesh->m_verbose)
{
cout << "\t adding noise to normals of amplitude "<< amp << " in range: ";
for(int i = 0; i < nvalues; ++i)
{
cout << values[2*i+1] <<"," << values[2*i+2] << " ";
}
cout << endl;
}
map<int,NodeSharedPtr>::iterator vIt;
map<int,NodeSharedPtr> surfverts;
// make a map of normal vertices to visit based on elements el
for (int i = 0; i < el.size(); ++i)
{
ElementSharedPtr e = el[i];
int nV = e->GetVertexCount();
for (int j = 0; j < nV; ++j)
{
int id = e->GetVertex(j)->m_id;
surfverts[id] = e->GetVertex(j);
}
}
for(vIt = surfverts.begin(); vIt != surfverts.end(); ++vIt)
{
bool AddNoise = false;
for(int i = 0; i < nvalues; ++i)
{
// check to see if point is in range
switch(nvalues)
{
case 1:
{
if(((vIt->second)->m_x > values[2*i+1])&&((vIt->second)->m_x < values[2*i+2]))
{
AddNoise = true;
}
}
break;
case 2:
{
if(((vIt->second)->m_x > values[2*i+1])&&
((vIt->second)->m_x < values[2*i+2])&&
((vIt->second)->m_y > values[2*i+3])&&
((vIt->second)->m_y < values[2*i+4]))
{
AddNoise = true;
}
}
break;
case 3:
{
if(((vIt->second)->m_x > values[2*i+1])&&
((vIt->second)->m_x < values[2*i+2])&&
((vIt->second)->m_y > values[2*i+3])&&
((vIt->second)->m_y < values[2*i+4])&&
((vIt->second)->m_z > values[2*i+5])&&
((vIt->second)->m_z < values[2*i+6]))
{
AddNoise = true;
}
}
break;
}
if(AddNoise)
{
// generate random unit vector;
Node rvec(0,rand(),rand(),rand());
rvec *= values[0]/sqrt(rvec.abs2());
Node normal = m_mesh->m_vertexNormals[vIt->first];
normal += rvec;
normal /= sqrt(normal.abs2());
m_mesh->m_vertexNormals[vIt->first] = normal;
}
}
}
}
// Allocate storage for interior points.
int nq = m_config["N"].as<int>();
int nquad = m_mesh->m_spaceDim == 3 ? nq*nq : nq;
Array<OneD, NekDouble> x(nq*nq);
Array<OneD, NekDouble> y(nq*nq);
Array<OneD, NekDouble> z(nq*nq);
Array<OneD, NekDouble> xc(nq*nq);
Array<OneD, NekDouble> yc(nq*nq);
Array<OneD, NekDouble> zc(nq*nq);
ASSERTL0(nq > 2, "Number of points must be greater than 2.");
LibUtilities::BasisKey B0(
LibUtilities::eOrtho_A, nq,
LibUtilities::PointsKey(
nq, LibUtilities::eGaussLobattoLegendre));
LibUtilities::BasisKey B1(
LibUtilities::eOrtho_B, nq,
LibUtilities::PointsKey(
nq, LibUtilities::eGaussRadauMAlpha1Beta0));
StdRegions::StdNodalTriExpSharedPtr stdtri =
MemoryManager<StdRegions::StdNodalTriExp>::AllocateSharedPtr(
B0, B1, LibUtilities::eNodalTriElec);
Array<OneD, NekDouble> xnodal(nq*(nq+1)/2), ynodal(nq*(nq+1)/2);
stdtri->GetNodalPoints(xnodal, ynodal);
int edgeMap[3][4][2] = {
{{0, 1}, {-1, -1}, {-1, -1 }, {-1, -1}}, // seg
{{0, 1}, {nq-1, nq}, {nq*(nq-1), -nq}, {-1, -1}}, // tri
{{0, 1}, {nq-1, nq}, {nq*nq-1, -1 }, {nq*(nq-1), -nq}} // quad
};
int vertMap[3][4][2] = {
{{0, 1}, {0, 0}, {0, 0}, {0, 0}}, // seg
{{0, 1}, {1, 2}, {2, 3}, {0, 0}}, // tri
{{0, 1}, {1, 2}, {2, 3}, {3, 0}}, // quad
};
for (int i = 0; i < el.size(); ++i)
{
// Construct a Nektar++ element to obtain coordinate points
// inside the element. TODO: Add options for various
// nodal/tensor point distributions + number of points to add.
ElementSharedPtr e = el[i];
LibUtilities::BasisKey B2(
LibUtilities::eModified_A, nq,
LibUtilities::PointsKey(
nq, LibUtilities::eGaussLobattoLegendre));
if (e->GetConf().m_e == LibUtilities::eSegment)
{
SpatialDomains::SegGeomSharedPtr geom =
boost::dynamic_pointer_cast<SpatialDomains::SegGeom>(
e->GetGeom(m_mesh->m_spaceDim));
LocalRegions::SegExpSharedPtr seg =
MemoryManager<LocalRegions::SegExp>::AllocateSharedPtr(
B2, geom);
seg->GetCoords(x,y,z);
nquad = nq;
}
else if (e->GetConf().m_e == LibUtilities::eTriangle)
{
SpatialDomains::TriGeomSharedPtr geom =
boost::dynamic_pointer_cast<SpatialDomains::TriGeom>(
e->GetGeom(3));
LocalRegions::NodalTriExpSharedPtr tri =
MemoryManager<LocalRegions::NodalTriExp>
::AllocateSharedPtr(
B0, B1, LibUtilities::eNodalTriElec, geom);
Array<OneD, NekDouble> coord(2);
tri->GetCoords(xc,yc,zc);
nquad = nq*(nq+1)/2;
for (int j = 0; j < nquad; ++j)
{
coord[0] = xnodal[j];
coord[1] = ynodal[j];
x[j] = stdtri->PhysEvaluate(coord, xc);
}
for (int j = 0; j < nquad; ++j)
{
coord[0] = xnodal[j];
coord[1] = ynodal[j];
y[j] = stdtri->PhysEvaluate(coord, yc);
}
for (int j = 0; j < nquad; ++j)
{
coord[0] = xnodal[j];
coord[1] = ynodal[j];
z[j] = stdtri->PhysEvaluate(coord, zc);
}
}
else if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
{
SpatialDomains::QuadGeomSharedPtr geom =
boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
e->GetGeom(3));
LocalRegions::QuadExpSharedPtr quad =
MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(
B2, B2, geom);
quad->GetCoords(x,y,z);
nquad = nq*nq;
}
else
{
ASSERTL0(false, "Unknown expansion type.");
}
// Zero z-coordinate in 2D.
if (m_mesh->m_spaceDim == 2)
{
Vmath::Zero(nquad, z, 1);
}
// Find vertex normals.
int nV = e->GetVertexCount();
vector<Node> v, vN;
for (int j = 0; j < nV; ++j)
{
v.push_back(*(e->GetVertex(j)));
ASSERTL1(m_mesh->m_vertexNormals.count(v[j].m_id) != 0,"Normal has not been defined");
vN.push_back(m_mesh->m_vertexNormals[v[j].m_id]);
}
vector<Node> tmp (nV);
vector<double> r (nV);
vector<Node> K (nV);
vector<Node> Q (nV);
vector<Node> Qp (nV);
vector<double> blend(nV);
vector<Node> out (nquad);
// Calculate segment length for 2D spherigon routine.
double segLength = sqrt((v[0] - v[1]).abs2());
// Perform Spherigon method to smooth manifold.
for (int j = 0; j < nquad; ++j)
{
Node P(0, x[j], y[j], z[j]);
Node N(0,0,0,0);
// Calculate generalised barycentric coordinates r[] and the
// Phong normal N = vN . r for this point of the element.
if (m_mesh->m_spaceDim == 2)
{
// In 2D the coordinates are given by a ratio of the
// segment length to the distance from one of the
// endpoints.
r[0] = sqrt((P - v[0]).abs2()) / segLength;
r[0] = max(min(1.0, r[0]), 0.0);
r[1] = 1.0 - r[0];
// Calculate Phong normal.
N = vN[0]*r[0] + vN[1]*r[1];
}
else if (m_mesh->m_spaceDim == 3)
{
for (int k = 0; k < nV; ++k)
{
tmp[k] = P - v[k];
}
// Calculate generalized barycentric coordinate system
// (see equation 6 of paper).
double weight = 0.0;
for (int k = 0; k < nV; ++k)
{
r[k] = 1.0;
for (int l = 0; l < nV-2; ++l)
{
r[k] *= CrossProdMag(tmp[(k+l+1) % nV],
tmp[(k+l+2) % nV]);
}
weight += r[k];
}
// Calculate Phong normal (equation 1).
for (int k = 0; k < nV; ++k)
{
r[k] /= weight;
N += vN[k]*r[k];
}
}
// Normalise Phong normal.
N /= sqrt(N.abs2());
for (int k = 0; k < nV; ++k)
{
// Perform steps denoted in equations 2, 3, 8 for C1
// smoothing.
double tmp1;
K [k] = P+N*((v[k]-P).dot(N));
tmp1 = (v[k]-K[k]).dot(vN[k]) / (1.0 + N.dot(vN[k]));
Q [k] = K[k] + N*tmp1;
Qp[k] = v[k] - N*((v[k]-P).dot(N));
}
// Apply C1 blending function to the surface. TODO: Add
// option to do (more efficient) C0 blending function.
SuperBlend(r, Qp, P, blend);
P.m_x = P.m_y = P.m_z = 0.0;
// Apply blending (equation 4).
for (int k = 0; k < nV; ++k)
{
P += Q[k]*blend[k];
}
out[j] = P;
}
// Push nodes into lines - TODO: face interior nodes.
// offset = 0 (seg), 1 (tri) or 2 (quad)
int offset = (int)e->GetConf().m_e-2;
for (int edge = 0; edge < e->GetEdgeCount(); ++edge)
{
eIt = visitedEdges.find(e->GetEdge(edge)->m_id);
if (eIt == visitedEdges.end())
{
bool reverseEdge = !(v[vertMap[offset][edge][0]] ==
*(e->GetEdge(edge)->m_n1));
// Clear existing curvature.
e->GetEdge(edge)->m_edgeNodes.clear();
if (e->GetConf().m_e != LibUtilities::eTriangle)
{
for (int j = 1; j < nq-1; ++j)
{
int v = edgeMap[offset][edge][0] +
j*edgeMap[offset][edge][1];
e->GetEdge(edge)->m_edgeNodes.push_back(
NodeSharedPtr(new Node(out[v])));
}
}
else
{
for (int j = 0; j < nq-2; ++j)
{
int v = 3 + edge*(nq-2) + j;
e->GetEdge(edge)->m_edgeNodes.push_back(
NodeSharedPtr(new Node(out[v])));
}
}
if (reverseEdge)
{
reverse(e->GetEdge(edge)->m_edgeNodes.begin(),
e->GetEdge(edge)->m_edgeNodes.end());
}
e->GetEdge(edge)->m_curveType =
LibUtilities::eGaussLobattoLegendre;
visitedEdges.insert(e->GetEdge(edge)->m_id);
}
}
// Add face nodes in manifold and full 3D case, but not for 2D.
if (m_mesh->m_spaceDim == 3)
{
vector<NodeSharedPtr> volNodes;
if (e->GetConf().m_e == LibUtilities::eQuadrilateral)
{
volNodes.resize((nq-2)*(nq-2));
for (int j = 1; j < nq-1; ++j)
{
for (int k = 1; k < nq-1; ++k)
{
int v = j*nq+k;
volNodes[(j-1)*(nq-2)+(k-1)] =
NodeSharedPtr(new Node(out[v]));
}
}
}
else
{
for (int j = 3+3*(nq-2); j < nquad; ++j)
{
volNodes.push_back(NodeSharedPtr(new Node(out[j])));
}
}
e->SetVolumeNodes(volNodes);
}
}
// Copy face nodes back into 3D element faces
if (m_mesh->m_expDim == 3)
{
set<pair<int,int> >::iterator it;
int elmt = 0;
for (it = m_mesh->m_spherigonSurfs.begin();
it != m_mesh->m_spherigonSurfs.end (); ++it, ++elmt)
{
FaceSharedPtr f = m_mesh->m_element[m_mesh->m_expDim][it->first]->
GetFace(it->second);
f->m_faceNodes = el[elmt]->GetVolumeNodes();
f->m_curveType = f->m_vertexList.size() == 3 ?
LibUtilities::eNodalTriElec :
LibUtilities::eGaussLobattoLegendre;
}
}
if (normalsGenerated)
{
m_mesh->m_vertexNormals.clear();
}
}
void Nektar::Utilities::ProcessSpherigon::SuperBlend ( vector< double > &  r,
vector< Node > &  Q,
Node &  P,
vector< double > &  blend 
)
protected

Calculate the $ C^1 $ blending function for spherigon implementation.

See equation (10) of the paper.

Parameters
rGeneralised barycentric coordinates of the point P.
QVector of vertices denoting this triangle/quad.
PPoint in the triangle to apply blending to.
blendThe resulting blending components for each vertex.

Definition at line 169 of file ProcessSpherigon.cpp.

References TOL_BLEND.

Referenced by Process().

{
vector<double> tmp(r.size());
double totBlend = 0.0;
int i;
int nV = r.size();
for (i = 0; i < nV; ++i)
{
blend[i] = 0.0;
tmp [i] = (Q[i]-P).abs2();
}
for (i = 0; i < nV; ++i)
{
int ip = (i+1) % nV, im = (i-1+nV) % nV;
if (r[i] > TOL_BLEND && r[i] < (1-TOL_BLEND))
{
blend[i] = r[i]*r[i]*(
r[im]*r[im]*tmp[im]/(tmp[im] + tmp[i]) +
r[ip]*r[ip]*tmp[ip]/(tmp[ip] + tmp[i]));
totBlend += blend[i];
}
}
for (i = 0; i < nV; ++i)
{
blend[i] /= totBlend;
if (r[i] >= (1-TOL_BLEND))
{
blend[i] = 1.0;
}
if (r[i] <= TOL_BLEND)
{
blend[i] = 0.0;
}
}
}
void Nektar::Utilities::ProcessSpherigon::UnitCrossProd ( Node &  a,
Node &  b,
Node &  c 
)
protected

Definition at line 118 of file ProcessSpherigon.cpp.

References Nektar::Utilities::Node::m_x, Nektar::Utilities::Node::m_y, and Nektar::Utilities::Node::m_z.

Referenced by GenerateNormals().

{
double inv_mag;
c.m_x = a.m_y*b.m_z - a.m_z*b.m_y;
c.m_y = a.m_z*b.m_x - a.m_x*b.m_z;
c.m_z = a.m_x*b.m_y - a.m_y*b.m_x;
inv_mag = 1.0/sqrt(c.m_x*c.m_x + c.m_y*c.m_y + c.m_z*c.m_z);
c.m_x = c.m_x*inv_mag;
c.m_y = c.m_y*inv_mag;
c.m_z = c.m_z*inv_mag;
}

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessSpherigon::className
static
Initial value:
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "spherigon"),
ProcessSpherigon::create)

Definition at line 55 of file ProcessSpherigon.h.

Generated on Wed Mar 4 2015 16:33:51 for Nektar++ by   doxygen 1.8.1.2