LCOV - code coverage report
Current view: top level - src/Integrator/Base - Mechanics.H (source / functions) Hit Total Coverage
Test: coverage_merged.info Lines: 158 226 69.9 %
Date: 2025-02-18 04:54:05 Functions: 87 192 45.3 %

          Line data    Source code
       1             : #ifndef INTEGRATOR_BASE_MECHANICS_H
       2             : #define INTEGRATOR_BASE_MECHANICS_H
       3             : 
       4             : #include "AMReX.H"
       5             : #include "BC/Operator/Elastic/Elastic.H"
       6             : #include "BC/Operator/Elastic/Constant.H"
       7             : #include "BC/Operator/Elastic/TensionTest.H"
       8             : #include "BC/Operator/Elastic/Expression.H"
       9             : #include "Integrator/Integrator.H"
      10             : #include "Numeric/Stencil.H"
      11             : #include "Model/Solid/Solid.H"
      12             : #include "Solver/Nonlocal/Linear.H"
      13             : #include "Solver/Nonlocal/Newton.H"
      14             : #include "Operator/Operator.H"
      15             : #include "IC/Constant.H"
      16             : #include "IC/Expression.H"
      17             : #include "IC/Trig.H"
      18             : 
      19             : namespace Integrator
      20             : {
      21             : namespace Base
      22             : {
      23             : template<class MODEL>
      24             : class Mechanics: virtual public Integrator
      25             : {
      26             : public:
      27             : 
      28             :     enum Type { Static, Dynamic, Disable };
      29             : 
      30         131 :     Mechanics() : Integrator()
      31          25 :     {}
      32             : 
      33          25 :     ~Mechanics()
      34             :     {
      35           4 :         delete ic_rhs;
      36          25 :         delete velocity_ic;
      37          25 :         delete bc;
      38          50 :     }
      39             : 
      40             :     // The mechanics integrator manages the solution of an elastic 
      41             :     // solve using the MLMG solver. 
      42          25 :     static void Parse(Mechanics& value, IO::ParmParse& pp)
      43             :     {
      44             :         BL_PROFILE("Integrator::Base::Mechanics::Parse");
      45          25 :         if (pp.contains("type"))
      46             :         {
      47          18 :             std::string type_str;
      48             :             // Type of mecahnics to use.
      49             :             // Static: do full implicit solve.
      50             :             // Dynamic: evolve dynamic equations with explicit dynamics
      51             :             // Disable: do nothing.
      52           9 :             pp_query_validate("type", type_str, {"disable","static","dynamic"}); 
      53           9 :             if (type_str == "static")        value.m_type = Type::Static;
      54           3 :             else if (type_str == "dynamic")  value.m_type = Type::Dynamic;
      55           3 :             else if (type_str == "disable")  value.m_type = Type::Disable;
      56           0 :             else Util::Abort(INFO, "Invalid type ", type_str, " specified");
      57             :         }
      58          25 :         if (value.m_type == Type::Disable) return;
      59             : 
      60             :         // Treat mechanics fields as changing in time. [false]
      61             :         // You should use this if you care about other physics driven by
      62             :         // the output of this integrator.
      63          19 :         pp_query_default("time_evolving", value.m_time_evolving,false);
      64             : 
      65          19 :         pp_query_default("plot_disp", value.plot_disp, true); // Include displacement field in output
      66          19 :         pp_query_default("plot_rhs", value.plot_rhs, true); // Include right-hand side in output
      67          19 :         pp_query_default("plot_psi", value.plot_psi, true); // Include :math:`\psi` field in output
      68          19 :         pp_query_default("plot_stress", value.plot_stress, true); // Include stress in output
      69          19 :         pp_query_default("plot_strain", value.plot_strain, true); // Include strain in output
      70             : 
      71          19 :         value.RegisterGeneralFab(value.disp_mf, 1, 2, value.plot_disp, "disp", value.m_time_evolving);
      72          19 :         value.RegisterGeneralFab(value.rhs_mf, 1, 2, value.plot_rhs, "rhs", value.m_time_evolving);
      73          19 :         value.RegisterGeneralFab(value.stress_mf, 1, 2, value.plot_stress, "stress", value.m_time_evolving);
      74          19 :         value.RegisterGeneralFab(value.strain_mf, 1, 2, value.plot_strain, "strain", value.m_time_evolving);
      75             : 
      76          19 :         if (value.m_type == Type::Static)
      77             :         {
      78             :             // Read parameters for :ref:`Solver::Nonlocal::Newton` solver
      79          19 :             pp_queryclass("solver", value.solver);
      80             :         }
      81          19 :         if (value.m_type == Type::Dynamic)
      82             :         {
      83           0 :             value.RegisterGeneralFab(value.vel_mf, 1, 2, "vel",true);
      84           0 :             value.RegisterGeneralFab(value.disp_old_mf, 1, 2, "dispold");
      85           0 :             value.RegisterGeneralFab(value.vel_old_mf, 1, 2, "velold");
      86           0 :             value.RegisterGeneralFab(value.ddw_mf, 1, 2);
      87           0 :             pp_forbid("viscous.mu","replaced with viscous.mu_dashpot");
      88           0 :             pp_query_default("viscous.mu_dashpot", value.mu_dashpot,0.0); // Dashpot damping (damps velocity)
      89           0 :             pp_forbid("viscous.mu2","replaced with viscous.mu_newton");
      90           0 :             pp_query_default("viscous.mu_newton", value.mu_newton,0.0); // Newtonian viscous damping (damps velocity gradient)
      91             : 
      92           0 :             std::string velocity_ic_str;
      93           0 :             pp_query_validate("velocity.ic.type",velocity_ic_str,{"none","expression"}); // Initializer for RHS
      94           0 :             if (velocity_ic_str == "expression") value.velocity_ic = new IC::Expression(value.geom, pp,"velocity.ic.expression");
      95             :         }
      96             : 
      97             :         // Select the mechanical boundary conditions
      98          19 :         pp.select<BC::Operator::Elastic::Constant,BC::Operator::Elastic::TensionTest,BC::Operator::Elastic::Expression>("bc",value.bc);
      99             : 
     100          19 :         pp_query_default("print_model", value.m_print_model, false); // Print out model variables (if enabled by model)
     101          19 :         if (value.m_print_model) value.RegisterGeneralFab(value.model_mf, 1, 2, "model", value.m_time_evolving);
     102          19 :         else value.RegisterGeneralFab(value.model_mf, 1, 2, value.m_time_evolving);
     103             : 
     104             :         // This if is here because the code currently allows for no
     105             :         // body force to be defined, in which case value.ic_rhs is set to
     106             :         // a nullptr. This should eventually be fixed - this is a stopgap.
     107          19 :         if (pp.contains("rhs.type"))
     108             :             // initial condition for right hand side (body force) 
     109           4 :             pp.select<IC::Constant,IC::Expression,IC::Trig>("rhs",value.ic_rhs,value.geom);
     110             : 
     111             :         // Timestep interval for elastic solves (default - solve every time)
     112          19 :         pp_query_default("interval", value.m_interval, 0);
     113             : 
     114          19 :         value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[0]), "disp_xhi_x");
     115          19 :         value.RegisterIntegratedVariable(&(value.disp_hi[0].data()[1]), "disp_xhi_y");
     116          19 :         value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[0]), "disp_yhi_x");
     117          19 :         value.RegisterIntegratedVariable(&(value.disp_hi[1].data()[1]), "disp_yhi_y");
     118          19 :         value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[0]), "trac_xhi_x");
     119          19 :         value.RegisterIntegratedVariable(&(value.trac_hi[0].data()[1]), "trac_xhi_y");
     120          19 :         value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[0]), "trac_yhi_x");
     121          19 :         value.RegisterIntegratedVariable(&(value.trac_hi[1].data()[1]), "trac_yhi_y");
     122             : 
     123             :         // Maximum multigrid coarsening level (default - none, maximum coarsening)
     124          19 :         pp_query_default("max_coarsening_level", value.m_max_coarsening_level,-1);
     125             : 
     126             :         // Whether to include residual output field
     127          19 :         pp_query_default("print_residual", value.m_print_residual,false); 
     128          19 :         if (value.m_print_residual) value.RegisterGeneralFab(value.res_mf, 1, 2, "res", false);
     129             : 
     130             :         // Whether to refine based on elastic solution
     131          19 :         pp_query_default("elastic_ref_threshold", value.m_elastic_ref_threshold,0.01); 
     132             : 
     133             :         // Set this to true to zero out the displacement before each solve.
     134             :         // (This is a temporary fix - we need to figure out why this is needed.)
     135          19 :         pp_query_default("zero_out_displacement", value.m_zero_out_displacement,false);
     136             : 
     137             :         // Time to start doing the elastic solve (by default, start immediately)
     138          19 :         pp_query_default("tstart", value.tstart,-1.0);
     139             :     }
     140             : 
     141             : protected:
     142             :     /// \brief Use the #ic object to initialize#Temp
     143         121 :     void Initialize(int lev) override
     144             :     {
     145             :         BL_PROFILE("Integrator::Base::Mechanics::Initialize");
     146         121 :         if (m_type == Mechanics<MODEL>::Type::Disable) return;
     147             : 
     148          74 :         disp_mf[lev]->setVal(Set::Vector::Zero());
     149             :         //disp_old_mf[lev]->setVal(Set::Vector::Zero());
     150             : 
     151          74 :         if (m_type == Type::Dynamic)
     152           0 :             if (velocity_ic)
     153             :             {
     154           0 :                 velocity_ic->Initialize(lev,vel_mf);
     155             :             }
     156             : 
     157          74 :         if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
     158          67 :         else rhs_mf[lev]->setVal(Set::Vector::Zero());
     159             :     }
     160             : 
     161             :     virtual void UpdateModel(int a_step, Set::Scalar a_time) = 0;
     162             : 
     163         682 :     virtual void TimeStepBegin(Set::Scalar a_time, int a_step) override
     164             :     {
     165             :         BL_PROFILE("Integrator::Base::Mechanics::TimeStepBegin");
     166         682 :         if (m_type == Mechanics<MODEL>::Type::Disable) return;
     167             : 
     168         898 :         for (int lev = 0; lev <= finest_level; ++lev)
     169             :         {
     170         472 :             if (m_zero_out_displacement) disp_mf[lev]->setVal(Set::Vector::Zero());
     171         472 :             rhs_mf[lev]->setVal(Set::Vector::Zero());
     172         472 :             if (ic_rhs) ic_rhs->Initialize(lev, rhs_mf);
     173             :         }
     174             : 
     175         426 :         UpdateModel(a_step, a_time);
     176             : 
     177         426 :         bc->SetTime(a_time);
     178         426 :         bc->Init(rhs_mf, geom);
     179             : 
     180         426 :         if (m_type != Mechanics<MODEL>::Type::Static) return;
     181         426 :         if (a_time < tstart) return;
     182         426 :         if (m_interval && a_step % m_interval) return;
     183             : 
     184         426 :         amrex::LPInfo info;
     185         426 :         if (m_max_coarsening_level >= 0)
     186           0 :             info.setMaxCoarseningLevel(m_max_coarsening_level);
     187        1278 :         Operator::Elastic<MODEL::sym> elastic_op(Geom(0, finest_level), grids, DistributionMap(0, finest_level), info);
     188         426 :         elastic_op.SetUniform(false);
     189         426 :         elastic_op.SetHomogeneous(false);
     190         426 :         elastic_op.SetBC(bc);
     191        1278 :         IO::ParmParse pp("elasticop");
     192         426 :         pp_queryclass(elastic_op);
     193             : 
     194         426 :         Set::Scalar tol_rel = 1E-8, tol_abs = 1E-8;
     195             : 
     196         426 :         solver.Define(elastic_op);
     197         426 :         if (psi_on) solver.setPsi(psi_mf);
     198         426 :         solver.solve(disp_mf, rhs_mf, model_mf, tol_rel, tol_abs);
     199         426 :         if (m_print_residual) solver.compLinearSolverResidual(res_mf, disp_mf, rhs_mf);
     200         426 :         solver.Clear();
     201             : 
     202         898 :         for (int lev = 0; lev <= disp_mf.finest_level; lev++)
     203             :         {
     204         472 :             amrex::Box domain = geom[lev].Domain();
     205         472 :             domain.convert(amrex::IntVect::TheNodeVector());
     206             : 
     207         472 :             const amrex::Real* DX = geom[lev].CellSize();
     208        1077 :             for (MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
     209             :             {
     210         605 :                 amrex::Box bx = mfi.nodaltilebox();
     211         605 :                 bx.grow(2);
     212         605 :                 bx = bx & domain;
     213         605 :                 amrex::Array4<MODEL>             const& model = model_mf[lev]->array(mfi);
     214         605 :                 amrex::Array4<Set::Matrix>       const& stress = stress_mf[lev]->array(mfi);
     215         605 :                 amrex::Array4<Set::Matrix>       const& strain = strain_mf[lev]->array(mfi);
     216         605 :                 amrex::Array4<const Set::Vector> const& disp = disp_mf[lev]->array(mfi);
     217             : 
     218             : 
     219      187862 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     220             :                 {
     221      187257 :                     auto sten = Numeric::GetStencil(i, j, k, bx);
     222      374514 :                     if (model(i, j, k).kinvar == Model::Solid::KinematicVariable::F)
     223             :                     {
     224       97068 :                         Set::Matrix F = Set::Matrix::Identity() + Numeric::Gradient(disp, i, j, k, DX, sten);
     225      145602 :                         stress(i, j, k) = model(i, j, k).DW(F);
     226       97068 :                         strain(i, j, k) = F;
     227             :                     }
     228             :                     else
     229             :                     {
     230      138723 :                         Set::Matrix gradu = Numeric::Gradient(disp, i, j, k, DX, sten);
     231      416169 :                         stress(i, j, k) = model(i, j, k).DW(gradu);
     232      277446 :                         strain(i, j, k) = 0.5 * (gradu + gradu.transpose());
     233             :                     }
     234             :                 });
     235             :             }
     236         472 :             Util::RealFillBoundary(*stress_mf[lev], geom[lev]);
     237         472 :             Util::RealFillBoundary(*strain_mf[lev], geom[lev]);
     238             :         }
     239             :     }
     240             : 
     241       55102 :     void Advance(int lev, Set::Scalar time, Set::Scalar dt) override
     242             :     {
     243             :         BL_PROFILE("Integrator::Base::Mechanics::Advance");
     244       55102 :         if (m_type == Mechanics<MODEL>::Type::Disable) return;
     245         718 :         const amrex::Real* DX = geom[lev].CellSize();
     246             : 
     247         718 :         amrex::Box domain = geom[lev].Domain();
     248        1436 :         domain.convert(amrex::IntVect::TheNodeVector());
     249         718 :         const amrex::Dim3 lo = amrex::lbound(domain), hi = amrex::ubound(domain);
     250             : 
     251         718 :         if (m_type == Type::Dynamic)
     252             :         {
     253             : 
     254           0 :             std::swap(*disp_mf[lev], *disp_old_mf[lev]);
     255           0 :             std::swap(*vel_mf[lev], *vel_old_mf[lev]);
     256             : 
     257             : 
     258           0 :             for (amrex::MFIter mfi(*disp_mf[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     259             :             {
     260           0 :                 Set::Patch<const Set::Vector> u     = disp_old_mf.Patch(lev,mfi);
     261           0 :                 Set::Patch<const Set::Vector> v     = vel_old_mf.Patch(lev,mfi);
     262           0 :                 Set::Patch<const Set::Vector> b     = rhs_mf.Patch(lev,mfi);
     263             : 
     264           0 :                 Set::Patch<Set::Vector>       unew  = disp_mf.Patch(lev,mfi);
     265           0 :                 Set::Patch<Set::Vector>       vnew  = vel_mf.Patch(lev,mfi);
     266           0 :                 Set::Patch<Set::Matrix>       eps   = strain_mf.Patch(lev,mfi);
     267           0 :                 Set::Patch<Set::Matrix>       sig   = stress_mf.Patch(lev,mfi);
     268             :                 //Set::Patch<MATRIX4>           ddw   = ddw_mf.Patch(lev,mfi);
     269           0 :                 Set::Patch<MODEL>             model = model_mf.Patch(lev,mfi);
     270             : 
     271           0 :                 amrex::Box bx = mfi.grownnodaltilebox() & domain;
     272           0 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     273             :                 {
     274           0 :                     auto sten = Numeric::GetStencil(i,j,k,bx);
     275           0 :                     eps(i, j, k) = Numeric::Gradient(u, i, j, k, DX,sten);
     276           0 :                     sig(i, j, k) = model(i, j, k).DW(eps(i, j, k));
     277             :                     //ddw(i,j,k) = model(i,j,k).DDW(eps(i,j,k));
     278             :                 });
     279             : 
     280           0 :                 bx = mfi.nodaltilebox() & domain;
     281           0 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     282             :                 {
     283             : 
     284           0 :                     bool AMREX_D_DECL(xmin = (i == lo.x), ymin = (j == lo.y), zmin = (k == lo.z));
     285           0 :                     bool AMREX_D_DECL(xmax = (i == hi.x), ymax = (j == hi.y), zmax = (k == hi.z));
     286             : 
     287           0 :                     if (AMREX_D_TERM(xmax || xmin, || ymax || ymin, || zmax || zmin))
     288             :                     {
     289           0 :                         auto sten = Numeric::GetStencil(i,j,k,domain);                        
     290             : 
     291           0 :                         auto bctype = bc->getType(i,j,k,domain);
     292             : 
     293           0 :                         for (int d = 0; d < AMREX_SPACEDIM; d++)
     294             :                         {
     295           0 :                             if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Displacement)
     296             :                             {
     297           0 :                                 unew(i,j,k)(d) = b(i,j,k)(d);
     298             :                             }
     299           0 :                             else if (bctype[d] == BC::Operator::Elastic::Elastic::Type::Traction)
     300             :                             {
     301             :                                 
     302           0 :                                 Set::Vector N = Set::Normal(AMREX_D_DECL(xmin,ymin,zmin),
     303             :                                                             AMREX_D_DECL(xmax,ymax,zmax));
     304             :                                 
     305           0 :                                 auto [phi, offdiag] = Numeric::GradientSplit(u,i,j,k,DX,sten);
     306             :                                 
     307           0 :                                 Set::Matrix A = Set::Matrix::Zero();
     308           0 :                                 Set::Vector rhs = b(i,j,k);
     309           0 :                                 Set::Matrix DW_F0 = model(i,j,k).DW(Set::Matrix::Zero());
     310           0 :                                 MATRIX4 ddw = model(i,j,k).DDW(eps(i,j,k));
     311             :                                 //Util::Message(INFO,b(i,j,k).transpose());
     312             : 
     313           0 :                                 for (int p = 0; p < AMREX_SPACEDIM; p++)
     314           0 :                                     for (int q = 0; q < AMREX_SPACEDIM; q++)
     315             :                                     {
     316           0 :                                         for (int r = 0; r < AMREX_SPACEDIM; r++)
     317           0 :                                             for (int s = 0; s < AMREX_SPACEDIM; s++)
     318             :                                             {
     319           0 :                                                 A(p,r) += ddw(p,q,r,s) * phi(s) * N(q);
     320             : 
     321           0 :                                                 rhs(p) -= ddw(p,q,r,s) * eps(i,j,k)(r,s) * N(q);
     322             :                                             }
     323             :                                         
     324           0 :                                         rhs(p) -= DW_F0(p,q)*N(q);
     325             :                                     }
     326           0 :                                 Set::Vector delta_u = (A.inverse() * rhs);
     327             : 
     328           0 :                                 unew(i,j,k)(d) = u(i,j,k)(d) + delta_u(d);
     329           0 :                                 vnew(i,j,k)(d) = (unew(i,j,k)(d) - u(i,j,k)(d))/dt;
     330             :                             }
     331             :                             else
     332             :                             {
     333           0 :                                 Util::Abort(INFO,"Elastic dynamics not supported for other BCs yet");
     334             :                             }
     335           0 :                         }
     336             :                     }
     337             :                     else
     338             :                     {
     339             :                         //Set::Matrix      gradu = Numeric::Gradient(u,i,j,k,DX);
     340             :                         //Set::Matrix3 gradgradu = Numeric::Hessian(u,i,j,k,DX);
     341             :                         
     342             :                         //Set::Vector f = ddw(i,j,k)*gradgradu;
     343           0 :                         Set::Vector f = Numeric::Divergence(sig, i, j, k, DX);
     344             : 
     345             :                         //MATRIX4 AMREX_D_DECL(
     346             :                         //    Cgrad1 = (Numeric::Stencil<MATRIX4, 1, 0, 0>::D(ddw, i, j, k, 0, DX)),
     347             :                         //    Cgrad2 = (Numeric::Stencil<MATRIX4, 0, 1, 0>::D(ddw, i, j, k, 0, DX)),
     348             :                         //    Cgrad3 = (Numeric::Stencil<MATRIX4, 0, 0, 1>::D(ddw, i, j, k, 0, DX)));
     349             : 
     350             :                         //f += AMREX_D_TERM(  ( Cgrad1*gradu).col(0),
     351             :                         //                    +(Cgrad2*gradu).col(1),
     352             :                         //                    +(Cgrad3*gradu).col(2));
     353             :                         
     354           0 :                         Set::Vector lapv = Numeric::Laplacian(v, i, j, k, DX);
     355           0 :                         f += mu_newton * lapv + b(i,j,k) - mu_dashpot*v(i,j,k);
     356             :                         
     357           0 :                         Set::Vector udotdot = f / rho;
     358           0 :                         vnew(i, j, k) = v(i, j, k) + dt * udotdot;
     359           0 :                         unew(i, j, k) = u(i, j, k) + dt * v(i, j, k);
     360             :                     }
     361             :                 });
     362             :             }
     363             :         }
     364             : 
     365        4424 :         for (amrex::MFIter mfi(*disp_mf[lev], false); mfi.isValid(); ++mfi)
     366             :         {
     367        3706 :             amrex::Box bx = mfi.grownnodaltilebox() & domain;
     368        3706 :             amrex::Array4<Set::Matrix>         const& eps = (*strain_mf[lev]).array(mfi);
     369        3706 :             amrex::Array4<Set::Matrix>         const& sig = (*stress_mf[lev]).array(mfi);
     370        3706 :             amrex::Array4<MODEL>               const& model = (*model_mf[lev]).array(mfi);
     371     7012752 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     372             :             {
     373     7059530 :                 model(i, j, k).Advance(dt, eps(i, j, k), sig(i, j, k),time);
     374             :             });
     375             :         }
     376         718 :         Util::RealFillBoundary(*model_mf[lev], geom[lev]);
     377             : 
     378             :     }
     379             : 
     380       12754 :     void Integrate(int amrlev, Set::Scalar /*time*/, int /*step*/,
     381             :         const amrex::MFIter& mfi, const amrex::Box& a_box) override
     382             :     {
     383             :         BL_PROFILE("Integrator::Base::Mechanics::Integrate");
     384       12754 :         if (m_type == Type::Disable) return;
     385             : 
     386         856 :         if (amrex::ParallelDescriptor::NProcs() > 1)
     387             :         {
     388           0 :             Util::Warning(INFO,"There is a known bug when calculating trac/disp in Base::Mechanics in parallel.");
     389           0 :             Util::Warning(INFO,"The thermo.dat values likely will not be correct; use the boxlib output instead.");
     390             :         }
     391             : 
     392         856 :         const amrex::Real* DX = geom[amrlev].CellSize();
     393         856 :         amrex::Box domain = geom[amrlev].Domain();
     394         856 :         domain.convert(amrex::IntVect::TheNodeVector());
     395             : 
     396         856 :         amrex::Box box = a_box;
     397         856 :         box.convert(amrex::IntVect::TheNodeVector());
     398             : 
     399             : 
     400             :         //Set::Scalar dv = AMREX_D_TERM(DX[0], *DX[1], *DX[2]);
     401             : #if AMREX_SPACEDIM == 2
     402         556 :         Set::Vector da0(DX[1], 0);
     403         556 :         Set::Vector da1(0, DX[0]);
     404             : #elif AMREX_SPACEDIM == 3
     405         300 :         Set::Vector da(DX[1] * DX[2], 0, 0);
     406             : #endif
     407             : 
     408         856 :         const Dim3 /*lo= amrex::lbound(domain),*/ hi = amrex::ubound(domain);
     409         856 :         const Dim3 /*boxlo= amrex::lbound(box),*/ boxhi = amrex::ubound(box);
     410             : 
     411         856 :         amrex::Array4<const Set::Matrix> const& stress = (*stress_mf[amrlev]).array(mfi);
     412         856 :         amrex::Array4<const Set::Vector> const& disp = (*disp_mf[amrlev]).array(mfi);
     413      136248 :         amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     414             :         {
     415             : #if AMREX_SPACEDIM == 2
     416       97892 :             if (i == hi.x && j < boxhi.y)
     417             :             {
     418        3072 :                 trac_hi[0] += (0.5 * (stress(i, j, k) + stress(i, j + 1, k)) * da0);
     419        2048 :                 disp_hi[0] = disp(i, j, k);
     420             :             }
     421       97892 :             if (j == hi.y && i < boxhi.x)
     422             :             {
     423        3264 :                 trac_hi[1] += (0.5 * (stress(i, j, k) + stress(i + 1, j, k)) * da1);
     424        2176 :                 disp_hi[1] = disp(i, j, k);
     425             :             }
     426             : #elif AMREX_SPACEDIM == 3
     427       37500 :             if (i == hi.x && (j < boxhi.y && k < boxhi.z))
     428             :             {
     429        9600 :                 trac_hi[0] += (0.25 * (stress(i, j, k) + stress(i, j + 1, k)
     430       14400 :                     + stress(i, j, k + 1) + stress(i, j + 1, k + 1)) * da);
     431        9600 :                 disp_hi[0] = disp(i, j, k);
     432             :             }
     433             : #endif
     434             :         });
     435             : 
     436             :     }
     437             : 
     438         528 :     void TagCellsForRefinement(int lev, amrex::TagBoxArray& a_tags, Set::Scalar /*time*/, int /*ngrow*/) override
     439             :     {
     440             :         BL_PROFILE("Integrator::Base::Mechanics::TagCellsForRefinement");
     441         528 :         if (m_type == Type::Disable) return;
     442             : 
     443         215 :         Set::Vector DX(geom[lev].CellSize());
     444         215 :         Set::Scalar DXnorm = DX.lpNorm<2>();
     445        1202 :         for (amrex::MFIter mfi(*strain_mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi)
     446             :         {
     447         987 :             amrex::Box bx = mfi.tilebox();
     448         987 :             bx.convert(amrex::IntVect::TheCellVector());
     449         987 :             amrex::Array4<char> const& tags = a_tags.array(mfi);
     450         987 :             amrex::Array4<Set::Matrix> const& eps = strain_mf[lev]->array(mfi);
     451      371675 :             amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
     452             :             {
     453      741376 :                 Set::Matrix3 grad = Numeric::NodeGradientOnCell(eps, i, j, k, DX.data());
     454      370688 :                 if (grad.norm() * DXnorm > m_elastic_ref_threshold)
     455         412 :                     tags(i, j, k) = amrex::TagBox::SET;
     456             :             });
     457             :         }
     458             :     }
     459             : 
     460             : protected:
     461             :     typedef Set::Matrix4<AMREX_SPACEDIM,MODEL::sym> MATRIX4;
     462             :     Set::Field<MODEL> model_mf;
     463             :     Set::Field<MATRIX4> ddw_mf;
     464             :     Set::Field<Set::Scalar> psi_mf;
     465             :     bool psi_on = false;
     466             : 
     467             :     int m_interval = 0;
     468             :     Type m_type = Type::Static;
     469             : 
     470             :     Set::Field<Set::Vector> disp_mf;
     471             :     Set::Field<Set::Vector> rhs_mf;
     472             :     Set::Field<Set::Vector> res_mf;
     473             :     Set::Field<Set::Matrix> stress_mf;
     474             :     Set::Field<Set::Matrix> strain_mf;
     475             : 
     476             :     // Only use these if using the "dynamics" option
     477             :     Set::Field<Set::Vector> disp_old_mf;
     478             :     Set::Field<Set::Vector> vel_mf;
     479             :     Set::Field<Set::Vector> vel_old_mf;
     480             :     //Set::Field<Set::Matrix4<AMREX_SPACEDIM,MODEL::sym>> ddw_mf;
     481             :     Set::Scalar rho = 1.0;
     482             :     Set::Scalar mu_dashpot = NAN;
     483             :     Set::Scalar mu_newton = NAN;
     484             : 
     485             :     //Set::Vector trac_lo[AMREX_SPACEDIM];
     486             :     Set::Vector trac_hi[AMREX_SPACEDIM];
     487             :     Set::Vector disp_hi[AMREX_SPACEDIM];
     488             : 
     489             : 
     490             :     IC::IC* ic_rhs = nullptr;
     491             : 
     492             :     IC::IC* velocity_ic = nullptr;
     493             : 
     494             :     Solver::Nonlocal::Newton<MODEL> solver;//(elastic.op);
     495             :     BC::Operator::Elastic::Elastic* bc = nullptr;
     496             : 
     497             :     Set::Scalar m_elastic_ref_threshold = 0.01;
     498             :     bool m_print_model = false;
     499             :     bool m_print_residual = false;
     500             :     bool m_time_evolving = false;
     501             :     int m_max_coarsening_level = -1;
     502             : 
     503             :     bool m_zero_out_displacement = false;
     504             : 
     505             :     bool plot_disp = true;
     506             :     bool plot_stress = true;
     507             :     bool plot_strain = true;
     508             :     bool plot_psi = true;
     509             :     bool plot_rhs = true;
     510             : 
     511             : 
     512             :     Set::Scalar tstart = -1;
     513             : };
     514             : }
     515             : }
     516             : #endif

Generated by: LCOV version 1.14