gmsh_io.cpp

#include <iostream>
#include <sstream>
#include <cassert>

#include "gmsh_io.h"
#include "entity.h"

std::string
quadrature_name(int order)
{
    std::stringstream ss;
    ss << "Gauss" << order;
    return ss.str();
}

std::string
basis_func_name(int order)
{
    (void) order;
    std::stringstream ss;
    //ss << "Lagrange";
    ss << "H1Legendre" << order;
    return ss.str();
}

std::string
basis_grad_name(int order)
{
    (void) order;
    std::stringstream ss;
    //ss << "GradLagrange";
    ss << "GradH1Legendre" << order;
    return ss.str();
}

model::model()
    : geometric_order(1), approximation_order(1)
{
    remesh(geometric_order);
}

model::model(int pgo, int pao)
    : geometric_order(pgo), approximation_order(pao)
{
    assert(geometric_order > 0);
    assert(approximation_order > 0);

    remesh(geometric_order);
}

model::model(const char *filename)
    : geometric_order(1), approximation_order(1)
{
    gmsh::open( std::string(filename) );
    remesh(geometric_order);
}

model::model(const char *filename, int pgo, int pao)
    : geometric_order(pgo), approximation_order(pao)
{
    assert(geometric_order > 0);
    assert(approximation_order > 0);

    gmsh::open( std::string(filename) );
    remesh(geometric_order);
}

model::~model()
{
    gmsh::clear();
}

void
model::remesh()
{
    gmm::generate( DIMENSION(3) );
    gmm::setOrder( geometric_order );
    entities.clear();
    import_gmsh_entities();
    map_boundaries();
}

void
model::remesh(int pgo)
{
    geometric_order = pgo;
    remesh();
}

std::vector<entity>::const_iterator
model::begin() const
{
    return entities.begin();
}

std::vector<entity>::const_iterator
model::end() const
{
    return entities.end();
}

std::vector<entity>::iterator
model::begin()
{
    return entities.begin();
}

std::vector<entity>::iterator
model::end()
{
    return entities.end();
}

entity&
model::entity_at(size_t pos)
{
    return entities.at(pos);
}

size_t
model::num_dofs(void) const
{
    size_t ret = 0;
    for (const auto& e : entities)
        ret += e.num_dofs();

    return ret;
}

size_t
model::num_fluxes(void) const
{
    size_t ret = 0;
    for (const auto& e : entities)
        ret += e.num_fluxes();

    return ret;
}

size_t
model::num_cells(void) const
{
    size_t ret = 0;
    for (const auto& e : entities)
        ret += e.num_cells();

    return ret;
}

void
model::update_connectivity(const entity& e, size_t entity_index)
{
    for (size_t iT = 0; iT < e.num_cells(); iT++)
    {
        for (size_t iF = 0; iF < e.num_faces_per_elem(); iF++)
        {
            size_t findex = e.num_faces_per_elem() * iT  + iF;
            neighbour_descriptor cd(entity_index, iT, iF);
            element_key fk( e.face(findex) );
            conn.connect( fk, cd );
        }
    }
}

void
model::import_gmsh_entities(void)
{
    gvp_t gmsh_entities;
    gm::getEntities(gmsh_entities, DIMENSION(3));

    size_t entity_index = 0;
    size_t dof_base = 0;
    size_t flux_base = 0;
    size_t index_base = 0;
    for (auto [dim, tag] : gmsh_entities)
    {
        std::vector<int> elemTypes;
        gmm::getElementTypes(elemTypes, dim, tag);
        
        assert(elemTypes.size() == 1);
        if (elemTypes.size() != 1)
            throw std::invalid_argument("Only one element type per entity is allowed");

        for (auto& elemType : elemTypes)
        {
            entity e({.dim = dim, .tag = tag, .etype = elemType,
                      .gorder = geometric_order, .aorder = approximation_order});
            
            e.sort_by_orientation();
            e.base(dof_base, flux_base, index_base);
            dof_base += e.num_dofs();
            flux_base += e.num_fluxes();
            index_base += e.num_cells();

            update_connectivity(e, entity_index);

            entities.push_back( std::move(e) );
            entity_index++;
        }
    }
}

void
model::map_boundaries(void)
{
    gvp_t entities;
    /* Create a table mapping the tag of a boundary to the list of
     * its face keys */
    std::map<size_t, std::vector<face_key>>   boundaries_elemTags;
    gm::getEntities(entities, DIMENSION(2));
    for (auto [dim, tag] : entities)
    {
        std::vector<int> elemTypes;
        gmm::getElementTypes(elemTypes, dim, tag);
        assert(elemTypes.size() == 1);
        for (auto& elemType : elemTypes)
        {
            element_key_factory ekf(DIMENSION(2), tag, elemType);
            boundary_map[tag] = ekf.get_keys();
        }
    }
}