b2domain_database.H

Go to the documentation of this file.

//------------------------------------------------------------------------
// b2domain_database.H --
//
//
// written by Mathias Doreille
//            Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
//            Thomas Blome <thomas.blome@dlr.de>
//
// (c) 2004-2012,2016,2017 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// (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 B2DOMAIN_DATABASE_H_
#define B2DOMAIN_DATABASE_H_
 
#include <cassert>
#include <map>
#include <set>
#include <sstream>
#include <string>
#include <vector>
 
#include "b2branch_property.H"
#include "b2component.H"
#include "b2domain_property.H"
#include "b2fortran_element_property.H"
#include "b2material_property.H"
#include "b2ppconfig.h"
#include "ip/b2ip_transformations.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model/b2node.H"
#include "utils/b2allocator.H"
#include "utils/b2object.H"
#include "utils/b2threading.H"
#include "utils/b2type.H"
 
namespace b2000::b2dbv3 {
 
class Domain : public b2000::Domain {
public:
    enum struct element_type_t {
        real_symmetric,
        real_unsymmetric,
        csda_symmetric,
        csda_unsymmetric,
        complex_symmetric,
        complex_unsymmetric,
        unknown
    };
 
    void set_model(b2000::Model& model) override;
 
    void set_case(b2000::Case& case_) override;
 
    void set_subcase_id(const int subcase_id_) override {
        for (auto& i : list_branch) { i.property.set_subcase_id(subcase_id_); }
    }
 
    bool have_temperature() const override {
        for (auto& i : list_branch) {
            if (i.property.have_temperature()) { return true; }
        }
        return false;
    }
 
    size_t get_number_of_dof() const override { return nb_dof; }
 
    size_t get_number_of_elements() const override {
        size_t nb = 0;
        for (list_branch_t::const_iterator i = list_branch.begin(); i != list_branch.end(); ++i) {
            nb += i->get_number_of_elements();
        }
        return nb;
    }
 
    size_t get_number_of_elements_and_subelements() const override {
        size_t nb = 0;
        for (list_branch_t::const_iterator i = list_branch.begin(); i != list_branch.end(); ++i) {
            nb += i->get_number_of_elements_and_subelements();
        }
        return nb;
    }
 
    size_t get_number_of_nodes() const override {
        size_t nb = 0;
        for (list_branch_t::const_iterator i = list_branch.begin(); i != list_branch.end(); ++i) {
            nb += i->get_number_of_nodes();
        }
        return nb;
    }
 
    const std::vector<Element::VariableInfo> get_value_info() const override { return value_info; }
 
    struct Branch {
        struct CompFortranElementProperty {
            bool operator()(
                  const FortranElementProperty* a, const FortranElementProperty* b) const {
                return *a < *b;
            }
        };
 
        struct ParentNodeMapping {
            size_t parent_internal_branch;
            size_t parent_internal_element;
            double global_icoor[3];
        };
 
        struct ParentElementMapping {
            size_t parent_internal_branch;
            size_t parent_internal_element;
            b2linalg::Matrix<double> global_icoor;
        };
 
        using list_nodes_t = std::vector<Node*>;
        using list_elements_t = std::vector<Element*>;
        using ElementProperties = std::set<FortranElementProperty*, CompFortranElementProperty>;
 
        // To iterate over all elements and sub-elements of this
        // branch.  The sub-elements are determined and iterated
        // over only after the parent element has been called.
        class ElementIterator {
        public:
            ElementIterator(Branch& branch)
                : estart(branch.list_all_elements.begin()), eend(branch.list_all_elements.end()) {}
 
            Element* next() {
                if (first) {
                    first = false;
                } else {
                    increment();
                }
                Element* e = current();
                return e;
            }
 
        private:
            Element* current() {
                if (in_sub) { return *sstart; }
                if (estart != eend) { return *estart; }
                return nullptr;
            }
 
            void increment() {
                if (!in_sub) {
                    std::pair<size_t, Element**> s = (*estart)->get_subelements();
                    if (s.second != 0) {
                        sstart = s.second;
                        send = s.second + s.first;
                        in_sub = true;
                    }
                } else if (++sstart == send) {
                    in_sub = false;
                }
                if (!in_sub) {
                    ++estart;
                    sstart = send = nullptr;
                }
            }
 
            Branch::list_elements_t::iterator estart;
            Branch::list_elements_t::iterator eend;
 
            Element** sstart{};
            Element** send{};
 
            bool in_sub{};
            bool first{true};
        };
 
        Branch(DomainProperty* dp) : property(*dp) {}
 
        size_t get_number_of_elements() const { return list_all_elements.size(); }
 
        size_t get_number_of_elements_and_subelements() const {
            size_t nb = 0;
            for (size_t i = 0; i != list_all_elements.size(); ++i) {
                nb += 1 + list_all_elements[i]->get_subelements().first;
            }
            return nb;
        }
 
        size_t get_number_of_nodes() const { return list_nodes.size(); }
 
        list_elements_t list_all_elements;
 
        list_nodes_t list_nodes;
        size_t nb_node_dof{};
        size_t nb_element_dof{};
        std::string id{};
 
        Transformation3D<double> btog_transformation;
        Transformation3D<b2000::csda<double>> btog_transformation_cs;
 
        std::vector<Transformation3D<double>> ltob;
        std::vector<Transformation3D<b2000::csda<double>>> ltob_cs;
        std::vector<ParentNodeMapping> parent_node_mapping;
        std::vector<ParentElementMapping> parent_element_mapping;
 
        BranchProperty property;
 
        ElementProperties element_properties;
 
        int min_number_of_dof_by_element{std::numeric_limits<int>::max()};
        int max_number_of_dof_by_element{};
        int min_number_of_dof_by_node{std::numeric_limits<int>::max()};
        int max_number_of_dof_by_node{};
    };
 
    std::vector<size_t> first_element_internal_id_of_branch;
    std::vector<size_t> first_node_internal_id_of_branch;
 
    using list_branch_t = std::vector<Branch>;
    list_branch_t list_branch;
 
    using branch_eid_to_branch_iid_t = std::map<size_t, size_t>;
    branch_eid_to_branch_iid_t branch_eid_to_branch_iid;
 
    using Materials = std::map<int, std::pair<MaterialProperty*, FortranElementProperty::type_t*>>;
    Materials materials;
    Materials properties;
 
    MaterialProperty* get_material_property_from_mid(size_t mid) {
        Materials::iterator i = materials.find(mid);
        if (i == materials.end()) { return nullptr; }
        return i->second.first;
    }
 
    MaterialProperty* get_material_property_from_pid(size_t pid) {
        Materials::iterator i = properties.find(pid);
        if (i == properties.end()) { return nullptr; }
        return i->second.first;
    }
 
    b2linalg::Matrix<double, b2linalg::Mcompressed_col> ltob_transformation;
    b2linalg::Matrix<b2000::csda<double>, b2linalg::Mcompressed_col> ltob_transformation_cs;
    b2linalg::Matrix<std::complex<double>, b2linalg::Mcompressed_col> ltob_transformation_c;
 
    std::vector<const double*> nodes_local_referential;
    std::vector<const b2000::csda<double>*> nodes_local_referential_cs;
 
    ip::Transformations transformations;
 
    Branch& get_branch_from_eid(size_t eid) {
        branch_eid_to_branch_iid_t::iterator i = branch_eid_to_branch_iid.find(eid);
        if (i == branch_eid_to_branch_iid.end()) {
            KeyError() << "Internal branch " << eid << " not found." << THROW;
        }
        return list_branch[i->second];
    }
 
    Branch& get_branch_from_eid(const std::string& eid) {
        std::istringstream s(eid);
        size_t seid;
        s >> seid;
        return get_branch_from_eid(seid);
    }
 
    class NodeIterator : public b2000::Domain::NodeIterator {
    public:
        Node* next() override {
            while (bstart != bend) {
                if (start != end) { return *start++; }
                ++bstart;
                if (bstart != bend) {
                    start = bstart->list_nodes.begin();
                    end = bstart->list_nodes.end();
                }
            }
            return nullptr;
        }
 
    private:
        NodeIterator(Domain& domain)
            : bstart(domain.list_branch.begin()), bend(domain.list_branch.end()) {
            while (bstart != bend && bstart->list_nodes.empty()) { ++bstart; }
            if (bstart != bend) {
                start = bstart->list_nodes.begin();
                end = bstart->list_nodes.end();
            }
        }
        friend class b2dbv3::Domain;
        list_branch_t::iterator bstart;
        list_branch_t::const_iterator bend;
        Branch::list_nodes_t::iterator start;
        Branch::list_nodes_t::const_iterator end;
    };
 
    node_iterator get_node_iterator() override { return node_iterator(new NodeIterator(*this)); }
 
    // To iterate over all elements and sub-elements.  The
    // sub-elements are determined and iterated over only after
    // the parent element has been called.
    class ElementIterator : public b2000::Domain::ElementIterator {
    public:
        Element* next() override {
            if (first) {
                first = false;
            } else {
                increment();
            }
            Element* e = current();
            return e;
        }
 
    protected:
        ElementIterator(Domain& domain)
            : bstart(domain.list_branch.begin()), bend(domain.list_branch.end()) {
            assert(bstart != bend);
            estart = bstart->list_all_elements.begin();
            eend = bstart->list_all_elements.end();
            in_sub = false;
            sstart = send = nullptr;
        }
 
        Element* current() {
            if (in_sub) { return *sstart; }
            if (estart != eend) { return *estart; }
            return nullptr;
        }
 
        void increment() {
            if (estart != eend) {
                if (!in_sub) {
                    std::pair<size_t, Element**> s = (*estart)->get_subelements();
                    if (s.second != 0) {
                        sstart = s.second;
                        send = s.second + s.first;
                        in_sub = true;
                    }
                } else if (++sstart == send) {
                    in_sub = false;
                }
            }
            if (!in_sub) {
                if (++estart == eend) {
                    if (++bstart != bend) {
                        estart = bstart->list_all_elements.begin();
                        eend = bstart->list_all_elements.end();
                    }
                }
                sstart = send = nullptr;
            }
        }
 
        friend class b2dbv3::Domain;
        bool first{true};
        list_branch_t::iterator bstart;
        list_branch_t::const_iterator bend;
        Branch::list_elements_t::iterator estart;
        Branch::list_elements_t::iterator eend;
        bool in_sub;
        Element** sstart;
        Element** send;
    };
 
    element_iterator get_element_iterator() override {
        return element_iterator(new ElementIterator(*this));
    }
 
    ~Domain() override {
        // Don't delete any sub-elements --> do not use ElementIterator.
        for (std::vector<Branch>::iterator branch = list_branch.begin();
             branch != list_branch.end(); ++branch) {
            for (size_t i = 0; i != branch->list_all_elements.size(); ++i) {
                Element* e = branch->list_all_elements[i];
                e->decref_or_free();
            }
        }
        {
            node_iterator i = get_node_iterator();
            for (Node* n = i->next(); n != nullptr; n = i->next()) { n->decref_or_free(); }
        }
        {
            MaterialIterator* i = get_material_iterator();
            MaterialProperty* m;
            while ((m = i->next()) != nullptr) { m->decref_or_free(); }
            delete i;
        }
        {
            MaterialIterator* i = get_property_iterator();
            MaterialProperty* m;
            while ((m = i->next()) != nullptr) { m->decref_or_free(); }
            delete i;
        }
    }
 
    template <typename T>
    void save_global_dof(
          SetName name, const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
          RTable desc = RTable());
 
    template <typename T>
    void save_global_dof_indexed(
          SetName name, const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof,
          RTable desc = RTable());
 
    void save_field(
          const std::string& name,
          const b2linalg::Vector<double, b2linalg::Vdense_constref>& dof) override {
        RTable desc;
        desc.set("TYPE", "NODE");
        desc.set("SYSTEM", "BRANCH");
        save_global_dof(name, dof, desc);
    }
 
    void save_field(
          const std::string& name,
          const b2linalg::Vector<b2000::csda<double>, b2linalg::Vdense_constref>& dof) override {
        RTable desc;
        desc.set("TYPE", "NODE");
        desc.set("SYSTEM", "BRANCH");
        save_global_dof(name, dof, desc);
    }
 
    void save_field(
          const std::string& name,
          const b2linalg::Vector<std::complex<double>, b2linalg::Vdense_constref>& dof) override {
        RTable desc;
        desc.set("TYPE", "NODE");
        desc.set("SYSTEM", "BRANCH");
        save_global_dof(name, dof, desc);
    }
 
    template <typename T>
    void load_global_dof(
          SetName name, b2linalg::Vector<T, b2linalg::Vdense_ref> dof, RTable* desc = nullptr);
 
    void load_node_set(SetName name, std::vector<Node*>& nodes_set) {
        size_t s = 0;
        for (size_t branch = 0; branch != list_branch.size(); ++branch) {
            name.branch(list_branch[branch].id);
            if (database->has_key(name)) { s += database->dims(name)[0]; }
        }
        nodes_set.clear();
        nodes_set.reserve(s);
        for (size_t branch = 0; branch != list_branch.size(); ++branch) {
            name.branch(list_branch[branch].id);
            if (database->has_key(name)) {
                std::vector<mcInt32> list_node;
                database->get(name, list_node);
                for (size_t n = 0; n != list_node.size(); ++n) {
                    const size_t nn = list_node[n];
                    if (nn < 1 || nn > list_branch[branch].list_nodes.size()) {
                        DBError() << "invalid internal node id " << nn << " in the set " << name
                                  << "." << THROW;
                    }
                    nodes_set.push_back(list_branch[branch].list_nodes[nn - 1]);
                }
            }
        }
    }
 
    template <typename T>
    void load_nodal_values(SetName name_, b2linalg::Matrix<T>& value) {
        SetName name(name_);
        bool init = false;
        size_t internal_node_id = 0;
        for (size_t branch = 0; branch != list_branch.size(); ++branch) {
            const size_t internal_node_id_end =
                  internal_node_id + list_branch[branch].get_number_of_nodes();
            name.branch(list_branch[branch].id);
            if (database->has_key(name)) {
                if (!init) {
                    if (value.size1() == 0) {
                        const std::vector<size_t> dims = database->dims(name);
                        if (dims.empty()) { DBError() << THROW; }
                        if (dims.size() == 1) {
                            value.resize(1, get_number_of_nodes());
                        } else {
                            value.resize(database->dims(name)[1], get_number_of_nodes());
                        }
                    } else {
                        value.resize(value.size1(), get_number_of_nodes());
                    }
                    init = true;
                }
                database->get(
                      name, value(b2linalg::Interval(0, value.size1()),
                                  b2linalg::Interval(internal_node_id, internal_node_id_end)));
            } else if (value.size1() != 0) {
                value(b2linalg::Interval(0, value.size1()),
                      b2linalg::Interval(internal_node_id, internal_node_id_end))
                      .set_zero();
            }
 
            internal_node_id = internal_node_id_end;
        }
    }
 
    const b2linalg::Matrix<double, b2linalg::Mrectangle>& get_nodal_values_double(SetName name) {
        b2linalg::Matrix<double, b2linalg::Mrectangle>& res = list_nodal_values[name];
        if (res.size2() == 0) { load_nodal_values(name, res); }
        return res;
    }
 
    const b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle>& get_nodal_values_csda(
          SetName name) {
        b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle>& res =
              list_nodal_values_cs[name];
        if (res.size2() == 0) { load_nodal_values(name, res); }
        return res;
    }
 
    template <typename T>
    const b2linalg::Matrix<T, b2linalg::Mrectangle>& get_nodal_values(SetName name);
 
    void save_state(const std::string& state_id) override {
        if (state_buffers_size.first == 0 && state_buffers_size.second == 0) { return; }
        {
            std::vector<int> vi(state_buffers_size.first);
            std::vector<double> vd(state_buffers_size.second);
            std::pair<int*, double*> bufptr(&vi[0], &vd[0]);
            for (list_branch_t::iterator b = list_branch.begin(); b != list_branch.end(); ++b) {
                for (Branch::list_elements_t::iterator e = b->list_all_elements.begin();
                     e != b->list_all_elements.end(); ++e) {
                    bufptr = (*e)->get_state_buffer(bufptr);
                }
            }
            if (!vi.empty()) { database->set(state_id + ".I", vi); }
            if (!vd.empty()) { database->set(state_id + ".F", vd); }
        }
    }
 
    void load_state(const std::string& state_id) override {
        if (state_buffers_size.first == 0 && state_buffers_size.second == 0) { return; }
        {
            std::vector<int> vi(state_buffers_size.first);
            std::vector<double> vd(state_buffers_size.second);
            if (!vi.empty()) { database->get(state_id + ".I", vi); }
            if (!vd.empty()) { database->get(state_id + ".F", vd); }
            std::pair<const int*, const double*> bufptr(&vi[0], &vd[0]);
            for (list_branch_t::iterator b = list_branch.begin(); b != list_branch.end(); ++b) {
                for (Branch::list_elements_t::iterator e = b->list_all_elements.begin();
                     e != b->list_all_elements.end(); ++e) {
                    bufptr = (*e)->set_state_buffer(bufptr);
                }
            }
        }
    }
 
    b2linalg::SparseMatrixConnectivityType get_dof_connectivity_type() const override {
        return connectivity_type;
    }
 
    element_type_t get_element_base_type_code() const {
        switch (static_cast<element_type_t>(element_base_type_code)) {
            case element_type_t::real_symmetric:
                return element_type_t::real_symmetric;
            case element_type_t::real_unsymmetric:
                return element_type_t::real_unsymmetric;
            case element_type_t::csda_symmetric:
                return element_type_t::csda_symmetric;
            case element_type_t::csda_unsymmetric:
                return element_type_t::csda_unsymmetric;
            case element_type_t::complex_symmetric:
                return element_type_t::complex_symmetric;
            case element_type_t::complex_unsymmetric:
                return element_type_t::complex_unsymmetric;
            default:
                return element_type_t::unknown;
        }
    }
 
    bool is_unsymmetric() const {
        if (element_base_type_code == static_cast<int>(element_type_t::unknown)) {
            Exception() << THROW;
        }
        return element_base_type_code & static_cast<int>(element_type_t::real_unsymmetric);
    }
 
    bool is_csda() const {
        if (element_base_type_code == static_cast<int>(element_type_t::unknown)) {
            Exception() << THROW;
        }
        return element_base_type_code & static_cast<int>(element_type_t::csda_symmetric);
    }
 
    bool is_complex() const {
        if (element_base_type_code == static_cast<int>(element_type_t::unknown)) {
            Exception() << THROW;
        }
        return element_base_type_code & static_cast<int>(element_type_t::complex_symmetric);
    }
 
    Node* get_node(const std::string& node_name) override {
        ObjectNameToNodeMap::iterator i = object_name_to_node_map.find(node_name);
        if (i == object_name_to_node_map.end()) { return nullptr; }
        return i->second;
    }
 
    Element* get_element(const std::string& element_name) override {
        ObjectNameToElementMap::iterator i = object_name_to_element_map.find(element_name);
        if (i == object_name_to_element_map.end()) { return nullptr; }
        return i->second;
    }
 
    void get_adjacent_elements(
          std::pair<int, const Node* const*> node_list,
          std::vector<AdjacentElement>& adjacent_elements) override {
        init_node_connectivity();
 
        std::map<Element*, std::vector<int>> tmp;
        for (const Node* const* n = node_list.second; n != node_list.second + node_list.first;
             ++n) {
            size_t nid = (*n)->get_id();
            std::vector<std::pair<Element*, int>>::const_iterator begin =
                  node_element_connectivity.begin() + node_element_connectivity_index[nid];
            std::vector<std::pair<Element*, int>>::const_iterator end =
                  node_element_connectivity.begin() + node_element_connectivity_index[nid + 1];
            for (; begin != end; ++begin) { tmp[begin->first].push_back(begin->second); }
        }
        size_t nb_adjacent_elements = 0;
        for (std::map<Element*, std::vector<int>>::iterator i = tmp.begin(); i != tmp.end(); ++i) {
            if (i->second.size() == size_t(node_list.first)) { nb_adjacent_elements += 1; }
        }
        adjacent_elements.resize(nb_adjacent_elements);
        nb_adjacent_elements = 0;
        for (std::map<Element*, std::vector<int>>::iterator i = tmp.begin(); i != tmp.end(); ++i) {
            if (i->second.size() == size_t(node_list.first)) {
                AdjacentElement& ad = adjacent_elements[nb_adjacent_elements++];
                ad.element = i->first;
                ad.internal_node_id_list.swap(i->second);
            }
        }
    }
 
    void get_elements_of_same_type_near(
          const Element* e, double dist, std::vector<Element*>& element_list) override;
 
    void get_slave_node(std::vector<std::pair<const Node*, const Node*>>& slave_node) override;
 
    void set_dof(const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof) override {
        global_dof_array_double = dof;
        global_dof_array_csda_double =
              b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>::null;
        global_dof_array_complex_double =
              b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>::null;
        global_dof_array_double_tmp.resize(0, dof.size2());
        list_element_aabb_intersection.clear();
    }
 
    void set_dof(const b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>& dof)
          override {
        global_dof_array_double = b2linalg::Matrix<double, b2linalg::Mrectangle_constref>::null;
        global_dof_array_csda_double = dof;
        global_dof_array_complex_double =
              b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>::null;
        global_dof_array_double_tmp.resize(0, 0);
        list_element_aabb_intersection.clear();
    }
 
    void set_dof(const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& dof)
          override {
        global_dof_array_double = b2linalg::Matrix<double, b2linalg::Mrectangle_constref>::null;
        global_dof_array_csda_double =
              b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>::null;
        global_dof_array_complex_double = dof;
        global_dof_array_double_tmp.resize(0, 0);
        list_element_aabb_intersection.clear();
    }
 
    const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& get_dof(double) override {
        if (global_dof_array_double.is_null()) {
            if (global_dof_array_csda_double.is_null()
                && global_dof_array_complex_double.is_null()) {
                return global_dof_array_double;
            } else {
                if (global_dof_array_double_tmp.size1()
                    != global_dof_array_complex_double.size1()) {
                    global_dof_array_double_tmp.copy_real(global_dof_array_complex_double);
                    global_dof_array_double_tmpref = global_dof_array_double_tmp;
                }
                return global_dof_array_double_tmpref;
            }
        } else {
            return global_dof_array_double;
        }
    }
 
    const b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>& get_dof(
          b2000::csda<double>) override {
        if (global_dof_array_csda_double.is_null()) {
            if (global_dof_array_double.is_null() && global_dof_array_complex_double.is_null()) {
                return global_dof_array_csda_double;
            } else {
                TypeError(
                      "Cannot obtain csda dof because the solver is "
                      "real-valued or complex-valued.")
                      << THROW;
            }
        } else {
            return global_dof_array_csda_double;
        }
        return b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>::null;
    }
 
    const b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>& get_dof(
          std::complex<double>) override {
        if (global_dof_array_complex_double.is_null()) {
            if (global_dof_array_double.is_null() && global_dof_array_csda_double.is_null()) {
                return global_dof_array_complex_double;
            } else {
                TypeError(
                      "Cannot obtain complex dof because the solver is "
                      "real-valued or csda-valued.")
                      << THROW;
            }
        } else {
            return global_dof_array_complex_double;
        }
        return b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>::null;
    }
 
    const std::pair<std::map<std::string, size_t>, b2linalg::Index>& get_dof_field() override {
        if (dof_field.second.size() != nb_dof) {
            dof_field.first.clear();
            dof_field.first["Unknown"] = 0;
            dof_field.second.resize(nb_dof);
            dof_field.second.set_zero();
            DegreesOfFreedom::ListDofField ldf = {
                  dof_field.first,
                  dof_field.second,
            };
            {
                node_iterator i = get_node_iterator();
                for (Node* node = i->next(); node != nullptr; node = i->next()) {
                    node->add_dof_field(ldf);
                }
            }
            {
                element_iterator i = get_element_iterator();
                for (Element* element = i->next(); element != nullptr; element = i->next()) {
                    element->add_dof_field(ldf);
                }
            }
        }
        return dof_field;
    }
 
    std::pair<size_t, size_t> get_internal_branch_and_element_id(const Element& e) {
        const size_t eid = e.get_id();
        const std::vector<size_t>::const_iterator i =
              std::upper_bound(
                    first_element_internal_id_of_branch.begin(),
                    first_element_internal_id_of_branch.end(), eid)
              - 1;
        return std::pair<size_t, size_t>(i - first_element_internal_id_of_branch.begin(), eid - *i);
    }
 
    std::pair<size_t, size_t> get_internal_branch_and_node_id(const Node& n) {
        const size_t eid = n.get_id();
        const std::vector<size_t>::const_iterator i =
              std::upper_bound(
                    first_node_internal_id_of_branch.begin(),
                    first_node_internal_id_of_branch.end(), eid)
              - 1;
        return std::pair<size_t, size_t>(i - first_node_internal_id_of_branch.begin(), eid - *i);
    }
 
    Element* get_parent_element_and_nodes_mapping(
          b2000::Domain& parent_domain, const b2000::Element& element,
          b2linalg::Matrix<double>& parent_nodes_internal_coor) override {
        std::pair<size_t, size_t> iid = get_internal_branch_and_element_id(element);
        if (list_branch[iid.first].parent_element_mapping.empty()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " is missing on the database." << THROW;
        }
        const Branch::ParentElementMapping& pem =
              list_branch[iid.first].parent_element_mapping[iid.second];
        if (!parent_nodes_internal_coor.is_null()) {
            parent_nodes_internal_coor = pem.global_icoor;
        }
        Domain& d = dynamic_cast<Domain&>(parent_domain);
        if (pem.parent_internal_branch >= d.list_branch.size()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " has some wrong internal branch id." << THROW;
        }
        if (pem.parent_internal_element
            >= d.list_branch[pem.parent_internal_branch].list_all_elements.size()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " has some wrong internal element id." << THROW;
        }
        return d.list_branch[pem.parent_internal_branch]
              .list_all_elements[pem.parent_internal_element];
    }
 
    Element* get_parent_element_mapping(
          b2000::Domain& parent_domain, const b2000::Node& node,
          b2linalg::Vector<double>& parent_internal_coor) override {
        std::pair<size_t, size_t> iid = get_internal_branch_and_node_id(node);
        if (list_branch[iid.first].parent_node_mapping.empty()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " is missing on the database." << THROW;
        }
        const Branch::ParentNodeMapping& pnm =
              list_branch[iid.first].parent_node_mapping[iid.second];
        if (!parent_internal_coor.is_null()) {
            parent_internal_coor.assign(pnm.global_icoor, pnm.global_icoor + 3);
        }
        Domain& d = dynamic_cast<Domain&>(parent_domain);
        if (pnm.parent_internal_branch >= d.list_branch.size()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " has some wrong internal branch id." << THROW;
        }
        if (pnm.parent_internal_element
            >= d.list_branch[pnm.parent_internal_branch].list_all_elements.size()) {
            DBError() << "The set MULTISCALE-NODE-MAPPINGS." << list_branch[iid.first].id
                      << " has some wrong internal element id." << THROW;
        }
        return d.list_branch[pnm.parent_internal_branch]
              .list_all_elements[pnm.parent_internal_element];
    }
 
    const double* get_node_local_referential_double(const Node& node) const override {
        if (nodes_local_referential.empty()) { return nullptr; }
        assert(node.get_id() < nodes_local_referential.size());
        return nodes_local_referential[node.get_id()];
    }
 
    const b2000::csda<double>* get_node_local_referential_csda(const Node& node) const override {
        if (nodes_local_referential_cs.empty()) { return nullptr; }
        assert(node.get_id() < nodes_local_referential_cs.size());
        return nodes_local_referential_cs[node.get_id()];
    }
 
protected:
    class MaterialIterator {
    public:
        MaterialProperty* next() {
            if (start == end) { return nullptr; }
            return (*start++).second.first;
        }
 
    private:
        MaterialIterator(Materials& materials) : start(materials.begin()), end(materials.end()) {}
        Materials::iterator start;
        Materials::const_iterator end;
        friend class b2dbv3::Domain;
    };
 
    MaterialIterator* get_material_iterator() { return new MaterialIterator(materials); }
 
    MaterialIterator* get_property_iterator() { return new MaterialIterator(properties); }
 
private:
    Allocator allocator;
    Model* model{};
    DataBase* database{};
    size_t nb_dof{};
    b2linalg::SparseMatrixConnectivityType connectivity_type{
          b2linalg::sparse_matrix_connectivity_unknown};
    DomainProperty property{*this};
    std::vector<Element::VariableInfo> value_info;
    std::pair<size_t, size_t> state_buffers_size;
    int element_base_type_code{static_cast<int>(element_type_t::real_symmetric)};
 
    b2linalg::Matrix<double, b2linalg::Mrectangle_constref> global_dof_array_double;
    b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle_constref>
          global_dof_array_csda_double;
    b2linalg::Matrix<std::complex<double>, b2linalg::Mrectangle_constref>
          global_dof_array_complex_double;
 
    b2linalg::Matrix<double, b2linalg::Mrectangle> global_dof_array_double_tmp;
    b2linalg::Matrix<double, b2linalg::Mrectangle_constref> global_dof_array_double_tmpref;
 
    std::vector<std::pair<Element*, int>> node_element_connectivity;
    std::vector<size_t> node_element_connectivity_index;
 
    std::map<std::string, b2linalg::Matrix<double, b2linalg::Mrectangle>> list_nodal_values;
    std::map<std::string, b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle>>
          list_nodal_values_cs;
 
    std::pair<std::map<std::string, size_t>, b2linalg::Index> dof_field;
 
    using merged_node_type_t = std::map<std::pair<Node::type_t*, Node::type_t*>, Node::type_t*>;
    merged_node_type_t merged_node_type;
 
    using list_one_node_pear_node_type_t = std::map<Node::type_t*, Node*>;
 
    using ObjectNameToNodeMap = std::map<std::string, Node*>;
    ObjectNameToNodeMap object_name_to_node_map;
    using ObjectNameToElementMap = std::map<std::string, Element*>;
    ObjectNameToElementMap object_name_to_element_map;
 
    list_one_node_pear_node_type_t list_one_node_pear_node_type;
    Node::type_t* get_merged_node_type(Node::type_t* a, Node::type_t* b) {
        if (a == b) { return a; }
        if (a > b) { std::swap(a, b); }
        Node::type_t* mnt = merged_node_type[std::pair<Node::type_t*, Node::type_t*>(a, b)];
        if (mnt != nullptr) { return mnt; }
        if (list_one_node_pear_node_type.empty()) {
            const std::map<std::string, ObjectType*>& node_type_list =
                  node::module.get_object_types_list();
            for (std::map<std::string, ObjectType*>::const_iterator i = node_type_list.begin();
                 i != node_type_list.end(); ++i) {
                Node::type_t* nt = dynamic_cast<Node::type_t*>(i->second);
                if (nt == nullptr) {
                    Exception() << "Internal error: Object identified by \"" << i->first
                                << "\" cannot be cast to Node::type_t*." << THROW;
                }
                list_one_node_pear_node_type[nt] = nt->new_object(allocator);
            }
        }
        Node* na = list_one_node_pear_node_type[a];
        Node* nb = list_one_node_pear_node_type[b];
        for (list_one_node_pear_node_type_t::iterator i = list_one_node_pear_node_type.begin();
             i != list_one_node_pear_node_type.end(); ++i) {
            if (na->is_subnode(i->second) && nb->is_subnode(i->second)) {
                // we have a candidate node
                if (mnt == nullptr) {
                    mnt = i->first;
                } else {
                    // we take the node with the minimum dof
                    if (list_one_node_pear_node_type[mnt]->get_number_of_dof()
                        > i->second->get_number_of_dof()) {
                        mnt = i->first;
                    }
                }
            }
        }
        if (mnt == nullptr) {
            IOError() << "Could not find a common node type for the nodes " << a->get_cpp_name()
                      << " and " << b->get_cpp_name() << "." << THROW;
        } else {
            merged_node_type[std::pair<Node::type_t*, Node::type_t*>(a, b)] = mnt;
        }
        // std::cout << "merge " <<  a->get_cpp_name() << " and " << b->get_cpp_name() << " get " <<
        // mnt->get_cpp_name() << std::endl;
        return mnt;
    }
 
    void init_node_connectivity();
 
    struct ElementIndex {
        const std::type_info& type;
        bool operator<(const ElementIndex& e) const { return type.before(e.type); }
    };
 
    struct ElementIntersectionIndex {
        const double dist;
        const std::type_info& type;
        bool operator<(const ElementIntersectionIndex& e) const {
            if (dist == e.dist) { return type.before(e.type); }
            return dist < e.dist;
        }
    };
 
    struct AABBox {
        static const int ndim = 3;
        double min[3];
        double max[3];
    };
 
    template <typename SET = std::set<Element*>>
    struct PairIntersectionList {
        void init(
              size_t size, const AABBox* start_aabb_,
              const std::vector<std::pair<Element*, size_t>>& list_element_) {
            start_aabb = start_aabb_;
            list_element = &list_element_;
            pair_list.assign(size, std::set<Element*>());
            nb_pair = 0;
        }
        void add(const AABBox* a, const AABBox* b) {
            size_t ai = a - start_aabb;
            size_t bi = b - start_aabb;
            if ((*list_element)[ai].second == (*list_element)[bi].second) { return; }
            if (ai > bi) { std::swap(ai, bi); }
            if (pair_list[(*list_element)[ai].first->get_id()]
                      .insert((*list_element)[bi].first)
                      .second) {
                ++nb_pair;
            }
            return;
        }
        size_t get_nb_pair() const { return nb_pair; }
        std::vector<SET> pair_list;
        size_t nb_pair;
        const AABBox* start_aabb;
        const std::vector<std::pair<Element*, size_t>>* list_element;
    };
 
    std::map<ElementIntersectionIndex, PairIntersectionList<std::set<Element*>>>
          list_element_aabb_intersection;
    std::map<ElementIndex, std::vector<std::pair<Element*, size_t>>> list_element_abbbb;
 
    b2spin_rw_mutex list_element_aabb_intersection_mutex;
 
public:
    using type_t = ObjectTypeComplete<Domain, b2000::Domain::type_t>;
    static type_t type;
};
 
template <>
inline const b2linalg::Matrix<double, b2linalg::Mrectangle>& Domain::get_nodal_values(
      SetName name) {
    return get_nodal_values_double(name);
}
 
template <>
inline const b2linalg::Matrix<b2000::csda<double>, b2linalg::Mrectangle>& Domain::get_nodal_values(
      SetName name) {
    return get_nodal_values_csda(name);
}
 
}  // namespace b2000::b2dbv3
 
template <typename T>
void b2000::b2dbv3::Domain::save_global_dof(
      SetName name, const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, RTable desc) {
    desc.set("SYSTEM", "BRANCH");
    desc.set("TYPE", "CELL");
    size_t pos = 0;
 
    SetName name_e = name;
    name_e.gname(name.get_gname() + "_E");
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        name_e.branch(list_branch[branch].id);
        if (list_branch[branch].max_number_of_dof_by_element == 0) { continue; }
        if (list_branch[branch].min_number_of_dof_by_element
            < list_branch[branch].max_number_of_dof_by_element) {
            b2linalg::Matrix<T, b2linalg::Mrectangle> m(
                  list_branch[branch].max_number_of_dof_by_element,
                  list_branch[branch].list_all_elements.size());
            Branch::list_elements_t& le = list_branch[branch].list_all_elements;
            for (size_t i = 0; i != m.size2(); ++i) {
                int s = le[i]->get_number_of_dof();
                std::copy(&dof[pos], &dof[pos + s], &m(0, i));
                pos += s;
            }
            database->set(name_e, m);
        } else {
            b2linalg::Interval interval(pos, pos + list_branch[branch].nb_element_dof);
            database->set(
                  name_e, b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(
                                list_branch[branch].max_number_of_dof_by_element,
                                list_branch[branch].list_all_elements.size(), dof(interval)));
            pos += list_branch[branch].nb_element_dof;
        }
        database->set_desc(name_e, desc);
    }
 
    desc.set("TYPE", "NODE");
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        name.branch(list_branch[branch].id);
        if (list_branch[branch].max_number_of_dof_by_node == 0) { continue; }
        if (list_branch[branch].min_number_of_dof_by_node
            < list_branch[branch].max_number_of_dof_by_node) {
            b2linalg::Matrix<T, b2linalg::Mrectangle> m(
                  list_branch[branch].max_number_of_dof_by_node,
                  list_branch[branch].list_nodes.size());
            Branch::list_nodes_t& ln = list_branch[branch].list_nodes;
            for (size_t i = 0; i != m.size2(); ++i) {
                int s = ln[i]->get_number_of_dof();
                std::copy(&dof[pos], &dof[pos + s], &m(0, i));
                pos += s;
            }
            database->set(name, m);
        } else {
            b2linalg::Interval interval(pos, pos + list_branch[branch].nb_node_dof);
            database->set(
                  name, b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(
                              list_branch[branch].max_number_of_dof_by_node,
                              list_branch[branch].list_nodes.size(), dof(interval)));
            pos += list_branch[branch].nb_node_dof;
        }
        database->set_desc(name, desc);
    }
}
 
template <typename T>
void b2000::b2dbv3::Domain::save_global_dof_indexed(
      SetName name, const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, RTable desc) {
    desc.set("SYSTEM", "BRANCH");
    desc.set("TYPE", "NODE");
    desc.set("INDEXED", true);
 
    SetName name_e = name;
    name_e.gname(name.get_gname() + "_E");
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        if (list_branch[branch].max_number_of_dof_by_element == 0) { continue; }
        Branch::list_elements_t& le = list_branch[branch].list_all_elements;
        std::vector<T> res;
        res.reserve(list_branch[branch].max_number_of_dof_by_element * le.size() * 3);
        for (size_t i = 0; i != le.size(); ++i) {
            size_t dof_id = 1;
            for (std::pair<size_t, size_t> dn = le[i]->get_global_dof_numbering();
                 dn.first != dn.second; ++dn.first, ++dof_id) {
                if (dof[dn.first] != T(0)) {
                    res.push_back(le[i]->get_id() + 1);
                    res.push_back(dof_id);
                    res.push_back(dof[dn.first]);
                }
            }
        }
        if (!res.empty()) {
            name_e.branch(list_branch[branch].id);
            database->set(
                  name_e,
                  b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(3, res.size() / 3, &res[0]));
            database->set_desc(name_e, desc);
        }
    }
 
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        name.branch(list_branch[branch].id);
        if (list_branch[branch].max_number_of_dof_by_node == 0) { continue; }
        Branch::list_nodes_t& ln = list_branch[branch].list_nodes;
        std::vector<T> res;
        res.reserve(list_branch[branch].max_number_of_dof_by_node * ln.size() * 3);
        for (size_t i = 0; i != ln.size(); ++i) {
            size_t dof_id = 1;
            for (std::pair<size_t, size_t> dn = ln[i]->get_global_dof_numbering();
                 dn.first != dn.second; ++dn.first, ++dof_id) {
                if (dof[dn.first] != T(0)) {
                    res.push_back(ln[i]->get_id() + 1);
                    res.push_back(dof_id);
                    res.push_back(dof[dn.first]);
                }
            }
        }
        if (!res.empty()) {
            name.branch(list_branch[branch].id);
            database->set(
                  name,
                  b2linalg::Matrix<T, b2linalg::Mrectangle_constref>(3, res.size() / 3, &res[0]));
            database->set_desc(name, desc);
        }
    }
}
 
template <typename T>
void b2000::b2dbv3::Domain::load_global_dof(
      SetName name, b2linalg::Vector<T, b2linalg::Vdense_ref> dof, RTable* rtable) {
    size_t pos = 0;
    bool zero_init = true;
 
    SetName name_e = name;
    name_e.gname(name.get_gname() + "_E");
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        if (list_branch[branch].max_number_of_dof_by_element == 0) { continue; }
        name_e.branch(list_branch[branch].id);
        if (!database->has_key(name_e)) {
            std::fill_n(&dof[pos], list_branch[branch].nb_element_dof, 0);
            pos += list_branch[branch].nb_element_dof;
        } else {
            if (rtable) { database->get_desc(name_e, *rtable); }
            if (list_branch[branch].min_number_of_dof_by_element
                < list_branch[branch].max_number_of_dof_by_element) {
                zero_init = false;
                b2linalg::Matrix<T, b2linalg::Mrectangle> m(
                      list_branch[branch].max_number_of_dof_by_element,
                      list_branch[branch].list_all_elements.size());
                database->get(name_e, m);
                Branch::list_elements_t& ln = list_branch[branch].list_all_elements;
                for (size_t i = 0; i != m.size2(); ++i) {
                    int s = ln[i]->get_number_of_dof();
                    std::copy(&m(0, i), &m(0, i) + s, &dof[pos]);
                    pos += s;
                }
            } else {
                zero_init = false;
                b2linalg::Interval interval(pos, pos + list_branch[branch].nb_element_dof);
                database->get(
                      (const std::string&)name_e,
                      b2linalg::Matrix<T, b2linalg::Mrectangle_increment_ref>(
                            b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(
                                  list_branch[branch].max_number_of_dof_by_element,
                                  list_branch[branch].list_all_elements.size(), dof(interval))));
                pos += list_branch[branch].nb_element_dof;
            }
        }
    }
 
    for (size_t branch = 0; branch != list_branch.size(); ++branch) {
        name.branch(list_branch[branch].id);
        if (!database->has_key(name)) {
            std::fill_n(&dof[pos], list_branch[branch].nb_node_dof, 0);
            pos += list_branch[branch].nb_node_dof;
            continue;
        }
        RTable desc_;
        RTable desc = rtable ? *rtable : desc_;
        bool has_desc = database->get_desc(name, desc, false);
        if (has_desc && desc.get_string("DOMAIN", "") == "DOF" && desc.get_bool("INDEXED", false)) {
            zero_init = false;
            b2linalg::Matrix<T, b2linalg::Mrectangle> m;
            database->get(name, m);
            std::fill_n(&dof[pos], list_branch[branch].nb_node_dof, 0);
            for (size_t i = 0; i != m.size2(); ++i) {
                size_t n = size_t(b2000::real(m(0, i))) - 1;
                const std::pair<size_t, size_t> node_dof =
                      list_branch[branch].list_nodes[n]->get_global_dof_numbering();
                ssize_t d = ssize_t(b2000::real(m(1, i))) - 1;
                if (d < 0) {
                    DBError() << "A dof number of a node cannot be a negative value."
                                 " In the dataset "
                              << name << ", row " << i + 1 << "." << THROW;
                }
                d += node_dof.first;
                if (d <= ssize_t(node_dof.second)) { dof[d] = m(2, i); }
            }
            pos += list_branch[branch].nb_node_dof;
        } else if (has_desc && desc.get_string("TYPE", "") == "NODE-INDEXED") {
            zero_init = false;
            b2linalg::Matrix<T, b2linalg::Mrectangle> m;
            database->get(name, m);
            std::fill_n(&dof[pos], list_branch[branch].nb_node_dof, 0);
            for (size_t i = 0; i != m.size2(); ++i) {
                size_t n = size_t(b2000::real(m(0, i))) - 1;
                const std::pair<size_t, size_t> node_dof =
                      list_branch[branch].list_nodes[n]->get_global_dof_numbering();
                const size_t s = std::min(node_dof.second - node_dof.first, m.size1());
                std::copy(&m(1, i), &m(1, i) + s, &dof[pos + node_dof.first]);
            }
            pos += list_branch[branch].nb_node_dof;
        } else {
            zero_init = false;
            if (list_branch[branch].min_number_of_dof_by_node
                < list_branch[branch].max_number_of_dof_by_node) {
                b2linalg::Matrix<T, b2linalg::Mrectangle> m(
                      list_branch[branch].max_number_of_dof_by_node,
                      list_branch[branch].list_nodes.size());
                database->get(name, m);
                Branch::list_nodes_t& ln = list_branch[branch].list_nodes;
                for (size_t i = 0; i != m.size2(); ++i) {
                    int s = ln[i]->get_number_of_dof();
                    std::copy(&m(0, i), &m(0, i) + s, &dof[pos]);
                    pos += s;
                }
            } else {
                b2linalg::Interval interval(pos, pos + list_branch[branch].nb_node_dof);
                database->get(
                      (const std::string&)name,
                      b2linalg::Matrix<T, b2linalg::Mrectangle_increment_ref>(
                            b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(
                                  list_branch[branch].max_number_of_dof_by_node,
                                  list_branch[branch].list_nodes.size(), dof(interval))));
                pos += list_branch[branch].nb_node_dof;
            }
        }
    }
    if (zero_init) {
        DBError() << "The set " << name << " could not be read or is empty." << THROW;
    }
}
 
#endif  // B2DOMAIN_DATABASE_H_