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 "Util/Util.H"
28 : #include "Base.H"
29 :
30 : namespace Set
31 : {
32 : template<>
33 : class Matrix4<AMREX_SPACEDIM,Sym::Isotropic>
34 : {
35 : Set::Scalar lambda=NAN, mu=NAN;
36 : public:
37 52221969 : AMREX_GPU_HOST_DEVICE Matrix4() {};
38 88305 : AMREX_GPU_HOST_DEVICE Matrix4(Set::Scalar a_lambda, Set::Scalar a_mu) : lambda(a_lambda), mu(a_mu) {};
39 :
40 : /// Note: for the Isotropic Matrix4 this routine works for **retrieval only**!
41 : /// If you try to assign a value using this with, say
42 : ///
43 : /// isotropicmatrix4(i,j,k,l) = 8.0
44 : ///
45 : /// you will get a `lvalue required as left operand of assignment` compile error.
46 : /// You should probably consider using a lower symmetry operator.
47 : AMREX_FORCE_INLINE
48 : Scalar operator () (const int i, const int j, const int k, const int l) const
49 : {
50 3880 : Set::Scalar ret = 0.0;
51 3880 : if (i==k && j==l) ret += mu;
52 3880 : if (i==l && j==k) ret += mu;
53 3880 : if (i==j && k==l) ret += lambda;
54 3880 : return ret;
55 : }
56 : void Randomize()
57 : {
58 : lambda = Util::Random();
59 : mu = Util::Random();
60 : }
61 0 : void Print (std::ostream& os )
62 : {
63 0 : os << "lambda = " << lambda << " mu = " << mu;
64 0 : }
65 450711 : Set::Scalar Lambda () const
66 : {
67 450711 : return lambda;
68 : }
69 450711 : Set::Scalar Mu () const
70 : {
71 450711 : return mu;
72 : }
73 562 : static Matrix4<AMREX_SPACEDIM,Sym::Isotropic> Zero()
74 : {
75 562 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> zero;
76 562 : zero.lambda = 0.0;
77 562 : zero.mu = 0.0;
78 562 : return zero;
79 : }
80 :
81 :
82 79 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> Inverse() const
83 : {
84 79 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> inv;
85 79 : inv.mu = 1./4./mu;
86 79 : inv.lambda = lambda / (2*mu*(3*lambda + 2*mu));
87 79 : return inv;
88 : }
89 : friend Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix &b);
90 : friend Set::Vector operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix3 &b);
91 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator - (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
92 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b);
93 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a);
94 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator / (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b);
95 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a);
96 : friend bool operator == (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
97 : friend Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator + (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b);
98 :
99 : //AMREX_GPU_HOST_DEVICE void operator = (Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda = a.lambda; mu = a.mu;}
100 126627 : AMREX_GPU_HOST_DEVICE void operator += (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda += a.lambda; mu += a.mu;}
101 : AMREX_GPU_HOST_DEVICE void operator -= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda -= a.lambda; mu -= a.mu;}
102 : AMREX_GPU_HOST_DEVICE void operator *= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda *= a.lambda; mu *= a.mu;}
103 : AMREX_GPU_HOST_DEVICE void operator /= (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda /= a.lambda; mu /= a.mu;}
104 : AMREX_GPU_HOST_DEVICE void operator *= (const Set::Scalar &alpha) {lambda *= alpha; mu *= alpha;}
105 : AMREX_GPU_HOST_DEVICE void operator /= (const Set::Scalar &alpha) {lambda /= alpha; mu /= alpha;}
106 :
107 : Set::Scalar Norm()
108 : {
109 : return std::sqrt(lambda*lambda + mu*mu);
110 : }
111 :
112 0 : bool contains_nan() const
113 : {
114 0 : if (std::isnan(lambda)) return true;
115 0 : if (std::isnan(mu)) return true;
116 0 : return false;
117 : }
118 : };
119 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
120 : Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix &b)
121 : {
122 50021357 : Set::Matrix ret;
123 :
124 : #if AMREX_SPACEDIM == 2
125 35046273 : ret(0,0) = (a.lambda + 2.*a.mu) * b(0,0) + a.lambda *b(1,1);
126 35046273 : ret(1,1) = a.lambda * b(0,0) + (a.lambda + 2.*a.mu)*b(1,1);
127 35046273 : ret(0,1) = a.mu*(b(0,1) + b(1,0)); ret(1,0) = ret(0,1);
128 :
129 : #elif AMREX_SPACEDIM == 3
130 14975084 : ret(0,0) = (a.lambda + 2.*a.mu) * b(0,0) + a.lambda *b(1,1) + a.lambda *b(2,2);
131 14975084 : ret(1,1) = a.lambda * b(0,0) + (a.lambda + 2.*a.mu)*b(1,1) + a.lambda *b(2,2);
132 14975084 : ret(2,2) = a.lambda * b(0,0) + a.lambda *b(1,1) + (a.lambda + 2.*a.mu)*b(2,2);
133 14975084 : ret(1,2) = a.mu*(b(1,2) + b(2,1)); ret(2,1) = ret(1,2);
134 14975084 : ret(2,0) = a.mu*(b(2,0) + b(0,2)); ret(0,2) = ret(2,0);
135 14975084 : ret(0,1) = a.mu*(b(0,1) + b(1,0)); ret(1,0) = ret(0,1);
136 :
137 : #endif
138 50021357 : return ret;
139 : }
140 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
141 : Set::Matrix operator * (const Set::Matrix &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a)
142 : {return a*b;}
143 :
144 :
145 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
146 : Set::Vector operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Matrix3 &b)
147 : {
148 11898313 : Set::Vector ret = Set::Vector::Zero();
149 :
150 37311639 : for (int i = 0; i < AMREX_SPACEDIM; i++)
151 81090078 : for (int j=0; j < AMREX_SPACEDIM; j++)
152 222707008 : ret(i) += a.mu*(b(i,j,j) + b(j,i,j)) + a.lambda*b(j,j,i);
153 :
154 11898313 : return ret;
155 : }
156 :
157 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
158 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b)
159 : {
160 25223672 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
161 25223672 : ret.mu = a.mu * b;
162 25223672 : ret.lambda = a.lambda * b;
163 25223672 : return ret;
164 : }
165 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
166 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator / (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Set::Scalar &b)
167 : {
168 25041668 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
169 25041668 : ret.mu = a.mu / b;
170 25041668 : ret.lambda = a.lambda / b;
171 25041668 : return ret;
172 : }
173 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
174 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator * (const Set::Scalar &b, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a)
175 : {
176 0 : return a*b;
177 : }
178 :
179 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
180 : bool operator == (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
181 : {
182 20 : if (a.mu != b.mu) return false;
183 20 : if (a.lambda != b.lambda) return false;
184 20 : return true;
185 : }
186 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
187 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator + (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a, const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
188 : {
189 70959 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret;// = Set::Vector::Zero();
190 70959 : ret.mu = a.mu + b.mu;
191 70959 : ret.lambda = a.lambda + b.lambda;
192 70959 : return ret;
193 : }
194 :
195 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
196 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> operator - (const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a,
197 : const Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &b)
198 : {
199 25042830 : Matrix4<AMREX_SPACEDIM,Sym::Isotropic> ret = a;
200 25042830 : ret.mu -= b.mu;
201 25042830 : ret.lambda -= b.lambda;
202 25042830 : return ret;
203 : }
204 :
205 :
206 : }
207 : #endif
|
