Line data Source code
1 : #ifndef OPERATOR_ELASTIC_H_
2 : #define OPERATOR_ELASTIC_H_
3 :
4 : #include <AMReX_MLCellLinOp.H>
5 : #include <AMReX_Array.H>
6 : #include <limits>
7 : #include "Set/Set.H"
8 : #include "Operator/Operator.H"
9 : #include "Model/Solid/Solid.H"
10 : #include "BC/Operator/Elastic/Elastic.H"
11 : #include "Numeric/Stencil.H"
12 : #include "IC/IC.H"
13 :
14 : using namespace amrex;
15 :
16 : namespace Operator
17 : {
18 : ///
19 : /// \todo Remove Elastic derived classes. They have been replaced with the Model construction.
20 : ///
21 : template<int SYM>
22 : class Elastic : public Operator<Grid::Node>
23 : {
24 : using MATRIX4 = Set::Matrix4<AMREX_SPACEDIM,SYM>;
25 : using TArrayBox = amrex::BaseFab<MATRIX4>;
26 : using MultiTab = amrex::FabArray<TArrayBox>;
27 : public:
28 : enum class BC {Displacement, Traction, Periodic, Neumann};
29 : enum class Boundary {Lo, Hi, None};
30 :
31 0 : Elastic () {}
32 : Elastic (const Vector<Geometry>& a_geom,
33 : const Vector<BoxArray>& a_grids,
34 : const Vector<DistributionMapping>& a_dmap,
35 : const LPInfo& a_info);
36 : virtual ~Elastic ();
37 : Elastic (const Elastic&) = delete;
38 : Elastic (Elastic&&) = delete;
39 : Elastic& operator= (const Elastic&) = delete;
40 : Elastic& operator= (Elastic&&) = delete;
41 :
42 : void define (const Vector<Geometry>& a_geom,
43 : const Vector<BoxArray>& a_grids,
44 : const Vector<DistributionMapping>& a_dmap,
45 : const LPInfo& a_info = LPInfo(),
46 : const Vector<FabFactory<FArrayBox> const*>& a_factory = {});
47 :
48 734 : virtual void SetHomogeneous (bool a_homogeneous) override {m_homogeneous = a_homogeneous;}
49 : void SetModel (Set::Matrix4<AMREX_SPACEDIM,SYM> &a_model);
50 : void SetModel (int amrlev, const MultiTab& a_model);
51 0 : void SetModel (const amrex::Vector<MultiTab> & a_model)
52 0 : { for (int ilev = 0; ilev < a_model.size(); ilev++) SetModel(ilev,a_model[ilev]);}
53 1634 : void SetModel (const Set::Field<Set::Matrix4<AMREX_SPACEDIM,SYM>> & a_model)
54 3322 : { for (int ilev = 0; ilev < a_model.size(); ilev++) SetModel(ilev,*a_model[ilev]);}
55 : void SetPsi (int amrlev, const amrex::MultiFab& a_psi);
56 0 : void SetPsi (int amrlev, const amrex::MultiFab& a_psi, const Set::Scalar &a_psi_small)
57 0 : {m_psi_small = a_psi_small; SetPsi(amrlev,a_psi);}
58 :
59 : /// The different types of Boundary Condtiions are listed in the `BC::Operator::Elastic` documentation
60 : ///
61 734 : void SetBC (::BC::Operator::Elastic::Elastic *a_bc)
62 : {
63 734 : m_bc = a_bc;
64 734 : m_bc_set = true;
65 734 : };
66 80368 : ::BC::Operator::Elastic::Elastic & GetBC()
67 : {
68 80368 : return *m_bc;
69 : }
70 :
71 : /// Compute strain \f$\mathbf{\epsilon}\f$ given the displacement field \f$\mathbf{u}\f$
72 : /// by
73 : /// \f[\mathbf{\epsilon}_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})\f]
74 : ///
75 : void Strain (int amrlev, amrex::MultiFab& epsfab, const amrex::MultiFab& ufab, bool voigt = false) const;
76 :
77 : /// Compute stress \f$\mathbf{\sigma}\f$ given the displacement field \f$\mathbf{u}\f$
78 : /// by
79 : /// \f[\mathbf{\sigma}_{ij} = \mathbb{C}_{ijkl}\,u_{k,l}\f]
80 : /// where, \f$\mathbb{C}\f$ is furnished by the templated `model`
81 : ///
82 : void Stress (int amrlev, amrex::MultiFab& sigmafab, const amrex::MultiFab& ufab, bool voigt = false, bool a_homogeneous=false);
83 :
84 : /// Compute energy density \f$\mathbb{W}\f$ given the displacement field \f$\mathbf{u}\f$
85 : /// by
86 : /// \f[\mathbb{W} = \mathbb{\Sigma}\mathbb{\Sigma}\mathbf{\sigma}_{ij}\mathbf{\epsilon}_{ij}\f]
87 : /// where \f$\mathbf{\sigma}\f$ and \f$\mathbf{\epsilon}\f$ are as defined above
88 : ///
89 : void Energy (int amrlev, amrex::MultiFab& energy, const amrex::MultiFab& u, bool a_homogeneous=false);
90 :
91 : //void Energy (int amrlev, amrex::MultiFab& energies, const amrex::MultiFab& u, std::vector<SOLID> models, bool a_homogeneous=false);
92 :
93 : /// This function is depricated and should not be used. Use the other `SetBC` function.
94 : ///
95 0 : void SetBC(const std::array<std::array<BC,AMREX_SPACEDIM>,AMREX_SPACEDIM> &a_bc_lo,
96 : const std::array<std::array<BC,AMREX_SPACEDIM>,AMREX_SPACEDIM> &a_bc_hi)
97 : {
98 0 : Util::Abort(INFO,"This is now depricated. Use the other SetBC function.");
99 0 : m_bc_lo = a_bc_lo;m_bc_hi = a_bc_hi;
100 0 : };
101 :
102 : void Error0x (int amrlev, int mglev, MultiFab& R0x, const MultiFab& x) const;
103 :
104 0 : void SetTesting(bool a_testing) {m_testing = a_testing;}
105 734 : void SetUniform(bool a_uniform) {m_uniform = a_uniform;}
106 0 : virtual void SetAverageDownCoeffs(bool a_average_down_coeffs) override
107 0 : {m_average_down_coeffs = a_average_down_coeffs;}
108 :
109 : using Operator::Reflux;
110 :
111 : protected:
112 :
113 : virtual void Diagonal (int amrlev, int mglev, amrex::MultiFab& diag) override;
114 :
115 : virtual void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const override final;
116 : virtual void FFlux (int amrlev, const MFIter& mfi,
117 : const std::array<FArrayBox*,AMREX_SPACEDIM>& flux,
118 : const FArrayBox& sol, const int face_only=0) const final;
119 :
120 :
121 :
122 :
123 262435 : virtual int getNComp() const override {return AMREX_SPACEDIM;};
124 0 : virtual bool isCrossStencil () const { return false; }
125 734 : virtual void prepareForSolve () override
126 : {
127 734 : if (!m_model_set) Util::Abort(INFO,"Attempting to use operator before calling SetModel!");
128 734 : if (!m_bc_set) Util::Warning(INFO,"Attempting to use operator before calling SetBC!");
129 734 : Operator<Grid::Node>::prepareForSolve();
130 734 : };
131 :
132 : private:
133 : /// Simple arrays storing boundary conditions for each component and each face.
134 : std::array<std::array<BC,AMREX_SPACEDIM>, AMREX_SPACEDIM> m_bc_lo; // m_bc_lo[face][dimension]
135 : std::array<std::array<BC,AMREX_SPACEDIM>, AMREX_SPACEDIM> m_bc_hi; // m_bc_hi[face][dimension]
136 :
137 : /// This is a multifab-type object containing objects of type
138 : /// Model::Solid::Elastic::Isotropic::Isotropic
139 : /// (or some other model type). T is the template argument.
140 : /// The models contain elastic constants and contain methods for converting strain to stress
141 : amrex::Vector<Set::Field<Set::Matrix4<AMREX_SPACEDIM,SYM>>> m_ddw_mf;
142 :
143 : // This is a mask variable.
144 : amrex::Vector<Set::Field<Set::Scalar>> m_psi_mf;
145 : Set::Scalar m_psi_small = 1E-8;
146 : bool m_psi_set = false;
147 :
148 : virtual void averageDownCoeffs () override;
149 : void averageDownCoeffsDifferentAmrLevels (int fine_amrlev);
150 :
151 : /// \fn averageDownCoeffsSameAmrLevel
152 : /// \brief Update coarse-level AMR coefficients with data from fine level
153 : ///
154 : /// This function is called before the solve, and it updates the Matrix4 coefficients
155 : /// on coarse levels with averaged versions on the fine levels.
156 : /// Typically, this is not necessary, since the coefficients are computed on each level
157 : /// already.
158 : /// However, there are some cases where this messes up the multigrid solver.
159 : /// (Phase field fracture is the primary example.)
160 : ///
161 : /// This is disabled by default. In order to enable, call
162 : ///
163 : /// elasticoperator.SetAverageDownCoeffs(true);
164 : ///
165 : void averageDownCoeffsSameAmrLevel (int amrlev);
166 :
167 : void FillBoundaryCoeff (MultiTab& sigma, const amrex::Periodicity& p);
168 : void FillBoundaryCoeff (MultiFab& psi, const amrex::Periodicity& p);
169 :
170 : bool m_testing = false;
171 : bool m_uniform = false;
172 : bool m_homogeneous = false;
173 :
174 : bool m_average_down_coeffs = false;
175 :
176 : ::BC::Operator::Elastic::Elastic *m_bc;
177 :
178 : bool m_model_set = false;
179 : bool m_bc_set = false;
180 :
181 : public:
182 734 : static void Parse(Elastic<SYM> & value, IO::ParmParse & pp)
183 : {
184 : // Regularization offset value used in near-singular elastic solves.
185 : // It should be small - if it is too large, you will get better convergence
186 : // but less correct values!
187 734 : pp_query("small",value.m_psi_small);
188 734 : }
189 :
190 : };
191 : }
192 :
193 : #endif
|
