b2element_klt_compute_shape_t.H

//------------------------------------------------------------------------
// b2element_klt_compute_shape_t.H --
//
//     Compute shape functions and derivatives for Bezier triangles.
//
// 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_COMPUTE_SHAPE_T_H_
#define _B2ELEMENT_KLT_COMPUTE_SHAPE_T_H_
 
#include <algorithm>
#include <cassert>
 
#include "elements/klt_shells/b2element_klt_generic_shape.H"
#include "elements/klt_shells/b2element_klt_util.H"
#include "elements/smr/b2element_continuum_integrate.H"
#include "elements/smr/b2element_continuum_shape.H"
#include "elements/smr/b2element_continuum_util.H"
#include "model/b2domain.H"
#include "model/b2element.H"
#include "model_imp/b2hnode.H"
#include "utils/b2tensor_calculus.H"
 
namespace b2000::klt {
 
template <int ORDER, typename DATA>
class TriangleBezier : public GenericShape {
public:
    static const int order = ORDER;
 
    static const size_t num_nodes = (ORDER + 2) * (ORDER + 1) / 2;
 
    static const size_t num_nodes_edge_c0 = ORDER + 1;
 
    static const size_t num_nodes_edge_c1 = ORDER * 2 + 1;
 
    static const size_t num_nodes_edge_c2 = ORDER * 3;
 
    TriangleBezier() : GenericShape((ORDER + 2) * (ORDER + 1) / 2, 3) {}
 
    IntegrationScheme create_ischeme() const {
        IS_T orig_ischeme;
        orig_ischeme.set_order((ORDER - 1) * 2);  // XXX change
        IntegrationScheme ischeme;
        for (int i = 0; i != orig_ischeme.get_num(); ++i) {
            double xi[2];
            double weight;
            orig_ischeme.get_point(i, xi, weight);
            ischeme.points.push_back(IntegrationScheme::Point(xi, weight));
        }
        return ischeme;
    }
 
    int get_edge_order(const int edge) const {
        assert(edge >= 0 && edge < 3);
        return 2 * ORDER;  // XXX change
    }
 
    int get_edge_shape_order(const int edge) const {
        assert(edge >= 0 && edge < 3);
        return ORDER;
    }
 
    void get_xi_all_from_xi_edge(const int edge, const double xi, double xi_all[2]) const {
        assert(edge >= 0 && edge < 3);
        const double mapping[3][2] = {{xi, 0.}, {1. - xi, xi}, {0, 1. - xi}};
        xi_all[0] = mapping[edge][0];
        xi_all[1] = mapping[edge][1];
    }
 
    IntegrationScheme create_ischeme_edge(const int edge) const {
        assert(edge >= 0 && edge < 3);
        IS_L orig_ischeme;
        orig_ischeme.set_order(get_edge_order(edge));
 
        IntegrationScheme ischeme;
        for (int i = 0; i != orig_ischeme.get_num(); ++i) {
            double xi;
            double weight;
            orig_ischeme.get_point(i, xi, weight);
            const double er = .5 * (1 + xi);
            const double mapping[3][2] = {{er, 0}, {1 - er, er}, {0, 1 - er}};
            ischeme.points.push_back(
                  IntegrationScheme::Point(mapping[edge][0], mapping[edge][1], .5 * weight));
        }
        return ischeme;
    }
 
    void get_edge_direction(const int edge, double edge_dir[2]) const {
        assert(edge >= 0 && edge < 3);
        switch (edge) {
            case 0:
                edge_dir[0] = 1;
                edge_dir[1] = 0;
                break;
            case 1:
                edge_dir[0] = -1;
                edge_dir[1] = 1;
                break;
            case 2:
                edge_dir[0] = 0;
                edge_dir[1] = -1;
                break;
        }
    }
 
    void get_node_indices(const int sub_type, const int sub_index, b2linalg::Index& indices) const {
        indices.clear();
        if (sub_type == SUB_ALL) {
            indices.reserve(num_nodes);
            for (size_t i = 0; i != num_nodes; ++i) { indices.push_back(i); }
        } else if (sub_type == SUB_EDGE_C0 || sub_type == SUB_EDGE_C1 || sub_type == SUB_EDGE_C2) {
            const int edge = sub_index;
            assert(edge >= 0 && edge < 3);
            const int continuity =
                  (sub_type == SUB_EDGE_C0 ? 0 : (sub_type == SUB_EDGE_C1 ? 1 : 2));
            const size_t n =
                  (continuity == 0 ? num_nodes_edge_c0
                                   : (continuity == 1 ? num_nodes_edge_c1 : num_nodes_edge_c2));
            indices.reserve(n);
            for (size_t i = 0; i != n; ++i) {
                const size_t ii = DATA::node_indices_edge_c2[edge][i];
                indices.push_back(ii);
            }
        } else if (
              sub_type == SUB_VERTEX_C0 || sub_type == SUB_VERTEX_C1 || sub_type == SUB_VERTEX_C2) {
            const int vertex = sub_index;
            assert(vertex >= 0 && vertex < 3);
            const int continuity =
                  (sub_type == SUB_VERTEX_C0 ? 0 : (sub_type == SUB_VERTEX_C1 ? 1 : 2));
            const size_t n = (continuity == 0 ? 1 : (continuity == 1 ? 3 : 6));
            indices.reserve(n);
            for (size_t i = 0; i != n; ++i) {
                const size_t ii = DATA::node_indices_vertex_c2[vertex][i];
                indices.push_back(ii);
            }
        }
    }
 
    void compute_nodes_interpolation_straight_d0(
          const double xi[2], b2linalg::Vector<double, b2linalg::Vdense>& N) const {
        N.resize(3);
        N[0] = 1 - xi[0] - xi[1];
        N[1] = xi[0];
        N[2] = xi[1];
    }
 
    void compute_nodes_interpolation_straight_d1(
          const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const {
        N.resize(3, 3);
 
        N[d00][0] = 1 - xi[0] - xi[1];
        N[d00][1] = xi[0];
        N[d00][2] = xi[1];
 
        N[d10][0] = -1;
        N[d10][1] = 1;
        N[d10][2] = 0;
 
        N[d01][0] = -1;
        N[d01][1] = 0;
        N[d01][2] = 1;
    }
 
    void compute_nodes_interpolation_straight_d2(
          const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const {
        N.resize(3, 6);
 
        N[d00][0] = 1 - xi[0] - xi[1];
        N[d00][1] = xi[0];
        N[d00][2] = xi[1];
 
        N[d10][0] = -1;
        N[d10][1] = 1;
        N[d10][2] = 0;
 
        N[d01][0] = -1;
        N[d01][1] = 0;
        N[d01][2] = 1;
 
        N[d20][0] = 0;
        N[d20][1] = 0;
        N[d20][2] = 0;
 
        N[d11][0] = 0;
        N[d11][1] = 0;
        N[d11][2] = 0;
 
        N[d02][0] = 0;
        N[d02][1] = 0;
        N[d02][2] = 0;
    }
};
 
template <int ORDER, typename DATA>
const size_t TriangleBezier<ORDER, DATA>::num_nodes_edge_c0;
 
template <int ORDER, typename DATA>
const size_t TriangleBezier<ORDER, DATA>::num_nodes_edge_c1;
 
template <int ORDER, typename DATA>
const size_t TriangleBezier<ORDER, DATA>::num_nodes_edge_c2;
 
}  // namespace b2000::klt
 
#endif  // B2ELEMENT_KLT_COMPUTE_SHAPE_T_H_