Commit 1bea063b authored by Dave Moxey's avatar Dave Moxey
Browse files

Fix a few cosmetic things

parent c3ae0ad6
////////////////////////////////////////////////////////////////////////////////
//
// File: SurfaceMeshing.cpp
// File: ProcessLoadCAD.cpp
//
// For more information, please see: http://www.nektar.info/
//
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: surfacemeshing object methods.
// Description: Load CAD module
//
////////////////////////////////////////////////////////////////////////////////
......
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessBL.h
// File: ProcessLoadCAD.h
//
// For more information, please see: http://www.nektar.info/
//
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Refine boundary layer of elements.
// Description: Load CAD module
//
////////////////////////////////////////////////////////////////////////////////
......
......@@ -53,7 +53,9 @@ void Octree::Process()
CompileSourcePointList();
if (m_mesh->m_verbose)
{
cout << "\tCurvature samples: " << m_SPList.size() << endl;
}
// make master octant based on the bounding box of the domain
m_dim = max((boundingBox[1] - boundingBox[0]) / 2.0,
......@@ -1021,14 +1023,15 @@ void Octree::CompileSourcePointList()
{
Array<OneD, NekDouble> x1(3), x2(3);
NekDouble r, d;
s >> x1[0] >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2] >> r >> d;
s >> x1[0] >> x1[1] >> x1[2] >> x2[0] >> x2[1] >> x2[2]
>> r >> d;
lsources.push_back(linesource(x1, x2, r, d));
}
}
fle.close();
//this takes any existing sourcepoints within the influence range
//and modifys them
// this takes any existing sourcepoints within the influence range
// and modifies them
for (int i = 0; i < m_SPList.size(); i++)
{
for (int j = 0; j < lsources.size(); j++)
......
......@@ -65,7 +65,7 @@ public:
: m_id(id), m_mesh(m)
{
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
};
}
/**
* @brief alternative constructor with mesh points already created
......@@ -74,7 +74,7 @@ public:
: m_id(id), m_mesh(m), m_meshpoints(n)
{
m_cadcurve = m_mesh->m_cad->GetCurve(m_id);
};
}
/**
* @brief execute meshing
......
......@@ -59,7 +59,9 @@ bool FaceMesh::ValidateCurves()
for(int i = 0; i < curvesInSurface.size(); i++)
{
vector<EdgeSharedPtr> es = m_curvemeshes[curvesInSurface[i]]->GetMeshEdges();
vector<EdgeSharedPtr> es =
m_curvemeshes[curvesInSurface[i]]->GetMeshEdges();
for(int j = i; j < curvesInSurface.size(); j++)
{
if(i == j)
......@@ -67,7 +69,8 @@ bool FaceMesh::ValidateCurves()
continue;
}
vector<EdgeSharedPtr> es2 = m_curvemeshes[curvesInSurface[j]]->GetMeshEdges();
vector<EdgeSharedPtr> es2 =
m_curvemeshes[curvesInSurface[j]]->GetMeshEdges();
for(int l = 0; l < es.size(); l++)
{
......@@ -83,16 +86,23 @@ bool FaceMesh::ValidateCurves()
continue;
}
Array<OneD, NekDouble> P3 = es2[k]->m_n1->GetCADSurfInfo(m_id);
Array<OneD, NekDouble> P4 = es2[k]->m_n2->GetCADSurfInfo(m_id);
Array<OneD, NekDouble> P3 =
es2[k]->m_n1->GetCADSurfInfo(m_id);
Array<OneD, NekDouble> P4 =
es2[k]->m_n2->GetCADSurfInfo(m_id);
NekDouble den = (P4[0]-P3[0])*(P2[1]-P1[1])
- (P2[0]-P1[0])*(P4[1]-P3[1]);
NekDouble den = (P4[0]-P3[0])*(P2[1]-P1[1]) - (P2[0]-P1[0])*(P4[1]-P3[1]);
if(fabs(den) < 1e-8)
{
continue;
}
NekDouble t = ((P1[0]-P3[0])*(P4[1]-P3[1]) - (P4[0]-P3[0])*(P1[1]-P3[1]))/den;
NekDouble u = (P1[0] - P3[0] + t*(P2[0]-P1[0])) / (P4[0] - P3[0]);
NekDouble t = ((P1[0]-P3[0])*(P4[1]-P3[1]) -
(P4[0]-P3[0])*(P1[1]-P3[1]))/den;
NekDouble u = (P1[0] - P3[0] + t*(P2[0]-P1[0])) /
(P4[0] - P3[0]);
if(t < 1.0 && t > 0.0 && u < 1.0 && u > 0.0)
{
......@@ -100,7 +110,9 @@ bool FaceMesh::ValidateCurves()
uv[0] = P1[0] + t * (P2[0] - P1[0]);
uv[1] = P1[1] + t * (P2[1] - P1[1]);
Array<OneD, NekDouble> loc = m_cadsurf->P(uv);
cout << endl << "Curve mesh error at " << loc[0] << " " << loc[1] << " " << loc[2] << " on face " << m_id << endl;
cout << endl << "Curve mesh error at " << loc[0] << " "
<< loc[1] << " " << loc[2] << " on face "
<< m_id << endl;
error = true;
}
}
......@@ -158,8 +170,8 @@ void FaceMesh::Mesh()
bool repeat = true;
int meshcounter = 1;
//continuously remesh until all triangles conform to the spacing in the
//octree
// continuously remesh until all triangles conform to the spacing in the
// octree
while (repeat)
{
repeat = Validate();
......@@ -174,8 +186,8 @@ void FaceMesh::Mesh()
meshcounter++;
}
//build a local version of the mesh (one set of triangles)
//this is done so edge connectivity infomration can be used for optimisation
// build a local version of the mesh (one set of triangles). this is done
// so edge connectivity infomration can be used for optimisation
BuildLocalMesh();
OptimiseLocalMesh();
......@@ -194,7 +206,8 @@ void FaceMesh::Mesh()
if(m_mesh->m_verbose)
{
cout << "\r "
cout << "\r "
" "
" ";
cout << scientific << "\r\t\tFace " << m_id << endl
<< "\t\t\tNodes: " << m_localNodes.size() << endl
......@@ -208,7 +221,7 @@ void FaceMesh::Mesh()
void FaceMesh::OptimiseLocalMesh()
{
//each optimisation algorithm is based on the work in chapeter 19
// each optimisation algorithm is based on the work in chapter 19
DiagonalSwap();
Smoothing();
......@@ -248,11 +261,15 @@ void FaceMesh::Smoothing()
for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
{
NodeSet::iterator f = m_inBoundary.find((*nit));
if (f != m_inBoundary.end()) // node is on curve so skip
// node is on curve so skip
if (f != m_inBoundary.end())
{
continue;
}
vector<NodeSharedPtr> connodes; // this can be real nodes or dummy
// nodes depending on the system
// this can be real nodes or dummy nodes depending on the system
vector<NodeSharedPtr> connodes;
vector<EdgeSharedPtr> edges = connectingedges[(*nit)->m_id];
vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
......@@ -266,11 +283,17 @@ void FaceMesh::Smoothing()
NodeSharedPtr J;
if (*nit == edges[i]->m_n1)
{
J = edges[i]->m_n2;
}
else if (*nit == edges[i]->m_n2)
{
J = edges[i]->m_n1;
}
else
{
ASSERTL0(false, "could not find node");
}
Array<OneD, NekDouble> ui = (*nit)->GetCADSurfInfo(m_id);
Array<OneD, NekDouble> uj = J->GetCADSurfInfo(m_id);
......@@ -278,9 +301,12 @@ void FaceMesh::Smoothing()
for (int j = 0; j < els.size(); j++)
{
vector<NodeSharedPtr> v = els[j]->GetVertexList();
// elememt is adjacent to J therefore no intersection on IJ
if (v[0] == J || v[1] == J || v[2] == J)
continue; // elememt is adjacent to J therefore no
// intersection on IJ
{
continue;
}
// need to find other edge
EdgeSharedPtr AtoB;
......@@ -314,8 +340,8 @@ void FaceMesh::Smoothing()
sort(lambda.begin(), lambda.end());
// make a new dummy node based on the system
Array<OneD, NekDouble> ud(2);
ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
Array<OneD, NekDouble> locd = m_cadsurf->P(ud);
NodeSharedPtr dn = boost::shared_ptr<Node>(
new Node(0, locd[0], locd[1], locd[2]));
......@@ -490,7 +516,8 @@ void FaceMesh::DiagonalSwap()
}
// determine signed area of alternate config
// //cout << A->GetNumCADSurf() << " " << B->GetNumCADSurf() << " " << C->GetNumCADSurf() << " " << D->GetNumCADSurf() << endl;
// //cout<< A->GetNumCADSurf() << " " << B->GetNumCADSurf() << " "
// // << C->GetNumCADSurf() << " " << D->GetNumCADSurf() << endl;
// ofstream file;
// file.open("pts.3D");
// file << "x y z value" << endl;
......@@ -619,7 +646,9 @@ void FaceMesh::DiagonalSwap()
for (int i = 0; i < links.size(); i++)
{
if (links[i].first->GetId() == tri1->GetId())
{
continue;
}
CA->m_elLink.push_back(links[i]);
}
......@@ -628,7 +657,9 @@ void FaceMesh::DiagonalSwap()
for (int i = 0; i < links.size(); i++)
{
if (links[i].first->GetId() == tri1->GetId())
{
continue;
}
BC->m_elLink.push_back(links[i]);
}
......@@ -637,7 +668,9 @@ void FaceMesh::DiagonalSwap()
for (int i = 0; i < links.size(); i++)
{
if (links[i].first->GetId() == tri2->GetId())
{
continue;
}
AD->m_elLink.push_back(links[i]);
}
......@@ -646,7 +679,9 @@ void FaceMesh::DiagonalSwap()
for (int i = 0; i < links.size(); i++)
{
if (links[i].first->GetId() == tri2->GetId())
{
continue;
}
DB->m_elLink.push_back(links[i]);
}
......@@ -756,7 +791,6 @@ void FaceMesh::BuildLocalMesh()
for (int i = 0; i < m_connec.size(); i++)
{
ElmtConfig conf(LibUtilities::eTriangle, 1, false, false);
vector<int> tags;
......@@ -834,9 +868,13 @@ void FaceMesh::Stretching()
uv[0] = bnds[0] + i * du;
uv[1] = bnds[2] + j * dv;
if (i == dxu - 1)
{
uv[0] = bnds[1];
}
if (j == dxv - 1)
uv[1] = bnds[3];
{
uv[1] = bnds[3];
}
Array<OneD, NekDouble> r = m_cadsurf->D1(uv);
NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
......@@ -846,7 +884,9 @@ void FaceMesh::Stretching()
rv *= dv;
if (rv < 1E-8)
{
continue;
}
m_str += ru / rv;
ct++;
......@@ -879,13 +919,19 @@ bool FaceMesh::Validate()
int numValid = 0;
if (r[0] < triDelta[0] && r[2] < triDelta[0])
{
numValid++;
}
if (r[1] < triDelta[1] && r[0] < triDelta[1])
{
numValid++;
}
if (r[2] < triDelta[2] && r[1] < triDelta[2])
{
numValid++;
}
if (numValid != 3)
{
......@@ -1063,7 +1109,9 @@ void FaceMesh::OrientateCurves()
nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
if (fabs(nt2[1] - nt1[1]) < 1e-30)
{
continue;
}
NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
......@@ -1079,7 +1127,9 @@ void FaceMesh::OrientateCurves()
nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
if (fabs(nt2[1] - nt1[1]) < 1e-30)
{
continue;
}
NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
......@@ -1101,7 +1151,9 @@ void FaceMesh::OrientateCurves()
nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
if (fabs(nt2[1] - nt1[1]) < 1e-30)
{
continue;
}
NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
......@@ -1117,7 +1169,9 @@ void FaceMesh::OrientateCurves()
nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
if (fabs(nt2[1] - nt1[1]) < 1e-30)
{
continue;
}
NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
......
////////////////////////////////////////////////////////////////////////////////
//
// File: SurfaceMeshing.cpp
// File: VolumeMesh.cpp
//
// For more information, please see: http://www.nektar.info/
//
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: surfacemeshing object methods.
// Description: Process volume meshing.
//
////////////////////////////////////////////////////////////////////////////////
......@@ -50,20 +50,16 @@ namespace NekMeshUtils
{
ModuleKey VolumeMesh::className = GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "volumemesh"),
VolumeMesh::create,
ModuleKey(eProcessModule, "volumemesh"), VolumeMesh::create,
"Generates a volume mesh");
VolumeMesh::VolumeMesh(MeshSharedPtr m) : ProcessModule(m)
{
m_config["blsurfs"] =
ConfigOption(false, "0", "Generate prisms on these surfs");
m_config["blthick"] =
ConfigOption(false, "0", "Prism layer thickness");
m_config["bllayers"] =
ConfigOption(false, "0", "Prism layers");
m_config["blprog"] =
ConfigOption(false, "0", "Prism progression");
m_config["blthick"] = ConfigOption(false, "0", "Prism layer thickness");
m_config["bllayers"] = ConfigOption(false, "0", "Prism layers");
m_config["blprog"] = ConfigOption(false, "0", "Prism progression");
}
VolumeMesh::~VolumeMesh()
......@@ -78,33 +74,31 @@ void VolumeMesh::Process()
bool makeBL;
vector<unsigned int> blSurfs;
if(m_config["blsurfs"].beenSet)
if (m_config["blsurfs"].beenSet)
{
makeBL = true;
makeBL = true;
m_mesh->m_numcomp = 2;
ParseUtils::GenerateSeqVector(m_config["blsurfs"].as<string>().c_str(),
blSurfs);
}
else
{
makeBL = false;
makeBL = false;
m_mesh->m_numcomp = 1;
}
NekDouble prefix = 100;
if(m_mesh->m_cad->GetNumSurf() > 100)
if (m_mesh->m_cad->GetNumSurf() > 100)
{
prefix*=10;
prefix *= 10;
}
TetMeshSharedPtr tet;
if (makeBL)
{
BLMeshSharedPtr blmesh = MemoryManager<BLMesh>::AllocateSharedPtr(
m_mesh, blSurfs,
m_config["blthick"].as<NekDouble>(),
m_config["bllayers"].as<int>(),
m_config["blprog"].as<NekDouble>());
m_mesh, blSurfs, m_config["blthick"].as<NekDouble>(),
m_config["bllayers"].as<int>(), m_config["blprog"].as<NekDouble>());
blmesh->Mesh();
/*ClearElementLinks();
......@@ -115,31 +109,30 @@ void VolumeMesh::Process()
ProcessComposites();
return;*/
//remesh the correct surfaces
// remesh the correct surfaces
vector<unsigned int> symsurfs = blmesh->GetSymSurfs();
vector<ElementSharedPtr> els = m_mesh->m_element[2];
vector<ElementSharedPtr> els = m_mesh->m_element[2];
m_mesh->m_element[2].clear();
for(int i = 0; i < els.size(); i++)
for (int i = 0; i < els.size(); i++)
{
vector<unsigned int>::iterator f = find(symsurfs.begin(),
symsurfs.end(),
els[i]->m_parentCAD->GetId());
vector<unsigned int>::iterator f = find(
symsurfs.begin(), symsurfs.end(), els[i]->m_parentCAD->GetId());
if(f == symsurfs.end())
if (f == symsurfs.end())
{
m_mesh->m_element[2].push_back(els[i]);
}
else
{
//remove element from links
// remove element from links
vector<EdgeSharedPtr> es = els[i]->GetEdgeList();
for(int j = 0; j < es.size(); j++)
for (int j = 0; j < es.size(); j++)
{
vector<pair<ElementSharedPtr,int> > lk = es[j]->m_elLink;
vector<pair<ElementSharedPtr, int> > lk = es[j]->m_elLink;
es[j]->m_elLink.clear();
for(int k = 0; k < lk.size(); k++)
for (int k = 0; k < lk.size(); k++)
{
if(lk[k].first == els[i])
if (lk[k].first == els[i])
{
continue;
}
......@@ -149,54 +142,56 @@ void VolumeMesh::Process()
}
}
for(int i = 0; i < symsurfs.size(); i++)
for (int i = 0; i < symsurfs.size(); i++)
{
set<int> cIds;
vector<EdgeLoop> e = m_mesh->m_cad->GetSurf(symsurfs[i])->GetEdges();
for(int i = 0; i < e.size(); i++)
vector<EdgeLoop> e =
m_mesh->m_cad->GetSurf(symsurfs[i])->GetEdges();
for (int i = 0; i < e.size(); i++)
{
for(int j = 0; j < e[i].edges.size(); j++)
for (int j = 0; j < e[i].edges.size(); j++)
{
cIds.insert(e[i].edges[j]->GetId());
}
}
//find the curve nodes which are on this symsurf
map<int,vector<NodeSharedPtr> > curveNodeMap;
// find the curve nodes which are on this symsurf
map<int, vector<NodeSharedPtr> > curveNodeMap;
NodeSet::iterator it;
for(it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end(); it++)
for (it = m_mesh->m_vertexSet.begin();
it != m_mesh->m_vertexSet.end(); it++)
{
vector<pair<int,CADCurveSharedPtr> > cc = (*it)->GetCADCurves();
for(int j = 0; j < cc.size(); j++)
vector<pair<int, CADCurveSharedPtr> > cc =
(*it)->GetCADCurves();
for (int j = 0; j < cc.size(); j++)
{
set<int>::iterator f = cIds.find(cc[j].first);
if(f != cIds.end())
if (f != cIds.end())
{
curveNodeMap[cc[j].first].push_back((*it));
}
}
}
//need to bubble sort the vectors
map<int,vector<NodeSharedPtr> >::iterator cit;
for(cit = curveNodeMap.begin(); cit != curveNodeMap.end(); cit++)
// need to bubble sort the vectors
map<int, vector<NodeSharedPtr> >::iterator cit;
for (cit = curveNodeMap.begin(); cit != curveNodeMap.end(); cit++)
{
vector<NekDouble> ts;
for(int i = 0; i < cit->second.size(); i++)
for (int i = 0; i < cit->second.size(); i++)
{
ts.push_back(cit->second[i]->GetCADCurveInfo(cit->first));
}
bool repeat = true;
while(repeat)
while (repeat)
{
repeat = false;
for(int i = 0; i < ts.size() - 1; i++)
for (int i = 0; i < ts.size() - 1; i++)
{
if(ts[i] > ts[i+1])
if (ts[i] > ts[i + 1])
{
swap(ts[i],ts[i+1]);
swap(cit->second[i],cit->second[i+1]);
swap(ts[i], ts[i + 1]);
swap(cit->second[i], cit->second[i + 1]);
repeat = true;
break;
}
......@@ -204,16 +199,18 @@ void VolumeMesh::Process()
}
}
//create quads
// create quads
map<NodeSharedPtr, NodeSharedPtr> nmap = blmesh->GetSymNodes();
for(cit = curveNodeMap.begin(); cit != curveNodeMap.end(); cit++)
for (cit = curveNodeMap.begin(); cit != curveNodeMap.end(); cit++)
{
for(int j = 0; j < cit->second.size() - 1; j++)
for (int j = 0; j < cit->second.size() - 1; j++)
{