Commit 2fc453a4 authored by Michael Turner's avatar Michael Turner

fix triangle face handling in face orientation of 3D elements

parent 15253780
......@@ -686,11 +686,41 @@ namespace Nektar
// belongs...
// Compute the length of edges on a base-face
if( f==1 || f==3 ) { // Face is a Triangle
for(i = 0; i < m_coordim; i++)
if( f==1 || f==3 )
{ // Face is a Triangle
if (baseVertex == m_verts[faceVerts[f][0]]->GetVid())
{
elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][0]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][0]])[i];
}
}
else if (baseVertex == m_verts[faceVerts[f][1]]->GetVid())
{
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][0]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][1]])[i];
}
}
else if (baseVertex == m_verts[faceVerts[f][2]]->GetVid())
{
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][2]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][0]])[i];
}
}
else
{
ASSERTL0(false, "Could not find matching vertex for the face");
}
}
else { // Face is a Quad
......@@ -779,12 +809,9 @@ namespace Nektar
norm = fabs(dotproduct2 / elementBaxis_length / faceBaxis_length);
// ASSERTL1(fabs(norm - 1.0) <
// NekConstants::kNekZeroTol,
// "These vectors should be parallel");
//This assert is commented out because tests fail, it has been for years
//there is a bug further up the code as this should not fail
//unsure what it is (MT 2017)
ASSERTL1(fabs(norm - 1.0) <
NekConstants::kNekZeroTol,
"These vectors should be parallel");
// if the inner product is negative, both B-axis point
// in reverse direction
......
......@@ -486,11 +486,41 @@ namespace Nektar
// Compute the length of edges on a base-face
if (f > 0)
{
for(i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[ faceVerts[f][1] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
elementBaxis[i] = (*m_verts[ faceVerts[f][2] ])[i] - (*m_verts[ faceVerts[f][0] ])[i];
}
if (baseVertex == m_verts[faceVerts[f][0]]->GetVid())
{
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][0]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][0]])[i];
}
}
else if (baseVertex == m_verts[faceVerts[f][1]]->GetVid())
{
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][0]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][1]])[i];
}
}
else if (baseVertex == m_verts[faceVerts[f][2]]->GetVid())
{
for (i = 0; i < m_coordim; i++)
{
elementAaxis[i] = (*m_verts[faceVerts[f][1]])[i] -
(*m_verts[faceVerts[f][2]])[i];
elementBaxis[i] = (*m_verts[faceVerts[f][2]])[i] -
(*m_verts[faceVerts[f][0]])[i];
}
}
else
{
ASSERTL0(false, "Could not find matching vertex for the face");
}
}
else
{
......@@ -581,12 +611,9 @@ namespace Nektar
norm = fabs(dotproduct2) / elementBaxis_length / faceBaxis_length;
// check that both these axis are indeed parallel
// ASSERTL1(fabs(norm - 1.0) <
// NekConstants::kNekZeroTol,
// "These vectors should be parallel");
//This assert is commented out because tests fail, it has been for years
//there is a bug further up the code as this should not fail
//unsure what it is (MT 2017)
ASSERTL1(fabs(norm - 1.0) <
NekConstants::kNekZeroTol,
"These vectors should be parallel");
// if the inner product is negative, both B-axis point
// in reverse direction
......
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