Skip to content

Commit

Permalink
Switch from global ids to point coordinates to compute orientation
Browse files Browse the repository at this point in the history 
See libMesh#3637 for a detailed discussion as to why
  • Loading branch information
@nmnobre
nmnobre committed Sep 7, 2023
1 parent 86558eb commit e0bcabc
Showing 1 changed file with 29 additions and 29 deletions.
58 changes: 29 additions & 29 deletions src/fe/fe_raviart_shape_3D.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ namespace libMesh
{


bool orientation(std::vector<dof_id_type> & arr)
bool orientation(std::vector<Point> & arr)
{
while (std::min_element(arr.begin(), arr.end()) != arr.begin())
std::rotate(arr.begin(), arr.begin() + 1, arr.end());
Expand Down Expand Up @@ -76,47 +76,47 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,
{
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 0.125*(zeta-1.0) );
else
return RealGradient( 0.0, 0.0, -0.125*(zeta-1.0) );
}
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(4), elem->node_id(0), elem->node_id(1), elem->node_id(5)};
std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
if (orientation(arr))
return RealGradient( 0.0, 0.125*(eta-1.0), 0.0 );
else
return RealGradient( 0.0, -0.125*(eta-1.0), 0.0 );
}
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(6), elem->node_id(5), elem->node_id(1), elem->node_id(2)};
std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.125*(xi+1.0), 0.0, 0.0 );
else
return RealGradient( -0.125*(xi+1.0), 0.0, 0.0 );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(7), elem->node_id(6), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 0.0, 0.125*(1.0+eta), 0.0 );
else
return RealGradient( 0.0, -0.125*(1.0+eta), 0.0 );
}
case 4:
{
std::vector<dof_id_type> arr = {elem->node_id(7), elem->node_id(3), elem->node_id(0), elem->node_id(4)};
std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
if (orientation(arr))
return RealGradient( 0.125*(xi-1.0), 0.0, 0.0 );
else
return RealGradient( -0.125*(xi-1.0), 0.0, 0.0 );
}
case 5:
{
std::vector<dof_id_type> arr = {elem->node_id(5), elem->node_id(6), elem->node_id(7), elem->node_id(4)};
std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 0.125*(1.0+zeta) );
else
Expand All @@ -141,31 +141,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,
{
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(2), elem->node_id(1)};
std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
if (orientation(arr))
return RealGradient( 2.0*xi, 2.0*eta, 2.0*zeta-2.0 );
else
return RealGradient( -2.0*xi, -2.0*eta, -2.0*zeta+2.0 );
}
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(3), elem->node_id(0)};
std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
if (orientation(arr))
return RealGradient( 2.0*xi, 2.0*eta-2.0, 2.0*zeta );
else
return RealGradient( -2.0*xi, -2.0*eta+2.0, -2.0*zeta );
}
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 2.0*xi, 2.0*eta, 2.0*zeta );
else
return RealGradient( -2.0*xi, -2.0*eta, -2.0*zeta );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 2.0*xi-2.0, 2.0*eta, 2.0*zeta );
else
Expand Down Expand Up @@ -256,15 +256,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
return RealGradient();
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(6), elem->node_id(5), elem->node_id(1), elem->node_id(2)};
std::vector<Point> arr = {elem->point(6), elem->point(5), elem->point(1), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.125, 0.0, 0.0 );
else
return RealGradient( -0.125, 0.0, 0.0 );
}
case 4:
{
std::vector<dof_id_type> arr = {elem->node_id(7), elem->node_id(3), elem->node_id(0), elem->node_id(4)};
std::vector<Point> arr = {elem->point(7), elem->point(3), elem->point(0), elem->point(4)};
if (orientation(arr))
return RealGradient( 0.125, 0.0, 0.0 );
else
Expand All @@ -288,15 +288,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
return RealGradient();
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(4), elem->node_id(0), elem->node_id(1), elem->node_id(5)};
std::vector<Point> arr = {elem->point(4), elem->point(0), elem->point(1), elem->point(5)};
if (orientation(arr))
return RealGradient( 0.0, 0.125, 0.0 );
else
return RealGradient( 0.0, -0.125, 0.0 );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(7), elem->node_id(6), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(7), elem->point(6), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 0.0, 0.125, 0.0 );
else
Expand All @@ -320,15 +320,15 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
return RealGradient();
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(1), elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 0.125 );
else
return RealGradient( 0.0, 0.0, -0.125 );
}
case 5:
{
std::vector<dof_id_type> arr = {elem->node_id(5), elem->node_id(6), elem->node_id(7), elem->node_id(4)};
std::vector<Point> arr = {elem->point(5), elem->point(6), elem->point(7), elem->point(4)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 0.125 );
else
Expand Down Expand Up @@ -360,31 +360,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
{
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(2), elem->node_id(1)};
std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
if (orientation(arr))
return RealGradient( 2.0, 0.0, 0.0 );
else
return RealGradient( -2.0, 0.0, 0.0 );
}
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(3), elem->node_id(0)};
std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
if (orientation(arr))
return RealGradient( 2.0, 0.0, 0.0 );
else
return RealGradient( -2.0, 0.0, 0.0 );
}
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 2.0, 0.0, 0.0 );
else
return RealGradient( -2.0, 0.0, 0.0 );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 2.0, 0.0, 0.0 );
else
Expand All @@ -403,31 +403,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
{
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(2), elem->node_id(1)};
std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
if (orientation(arr))
return RealGradient( 0.0, 2.0, 0.0 );
else
return RealGradient( 0.0, -2.0, 0.0 );
}
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(3), elem->node_id(0)};
std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
if (orientation(arr))
return RealGradient( 0.0, 2.0, 0.0 );
else
return RealGradient( 0.0, -2.0, 0.0 );
}
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 0.0, 2.0, 0.0 );
else
return RealGradient( 0.0, -2.0, 0.0 );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.0, 2.0, 0.0 );
else
Expand All @@ -446,31 +446,31 @@ RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
{
case 0:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(2), elem->node_id(1)};
std::vector<Point> arr = {elem->point(0), elem->point(2), elem->point(1)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 2.0 );
else
return RealGradient( 0.0, 0.0, -2.0 );
}
case 1:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(3), elem->node_id(0)};
std::vector<Point> arr = {elem->point(1), elem->point(3), elem->point(0)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 2.0 );
else
return RealGradient( 0.0, 0.0, -2.0 );
}
case 2:
{
std::vector<dof_id_type> arr = {elem->node_id(1), elem->node_id(2), elem->node_id(3)};
std::vector<Point> arr = {elem->point(1), elem->point(2), elem->point(3)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 2.0 );
else
return RealGradient( 0.0, 0.0, -2.0 );
}
case 3:
{
std::vector<dof_id_type> arr = {elem->node_id(0), elem->node_id(3), elem->node_id(2)};
std::vector<Point> arr = {elem->point(0), elem->point(3), elem->point(2)};
if (orientation(arr))
return RealGradient( 0.0, 0.0, 2.0 );
else
Expand Down

0 comments on commit e0bcabc

Please sign in to comment.