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 749 : 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 1803 : void SetModel (const Set::Field<Set::Matrix4<AMREX_SPACEDIM,SYM>> & a_model)
54 3799 : { 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 749 : void SetBC (::BC::Operator::Elastic::Elastic *a_bc)
62 : {
63 749 : if (this->Geom(0).isAnyPeriodic())
64 : {
65 0 : Util::Warning(INFO,"Looks like you're using a periodic domain. \nThat is currently VERY DANGEROUS when using linear elastic solver!!!");
66 : }
67 749 : m_bc = a_bc;
68 749 : m_bc_set = true;
69 749 : };
70 130596 : ::BC::Operator::Elastic::Elastic & GetBC()
71 : {
72 130596 : return *m_bc;
73 : }
74 :
75 : /// Compute strain \f$\mathbf{\epsilon}\f$ given the displacement field \f$\mathbf{u}\f$
76 : /// by
77 : /// \f[\mathbf{\epsilon}_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i})\f]
78 : ///
79 : void Strain (int amrlev, amrex::MultiFab& epsfab, const amrex::MultiFab& ufab, bool voigt = false) const;
80 :
81 : /// Compute stress \f$\mathbf{\sigma}\f$ given the displacement field \f$\mathbf{u}\f$
82 : /// by
83 : /// \f[\mathbf{\sigma}_{ij} = \mathbb{C}_{ijkl}\,u_{k,l}\f]
84 : /// where, \f$\mathbb{C}\f$ is furnished by the templated `model`
85 : ///
86 : void Stress (int amrlev, amrex::MultiFab& sigmafab, const amrex::MultiFab& ufab, bool voigt = false, bool a_homogeneous=false);
87 :
88 : /// Compute energy density \f$\mathbb{W}\f$ given the displacement field \f$\mathbf{u}\f$
89 : /// by
90 : /// \f[\mathbb{W} = \mathbb{\Sigma}\mathbb{\Sigma}\mathbf{\sigma}_{ij}\mathbf{\epsilon}_{ij}\f]
91 : /// where \f$\mathbf{\sigma}\f$ and \f$\mathbf{\epsilon}\f$ are as defined above
92 : ///
93 : void Energy (int amrlev, amrex::MultiFab& energy, const amrex::MultiFab& u, bool a_homogeneous=false);
94 :
95 : //void Energy (int amrlev, amrex::MultiFab& energies, const amrex::MultiFab& u, std::vector<SOLID> models, bool a_homogeneous=false);
96 :
97 : /// This function is depricated and should not be used. Use the other `SetBC` function.
98 : ///
99 0 : void SetBC(const std::array<std::array<BC,AMREX_SPACEDIM>,AMREX_SPACEDIM> &a_bc_lo,
100 : const std::array<std::array<BC,AMREX_SPACEDIM>,AMREX_SPACEDIM> &a_bc_hi)
101 : {
102 0 : Util::Abort(INFO,"This is now depricated. Use the other SetBC function.");
103 0 : m_bc_lo = a_bc_lo;m_bc_hi = a_bc_hi;
104 0 : };
105 :
106 : void Error0x (int amrlev, int mglev, MultiFab& R0x, const MultiFab& x) const;
107 :
108 0 : void SetTesting(bool a_testing) {m_testing = a_testing;}
109 749 : void SetUniform(bool a_uniform) {m_uniform = a_uniform;}
110 0 : virtual void SetAverageDownCoeffs(bool a_average_down_coeffs) override
111 0 : {m_average_down_coeffs = a_average_down_coeffs;}
112 :
113 : using Operator::Reflux;
114 :
115 : protected:
116 :
117 : virtual void Diagonal (int amrlev, int mglev, amrex::MultiFab& diag) override;
118 :
119 : virtual void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const override final;
120 : virtual void FFlux (int amrlev, const MFIter& mfi,
121 : const std::array<FArrayBox*,AMREX_SPACEDIM>& flux,
122 : const FArrayBox& sol, const int face_only=0) const final;
123 :
124 :
125 :
126 :
127 920882 : virtual int getNComp() const override {return AMREX_SPACEDIM;};
128 0 : virtual bool isCrossStencil () const { return false; }
129 748 : virtual void prepareForSolve () override
130 : {
131 748 : if (!m_model_set) Util::Abort(INFO,"Attempting to use operator before calling SetModel!");
132 748 : if (!m_bc_set) Util::Warning(INFO,"Attempting to use operator before calling SetBC!");
133 748 : Operator<Grid::Node>::prepareForSolve();
134 747 : };
135 :
136 : private:
137 : /// Simple arrays storing boundary conditions for each component and each face.
138 : std::array<std::array<BC,AMREX_SPACEDIM>, AMREX_SPACEDIM> m_bc_lo; // m_bc_lo[face][dimension]
139 : std::array<std::array<BC,AMREX_SPACEDIM>, AMREX_SPACEDIM> m_bc_hi; // m_bc_hi[face][dimension]
140 :
141 : /// This is a multifab-type object containing objects of type
142 : /// Model::Solid::Elastic::Isotropic::Isotropic
143 : /// (or some other model type). T is the template argument.
144 : /// The models contain elastic constants and contain methods for converting strain to stress
145 : amrex::Vector<Set::Field<Set::Matrix4<AMREX_SPACEDIM,SYM>>> m_ddw_mf;
146 :
147 : // This is a mask variable.
148 : amrex::Vector<Set::Field<Set::Scalar>> m_psi_mf;
149 : Set::Scalar m_psi_small = 1E-8;
150 : bool m_psi_set = false;
151 :
152 : virtual void averageDownCoeffs () override;
153 : void averageDownCoeffsDifferentAmrLevels (int fine_amrlev);
154 :
155 : /// \fn averageDownCoeffsSameAmrLevel
156 : /// \brief Update coarse-level AMR coefficients with data from fine level
157 : ///
158 : /// This function is called before the solve, and it updates the Matrix4 coefficients
159 : /// on coarse levels with averaged versions on the fine levels.
160 : /// Typically, this is not necessary, since the coefficients are computed on each level
161 : /// already.
162 : /// However, there are some cases where this messes up the multigrid solver.
163 : /// (Phase field fracture is the primary example.)
164 : ///
165 : /// This is disabled by default. In order to enable, call
166 : ///
167 : /// elasticoperator.SetAverageDownCoeffs(true);
168 : ///
169 : void averageDownCoeffsSameAmrLevel (int amrlev);
170 :
171 : void FillBoundaryCoeff (MultiTab& sigma, const Geometry& geom);
172 : void FillBoundaryCoeff (MultiFab& psi, const Geometry& geom);
173 :
174 : bool m_testing = false;
175 : bool m_uniform = false;
176 : bool m_homogeneous = false;
177 :
178 : bool m_average_down_coeffs = false;
179 :
180 : ::BC::Operator::Elastic::Elastic *m_bc;
181 :
182 : bool m_model_set = false;
183 : bool m_bc_set = false;
184 :
185 : public:
186 749 : static void Parse(Elastic<SYM> & value, IO::ParmParse & pp)
187 : {
188 : // Regularization offset value used in near-singular elastic solves.
189 : // It should be small - if it is too large, you will get better convergence
190 : // but less correct values!
191 749 : pp_query("small",value.m_psi_small);
192 749 : }
193 :
194 : };
195 : }
196 :
197 : #endif
|
