Commit ce8804d4 authored by Dave Moxey's avatar Dave Moxey

Apply clang-format to tidy things up

parent 382be47b
......@@ -296,7 +296,7 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
}
}
}
else if(geom->GetShapeType() == LibUtilities::eHexahedron)
else if (geom->GetShapeType() == LibUtilities::eHexahedron)
{
vector<Array<OneD, NekDouble> > xyz;
for (int i = 0; i < geom->GetNumVerts(); i++)
......@@ -318,53 +318,62 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
{
for (int i = 0; i < b[0]->GetNumPoints(); i++)
{
NekDouble a1 = 0.5 * (1 - eta1[i]);
NekDouble a2 = 0.5 * (1 + eta1[i]);
NekDouble b1 = 0.5 * (1 - eta2[j]),
b2 = 0.5 * (1 + eta2[j]);
NekDouble c1 = 0.5 * (1 - eta3[k]),
c2 = 0.5 * (1 + eta3[k]);
NekDouble a1 = 0.5 * (1 - eta1[i]);
NekDouble a2 = 0.5 * (1 + eta1[i]);
NekDouble b1 = 0.5 * (1 - eta2[j]),
b2 = 0.5 * (1 + eta2[j]);
NekDouble c1 = 0.5 * (1 - eta3[k]),
c2 = 0.5 * (1 + eta3[k]);
DNekMat dxdz(3, 3, 1.0, eFULL);
dxdz(0,0) = - 0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0]
+ 0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0]
- 0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0]
+ 0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
dxdz(1,0) = - 0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1]
+ 0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1]
- 0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1]
+ 0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
dxdz(2,0) = - 0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2]
+ 0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2]
- 0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2]
+ 0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
dxdz(0,1) = - 0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0]
+ 0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0]
- 0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0]
+ 0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
dxdz(1,1) = - 0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1]
+ 0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1]
- 0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1]
+ 0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
dxdz(2,1) = - 0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2]
+ 0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2]
- 0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2]
+ 0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
dxdz(0,0) = - 0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0]
- 0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0]
+ 0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0]
+ 0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
dxdz(1,0) = - 0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1]
- 0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1]
+ 0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1]
+ 0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
dxdz(2,0) = - 0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2]
- 0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2]
+ 0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2]
+ 0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
dxdz(0, 0) =
-0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0] +
0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0] -
0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0] +
0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
dxdz(1, 0) =
-0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1] +
0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1] -
0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1] +
0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
dxdz(2, 0) =
-0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2] +
0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2] -
0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2] +
0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
dxdz(0, 1) =
-0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0] +
0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0] -
0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0] +
0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
dxdz(1, 1) =
-0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1] +
0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1] -
0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1] +
0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
dxdz(2, 1) =
-0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2] +
0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2] -
0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2] +
0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
dxdz(0, 0) =
-0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0] -
0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0] +
0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0] +
0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
dxdz(1, 0) =
-0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1] -
0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1] +
0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1] +
0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
dxdz(2, 0) =
-0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2] -
0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2] +
0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2] +
0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
dxdz.Invert();
ret.push_back(dxdz);
......@@ -380,7 +389,9 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
return ret;
}
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(LocalRegions::ExpansionSharedPtr e, bool s)
Array<OneD, NekDouble> ProcessQualityMetric::GetQ(
LocalRegions::ExpansionSharedPtr e,
bool s)
{
SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
StdRegions::StdExpansionSharedPtr chi = e->GetGeom()->GetXmap();
......
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: NodalHexElec
// Description: Nodal hexahedron with 3D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
......@@ -55,7 +55,7 @@ void NodalHexElec::CalculatePoints()
for (int i = 0; i < 3; i++)
{
m_points[i] = Array<OneD, DataType>(numPoints*numPoints*numPoints);
m_points[i] = Array<OneD, DataType>(numPoints * numPoints * numPoints);
}
for (int k = 0, ct = 0; k < numPoints; k++)
......@@ -76,7 +76,7 @@ void NodalHexElec::CalculateWeights()
{
unsigned int numPoints = GetNumPoints();
m_weights = Array<OneD, DataType>(numPoints*numPoints*numPoints);
m_weights = Array<OneD, DataType>(numPoints * numPoints * numPoints);
for (int k = 0, ct = 0; k < numPoints; k++)
{
......@@ -101,6 +101,5 @@ boost::shared_ptr<PointsBaseType> NodalHexElec::Create(const PointsKey &key)
returnval->Initialize();
return returnval;
}
}
}
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: NodalHexElec
// Description: Nodal hexahedron with 3D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
......@@ -67,7 +67,10 @@ private:
{
}
Array<OneD, NekDouble> m_e0, m_ew;
/// 1D GLL points.
Array<OneD, NekDouble> m_e0;
/// 1D GLL weights.
Array<OneD, NekDouble> m_ew;
void CalculatePoints();
void CalculateWeights();
......
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: NodalQuadElec
// Description: Nodal quadrilateral with 2D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
......@@ -55,7 +55,7 @@ void NodalQuadElec::CalculatePoints()
for (int i = 0; i < 2; i++)
{
m_points[i] = Array<OneD, DataType>(numPoints*numPoints);
m_points[i] = Array<OneD, DataType>(numPoints * numPoints);
}
for (int j = 0, ct = 0; j < numPoints; j++)
......@@ -72,7 +72,7 @@ void NodalQuadElec::CalculateWeights()
{
unsigned int numPoints = GetNumPoints();
m_weights = Array<OneD, DataType>(numPoints*numPoints);
m_weights = Array<OneD, DataType>(numPoints * numPoints);
for (int j = 0, ct = 0; j < numPoints; j++)
{
......@@ -94,6 +94,5 @@ boost::shared_ptr<PointsBaseType> NodalQuadElec::Create(const PointsKey &key)
returnval->Initialize();
return returnval;
}
}
}
......@@ -29,7 +29,7 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: NodalQuadElec
// Description: Nodal quadrilateral with 2D GLL distribution
//
///////////////////////////////////////////////////////////////////////////////
......@@ -67,7 +67,10 @@ private:
{
}
Array<OneD, NekDouble> m_e0, m_ew;
/// 1D GLL points
Array<OneD, NekDouble> m_e0;
/// 1D GLL weights
Array<OneD, NekDouble> m_ew;
void CalculatePoints();
void CalculateWeights();
......
......@@ -125,7 +125,8 @@ void Mesh::MakeOrder(int order,
{
pTypes[LibUtilities::eSegment] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eTriangle] = LibUtilities::eNodalTriElec;
pTypes[LibUtilities::eQuadrilateral] = LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::eQuadrilateral] =
LibUtilities::eGaussLobattoLegendre;
pTypes[LibUtilities::ePrism] = LibUtilities::eNodalPrismElec;
pTypes[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetElec;
pTypes[LibUtilities::eHexahedron] = LibUtilities::eGaussLobattoLegendre;
......@@ -165,15 +166,15 @@ void Mesh::MakeOrder(int order,
boost::unordered_set<int> processedEdges, processedFaces, processedVolumes;
//note if CAD previously existed on the face or edge, the new points need
//to be projected onto the CAD entity.
// note if CAD previously existed on the face or edge, the new points need
// to be projected onto the CAD entity.
// Call MakeOrder with our generated geometries on each edge to fill in edge
// interior nodes.
int ct = 0;
for(eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++, ct++)
for (eit = m_edgeSet.begin(); eit != m_edgeSet.end(); eit++, ct++)
{
if(m_verbose)
if (m_verbose)
{
LibUtilities::PrintProgressbar(
ct, m_edgeSet.size(), "MakeOrder: Edges: ");
......@@ -193,9 +194,9 @@ void Mesh::MakeOrder(int order,
// Call MakeOrder with our generated geometries on each face to fill in face
// interior nodes.
ct = 0;
for(fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++, ct++)
for (fit = m_faceSet.begin(); fit != m_faceSet.end(); fit++, ct++)
{
if(m_verbose)
if (m_verbose)
{
LibUtilities::PrintProgressbar(
ct, m_faceSet.size(), "MakeOrder: Faces: ");
......@@ -251,18 +252,19 @@ void Mesh::MakeOrder(int order,
const int nElmt = m_element[m_expDim].size();
for (int i = 0; i < nElmt; ++i)
{
if(m_verbose)
if (m_verbose)
{
LibUtilities::PrintProgressbar(
i, nElmt, "MakeOrder: Elements: ");
LibUtilities::PrintProgressbar(i, nElmt, "MakeOrder: Elements: ");
}
ElementSharedPtr el = m_element[m_expDim][i];
el->MakeOrder(order, volGeoms[el->GetId()], pTypes[el->GetConf().m_e],
m_spaceDim, id);
}
if(m_verbose)
if (m_verbose)
{
cout << endl;
}
}
}
......
......@@ -150,5 +150,5 @@ IF(NEKTAR_USE_MESHGEN)
ADD_NEKTAR_TEST (MeshGen/2d-cad)
ADD_NEKTAR_TEST (MeshGen/2d-naca)
# varopti tests
ADD_NEKTAR_TEST (MeshGen/varopti_cubesphere)
ADD_NEKTAR_TEST_LENGTHY(MeshGen/varopti_cubesphere)
ENDIF()
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessJac.h
// File: Hessian.hxx
//
// 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: Calculate jacobians of elements.
// Description: Utility functions for Hessian matrices
//
////////////////////////////////////////////////////////////////////////////////
......@@ -41,45 +41,47 @@ namespace Nektar
namespace Utilities
{
//returns zero if hessian is good, 1 if indefinite
template<int DIM> int NodeOpti::IsIndefinite()
/**
* @brief Returns 1 if Hessian matrix is indefinite and 0 otherwise.
*
* Specialised versions of this function exist only for 2x2 and 3x3 matrices.
*/
template <int DIM> int NodeOpti::IsIndefinite()
{
ASSERTL0(false,"DIM error");
ASSERTL0(false, "DIM error");
return 0;
}
template<> int NodeOpti::IsIndefinite<2>()
template <> int NodeOpti::IsIndefinite<2>()
{
Array<OneD, NekDouble> eigR(2);
Array<OneD, NekDouble> eigI(2);
NekMatrix<NekDouble> H(2,2);
H(0,0) = m_grad[2];
H(1,0) = m_grad[3];
H(0,1) = H(1,0);
H(1,1) = m_grad[4];
NekMatrix<NekDouble> H(2, 2);
H(0, 0) = m_grad[2];
H(1, 0) = m_grad[3];
H(0, 1) = H(1, 0);
H(1, 1) = m_grad[4];
//cout << H << endl << endl;
// cout << H << endl << endl;
int nVel = 2;
int nVel = 2;
char jobvl = 'N', jobvr = 'N';
int worklen = 8*nVel, info;
int worklen = 8 * nVel, info;
NekDouble dum;
DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
Array<OneD, NekDouble> vl (nVel*nVel);
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi (nVel);
Array<OneD, NekDouble> wi(nVel);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel,
&(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
&dum, nVel,
&work[0], worklen, info);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
ASSERTL0(!info,"dgeev failed");
ASSERTL0(!info, "dgeev failed");
if(eval(0,0) < 0.0 || eval(1,1) < 0.0)
if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0)
{
if(eval(0,0) < 0.0 && eval(1,1) < 0.0)
if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0)
{
return 2;
}
......@@ -92,41 +94,39 @@ template<> int NodeOpti::IsIndefinite<2>()
return 0;
}
template<> int NodeOpti::IsIndefinite<3>()
template <> int NodeOpti::IsIndefinite<3>()
{
Array<OneD, NekDouble> eigR(3);
Array<OneD, NekDouble> eigI(3);
NekMatrix<NekDouble> H(3,3);
H(0,0) = m_grad[3];
H(1,0) = m_grad[4];
H(0,1) = H(1,0);
H(2,0) = m_grad[5];
H(0,2) = H(2,0);
H(1,1) = m_grad[6];
H(2,1) = m_grad[7];
H(1,2) = H(2,1);
H(2,2) = m_grad[8];
int nVel = 3;
NekMatrix<NekDouble> H(3, 3);
H(0, 0) = m_grad[3];
H(1, 0) = m_grad[4];
H(0, 1) = H(1, 0);
H(2, 0) = m_grad[5];
H(0, 2) = H(2, 0);
H(1, 1) = m_grad[6];
H(2, 1) = m_grad[7];
H(1, 2) = H(2, 1);
H(2, 2) = m_grad[8];
int nVel = 3;
char jobvl = 'N', jobvr = 'N';
int worklen = 8*nVel, info;
int worklen = 8 * nVel, info;
NekDouble dum;
DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
Array<OneD, NekDouble> vl (nVel*nVel);
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi (nVel);
Array<OneD, NekDouble> wi(nVel);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel,
&(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
&dum, nVel,
&work[0], worklen, info);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
ASSERTL0(!info,"dgeev failed");
ASSERTL0(!info, "dgeev failed");
if(eval(0,0) < 0.0 || eval(1,1) < 0.0 || eval(2,2) < 0.0)
if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0 || eval(2, 2) < 0.0)
{
if(eval(0,0) < 0.0 && eval(1,1) < 0.0 && eval(2,2))
if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0 && eval(2, 2))
{
return 2;
}
......@@ -139,96 +139,99 @@ template<> int NodeOpti::IsIndefinite<3>()
return 0;
}
template<int DIM> void NodeOpti::MinEigen(NekDouble &val)
/**
* @brief Calculates minimum eigenvalue of Hessian matrix.
*
* Specialised versions of this function exist only for 2x2 and 3x3 matrices.
*/
template <int DIM> void NodeOpti::MinEigen(NekDouble &val)
{
ASSERTL0(false,"DIM error");
ASSERTL0(false, "DIM error");
}
template<> void NodeOpti::MinEigen<2>(NekDouble &val)
template <> void NodeOpti::MinEigen<2>(NekDouble &val)
{
Array<OneD, NekDouble> eigR(2);
Array<OneD, NekDouble> eigI(2);
NekMatrix<NekDouble> H(2,2);
H(0,0) = m_grad[2];
H(1,0) = m_grad[3];
H(0,1) = H(1,0);
H(1,1) = m_grad[4];
NekMatrix<NekDouble> H(2, 2);
H(0, 0) = m_grad[2];
H(1, 0) = m_grad[3];
H(0, 1) = H(1, 0);
H(1, 1) = m_grad[4];
int nVel = 2;
int nVel = 2;
char jobvl = 'N', jobvr = 'V';
int worklen = 8*nVel, info;
int worklen = 8 * nVel, info;
DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec (nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl (nVel*nVel);
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec(nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi (nVel);
Array<OneD, NekDouble> wi(nVel);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel,
&(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
&(evec.GetPtr())[0], nVel,
&work[0], worklen, info);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &(evec.GetPtr())[0], nVel, &work[0],
worklen, info);
ASSERTL0(!info,"dgeev failed");
ASSERTL0(!info, "dgeev failed");
int minI;
NekDouble tmp = std::numeric_limits<double>::max();
for(int i = 0; i < 2; i++)
for (int i = 0; i < 2; i++)
{
if(eval(i,i) < tmp)
if (eval(i, i) < tmp)
{
minI = i;
tmp = eval(i,i);
tmp = eval(i, i);
}
}
val = eval(minI,minI);
val = eval(minI, minI);
}
template<> void NodeOpti::MinEigen<3>(NekDouble &val)
template <> void NodeOpti::MinEigen<3>(NekDouble &val)
{
Array<OneD, NekDouble> eigR(3);
Array<OneD, NekDouble> eigI(3);
NekMatrix<NekDouble> H(3,3);
H(0,0) = m_grad[3];
H(1,0) = m_grad[4];
H(0,1) = H(1,0);
H(2,0) = m_grad[5];
H(0,2) = H(2,0);
H(1,1) = m_grad[6];
H(2,1) = m_grad[7];
H(1,2) = H(2,1);
H(2,2) = m_grad[8];
int nVel = 3;
NekMatrix<NekDouble> H(3, 3);
H(0, 0) = m_grad[3];
H(1, 0) = m_grad[4];
H(0, 1) = H(1, 0);
H(2, 0) = m_grad[5];
H(0, 2) = H(2, 0);
H(1, 1) = m_grad[6];
H(2, 1) = m_grad[7];
H(1, 2) = H(2, 1);
H(2, 2) = m_grad[8];
int nVel = 3;
char jobvl = 'N', jobvr = 'V';
int worklen = 8*nVel, info;
int worklen = 8 * nVel, info;
DNekMat eval (nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec (nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl (nVel*nVel);
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec(nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi (nVel);
Array<OneD, NekDouble> wi(nVel);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel,
&(eval.GetPtr())[0], &wi[0], &vl[0], nVel,
&(evec.GetPtr())[0], nVel,
&work[0], worklen, info);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &(evec.GetPtr())[0], nVel, &work[0],
worklen, info);
ASSERTL0(!info,"dgeev failed");
ASSERTL0(!info, "dgeev failed");
int minI;
NekDouble tmp = std::numeric_limits<double>::max();
for(int i = 0; i < 3; i++)
for (int i = 0; i < 3; i++)
{
if(eval(i,i) < tmp)
if (eval(i, i) < tmp)
{
minI = i;
tmp = eval(i,i);
tmp = eval(i, i);
}
}
val = eval(minI,minI);
val = eval(minI, minI);
}
}
......
......@@ -55,34 +55,19 @@ NodeOptiFactory &GetNodeOptiFactory()
return asd;
}
const NekDouble NodeOpti::gam = numeric_limits<float>::epsilon()*10;
void NodeOpti::CalcMinJac()
{
m_minJac = numeric_limits<double>::max();
map<LibUtilities::ShapeType,vector<ElUtilSharedPtr> >::iterator typeIt;
for(typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
map<LibUtilities::ShapeType, vector<ElUtilSharedPtr> >::iterator typeIt;
for (typeIt = m_data.begin(); typeIt != m_data.end(); typeIt++)
{
for(int i = 0; i < typeIt->second.size(); i++)
for (int i = 0; i < typeIt->second.size(); i++)
{
m_minJac = min(m_minJac,typeIt->second[i]->GetMinJac());
m_minJac = min(m_minJac, typeIt->second[i]->GetMinJac());
}
}
}
NekDouble dir[12][3] = {{1,0,0},
{-1,0,0},
{0,1,0},
{0,-1,0},
{0,0,1},
{0,0,-1},