b2element.H

//------------------------------------------------------------------------
// b2element.H --
//
//
// written by Mathias Doreille
//            Thomas Blome <thomas.blome@dlr.de>
//
// Copyright (c) 2004-2012,2016,2017 SMR Engineering & Development SA
//                                   2502 Bienne, Switzerland
//
// Copyright (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
//                    Linder Höhe, 51147 Köln
//
// 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_H_
#define B2ELEMENT_H_
 
#include <type_traits>
 
#include "b2ppconfig.h"
#include "model/b2degrees_of_freedom.H"
#include "utils/b2csda.H"
#include "utils/b2exception.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2object.H"
#include "utils/b2type.H"
 
namespace b2000 {
 
class Node;
class Case;
class GradientContainer;
class SolverHints;
class Model;
class Domain;
class RTable;
 
namespace solver {
template <typename T, typename MATRIX_FORMAT>
class TypedSolver;
}
 
namespace element {
extern Module module;
}
 
class ElementProperty : public Object {
public:
    using type_t = ObjectTypeIncomplete<ElementProperty>;
    static type_t type;
 
    virtual void set_model_and_case(Model& model, Case& case_) {}
 
    virtual ElementProperty* copy() {
        UnimplementedError() << THROW;
        return nullptr;
    }
 
    const std::string& get_object_name() const override {
        static const std::string defaut = "<unknown object name>";
        return defaut;
    }
 
    virtual void set_object_name(const std::string& name) {}
};
 
template <typename T>
class Field : public Object {
public:
    virtual void get_value(
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& lcoor,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& coor,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& coor_deformed,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& d_coor_d_lcoor,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& d_coor_deformed_d_lcoor,
          b2linalg::Vector<T, b2linalg::Vdense>& value) const {
        UnimplementedError() << THROW;
    }
 
    virtual bool get_multiply_with_density() const { return false; }
};
 
class Element : public DegreesOfFreedom {
public:
    using type_t = ObjectTypeIncomplete<Element>;
    Element() = default;
 
    explicit Element(const bool symmetric) : symmetric_{symmetric} {}
 
 
    ~Element() override = default;
 
    Element(const Element&) = delete;
    Element(const Element&&) = delete;
    Element& operator=(const Element&) = delete;
    Element& operator=(const Element&&) = delete;
 
    size_t get_id() const { return id; }
 
    void set_id(size_t id_) { id = id_; }
 
    int get_elem_type_number() const { return elno_; }
 
    void set_elem_type_number(int elno) { elno_ = elno; }
 
    const std::string& get_object_name() const override { return object_name; }
 
    void set_object_name(const std::string& name) { object_name = name; }
 
 
    int get_number_of_dof() const override { return 0; }
 
    size_t set_global_dof_numbering(size_t index) override { return index; }
 
    std::pair<size_t, size_t> get_global_dof_numbering() const override {
        return std::pair<size_t, size_t>(0, 0);
    }
 
    [[deprecated]] virtual void get_dof_numbering(b2linalg::Index& dof_numbering) {
        dof_numbering.resize(0);
    }
 
 
    virtual std::pair<int, Node* const*> get_nodes() const {
        return std::pair<int, Node* const*>(0, nullptr);
    }
 
    virtual void set_nodes(std::pair<int, Node* const*> nodes) {}
 
    virtual const ElementProperty* get_property() const { return nullptr; }
 
    virtual void set_property(ElementProperty* property) {}
 
    virtual void set_additional_properties(const RTable& rtable) {}
 
    virtual void init(Model& model) {}
 
 
    virtual std::pair<size_t, Element**> get_subelements() {
        return std::pair<size_t, Element**>(0, nullptr);
    }
 
    class iterator {
        friend class Element;
 
    public:
        Element& operator*() {
            if (element != nullptr) { return *element; }
            return **sstart;
        }
 
        const Element& operator*() const {
            if (element != nullptr) { return *element; }
            return **sstart;
        }
 
        void operator++() {
            if (element != nullptr) {
                std::pair<size_t, Element**> s = element->get_subelements();
                if (s.second != 0) {
                    sstart = s.second;
                    send = s.second + s.first;
                }
                element = nullptr;
            } else if (sstart < send) {
                if (++sstart == send) { sstart = send = nullptr; }
            }
        }
 
        bool operator==(const iterator& o) const {
            return element == o.element && sstart == o.sstart;
        }
 
        bool operator!=(const iterator& o) const {
            return element != o.element || sstart != o.sstart;
        }
 
    private:
        iterator() : element(nullptr), sstart(nullptr), send(nullptr) {}
 
        iterator(Element& element_) : element(&element_), sstart(nullptr), send(nullptr) {}
 
        Element* element;
        Element** sstart;
        Element** send;
    };
 
    iterator begin() { return iterator(*this); }
 
    iterator end() const { return iterator(); }
 
 
    enum VariableInfo {
        zero = 0,
 
        constant = 1,
 
        linear = 2,
 
        nonlinear = 3
    };
 
    virtual const std::vector<VariableInfo> get_value_info() const {
        return std::vector<VariableInfo>();
    }
 
    virtual std::pair<int, VariableInfo> get_constraint_info() {
        return std::pair<int, VariableInfo>(0, zero);
    }
 
    virtual int edge_field_order(const int edge_id, const std::string& field_name) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual bool edge_field_linear_on_dof(const int edge_id, const std::string& field_name) {
        UnimplementedError() << THROW;
        return false;
    }
 
    virtual int edge_field_polynomial_sub_edge(
          const int edge_id, const std::string& field_name,
          b2linalg::Vector<double, b2linalg::Vdense>& sub_nodes) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual void edge_geom(
          const int edge_id, const double internal_coor, b2linalg::Vector<double>& geom,
          b2linalg::Vector<double>& d_geom_d_icoor) {
        UnimplementedError() << THROW;
    }
 
    virtual int face_field_order(const int face_id, const std::string& field_name) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual bool face_field_linear_on_dof(const int face_id, const std::string& field_name) {
        UnimplementedError() << THROW;
        return false;
    }
 
    struct Triangle {
        Triangle(int n1, int n2, int n3) {
            node[0] = n1;
            node[1] = n2;
            node[2] = n3;
        }
        size_t node[3];
    };
 
    virtual int face_field_polynomial_sub_face(
          const int face_id, const std::string& field_name,
          b2linalg::Matrix<double, b2linalg::Mrectangle>& sub_nodes,
          std::vector<Triangle>& sub_faces) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual void face_geom(
          const int face_id,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
          b2linalg::Vector<double>& geom,
          b2linalg::Matrix<double, b2linalg::Mrectangle>& d_geom_d_icoor) {
        UnimplementedError() << THROW;
    }
 
 
    virtual int body_field_order(const std::string& field_name) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual bool body_field_linear_on_dof(const std::string& field_name) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    struct Tetrahedral {
        int node[4];
    };
 
    virtual int body_field_polynomial_sub_volume(
          const std::string& field_name, b2linalg::Matrix<double, b2linalg::Mrectangle>& sub_nodes,
          std::vector<Tetrahedral>& sub_volumes) {
        UnimplementedError() << THROW;
        return 0;
    }
 
    virtual void body_geom(
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
          b2linalg::Vector<double>& geom,
          b2linalg::Matrix<double, b2linalg::Mrectangle>& d_geom_d_icoor) {
        UnimplementedError() << THROW;
    }
 
    virtual std::pair<int, int> get_state_buffer_size() const { return std::pair<int, int>(0, 0); }
 
    virtual std::pair<int*, double*> get_state_buffer(std::pair<int*, double*> buffers) const {
        return buffers;
    }
 
    virtual std::pair<const int*, const double*> set_state_buffer(
          std::pair<const int*, const double*> buffers) {
        return buffers;
    }
 
    bool is_symmetric() const { return symmetric_; }
 
 
    static std::vector<bool> d_value_d_dof_flags_null;
 
    static type_t type;
 
protected:
    const bool symmetric_{true};
 
private:
    size_t id{};
    std::string object_name{"???"};
    int elno_;
};
 
template <typename T>
class TypedElement : public Element {
public:
    using type_t = ObjectTypeIncomplete<TypedElement<T>, Element::type_t>;
 
    TypedElement() = default;
 
    explicit TypedElement(const bool symmetric) : Element{symmetric} {}
 
    virtual void get_constraint(
          Model& model, const bool linear, const double time,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
          b2linalg::Index& dof_numbering, b2linalg::Vector<T, b2linalg::Vdense>& constraint,
          b2linalg::Matrix<T, b2linalg::Mcompressed_col>& trans_d_constraint_d_dof,
          b2linalg::Vector<T, b2linalg::Vdense>& d_constraint_d_time) {
        UnimplementedError() << THROW;
    }
 
    virtual void edge_field_value(
          const int edge_id, const std::string& field_name, const double internal_coor,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
          b2linalg::Vector<T, b2linalg::Vdense>& value,
          b2linalg::Vector<T, b2linalg::Vdense>& d_value_d_icoor, b2linalg::Index& dof_numbering,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& d_value_d_dof,
          b2linalg::Index& d_value_d_dof_dep_col = b2linalg::Index::null) {
        UnimplementedError() << THROW;
    }
 
    virtual void face_field_value(
          const int face_id, const std::string& field_name,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
          b2linalg::Vector<T, b2linalg::Vdense>& value,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& d_value_d_icoor,
          b2linalg::Index& dof_numbering, b2linalg::Matrix<T, b2linalg::Mrectangle>& d_value_d_dof,
          b2linalg::Index& d_value_d_dof_dep_col) {
        UnimplementedError() << THROW;
    }
 
    virtual void body_field_value(
          const std::string& field_name,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& internal_coor,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, const double time,
          b2linalg::Vector<T, b2linalg::Vdense>& value,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& d_value_d_icoor,
          b2linalg::Index& dof_numbering,
          b2linalg::Matrix<T, b2linalg::Mrectangle>& d_value_d_dof) {
        UnimplementedError() << THROW;
    }
 
    virtual void field_volume_integration(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdot,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdotdot, const Field<T>& f,
          b2linalg::Index& dof_numbering, b2linalg::Vector<T, b2linalg::Vdense>& discretised_field,
          b2linalg::Matrix<T, b2linalg::Mpacked>& d_discretised_field_d_dof) {
        if (!dof_numbering.is_null()) { dof_numbering.resize(0); }
        if (!discretised_field.is_null()) { discretised_field.resize(0); }
        if (!d_discretised_field_d_dof.is_null()) { d_discretised_field_d_dof.resize(0); }
    }
 
    virtual void field_face_integration(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdot,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdotdot, const Field<T>& field,
          int face, b2linalg::Index& dof_numbering,
          b2linalg::Vector<T, b2linalg::Vdense>& discretised_field,
          b2linalg::Matrix<T, b2linalg::Mpacked>& d_discretised_field_d_dof) {
        if (!dof_numbering.is_null()) { dof_numbering.resize(0); }
        if (!discretised_field.is_null()) { discretised_field.resize(0); }
        if (!d_discretised_field_d_dof.is_null()) { d_discretised_field_d_dof.resize(0); }
    }
 
    virtual void field_edge_integration(
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdot,
          const b2linalg::Vector<T, b2linalg::Vdense_constref>& dofdotdot, const Field<T>& f,
          int edge, b2linalg::Index& dof_numbering,
          b2linalg::Vector<T, b2linalg::Vdense>& discretised_field,
          b2linalg::Matrix<T, b2linalg::Mpacked>& d_discretised_field_d_dof) {
        UnimplementedError() << THROW;
    }
 
    virtual void get_axe_aligned_bounding_box(
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof, int dim, double* min,
          double* max) {
        UnimplementedError() << THROW;
    }
 
    template <typename MATRIX_FORMAT>
    void AssembleElementEffectiveSystem(
          const solver::TypedSolver<T, MATRIX_FORMAT>& solver,
          b2linalg::Matrix<T, typename MATRIX_FORMAT::dense>& k_eff,
          b2linalg::Vector<T, b2linalg::Vdense>& f_eff, b2linalg::Index& index) {
        // Define variables to use deprecated interface get_value().
        b2000::Model& model{solver.GetModel()};  // Make this constant (TB).
        const bool is_linear{solver.IsLinear()};
        const EquilibriumSolution equilibrium_solution{false};  // Über solver abfragen (TB).
        const double time{solver.GetTime()};
        const double step_size{solver.GetStepSize()};
        const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof{solver.GetDof()};
 
        std::vector<bool> matrix_flags{
              true /*k*/, false /*d*/, false /*m*/};  
 
        if (solver.IsDamped()) { matrix_flags[1] = true; }
        if (solver.IsDynamic()) { matrix_flags[2] = true; }
 
        [[likely]] if (symmetric_) {
            std::vector<b2linalg::Matrix<T, b2linalg::Mpacked>> matrix_sym{matrix_flags.size()};
 
            this->get_value(
                  model, is_linear, equilibrium_solution, time, step_size, dof,
                  nullptr /*GradientContainer*/, nullptr /*SolverHints**/, index, f_eff,
                  matrix_flags, matrix_sym, b2linalg::Vector<T, b2linalg::Vdense>::null);
 
            if (solver.IsDamped() || solver.IsDynamic()) {
                auto [mat_dyn, vec_dyn] = solver.GetTimeIntegrator()->GetEffectiveElementSystem(
                      matrix_sym[0], matrix_sym[1], matrix_sym[2], dof, index);
 
                k_eff = mat_dyn;
                // f_eff += vec_dyn;
            } else {
                k_eff = matrix_sym[0];  
                // k_eff += matrix_sym[0];   //!< Return proper stiffness matrix.
            }
        } else {
            std::vector<b2linalg::Matrix<T, b2linalg::Mrectangle>> matrix_unsym{
                  matrix_flags.size()};
 
            this->get_value(
                  model, is_linear, equilibrium_solution, time, step_size, dof,
                  nullptr /*GradientContainer*/, nullptr /*SolverHints**/, index, f_eff,
                  matrix_flags, matrix_unsym, b2linalg::Vector<T, b2linalg::Vdense>::null);
 
            if (solver.IsDamped() || solver.IsDynamic()) {
                auto [mat_dyn, vec_dyn] = solver.GetTimeIntegrator()->GetEffectiveElementSystem(
                      matrix_unsym[0], matrix_unsym[1], matrix_unsym[2], dof, index);
 
                k_eff = mat_dyn;
                // f_eff += vec_dyn;
            } else {
                k_eff = matrix_unsym[0];  
            }
        }
    };
 
    template <typename MATRIX_FORMAT>
    b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> AssembleElementEffectiveMatrix(
          const solver::TypedSolver<T, MATRIX_FORMAT>& solver, b2linalg::Index& index) {
        // Define variables to use deprecated interface get_value().
        b2000::Model& model{solver.GetModel()};  // Make this constant (TB).
        const bool is_linear{solver.IsLinear()};
        const EquilibriumSolution equilibrium_solution{false};  // Über solver abfragen (TB).
        const double time{solver.GetTime()};
        const double step_size{solver.GetStepSize()};
        const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof{solver.GetDof()};
 
        b2linalg::Matrix<T, typename MATRIX_FORMAT::dense> k_eff;
 
        std::vector<bool> stiffness{true};
 
        // If the element is symmetric use local matrix_sym matrix.
        [[likely]] if (symmetric_) {
            std::vector<b2linalg::Matrix<T, b2linalg::Mpacked>> matrix_sym{stiffness.size()};
 
            this->get_value(
                  model, is_linear, equilibrium_solution, time, step_size, dof,
                  nullptr /*GradientContainer*/, nullptr /*SolverHints**/, index,
                  b2linalg::Vector<T, b2linalg::Vdense>::null,
                  TypedElement<T>::d_value_d_dof_flags_null, matrix_sym,
                  b2linalg::Vector<T, b2linalg::Vdense>::null);
 
            k_eff = matrix_sym[0];
        } else {
            // If the element is unsymmetric use matrix_unsym.
            std::vector<b2linalg::Matrix<T, b2linalg::Mrectangle>> matrix_unsym{stiffness.size()};
 
            this->get_value(
                  model, is_linear, equilibrium_solution, time, step_size, dof,
                  nullptr /*GradientContainer*/, nullptr /*SolverHints**/, index,
                  b2linalg::Vector<T, b2linalg::Vdense>::null,
                  TypedElement<T>::d_value_d_dof_flags_null, matrix_unsym,
                  b2linalg::Vector<T, b2linalg::Vdense>::null);
 
            k_eff = matrix_unsym[0];
        }
 
        return k_eff;
    };
 
    template <typename MATRIX_FORMAT>
    T ComputeElementError(const solver::TypedSolver<T, MATRIX_FORMAT>& solver) {
        T result{};
        return result;
    };
 
    template <typename MATRIX_FORMAT>
    b2linalg::Vector<T, b2linalg::Vdense> AssembleElementEffectiveVector(
          const solver::TypedSolver<T, MATRIX_FORMAT>& solver, b2linalg::Index index) {
        // Define variables to use deprecated interface get_value().
        b2000::Model& model{solver.GetModel()};  // Make this constant (TB).
        const bool is_linear{solver.IsLinear()};
        const EquilibriumSolution equilibrium_solution{false};  // Über solver abfragen (TB).
        const double time{solver.GetTime()};
        const double step_size{solver.GetStepSize()};
        const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof{solver.GetDof()};
 
        b2linalg::Vector<T, b2linalg::Vdense> f_eff{};
 
        std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense>>
              dummy{};  // Größe richtig? (TB)
 
        this->get_value(
              model, is_linear, equilibrium_solution, time, step_size, dof,
              nullptr /*GradientContainer*/, nullptr /*SolverHints**/, index, f_eff,
              TypedElement<T>::d_value_d_dof_flags_null, dummy,
              b2linalg::Vector<T, b2linalg::Vdense>::null);
 
        return f_eff;
    };
 
    template <typename MATRIX_FORMAT>
    void ComputeElementGradient(
          const solver::TypedSolver<T, MATRIX_FORMAT>& solver,
          GradientContainer* const gradient_container) {
        b2000::Model& model{solver.GetModel()};  // Make this constant (TB).
        const bool is_linear{solver.IsLinear()};
        const EquilibriumSolution equilibrium_solution{false};  // Über solver abfragen (TB).
        const double time{solver.GetTime()};
        const double step_size{solver.GetStepSize()};
 
        const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof{solver.GetDof()};
 
        std::vector<b2linalg::Matrix<T, typename MATRIX_FORMAT::dense>> dummy{};
        b2linalg::Index index;
 
        this->get_value(
              model, is_linear, equilibrium_solution, time, step_size, dof, gradient_container,
              nullptr /*SolverHints**/, index, b2linalg::Vector<T, b2linalg::Vdense>::null,
              TypedElement<T>::d_value_d_dof_flags_null, dummy,
              b2linalg::Vector<T, b2linalg::Vdense>::null);
    };
 
    virtual void get_value(
          Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
          const double time, const double delta_time,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
          GradientContainer* gradient_container, SolverHints* solver_hints,
          b2linalg::Index& dof_numbering, b2linalg::Vector<T, b2linalg::Vdense>& value,
          const std::vector<bool>& d_value_d_dof_flags,
          std::vector<b2linalg::Matrix<T, b2linalg::Mrectangle>>& d_value_d_dof,
          b2linalg::Vector<T, b2linalg::Vdense>& d_value_d_time) {
        dof_numbering.clear();
        if (!value.is_null()) { value.resize(0); }
        for (size_t i = 0; i != d_value_d_dof.size(); ++i) {
            d_value_d_dof[i].resize(size_t(0), size_t(0));
        }
        if (!d_value_d_time.is_null()) { d_value_d_time.resize(0); }
    }
 
    virtual void get_value(
          Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
          const double time, const double delta_time,
          const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& dof,
          GradientContainer* gradient_container, SolverHints* solver_hints,
          b2linalg::Index& dof_numbering, b2linalg::Vector<T, b2linalg::Vdense>& value,
          const std::vector<bool>& d_value_d_dof_flags,
          std::vector<b2linalg::Matrix<T, b2linalg::Mpacked>>& d_value_d_dof,
          b2linalg::Vector<T, b2linalg::Vdense>& d_value_d_time) {
        dof_numbering.clear();
        if constexpr (std::is_same<T, b2000::csda<double>>::value) {
            b2linalg::Matrix<double, b2linalg::Mrectangle> r_dof;
            copy_m(dof, r_dof);
            b2linalg::Vector<double, b2linalg::Vdense> r_value;
            std::vector<b2linalg::Matrix<T, b2linalg::Mpacked>> r_d_value_d_dof;
            b2linalg::Vector<double, b2linalg::Vdense> r_d_value_d_time;
            if (!value.is_null()) { r_value.resize(0); }
            for (size_t i = 0; i != r_d_value_d_dof.size(); ++i) {
                r_d_value_d_dof[i].resize(size_t(0), size_t(0));
            }
            if (!r_d_value_d_time.is_null()) { r_d_value_d_time.resize(0); }
            copy_v(r_value, value);
            copy_vm(r_d_value_d_dof, d_value_d_dof);
            copy_v(r_d_value_d_time, d_value_d_time);
        } else {
            if (!value.is_null()) { value.resize(0); }
            for (size_t i = 0; i != d_value_d_dof.size(); ++i) {
                d_value_d_dof[i].resize(size_t(0), size_t(0));
            }
            if (!d_value_d_time.is_null()) { d_value_d_time.resize(0); }
        }
    }
 
    static type_t type;
    static std::vector<b2linalg::Matrix<T, b2linalg::Mrectangle>> d_value_d_dof_null;
};
 
template <typename T>
std::vector<b2linalg::Matrix<T>> TypedElement<T>::d_value_d_dof_null;
 
template <typename T>
typename TypedElement<T>::type_t TypedElement<T>::type(
      "TYPED_VALUE_ELEMENT", type_name<T>(), StringList(), element::module, &Element::type);
 
}  // namespace b2000
 
#endif  // B2_ELEMENT_H_