5: SBGAT User's Manual

SbgatCore

This page provides a breakdown of SbgatCore's capabilities. It is not nearly as exhaustive as the technical documentation, but more geared towards helping users to embed SBGAT in their own software.

SBGATRefFrame

SBGATRefFrame is used to represent a reference frame : the conjunction of an origin and a set of three orthogonal vectors defining the orientation of this frame with respect to another. This class is not used on its own, but is wrapped within an instance of SBGATFrameGraph.

In SBGAT, frames are simply identified by their names (ex: "N" , "B", "my_frame").

SBGATFrameGraph

SBGATFrameGraph is used to construct a frame conversion graph. Let us imagine that a given problem features the following frames:

N : the inertial frame of reference
B: a body-fixed frame attached to some small-body center-of-mass
L: a frame associated to some local feature on the small body
S : a spacecraft frame

A typical frame graph would then look something like

We can create this graph through

        // Creating the frame graph
        SBGATFrameGraph frame_graph;
        

        // Adding frames to the graph
	frame_graph.add_frame("N");
	frame_graph.add_frame("B");
        frame_graph.add_frame("L");
        frame_graph.add_frame("S");

        // Adding edges to the graph
	frame_graph.add_transform("N","B");
        frame_graph.add_transform("N","S");
        frame_graph.add_transform("B","L");

Notice the order of the arguments in the add_transform method. The first argument refers to the parent frame, while the second one denotes the child frame. By default, all frame transforms have their origin and their attitude MRP set to (0,0,0).

Specifying non-trivial frame transforms can be done like so (for the specific case of the B,N transform)

arma::vec::fixed<3> mrp = {0.1,0.2,0.3};
arma::vec::fixed<3> origin = {-0.1,0.2,-0.3};

frame_graph.set_transform_mrp("N","B",mrp);
frame_graph.set_transform_origin("N","B",origin);

At this stage, the origin of the B frame expressed in the N frame is thus origin, and the attitude of the B frame with respect to the N frame is represented by mrp. Specifically, mrp instantiates the direction cosine matrix [BN] as defined by Schaub and Junkins. The frame graph can then be used to convert vectors from one frame to the other.

Here is the conversion of a vector expressed in the S frame, r_in_S, into the L frame:

arma::vec::fixed<3> r_in_L = frame_graph.convert(r_in_S,"S","L");

This conversion will not conserve the vector's magnitude. However, if the converted vector is a unit vector (aka a direction), then the norm will be conserved only under the rotational component of the conversion. This can be achieved by setting the last argument of convert, conserve_norm, to false. In so many words,

arma::vec::fixed<3> u_in_L = frame_graph.convert(u_in_S,"S","L",true);

The last argument of convert is set to false by default.

SBGATFilter

A generic class serving a a shape container and providing a few useful methods to access the shape's underlying facets, vertices and edges, inheriting from vtk's vtkPolyDataAlgorithm. After being instantiated,SBGATFilter objects store the vertices coordinates and the connectivity table of the input vtkPolyData. This data duplication allows for much quicker access to the underlying coordinates than using the vtkPolyData access methods. Besides providing accessors to the shape, this class provides unit consistency method (SetScaleMeters and SetScaleKiloMeters) that must be called prior to calling Update so as to ensure consistency with the provided shape. By default, it is expected that the provided shape will have its coordinates expressed in meters.

SBGATMassProperties

This class provides methods allowing computation of volume, area, shape index, center of mass, inertia tensor and principal axes of a polyhedral mesh of constant density. It inherits from SBGATFilter. A typical use case would look like:

vtkSmartPointer<SBGATMassProperties> mass_filter = vtkSmartPointer<SBGATMassProperties>::New();
	mass_filter -> SetInputData(my_polydata);
	mass_filter -> Update();
        
        // Evaluates the oriented sum of surface normals, which should evaluate to zero 
        // if the provided data is closed.
	assert(mass_filter -> CheckClosed());

	// Getting the center of mass, computed under the assumption of constant density
	auto com_sbgat = mass_filter -> GetCenterOfMass();
        ...

SBGATPolyhedronGravityModel

SBGATPolyhedronGravityModel is a module of SBGAT allowing the construction and evaluation of a constant-density Polyhedron Gravity Model. It inherits from `SBGATMassProperties'. Analytical insight into the foundations of this gravity model can be found in Werner and Scheeres' seminal paper

Werner, R. A., & Scheeres, D. J. (1997). Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia. Celestial Mechanics and Dynamical Astronomy, 65(3), 313–344. https://doi.org/10.1007/BF00053511

The surface polyhedron gravity model implemented in SBGAT provides static methods to compute the following quantities over a subset of a queried shape model:

gravitational slopes
inertial gravity potentials
body-fixed gravity potentials
inertial gravity acceleration magnitudes
body-fixed gravity acceleration magnitudes

Construction

The construction of the Polyhedron Gravity Model follows the same steps as in SBGATMassProperties:

        // Reading the shape model from file
	vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
	reader -> SetFileName("../input/KW4Alpha.obj");
	reader -> Update(); 

	// Setting the density
	double density = 2000.0;//Density of KW4 (kg/m^3). 

	// An instance of SBGATPolyhedronGravityModel is created to evaluate the PGM of 
	// the considered polytdata
	vtkSmartPointer<SBGATPolyhedronGravityModel> pgm_filter = vtkSmartPointer<SBGATPolyhedronGravityModel>::New();
	pgm_filter -> SetInputConnection(reader -> GetOutputPort());
	pgm_filter -> SetDensity(density);
	pgm_filter -> SetScaleKiloMeters();// The input shape model has its coordinates expressed in km 
	pgm_filter -> Update();

We can now break down the above code step by step:

Reading the shape model from file

vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
reader -> SetFileName("../input/KW4Alpha.obj");
reader -> Update();

This three lines simply create an instance of VTK's vtkOBJReader allowing us to parse and load .obj files. The Update() call is critical and should not be left-out. The output of this filter is effectively an instance of vtkSmartPointer<vtkPolyData>, representing the polyhedron. The creation of this object can be delayed through the use of GetOutputPort(). Otherwise, the object could have been directly retrieved by calling reader -> GetOutputData() instead. Refer to the VTK documentation for more insight into VTK's filters.

Setting the density

double density = 2000.0; // Density of KW4 (kg/m^3).

Note that SI units are always used when asked for density or other intrinsic parameters. But one could ask "What if my shape model coordinates are expressed in kilometers?". This is addressed in the next block of code

An instance of SBGATPolyhedronGravityModel is created to evaluate the PGM of the considered polytdata

An instance of SBGATPolyhedronGravityModel is created

vtkSmartPointer<SBGATPolyhedronGravityModel> pgm_filter = vtkSmartPointer<SBGATPolyhedronGravityModel>::New();

The input port of this filter is connected to the output port of the reader. Again, this could have been replaced with SetInputData if a smart pointer to a vtkPolydata was readily available

pgm_filter -> SetInputConnection(reader -> GetOutputPort());

The density is set

pgm_filter -> SetDensity(density);

This is were the units ambiguity is resolved. Since the provided shape model is expressed in kilometers, SetScaleKiloMeters will automatically rescale it so as to have its coordinates expressed in meters in the vtkPolyData.

pgm_filter -> SetScaleKiloMeters();// The input shape model has its coordinates expressed in km

The dyads of the polyhedron gravity model are finally constructed and stored.

pgm_filter -> Update();

Evaluation

Once the Polyhedron Gravity Model has been constructed, it can be evaluated at an arbitrary position to obtain the resulting inertial gravity acceleration and/or potential. Notice the use of GetPotentialAcceleration to compute both the acceleration and potential together. This method is much faster than successive calls to GetPotential and GetAcceleration. All PGM methods take position arguments expressed in meters.

arma::vec::fixed<3>  pos = {3e3,5e3,-2e3};
arma::vec::fixed<3> pgm_acc = pgm_filter -> GetAcceleration(pos);// acceleration expressed in m/s^2
double pgm_pot = pgm_filter -> GetPotential(pos);//potential expressed in m^2/s^2
pgm_filter -> GetPotentialAcceleration(pos,pgm_pot, pgm_acc);// compute both

The coordinates of the query point and the gravity acceleration are expressed in whichever frame the shape model came in. Ensuring that the shape model is centered at its center-of-mass and aligned with its principal axes is recommended. SbgatGui provides an alignment function that ensures this is the case.

Surface gravity

Although the surface gravity acceleration and corresponding gravity potentials & slopes can perfectly be computed manually, using the previously presented methods, a self-contained static method SBGATPolyhedronGravityModel::ComputeSurfacePGM is available to greatly facilitate this process. The method's synopsis reads

static void ComputeSurfacePGM(
  vtkSmartPointer<vtkPolyData> selected_shape,
  const std::vector<unsigned int> & queried_elements,
  bool is_in_meters,
  double density,
  const arma::vec::fixed<3> & omega,
  std::vector<double> & slopes,
  std::vector<double> & potentials,
  std::vector<double> & acc_magnitudes,
  std::vector<double> & acc_body_fixed_magnitudes);

This method will populate the output vectors slopes, potentials, acc_magnitudes and acc_body_fixed_magnitudes at the queried surface elements.

Below is a simple use case demonstrating this method in conjunction with SBGATPolyhedronGravityModel::SaveSurfacePGM.

        // polydata == KW4 BARYCENTERED
 
	std::vector<unsigned int> queried_elements;
	for (int i = 0; i < polydata -> GetNumberOfCells(); ++i){
		queried_elements.push_back(i);
	}

	// Rotation period of KW4
	arma::vec omega = {0,0,0.0006312};

	std::vector<double> slopes,inertial_potentials,
	body_fixed_potentials,
	inertial_acc_magnitudes,
	body_fixed_acc_magnitudes;

	SBGATPolyhedronGravityModel::ComputeSurfacePGM(polydata,
		queried_elements,
		false,
		KW4_density_in_kg_m3,
		omega,
		slopes,
		inertial_acc_magnitudes,
		body_fixed_acc_magnitudes,
		inertial_acc_magnitudes,
		body_fixed_acc_magnitudes);


	SBGATPolyhedronGravityModel::SaveSurfacePGM(polydata,
		queried_elements,
		false,
		KW4_density_in_kg_m3 * KW4_volume_in_m3 ,
		omega,
		slopes,
		inertial_acc_magnitudes,
		body_fixed_acc_magnitudes,
		inertial_acc_magnitudes,
		body_fixed_acc_magnitudes,
		"surface_pgm.json");

SBGATSrpYorp

SBGATSphericalHarmo

SBGATSphericalHarmo provides an interface to the algorithms used in the computation and evaluation of the spherical harmonics expension of an exterior gravity field about a constant-density small body. A typical use case is shown below for a polydata representative of a small body of constant density equal to 2000 kg/m^3, whose coordinates are expressed in kilometers. The small body's enclosing radius is equal to 1500 m. The coefficients up to degree and order 10 are computed.

vtkSmartPointer<SBGATSphericalHarmo> spherical_harmonics = vtkSmartPointer<SBGATSphericalHarmo>::New();
spherical_harmonics -> SetInputData(my_polydata);
spherical_harmonics -> SetDensity(2000);
spherical_harmonics -> SetScaleKiloMeters();
spherical_harmonics -> SetReferenceRadius(1500);
spherical_harmonics -> IsNormalized(); // can be skipped as normalized coefficients is the default parameter
spherical_harmonics -> SetDegree(10);
spherical_harmonics -> Update();
...
arma::vec::fixed<3>  pos = {3e3,5e3,-2e3};// coordinates of the query point in meters
arma::vec::fixed<3> sharm_acc = spherical_harmonics -> GetAcceleration(pos);
...

Once properly defined, the instance of SBGATSphericalHarmo can be saved to a .json file storing the computed coefficients along with the rest of the specified inputs.

// The spherical harmonics are saved to a file
	spherical_harmonics -> SaveToJson("harmo.json");

An instance of SBGATSphericalHarmo can also be initialized directly from a .json file conforming with the format defined in SBGATSphericalHarmo.hpp.

vtkSmartPointer<SBGATSphericalHarmo> spherical_harmonics_from_file = vtkSmartPointer<SBGATSphericalHarmo>::New();
spherical_harmonics_from_file -> LoadFromJson("harmo.json");
// spherical_harmonics_from_file is now ready for use

SBGATObjWriter

SBGATObjWriter is no less than the exact reproduction of David Doria's vtkPolydataWriter. Its usage is fairly trivial as show below:

vtkSmartPointer<SBGATObjWriter> writer = vtkSmartPointer<SBGATObjWriter>::New();
writer -> SetInputData(polydata_to_save);
writer -> SetFileName("saved_polydata.obj");
writer -> Update();//effectively writes the polydata to the obj file

SBGATFilterUQ

Virtual class serving as the base of dedicated uncertainty quantification classes in SBGAT. This class provides setters/getters to the considered shape's covariance matrix, enabling global or local uncertainty regions to be created. It also provides a set of IO methods for saving and/or loading the covariance matrix to a text file/JSON file (non-zero partitions of the covariance only).

SBGATMassPropertiesUQ

Evaluate the formal uncertainty in the volume, center of mass, inertia tensor parametrization from a topologically-closed, constant density polyhedron. Inherits from SBGATFilterUQ. A typical use case is provided below:

// Reading shape model
	vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
	reader -> SetFileName(PATH_SHAPE.c_str());
	reader -> Update(); 

	// An instance of SBGATMassProperties is created to evaluate the mass properties
	// of the considered shape
	vtkSmartPointer<SBGATMassProperties> mass_filter = vtkSmartPointer<SBGATMassProperties>::New();
	mass_filter -> SetInputConnection(reader -> GetOutputPort());
	mass_filter -> SetDensity(DENSITY);

	// Sets the scaling factor of the input shape
	if(UNIT_IN_METERS){
		mass_filter -> SetScaleMeters();
	} else{
		mass_filter -> SetScaleKiloMeters();
	}

	// Computes the volume, center-of-mass,...
	mass_filter -> Update();

	// An instance of SBGATMassPropertiesUQ is created to perform
	// uncertainty quantification in the inertias of the considered shape
	SBGATMassPropertiesUQ mass_uq;

	// The SBGATMassPropertiesUQ is provided with the reference inertias
	mass_uq.SetModel(mass_filter);

	// Populate the shape vertices covariance
	mass_uq.ComputeVerticesCovarianceGlobal(ERROR_STANDARD_DEV,CORRELATION_DISTANCE);
	
	// Regularizing the covariance
	int regularized_eigen_values = mass_uq.RegularizeCovariance();

	std::cout << regularized_eigen_values << " eigenvalues were regularized\n";
	
	arma::mat C_CC = mass_uq.GetCovarianceSquareRoot();
	arma::mat P_CC = mass_uq.GetVerticesCovariance();
	std::cout << "Maximum absolute error in covariance square root: " << arma::abs(P_CC - C_CC * C_CC.t()).max() << std::endl;
	
	std::cout << "Saving shape covariance ...\n";
	P_CC.save(OUTPUT_DIR + "full_covariance.txt",arma::raw_ascii);

	// Saving baseline slices
	mass_uq.TakeAndSaveSlice(0,OUTPUT_DIR + "baseline_slice_x.txt",1e-6);
	mass_uq.TakeAndSaveSlice(1,OUTPUT_DIR + "baseline_slice_y.txt",1e-6);
	mass_uq.TakeAndSaveSlice(2,OUTPUT_DIR + "baseline_slice_z.txt",1e-6);



	// The analytical uncertainties are computed
	auto start = std::chrono::system_clock::now();
	mass_uq.PrecomputeMassPropertiesPartials();
	const arma::mat & partialIPartialC = mass_uq.GetPartialIPartialC();
	const arma::mat & partialComPartialC = mass_uq.GetPartialComPartialC();
	const arma::rowvec & partialVolumePartialC = mass_uq.GetPartialVolumePartialC();
	double sigma_vol = std::sqrt(arma::dot(partialVolumePartialC.t(),P_CC * partialVolumePartialC.t()));
	arma::mat cov_com = partialComPartialC * P_CC * partialComPartialC.t();
	arma::mat cov_I = partialIPartialC * P_CC * partialIPartialC.t();

	std::cout << "Standard deviation on volume (m^3): " << sigma_vol << " \n";
	cov_com.print("Covariance in center-of-mass (m^4): ");
	cov_I.print("Covariance in unit-density inertia tensor (m^10): ");

	auto end = std::chrono::system_clock::now();
	std::chrono::duration<double> elapsed_seconds = end - start;
	std::cout << "Time elapsed computing analytical uncertainties: " << elapsed_seconds.count() << " s" << std::endl;

SBGATPolyhedronGravityModelUQ

This class can be used to compute the gravity potential variance, acceleration covariance deriving from the constant density polyhedron gravity model assuming that the underlying shape vertices are outcomes of a Gaussian distribution of known mean and covariance. A typical example is shown below:

// Reading
	vtkSmartPointer<vtkOBJReader> reader = vtkSmartPointer<vtkOBJReader>::New();
	reader -> SetFileName(PATH_SHAPE.c_str());
	reader -> Update(); 

	// An instance of SBGATPolyhedronGravityModel is created to evaluate the PGM of 
	// the considered polytdata
	vtkSmartPointer<SBGATPolyhedronGravityModel> pgm_filter = vtkSmartPointer<SBGATPolyhedronGravityModel>::New();
	pgm_filter -> SetInputConnection(reader -> GetOutputPort());
	pgm_filter -> SetDensity(DENSITY);

	
	if(UNIT_IN_METERS){
		pgm_filter -> SetScaleMeters();
	} else{
		pgm_filter -> SetScaleKiloMeters();
	}

	std::cout << "Building pgm ...\n";
	pgm_filter -> Update();

	// An instance of SBGATPolyhedronGravityModelUQ is created to perform
	// uncertainty quantification from the PGM associated to the shape
	SBGATPolyhedronGravityModelUQ pgm_uq;
	pgm_uq.SetModel(pgm_filter);
	pgm_uq.PrecomputeMassPropertiesPartials();
	std::cout << "Populating shape covariance ...\n";

	// Populate the shape vertices covariance
	pgm_uq.ComputeVerticesCovarianceGlobal(ERROR_STANDARD_DEV,CORRELATION_DISTANCE);
	// Regularizing the covariance
	int regularized_eigen_values = pgm_uq.RegularizeCovariance();

	std::cout << regularized_eigen_values << " eigenvalues were regularized\n";

	std::cout << "Saving non-zero partition of shape covariance ...\n";

	// Save the covariance
	pgm_uq.SaveNonZeroVerticesCovariance(OUTPUT_DIR + "shape_covariance.json");

	// Saving baseline slices
	pgm_uq.TakeAndSaveSlice(0,OUTPUT_DIR + "baseline_slice_x.txt",0);
	pgm_uq.TakeAndSaveSlice(1,OUTPUT_DIR + "baseline_slice_y.txt",0);
	pgm_uq.TakeAndSaveSlice(2,OUTPUT_DIR + "baseline_slice_z.txt",0);
        
        // Positions where to evaluate the uncertainty model
	std::vector<arma::vec::fixed<3> > all_positions = {
		arma::vec::fixed<3>({300,0,0}),
		arma::vec::fixed<3>({400,0,0}),
		arma::vec::fixed<3>({500,0,0}),
		arma::vec::fixed<3>({-300,0,0}),
		arma::vec::fixed<3>({-400,0,0}),
		arma::vec::fixed<3>({-500,0,0}),
		arma::vec::fixed<3>({0,300,0}),
		arma::vec::fixed<3>({0,400,0}),
		arma::vec::fixed<3>({0,500,0}),
		arma::vec::fixed<3>({0,-300,0}),
		arma::vec::fixed<3>({0,-400,0}),
		arma::vec::fixed<3>({0,-500,0}),
		arma::vec::fixed<3>({0,0,300}),
		arma::vec::fixed<3>({0,0,400}),
		arma::vec::fixed<3>({0,0,500}),
		arma::vec::fixed<3>({0,0,-300}),
		arma::vec::fixed<3>({0,0,-400}),
		arma::vec::fixed<3>({0,0,-500}),
	};


	// Analytical UQ

	std::vector<arma::mat::fixed<3,3> > analytical_covariances_acc(all_positions.size());
	std::vector<double> analytical_variances_pot(all_positions.size());

	std::cout << "Computing analytical uncertainties ... ";
	auto start = std::chrono::system_clock::now();
	#pragma omp parallel for
	for (int e = 0; e < all_positions.size(); ++e){
		analytical_variances_pot[e] = pgm_uq.GetVariancePotential(all_positions[e],HOLD_MASS_CONSTANT);
		analytical_covariances_acc[e] = pgm_uq.GetCovarianceAcceleration(all_positions[e],HOLD_MASS_CONSTANT);
	}
	auto end = std::chrono::system_clock::now();

	std::chrono::duration<double> elapsed_seconds = end-start;
	std::cout << "Done computing analytical uncertainties in " << elapsed_seconds.count() << " s\n";

Provide feedback

Saved searches

Use saved searches to filter your results more quickly