FMath - Efficient and expressive use of abstract Fields¶
General Concepts¶
The library revolves around the Field
template class, which is a convenience
wrapper around std::vector
with added operators and functions. The Eigen library
is used for vector operations.
Expression templates ensure that mathematical operations are used efficiently, avoiding temporaries.
The library provides convenience typedef
s:
FMath::scalar
= default isdouble
, can be defined differently byFMATH_SCALAR_TYPE
FMath::ScalarField
=FMath::Field<FMath::scalar>
FMath::Vector3
=Eigen::Vector3<FMath::scalar>
FMath::VectorX
=Eigen::VectorX<FMath::scalar>
FMath::VectorField
=FMath::Field<FMath::Vector3>
Using¶
Requirements:
- C++17
- OpenMP 4.5 (optional)
- (later maybe a CUDA version which supports C++17 features)
Adding this library to your project should be trivial. You can do it manually:
- copy the
FMath
folder into your directory of choice - copy the
thirdparty/Eigen
folder or provide your own - make sure the
FMath
andEigen
folders are in your include-directories - optionally define
FMATH_SCALAR_TYPE
- optionally add OpenMP compiler flags to activate parallelisation
or using CMake:
- copy the entire repository folder into your directory of choice
- TODO…
Basics¶
Reductions¶
#include <FMath/Core>
using FMath::Field;
using FMath::scalar;
using FMath::Vector3;
// Single field reduction
int N = 10;
Field<Vector3> vf(N, Vector3{0,0,1});
Vector3 mean = vf.mean();
// N-dimensional dot-product
FMath::ScalarField sf1(N), sf2(N);
sf1 *= sf2;
FMath::scalar dot = sf1.sum();
// More efficient version of the dot product:
dot = (sf1*sf2).sum();
Operators¶
#include <FMath/Core>
FMath::VectorField vf1(N), vf2(N);
FMath::ScalarField sf1(N);
// This will produce an expression object, due to auto
auto vf = vf1 + vf2*vf2;
// This will actually evaluate the expression
FMath::ScalarField sf_result = vf.dot(vf1);
Convenience math functions¶
#include <FMath/Core>
FMath::VectorField vf1(N), vf2(N);
// Element-wise dot product
FMath::ScalarField sf_dot = vf1.dot(vf2);
// Element-wise cross product
FMath::VectorField vf_cross = vf1.cross(vf2);
Other convenience functions¶
Copying or re-interpreting a Field
as an Eigen::VectorX
#include <FMath/Core>
FMath::ScalarField sf(N);
FMath::VectorField vf(N);
// Copy a scalar field to a new N-dimensional vector
FMath::VectorX vec1 = sf.asRef<VectorX>();
// Copy a vector field to a new 3N-dimensional vector
FMath::VectorX vec2 = vf.asRef<VectorX>();
// Interpret a scalar field as a N-dimensional vector without copying
Eigen::Ref<VectorX> vecRef1 = sf.asRef<VectorX>();
// Interpret a vector field as a 3N-dimensional vector without copying
Eigen::Ref<VectorX> vecRef2 = vf.asRef<VectorX>();
Subsets¶
Indexed subset¶
Extracting and operating on an indexed subset of a Field
can be performed by passing a std::vector<std::size_t>
or anything compatible, such as a Field<int>
into the
access operator.
#include <FMath/Core>
using FMath::Field;
using FMath::scalar;
// A Field of size N1 and an index list of size N2<N1
Field<scalar> sf1(N1);
Field<int> index_list1(N2);
// Set the indices of the Field entries you wish to extract...
// (this can also be used to re-order a Field)
// Extract the indexed set
Field<scalar> sf_subset1 = sf1[index_list1];
// Extract a small set via an initializer list
Field<scalar> sf_subset2 = sf1[{0,3,22}];
// Operate on subsets, combining different index lists
Field<scalar> sf2(N1);
Field<scalar> sf3(N3);
Field<int> index_list2(N2);
sf1[index_list1] = sf2[index_list1] + sf3[index_list2];
Contiguous slices¶
Field.slice()
takes:
std::size_t begin = 0
std::optional<std::size_t> end = {}
std::size_t stride = 1
#include <FMath/Core>
using FMath::Field;
using FMath::scalar;
Field<scalar> sf1{1,2,3,4,5};
// Create field from strided slice of other field
Field<scalar> sf2 = sf1.slice(0, {}, 2);
// Expected contents
// sf2.size() == 3
// sf2[0] == 1
// sf2[1] == 3
// sf2[2] == 5
// Assign slice from slice
sf2.slice(0,1) = sf1.slice(2,3);
// Expected contents
// sf2[0] == 3
// sf2[1] == 4
// Resize to size `5` and set entire field to value `2`
sf2.resize(5);
sf2 = 2;
// Assign value `3` to strided slice
sf2.slice(0, {}, 2) = 3;
// Expected contents
// sf2[0] == 3
// sf2[1] == 2
// sf2[2] == 3
// sf2[3] == 2
// sf2[4] == 3
Lambdas¶
Simple example¶
#include <FMath/Core>
using FMath::Field;
using FMath::scalar;
Field<scalar> sf(5, 0.0);
auto lambda = [](std::size_t i, scalar& val)
{
// Every value of the field is set to a calculated value
val = i+5;
};
sf.apply_lambda(lambda);
// Example entries of `sf`
// sf[0] == 5
// sf[3] == 8
Complex example¶
This example shows a reasonably complex calculation:
- each entry of
intermediate
should contain the sum of its “neighbours” inorientations
- each entry of
res
should contain the scalar product of said sum withorientations
As should be, the first operation (in form of a lambda) and the second will be performed in
parallel over the elements of the result Field
.
#include <FMath/Core>
using FMath::Field;
using FMath::scalar;
using FMath::Vector3;
int N = 5;
std::vector<int> relative_indices({-1,+1});
Field<Vector3> orientations(N, Vector3{0,0,1});
Field<Vector3> intermediate(N, Vector3{0,0,0});
Field<scalar> res(N, 0.0);
auto lambda = [&](std::size_t idx, Vector3& val)
{
// Inside entries get both directions
if( idx > 0 && idx < orientations.size()-1 )
{
for( auto& rel : relative_indices )
{
val += orientations[idx+rel];
}
}
// Left end: right-hand side
else if( idx == 0 )
val += orientations[idx+1];
// Right end: left-hand side
else if( idx == orientations.size()-1 )
val += orientations[idx-1];
};
// Result: scalar field of products of vectors in `orientations` and `intermediate`
res = orientations.dot(intermediate.applied_lambda(lambda));
// Contents of result:
// res[0] == 1
// res[1] == 2
// res[2] == 2
// res[3] == 2
// res[4] == 1