Skip to content
Snippets Groups Projects
physical_element.cpp 8.94 KiB
#include <iostream>
#include <cassert>

#include "physical_element.h"
#include "gmsh_io.h"
#include "sgr.hpp"

physical_element::physical_element()
{}

size_t
physical_element::original_position(void) const
{
    return m_original_position;
}

int
physical_element::dimension(void) const
{
    return m_dimension;
}

int
physical_element::orientation(void) const
{
    return m_orientation;
}

double
physical_element::measure(void) const
{
    return m_measure;
}

vec_quadpt_3d
physical_element::integration_points(void) const
{
    return m_phys_quadpoints;
}

int
physical_element::element_tag() const
{
    return m_element_tag;
}

int
physical_element::element_type() const
{
    return m_gmsh_elemtype;
}

vec_size_t
physical_element::node_tags(void) const
{
    return m_node_tags;
}

/* Return only one determinant: this is used for geometric order 1
 * where the determinants and the jacobians are constant */
double
physical_element::determinant(void) const
{
    assert(m_geometric_order == 1);
    return m_determinants[0];
}

double
physical_element::determinant(size_t iQp) const
{
    assert(iQp < m_determinants.size());
    return m_determinants[iQp];
}

/* Return all the determinants in all the integration points: this
 * is for elements with geometric order > 1 */
vec_double
physical_element::determinants(void) const
{
    assert(m_geometric_order >= 1);
    return m_determinants;
}

/* Return only one jacobian: this is used for geometric order 1
 * where the determinants and the jacobians are constant */
mat3d
physical_element::jacobian(void) const
{
    assert(m_geometric_order == 1);
    return m_jacobians[0];
}

/* Return the iQp-th jacobian: this is used for all
 * geometric orders 1 */
mat3d
physical_element::jacobian(size_t iQp) const
{
    assert(iQp < m_jacobians.size());
    return m_jacobians[iQp];
}

/* Return all the jacobians in all the integration points: this
 * is for elements with geometric order > 1 */
vec_mat3d
physical_element::jacobians(void) const
{
    assert(m_geometric_order >= 1);
    return m_jacobians;
}

point_3d
physical_element::barycenter(void) const
{
    return m_barycenter;
}

physical_elements_factory::physical_elements_factory(const entity_params& ep)
    : dim(ep.dim), tag(ep.tag), elemType(ep.etype), geom_order(ep.gorder),
      approx_order(ep.aorder)
{
    gmm::getElementsByType(elemType, cellTags, cellNodesTags, tag);
    gmm::getIntegrationPoints(elemType, quadrature_name(2*approx_order), cellIps, cellIws);
    gmm::getJacobians(elemType, cellIps, cellJacs, cellDets, cellPpts, tag);
    gmm::getBasisFunctionsOrientationForElements(elemType, basis_func_name(approx_order),
        cellOrientations, tag);
    assert(cellOrientations.size() == cellTags.size());
    gmm::getBarycenters(elemType, tag, false, false, barycenters);

    std::vector<double>     coord;

#ifdef USE_INITIAL_4_8_4_API
    gmm::getKeysForElements(elemType, basis_func_name(approx_order),
        keypairs, coord, tag, false);
#else
    gmm::getKeysForElements(elemType, basis_func_name(approx_order),
        tagKeys, entityKeys, coord, tag, false);
    assert(tagKeys.size() == entityKeys.size());
#endif

    keys_per_elem = gmm::getNumberOfKeysForElements(elemType, basis_func_name(approx_order));
#ifdef USE_INITIAL_4_8_4_API
    assert(keys_per_elem*cellTags.size() == keypairs.size());
#else
    assert(keys_per_elem*cellTags.size() == tagKeys.size());
#endif
}

std::vector<physical_element>
physical_elements_factory::get_elements()
{
    auto num_elems = cellTags.size();
    auto num_nodes = cellNodesTags.size()/cellTags.size();
    assert( cellTags.size()*num_nodes == cellNodesTags.size() );

    auto num_gf = cellIws.size();
    assert( num_gf * num_elems == cellDets.size() );

    std::vector<physical_element> ret;
    ret.reserve(num_elems);

    for (size_t elem = 0; elem < num_elems; elem++)
    {
        physical_element new_pe;
        new_pe.m_geometric_order = geom_order;
        new_pe.m_approximation_order = approx_order;
        new_pe.m_original_position = elem;
        new_pe.m_dimension = dim;
        new_pe.m_parent_entity_tag = tag;
        new_pe.m_orientation = cellOrientations[elem];
        new_pe.m_element_tag = cellTags[elem];
        new_pe.m_measure = 0.0;
        new_pe.m_gmsh_elemtype = elemType;
        new_pe.m_barycenter = point_3d(barycenters[3*elem+0],
                                       barycenters[3*elem+1],
                                       barycenters[3*elem+2]);

        new_pe.m_bf_keys.resize(keys_per_elem);
        for (size_t i = 0; i < keys_per_elem; i++)
        {
#ifdef USE_INITIAL_4_8_4_API
            auto [vi, vu] = keypairs[keys_per_elem*elem + i];
#else
            auto vi = tagKeys[keys_per_elem*elem + i];
            auto vu = entityKeys[keys_per_elem*elem + i];
#endif
            new_pe.m_bf_keys[i] = bf_key(vi,vu);
        }

        new_pe.m_node_tags.resize(num_nodes);
        for (size_t i = 0; i < num_nodes; i++)
            new_pe.m_node_tags[i] = cellNodesTags[num_nodes*elem+i];

        auto elem_base = elem*num_gf;

        for (size_t gf = 0; gf < num_gf; gf++)
        {
            auto gf_offset = elem_base + gf;

            const auto JSIZE = 3*3; /* Jacobian matrix size */
            auto jacs_base = gf_offset*JSIZE;
            assert(jacs_base+8 < cellJacs.size());

            /* GMSH returns jacobians by column */
            Eigen::Matrix3d jac;
            jac(0,0) = cellJacs[jacs_base+0];   /* dx/du */
            jac(1,0) = cellJacs[jacs_base+1];   /* dy/du */
            jac(2,0) = cellJacs[jacs_base+2];   /* dz/du */
            jac(0,1) = cellJacs[jacs_base+3];   /* dx/dv */
            jac(1,1) = cellJacs[jacs_base+4];   /* dy/dv */
            jac(2,1) = cellJacs[jacs_base+5];   /* dz/dv */
            jac(0,2) = cellJacs[jacs_base+6];   /* dx/dw */
            jac(1,2) = cellJacs[jacs_base+7];   /* dy/dw */
            jac(2,2) = cellJacs[jacs_base+8];   /* dz/dw */
            new_pe.m_jacobians.push_back(jac);

            auto dets_ofs = gf_offset;
            assert( dets_ofs < cellDets.size() );
            new_pe.m_determinants.push_back(cellDets[dets_ofs]);
            new_pe.m_measure += cellDets[dets_ofs]*cellIws[gf];

            const auto PSIZE = 3;
            auto pts_base = gf_offset*PSIZE;
            assert(pts_base + 2 < cellPpts.size());
            point_3d pt(cellPpts[pts_base+0], cellPpts[pts_base+1], cellPpts[pts_base+2]);
            quadrature_point_3d qpt(pt, cellIws[gf]*cellDets[dets_ofs]);
            new_pe.m_phys_quadpoints.push_back(qpt);
        }

        ret.push_back( std::move(new_pe) );
    }

    return ret;
}

element_key::element_key()
{}

element_key::element_key(const physical_element& pe)
{
    m_dim       = pe.dimension();
    m_elemType  = pe.element_type();

    auto nodes = pe.node_tags();

    if (m_dim == 2)
    {
        m_key_data[0] = 3;
        m_key_data[1] = nodes[0];
        m_key_data[2] = nodes[1];
        m_key_data[3] = nodes[2];
        std::sort(std::next(m_key_data.begin()), m_key_data.end());
        return;
    }

    std::stringstream ss;
    ss << "Key not implemented for element type " << m_elemType;
    throw std::invalid_argument(ss.str());
}

bool
element_key::operator<(const element_key& other) const
{
    if (m_dim != other.m_dim or m_elemType != other.m_elemType)
        return false;

    if (m_key_data[0] != other.m_key_data[0])
        return false;

    auto mb = std::next(m_key_data.begin());
    auto me = m_key_data.end();
    auto ob = std::next(other.m_key_data.begin());
    auto oe = other.m_key_data.end();

    return std::lexicographical_compare(mb, me, ob, oe);
}

bool
element_key::operator==(const element_key& other) const
{
    if (m_dim != other.m_dim or m_elemType != other.m_elemType)
        return false;

    if (m_key_data[0] != other.m_key_data[0])
        return false;

    auto mb = std::next(m_key_data.begin());
    auto me = m_key_data.end();
    auto ob = std::next(other.m_key_data.begin());

    return std::equal(mb, me, ob);
}


std::ostream& operator<<(std::ostream& os, const element_key& ek)
{
    os << "ek: " << ek.m_dim << " " << ek.m_elemType << " (";
    for (size_t i = 0; i < 8; i++)
        os << ek.m_key_data[i] << " ";
    os << ")";
    return os;
}



element_key_factory::element_key_factory(int dim, int tag, int etype)
{
    assert(dim == 2);

    std::vector<size_t> nTags;
    gmm::getElementFaceNodes(etype, 3, nTags, tag, true);

    for (size_t i = 0; i < nTags.size(); i+= 3)
    {
        element_key ek;
        ek.m_dim = dim;
        ek.m_elemType = etype;
        ek.m_key_data[0] = 3;
        ek.m_key_data[1] = nTags[i+0];
        ek.m_key_data[2] = nTags[i+1];
        ek.m_key_data[3] = nTags[i+2];

        std::sort(std::next(ek.m_key_data.begin()), ek.m_key_data.end());

        ekeys.push_back( ek );
    }
    std::sort(ekeys.begin(), ekeys.end());

    /* Not sure why there are duplicates */
    ekeys.erase( std::unique(ekeys.begin(), ekeys.end()), ekeys.end() );
}

std::vector<element_key>
element_key_factory::get_keys(void) const
{
    return ekeys;
}