-
Notifications
You must be signed in to change notification settings - Fork 2
5: SBGAT User's Manual
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
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
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
N
/ \
S B
\
L
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.
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.
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
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
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();
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.
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");
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
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
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).
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;
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";