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