# Finding paths between triangles efficiently in 3D geometry

An update to this post can be found here

I've been writing some functions used to find paths between two types of triangles - alphas and betas. Alphas are triangles that have been in a zone we consider important, have an "interesting" value above a given threshold, and are "active". Betas are essentially anything that isn't an Alpha.

The position of the zone and the geometry of the model can change between invocations.

This is written in C++03, compiled into a MEX file (.mexw64) to be executed by MATLAB R2016B on a Linux machine. Those are all hard limits.

This code uses a good deal of functions and data from an external libraries and objects. However, most of the methods used are very simple array lookups, nothing performance-hindering.

Everything works correctly so far in testing, but performance has become a significant problem.

The code:

// Doxygen block exists here

// Various includes go here

// Only needed because ultimately MATLAB needs an error code, not a C++
// exception
#define SUCCESS 0
#define DYN_ALLOC_ERR 1

/*
* Design notes: I considered having a modified version of the geometry checking
* section that just pushed _new_ alphas onto the global vector of alphas
* and set a flag to cull old ones, but that actually seemed significantly less
* efficient than just making a new one each time since all the same checks need
* to be made regardless however, if validAlphas gets very large, this could
* be highly inefficient due to push_back, so the idea of pushing on
* just new ones and culling may actually be the best solution.
*
* Also: Could maintain a sum of validAlphas and use that to re-size it upon
* re-creation to save some overhead from resizing, maybe?
*/

// The indices of two triangles who have a valid alpha and a beta that can
// be seen from it. 120k paths * 8 bytes per pair = ~1GB on the heap.
// No, ushort won't be enough here.
struct ABPair
{
unsigned int alphaTriIndex;
unsigned int betaTriIndex;
};

// Useful for multithreading, stolen and modified from craytracer.h
{
CRayTracer* rt;
unsigned uThreadID;                     // Index
std::vector<ABPair*> validPathsThread;  // valid pairs that this thread
// found
unsigned int numTris;                   // Input argument, the number of
// triangles in the mesh
double distThreshold;                   // Input argument, the maximum
// distance between triangles
};

// Exceptions for experimentation
class PThreadException: public std::exception
{
virtual const char* what() const throw()
{
return "Exception occured in a pthread_attr_init or pthread_create\n";
}
};

class DynAllocationException: public std::exception
{
virtual const char* what() const throw()
{
return "Exception occured when attempting to malloc or calloc\n";
}
};

// Note: Globals must exist here so that when the MEX file exits and goes
// Back to MATLAB, this information is maintained.
// (AFAIK, mexMakeMemoryPersistant() wouldn't make sense here. I can link
// to discussions on that topic)

// An indicator for every tri to tell if it has been removed (neccesary for
// maintaining previous state to check if we need to call findPaths())
static bool* triActive = NULL;

// A map from a given face to the element it resides in
static unsigned int* faceToElementMap = NULL;

// All valid paths. Must be maintained, because we don't know if
// findPaths() will be called. It may not be if geometry hasnt changed.
static std::vector<ABPair*> validPaths;

// The previous state of what alphas were considered valid. Neccesary to see
// if a change has occured that isn't geometry-based.
static unsigned int* validAlphaIndex;

// Not neccesary as a global if it's just getting re-generated each time.
// However, if we do decide to append and cull instead of regenerating, this
// needs to stay.
static std::vector<unsigned int> validAlphas;

// Needed so we can accurately determine alphas.
// I removed this in the past, thinking it wasn't needed. As it turns out,
// it's absolutely needed and very helpful.
static bool* hasBeenInZone;

// Useful everywhere
CRayTracerClass* rayTracer = NULL;
NanoRTWrapperClass nanoRTWrapper = NULL;

// Function declarations
// Not required, but prevents warnings depending on how functions are ordered
// and call each other
// (Including the mexFunction here would be redundant, as it exists in mex.h)
void exitFcn();
bool isTriInZoneRadius(const unsigned int itri);
bool checkForModelChanges(const unsigned int numTris,
const float* iValues,
const double iThreshold
);
void initialize(const float* elemFace,
const unsigned int numElems,
const unsigned int facePerElMax,
unsigned int* numTri,
unsigned int* numFace
);
void findPathsThreadSpooler(const unsigned int numTris,
const double distThreshold
);
void mapFacesToElements(const float* elemFace,
const unsigned int numElems,
const unsigned int facePerElMax
);
bool checkPairValid(const unsigned int i,
const unsigned int j,
const double distThreshold
);
bool isTriAlpha(const unsigned int itri,
const float* iValues,
const double iThreshold
);
void findPaths(const unsigned int numTris,
const double distThreshold
);
//mainfunc declaration goes here

/**
* @brief exitFcn - Cleans up malloc'd or calloc'd memory if someone in the
* MATLAB script calls "clear mexfilename" or "clear all".
*/
void exitFcn()
{
//mexPrintf("exitFcn() called\n");

if(triActive)
{
free(triActive);
}
if(faceToElementMap)
{
free(faceToElementMap);
}
if(validAlphaIndex)
{
free(validAlphaIndex);
}
if(hasBeenInZone)
{
free(hasBeenInZone);
}
for(unsigned int i = 0; i < validPaths.size(); i++)
{
free(validPaths[i]);
}
}

/**
* @brief Checks if a given tri is currently in the zone's external radius.
* Implementation stolen from CRayTracerClass::trace_inverseray_tri
* Not sure if we need to raytrace, so I omitted it
* @param itri - The index of the triangle to check
* @return True if in the radius, false if not
*/
bool isTriInZoneRadius(const unsigned int itri)
{
//Omitted
}

/**
* @brief Checks if the model has changed (either in terms of alphas or
* geometry) and re-generates the vector of alphas
* @param numTris -     The number of triangles in the finite mesh
* @param iValues -     The ivalue at each node
* @param iThreshold -  The interesting value threshold beyond which an alpha
* is interesting enough to be valid
* @return True if the list of alphas or the geometry has changed, false if
* neither have
*/
bool checkForModelChanges(const unsigned int numTris,
const float* iValues,
const double iThreshold
)
{
bool modelChanged = false;
bool isAlpha;
bool currentlyActive;

// Two checks need to happen - geometry changes and for the list of valid
// alphas to change
// Also regenerates the vector of valid alphas from scratch as it goes

for(unsigned int i = 0; i < numTris; i++)
{
// Active means it has 1 exposed face, not 2 (internal) or 0 (gone)
currentlyActive = nanoRTWrapper->getTriActive(i);

// Has the geometry changed?
if(currentlyActive != triActive[i])
{
modelChanged = true;
triActive[i] = currentlyActive;
}

// Get whether this triangle is an alpha:
isAlpha = isTriAlpha(i, iValues, iThreshold);

// Triangle is a valid alpha now, but wasn't before
if((isAlpha == true) && (validAlphaIndex[i] == false))
{
validAlphaIndex[i] = true;
modelChanged = true;
}
// Was valid before, is no longer valid now
else if((isAlpha == false) && (validAlphaIndex[i] == true))
{
validAlphaIndex[i] = false;
modelChanged = true;
//cullalphasFlag = true;
}

// Generating the set of all valid alphas
if(isAlpha)
{
validAlphas.push_back(i);
}
}

return modelChanged;
}

/**
* @brief Initializes this MEX file for its first run
* @param rt -            A pointer to the raytracer object
* @param numTris -       The total number of triangles in the finite mesh
* @param numFaces -      The total number of faces in the finite mesh
* @param elemFace -      The map of elements to the faces that they have
* @param numElems -      The number of elements in the finite mesh
* @param facePerElMax -  The maximum number of faces per element
*/
void initialize(const float* elemFace,
const unsigned int numElems,
const unsigned int facePerElMax,
unsigned int* numTri,
unsigned int* numFace
)
{
DynAllocationException e;

// Fetch number of tris and faces
// Must be done every time, since we're storing locally and not globally
// However:
// They're never modified
// They never change between calls from the MATLAB script
// They're used frequently in many functions
// I think that's a strong candidate for being a global

unsigned int numTris = nanoRTWrapper->getTriCount();
*numTri = numTris;

unsigned int numFaces = nanoRTWrapper->getFaceCount();
*numFace = numFaces;

/*
* Allocate some space for things we need to be persistent between runs of
* this MEX file. And check that the allocation succeeded, of course.
*/
if(triActive == NULL)
{
if(NULL ==
(triActive =
static_cast<bool*>(calloc(numTris, sizeof(bool))))
)
{
throw e;
}
}
if(hasBeenInZone == NULL)
{
if(NULL ==
(hasBeenInZone =
static_cast<bool*>(calloc(numTris, sizeof(bool))))
)
{
throw e;
}
}
if(validAlphaIndex == NULL)
{
if(NULL ==
(validAlphaIndex =
static_cast<unsigned int*>(calloc(numTris, sizeof(unsigned int))))
)
{
throw e;
}
}
if(faceToElementMap == NULL)
{
if(NULL ==
(faceToElementMap =
static_cast<unsigned int*>(calloc(numFaces, sizeof(unsigned int))))
)
{
throw e;
}
mapFacesToElements(elemFace, numElems, facePerElMax);
}

return;
}

/**
* @brief Is something that can be used by pthread_create(). Threads will skip
* over some of the work, and do isValidPair on others. Thus...multithreading.
* @param data - The data structure that will hold the results and arguments
*/
{
const double distThreshold = thisThreadsData->distThreshold;
const unsigned int numTris = thisThreadsData->numTris;

// Loop over all valid alphas
for(unsigned int i = 0; i < validAlphas.size(); i++)
{
{
// Loop over all triangles (potential betas)
for(unsigned int j = 0; j < numTris; j++)
{
if(checkPairValid(i, j, distThreshold))
{
ABPair* temp =
static_cast<ABPair*>(malloc(sizeof(ABPair)));

temp->alphaTriIndex = validAlphas[i];
temp->betaTriIndex = j;

}
}
}
}
return NULL;
}

/**
* @brief Creates as many threads as the system has available, and then uses
* pthread_create() to dish out the work of findPaths()
* @param numTris - The number of triangles in the finite mesh
* @param distThreshold - The maximum distance an alpha and beta can be
* apart
*/
void findPathsThreadSpooler(const unsigned int numTris,
const double distThreshold
)
{
int rc;

// I think this is checking to make sure something doesn't already exist,
// not sure what though
{
throw e;
}

// We know how many threads the system supports
// So all this does is walk through an array of them and start them up
{

data.rt = rayTracer;
data.numTris = numTris;
data.distThreshold = distThreshold;

{
{
throw e;
}
}
else
{
}
}

// Join all threads
{

{
void* res;

{
throw e;
}
}

// validPathsThread is the set of ABPairs that this thread found
// while validPaths is the globally maintained set of valid paths
// Take each thread's results and merge it into the overall results
validPaths.insert(validPaths.end(),
}

return;
}

/*
void cullAlphas()
{
for(unsigned int i = 0; i < validAlphas.size(); i++)
{
if(!isValidAlpha(validAlphas[i]))
{
validAlphas.erase(i);
}
}
}
*/

/**
* @brief Determines the elements that each face belongs to
* @details the MATLAB script maintains a map of all faces per element.
* This is the opposite of what we want. Accessing it linearly
* walks by column, not by row. Also, MATLAB stores everything 1-indexed.
* Finally, the MATLAB script left them stored as the default, which are singles.
* @param elemFace - A MATLAB facePerElMax by numElems array, storing which
* faces belong to each element (elements being the row number)
* @param numElems - The total number of elements (rows) in the array
* @param facePerElMax - The max number of faces per element (the number of
* columns)
*/
void mapFacesToElements(const float* elemFace,
const unsigned int numElems,
const unsigned int facePerElMax
)
{
unsigned int i;

// elemFace[0] = 1. We don't know how elemFace will be structured precisely,
// so we need to keep going until we find a face in it that equals our number
// of faces, since it's 1-indexed.
for(i = 0; i < (numElems * facePerElMax); i++)
{
faceToElementMap[static_cast<unsigned int>(elemFace[i]) - 1] =
(i % numElems);

// Is the next face for that element a NaN? If so, we can skip it. Keep
// skipping until the next element WON'T be NaN.
// Don't cast here, as NaN only exists for floating point numbers,
// not integers.
while(isnan(elemFace[i + 1]) && ((i + 1) < (numElems * facePerElMax)))
{
i++;
}
}
}

/**
* @brief checkPairValid - Checks if a pair of an alpha index (of validAlphas),
* beta index form a valid path
* @param i -             Index into validAlphas
* @param j -             Index into all tris (potential beta)
* @param distThreshold - The max distance the tri's centers can be apart
* @return Whether the pair forms a valid path
*/
bool checkPairValid(const unsigned int i,
const unsigned int j,
const double distThreshold
)
{
double path_dist_sqrd;
double path_dist;
double alphaCoords[3];
double betaCoords[3];
nanort::Ray<double> ray;

// If they're not an alpha currently, they must be a potential beta,
// must also be alive
if(!validAlphaIndex[j] && triActive[j])
{
alphaCoords[0] =
rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.x();
alphaCoords[1] =
rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.y();
alphaCoords[2] =
rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.z();

betaCoords[0] = rayTracer->m_vecTriFixedInfo[j].center.x();
betaCoords[1] = rayTracer->m_vecTriFixedInfo[j].center.y();
betaCoords[2] = rayTracer->m_vecTriFixedInfo[j].center.z();

// Determine distance squared between alpha and beta
// (x2-x1)^2 + (y2-y1)^2 +(z2-z1)^2
path_dist_sqrd = pow((betaCoords[0] - alphaCoords[0]), 2)
+ pow((betaCoords[1] - alphaCoords[1]), 2)
+ pow((betaCoords[2] - alphaCoords[2]), 2);

// Doing this instead of doing the sqrt to save doing the sqrt when not
// needed for performance
if(path_dist_sqrd <= pow(distThreshold, 2))
{
path_dist = sqrt(path_dist_sqrd);

// Set up a nanort::Ray's origin, direction, and max distance
ray.org[0] = alphaCoords[0]; // x
ray.org[1] = alphaCoords[1]; // y
ray.org[2] = alphaCoords[2]; // z

ray.dir[0] = (betaCoords[0] - alphaCoords[0]) / path_dist;
ray.dir[1] = (betaCoords[1] - alphaCoords[1]) / path_dist;
ray.dir[2] = (betaCoords[2] - alphaCoords[2]) / path_dist;

// TODO: Subtract some EPSILON here so it doesn't report a hit because it
// hit the beta itself (assuming that's how it works)
ray.max_t = path_dist;

// Call ShootRay() to check if there is a path (calls nanoRT)
if(!(nanoRTWrapper->shootRay(ray)))
{
return true;
}
else
{
// There's no path
return false;
}
}
else
{
// The distance is too far between alpha and beta
return false;
}
}
else
{
// The beta is either dead or currently an alpha
return false;
}
}

/**
* @brief Determines if a given triangle is a valid alpha.
* @param itri - The triangle index to check
* @return True if it is an alpha, false if it is not
*/
bool isTriAlpha(const unsigned int itri,
const float* iValues,
const double iThreshold
)
{
double tri_avg_interesting;
const unsigned int* tri_nodes;

// Do the simple checks first, as it's more performant to do so
// alternate consideration (acccuracy, wouldn't affect performance)
//if(triActive[itri] && (hasBeenAlpha[itri] || isTriInZoneRadius(itri)))
if(triActive[itri] && (hasBeenInZone[itri] || isTriInZoneRadius(itri)))
{
// Retrieve the average iValue of this triangle
tri_nodes = nanoRTWrapper->getTriNodes(itri);

tri_avg_interesting = (iValues[tri_nodes[0]]
+ iValues[tri_nodes[1]]
+ iValues[tri_nodes[2]]) / 3;

if(tri_avg_interesting > iThreshold)
{
return true;
}
}

return false;
}

/**
* @brief Uses the raytracer object's current state as well as arguments to
* generate pairs of unobstructed paths between alphas and betas.
* @param numTris - The number of triangles in the finite element mesh
* @param distThreshold - The max distance an alpha and beta pair can be
* apart before not being considered in calculations
*/
void findPaths(const unsigned int numTris,
const double distThreshold
)
{
// This function once held more importance, but yes, it could be omitted
// at this point.

// Spool up some threads to take care of the work
try
{
}
catch(DynAllocationException& e)
{
throw e;
}

return;
}

// Doxygen header (omitted)
int mainFunc(args)
{
// Initialize the program if we're on a first run
try
{
initialize(elemFace, numElems, facePerElMax, &numTris, &numFaces);
}
catch(DynAllocationException& e)
{
return DYN_ALLOC_ERR;
}

// Need to see if we need to call findPaths
if(checkForModelChanges(numTris, iValues, iThreshold))
{
// Remove old list of valid paths
for(unsigned int i = 0; i < validPaths.size(); i++)
{
free(validPaths[i]);
}
validPaths.clear();

try
{
findPaths(numTris,
distThreshold);
}
{
}
}

//mexPrintf("Number of valid paths: %d\n", validPaths.size());

/*
* Walk over all paths and do some calculations
*/

return SUCCESS;
}

// Doxygen block goes here, omitted
void mexFunction(int nlhs,
mxArray *plhs[],
int nrhs,
const mxArray *prhs[]
)
{
// Register exit function
mexAtExit(exitFcn);

// Prep for writing out results

// Checking to make sure # of arguments was right from MATLAB

// Input argument handling to convert from mxArrays to double*, float*, etc

//*errcode = mainFunc(some args)

// retrieve execution time in clock cycles, convert to seconds, print

// Put the outputs in plhs
}


Callgraph(?):

This isn't exactly a callgraph, but it might be useful to get an idea of the flow of the program.

The Problem: Performance

For medium-size models (104k tris, 204k faces, 51k elems) it can take up to a couple seconds for this to complete, even though the worst of it is multi-threaded on a powerful 4C/8T machine. (roughly 100*104k size loop)

For any models where the number of alphas is very large (50K) it can take up to three minutes for a single execution to complete because of how large that double-nested for loop must become. (50k^2 size loop)

Possible optimizations:

An optimization that may be worth considering is creating a sphere around all alphas, and use its center and radius to cull the remaining triangles down to a smaller size based on the threshold distance. However, the benefit of this is extremely variable and may actually be zero for smaller meshes. And to create the sphere, we need to find the two triangle centers that are the furthest apart, an O(alphas^2) operation...very slow if there are many. And there's no easy way to tell in advance if using this technique would be beneficial for the given state of the mesh, so it's not easy to toggle on and off as needed.

As well, it's possible something can be done with nanoRT's BVHs, but I'm not very familiar with BVH's or what they'd let me do in this

Note: How it's being used:

The MATLAB script will likely call this many times. In small models, it may finish its own loop within tenths of a second and then call ours again. In larger ones, there may be half a second between calls. In total, this may be called hundreds of times.

Note: How it's being built:

This isn't built using the MEX command in MATLAB, nor by using Visual Studio. Instead, g++ is used to create an object file (.o) and then g++ is used again to create the .mexw64 file in a method I'm not entirely familiar with. (This is also a hard limit I cannot touch)

I occasionally compiled with very aggressive warnings enabled to catch things like sign conversion, promotions, bad casts, etc.

Profiling:

I would love to be able to profile this code. However, it seems impossible. MEX files built using MEX command in MATLAB can be done. MEX files compiled in Visual Studio can be profiled. But we're not doing either of those, and so when I try to profile with either MATLAB or Visual Studio, it just doesn't work.

Even if I could, I don't think it would reveal anything surprising. The numbers we're working with are large, so the double-nested loop at the core of it grows very large.

If need be, I could stick a ton of clock() calls around the start/end of every method to get some more info, though that's not precise. (Line-by-line)

Final Notes:

I'm fresh out of college and this is my first real production code. I'm more familiar with C than C++, so I'm sure that bled into the code. Suggestions for style are always welcome, though performance is the most major concern here.

Update:

Thanks to all for the help. I wish I could mark all 3 answers as solutions, as each were helpful. In total, it's managed to shave 8.2% off of the total runtime. And it's just generally less C-like and more C++ like. I may post an updated question now that this has been taken care of.

• Is C++11 not usable for you for some reason? – Edward Jul 3 '19 at 19:35
• @Edward Yes, due the environment this will ultimately be run on C++11 is impossible. – Tyler Shellberg Jul 3 '19 at 19:37
• Why is mexFunction mostly empty? It seems your MEX-file doesn’t do anything at all except register the at-exit function. – Cris Luengo Jul 4 '19 at 2:09
• @CrisLuengo I omitted that code, as it wasn't relevant to the performance problem and would've just cluttered the post. I just wanted to show that it existed and what it generally did with the comments. – Tyler Shellberg Jul 4 '19 at 2:14

  const double distThreshold = thisThreadsData->distThreshold;

for(unsigned int i = 0; i < validAlphas.size(); i++)
{
// Loop over all triangles (potential betas)
for(unsigned int j = 0; j < numTris; j++)
{
if(checkPairValid(i, j, distThreshold))


calls checkPairValid():

bool checkPairValid(const unsigned int i,
const unsigned int j,
const double distThreshold
)
{
double path_dist_sqrd;

if(!validAlphaIndex[j] && triActive[j])
{
alphaCoords[0] = rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.x();
//...

// Determine distance squared between alpha and beta
// (x2-x1)^2 + (y2-y1)^2 +(z2-z1)^2
path_dist_sqrd = pow((betaCoords[0] - alphaCoords[0]), 2)
+ pow((betaCoords[1] - alphaCoords[1]), 2)
+ pow((betaCoords[2] - alphaCoords[2]), 2);

// Doing this instead of doing the sqrt to save doing the sqrt when not
// needed for performance
if(path_dist_sqrd <= pow(distThreshold, 2))
{


If you moved the !validAlphaIndex[j] && triActive[j] check out of checkPairValid(), to the caller (the double loop), you could avoid the overhead of the function call when those indices are not valid. Alternately, generate a list of the j indices once, and just loop over the valid indices, replacing $$\O(n^2)\$$ validity/active checks with a $$\O(n)\$$ validity/active check.

Also, you are doing calculations on path_dist_sqrd to avoid the square-root. Yet, you are still calling pow(distThreshold, 2) on the order of $$\50,000^2\$$ times, when that value is constant. pow() is not a fast function, or sqrt(x) would simply be defined as pow(x, 0.5).

Additionally, each call of checkPairValid() has to look up validAlphas[i], to get the correct index into rayTracer->m_vecTriFixedInfo[__]. This value is needed three times, and the compiler will likely cache the value. But it will not be cached across the 50,000 calls of the inner loop, where i remains constant. Instead of passing in an index to an index, look up validAlpha[i] in the outer loop, and pass that triangle index to checkPairValid(). Or lookup and cache rayTracer->m_vecTriFixedInfo[validAlphas[i]].center in the outer loop, and pass a reference to that coordinate to checkPairValid().

• Thanks for the tip on pow(), that makes good sense. I think I was mixing up the speed of doubling (equivalent to a shift) with that pow(x, 2). I've also moved the checks out of checkPairValid to avoid the function call overhead, that's a good catch. The validAlpha[i] as well was a great point, and I've fixed that now. (By the way, what's the practice on this site? Do I update my post with my fixes as I make them?) – Tyler Shellberg Jul 3 '19 at 21:44
• No, don't update your post. See What should I do when someone answers my question? – AJNeufeld Jul 3 '19 at 22:08

## Efficiency

Things that jump out.

std::vector<ABPair*> validPathsThread;


An array of pointer looks odd. Especially since ABPair is simply a pair of integers. I notice in the code you are using malloc() to allocate them which itself is a very expensive operation.

I would just remove the pointer here:

Then this code:

      ABPair* temp =
static_cast<ABPair*>(malloc(sizeof(ABPair)));

temp->alphaTriIndex = validAlphas[i];
temp->betaTriIndex = j;



simplifies to:

      validPathsThread.emplace_back(validAlphas[i], j);


## Simplify Types

You don't need to define your own what() on exceptions.

class PThreadException: public std::exception
{
virtual const char* what() const throw()
{
return "Exception occured in a pthread_attr_init or pthread_create\n";
}
};

class DynAllocationException: public std::exception
{
virtual const char* what() const throw()
{
return "Exception occured when attempting to malloc or calloc\n";
}
};


Simplify to:

// Exceptions for experimentation
struct PThreadException: public std::runtime_exception
{
};

struct DynAllocationException: public std::runtime_exception
{
DynAllocationException(): std::exception("Exception occured when attempting to malloc or calloc\n") {}
};


Simple bag types like this:

struct ABPair
{
unsigned int alphaTriIndex;
unsigned int betaTriIndex;
};


Can be simpler defined using standard types.

using ABPair = std::pair<unsigned int, unsigned int>;


## Memory Management

There is a lot of memory management code that looks like C. You should remove all calls to malloc/calloc/free and replace with new/delete then you should remove all references to new/delete with containers (using reserve() or resize() to make sure you have pre-allocated enough space).

## Throwing Exceptions

Don't do this:

  DynAllocationException e;
...
throw e;


Just do this:

            throw DynAllocationException();


## Uneeded Loop

  // Loop over all valid alphas
for(unsigned int i = 0; i < validAlphas.size(); i++)
{
{
// Do Work
}
}


This is the same as writing:

  for(unsigned int i = uThreadID; i < validAlphas.size(); i += uNumThreads) {
// Do Work
}

• Are you sure on those exceptions? They're not the same, and both fail for me. The pthread one fails with no matching function for call to 'std::exception::exception(const char [60]) and the DynAllocationException fails with error: expected unqualified-id before string constant – Tyler Shellberg Jul 3 '19 at 21:22
• I did originally have them as a pair, but converted to a struct before we realized storing some extra information about the paths would be too costly in memory. I do have one question though, could I just revise to validPaths.push_back(std::make_pair(validAlphas[i], j))? Or is there a reason to not go that route? – Tyler Shellberg Jul 3 '19 at 21:26
• Fixed the exceptions. Yes you can do: validPaths.push_back(std::make_pair(validAlphas[i], j)) But using validPaths.emplace_back(validAlphas[i], j) would be even better. – Martin York Jul 3 '19 at 22:06
• Those exceptions still don't seem to work: error: expected class-name before '{' token right after the inheritance. Also, emplace_back doesn't work either on a vector of pairs,I get 'class std::vector<std::pair<unsigned int, unsigned int> >' has no member named 'emplace_back'; did you mean 'push_back'? – Tyler Shellberg Jul 3 '19 at 22:44
• std::runtime_exception is defined in <stdexcept>. The method emplace_back() comes in C++11 – Martin York Jul 4 '19 at 5:53

Just one quick comment to add to the previous answers. You have static variables like this:

static bool* triActive = NULL;


These are assigned the pointer to a memory block, and you need an at-exit function to free these memory blocks. In C++ it is considered good practice to avoid “naked pointers” like these. If you follow that advice, you won’t need the at-exit function any more.

For example, you could instead do:

std::vector<bool> triActive;


Your initialize function then just needs to resize the array to the appropriate size:

if(triActive.empty())
triActive.resize(numTris, false);


Note that you also don’t need the static keyword here. In file scope, all variables are global, and will live for the duration of the program (since this is a DLL/SO, for the duration of it being loaded in memory). static here means to not make it visible outside the file.

Inside a function, the static keyword means to keep the variable alive in between function calls.

• I suppose it's not obvious, but that's a C-style way of showing an array. So faceToElementMap, validAlphaIndex, hasBeenInZone, etc, are all arrays, which get malloc'd for many, many values. (hundreds of thousands depending on the size of the mesh) – Tyler Shellberg Jul 4 '19 at 2:17
• @Tyler: You’re right, sorry, I didn’t read the code carefully enough. I’ve changed the answer to suggest something different with these static variables. – Cris Luengo Jul 4 '19 at 2:30