// Parasolid.cpp: implementation of the CParasolid class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------------
CParasolid::CParasolid()
//---------------------------------------------------------------------------
{// initialize

	RegisterFrustrum();

#ifdef PARASOLID
	PK_DELTA_frustrum_t frustrum = {
		FRU_delta_open_for_write,
		FRU_delta_open_for_read,
		FRU_delta_close,
		FRU_delta_write,
		FRU_delta_read,
		FRU_delta_delete,
	};

	PK_ERROR_code_t result1 = PK_DELTA_register_callbacks(frustrum);
	ASSERT(result1 == PK_ERROR_no_errors);

	// startup modeler
	
	PK_SESSION_start_o_t options;
	PK_SESSION_start_o_m(options);
	PK_ERROR_code_t result = PK_SESSION_start(&options);
	ASSERT(result == PK_ERROR_no_errors);

#else

	AfxMessageBox("This executable has been compiled without Parasolid library.\n"
		"No STL geometry will be produced.");

#endif

}

//---------------------------------------------------------------------------
CParasolid::~CParasolid()
//---------------------------------------------------------------------------
{// shutdown modeler

#ifdef PARASOLID
	PK_ERROR_code_t result = PK_SESSION_stop();
	ASSERT(result == PK_ERROR_no_errors);
#endif

}

//---------------------------------------------------------------------------
void CParasolid::Test()
//---------------------------------------------------------------------------
{

	BODY sheet = CreateSolidTriangle(CVec(0,0,1), CVec(1,1,1), CVec(1,2,3), 0.1);
//	BODY block = CreateSolidBlock(CVec(0.1,0,0), CVec(0,1,1).Normalized(), CVec(0,1,-1).Normalized(), CVec(1,2,3));
//	BODY cylinder = CreateSolidCylinder(CVec(0,0,0), CVec(0,0,1), CVec(0,1,0), 0.6, 1.1);

//	BODY cylinder2 = CreateSolidCylinder(CVec(3.2,0,0), CVec(0,1,0), 0.5, 3);
//	BODY sphere = CreateSolidSphere(CVec(0,0,0), 1.2);
//	Unite(block, cylinder);
//	BODY ass2 = Unite(cylinder2, ass);

//	BODY ass = Subtract(block, cylinder);

//	BlendAll(ass, 0.05);
	SaveModel(sheet, "ptest11"); 

}

//---------------------------------------------------------------------------
int CParasolid::report_ifail(int ifail, int expect)
//---------------------------------------------------------------------------
{
	if ( ifail != expect ) {
		TRACE("Non-zero ifail %d. Expected %d\n",ifail,expect);
		return 1;
    } else {
		return 0;
    }

}

//---------------------------------------------------------------------------
void CParasolid::wrong_result(int nwrong)
//---------------------------------------------------------------------------
{
	TRACE("*********Unexpected number %d *********\n",nwrong);
}

//---------------------------------------------------------------------------
void CParasolid::AcceptanceTest()
//---------------------------------------------------------------------------
{
	// declarations

#ifdef PARASOLID
	double vec1[3], vec2[3];
	double cyl_rad  = 2.0;
	double cyl_ht   = 12.0;
	double distance = 99.0;
	double width    =  6.0;
	double height   =  8.0;
	double depth    = 10.0;
	int ifail,nfaces,nbod,vrsion,nchars;
	int pos,nitems;
	int tworld,tass1,tget,tass2,ttran,tbox1,tbox2;
	int tcyl,tfl1,tfl2,tbyli,tbody[1];
	int        usrfld = 0;
	int        pcode  = 5;
	int  pival  = 2000000;
	double   prval  = 0.0;
	int      fault_tokens;
	int        fault_tags;
	int chcken_mxflts = 0;
	int chcken_nopts  = 0;
	int chcken_option[1];
	int             pdata;
	int            nfault;
	int    fg_key_len = 5;
	char        fg_key[6];
	int     fg_nspace = 0;
	int     fg_nints  = 0;
	int        fg_ints[1];
	int     fg_nreals = 0;
	double    fg_reals[1];
	int       fg_tag_surf;
	int     token1, token2;

	strcpy(fg_key, "dummy");

	TRACE("       Parasolid Acceptance Test\n");
	TRACE("       =========================\n");
	TRACE("\n");
	TRACE(" This acceptance test will try out some of the features which\n");
	TRACE("are available through the Kernel Interface of Parasolid.\n");
	TRACE(" Although the tests are not exhaustive, they make use\n");
	TRACE("use of some of the major parts of the system, including:\n");
	TRACE(" - starting and stopping the modeller\n");
	TRACE(" - model creation and deletion\n");
	TRACE(" - Boolean operations\n");
	TRACE(" - archived model files\n");
	TRACE(" - snapshot files of tag memory\n");
	TRACE(" - error handling\n");
	TRACE(" - rollback.\n");
	TRACE("\n");
	TRACE(" The following Parasolid routines are called:\n");
	TRACE("\n");
	TRACE("APPTRA CHCKEN CLABYS CRBXSO CRCYSO CRETRA DELENT GETMOD GETSNP\n");
	TRACE("GTTGLI IDCOEN INTBYS MERGEN ROLBLM ROLSMK SAVMOD SAVSNP SEINTP\n");
	TRACE("STAMOD STOMOD UNIBYS\n");
	TRACE("\n");
	TRACE(" The following files will be generated by the test, and can be\n");
	TRACE("deleted on completion:\n");
	TRACE(" - jfile.jnl_txt (journal file)\n");
	TRACE(" - xfile.xmt_txt (archived body)\n");
	TRACE(" - sfile.snp_bin (snapshot file)\n");
	TRACE("\n");
	TRACE(" As the test proceeds, messages are output which give the name\n");
	TRACE("of the Parasolid routine about to  be called, and a brief\n");
	TRACE("decription of what it will do. If an unexpected result is \n");
	TRACE("detected, the tests will stop.  Hence, the routine causing an\n");
	TRACE("error can be found from the messages.\n");
	TRACE("\n");
	TRACE("\n");

	/*
    From V7.0, Parasolid is `DLL-capable', ie, the Parasolid library contains
    stub functions which transfer the Downward Interface calls to registered
    callback functions, allowing Parasolid to be delivered as a self-contained
    image as well as an archive/object library. This does affect how the
    application is linked against Parasolid in order to obtain a correct
    runtime match between the application's Frustrum and the Parasolid stubs,
    ie, for the parasolid_test.c program, the behaviour of the example
    frustrum.lib and fg.lib.
	
	  We choose to always register all our example routines: when linking
	  against the archive library
	  cc parasolid_test.c frustrum.lib fg.lib parasolid.lib
	  or      cc parasolid_test.c frustrum.o   fg.o   parasolid.lib
	  then our routines are preloaded into the image and the stubs in
	  parasolid.lib are never required - the registry does occur but nothing in
	  Parasolid then uses the registration. Note that the V6 equivalent
	  cc parasolid_test.c parasolid.lib frustrum.lib fg.lib
	  will not work: only the stub routines are loaded from the Parasolid library
	  and STAMOD will fail. When linking against a standalone image
	  cc parasolid_test.c frustrum.lib fg.lib parasolid.so
	  or      cc parasolid_test.c frustrum.o   fg.o   parasolid.so
	  then both the example routines and the Parasolid stubs are loaded, the
	  stubs being hidden in the image and performing immediate transfer to the
	  active code from the example libraries.
	*/
	
	PK_SESSION_frustrum_t fru;
	TRACE("Registering Example Frustrum ...\n");
	PK_SESSION_frustrum_o_m( fru );
	fru.fstart = FSTART;
	fru.fabort = FABORT;
	fru.fstop  = FSTOP;
	fru.fmallo = FMALLO;
	fru.fmfree = FMFREE;
	fru.gosgmt = GOSGMT;
	fru.goopsg = GOOPSG;
	fru.goclsg = GOCLSG;
	fru.gopixl = GOPIXL;
	fru.gooppx = GOOPPX;
	fru.goclpx = GOCLPX;
	fru.ffoprd = FFOPRD;
	fru.ffopwr = FFOPWR;
	fru.ffclos = FFCLOS;
	fru.ffread = FFREAD;
	fru.ffwrit = FFWRIT;
	fru.ffoprb = FFOPRB;
	fru.ffseek = FFSEEK;
	fru.fftell = FFTELL;
	fru.fgcrcu = FGCRCU;
	fru.fgcrsu = FGCRSU;
	fru.fgevcu = FGEVCU;
	fru.fgevsu = FGEVSU;
	fru.fgprcu = FGPRCU;
	fru.fgprsu = FGPRSU;
	PK_SESSION_register_frustrum( &fru );

	// begin test
	
	TRACE("STAMOD: Start the modeller and check that the version is suitable.\n");
	TRACE("        A journal file jfile.jnl_txt will be written.\n");
	nchars = 5;
	token1 = KI_TRUE;
	STAMOD(&token1,&nchars,"jfile",&usrfld, &tworld, &vrsion, &ifail);
	if ( report_ifail(ifail,0) != 0 ) return;
	TRACE("STAMOD: returns version %d\n",vrsion);

	TRACE("SEINTP: Set interface parameter to allow rollback and rollforward\n");
	SEINTP( &pcode, &pival, &prval, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("DELENT: Attempt to delete invalid entity (the null tag)\n");
	TRACE("        This should result in the ifail KI_not_a_tag  22\n");
	token1 = NULTAG;
	DELENT( &token1, &ifail );
	if ( report_ifail(ifail,KI_not_a_tag) != 0 ) return;

	TRACE("CRBXSO: Create a block\n");
	vec1[0] = 0.0;
	vec1[1] = 0.0;
	vec1[2] = 0.0;
	vec2[0] = 0.0;
	vec2[1] = 0.0;
	vec2[2] = 1.0;
	CRBXSO( vec1, vec2, &width, &height, &depth, &tbox1, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("IDCOEN: Get a list of all the faces and check that there are\n");
	TRACE("        the correct number of them (6).\n");
	token1 = TYTOFA;
	IDCOEN( &tbox1, &token1, &tfl1, &nfaces, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if( nfaces != 6 ) wrong_result( nfaces);

	TRACE("SAVMOD: Save the block in a file XFILE. This will be read in later.\n");
	nchars = 5;
	SAVMOD( &tbox1, &nchars, "xfile", &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("DELENT: Delete the block.\n");
	DELENT( &tbox1, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CRBXSO: Create another solid block\n");
	vec1[0] = 0.0;
	vec1[1] = 0.0;
	vec1[2] = 0.0;
	vec2[0] = 0.0;
	vec2[1] = 0.0;
	vec2[2] = 1.0;
	CRBXSO( vec1, vec2, &width, &height, &depth, &tbox2, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CRCYSO: Create a solid cylinder\n");
	vec1[0] = 0.0;
	vec1[1] = 0.0;
	vec1[2] = -1.0;
	vec2[0] = 0.0;
	vec2[1] = 0.0;
	vec2[2] = 1.0;
	cyl_rad = 2.0;
	cyl_ht  = 12.0;
	CRCYSO( vec1, vec2, &cyl_rad, &cyl_ht, &tcyl, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CRETRA: Create a translation transformation\n");
	vec1[0] = 1.0;
	vec1[1] = 0.0;
	vec1[2] = 0.0;
	distance = 99.0;
	CRETRA( vec1, &distance, &ttran, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("ROLSMK: Set rollback mark\n");
	ROLSMK( &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("APPTRA: Apply transformation to the cylinder to move it away from\n");
	TRACE("        the block.\n");
	APPTRA( &tcyl, &ttran, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CLABYS: Test if the bodies clash, which they should not do.\n");
	token1 = NULTAG;
	token2 = KI_FALSE;
	CLABYS( &tbox2, &token1, &tcyl,&token1, &token2,
		&tfl1, &tfl2, &nfaces, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	TRACE("        Number of faces returned by CLABYS %d\n", nfaces);
	if( nfaces != 0 ) wrong_result( nfaces );

	TRACE("ROLBLM: Roll back to undo the transformation.\n");
	ROLBLM( &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CLABYS: Test that the bodies now do clash.\n");
	token1 = NULTAG;
	token2 = KI_TRUE;
	CLABYS( &tbox2, &token1, &tcyl, &token1, &token2,
		&tfl1, &tfl2, &nfaces, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	TRACE("        Number of faces now returned by CLABYS %d\n", nfaces);
	if( nfaces != 2 ) wrong_result( nfaces);

	TRACE("INTBYS: Intersect the block and the cylinder\n");
	TRACE("        A single new body should be created\n");
	INTBYS( &tcyl, &tbox2, &tass1, &nbod, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nbod != 1 ) wrong_result( nbod );

	TRACE("IDCOEN: Get a list of the bodies in the resulting assembly\n");
	token1 = TYTOBY;
	IDCOEN( &tass1, &token1, &tbyli, &nbod, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nbod != 1 ) wrong_result( nbod );

	TRACE("GTTGLI: Get the body from the list\n");
	pos = 1;
	nitems = 1;
	GTTGLI( &tbyli, &pos, &nitems, tbody, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	TRACE("CHCKEN: Check the body\n");
	CHCKEN( &tbody[0], &chcken_mxflts, &chcken_nopts, chcken_option,
		&fault_tokens, &fault_tags,
		&pdata, &nfault, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("SAVSNP: Save the state of tag memory to a file\n");
	nchars = 5;
	token1 = KI_FALSE;
	SAVSNP( &nchars, "sfile", &token1, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("DELENT: Delete the assembly\n");
	DELENT( &tass1, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("STOMOD: Stop the modeller.\n");
	STOMOD(&ifail);
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("STAMOD: Restart the modeller, this time without a journal file.\n");
	nchars = 0;
	token1 = KI_FALSE;
	STAMOD(&token1,&nchars,"jfile",&usrfld, &tworld, &vrsion, &ifail);
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("GETSNP: Restore the saved state of tag memory\n");
	nchars = 5;
	token1 = KI_FALSE;
	GETSNP( &nchars, "sfile", &token1, &token1, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("GETMOD: Try to get non-existent model. This should result in\n");
	TRACE("        ifail KI_key_not_found  58\n");
	nchars = 6;
	GETMOD( &nchars, "xyz123", &tget, &ifail );
	if ( report_ifail(ifail,KI_key_not_found) != 0 ) return;

	TRACE("GETMOD: Get the saved body.\n");
	nchars = 5;
	GETMOD( &nchars, "xfile", &tget, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("CRCYSO: Create another solid cylinder\n");
	vec1[0] = 0.0;
	vec1[1] = 0.0;
	vec1[2] = 0.0;
	vec2[0] = 0.0;
	vec2[1] = 0.0;
	vec2[2] = 1.0;
	cyl_rad = 2.0;
	cyl_ht  = 12.0;
	CRCYSO( vec1, vec2, &cyl_rad, &cyl_ht, &tcyl, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("UNIBYS: Unite the new cylinder and saved block\n");
	UNIBYS( &tget, &tcyl, &tass2, &nbod, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nbod != 1 ) wrong_result( nbod );

	TRACE("IDCOEN: Get a list of the bodies from the resulting assembly\n");
	token1 = TYTOBY;
	IDCOEN( &tass2, &token1, &tbyli, &nbod, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nbod != 1 ) wrong_result( nbod);

	TRACE("GTTGLI: Get the body from the list\n");
	pos = 1;
	nitems = 1;
	GTTGLI( &tbyli, &pos, &nitems, tbody, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	TRACE("CHCKEN: Check the body\n");
	CHCKEN( &tbody[0], &chcken_mxflts, &chcken_nopts, chcken_option,
		&fault_tokens, &fault_tags,
		&pdata, &nfault, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("IDCOEN: Get a list of all the faces and check that there are 9\n");
	token1 = TYTOFA;
	IDCOEN( &tbody[0], &token1, &tfl1, &nfaces, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nfaces != 9 ) wrong_result( nfaces );

	TRACE("MERGEN: Merge the superfluous faces away.\n");
	MERGEN( &tbody[0], &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("IDCOEN: Get a list of all the faces and check that there are 8\n");
	token1 = TYTOFA;
	IDCOEN( &tbody[0], &token1, &tfl1, &nfaces, &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;
	if(nfaces != 8 ) wrong_result( nfaces );

	TRACE("CRFGSU: Cause a dummy call to be made to a FG routine (FGCRSU)\n");
	CRFGSU( &fg_key_len, fg_key,
		&fg_nspace, &fg_nints, fg_ints, &fg_nreals, fg_reals,
		&fg_tag_surf, &ifail );

	TRACE("STOMOD: Close down the modeller\n");
	STOMOD( &ifail );
	if ( report_ifail(ifail,0) != 0 ) return;

	TRACE("*** TEST SUCCEEDED ***\n");
#endif

}


//---------------------------------------------------------------------------
void CParasolid::RegisterFrustrum()
//---------------------------------------------------------------------------
{// register frustrum into parasolid

#ifdef PARASOLID

	PK_SESSION_frustrum_t fru;

	// load structure with default values

	PK_SESSION_frustrum_o_m( fru );

	// load function addresses

	fru.fstart = FSTART;
	fru.fabort = FABORT;
	fru.fstop  = FSTOP;
	fru.fmallo = FMALLO;
	fru.fmfree = FMFREE;
	fru.gosgmt = GOSGMT;
	fru.goopsg = GOOPSG;
	fru.goclsg = GOCLSG;
	fru.gopixl = GOPIXL;
	fru.gooppx = GOOPPX;
	fru.goclpx = GOCLPX;
	fru.ffoprd = FFOPRD;
	fru.ffopwr = FFOPWR;
	fru.ffclos = FFCLOS;
	fru.ffread = FFREAD;
	fru.ffwrit = FFWRIT;
	fru.ffoprb = FFOPRB;
	fru.ffseek = FFSEEK;
	fru.fftell = FFTELL;
	fru.fgcrcu = FGCRCU;
	fru.fgcrsu = FGCRSU;
	fru.fgevcu = FGEVCU;
	fru.fgevsu = FGEVSU;
	fru.fgprcu = FGPRCU;
	fru.fgprsu = FGPRSU;

	// register

	PK_SESSION_register_frustrum( &fru );
#endif

}


//---------------------------------------------------------------------------
BODY CParasolid::CreateSolidBlock(CVec &origin, CVec &xaxis, CVec &yaxis, CVec &size)
//---------------------------------------------------------------------------
{// create solid block

	BODY box=0;

#ifdef PARASOLID

	ASSERT(fabs(xaxis.Length()-1) < 1E-5);
	ASSERT(fabs(yaxis.Length()-1) < 1E-5);
	ASSERT(fabs(xaxis.Dot(yaxis)) < 1E-5);

	// create local coord system

	PK_AXIS2_sf_t basis_set;

	basis_set.location.coord[0] = origin.x;
	basis_set.location.coord[1] = origin.y;
	basis_set.location.coord[2] = origin.z;

	basis_set.axis.coord[0] = xaxis.x;
	basis_set.axis.coord[1] = xaxis.y;
	basis_set.axis.coord[2] = xaxis.z;

	basis_set.ref_direction.coord[0] = yaxis.x;
	basis_set.ref_direction.coord[1] = yaxis.y;
	basis_set.ref_direction.coord[2] = yaxis.z;

	// create block

	PK_ERROR_code_t result = PK_BODY_create_solid_block(size.x, size.y, size.z, &basis_set, &box);
	ASSERT(result == PK_ERROR_no_errors);

#endif

	return box;

}

//---------------------------------------------------------------------------
void CParasolid::SaveModel(BODY tag, CString filename)
//---------------------------------------------------------------------------
{// save a model to file

#ifdef PARASOLID

	PK_PART_transmit_o_t options;
	PK_PART_transmit_o_m(options);
	options.transmit_format = PK_transmit_format_text_c;

	PK_ERROR_code_t result = PK_PART_transmit(1, &tag, filename.GetBuffer(255), &options);
	ASSERT(result == PK_ERROR_no_errors);

#endif

}

//---------------------------------------------------------------------------
BODY CParasolid::CreateSolidCylinder(CVec &origin, CVec &xaxis, CVec &yaxis, double radius, double height)
//---------------------------------------------------------------------------
{// create a solid cylinder

	BODY cylinder=0;

	ASSERT(fabs(xaxis.Length()-1) < 1E-5);
	ASSERT(fabs(yaxis.Length()-1) < 1E-5);
	ASSERT(fabs(xaxis.Dot(yaxis)) < 1E-5);

#ifdef PARASOLID

	// create local coord system

	PK_AXIS2_sf_t basis_set;

	basis_set.location.coord[0] = origin.x;
	basis_set.location.coord[1] = origin.y;
	basis_set.location.coord[2] = origin.z;

	basis_set.axis.coord[0] = xaxis.x;
	basis_set.axis.coord[1] = xaxis.y;
	basis_set.axis.coord[2] = xaxis.z;

	basis_set.ref_direction.coord[0] = yaxis.x;
	basis_set.ref_direction.coord[1] = yaxis.y;
	basis_set.ref_direction.coord[2] = yaxis.z;
	
	PK_ERROR_code_t result = PK_BODY_create_solid_cyl(radius,height,&basis_set,&cylinder);
	ASSERT(result == PK_ERROR_no_errors);

#endif

	return cylinder;

}

//---------------------------------------------------------------------------
BODY CParasolid::Unite(BODY& target, BODY tool)
//---------------------------------------------------------------------------
{// unite two bodies

	if (target == 0) {
		target = tool;
		return tool;
	}

	if (tool == 0)
		return target;

	BODY* assembly=0;

#ifdef PARASOLID

	int bodies;

	// set options

	PK_BODY_boolean_o_t options;
	PK_BODY_boolean_o_m(options);
	options.allow_disjoint = PK_LOGICAL_true;

	// perfotm subtraction

	PK_ERROR_code_t res = PK_BODY_boolean(target, 1, &tool, &options, &bodies, &assembly);
	ASSERT(res == PK_ERROR_no_errors); // 547=Nonmanifold
	ASSERT(bodies == 1);

#endif

	return *assembly;

}

//---------------------------------------------------------------------------
BODY CParasolid::Subtract(BODY target, BODY tool)
//---------------------------------------------------------------------------
{// subtract two bodies

	if (tool == 0)
		return target;

	BODY* assembly=0;
	ASSERT(target != 0);

#ifdef PARASOLID

	int bodies;

	// set options

	PK_BODY_boolean_o_t options;
	PK_BODY_boolean_o_m(options);
	options.allow_disjoint = PK_LOGICAL_true;
	options.function = PK_boolean_subtract_c;

	// perfotm subtraction

	PK_ERROR_code_t res = PK_BODY_boolean(target, 1, &tool, &options, &bodies, &assembly);
	ASSERT(res == PK_ERROR_no_errors); // 547=Nonmanifold
	ASSERT(bodies == 1);

#endif

	return *assembly;

}

//---------------------------------------------------------------------------
BODY CParasolid::Subtract(BODY target, int n, BODY* tool)
//---------------------------------------------------------------------------
{// subtract two bodies

	BODY* assembly=0;
	ASSERT(target != 0);

#ifdef PARASOLID

	int bodies;

	// set options

	PK_BODY_boolean_o_t options;
	PK_BODY_boolean_o_m(options);
	options.allow_disjoint = PK_LOGICAL_true;
	options.function = PK_boolean_subtract_c;

	// perfotm subtraction

	PK_ERROR_code_t res = PK_BODY_boolean(target, n, tool, &options, &bodies, &assembly);
	ASSERT(res == PK_ERROR_no_errors); // 547=Nonmanifold
	ASSERT(bodies == 1);

#endif

	return *assembly;

}

//---------------------------------------------------------------------------
BODY CParasolid::Unite(BODY target, int n, BODY* tool)
//---------------------------------------------------------------------------
{// Unite many bodies

	BODY* assembly=0;

	// set options

#ifdef PARASOLID

	int bodies;

	PK_BODY_boolean_o_t options;
	PK_BODY_boolean_o_m(options);
	options.allow_disjoint = PK_LOGICAL_true;

	// perfotm subtraction

	PK_ERROR_code_t res = PK_BODY_boolean(target, n, tool, &options, &bodies, &assembly);
	ASSERT(res == PK_ERROR_no_errors); // 547=Nonmanifold 22=PK_ERROR_not_a_tag
	ASSERT(bodies == 1);

#endif

	return *assembly;

}

//---------------------------------------------------------------------------
BODY CParasolid::CreateSolidSphere(CVec &origin, CVec &xaxis, CVec &yaxis, double radius)
//---------------------------------------------------------------------------
{// create a solid sphere

	ASSERT(fabs(xaxis.Length()-1) < 1E-5);
	ASSERT(fabs(yaxis.Length()-1) < 1E-5);
	ASSERT(fabs(xaxis.Dot(yaxis)) < 1E-5);

	BODY sphere=0;

#ifdef PARASOLID

	// create local coord system

	PK_AXIS2_sf_t basis_set;

	basis_set.location.coord[0] = origin.x;
	basis_set.location.coord[1] = origin.y;
	basis_set.location.coord[2] = origin.z;

	basis_set.axis.coord[0] = xaxis.x;
	basis_set.axis.coord[1] = xaxis.y;
	basis_set.axis.coord[2] = xaxis.z;

	basis_set.ref_direction.coord[0] = yaxis.x;
	basis_set.ref_direction.coord[1] = yaxis.y;
	basis_set.ref_direction.coord[2] = yaxis.z;

	// create sphere

	PK_ERROR_code_t result = PK_BODY_create_solid_sphere(radius, &basis_set, &sphere);
	ASSERT(result == PK_ERROR_no_errors);

#endif

	return sphere;

}
 
//---------------------------------------------------------------------------
BODY CParasolid::CreateSolidCone(CVec &origin, CVec &xaxis, CVec& yaxis, double apexradius, double semiang, double height)
//---------------------------------------------------------------------------
{// create a solid cone

	ASSERT(fabs(xaxis.Length()-1) < 1E-5);
	ASSERT(fabs(yaxis.Length()-1) < 1E-5);
	ASSERT(fabs(xaxis.Dot(yaxis)) < 1E-5);
	ASSERT(semiang > 0 && semiang < atan(1)*2);

	BODY cone=0;

#ifdef PARASOLID

	// create local coord system

	PK_AXIS2_sf_t basis_set;

	basis_set.location.coord[0] = origin.x;
	basis_set.location.coord[1] = origin.y;
	basis_set.location.coord[2] = origin.z;

	basis_set.axis.coord[0] = xaxis.x;
	basis_set.axis.coord[1] = xaxis.y;
	basis_set.axis.coord[2] = xaxis.z;

	basis_set.ref_direction.coord[0] = yaxis.x;
	basis_set.ref_direction.coord[1] = yaxis.y;
	basis_set.ref_direction.coord[2] = yaxis.z;

	// create cone

	PK_ERROR_code_t result = PK_BODY_create_solid_cone(apexradius, height, semiang, &basis_set, &cone);
	ASSERT(result == PK_ERROR_no_errors);

#endif

	return cone;

}

//---------------------------------------------------------------------------
void CParasolid::Blend(BODY tag, double radius)
//---------------------------------------------------------------------------
{ // round all edges of a body

#ifdef PARASOLID

	int edges;
	BODY* edge;
	int totfaces;

	do {
		
		totfaces = 0;
		
		// get list of edges
		
		PK_ERROR_code_t res = PK_BODY_ask_edges(tag, &edges, &edge);
		ASSERT(res == PK_ERROR_no_errors);
		
		for (int i = 0; i<edges; i++) {
			
			PK_LOGICAL_t is_entity;
			res = PK_ENTITY_is(edge[i], &is_entity);
			ASSERT(res == PK_ERROR_no_errors);
			
			if (is_entity) {
				
				// set blend options
				
				PK_EDGE_set_blend_constant_o_t options;
				PK_EDGE_set_blend_constant_o_m(options);
				
				// attach blend
				
				int edges1;
				BODY* edge1;
				//		res = PK_EDGE_set_blend_constant(edges, edge, radius, &options, &edges1, &edge1);
				res = PK_EDGE_set_blend_constant(1, &edge[i], radius, &options, &edges1, &edge1);
				ASSERT(res == PK_ERROR_no_errors);
				
				// fix blends
				
				PK_BODY_fix_blends_o_t foptions;
				PK_BODY_fix_blends_o_m(foptions);
				
				int nfaces;
				PK_FACE_t* facelist;
				PK_FACE_array_t* underfacelist;
				int* edgelist;
				PK_blend_fault_t fault;
				PK_EDGE_t faultedge;
				PK_TOPOL_t fault_topol;
				
				res = PK_BODY_fix_blends(tag, &foptions, &nfaces, &facelist, &underfacelist, 
					&edgelist, &fault, &faultedge, &fault_topol);
				ASSERT(res == PK_ERROR_no_errors);
				totfaces += nfaces;
				
			}
			
		}
		
	} while (totfaces > 0);

#endif

}

//---------------------------------------------------------------------------
void CParasolid::Blend(BODY tag, CVec pos, double radius)
//---------------------------------------------------------------------------
{ // round all edges of a body that poass through vector 'pos'

#ifdef PARASOLID

	int edges;
	BODY* edge;

	PK_VECTOR_t vector;
	vector.coord[0] = pos.x;
	vector.coord[1] = pos.y;
	vector.coord[2] = pos.z;
	
	// get list of edges
	
	PK_ERROR_code_t res = PK_BODY_ask_edges(tag, &edges, &edge);
	ASSERT(res == PK_ERROR_no_errors);
	
	for (int i = 0; i<edges; i++) {
		
		PK_LOGICAL_t is_entity;
		res = PK_ENTITY_is(edge[i], &is_entity);
		ASSERT(res == PK_ERROR_no_errors);
		
		if (is_entity) {
			
			PK_TOPOL_t topol;
			PK_ERROR_code_t res1 = PK_EDGE_contains_vector(edge[i], vector, &topol);
			ASSERT(res1 == PK_ERROR_no_errors);
			
			if (topol != PK_ENTITY_null) {
				
				// set blend options
				
				PK_EDGE_set_blend_constant_o_t options;
				PK_EDGE_set_blend_constant_o_m(options);
				
				// attach blend
				
				int edges1;
				BODY* edge1;
				//		res = PK_EDGE_set_blend_constant(edges, edge, radius, &options, &edges1, &edge1);
				res = PK_EDGE_set_blend_constant(1, &edge[i], radius, &options, &edges1, &edge1);
				ASSERT(res == PK_ERROR_no_errors);
				
				// fix blends
				
				PK_BODY_fix_blends_o_t foptions;
				PK_BODY_fix_blends_o_m(foptions);
				
				int nfaces;
				PK_FACE_t* facelist;
				PK_FACE_array_t* underfacelist;
				int* edgelist;
				PK_blend_fault_t fault;
				PK_EDGE_t faultedge;
				PK_TOPOL_t fault_topol;
				
				res = PK_BODY_fix_blends(tag, &foptions, &nfaces, &facelist, &underfacelist, 
					&edgelist, &fault, &faultedge, &fault_topol);
				ASSERT(res == PK_ERROR_no_errors);
				
			}
			
		}
		
	}

#endif

}

//---------------------------------------------------------------------------
BODY CParasolid::CreateSolidTriangle(CVec p0, CVec p1, CVec p2, double thickness)
//---------------------------------------------------------------------------
{// create a solid triangle in space

	CVec norm = ((p1-p0).Cross(p2-p0)).Normalized();
	BODY sheet=0;

#ifdef PARASOLID

	// create a sheet profile

	PK_VECTOR_t v[3];

	v[0].coord[0] = p0.x;
	v[0].coord[1] = p0.y;
	v[0].coord[2] = p0.z;
	v[1].coord[0] = p1.x;
	v[1].coord[1] = p1.y;
	v[1].coord[2] = p1.z;
	v[2].coord[0] = p2.x;
	v[2].coord[1] = p2.y;
	v[2].coord[2] = p2.z;

	PK_BODY_create_sheet_planar_o_t options;
	PK_BODY_create_sheet_planar_o_m(options);

	options.plane.location.coord[0] = p0.x;
	options.plane.location.coord[1] = p0.y;
	options.plane.location.coord[2] = p0.z;
	options.plane.axis.coord[0] = norm.x;
	options.plane.axis.coord[1] = norm.y;
	options.plane.axis.coord[2] = norm.z;

	PK_ERROR_code_t res = PK_BODY_create_sheet_planar(3, v, &options, &sheet);
	ASSERT(res == PK_ERROR_no_errors);

	// create an extrusion
	
	PK_BODY_thicken_o_t options1;
	PK_BODY_thicken_o_m(options1);

	PK_TOPOL_track_r_t tracking;
	PK_TOPOL_local_r_t results;

	res = PK_BODY_thicken_3(sheet, thickness/2, thickness/2, 0.001, &options1, &tracking, &results);
	ASSERT(res == PK_ERROR_no_errors);

#endif

	return sheet;

}

//---------------------------------------------------------------------------
void CParasolid::Rotate(BODY &tag, CVec pos, CVec ax, double ang)
//---------------------------------------------------------------------------
{ // rotate body about an axis

#ifdef PARASOLID

	PK_VECTOR_t position;
	position.coord[0] = pos.x;
	position.coord[1] = pos.y;
	position.coord[2] = pos.z;

	PK_VECTOR1_t axis;
	axis.coord[0] = ax.x;
	axis.coord[1] = ax.y;
	axis.coord[2] = ax.z;

	// create transformation

	PK_TRANSF_t transf;

	PK_ERROR_code_t res = PK_TRANSF_create_rotation(position, axis, ang, &transf);
	ASSERT(res == PK_ERROR_no_errors);

	// apply transformation

	int n_replaces;
	PK_GEOM_t* replaces;
	PK_LOGICAL_t* exact;

	res = PK_BODY_transform(tag, transf, 0.001, &n_replaces, &replaces, &exact);
	ASSERT(res == PK_ERROR_no_errors);

#endif

}

//---------------------------------------------------------------------------
bool CParasolid::SaveModelSTL(BODY tag, double tol, CString filename)
//---------------------------------------------------------------------------
{ // generate STL file. Return success

  // BINARY STL file format is accessed by byte. The format is as follows: the first 80 bytes
  // are used for description, and the next 4 bytes represents the total number of the facets
  // (Long Int), followed by the facet information (normal and 3 vertices), the normal and vertices
  // are stored in floating point format, each occupying 4 bytes. At the end of each facet information
  // section, there are two bytes spaces, then the next facet is repeated till the end of the file. 
  // When BINARY format is used to describe STL file, the data size is much smaller than ASCII format, 
  // so most STL files available now use BINARY format. 

  // <BINARY STL file format>::=<STL file entity name><facet number N><facet info> 
  // <STL file entity name >::=<80 bytes entity name, spaces are used to fill the blank> 
  // <facet number N>::=<4 bytes long integer> 
  // <facet info>::=<facet normal><facet vertices><2 bytes spaces><facet normal ><facet vertices><2 bytes spaces> ...  ... 
  // <facet normal>::=<lx,ly,lz, float, 12 bytes> 
  // <facet vertex coordinates>::=<x1,y1,z1,x2,y2,z2,x3,y3,z3, float, 36 bytes>  

  // (see http://www.vr.clemson.edu/rp/rp_stlfile.htm)

#ifdef PARASOLID

	// get facet decription

	PK_TOPOL_facet_o_t options;
	PK_TOPOL_facet_mesh_o_m(options.control);
	PK_TOPOL_facet_choice_o_m(options.choice);

	// set tolerance

	options.control.is_curve_chord_tol = PK_LOGICAL_true;
	options.control.curve_chord_tol = tol;
	options.control.is_surface_plane_tol = PK_LOGICAL_true;
	options.control.surface_plane_tol = tol;
	options.control.is_facet_plane_tol = PK_LOGICAL_true;
	options.control.facet_plane_tol = tol;
	
	// ask for facet normals as well

	options.choice.vertex_normal = PK_LOGICAL_true;
	options.choice.normal_vec = PK_LOGICAL_true;

	// compute

	PK_TOPOL_facet_r_t tables;

	PK_ERROR_code_t res = PK_TOPOL_facet(1, &tag, NULL, NULL, &options, &tables);
	ASSERT(res == PK_ERROR_no_errors);
	ASSERT(tables.number_of_facets > 3);

	// generate file

	CFile file;
	if (file.Open(filename, CFile::modeCreate | CFile::modeWrite) == 0) {
		return false;
	}

	STLHEADER header;
	STLFACETINFO facetinfo;

	header.nfacets = tables.number_of_facets;
	strcpy(header.description, "Livetruss / mm units. (C) Lipson/Pollack, Brandeis Univeristy");
	file.Write(&header, sizeof(header));

	int i;
	int facets=0;
	int lastfacet = -1;
	int vcnt = 0;
	CVec n1,n2,n3;
	int fcnt = 0;

	for (i=0; i<tables.facet_fin.length; i++) {
		int facet = tables.facet_fin.data[i].facet;
		int fin = tables.facet_fin.data[i].fin;
		if (facet != lastfacet) {
			if (lastfacet != -1) { // write facet
				ASSERT(vcnt = 3);
				CVec v1(facetinfo.x1,facetinfo.y1,facetinfo.z1);
				CVec v2(facetinfo.x2,facetinfo.y2,facetinfo.z2);
				CVec v3(facetinfo.x3,facetinfo.y3,facetinfo.z3);
//				CVec norm = -((v3-v2).Cross(v2-v1)).Normalized();
				CVec norm2 = (n1+n2+n3).Normalized();
//				if (norm.Dot(norm2) < 0) { // ensure normals have correct sign
//					norm = -norm;
//				} 
				facetinfo.nx = (float) norm2.x;
				facetinfo.ny = (float) norm2.y;
				facetinfo.nz = (float) norm2.z;
				fcnt++;
				file.Write(&facetinfo, 50);
			}
			lastfacet = facet;
			vcnt = 0;
		}
		ASSERT(fin >= 0 && fin < tables.fin_vertex.length);
		int vertex = tables.fin_vertex.vertex[fin];
		ASSERT(vertex >= 0 && vertex < tables.vertex_point.length);
		int point = tables.vertex_point.point[vertex];
		ASSERT(point >= 0 && point < tables.point_vec.length);
		int normal = tables.vertex_normal.normal[vertex];
		ASSERT(normal >= 0 && normal < tables.normal_vec.length);
		if (vcnt == 0) {
			facetinfo.x1 = (float) tables.point_vec.vec[point].coord[0]*1000;
			facetinfo.y1 = (float) tables.point_vec.vec[point].coord[1]*1000;
			facetinfo.z1 = (float) tables.point_vec.vec[point].coord[2]*1000;
			n1 = CVec(tables.normal_vec.vec[normal].coord[0],tables.normal_vec.vec[normal].coord[1],tables.normal_vec.vec[normal].coord[2]);
			vcnt++;
		} else if (vcnt == 1) {
			facetinfo.x2 = (float) tables.point_vec.vec[point].coord[0]*1000;
			facetinfo.y2 = (float) tables.point_vec.vec[point].coord[1]*1000;
			facetinfo.z2 = (float) tables.point_vec.vec[point].coord[2]*1000;
			n2 = CVec(tables.normal_vec.vec[normal].coord[0],tables.normal_vec.vec[normal].coord[1],tables.normal_vec.vec[normal].coord[2]);
			vcnt++;
		} else if (vcnt == 2) {
			facetinfo.x3 = (float) tables.point_vec.vec[point].coord[0]*1000;
			facetinfo.y3 = (float) tables.point_vec.vec[point].coord[1]*1000;
			facetinfo.z3 = (float) tables.point_vec.vec[point].coord[2]*1000;
			n3 = CVec(tables.normal_vec.vec[normal].coord[0],tables.normal_vec.vec[normal].coord[1],tables.normal_vec.vec[normal].coord[2]);
			vcnt++;
		} else {
			ASSERT(false);
		}
	}

				ASSERT(vcnt = 3);
				CVec v1(facetinfo.x1,facetinfo.y1,facetinfo.z1);
				CVec v2(facetinfo.x2,facetinfo.y2,facetinfo.z2);
				CVec v3(facetinfo.x3,facetinfo.y3,facetinfo.z3);
//				CVec norm = -((v3-v2).Cross(v2-v1)).Normalized();
				CVec norm2 = (n1+n2+n3).Normalized();
//				if (norm.Dot(norm2) < 0) { // ensure normals have correct sign
//					norm = -norm;
//				} 
				facetinfo.nx = (float) norm2.x;
				facetinfo.ny = (float) norm2.y;
				facetinfo.nz = (float) norm2.z;
				fcnt++;
				file.Write(&facetinfo, 50);

	ASSERT(fcnt == header.nfacets);

	file.Close();

	TRACE("Generated %d facets\n", header.nfacets);
	// free alocated tables

#if 0
	if (tables.facet_fin.data != 0)
		free(tables.facet_fin.data);
	if (tables.fin_vertex.vertex != 0)
		free(tables.fin_vertex.vertex);
	if (tables.vertex_point.point != 0)
		free(tables.vertex_point.point);
	if (tables.point_vec.vec != 0)
		free(tables.point_vec.vec);
	if (tables.fin_fin.fin != 0)
		free(tables.fin_fin.fin);
#endif

#endif

	return true;

}

//---------------------------------------------------------------------------
void CParasolid::Disjoin(BODY tag, int& bodies, BODY **body)
//---------------------------------------------------------------------------
{

#ifdef PARASOLID

	PK_ERROR_code_t res = PK_BODY_disjoin(tag, &bodies, body);
	ASSERT(res == PK_ERROR_no_errors);

#endif

}

//---------------------------------------------------------------------------
void CParasolid::ComputeUnalignedBox(BODY tag, CVec &origin, CVec &xaxis, CVec &yaxis, CVec &zaxis)
//---------------------------------------------------------------------------
{// compute unaligned bonding box

#ifdef PARASOLID

	int i,j,k,m;

	// get facet decription

	PK_TOPOL_facet_o_t options;
	PK_TOPOL_facet_mesh_o_m(options.control);
	PK_TOPOL_facet_choice_o_m(options.choice);
	PK_TOPOL_facet_r_t tables;

	// set low tolerance

	double tol = 0.01;
	options.control.is_curve_chord_tol = PK_LOGICAL_true;
	options.control.curve_chord_tol = tol;
	options.control.is_surface_plane_tol = PK_LOGICAL_true;
	options.control.surface_plane_tol = tol;
	options.control.is_facet_plane_tol = PK_LOGICAL_true;
	options.control.facet_plane_tol = tol;

	PK_ERROR_code_t res = PK_TOPOL_facet(1, &tag, NULL, NULL, &options, &tables);
	ASSERT(res == PK_ERROR_no_errors);

	CVec* p = new CVec[tables.point_vec.length+1];
	for (i=0; i<tables.point_vec.length; i++) {
		p[i] = CVec(tables.point_vec.vec[i].coord[0], tables.point_vec.vec[i].coord[1], tables.point_vec.vec[i].coord[2]);
	}

	// try permutations

	double volmin = 1E30;
	int skip = tables.point_vec.length/75;
	if (skip < 1) 
		skip = 1;

	TRACE("Points to consider = %d (skip=%d)\n", tables.point_vec.length/skip, skip);

	for (i=0; i<tables.point_vec.length; i += skip) {
		for (j=0; j<i; j += skip) {
			for (k=0; k<j; k += skip) {

				// compute candidate coord system 

				CVec cxaxis = (p[j]-p[i]).Normalized();
				CVec cyaxis = (p[k]-p[i]);
				CVec czaxis = (cxaxis.Cross(cyaxis)).Normalized();
				cyaxis = (czaxis.Cross(cxaxis)).Normalized();
				if (cxaxis.Length() == 0 || cyaxis.Length() == 0 || czaxis.Length() == 0)
					continue;

				// compute bound box

				double x,y,z;
				double zmax = -1E30, xmax = -1E30, ymax = -1E30;
				double zmin = 1E30, xmin = 1E30, ymin = 1E30;

				for (m=0; m<tables.point_vec.length; m++) {
					x = p[m].Dot(cxaxis);
					y = p[m].Dot(cyaxis);
					z = p[m].Dot(czaxis);
					if (x > xmax) 
						xmax = x;
					if (x < xmin) 
						xmin = x;
					if (y > ymax) 
						ymax = y;
					if (y < ymin) 
						ymin = y;
					if (z > zmax) 
						zmax = z;
					if (z < zmin) 
						zmin = z;
					if (m == 100 && (xmax-xmin)*(ymax-ymin)*(zmax-zmin) > volmin)
						break;
				}

				double vol = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
				if (vol < volmin && vol > 0) {
					volmin = vol;
					origin = xmin*cxaxis + ymin*cyaxis + zmin*czaxis;
					xaxis = cxaxis*(xmax-xmin);
					yaxis = cyaxis*(ymax-ymin);
					zaxis = czaxis*(zmax-zmin);
				}
					
			}
		}
	}

	delete [] p;

	// rearrange axes so that xaxis is longest and zaxis is shortest

	for (i=0; i<3; i++) {
		if (xaxis.Length() < zaxis.Length()) {
			CVec tmp = xaxis;
			xaxis = zaxis;
			zaxis = tmp;
		}
		if (xaxis.Length() < yaxis.Length()) {
			CVec tmp = xaxis;
			xaxis = yaxis;
			yaxis = tmp;
		}
		if (yaxis.Length() < zaxis.Length()) {
			CVec tmp = yaxis;
			yaxis = zaxis;
			zaxis = tmp;
		}
	}

	// rearrange system so that it is right-handed

	if ((xaxis.Cross(yaxis)).Dot(zaxis)< 0) {
		// switch z axis
		origin += zaxis;
		zaxis = -zaxis;
	}

#endif

}

//---------------------------------------------------------------------------
void CParasolid::Translate(BODY tag, CVec delta)
//---------------------------------------------------------------------------
{ // translate a body by delta

#ifdef PARASOLID

	PK_VECTOR_t translation;
	translation.coord[0] = delta.x;
	translation.coord[1] = delta.y;
	translation.coord[2] = delta.z;

	// create transformation

	PK_TRANSF_t transf;

	PK_ERROR_code_t res = PK_TRANSF_create_translation(translation, &transf);
	ASSERT(res == PK_ERROR_no_errors);

	// apply transformation

	int n_replaces;
	PK_GEOM_t* replaces;
	PK_LOGICAL_t* exact;

	res = PK_BODY_transform(tag, transf, 0.001, &n_replaces, &replaces, &exact);
	ASSERT(res == PK_ERROR_no_errors);

#endif

}

//---------------------------------------------------------------------------
CVec CParasolid::ComputeBoundaryAverage(BODY tag)
//---------------------------------------------------------------------------
{ // compute pseudo centroid by averaging boundary points

	CVec average(0,0,0);

#ifdef PARASOLID

	int i;

	// get facet decription

	PK_TOPOL_facet_o_t options;
	PK_TOPOL_facet_mesh_o_m(options.control);
	PK_TOPOL_facet_choice_o_m(options.choice);
	PK_TOPOL_facet_r_t tables;

	PK_ERROR_code_t res = PK_TOPOL_facet(1, &tag, NULL, NULL, &options, &tables);
	ASSERT(res == PK_ERROR_no_errors);

	for (i=0; i<tables.point_vec.length; i++) {
		average += CVec(tables.point_vec.vec[i].coord[0], tables.point_vec.vec[i].coord[1], tables.point_vec.vec[i].coord[2]);
	}
	average /= tables.point_vec.length;

#endif

	return average;

}

//---------------------------------------------------------------------------
BODY CParasolid::CreateEmptyAssembly()
//---------------------------------------------------------------------------
{

#ifdef PARASOLID

	PK_ASSEMBLY_t assembly;
	PK_ERROR_code_t res = PK_ASSEMBLY_create_empty(&assembly);
	ASSERT(res == PK_ERROR_no_errors);
	return assembly;

#else

	return 0;

#endif

}

//---------------------------------------------------------------------------
void CParasolid::AddToAssembly(BODY assembly, BODY tag)
//---------------------------------------------------------------------------
{// Add instance of body to an assembly "as is" (without transformation);

#ifdef PARASOLID

	PK_INSTANCE_sf_t instance_sf;
	instance_sf.assembly = assembly;
	instance_sf.part = tag;
	instance_sf.transf = PK_ENTITY_null;

	PK_INSTANCE_t instance;
	PK_ERROR_code_t res = PK_INSTANCE_create(&instance_sf, &instance);
	ASSERT(res == PK_ERROR_no_errors);

#endif

}

//---------------------------------------------------------------------------
BODY CParasolid::Copy(BODY tag)
//---------------------------------------------------------------------------
{

	BODY newtag = 0;

#ifdef PARASOLID
	PK_ERROR_code_t res = PK_ENTITY_copy(tag, &newtag);
	ASSERT(res == PK_ERROR_no_errors);
#endif

	return newtag;
}