helpers.h

Go to the documentation of this file.

 
#pragma once
 
#include "utils/testing.h"
 
#include <jem/base/Array.h>
#include <jem/base/CString.h>
#include <jem/base/tuple/Tuple.h>
#include <jem/numeric/Quaternion.h>
#include <jem/numeric/algebra.h>
#include <jem/numeric/func/Function.h>
 
#include <jive/Array.h>
#include <jive/util/Constraints.h>
#include <jive/util/FuncUtils.h>
 
#define WARN(key) jem::System::warn() << key << "\n";
#define ERR(key) jem::System::err() << key << "\n";
 
using jem::ALL;
using jem::Array;
using jem::idx_t;
using jem::SliceFromTo;
using jem::String;
using jem::numeric::dotProduct;
using jem::numeric::Function;
using jem::numeric::matmul;
using jem::numeric::MatmulChain;
using jem::numeric::norm2;
using jem::numeric::Quaternion;
 
using jive::Cubix;
using jive::IdxVector;
using jive::Matrix;
using jive::Ref;
using jive::Vector;
using jive::util::FuncUtils;
 
namespace jive_helpers
{
typedef Array<double, 4> Quadix;
typedef Array<IdxVector, 2> IdxVectorMatrix;
 
const Vector e1 = {1., 0., 0.};
const Vector e2 = {0., 1., 0.};
const Vector e3 = {0., 0., 1.};
const Vector vec0 = {0., 0., 0.};
const double TINY = __DBL_EPSILON__ * 1e6;
 
Vector funcGrad(const Ref<Function> func, const Vector &args);
 
Vector evalFuncs(const FuncUtils::FuncArray funcs, const Vector &args);
 
Matrix funcHessian(const Ref<Function> func, const Vector &args);
 
Matrix gradFuncs(const FuncUtils::FuncArray funcs, const Vector &args);
 
Matrix eye(const idx_t dim = 3);
 
void logMat(const Vector &rv, const Matrix &R);
 
void vec2mat(const Matrix &mat, const Vector &vec);
 
void mat2vec(const Vector &vec, const Matrix &mat);
 
// -----------------------------------------------------------------------
// Temporary arrays: why and how we avoid them
//
// A JEM array (Vector, Matrix, Cubix) is a small handle that points to a
// block of numbers. Copying a handle ("Vector a = b;") is cheap and does NOT
// copy the numbers: a and b then show the same numbers. Creating a new array,
// or writing "x = matmul(A, v)", asks the memory manager for a fresh block of
// numbers and gives it back afterwards. For the 3 to 6 numbers of a rod
// element this bookkeeping costs more than the arithmetic itself, and in the
// force assembly it happens millions of times per run. Therefore the code that
// is executed for every element and every time step ("hot code") follows
// three rules; everywhere else (set-up, output, rarely used code) the short
// forms matmul(A, v) and skew(v) stay, because they are easier to read.
//
//  1. Work arrays: arrays that live as (mutable) members of the class, are
//     created once in the constructor/configure and overwritten in every
//     call. A work array is always written completely before it is read
//     within the same call, and it is never used to carry values from one
//     call to the next. A class with work arrays must not be used by several
//     threads at once.
//  2. 3 x 3 products: mul3(C, A, B) and mul3(y, A, x). The result may be one
//     of the factors.
//  3. All other matrix products: a MatmulChain member (jem::numeric), which
//     keeps its result array and reuses it. The result is only valid until
//     the next call of the same chain, so use it directly or copy it first.
// -----------------------------------------------------------------------
 
void mul3(const Matrix &C, const Matrix &A, const Matrix &B);
 
void mul3(const Vector &y, const Matrix &A, const Vector &x);
 
void expVec(const Matrix &Exp, const Vector &v);
 
void expVecP(const Matrix &ExpP, const Vector &v, const Vector &vP);
 
double trace(const Matrix &mat);
 
idx_t factorial(const idx_t n);
 
idx_t binom(const idx_t n, const idx_t k);
 
double matrixNorm2(const Matrix &mat);
 
Matrix skew(const Vector &vec);
 
void skew(const Matrix &S, const Vector &vec);
 
Vector unskew(const Matrix &mat);
 
void setDofValue(jive::util::Constraints &cons, const idx_t idof, const double rval);
 
void linkDofs(jive::util::Constraints &cons, const idx_t idof, const idx_t jdof, const double coeff);
}; // namespace jive_helpers