Line data Source code
1 : //
2 : // The isotropic tensor is defined to be
3 : //
4 : // .. math::
5 : //
6 : // C_{ijkl} = \mu(d_{il}d_{jk} + d_{ik}d_{jl}) + \lambda d_{ij} d_{kl}
7 : //
8 : // The inverse ("compliance") tensor is
9 : //
10 : // .. math::
11 : //
12 : // S_{ijkl} = ((1+\nu)/2E)(d_{il}d_{jl} + d_{ik}d_{jl}) - (\nu/E)d_{ij}d_{kl}
13 : //
14 : // Replacing E, nu with Lame constants gives:
15 : //
16 : // .. math::
17 : //
18 : // S_{ijkl} = (1/(4 \mu))(d_{il}d_{jl} + d_{ik}d_{jl}) + (\lambda/(2*\mu*(3\lambda+2\mu))) * d_{ij} * d_{kl}
19 : //
20 : // For reference: http://solidmechanics.org/text/Chapter3_2/Chapter3_2.htm
21 : // https://en.wikipedia.org/wiki/Lam%C3%A9_parameters
22 : //
23 :
24 : #ifndef SET_MATRIX4_ISOTROPIC_H
25 : #define SET_MATRIX4_ISOTROPIC_H
26 :
27 : #include "AMReX_Config.H"
28 : #include "Util/Util.H"
29 : #include "Set/Set.H"
30 : #include "Set/Base.H"
31 : #include "Matrix3.H"
32 : #include "Matrix4.H"
33 :
34 : namespace Set
35 : {
36 : template<>
37 : class Matrix4<AMREX_SPACEDIM,Sym::Isotropic>
38 : {
39 : Set::Scalar lambda=NAN, mu=NAN;
40 : public:
41 52222959 : AMREX_GPU_HOST_DEVICE Matrix4() {};
42 88301 : AMREX_GPU_HOST_DEVICE Matrix4(Set::Scalar a_lambda, Set::Scalar a_mu) : lambda(a_lambda), mu(a_mu) {};
43 :
44 : /// Note: for the Isotropic Matrix4 this routine works for **retrieval only**!
45 : /// If you try to assign a value using this with, say
46 : ///
47 : /// isotropicmatrix4(i,j,k,l) = 8.0
48 : ///
49 : /// you will get a `lvalue required as left operand of assignment` compile error.
50 : /// You should probably consider using a lower symmetry operator.
51 : AMREX_FORCE_INLINE
52 : Scalar operator () (const int i, const int j, const int k, const int l) const
53 : {
54 3880 : Set::Scalar ret = 0.0;
55 3880 : if (i==k && j==l) ret += mu;
56 3880 : if (i==l && j==k) ret += mu;
57 3880 : if (i==j && k==l) ret += lambda;
58 3880 : return ret;
59 : }
60 2 : void Randomize()
61 : {
62 2 : lambda = Util::Random();
63 2 : mu = Util::Random();
64 2 : }
65 2 : static Matrix4<AMREX_SPACEDIM,Sym::Isotropic> Random()
66 : {
67 2 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;
68 2 : ret.Randomize();
69 2 : return ret;
70 : }
71 : void Print (std::ostream& os )
72 : {
73 : os << "lambda = " << lambda << " mu = " << mu;
74 : }
75 1071238 : Set::Scalar Lambda () const
76 : {
77 1071238 : return lambda;
78 : }
79 1071238 : Set::Scalar Mu () const
80 : {
81 1071238 : return mu;
82 : }
83 : Set::Scalar Youngs() const
84 : {
85 : return mu * ( 3.0*lambda + 2.0*mu ) / (lambda + mu);
86 : }
87 : Set::Scalar Nu() const
88 : {
89 : return 0.5 * lambda / (lambda + mu);
90 : }
91 562 : static Matrix4<AMREX_SPACEDIM,Sym::Isotropic> Zero()
92 : {
93 562 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> zero;
94 562 : zero.lambda = 0.0;
95 562 : zero.mu = 0.0;
96 562 : return zero;
97 : }
98 :
99 :
100 79 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> Inverse() const
101 : {
102 79 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> inv;
103 79 : inv.mu = 1./4./mu;
104 79 : inv.lambda = lambda / (2*mu*(3*lambda + 2*mu));
105 79 : return inv;
106 : }
107 : friend Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix &b);
108 : friend Set::Vector operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix3 &b);
109 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator - (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
110 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b);
111 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a);
112 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator / (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b);
113 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a);
114 : friend bool operator == (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
115 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator + (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
116 :
117 : //AMREX_GPU_HOST_DEVICE void operator = (Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda = a.lambda; mu = a.mu;}
118 126627 : AMREX_GPU_HOST_DEVICE void operator += (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda += a.lambda; mu += a.mu;}
119 : AMREX_GPU_HOST_DEVICE void operator -= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda -= a.lambda; mu -= a.mu;}
120 : AMREX_GPU_HOST_DEVICE void operator *= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda *= a.lambda; mu *= a.mu;}
121 : AMREX_GPU_HOST_DEVICE void operator /= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda /= a.lambda; mu /= a.mu;}
122 : AMREX_GPU_HOST_DEVICE void operator *= (const Set::Scalar &alpha) {lambda *= alpha; mu *= alpha;}
123 : AMREX_GPU_HOST_DEVICE void operator /= (const Set::Scalar &alpha) {lambda /= alpha; mu /= alpha;}
124 :
125 : Set::Scalar Norm()
126 : {
127 : return std::sqrt(lambda*lambda + mu*mu);
128 : }
129 :
130 : bool contains_nan() const
131 : {
132 : if (std::isnan(lambda)) return true;
133 : if (std::isnan(mu)) return true;
134 : return false;
135 : }
136 : };
137 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
138 : Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix &b)
139 : {
140 50024085 : Set::Matrix ret;
141 :
142 : #if AMREX_SPACEDIM == 2
143 35049001 : ret(0,0) = (a.lambda + 2.*a.mu) * b(0,0) + a.lambda *b(1,1);
144 35049001 : ret(1,1) = a.lambda * b(0,0) + (a.lambda + 2.*a.mu)*b(1,1);
145 35049001 : ret(0,1) = a.mu*(b(0,1) + b(1,0)); ret(1,0) = ret(0,1);
146 :
147 : #elif AMREX_SPACEDIM == 3
148 14975084 : ret(0,0) = (a.lambda + 2.*a.mu) * b(0,0) + a.lambda *b(1,1) + a.lambda *b(2,2);
149 14975084 : ret(1,1) = a.lambda * b(0,0) + (a.lambda + 2.*a.mu)*b(1,1) + a.lambda *b(2,2);
150 14975084 : ret(2,2) = a.lambda * b(0,0) + a.lambda *b(1,1) + (a.lambda + 2.*a.mu)*b(2,2);
151 14975084 : ret(1,2) = a.mu*(b(1,2) + b(2,1)); ret(2,1) = ret(1,2);
152 14975084 : ret(2,0) = a.mu*(b(2,0) + b(0,2)); ret(0,2) = ret(2,0);
153 14975084 : ret(0,1) = a.mu*(b(0,1) + b(1,0)); ret(1,0) = ret(0,1);
154 :
155 : #endif
156 50024085 : return ret;
157 : }
158 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
159 : Set::Matrix operator * (const Set::Matrix &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a)
160 : {return a*b;}
161 :
162 :
163 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
164 : Set::Vector operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix3 &b)
165 : {
166 11898561 : Set::Vector ret = Set::Vector::Zero();
167 :
168 37312383 : for (int i = 0; i < AMREX_SPACEDIM; i++)
169 81091566 : for (int j=0; j < AMREX_SPACEDIM; j++)
170 222710976 : ret(i) += a.mu*(b(i,j,j) + b(j,i,j)) + a.lambda*b(j,j,i);
171 :
172 11898561 : return ret;
173 : }
174 :
175 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
176 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b)
177 : {
178 25224168 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
179 25224168 : ret.mu = a.mu * b;
180 25224168 : ret.lambda = a.lambda * b;
181 25224168 : return ret;
182 : }
183 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
184 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator / (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b)
185 : {
186 25042164 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
187 25042164 : ret.mu = a.mu / b;
188 25042164 : ret.lambda = a.lambda / b;
189 25042164 : return ret;
190 : }
191 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
192 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a)
193 : {
194 0 : return a*b;
195 : }
196 :
197 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
198 : bool operator == (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
199 : {
200 20 : if (a.mu != b.mu) return false;
201 20 : if (a.lambda != b.lambda) return false;
202 20 : return true;
203 : }
204 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
205 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator + (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
206 : {
207 70959 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
208 70959 : ret.mu = a.mu + b.mu;
209 70959 : ret.lambda = a.lambda + b.lambda;
210 70959 : return ret;
211 : }
212 :
213 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
214 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator - (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,
215 : const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
216 : {
217 25043326 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret = a;
218 25043326 : ret.mu -= b.mu;
219 25043326 : ret.lambda -= b.lambda;
220 25043326 : return ret;
221 : }
222 :
223 :
224 : }
225 : #endif
|
