b2eidw_mesh_deformation.H

//------------------------------------------------------------------------
// b2eidw_mesh_deformation.H --
//
// written by Mathias Doreille
//
// Copyright (c) 2011-2012 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 B2EIDW_MESH_DEFORMATION_H_
#define B2EIDW_MESH_DEFORMATION_H_
 
#include <algorithm>
#include <cmath>
#include <set>
#include <vector>
 
#include "utils/b2database.H"
#include "utils/b2finite_rotation.H"
#include "utils/b2linear_algebra.H"
#include "utils/b2tensor_calculus.H"
 
namespace b2000 {
 
class EIDW_Mesh_deformation {
public:
    EIDW_Mesh_deformation() : power_disp(0), power_rot(0), power_rot_dist(0), boundary_size(0) {}
 
    void set_parameter(
          int power_disp_, int power_rot_, int power_rot_dist_, double boundary_size_) {
        power_disp = power_disp_;
        power_rot = power_rot_;
        power_rot_dist = power_rot_dist_;
        boundary_size = boundary_size_;
    }
 
    template <int ndim>
    void compute_deformation(
          const b2linalg::Matrix<double>& node_coor, b2linalg::Matrix<double>& node_disp,
          const b2linalg::Index& imposed_node_id, const b2linalg::Vector<double>& node_area,
          const b2linalg::Matrix<double>& node_normal,
          const b2linalg::Matrix<double>& node_rotation) {
        size_t imposed_node_ptr = 0;
        for (size_t i = 0; i != node_disp.size2(); ++i) {
            if (i == imposed_node_id[imposed_node_ptr]) {
                ++imposed_node_ptr;
            } else {
                double min_dist = 1e50;
                double p = 0;
                double pr = 0;
                double disp[ndim];
                std::fill_n(disp, ndim, 0);
                double disp_r[ndim];
                std::fill_n(disp_r, ndim, 0);
                for (size_t j = 0; j != imposed_node_id.size(); ++j) {
                    const size_t jj = imposed_node_id[j];
                    double v[ndim];
                    double nv = 0;
                    for (int k = 0; k != ndim; ++k) {
                        const double tmp = node_coor(k, i) - node_coor(k, jj);
                        v[k] = tmp;
                        nv += tmp * tmp;
                    }
                    nv = std::sqrt(nv);
                    min_dist = std::min(min_dist, nv);
                    double vr[ndim];
                    std::fill_n(vr, ndim, 0);
                    {
                        const double* rotator = &node_rotation(0, j);
                        for (int lj = 0; lj != ndim; ++lj) {
                            for (int li = 0; li != ndim; ++li, ++rotator) {
                                vr[li] += *rotator * v[lj];
                            }
                        }
                    }
                    const double omega = node_area[j] * std::pow(nv, -power_disp);
                    const double omega_r = node_area[j] * std::pow(nv, -power_rot);
                    p += omega;
                    pr += omega_r;
                    for (int k = 0; k != ndim; ++k) {
                        disp[k] += omega * node_disp(k, jj);
                        disp_r[k] += omega_r * (vr[k] - v[k]);
                    }
                }
                p = 1.0 / p;
                pr = std::pow(std::max(1.0, min_dist - boundary_size + 1), -power_rot_dist) / pr;
                for (int k = 0; k != ndim; ++k) { node_disp(k, i) = disp[k] * p + disp_r[k] * pr; }
            }
        }
    }
 
private:
    int power_disp;
    int power_rot;
    int power_rot_dist;
    double boundary_size;
};
 
class EIDW_Mesh_deformationDB : public EIDW_Mesh_deformation {
public:
    EIDW_Mesh_deformationDB() : database(nullptr) {}
 
    void init(DataBase& db) { database = &db; }
 
    void compute_deformation(const int ndim) {
        b2linalg::Matrix<double> node_coor;
        database->get(SetName("COOR", 1), node_coor);
 
        b2linalg::Matrix<double> node_disp(ndim, node_coor.size2());
        b2linalg::Index imposed_node_id;
        {
            std::set<size_t> imposed_node_id_set;
            b2linalg::Matrix<double> ebc;
            database->get(SetName("EBC", 1, 0, 0, 1), ebc);
            for (size_t i = 0; i != ebc.size2(); ++i) {
                const size_t n = int(ebc(0, i)) - 1;
                imposed_node_id_set.insert(n);
                node_disp(int(ebc(1, i)) - 1, n) = ebc(2, i);
            }
            imposed_node_id.assign(imposed_node_id_set.begin(), imposed_node_id_set.end());
        }
 
        b2linalg::Vector<double> node_area(imposed_node_id.size());
        b2linalg::Matrix<double> node_normal(ndim, imposed_node_id.size());
        b2linalg::Matrix<double> node_rotation(ndim * ndim, imposed_node_id.size());
 
        switch (ndim) {
            case 2:
                compute_surface_node_2(
                      node_coor, node_disp, imposed_node_id, node_area, node_normal, node_rotation);
                EIDW_Mesh_deformation::compute_deformation<2>(
                      node_coor, node_disp, imposed_node_id, node_area, node_normal, node_rotation);
                break;
            case 3:
                compute_surface_node_3(
                      node_coor, node_disp, imposed_node_id, node_area, node_normal, node_rotation);
                EIDW_Mesh_deformation::compute_deformation<3>(
                      node_coor, node_disp, imposed_node_id, node_area, node_normal, node_rotation);
                break;
            default:
                Exception() << THROW;
                break;
        }
 
        database->set(SetName("DISP", 1, 0, 0, 1), node_disp);
    }
 
private:
    DataBase* database;
 
    void compute_surface_node_2(
          const b2linalg::Matrix<double>& node_coor, const b2linalg::Matrix<double>& node_disp,
          const b2linalg::Index& imposed_node_id, b2linalg::Vector<double>& node_area,
          b2linalg::Matrix<double>& node_normal, b2linalg::Matrix<double>& node_rotation) {
        ArrayRTable etab;
        database->get(SetName("ETAB", 1), etab);
        std::vector<int> nb_contrib(imposed_node_id.size());
        for (size_t j = 0; j != etab.size(); ++j) {
            const RTable& rtable = etab[j];
            if (rtable.get_int("ITYP", -1) != 266) { continue; }
            std::vector<int> nodes;
            rtable.get("NODES", nodes);
            if (nodes.size() != 2) { Exception() << THROW; }
            nodes[0] -= 1;
            nodes[1] -= 1;
            double v[2];
            double dv[2];
            double nv = 0;
            double ndv = 0;
            for (int i = 0; i != 2; ++i) {
                double tmp = node_coor(i, nodes[1]) - node_coor(i, nodes[0]);
                v[i] = tmp;
                nv += tmp * tmp;
                tmp += node_disp(i, nodes[1]) - node_disp(i, nodes[0]);
                dv[i] = tmp;
                ndv += tmp * tmp;
            }
            nv = std::sqrt(nv);
            ndv = std::sqrt(ndv);
            const double n[2] = {v[1] / nv, -v[0] / nv};
            const double omega = (v[0] * dv[1] - v[1] * dv[0]) / (nv * ndv);
            nv *= 0.5;
            for (int i = 0; i != 2; ++i) {
                const size_t node_id =
                      std::lower_bound(
                            imposed_node_id.begin(), imposed_node_id.end(), size_t(nodes[i]))
                      - imposed_node_id.begin();
                ++nb_contrib[node_id];
                node_area[node_id] += nv;
                for (int k = 0; k != 2; ++k) { node_normal(k, node_id) += n[k]; }
                node_rotation(0, node_id) += omega;
            }
        }
        for (size_t j = 0; j != imposed_node_id.size(); ++j) {
            const double omega = node_rotation(0, j) / nb_contrib[j];
            const double c = std::cos(omega);
            const double s = std::sin(omega);
            node_rotation(0, j) = c;
            node_rotation(1, j) = s;
            node_rotation(2, j) = -s;
            node_rotation(3, j) = c;
            normalise_2(&node_normal(0, j));
        }
    }
 
    void compute_surface_node_3(
          const b2linalg::Matrix<double>& node_coor, const b2linalg::Matrix<double>& node_disp,
          const b2linalg::Index& imposed_node_id, b2linalg::Vector<double>& node_area,
          b2linalg::Matrix<double>& node_normal, b2linalg::Matrix<double>& node_rotation) {
        ArrayRTable etab;
        database->get(SetName("ETAB", 1), etab);
        std::vector<int> nb_contrib(imposed_node_id.size());
        for (size_t j = 0; j != etab.size(); ++j) {
            const RTable& rtable = etab[j];
            switch (rtable.get_int("ITYP", -1)) {
                case 268:  // T3
                {
                    std::vector<int> nodes;
                    rtable.get("NODES", nodes);
                    if (nodes.size() != 3) { Exception() << THROW; }
                    for (int k = 0; k != 3; ++k) { nodes[k] -= 1; }
                    double v[2][3];
                    double dv[2][3];
                    for (int k = 0; k != 2; ++k) {
                        for (int i = 0; i != 3; ++i) {
                            const double tmp = node_coor(i, nodes[k + 1]) - node_coor(i, nodes[0]);
                            v[k][i] = tmp;
                            dv[k][i] = tmp + node_disp(i, nodes[k + 1]) - node_disp(i, nodes[0]);
                        }
                    }
                    double n[3];
                    outer_product_3(v[0], v[1], n);
                    const double s = normalise_3(n) / 3;
                    double dn[3];
                    outer_product_3(dv[0], dv[1], dn);
                    normalise_3(dn);
                    double rot[3];
                    outer_product_3(n, dn, rot);
                    for (int k = 0; k != 3; ++k) {
                        const size_t node_id = std::lower_bound(
                                                     imposed_node_id.begin(), imposed_node_id.end(),
                                                     size_t(nodes[k]))
                                               - imposed_node_id.begin();
                        ++nb_contrib[node_id];
                        node_area[node_id] += s;
                        for (int i = 0; i != 3; ++i) {
                            node_normal(i, node_id) += n[i];
                            node_rotation(i, node_id) += rot[i];
                        }
                    }
                    break;
                }
                case 267:  // Q4
                {
                    std::vector<int> nodesq;
                    rtable.get("NODES", nodesq);
                    if (nodesq.size() != 4) { Exception() << THROW; }
                    for (int k = 0; k != 4; ++k) { nodesq[k] -= 1; }
                    for (int kk = 0; kk != 4; ++kk) {
                        const int nodes[3] = {
                              nodesq[kk], nodesq[(kk + 1) % 4], nodesq[(kk + 2) % 4]};
                        double v[2][3];
                        double dv[2][3];
                        for (int k = 0; k != 2; ++k) {
                            for (int i = 0; i != 3; ++i) {
                                const double tmp =
                                      node_coor(i, nodes[k + 1]) - node_coor(i, nodes[0]);
                                v[k][i] = tmp;
                                dv[k][i] =
                                      tmp + node_disp(i, nodes[k + 1]) - node_disp(i, nodes[0]);
                            }
                        }
                        double n[3];
                        outer_product_3(v[0], v[1], n);
                        const double s = normalise_3(n) / 2;
                        double dn[3];
                        outer_product_3(dv[0], dv[1], dn);
                        normalise_3(dn);
                        double rot[3];
                        outer_product_3(n, dn, rot);
                        const size_t node_id =
                              std::lower_bound(
                                    imposed_node_id.begin(), imposed_node_id.end(), nodes[1])
                              - imposed_node_id.begin();
                        ++nb_contrib[node_id];
                        node_area[node_id] += s;
                        for (int i = 0; i != 3; ++i) {
                            node_normal(i, node_id) += n[i];
                            node_rotation(i, node_id) += rot[i];
                        }
                    }
                    break;
                }
            }
        }
 
        for (size_t j = 0; j != imposed_node_id.size(); ++j) {
            double omega[3];
            for (int i = 0; i != 3; ++i) { omega[i] = node_rotation(i, j) / nb_contrib[j]; }
            FiniteRotation<double>(omega).get_rotator(
                  reinterpret_cast<double(*)[3]>(&node_rotation(0, j)));
            normalise_3(&node_normal(0, j));
        }
    }
};
 
}  // namespace b2000
 
#endif /* B2EIDW_MESH_DEFORMATION_H_ */