body.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/hybridutilities.hh>
#include "body.hh"

#include <dune/tectonic/utils/debugutils.hh>

// -------------------------------
// --- LeafBody Implementation ---
// -------------------------------
template <class GridTEMPLATE, class VectorType>
LeafBody<GridTEMPLATE, VectorType>::LeafBody(
        const std::shared_ptr<BodyData<dim>>& bodyData,
        const std::shared_ptr<GridTEMPLATE>& hostGrid) :
    bodyData_(bodyData),
    hostGrid_(hostGrid),
    deformation_(std::make_shared<DeformationFunction>(hostGrid_->leafGridView())),
    grid_(std::make_shared<GridType>(*hostGrid_, *deformation_)),
    gridView_(grid_->leafGridView()),
    assembler_(std::make_shared<Assembler>(gridView_)),
    matrices_()
{

    boundaryConditions_.resize(0);

    externalForce_ = [&](double relativeTime, VectorType &ell) {
        // Problem formulation: right-hand side
        std::vector<std::shared_ptr<BoundaryCondition>> leafNeumannConditions;
        boundaryConditions("neumann", leafNeumannConditions);
        ell.resize(gravityFunctional_.size());
        ell = gravityFunctional_;

        for (size_t i=0; i<leafNeumannConditions.size(); i++) {
            const auto& leafNeumannCondition = leafNeumannConditions[i];

            VectorType b;
            assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(),
                                        *leafNeumannCondition->boundaryNodes(),
                                        b,
                                        *leafNeumannCondition->boundaryFunction(),
                                        relativeTime);
            ell += b;
        }
    };
}


template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::assemble() {

    // assemble matrices
    assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
    //print(*matrices_.elasticity, "elasticity");

    assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
    //print(*matrices_.damping, "viscosity");

    assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
    //print(*matrices_.mass, "mass");

    // assemble forces
    assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
    //print(gravityFunctional_, "gravity");
}

// setter and getter

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::data() const
-> const std::shared_ptr<BodyData<dim>>& {

    return bodyData_;
}


template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::setDeformation(const VectorType& def) {
    deformation_->setDeformation(def);
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::deformation() const
-> DeformationFunction& {

    return *deformation_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::grid() const
-> std::shared_ptr<GridType> {

    return grid_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::hostGrid() const
-> std::shared_ptr<GridTEMPLATE> {

    return hostGrid_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::gridView() const
-> const GridView& {

    return gridView_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::nVertices() const
-> size_t {
    return gridView_.size(dim);
}


template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::assembler() const
-> const std::shared_ptr<Assembler>& {

    return assembler_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::matrices() const
-> const Matrices<Matrix,1>& {

    return matrices_;
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::externalForce() const
-> const ExternalForce& {

    return externalForce_;
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
    boundaryConditions_.emplace_back(boundaryCondition);
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::boundaryConditions(
        const std::string& tag,
        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {

    selectedConditions.resize(0);
    for (size_t i=0; i<boundaryConditions_.size(); i++) {
        if (boundaryConditions_[i]->tag() == tag)
            selectedConditions.emplace_back(boundaryConditions_[i]);
    }
}

template <class GridTEMPLATE, class VectorType>
auto LeafBody<GridTEMPLATE, VectorType>::boundaryConditions() const
-> const std::vector<std::shared_ptr<BoundaryCondition>>& {

    return boundaryConditions_;
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::boundaryPatchNodes(
        const std::string& tag,
        BoundaryPatchNodes& nodes) const {

    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    nodes.resize(_boundaryConditions.size());

    for (size_t i=0; i<nodes.size(); i++) {
        nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
    }
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::boundaryNodes(
        const std::string& tag,
        BoundaryNodes& nodes) const {

    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    nodes.resize(_boundaryConditions.size());

    for (size_t i=0; i<nodes.size(); i++) {
        nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
    }
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::boundaryFunctions(
        const std::string& tag,
        BoundaryFunctions& functions) const {

    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    functions.resize(_boundaryConditions.size());

    for (size_t i=0; i<functions.size(); i++) {
        functions[i] = _boundaryConditions[i]->boundaryFunction().get();
    }
}

template <class GridTEMPLATE, class VectorType>
void LeafBody<GridTEMPLATE, VectorType>::boundaryPatches(
        const std::string& tag,
        BoundaryPatches& patches) const {

    std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    patches.resize(_boundaryConditions.size());

    for (size_t i=0; i<patches.size(); i++) {
        patches[i] = _boundaryConditions[i]->boundaryPatch().get();
    }
}


// ---------------------------
// --- Body Implementation ---
// ---------------------------

// setter and getter

template <class GridView>
auto Body<GridView>::data() const
-> const std::shared_ptr<BodyData<dim>>& {

    return bodyData_;
}


template <class GridView>
auto Body<GridView>::grid() const
-> std::shared_ptr<GridType> {

    return grid_;
}

template <class GridView>
auto Body<GridView>::gridView() const
-> const GridView& {

    return gridView_;
}

template <class GridView>
auto Body<GridView>::nVertices() const
-> size_t {
    return gridView_.size(dim);
}


template <class GridView>
void Body<GridView>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {

    boundaryConditions_.push_back(boundaryCondition);
}


template <class GridView>
void Body<GridView>::boundaryConditions(
        const std::string& tag,
        std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {

    selectedConditions.resize(0);
    for (size_t i=0; i<boundaryConditions_.size(); i++) {
        if (boundaryConditions_[i]->tag() == tag)
            selectedConditions.push_back(boundaryConditions_[i]);
    }
}


template <class GridView>
auto Body<GridView>::boundaryConditions() const
-> const std::vector<std::shared_ptr<BoundaryCondition>>& {

    return boundaryConditions_;
}


template <class GridView>
void Body<GridView>::boundaryPatchNodes(
        const std::string& tag,
        BoundaryPatchNodes& nodes) const {

    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    nodes.resize(_boundaryConditions.size());

    for (size_t i=0; i<nodes.size(); i++) {
        nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
    }
}

template <class GridView>
void Body<GridView>::boundaryNodes(
        const std::string& tag,
        BoundaryNodes& nodes) const {

    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    nodes.resize(_boundaryConditions.size());

    for (size_t i=0; i<nodes.size(); i++) {
        nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
    }
}

template <class GridView>
void Body<GridView>::boundaryFunctions(
        const std::string& tag,
        BoundaryFunctions& functions) const {

    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    functions.resize(_boundaryConditions.size());

    for (size_t i=0; i<functions.size(); i++) {
        functions[i] = _boundaryConditions[i]->boundaryFunction().get();
    }
}

template <class GridView>
void Body<GridView>::boundaryPatches(
        const std::string& tag,
        BoundaryPatches& patches) const {

    std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
    this->boundaryConditions(tag, _boundaryConditions);

    patches.resize(_boundaryConditions.size());

    for (size_t i=0; i<patches.size(); i++) {
        patches[i] = _boundaryConditions[i]->boundaryPatch().get();
    }
}

#include "body_tmpl.cc"