chull2d.h

crack modeling with xfem

/* 2008 (c) Dorival M. Pedroso */

#ifndef MPM_CHULL2D_H
#define MPM_CHULL2D_H

// STL
#include<algorithm> // for std::sort

// Local
#include "defs.h"
#include "array.h"

/** Lexicographically comparision between points. */
inline bool CompareVec2D (Vector3D const & A, Vector3D const & B)
{
	return (A(0)<B(0) || (A(0)==B(0) && A(1)<B(1)));
}

/** 2D cross product.
 * Return a positive value, if OAB makes a counter-clockwise turn,
 * negative for clockwise turn, and zero if the points are collinear. */
inline double Cross2D (Vector3D const & O, Vector3D const & A, Vector3D const & B)
{
	return (A(0)-O(0)) * (B(1)-O(1)) - (A(1)-O(1)) * (B(0)-O(0));
}

/** Convex hull using Andrew's monotone chain 2D algorithm.
 * Returns a list of points on the convex hull in counter-clockwise order. */
inline void CHull2D (Array<Vector3D> & P, Array<size_t> & H)
{
	// Input
	int n = P.Size();
	int k = 0;
	Array<Vector3D> h(2*n); // points on the hull
	Array<size_t>   j(2*n); // indexes of the points in the hull

	// Sort points lexicographically
	std::sort (P.GetPtr(), P.GetPtr()+P.Size(), CompareVec2D);

	// Build lower hull
	for (int i=0; i<n; i++)
	{
		while (k>=2 && Cross2D (h[k-2], h[k-1], P[i])<=0) k--;
		h[k++] = P[i];
		j[k-1] = i;
	}

	// Build upper hull
	for (int i=n-2, t=k+1; i>=0; i--)
	{
		while (k >= t && Cross2D (h[k-2], h[k-1], P[i])<=0) k--;
		h[k++] = P[i];
		j[k-1] = i;
	}

	// Results
	if (k>1)
	{
		H.Resize(k-1);
		for (int i=0; i<k-1; ++i) H[i] = j[i];
	}
}

#endif // MPM_CHULL2D_H

/* 2008 (c) Dorival M. Pedroso */