Mechanics.H

Go to the documentation of this file.

//
// This is a general purpose integrator that focuses on solving 
// elasticity/mechanics equations in the absence of other multiphysics
// simulations.
// It is enabled by :code:`alamo.program=mechanics` if used on its own, in
// which case there is no prefix.
// If it is being used by another integrator, see that integrator to determine
// the value of :code:`[prefix]` (often equal to :code:`elastic`).
//
// This integrator inherits from :ref:`Integrator::Base::Mechanics`; see 
// documentation for that integrator for additional parameters.
//
// :bdg-primary-line:`Model setup`
// There are two basic tools for setting up a mechanics problem in this
// integrator.
//
// 1. The :code:`eta` field: this is used to mix models of different types.
//    Use :code:`nmodels` to specify how many material models to use, and then
//    specify each model as :code:`model1`, :code:`model2`, etc.
//    The type of moel is specified using the :code:`alamo.program.mechanics.model`
//    input.
// 
//    Once you have done this, you must determine the spatial distribution of each
//    model. This is done by choosing an IC for the eta field with :code:`ic.type`.
//    The :code:`IC::Expression` is the most general and recommended.
//    The models are then mixed linearly, i.e.
//
//    .. math::
//
//       W_{\textrm{eff}} = \sum_{i=1}^N W_i\,\eta_i(\mathbf{x})
// 
//    See the :ref:`Eshelby` test for an example of model mixing.
// 
// 2. The :code:`psi` field: this is used specifically for cases where a "void" region
//    is desired. Its usage is similar to the :code:`eta` case, and is conceptually
//    similar in that it scales the model field to near-zero in order to mimic the
//    (lack of) mechanical behavior in an empty region.
//    It is important to use :code:`psi` here, for reasons that are discussed in detail
//    in
//    `this paper <https://doi.org/10.1007/s00466-023-02325-8>`_.
//    The initialization of :code:`psi` is similar to that for :code:`eta`.
//    
//    See the :ref:`PlateHole` and :ref:`RubberPlateHole` for canonical exmaples.
//    The :ref:`Integrator::Fracture` and :ref:`Integrator::TopOp` integrators are examples
//    of integrators that leverage the psi property.
//
// :bdg-primary-line:`Body forces` 
// currently have limited support due to the relatively low number of 
// times they are needed. See the :ref:`Integrator::Base::Mechanics` documentation for 
// detail. 
// See the :ref:`TrigTest` test for examples of current body force implementation.
//
// :bdg-primary-line:`Boundary conditions` 
// are implemented using the :ref:`BC::Operator::Elastic` classes.
// See the documentation on these classes for more detail.
// See any of the mechanics-based tests for examples of boundary condition application.
//
 
#ifndef INTEGRATOR_MECHANICS_H
#define INTEGRATOR_MECHANICS_H
#include <iostream>
#include <fstream>
#include <iomanip>
#include <numeric>
 
#include "AMReX.H"
#include "AMReX_ParallelDescriptor.H"
#include "AMReX_ParmParse.H"
 
#include "IO/ParmParse.H"
#include "Integrator/Base/Mechanics.H"
 
 
#include "IC/IC.H"
#include "BC/BC.H"
#include "BC/Operator/Elastic/Constant.H"
#include "BC/Operator/Elastic/TensionTest.H"
#include "BC/Operator/Elastic/Expression.H"
 
#include "IC/Ellipse.H"
#include "IC/Voronoi.H"
#include "IC/Constant.H"
#include "IC/Expression.H"
#include "IC/BMP.H"
#include "IC/PNG.H"
#include "IC/PSRead.H"
#include "IC/Constant.H"
#include "BC/Constant.H"
#include "Numeric/Stencil.H"
 
#include "Model/Solid/Solid.H"
#include "Solver/Nonlocal/Linear.H"
#include "Solver/Nonlocal/Newton.H"
 
#include "Operator/Operator.H"
#include "Base/Mechanics.H"
 
namespace Integrator
{
template<class MODEL>
class Mechanics : virtual public Base::Mechanics<MODEL>
{
public:
    static constexpr const char *name = "mechanics";
 
    Mechanics() : Base::Mechanics<MODEL>() { }
    Mechanics(IO::ParmParse& pp) : Base::Mechanics<MODEL>()
    {
        Parse(*this, pp);
    }
 
    ~Mechanics()
    {
        delete ic_eta;
        delete ic_psi;
        delete ic_trac_normal;
        delete bc_psi;
        delete bc_trac_normal;
    }
 
    // Mechanics inputs. See also :ref:`Integrator::Base::Mechanics`
    static void Parse(Mechanics& value, IO::ParmParse& pp)
    {
        pp.queryclass<Base::Mechanics<MODEL>>(value);
 
        int nmodels;
        pp_query_default("nmodels", nmodels, 1); // Number of elastic model varieties
 
        pp.queryclass_enumerate<MODEL>("model",value.models,nmodels);
 
        Util::Assert(INFO, TEST(value.models.size() > 0));
        value.RegisterNodalFab(value.eta_mf, value.models.size(), 2, "eta", true);
        // Refinement threshold for eta field
        pp_query_default("eta_ref_threshold", value.m_eta_ref_threshold, 0.01);
        // Refinement threshold for strain gradient
        pp_query_default("ref_threshold", value.m_elastic_ref_threshold, 0.01);
        // Explicity impose neumann condition on model at domain boundaries (2d only)
        pp_query_default("model_neumann_boundary",value.model_neumann_boundary,false);
 
        // Read in IC for eta
        if (nmodels > 1 && pp.contains("ic.type"))
        {
            // Select the initial condition for eta
            pp.select<IC::Constant,IC::Ellipse,IC::Voronoi,IC::BMP,IC::PNG,IC::Expression,IC::PSRead>("ic",value.ic_eta,value.geom);
 
            value.eta_reset_on_regrid = true;
            // Whether to re-initialize eta when re-gridding occurs.
            // Default is false unless eta ic is set, then default is.
            // true.
            pp_query_default("eta.reset_on_regrid", value.eta_reset_on_regrid, true);
        }
 
        // Read in IC for psi
        if (pp.contains("psi.ic.type"))
        {
            // Select initial condition for psi field
            pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead,IC::PNG>("psi.ic",value.ic_psi,value.geom);
 
            value.bc_psi = new BC::Nothing();
            value.RegisterNewFab(value.psi_mf, value.bc_psi, 1, 2, "psi", value.plot_psi);
            value.psi_on = true;
 
            value.psi_reset_on_regrid = true;
            // Whether to re-initialize psi when re-gridding occurs.
            // Default is false unless a psi ic is set, then default is
            // true.
            pp_query_default("psi.reset_on_regrid", value.psi_reset_on_regrid, true);
        }
 
        if (pp.contains("trac_normal.ic.type"))
        {
            // Read in IC for the "normal traction" field (applied at diffuse boundaries)
            pp.select<IC::Ellipse,IC::Constant,IC::Expression,IC::PSRead>("trac_normal.ic",value.ic_trac_normal,value.geom); 
 
            if (value.ic_trac_normal)
            {
                value.bc_trac_normal = new BC::Nothing();
                value.RegisterNewFab(value.trac_normal_mf, value.bc_trac_normal, 1, 2, "trac_normal", true);
            }
        }
        Util::Message(INFO);
 
    }
 
    void Initialize(int lev) override
    {
        Base::Mechanics<MODEL>::Initialize(lev);
        eta_mf[lev]->setVal(0.0);
        if (models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf);
        else eta_mf[lev]->setVal(1.0);
 
        if (psi_on) ic_psi->Initialize(lev, psi_mf);
        if (ic_trac_normal) ic_trac_normal->Initialize(lev, trac_normal_mf);
    }
 
    virtual void UpdateModel(int a_step, Set::Scalar time) override
    {
        if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
 
        if (ic_trac_normal)
        {
            for (int lev = 0; lev <= finest_level; ++lev)
            {
                ic_trac_normal->Initialize(lev, trac_normal_mf, time);
                psi_mf[lev]->FillBoundary();
                amrex::Box domain = this->geom[lev].Domain();
                domain.convert(amrex::IntVect::TheNodeVector());
                Set::Vector DX(geom[lev].CellSize());
 
                for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
                {
                    amrex::Box bx = mfi.nodaltilebox();
                    bx = bx & domain;
 
                    amrex::Array4<Set::Vector> const& rhs = rhs_mf[lev]->array(mfi);
                    amrex::Array4<const Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
                    amrex::Array4<const Set::Scalar> const& trac_normal = trac_normal_mf[lev]->array(mfi);
 
                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                        Set::Vector grad = Numeric::CellGradientOnNode(psi, i, j, k, 0, DX.data());
                        rhs(i,j,k) = trac_normal(i,j,k) * grad;
                    });
                }
                Util::RealFillBoundary(*rhs_mf[lev], geom[lev]);
            }
        }
 
        if (a_step > 0) return;
 
        for (int lev = 0; lev <= finest_level; ++lev)
        {
            eta_mf[lev]->FillBoundary();
 
            amrex::Box domain = this->geom[lev].Domain();
            domain.convert(amrex::IntVect::TheNodeVector());
 
            Set::Vector DX(geom[lev].CellSize());
 
            for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
            {
                amrex::Box bx = mfi.grownnodaltilebox();
 
                amrex::Array4<MODEL> const& model = model_mf[lev]->array(mfi);
                amrex::Array4<const Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
 
                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
                    model(i, j, k) = MODEL::Zero();
                    for (unsigned int n = 0; n < models.size(); n++)
                        model(i, j, k) += eta(i, j, k, n) * models[n];
                });
            }
 
 
            if (model_neumann_boundary)
            {
                Util::AssertException(INFO,TEST(AMREX_SPACEDIM==2),"neumann boundary works in 2d only");
                for (MFIter mfi(*model_mf[lev], false); mfi.isValid(); ++mfi)
                {
                    amrex::Box bx = mfi.grownnodaltilebox() & domain;
                    amrex::Array4<MODEL> const& model = this->model_mf[lev]->array(mfi);
                    const Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
 
                    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
                    {
                        if      (i==lo.x && j==lo.y)
                            model(i,j,k) = 0.5*(model(i+1,j,k)+model(i,j+1,k));
                        else if (i==lo.x && j==hi.y)
                            model(i,j,k) = 0.5*(model(i+1,j,k)+model(i,j-1,k));
                        else if (i==hi.x && j==lo.y)
                            model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j+1,k));
                        else if (i==hi.x && j==hi.y)
                            model(i,j,k) = 0.5*(model(i-1,j,k)+model(i,j-1,k));
 
                        else if (i==lo.x)
                            model(i,j,k) = model(i+1,j,k);
                        else if (i==hi.x)
                            model(i,j,k) = model(i-1,j,k);
                        else if (j==lo.y)
                            model(i,j,k) = model(i,j+1,k);
                        else if (j==hi.y)
                            model(i,j,k) = model(i,j-1,k);
                    });
                }
            }
 
            Util::RealFillBoundary(*model_mf[lev], geom[lev]);
        }
    }
 
    void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar a_time, int a_ngrow) override
    {
        if (m_type == Base::Mechanics<MODEL>::Type::Disable) return;
        Base::Mechanics<MODEL>::TagCellsForRefinement(lev, a_tags, a_time, a_ngrow);
 
        Set::Vector DX(geom[lev].CellSize());
        Set::Scalar DXnorm = DX.lpNorm<2>();
        for (amrex::MFIter mfi(*eta_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
        {
            amrex::Box bx = mfi.nodaltilebox();
            amrex::Array4<char> const& tags = a_tags.array(mfi);
            amrex::Array4<Set::Scalar> const& eta = eta_mf[lev]->array(mfi);
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
            {
                auto sten = Numeric::GetStencil(i, j, k, bx);
                for (int n = 0; n < eta_mf[lev]->nComp(); n++)
                {
                    Set::Vector grad = Numeric::Gradient(eta, i, j, k, n, DX.data(), sten);
                    if (grad.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
                        tags(i, j, k) = amrex::TagBox::SET;
                }
            });
            if (psi_on)
            {
                amrex::Array4<Set::Scalar> const& psi = psi_mf[lev]->array(mfi);
                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
                {
                    auto sten = Numeric::GetStencil(i, j, k, bx);
                    {
                        Set::Vector gradpsi = Numeric::Gradient(psi, i, j, k, 0, DX.data(), sten);
                        if (gradpsi.lpNorm<2>() * DXnorm > m_eta_ref_threshold)
                            tags(i, j, k) = amrex::TagBox::SET;
                    }
                });
            }
        }
    }
 
    void Regrid(int lev, Set::Scalar time) override
    {
        if (eta_reset_on_regrid && models.size() > 1 && ic_eta) ic_eta->Initialize(lev, eta_mf, time);
        if (psi_reset_on_regrid) ic_psi->Initialize(lev, psi_mf, time);
    }
 
protected:
    Set::Field<Set::Scalar> eta_mf;
    Set::Scalar m_eta_ref_threshold = NAN;
    std::vector<MODEL> models;
    IC::IC<Set::Scalar>* ic_eta = nullptr;
    IC::IC<Set::Scalar>* ic_psi = nullptr;
    IC::IC<Set::Scalar>* ic_trac_normal = nullptr;
    BC::BC<Set::Scalar>* bc_psi = nullptr;
    BC::BC<Set::Scalar>* bc_trac_normal = nullptr;
    bool psi_reset_on_regrid = false;
    bool eta_reset_on_regrid = false;
 
    bool model_neumann_boundary = false;
 
    Set::Field<Set::Scalar> trac_normal_mf;
 
    using Base::Mechanics<MODEL>::m_type;
    using Base::Mechanics<MODEL>::finest_level;
    using Base::Mechanics<MODEL>::geom;
    using Base::Mechanics<MODEL>::model_mf;
    using Base::Mechanics<MODEL>::psi_mf;
    using Base::Mechanics<MODEL>::psi_on;
    using Base::Mechanics<MODEL>::rhs_mf;
};
 
 
 
 
 
 
 
 
 
 
} // namespace Integrator
#endif