b2element_klt_util.H

//------------------------------------------------------------------------
// b2element_klt_util.H --
//
//     Data structures for scalar and vector quantities with
//     derivatives by the dof and by the element-internal directions.
//
// written by Thomas Ludwig
//            Harald Klimach <harald.klimach@dlr.de>
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//
// Copyright (c) 2015,2016,2017 SMR Engineering & Development SA
//                              2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2ELEMENT_KLT_UTIL_H__
#define __B2ELEMENT_KLT_UTIL_H__
 
// #define KLT_WITH_HIGH_ORDER_ELEMENTS 1
 
#include <iostream>
#include <vector>
 
#include "b2ppconfig.h"
#include "utils/b2linear_algebra.H"
#include "utils/b2tensor_calculus.H"
#include "utils/b2vec3.H"
 
namespace b2000::klt {
 
enum {
    DOF_TYPE_UNKNOWN = 0,
    DOF_TYPE_DX = 1,
    DOF_TYPE_DY = 2,
    DOF_TYPE_DZ = 3,
    DOF_TYPE_DW = 4,
    DOF_TYPE_RX = 5,
    DOF_TYPE_RY = 6,
    DOF_TYPE_RZ = 7
};
 
enum {
    SUB_ALL = 0,
    SUB_EDGE_C0 = 1,
    SUB_EDGE_C1 = 2,
    SUB_EDGE_C2 = 3,
    SUB_VERTEX_C0 = 4,
    SUB_VERTEX_C1 = 5,
    SUB_VERTEX_C2 = 6
};
 
enum {
    GEOM_ISOPARAMETRIC,
    GEOM_ANALYTIC_CYLINDER,
    GEOM_ANALYTIC_HYPERBOLIC_PARABOLOID,
    GEOM_ANALYTIC_STRAIGHT_EDGES
};
 
class IntegrationScheme {
public:
    size_t get_num() const { return int(points.size()); }
 
    void get_point(const size_t ip, double xi[2], double& weight) const {
        const Point& p = points[ip];
        std::copy(p.xi, p.xi + 2, xi);
        weight = p.weight;
    }
 
    struct Point {
        double xi[2];
        double weight;
 
        Point(const double xi_[2], const double weight_) : weight(weight_) {
            std::copy(xi_, xi_ + 2, xi);
        }
 
        Point(const double r, const double s, const double weight_) : weight(weight_) {
            xi[0] = r;
            xi[1] = s;
        }
    };
 
    typedef std::vector<Point> PointVector;
 
    PointVector points;
};
 
template <typename T>
class Array2 {
public:
    Array2() : size1_(0), size2_(0), m2_(0), m3_(0) {}
 
    Array2(size_t size1, size_t size2, const T& t = T())
        : size1_(size1), size2_(size2), m2_(size1_), m3_(m2_ * size2_), data_(m3_, t) {}
 
    void set_zero() { std::fill(data_.begin(), data_.end(), T()); }
 
    void swap(Array2& o) {
        std::swap(size1_, o.size1_);
        std::swap(size2_, o.size2_);
        std::swap(m2_, o.m2_);
        std::swap(m3_, o.m3_);
        data_.swap(o.data_);
    }
 
    template <typename T1>
    Array2& operator=(const Array2<T1>& o) {
        size1_ = o.size1_;
        size2_ = o.size2_;
        m2_ = o.m2_;
        m3_ = o.m3_;
        data_.resize(m3_);
        for (size_t i = 0; i != m3_; ++i) { data_[i] = o.data_[i]; }
        return *this;
    }
 
    Array2& operator=(const Array2& o) {
        size1_ = o.size1_;
        size2_ = o.size2_;
        m2_ = o.m2_;
        m3_ = o.m3_;
        data_ = o.data_;
        return *this;
    }
 
    size_t size1() const { return size1_; }
 
    size_t size2() const { return size2_; }
 
    size_t msize1() const { return 1; }
 
    size_t msize2() const { return m2_; }
 
    size_t msize3() const { return m3_; }
 
    size_t esize() const { return m3_; }
 
    void resize(size_t size1, size_t size2, const T& t = T()) {
        size1_ = size1;
        size2_ = size2;
        m2_ = size1_;
        m3_ = m2_ * size2_;
        data_.resize(m3_, t);
    }
 
    size_t index(size_t i1, size_t i2) const { return i1 * m2_ + i2; }
 
    const T& operator()(size_t i1, size_t i2) const { return data_[index(i1, i2)]; }
 
    T& operator()(size_t i1, size_t i2) { return data_[index(i1, i2)]; }
 
    const T& operator[](size_t i) const { return data_[i]; }
 
    T& operator[](size_t i) { return data_[i]; }
 
private:
    size_t size1_;
    size_t size2_;
    size_t m2_;
    size_t m3_;
    std::vector<T> data_;
};
 
struct DofScalar {
    size_t ndof;
    bool first;
    bool second;
    double d0;               // value
    std::vector<double> d1;  // first derivatives
    Array2<double> d2;       // second derivatives
 
    DofScalar() : ndof(0), first(false), second(false), d0(0.) {}
 
    DofScalar(const size_t ndof_, const bool first_, const bool second_)
        : ndof(ndof_),
          first(first_),
          second(second_),
          d0(0.),
          d1((first ? ndof : 0), 0.),
          d2((second ? ndof : 0), (second ? ndof : 0), 0.) {}
 
    void resize(const size_t ndof_, const bool first_, const bool second_) {
        first = first_;
        second = second_;
        d1.resize((first ? ndof_ : 0), 0.);
        d2.resize((second ? ndof_ : 0), (second ? ndof_ : 0), 0.);
        ndof = ndof_;
    }
 
    size_t size() const { return ndof; }
 
    DofScalar& set_zero() {
        d0 = 0.;
        if (first) { std::fill(d1.begin(), d1.end(), 0); }
        if (second) { d2.set_zero(); }
        return *this;
    }
 
    DofScalar& operator*=(const double o) {
        d0 *= o;
        for (size_t i = 0; i != d1.size(); ++i) { d1[i] *= o; }
        for (size_t i = 0; i != d2.esize(); ++i) { d2[i] *= o; }
        return *this;
    }
 
    DofScalar& operator+=(const DofScalar& o) {
        assert(size() == o.size());
        d0 += o.d0;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] += o.d1[i]; }
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] += o.d2[i]; }
        }
        return *this;
    }
 
    DofScalar& assign(const DofScalar& o) {
        assert(ndof == o.ndof);
        d0 = o.d0;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] = o.d1[i]; }
        } else if (first) {
            std::fill(d1.begin(), d1.end(), 0);
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] = o.d2[i]; }
        } else if (second) {
            d2.set_zero();
        }
        return *this;
    }
 
    DofScalar& muladd(const DofScalar& o, const double f) {
        assert(ndof == o.ndof);
        d0 += o.d0 * f;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] += f * o.d1[i]; }
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] += f * o.d2[i]; }
        }
        return *this;
    }
 
    DofScalar& add_at_offset(const size_t offset, const DofScalar& o) {
        if (ndof < offset + o.ndof) { resize(offset + o.ndof, first, second); }
        d0 += o.d0;
        for (size_t k = 0; k != o.ndof; ++k) { d1[offset + k] = o.d1[k]; }
        for (size_t k = 0; k != o.ndof; ++k) {
            for (size_t l = 0; l <= k; ++l) { d2(offset + k, offset + l) = o.d2(k, l); }
        }
        return *this;
    }
};
 
//
// Vector quantities with derivatives by the dof (up to 2x).
//
template <int NUMD>
struct DofVector {
    DofScalar e[NUMD];
 
    DofVector() {}
 
    DofVector(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int j = 0; j != NUMD; ++j) { e[j].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    DofVector& set_zero() {
        for (int j = 0; j != NUMD; ++j) { e[j].set_zero(); }
        return *this;
    }
 
    const DofScalar& operator[](const size_t i) const { return e[i]; }
 
    DofScalar& operator[](const size_t i) { return e[i]; }
 
    DofVector& operator*=(const double o) {
        for (int j = 0; j != NUMD; ++j) { e[j] *= o; }
        return *this;
    }
 
    DofVector& operator+=(const DofVector& o) {
        for (int j = 0; j != NUMD; ++j) { e[j] += o[j]; }
        return *this;
    }
 
    DofVector& add_at_offset(const size_t offset, const DofVector& o) {
        for (int j = 0; j != NUMD; ++j) { e[j].add_at_offset(offset, o[j]); }
        return *this;
    }
};
 
//
// Matrix (3x3) quantities with derivatives by the dof (up to 2x).
//
struct DofMatrix {
    DofVector<3> e[3];
 
    DofMatrix() {}
 
    DofMatrix(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int j = 0; j != 3; ++j) { e[j].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    DofMatrix& set_zero() {
        for (int j = 0; j != 3; ++j) { e[j].set_zero(); }
        return *this;
    }
 
    const DofVector<3>& operator[](const size_t i) const { return e[i]; }
 
    DofVector<3>& operator[](const size_t i) { return e[i]; }
 
    DofMatrix& operator*=(const double o) {
        for (int j = 0; j != 3; ++j) { e[j] *= o; }
        return *this;
    }
 
    DofMatrix& operator+=(const DofMatrix& o) {
        for (int j = 0; j != 3; ++j) { e[j] += o[j]; }
        return *this;
    }
 
    DofMatrix& add_at_offset(const size_t offset, const DofMatrix& o) {
        for (int j = 0; j != 3; ++j) { e[j].add_at_offset(offset, o[j]); }
        return *this;
    }
};
 
// Indices for derivatives by the element-local system directions.
const int d00 = 0;
const int d10 = 1;
const int d01 = 2;
const int d20 = 3;
const int d11 = 4;
const int d02 = 5;
const int d30 = 6;
const int d21 = 7;
const int d12 = 8;
const int d03 = 9;
 
//
// Scalar quantities with derivatives by the element-local system
// directions (up to 3x), and derivatives by the dof (up to 2x).
//
template <int NUMD>
struct ScalarXiDof {
    DofScalar e[NUMD];
 
    ScalarXiDof() {}
 
    ScalarXiDof(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int i = 0; i != NUMD; ++i) { e[i].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    ScalarXiDof& set_zero() {
        for (int i = 0; i != NUMD; ++i) { e[i].set_zero(); }
        return *this;
    }
 
    const DofScalar& operator[](const size_t i) const {
        // assert(i < NUMD);
        return e[i];
    }
 
    DofScalar& operator[](const size_t i) {
        // assert(i < NUMD);
        return e[i];
    }
 
    ScalarXiDof& operator*=(const double o) {
        for (int i = 0; i != NUMD; ++i) { e[i] *= o; }
        return *this;
    }
 
    ScalarXiDof& operator+=(const ScalarXiDof& o) {
        for (int i = 0; i != NUMD; ++i) { e[i] += o[i]; }
        return *this;
    }
 
    ScalarXiDof& muladd(const ScalarXiDof& o, const double f) {
        for (int i = 0; i != NUMD; ++i) { e[i].muladd(o[i], f); }
        return *this;
    }
 
    ScalarXiDof& add_at_offset(const size_t offset, const ScalarXiDof& o) {
        for (int i = 0; i != NUMD; ++i) { e[i].add_at_offset(offset, o[i]); }
        return *this;
    }
};
 
//
// Vector quantities with derivatives by the element-local system
// directions (up to 3x), and derivatives by the dof (up to 2x).
//
template <int NUMD>
struct VectorXiDof {
    ScalarXiDof<NUMD> e[3];
 
    VectorXiDof() {}
 
    VectorXiDof(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int j = 0; j != 3; ++j) { e[j].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    VectorXiDof& set_zero() {
        for (int j = 0; j != 3; ++j) { e[j].set_zero(); }
        return *this;
    }
 
    const ScalarXiDof<NUMD>& operator[](const size_t i) const {
        // assert(i < 3);
        return e[i];
    }
 
    ScalarXiDof<NUMD>& operator[](const size_t i) {
        // assert(i < 3);
        return e[i];
    }
 
    VectorXiDof& operator*=(const double o) {
        for (int j = 0; j != 3; ++j) { e[j] *= o; }
        return *this;
    }
 
    VectorXiDof& operator+=(const VectorXiDof& o) {
        for (int j = 0; j != 3; ++j) { e[j] += o[j]; }
        return *this;
    }
 
    VectorXiDof& add_at_offset(const size_t offset, const VectorXiDof& o) {
        for (int j = 0; j != 3; ++j) { e[j].add_at_offset(offset, o[j]); }
        return *this;
    }
};
 
//
// Matrix (3x3) quantities with derivatives by the element-local
// system directions (up to 3x), and derivatives by the dof (up to
// 2x).
//
template <int NUMD>
struct MatrixXiDof {
    VectorXiDof<NUMD> e[3];
 
    MatrixXiDof() {}
 
    MatrixXiDof(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int j = 0; j != 3; ++j) { e[j].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    MatrixXiDof& set_zero() {
        for (int j = 0; j != 3; ++j) { e[j].set_zero(); }
        return *this;
    }
 
    const VectorXiDof<NUMD>& operator[](const size_t i) const { return e[i]; }
 
    VectorXiDof<NUMD>& operator[](const size_t i) { return e[i]; }
 
    MatrixXiDof& operator*=(const double o) {
        for (int j = 0; j != 3; ++j) { e[j] *= o; }
        return *this;
    }
 
    MatrixXiDof& operator+=(const MatrixXiDof& o) {
        for (int j = 0; j != 3; ++j) { e[j] += o[j]; }
        return *this;
    }
 
    MatrixXiDof& add_at_offset(const size_t offset, const MatrixXiDof& o) {
        for (int j = 0; j != 3; ++j) { e[j].add_at_offset(offset, o[j]); }
        return *this;
    }
};
 
typedef Vec3<double> Vec3d;
 
//
// Value and up to 2nd derivatives of 3d vectors.
//
struct DofVec3d {
    size_t ndof;
    bool first;
    bool second;
    Vec3d d0;               // value
    std::vector<Vec3d> d1;  // first derivatives
    Array2<Vec3d> d2;       // second derivatives
 
    DofVec3d() : ndof(0), first(false), second(false) {}
 
    DofVec3d(const size_t ndof_, const bool first_, const bool second_)
        : ndof(ndof_),
          first(first_),
          second(second_),
          d1((first ? ndof : 0), Vec3d()),
          d2((second ? ndof : 0), (second ? ndof : 0), Vec3d()) {}
 
    void resize(const size_t ndof_, const bool first_, const bool second_) {
        first = first_;
        second = second_;
        d1.resize((first ? ndof_ : 0), Vec3d());
        d2.resize((second ? ndof_ : 0), (second ? ndof_ : 0), Vec3d());
        ndof = ndof_;
    }
 
    size_t size() const { return ndof; }
 
    DofVec3d& set_zero() {
        d0 = Vec3d();
        if (first) { std::fill(d1.begin(), d1.end(), Vec3d()); }
        if (second) { d2.set_zero(); }
        return *this;
    }
 
    DofVec3d& operator*=(const double o) {
        d0 *= o;
        for (size_t i = 0; i != d1.size(); ++i) { d1[i] *= o; }
        for (size_t i = 0; i != d2.esize(); ++i) { d2[i] *= o; }
        return *this;
    }
 
    DofVec3d& operator+=(const DofVec3d& o) {
        assert(size() == o.size());
        d0 += o.d0;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] += o.d1[i]; }
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] += o.d2[i]; }
        }
        return *this;
    }
 
    DofVec3d& assign(const DofVec3d& o) {
        assert(ndof == o.ndof);
        d0 = o.d0;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] = o.d1[i]; }
        } else if (first) {
            for (size_t i = 0; i != d1.size(); ++i) { d1[i].set_zero(); }
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] = o.d2[i]; }
        } else if (second) {
            d2.set_zero();
        }
        return *this;
    }
 
    DofVec3d& muladd(const DofVec3d& o, const double f) {
        assert(ndof == o.ndof);
        d0 += o.d0 * f;
        if (first && o.first) {
            assert(d1.size() == o.d1.size());
            for (size_t i = 0; i != d1.size(); ++i) { d1[i] += f * o.d1[i]; }
        }
        if (second && o.second) {
            assert(d2.esize() == o.d2.esize());
            for (size_t i = 0; i != d2.esize(); ++i) { d2[i] += f * o.d2[i]; }
        }
        return *this;
    }
};
 
//
// Vec3 quantities with derivatives by the element-local system
// directions (up to 3x), and derivatives by the dof (up to 2x).
//
template <int NUMD>
struct XiDofVec3d {
    DofVec3d e[NUMD];
 
    XiDofVec3d() {}
 
    XiDofVec3d(const size_t ndof, const bool first, const bool second) {
        resize(ndof, first, second);
    }
 
    void resize(const size_t ndof, const bool first, const bool second) {
        for (int i = 0; i != NUMD; ++i) { e[i].resize(ndof, first, second); }
    }
 
    size_t size() const { return e[0].size(); }
 
    XiDofVec3d& set_zero() {
        for (int i = 0; i != NUMD; ++i) { e[i].set_zero(); }
        return *this;
    }
 
    const DofVec3d& operator[](const size_t i) const {
        // assert(i < NUMD);
        return e[i];
    }
 
    DofVec3d& operator[](const size_t i) {
        // assert(i < NUMD);
        return e[i];
    }
 
    XiDofVec3d& operator*=(const double o) {
        for (int i = 0; i != NUMD; ++i) { e[i] *= o; }
        return *this;
    }
 
    XiDofVec3d& operator+=(const XiDofVec3d& o) {
        for (int i = 0; i != NUMD; ++i) { e[i] += o[i]; }
        return *this;
    }
 
    XiDofVec3d& muladd(const XiDofVec3d& o, const double f) {
        for (int i = 0; i != NUMD; ++i) { e[i].muladd(o[i], f); }
        return *this;
    }
 
    XiDofVec3d& add_at_offset(const size_t offset, const XiDofVec3d& o) {
        for (int i = 0; i != NUMD; ++i) { e[i].add_at_offset(offset, o[i]); }
        return *this;
    }
};
 
}  // namespace b2000::klt
 
#endif  // B2_ELEMENT_KLT_UTIL_H_