Commit fb736927 authored by Michael Turner's avatar Michael Turner

clean up

parent aa368f8d
......@@ -46,8 +46,10 @@ void BLMesh::Mesh()
//At this stage the surface mesh is complete and the elements know their
//neigbours through element links in the edges, this includes quads,
//here elements are made for the boundary layer they will need to know links (maybe facelinks)
//so that the tetmeshing module can extract the surface upon which it needs to mesh (top of the bl and the rest of the surface)
//here elements are made for the boundary layer they will need to know
//links (maybe facelinks)
//so that the tetmeshing module can extract the surface upon which it needs
// to mesh (top of the bl and the rest of the surface)
vector<ElementSharedPtr> quad;
vector<ElementSharedPtr> ptri; //triangles to grow prisms onto
......@@ -63,7 +65,8 @@ void BLMesh::Mesh()
break;
}
}
if(m_mesh->m_element[2][i]->GetConf().m_e == LibUtilities::eQuadrilateral)
if(m_mesh->m_element[2][i]->GetConf().m_e ==
LibUtilities::eQuadrilateral)
{
quad.push_back(m_mesh->m_element[2][i]);
}
......@@ -79,8 +82,11 @@ void BLMesh::Mesh()
vector<EdgeSharedPtr> e = quad[i]->GetEdgeList();
for(int j = 0; j < e.size(); j++)
{
if((e[j]->m_n1->GetNumCadCurve() > 0 && e[j]->m_n2->GetNumCadCurve() > 0) ||
(!(e[j]->m_n1->GetNumCadCurve() > 0) && !(e[j]->m_n2->GetNumCadCurve() > 0))) //if both or none are on curve skip
//if both or none are on curve skip
if((e[j]->m_n1->GetNumCadCurve() > 0 &&
e[j]->m_n2->GetNumCadCurve() > 0) ||
(!(e[j]->m_n1->GetNumCadCurve() > 0) &&
!(e[j]->m_n2->GetNumCadCurve() > 0)))
{
continue;
}
......@@ -150,17 +156,24 @@ void BLMesh::Mesh()
vector<pair<int, CADSurfSharedPtr> > surfs = n[j]->GetCADSurfs();
for(int s = 0; s < surfs.size(); s++)
{
Array<OneD, NekDouble> N = surfs[s].second->N(n[j]->GetCADSurfInfo(surfs[s].first));
for(int k = 0; k < 3; k++) AN[k] += N[k];
Array<OneD, NekDouble> N =
surfs[s].second->N(n[j]->GetCADSurfInfo(surfs[s].first));
for(int k = 0; k < 3; k++)
{
AN[k] += N[k];
}
}
NekDouble mag = sqrt(AN[0]*AN[0] + AN[1]*AN[1] + AN[2]*AN[2]);
for(int k = 0; k < 3; k++) AN[k] /= mag;
for(int k = 0; k < 3; k++)
{
AN[k] /= mag;
}
Array<OneD, NekDouble> loc = n[j]->GetLoc();
Array<OneD, NekDouble> np(3);
for(int k = 0; k < 3; k++) np[k] = loc[k] + AN[k]*m_bl;
NodeSharedPtr nn = boost::shared_ptr<Node>(new Node(m_mesh->m_numNodes++,
np[0],np[1],np[2]));
NodeSharedPtr nn = boost::shared_ptr<Node>(new Node(
m_mesh->m_numNodes++, np[0], np[1], np[2]));
pn[nm[j*2+1]] = nn;
blpair[n[j]->m_id] = nn;
}
......@@ -174,7 +187,8 @@ void BLMesh::Mesh()
m_mesh->m_element[3].push_back(E);
//need to give the surface element some information about which prism is above it
//need to give the surface element some information about
//which prism is above it
//so that tetmesh can infer the psudo surface
vector<NodeSharedPtr> faceNodes;
vector<EdgeSharedPtr> edgeList = ptri[i]->GetEdgeList();
......@@ -187,7 +201,9 @@ void BLMesh::Mesh()
if(f[j]->m_vertexList.size() != 3) //quad
continue;
if(!(F == f[j])) //only two triangle faces so if its not this one, this is the psudo surfaces
//only two triangle faces so if its not this one,
//this is the psudo surfaces
if(!(F == f[j]))
{
m_surftopriface[ptri[i]->GetId()] = f[j];
}
......
////////////////////////////////////////////////////////////////////////////////
//
// File: TetMesh.h
// File: BLMesh.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: class for tet meshing
// Description: class for boundary layer meshing
//
////////////////////////////////////////////////////////////////////////////////
......@@ -58,14 +58,14 @@ public:
*@brief default constructor
*/
BLMesh(MeshSharedPtr m, std::vector<unsigned int> bls,
std::vector<unsigned int> syms, NekDouble b) :
m_mesh(m), m_blsurfs(bls), m_symsurfs(syms), m_bl(b)
std::vector<unsigned int> syms, NekDouble b) :
m_mesh(m), m_blsurfs(bls), m_symsurfs(syms), m_bl(b)
{
};
/**
*@brief execute tet meshing
*@brief execute bl meshing
*/
void Mesh();
......
......@@ -33,8 +33,8 @@
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NekMeshUtils_CADSYSTEM_CADCURVE
#define NekMeshUtils_CADSYSTEM_CADCURVE
#ifndef NEKMESHUTILS_CADSYSTEM_CADCURVE
#define NEKMESHUTILS_CADSYSTEM_CADCURVE
#include <boost/shared_ptr.hpp>
......@@ -88,17 +88,19 @@ public:
NekDouble Length(NekDouble ti, NekDouble tf);
/**
* @brief Gets the location (x,y,z) in an array out of the curve at point \p t.
* @brief Gets the location (x,y,z) in an array out of the curve at
* point \p t.
*
* @param t Parametric coordinate
* @return Array of x,y,z
*/
Array<OneD, NekDouble> P(NekDouble t);
/**
* @brief Gets the second derivatives at t
*/
Array<OneD, NekDouble> D2(NekDouble t);
/**
* @brief Calculates the parametric coordinate and arclength location
* defined by \p s.
......@@ -140,7 +142,8 @@ public:
m_mainVerts = falVert;
}
/// get the ids of the vertices that are the ends of the curve, which are in the main cad list
/// get the ids of the vertices that are the ends of the curve,
/// which are in the main cad list
std::vector<CADVertSharedPtr> GetVertex()
{
return m_mainVerts;
......
......@@ -83,7 +83,7 @@ public:
protected:
/// ID of the vert.
int m_id;
///
/// type of the cad object
cadType m_type;
};
......
......@@ -54,7 +54,8 @@ class CADCurve;
typedef boost::shared_ptr<CADCurve> CADCurveSharedPtr;
/**
* @brief struct which descibes a collection of cad edges which for a loop on the cad surface
* @brief struct which descibes a collection of cad edges which for a
* loop on the cad surface
*/
struct EdgeLoop
{
......@@ -202,7 +203,8 @@ private:
void Test(Array<OneD, NekDouble> uv);
/// normal
bool m_correctNormal;
/// flag to alert the mesh generation to a potential problem is both curves have only two points in the mesh
/// flag to alert the mesh generation to a potential problem is both
/// curves have only two points in the mesh
bool m_hasTwoCurves;
/// OpenCascade object for surface.
BRepAdaptor_Surface m_occSurface;
......
......@@ -10,7 +10,6 @@ SET(NEKMESHUTILS_SOURCES
MeshElements/Pyramid.cpp
MeshElements/Prism.cpp
MeshElements/Hexahedron.cpp
Optimisation/Guass-Seidel.cpp
Optimisation/BGFS-B.cpp
)
IF(NEKTAR_USE_MESH)
......@@ -49,7 +48,6 @@ SET(NEKMESHUTILS_HEADERS
MeshElements/Pyramid.h
MeshElements/Prism.h
MeshElements/Hexahedron.h
Optimisation/Guass-Seidel.h
Optimisation/BGFS-B.h
Optimisation/OptimiseObj.h
)
......
......@@ -78,7 +78,8 @@ class TetGenInterface
void GetNewPoints(int num, std::vector<Array<OneD, NekDouble> > &newp);
/**
* @brief refines a previously made tetmesh with node delta information from the Octree
* @brief refines a previously made tetmesh with node delta information
* from the Octree
*/
void RefineMesh(std::map<int, NekDouble> delta);
......
......@@ -136,7 +136,8 @@ void TriangleInterface::Mesh(bool Quiet, bool Quality)
if(out.numberofpoints-out.numberofedges+out.numberoftriangles != 2 - m_centers.size())
{
cout << endl << "epc wrong" << endl;
cout << out.numberofpoints-out.numberofedges+out.numberoftriangles << " " << m_centers.size() << " " << sid << endl;
cout << out.numberofpoints-out.numberofedges+out.numberoftriangles
<< " " << m_centers.size() << " " << sid << endl;
}
}
......@@ -212,7 +213,8 @@ void TriangleInterface::Extract(std::vector<std::vector<NodeSharedPtr> > &Connec
n2 = nodemap.find(out.trianglelist[i*3+1]);
n3 = nodemap.find(out.trianglelist[i*3+2]);
ASSERTL0(n1!=nodemap.end() && n2!=nodemap.end() && n3!=nodemap.end(),"node index error");
ASSERTL0(n1!=nodemap.end() && n2!=nodemap.end() && n3!=nodemap.end(),
"node index error");
vector<NodeSharedPtr> tri(3);
tri[0] = n1->second;
......
......@@ -83,7 +83,8 @@ namespace NekMeshUtils
return m_edgeNodes.size() + 2;
}
NEKMESHUTILS_EXPORT void GetCurvedNodes(std::vector<NodeSharedPtr> &nodeList) const
NEKMESHUTILS_EXPORT void GetCurvedNodes(
std::vector<NodeSharedPtr> &nodeList) const
{
nodeList.push_back(m_n1);
for (int k = 0; k < m_edgeNodes.size(); ++k)
......@@ -117,7 +118,8 @@ namespace NekMeshUtils
}
/// Generate a SpatialDomains::SegGeom object for this edge.
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(int coordDim)
NEKMESHUTILS_EXPORT SpatialDomains::SegGeomSharedPtr GetGeom(
int coordDim)
{
// Create edge vertices.
SpatialDomains::PointGeomSharedPtr p[2];
......
......@@ -240,9 +240,8 @@ namespace NekMeshUtils
l += jn(i,0) * dk[i];
}
//cout << lam << ": " << fo << " " << fn << " : " << 0.9*fabs(r) << " " << fabs(l) << endl;
}while(fn > fo + c || fabs(l) > 0.9*fabs(r));
// wolfe conditions
//tst at this point is the new all vector
//now need to update hessians
......@@ -272,7 +271,8 @@ namespace NekMeshUtils
if(d3(0,0) > 2.2E-16 * ynorm)
{
B = B + y * yT * (1.0 / d1(0,0)) - B * s * sT * B * (1.0 / d2(0,0));
H = H + (d3(0,0) + n1(0,0)) / d3(0,0) / d3(0,0) * s * sT - 1.0/d3(0,0) * (H * y * sT + s * yT * H);
H = H + (d3(0,0) + n1(0,0)) / d3(0,0) / d3(0,0) * s * sT -
1.0/d3(0,0) * (H * y * sT + s * yT * H);
}
J = Jn;
......
////////////////////////////////////////////////////////////////////////////////
//
// File: SurfaceMeshing.cpp
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: surfacemeshing object methods.
//
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/Optimisation/Guass-Seidel.h>
using namespace std;
namespace Nektar
{
namespace NekMeshUtils
{
Array<OneD, NekDouble> gsOptimise(NekDouble alpha, Array<OneD, NekDouble> x, DNekMat H, DNekMat J)
{
//Array<OneD, NekDouble> eig_r(x.num_elements());
//Array<OneD, NekDouble> eig_i(x.num_elements());
//H.EigenSolve(eig_r, eig_i);
//cout << eig_r[0] << " " << eig_r[1] << " " << eig_r[2] << endl;
Array<OneD, NekDouble> dx(x.num_elements(),0.0);
NekDouble Diff;
int ct = 0;
do
{
Array<OneD, NekDouble> dxn(x.num_elements(),0.0);
for(int i = 0; i < x.num_elements(); i++)
{
NekDouble presum = 0.0;
for(int j = 0; j < i; j++)
{
presum += H(i,j)*dxn[j];
}
NekDouble postsum = 0.0;
for(int j = i + 1; j < x.num_elements(); j++)
{
postsum += H(i,j)*dx[j];
}
dxn[i] = 1.0/H(i,i) * (-J(i,0) - presum - postsum);
}
Diff = 0.0;
for(int i = 0; i < x.num_elements(); i++)
{
Diff += (dxn[i] - dx[i]) * (dxn[i] - dx[i]);
}
Diff = sqrt(Diff);
for(int i = 0; i < x.num_elements(); i++)
{
dx[i] = dxn[i];
}
if(ct>1000000)
{
cout << endl << endl << J << endl << endl << H << endl;
cout << "Failed to converge" << endl;
abort();
}
ct++;
}
while(Diff > 1E-6);
return dx;
}
}
}
////////////////////////////////////////////////////////////////////////////////
//
// File: Curvemesh.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: object for individual curve meshes.
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_MESHUTILS_OPTIMISATION_GUASS_SEIDEL_H
#define NEKTAR_MESHUTILS_OPTIMISATION_GUASS_SEIDEL_H
#include <LocalRegions/MatrixKey.h>
namespace Nektar
{
namespace NekMeshUtils
{
Array<OneD, NekDouble> gsOptimise(NekDouble alpha, Array<OneD, NekDouble> x, DNekMat H, DNekMat J);
}
}
#endif
......@@ -124,7 +124,8 @@ void InputCAD::Process()
ParseUtils::GenerateSeqVector(bl.c_str(), blsurfs);
sort(symsurfs.begin(), symsurfs.end());
sort(blsurfs.begin(), blsurfs.end());
ASSERTL0(blsurfs.size() > 0, "No surfaces selected to make boundary layer on");
ASSERTL0(blsurfs.size() > 0,
"No surfaces selected to make boundary layer on");
}
if(pSession->DefinesSolverInfo("UserDefinedSpacing"))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment