Line data Source code
1 : #ifndef MODEL_SOLID_FINITE_PSEUDOAFFINE_CUBIC_H_
2 : #define MODEL_SOLID_FINITE_PSEUDOAFFINE_CUBIC_H_
3 :
4 : #include "IO/ParmParse.H"
5 : #include "Model/Solid/Finite/PseudoLinear/Cubic.H"
6 :
7 : namespace Model
8 : {
9 : namespace Solid
10 : {
11 : namespace Finite
12 : {
13 : namespace PseudoAffine
14 : {
15 : class Cubic : public PseudoLinear::Cubic
16 : {
17 : public:
18 : static constexpr const char *name = "cubic";
19 :
20 4 : Cubic() {};
21 66 : Cubic(PseudoLinear::Cubic base) : PseudoLinear::Cubic(base) {};
22 96 : virtual ~Cubic() {};
23 :
24 480 : Set::Scalar W(const Set::Matrix& F) const override
25 : {
26 480 : Set::Matrix F0inv = F0.inverse();
27 480 : return PseudoLinear::Cubic::W(F * F0inv);
28 : }
29 1960 : Set::Matrix DW(const Set::Matrix& F) const override
30 : {
31 1960 : Set::Matrix F0inv = F0.inverse();
32 1960 : return PseudoLinear::Cubic::DW(F * F0inv) * F0inv.transpose();
33 : }
34 20 : Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Major> DDW(const Set::Matrix& F) const override
35 : {
36 20 : Set::Matrix F0inv = F0.inverse();
37 20 : auto DDW = PseudoLinear::Cubic::DDW(F * F0inv);
38 20 : auto ret = Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Major>::Zero();
39 70 : for (int i = 0; i < AMREX_SPACEDIM; i++)
40 180 : for (int J = 0; J < AMREX_SPACEDIM; J++)
41 480 : for (int k = 0; k < AMREX_SPACEDIM; k++)
42 1320 : for (int L = 0; L < AMREX_SPACEDIM; L++)
43 : {
44 970 : ret(i,J,k,L) = 0.0;
45 3720 : for (int Q = 0; Q < AMREX_SPACEDIM; Q++)
46 10680 : for (int R = 0; R < AMREX_SPACEDIM; R++)
47 15860 : ret(i,J,k,L) += DDW(i,Q,k,R)*F0inv(J,Q)*F0inv(L,R);
48 : }
49 40 : return ret;
50 : }
51 :
52 : public:
53 : Set::Matrix F0 = Set::Matrix::Identity();
54 :
55 : public:
56 2 : static Cubic Zero()
57 : {
58 2 : Cubic ret = PseudoLinear::Cubic::Zero();
59 2 : ret.F0 = Set::Matrix::Zero();
60 2 : return ret;
61 0 : }
62 64 : static Cubic Random()
63 : {
64 64 : Cubic ret = PseudoLinear::Cubic::Random();
65 64 : ret.F0 = Set::Matrix::Random();
66 64 : return ret;
67 0 : }
68 :
69 : AMREX_FORCE_INLINE
70 : static Cubic Combine(const std::vector<Cubic> &models, const std::vector<Set::Scalar> &eta, int order)
71 : {
72 0 : Cubic ret = Cubic::Zero();
73 0 : if (order == 1)
74 : {
75 0 : Set::Scalar etasum = 0.;
76 0 : for (unsigned int n = 0 ; n < models.size(); n++) etasum += eta[n];
77 0 : for (unsigned int n = 0 ; n < models.size(); n++)
78 : {
79 0 : ret.C11 += models[n].C11 * (eta[n] / etasum);
80 0 : ret.C12 += models[n].C12 * (eta[n] / etasum);
81 0 : ret.C44 += models[n].C44 * (eta[n] / etasum);
82 0 : ret.q += models[n].q * (eta[n] / etasum);
83 0 : ret.F0 += models[n].F0 * (eta[n] / etasum);
84 : }
85 0 : return ret;
86 : }
87 0 : else if (order == 2)
88 : {
89 0 : Set::Scalar etasum = 0.;
90 0 : for (unsigned int n = 0 ; n < models.size(); n++) etasum += eta[n]*eta[n];
91 0 : for (unsigned int n = 0 ; n < models.size(); n++)
92 : {
93 0 : ret.C11 += models[n].C11 * (eta[n]*eta[n] / etasum);
94 0 : ret.C12 += models[n].C12 * (eta[n]*eta[n] / etasum);
95 0 : ret.C44 += models[n].C44 * (eta[n]*eta[n] / etasum);
96 0 : ret.q += models[n].q * (eta[n]*eta[n] / etasum);
97 0 : ret.F0 += models[n].F0 * (eta[n]*eta[n] / etasum);
98 : }
99 0 : return ret;
100 : }
101 :
102 0 : Util::Exception(INFO,"invalid value used for order, ",order);
103 0 : return ret; // should never happen
104 0 : }
105 :
106 :
107 0 : static void Parse(Cubic& value, IO::ParmParse& pp)
108 : {
109 0 : PseudoLinear::Cubic::Parse(value, pp);
110 0 : if (pp.contains("eps0") && pp.contains("F0"))
111 : {
112 0 : Util::Abort("Cannot specify both F0 and eps0");
113 : }
114 0 : else if (pp.contains("F0"))
115 : {
116 0 : Set::Matrix F0;
117 0 : pp_queryarr("F0", F0); // Large-strain eigendeformation (Identity = no deformation)
118 0 : value.F0 = F0;
119 : }
120 0 : else if (pp.contains("eps0"))
121 : {
122 0 : Set::Matrix eps0;
123 0 : pp_queryarr("eps0",eps0); // Small-strain eigendeformation (Zero = no deformation)
124 0 : value.F0 = eps0 + Set::Matrix::Identity();
125 : }
126 : else
127 : {
128 0 : value.F0 = Set::Matrix::Identity();
129 : }
130 0 : Util::Assert(INFO,TEST(fabs(value.F0.determinant()) > 1E-8 ),"F0 must be non-singular");
131 0 : }
132 :
133 : #define OP_CLASS Cubic
134 : #define OP_VARS X(C11) X(C12) X(C44) X(q) X(F0)
135 : #include "Model/Solid/InClassOperators.H"
136 : };
137 : #include "Model/Solid/ExtClassOperators.H"
138 :
139 : }
140 : }
141 : }
142 : }
143 :
144 :
145 : template<>
146 : ALAMO_SINGLE_DEFINITION
147 0 : int Set::Field<Model::Solid::Finite::PseudoAffine::Cubic>::NComp() const
148 : {
149 0 : return 4;
150 : }
151 :
152 : template<>
153 : ALAMO_SINGLE_DEFINITION
154 0 : std::string Set::Field<Model::Solid::Finite::PseudoAffine::Cubic>::Name(int i) const
155 : {
156 0 : if (i == 0) return name + "_F0xx";
157 0 : if (i == 1) return name + "_F0xy";
158 0 : if (i == 2) return name + "_F0yx";
159 0 : if (i == 3) return name + "_F0yy";
160 0 : return name;
161 : }
162 :
163 : template<>
164 : ALAMO_SINGLE_DEFINITION
165 0 : void Set::Field<Model::Solid::Finite::PseudoAffine::Cubic>::Copy(int a_lev, amrex::MultiFab& a_dst, int a_dstcomp, int a_nghost) const
166 : {
167 0 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
168 : {
169 0 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
170 0 : if (bx.ok())
171 : {
172 0 : amrex::Array4<const Model::Solid::Finite::PseudoAffine::Cubic> const& src = ((*this)[a_lev])->array(mfi);
173 0 : amrex::Array4<Set::Scalar> const& dst = a_dst.array(mfi);
174 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
175 : {
176 0 : dst(i, j, k, a_dstcomp + 0) = src(i, j, k).F0(0, 0);
177 0 : dst(i, j, k, a_dstcomp + 1) = src(i, j, k).F0(0, 1);
178 0 : dst(i, j, k, a_dstcomp + 2) = src(i, j, k).F0(1, 0);
179 0 : dst(i, j, k, a_dstcomp + 3) = src(i, j, k).F0(1, 1);
180 0 : });
181 : }
182 0 : }
183 0 : }
184 :
185 :
186 : #endif
|
