Skip to content

Commit

Permalink
Refactored to eliminate the need for the ESEAS-ordering lookup methods.
Browse files Browse the repository at this point in the history
  • Loading branch information
@CamelliaDPG
CamelliaDPG committed Apr 2, 2020
1 parent fb98ff9 commit c8c0461
Showing 1 changed file with 15 additions and 101 deletions.
116 changes: 15 additions & 101 deletions packages/intrepid2/src/Discretization/Basis/Intrepid2_IntegratedLegendreBasis_HGRAD_TET.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,88 +59,6 @@

namespace Intrepid2
{
KOKKOS_INLINE_FUNCTION
int localTetFaceBasisOrdinalESEAS(const int polyOrder, const int i, const int j)
{
/*
ESEAS numbers the (i,j)'s here somewhat differently from the loops
we have in the face basis computations below. Because we want to
reuse the jacobi values for each i value, we prefer our loop structure,
but we want to number the basis the same way ESEAS does. This lambda
allows that.
For faster execution, we could derive a formula in terms of i,j.
This is not likely to be a performance bottleneck, however, and we prefer
the clarity of making the loop ordering explicit.
The loop ordering below is derived from the code in shape3DHTet in ESEAS.
*/
int localFaceBasisOrdinal = 0;
const int min_j = 1;
const int min_i = 2;
const int min_ij = min_i + min_j;
const int max_ij = polyOrder;
for (int totalPolyOrder=min_ij; totalPolyOrder <= max_ij; totalPolyOrder++)
{
for (int ii=min_i; ii<=totalPolyOrder-min_j; ii++)
{
const int jj = totalPolyOrder - i;

if ((i==ii) && (j==jj))
{
return localFaceBasisOrdinal;
}

localFaceBasisOrdinal++;
}
}
return -1;
}

KOKKOS_INLINE_FUNCTION
int localTetInteriorBasisOrdinalESEAS(const int polyOrder, const int i, const int j, const int k)
{
/*
ESEAS numbers the (i,j,k)'s here somewhat differently from the loops
we have in the interior basis computations below. Because we want to
reuse the jacobi values for each (i+j) value, we prefer our loop structure,
but we want to number the basis the same way ESEAS does. This lambda
allows that.
For faster execution, we could derive a formula in terms of i,j.
This is not likely to be a performance bottleneck, however, and we prefer
the clarity of making the loop ordering explicit.
The loop ordering below is derived from the code in shape3DHTet in ESEAS.
*/
int localInteriorBasisOrdinal = 0;

const int min_i = 2;
const int min_j = 1;
const int min_k = 1;
const int min_ij = min_i + min_j;
const int min_ijk = min_ij + min_k;
const int max_ijk = polyOrder;

for (int totalPolyOrder_ijk=min_ijk; totalPolyOrder_ijk <= max_ijk; totalPolyOrder_ijk++)
{
for (int totalPolyOrder_ij=min_ij; totalPolyOrder_ij <= totalPolyOrder_ijk-min_k; totalPolyOrder_ij++)
{
for (int ii=min_i; ii <= totalPolyOrder_ij-min_j; ii++)
{
const int jj = totalPolyOrder_ij - i;
const int kk = totalPolyOrder_ijk - totalPolyOrder_ij;
if ((i==ii) && (j==jj) && (k==kk))
{
return localInteriorBasisOrdinal;
}
localInteriorBasisOrdinal++;
}
}
}
return -1;
}

/** \class Intrepid2::Hierarchical_HGRAD_TET_Functor
\brief Functor for computing values for the IntegratedLegendreBasis_HGRAD_TET class.
Expand Down Expand Up @@ -757,43 +675,39 @@ namespace Intrepid2

// **** face functions **** //
const int numFunctionsPerFace = ((polyOrder-1)*(polyOrder-2))/2;
const int max_ij_sum = polyOrder;
const int numFaces = 4;
for (int faceOrdinal=0; faceOrdinal<numFaces; faceOrdinal++)
{
for (int i=2; i<max_ij_sum; i++)
for (int totalPolyOrder=3; totalPolyOrder<=polyOrder_; totalPolyOrder++)
{
for (int j=1; i+j<=max_ij_sum; j++)
const int totalFaceDofs = (totalPolyOrder-2)*(totalPolyOrder-1)/2;
const int totalFaceDofsPrevious = (totalPolyOrder-3)*(totalPolyOrder-2)/2;
const int faceDofsForPolyOrder = totalFaceDofs - totalFaceDofsPrevious;
for (int i=0; i<faceDofsForPolyOrder; i++)
{
const int faceLocalBasisOrdinal = localTetFaceBasisOrdinalESEAS(polyOrder_, i, j);
const int fieldOrdinal = fieldOrdinalOffset + faceLocalBasisOrdinal;
this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = i+j;

TEUCHOS_TEST_FOR_EXCEPTION(faceLocalBasisOrdinal < 0, std::invalid_argument, "Internal error: negative faceLocalBasisOrdinal");
this->fieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder;
fieldOrdinalOffset++;
}
}
fieldOrdinalOffset += numFunctionsPerFace;
}

// **** interior functions **** //
const int numFunctionsPerVolume = ((polyOrder-1)*(polyOrder-2)*(polyOrder-3))/6;
const int max_ijk_sum = polyOrder;
const int numVolumes = 1; // interior
for (int volumeOrdinal=0; volumeOrdinal<numVolumes; volumeOrdinal++)
{
for (int i=2; i<=max_ijk_sum-2; i++)
for (int totalPolyOrder=4; totalPolyOrder<=polyOrder_; totalPolyOrder++)
{
for (int j=1; i+j<=max_ijk_sum-1; j++)
const int totalInteriorDofs = (totalPolyOrder-3)*(totalPolyOrder-2)*(totalPolyOrder-1)/6;
const int totalInteriorDofsPrevious = (totalPolyOrder-4)*(totalPolyOrder-3)*(totalPolyOrder-2)/6;
const int interiorDofsForPolyOrder = totalInteriorDofs - totalInteriorDofsPrevious;

for (int i=0; i<interiorDofsForPolyOrder; i++)
{
for (int k=1; i+j+k<=max_ijk_sum; k++)
{
const int interiorLocalBasisOrdinal = localTetInteriorBasisOrdinalESEAS(polyOrder_, i, j, k);
const int fieldOrdinal = fieldOrdinalOffset + interiorLocalBasisOrdinal;
this->fieldOrdinalPolynomialDegree_(fieldOrdinal,0) = i+j+k;
}
this->fieldOrdinalPolynomialDegree_(fieldOrdinalOffset,0) = totalPolyOrder;
fieldOrdinalOffset++;
}
}
fieldOrdinalOffset += numFunctionsPerVolume;
}

INTREPID2_TEST_FOR_EXCEPTION(fieldOrdinalOffset != this->basisCardinality_, std::invalid_argument, "Internal error: basis enumeration is incorrect");
Expand Down

0 comments on commit c8c0461

Please sign in to comment.