generic_field.impl.hh

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 *
 * 
 * @file    generic_field.impl.hh
 * @brief   
 */
 
#ifndef GENERIC_FIELD_IMPL_HH_
#define GENERIC_FIELD_IMPL_HH_
 
#include <memory>
#include "mesh/mesh.h"
#include "mesh/partitioning.hh"
#include "mesh/accessors.hh"
 
#include "fields/generic_field.hh"
#include "fields/field_fe.hh"
#include "la/vector_mpi.hh"
#include "fields/fe_value_handler.hh"
#include "fields/field_constant.hh"
 
#include "fem/mapping_p1.hh"
#include "fem/fe_p.hh"
#include "fem/dofhandler.hh"
#include "fem/discrete_space.hh"
 
 
template <int spacedim>
auto GenericField<spacedim>::region_id(Mesh &mesh) -> IndexField {
    IndexField region_id;
    region_id.name("region_id");
    region_id.units( UnitSI::dimensionless() );
    region_id.set_mesh(mesh);
 
    RegionSet all_regions=mesh.region_db().get_region_set("ALL");
    for(Region reg : all_regions) {
        auto field_algo=std::make_shared<FieldConstant<spacedim, DoubleScalar>>();
        field_algo->set_value(reg.id());
        region_id.set(
                field_algo,
                0.0,              // time=0.0
                { reg.label() }); // set to region by its name
    }
    return region_id;
}
 
template <int spacedim>
auto GenericField<spacedim>::subdomain(Mesh &mesh) -> IndexField {
    MixedPtr<FE_P_disc> fe(0);
    std::shared_ptr<DOFHandlerMultiDim> dh = std::make_shared<DOFHandlerMultiDim>(mesh);
    std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &mesh, fe);
    dh->distribute_dofs(ds);
 
    auto field_subdomain_data = mesh.get_part()->subdomain_id_field_data();
    unsigned int data_size = field_subdomain_data->size();
    VectorMPI data_vec(data_size);
    ASSERT_EQ(dh->max_elem_dofs(), 1);
    unsigned int i_ele=0;
    for (auto cell : dh->own_range()) {
        data_vec.set( cell.get_loc_dof_indices()(0), (*field_subdomain_data)[i_ele] );
        ++i_ele;
    }
    std::shared_ptr< FieldFE<spacedim, DoubleScalar> > field_ptr = std::make_shared< FieldFE<spacedim, DoubleScalar> >();
    field_ptr->set_fe_data(dh, data_vec);
 
    IndexField subdomain;
    subdomain.name("subdomain");
    subdomain.units( UnitSI::dimensionless() );
    subdomain.set_mesh(mesh);
 
    subdomain.set(
        field_ptr,
        0.0,        // time=0.0
        { "ALL" }); // ALL regions
 
    return subdomain;
}
 
 
template <int spacedim>
auto GenericField<spacedim>::element_measure(Mesh &mesh) -> IndexField {
    MixedPtr<FE_P_disc> fe(0);
    std::shared_ptr<DOFHandlerMultiDim> dh = std::make_shared<DOFHandlerMultiDim>(mesh);
    std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &mesh, fe);
    dh->distribute_dofs(ds);
 
    VectorMPI data_vec(dh->lsize());
    ASSERT_EQ(dh->max_elem_dofs(), 1);
    unsigned int i_ele=0;
    for (auto cell : dh->own_range()) {
        data_vec.set( cell.get_loc_dof_indices()(0), cell.elm().measure() );
        ++i_ele;
    }
    std::shared_ptr< FieldFE<spacedim, DoubleScalar> > field_ptr = std::make_shared< FieldFE<spacedim, DoubleScalar> >();
    field_ptr->set_fe_data(dh, data_vec);
 
    IndexField measure;
    measure.name("element_measure");
    measure.units( UnitSI().md() );
    measure.set_mesh(mesh);
 
    measure.set(
        field_ptr,
        0.0,        // time=0.0
        { "ALL" }); // ALL regions
 
    return measure;
}
 
 
template <int spacedim>
auto GenericField<spacedim>::element_diameter(Mesh &mesh) -> IndexField {
    MixedPtr<FE_P_disc> fe(0);
    std::shared_ptr<DOFHandlerMultiDim> dh = std::make_shared<DOFHandlerMultiDim>(mesh);
    std::shared_ptr<DiscreteSpace> ds = std::make_shared<EqualOrderDiscreteSpace>( &mesh, fe);
    dh->distribute_dofs(ds);
 
    VectorMPI data_vec(dh->lsize());
    ASSERT_EQ(dh->max_elem_dofs(), 1);
    unsigned int i_ele=0;
    for (auto cell : dh->own_range()) {
        double h = 0;
        for (auto side : cell.side_range())
            h = max(h, side.diameter());
        data_vec.set( cell.get_loc_dof_indices()(0), h );
        ++i_ele;
    }
    std::shared_ptr< FieldFE<spacedim, DoubleScalar> > field_ptr = std::make_shared< FieldFE<spacedim, DoubleScalar> >();
    field_ptr->set_fe_data(dh, data_vec);
 
    IndexField diam;
    diam.name("element_diameter");
    diam.units( UnitSI().m() );
    diam.set_mesh(mesh);
 
    diam.set(
        field_ptr,
        0.0,        // time=0.0
        { "ALL" }); // ALL regions
 
    return diam;
}
 
 
#endif /* GENERAL_FIELD_IMPL_HH_ */