Removed variable index from ForwardBlock class.

Feature was in the way of having (for example) face variables depending on cell variables.
This commit is contained in:
Atgeirr Flø Rasmussen 2013-05-01 09:03:56 +02:00
parent 633de102c2
commit 92f56f6d92

View File

@ -39,7 +39,7 @@ namespace AutoDiff
typedef Eigen::SparseMatrix<Scalar> M; typedef Eigen::SparseMatrix<Scalar> M;
/// Named constructor pattern used here. /// Named constructor pattern used here.
static ForwardBlock constant(const int index, const V& val, const std::vector<int>& blocksizes) static ForwardBlock constant(const V& val, const std::vector<int>& blocksizes)
{ {
std::vector<M> jac; std::vector<M> jac;
const int num_elem = val.size(); const int num_elem = val.size();
@ -49,12 +49,11 @@ namespace AutoDiff
for (int i = 0; i < num_blocks; ++i) { for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]); jac[i] = M(num_elem, blocksizes[i]);
} }
return ForwardBlock(index, val, jac); return ForwardBlock(val, jac);
} }
static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes) static ForwardBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
{ {
std::vector<M> jac; std::vector<M> jac;
const int num_elem = val.size(); const int num_elem = val.size();
const int num_blocks = blocksizes.size(); const int num_blocks = blocksizes.size();
@ -64,22 +63,22 @@ namespace AutoDiff
jac[i] = M(num_elem, blocksizes[i]); jac[i] = M(num_elem, blocksizes[i]);
} }
// ... then set the one corrresponding to this variable to identity. // ... then set the one corrresponding to this variable to identity.
assert(blocksizes[index] == num_elem);
jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1)); jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1));
for (typename M::Index row = 0; row < val.size(); ++row) { for (typename M::Index row = 0; row < val.size(); ++row) {
jac[index].insert(row, row) = Scalar(1.0); jac[index].insert(row, row) = Scalar(1.0);
} }
return ForwardBlock(index, val, jac); return ForwardBlock(val, jac);
} }
static ForwardBlock function(const int index, const V& val, const std::vector<M>& jac) static ForwardBlock function(const V& val, const std::vector<M>& jac)
{ {
return ForwardBlock(index, val, jac); return ForwardBlock(val, jac);
} }
/// Operators. /// Operator +
ForwardBlock operator+(const ForwardBlock& rhs) ForwardBlock operator+(const ForwardBlock& rhs)
{ {
assert(index() == rhs.index());
std::vector<M> jac = jac_; std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks(); int num_blocks = numBlocks();
@ -88,13 +87,12 @@ namespace AutoDiff
assert(jac[block].cols() == rhs.jac_[block].cols()); assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] += rhs.jac_[block]; jac[block] += rhs.jac_[block];
} }
return function(index(), val_ + rhs.val_, jac); return function(val_ + rhs.val_, jac);
} }
/// Operators. /// Operator -
ForwardBlock operator-(const ForwardBlock& rhs) ForwardBlock operator-(const ForwardBlock& rhs)
{ {
assert(index() == rhs.index());
std::vector<M> jac = jac_; std::vector<M> jac = jac_;
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
int num_blocks = numBlocks(); int num_blocks = numBlocks();
@ -103,13 +101,12 @@ namespace AutoDiff
assert(jac[block].cols() == rhs.jac_[block].cols()); assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] -= rhs.jac_[block]; jac[block] -= rhs.jac_[block];
} }
return function(index(), val_ - rhs.val_, jac); return function(val_ - rhs.val_, jac);
} }
/// Operators. /// Operator *
ForwardBlock operator*(const ForwardBlock& rhs) ForwardBlock operator*(const ForwardBlock& rhs)
{ {
assert(index() == rhs.index());
int num_blocks = numBlocks(); int num_blocks = numBlocks();
std::vector<M> jac(num_blocks); std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
@ -122,13 +119,12 @@ namespace AutoDiff
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
jac[block] = D2*jac_[block] + D1*rhs.jac_[block]; jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
} }
return function(index(), val_ * rhs.val_, jac); return function(val_ * rhs.val_, jac);
} }
/// Operators. /// Operator /
ForwardBlock operator/(const ForwardBlock& rhs) ForwardBlock operator/(const ForwardBlock& rhs)
{ {
assert(index() == rhs.index());
int num_blocks = numBlocks(); int num_blocks = numBlocks();
std::vector<M> jac(num_blocks); std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
@ -137,12 +133,11 @@ namespace AutoDiff
D D2 = rhs.val_.matrix().asDiagonal(); D D2 = rhs.val_.matrix().asDiagonal();
D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal(); D D3 = std::pow(rhs.val_, -2).matrix().asDiagonal();
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
std::cout << jac[block].rows() << ' ' << rhs.jac_[block].rows() << std::endl;
assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]); jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
} }
return function(index(), val_ / rhs.val_, jac); return function(val_ / rhs.val_, jac);
} }
/// I/O. /// I/O.
template <class Ostream> template <class Ostream>
@ -150,41 +145,43 @@ namespace AutoDiff
print(Ostream& os) const print(Ostream& os) const
{ {
int num_blocks = jac_.size(); int num_blocks = jac_.size();
os << "Variable index " << index_ << ":\nValue =\n" << val_ << "\n\nJacobian =\n"; os << "Value =\n" << val_ << "\n\nJacobian =\n";
for (int i = 0; i < num_blocks; ++i) { for (int i = 0; i < num_blocks; ++i) {
os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n"; os << "Sub Jacobian #" << i << '\n' << jac_[i] << "\n";
} }
return os; return os;
} }
/// Index of this variable.
int index() const
{
return index_;
}
/// Number of variables or Jacobian blocks. /// Number of variables or Jacobian blocks.
int numBlocks() const int numBlocks() const
{ {
return jac_.size(); return jac_.size();
} }
/// Function value
const V& value() const
{
return val_;
}
/// Function derivatives
const std::vector<M>& derivative() const
{
return jac_;
}
private: private:
ForwardBlock(const int index, ForwardBlock(const V& val,
const V& val,
const std::vector<M>& jac) const std::vector<M>& jac)
: index_(index), val_(val), jac_(jac) : val_(val), jac_(jac)
{ {
const int num_elem = val_.size(); const int num_elem = val_.size();
const int num_blocks = jac_.size(); const int num_blocks = jac_.size();
assert(index_ < num_blocks);
assert(num_elem == jac_[index].cols());
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
assert(num_elem == jac_[block].rows()); assert(num_elem == jac_[block].rows());
} }
} }
int index_;
V val_; V val_;
std::vector<M> jac_; std::vector<M> jac_;
}; };