OOFEM 3.0
Loading...
Searching...
No Matches
lsmastermat.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "lsmastermat.h"
36#include "sm/Materials/isolinearelasticmaterial.h"
37#include "gausspoint.h"
38#include "floatmatrix.h"
39#include "floatarray.h"
40#include "floatmatrixf.h"
41#include "floatarrayf.h"
42#include "intarray.h"
43#include "stressvector.h"
44#include "strainvector.h"
45#include "mathfem.h"
46#include "contextioerr.h"
47#include "datastream.h"
48#include "classfactory.h"
49
50namespace oofem {
51REGISTER_Material(LargeStrainMasterMaterial);
52
53LargeStrainMasterMaterial :: LargeStrainMasterMaterial(int n, Domain *d) : StructuralMaterial(n, d)
54{
55}
56
57void
58LargeStrainMasterMaterial :: initializeFrom(InputRecord &ir)
59{
60 IR_GIVE_OPTIONAL_FIELD(ir, slaveMat, _IFT_LargeStrainMasterMaterial_slaveMat); // number of slave material
61 IR_GIVE_OPTIONAL_FIELD(ir, m, _IFT_LargeStrainMasterMaterial_m); // type of Set-Hill strain tensor
62}
63
64std::unique_ptr<MaterialStatus>
65LargeStrainMasterMaterial :: CreateStatus(GaussPoint *gp) const
66{
67 return std::make_unique<LargeStrainMasterMaterialStatus>(gp, domain, slaveMat);
68}
69
70
71FloatArrayF<9>
72LargeStrainMasterMaterial :: giveFirstPKStressVector_3d(const FloatArrayF<9> &vF, GaussPoint *gp, TimeStep *tStep) const
73{
74 auto status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
75 this->initTempStatus(gp);
76
77 auto sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
78
79 //store of deformation gradient into 3x3 matrix
80 auto F = from_voigt_form_9(vF);
81 //compute right Cauchy-Green tensor(C), its eigenvalues and eigenvectors
82 auto C = Tdot(F, F);
83 // compute eigen values and eigen vectors of C
84 FloatArray eVals;
85 FloatMatrix eVecs;
86 FloatMatrix(C).jaco_(eVals, eVecs, 15);
87
88 // compute Seth - Hill's strain measure, it depends on mParameter
89 double lambda1 = eVals.at(1);
90 double lambda2 = eVals.at(2);
91 double lambda3 = eVals.at(3);
92 double E1, E2, E3;
93 if ( m == 0 ) {
94 E1 = 1. / 2. * log(lambda1);
95 E2 = 1. / 2. * log(lambda2);
96 E3 = 1. / 2. * log(lambda3);
97 } else {
98 E1 = 1. / ( 2. * m ) * ( pow(lambda1, m) - 1. );
99 E2 = 1. / ( 2. * m ) * ( pow(lambda2, m) - 1. );
100 E3 = 1. / ( 2. * m ) * ( pow(lambda3, m) - 1. );
101 }
102
103 FloatMatrixF<3,3> SethHillStrain;
104 for ( int i = 1; i < 4; i++ ) {
105 for ( int j = 1; j < 4; j++ ) {
106 SethHillStrain.at(i, j) = E1 * eVecs.at(i, 1) * eVecs.at(j, 1) + E2 *eVecs.at(i, 2) * eVecs.at(j, 2) + E3 *eVecs.at(i, 3) * eVecs.at(j, 3);
107 }
108 }
109
110 auto SethHillStrainVector = to_voigt_strain_33(SethHillStrain);
111 auto stressVector = sMat->giveRealStressVector_3d(SethHillStrainVector, gp, tStep);
112
113 auto T = this->constructTransformationMatrix(eVecs);
114 stressVector.at(4) = 2;
115 stressVector.at(5) = 2;
116 stressVector.at(6) = 2;
117
118 auto stressM = dot(T, FloatArrayF<6>(stressVector));
119 stressM.at(4) *= 1. / 2.;
120 stressM.at(5) *= 1. / 2.;
121 stressM.at(6) *= 1. / 2.;
122
123 //auto [L1, L2] = this->constructL1L2TransformationMatrices(eVals, stressM, E1, E2, E3); // c++17
124 auto tmp = this->constructL1L2TransformationMatrices(eVals, stressM, E1, E2, E3);
125 auto L1 = tmp.first;
126 auto L2 = tmp.second;
127
128 auto P = Tdot(T, dot(L1, T));
129
130 //transformation of the stress to the 2PK stress and then to 1PK
131 stressVector.at(4) *= 0.5;
132 stressVector.at(5) *= 0.5;
133 stressVector.at(6) *= 0.5;
134 auto secondPK = dot(P, FloatArrayF<6>(stressVector));
135
136 auto firstPK = dot(F, from_voigt_stress_6(secondPK)); // P = F*S
137 auto firstPKv = to_voigt_form_33(firstPK);
138 auto TL = Tdot(T, dot(L2, T));
139 status->setPmatrix(P);
140 status->setTLmatrix(TL);
141 status->letTempStressVectorBe(firstPKv);
142 return firstPKv;
143}
144
145
146FloatMatrixF<6,6>
147LargeStrainMasterMaterial :: constructTransformationMatrix(const FloatMatrixF<3,3> &eVecs) const
148{
149 FloatMatrixF<6,6> answer;
150 answer.at(1, 1) = eVecs.at(1, 1) * eVecs.at(1, 1);
151 answer.at(1, 2) = eVecs.at(2, 1) * eVecs.at(2, 1);
152 answer.at(1, 3) = eVecs.at(3, 1) * eVecs.at(3, 1);
153 answer.at(1, 4) = eVecs.at(2, 1) * eVecs.at(3, 1);
154 answer.at(1, 5) = eVecs.at(1, 1) * eVecs.at(3, 1);
155 answer.at(1, 6) = eVecs.at(1, 1) * eVecs.at(2, 1);
156
157 answer.at(2, 1) = eVecs.at(1, 2) * eVecs.at(1, 2);
158 answer.at(2, 2) = eVecs.at(2, 2) * eVecs.at(2, 2);
159 answer.at(2, 3) = eVecs.at(3, 2) * eVecs.at(3, 2);
160 answer.at(2, 4) = eVecs.at(2, 2) * eVecs.at(3, 2);
161 answer.at(2, 5) = eVecs.at(1, 2) * eVecs.at(3, 2);
162 answer.at(2, 6) = eVecs.at(1, 2) * eVecs.at(2, 2);
163
164 answer.at(3, 1) = eVecs.at(1, 3) * eVecs.at(1, 3);
165 answer.at(3, 2) = eVecs.at(2, 3) * eVecs.at(2, 3);
166 answer.at(3, 3) = eVecs.at(3, 3) * eVecs.at(3, 3);
167 answer.at(3, 4) = eVecs.at(2, 3) * eVecs.at(3, 3);
168 answer.at(3, 5) = eVecs.at(1, 3) * eVecs.at(3, 3);
169 answer.at(3, 6) = eVecs.at(1, 3) * eVecs.at(2, 3);
170
171 answer.at(4, 1) = 2 * eVecs.at(1, 2) * eVecs.at(1, 3);
172 answer.at(4, 2) = 2 * eVecs.at(2, 2) * eVecs.at(2, 3);
173 answer.at(4, 3) = 2 * eVecs.at(3, 2) * eVecs.at(3, 3);
174 answer.at(4, 4) = eVecs.at(2, 2) * eVecs.at(3, 3) + eVecs.at(3, 2) * eVecs.at(2, 3);
175 answer.at(4, 5) = eVecs.at(1, 2) * eVecs.at(3, 3) + eVecs.at(3, 2) * eVecs.at(1, 3);
176 answer.at(4, 6) = eVecs.at(1, 2) * eVecs.at(2, 3) + eVecs.at(2, 2) * eVecs.at(1, 3);
177
178 answer.at(5, 1) = 2 * eVecs.at(1, 1) * eVecs.at(1, 3);
179 answer.at(5, 2) = 2 * eVecs.at(2, 1) * eVecs.at(2, 3);
180 answer.at(5, 3) = 2 * eVecs.at(3, 1) * eVecs.at(3, 3);
181 answer.at(5, 4) = eVecs.at(2, 1) * eVecs.at(3, 3) + eVecs.at(3, 1) * eVecs.at(2, 3);
182 answer.at(5, 5) = eVecs.at(1, 1) * eVecs.at(3, 3) + eVecs.at(3, 1) * eVecs.at(1, 3);
183 answer.at(5, 6) = eVecs.at(1, 1) * eVecs.at(2, 3) + eVecs.at(2, 1) * eVecs.at(1, 3);
184
185 answer.at(6, 1) = 2 * eVecs.at(1, 1) * eVecs.at(1, 2);
186 answer.at(6, 2) = 2 * eVecs.at(2, 1) * eVecs.at(2, 2);
187 answer.at(6, 3) = 2 * eVecs.at(3, 1) * eVecs.at(3, 2);
188 answer.at(6, 4) = eVecs.at(2, 1) * eVecs.at(3, 2) + eVecs.at(3, 1) * eVecs.at(2, 2);
189 answer.at(6, 5) = eVecs.at(1, 1) * eVecs.at(3, 2) + eVecs.at(3, 1) * eVecs.at(1, 2);
190 answer.at(6, 6) = eVecs.at(1, 1) * eVecs.at(2, 2) + eVecs.at(2, 1) * eVecs.at(1, 2);
191 return answer;
192}
193
194
195std::pair<FloatMatrixF<6,6>, FloatMatrixF<6,6>>
196LargeStrainMasterMaterial :: constructL1L2TransformationMatrices(const FloatArrayF<3> &eigenValues, const FloatArrayF<6> &stressM, double E1, double E2, double E3) const
197{
198 FloatMatrixF<6,6> answer1, answer2;
199 double gamma12, gamma13, gamma23, gamma112, gamma221, gamma113, gamma331, gamma223, gamma332, gamma;
200 double lambda1 = eigenValues.at(1);
201 double lambda2 = eigenValues.at(2);
202 double lambda3 = eigenValues.at(3);
203 double lambda1P = pow(lambda1, m - 1);
204 double lambda2P = pow(lambda2, m - 1);
205 double lambda3P = pow(lambda3, m - 1);
206 if ( ( lambda1 == lambda2 ) && ( lambda2 == lambda3 ) ) { // three equal eigenvalues
207 gamma12 = gamma13 = gamma23 = 1. / 2. * lambda1P;
208
209 answer2.at(1, 1) = 2 * stressM.at(1) * ( m - 1 ) * pow(lambda1, m - 2);
210 answer2.at(2, 2) = 2 * stressM.at(2) * ( m - 1 ) * pow(lambda2, m - 2);
211 answer2.at(3, 3) = 2 * stressM.at(3) * ( m - 1 ) * pow(lambda3, m - 2);
212 answer2.at(4, 4) = 1. / 2. * ( stressM.at(2) + stressM.at(3) ) * ( m - 1 ) * pow(lambda2, m - 2);
213 answer2.at(5, 5) = 1. / 2. * ( stressM.at(1) + stressM.at(3) ) * ( m - 1 ) * pow(lambda3, m - 2);
214 answer2.at(6, 6) = 1. / 2. * ( stressM.at(1) + stressM.at(2) ) * ( m - 1 ) * pow(lambda1, m - 2);
215
216 answer2.at(1, 5) = answer2.at(5, 1) = stressM.at(5) * ( m - 1 ) * pow(lambda3, m - 2);
217 answer2.at(1, 6) = answer2.at(6, 1) = stressM.at(6) * ( m - 1 ) * pow(lambda1, m - 2);
218 answer2.at(2, 4) = answer2.at(4, 2) = stressM.at(4) * ( m - 1 ) * pow(lambda2, m - 2);
219 answer2.at(2, 6) = answer2.at(6, 2) = stressM.at(6) * ( m - 1 ) * pow(lambda1, m - 2);
220 answer2.at(3, 4) = answer2.at(4, 3) = stressM.at(4) * ( m - 1 ) * pow(lambda2, m - 2);
221 answer2.at(3, 5) = answer2.at(5, 3) = stressM.at(5) * ( m - 1 ) * pow(lambda3, m - 2);
222
223 answer2.at(4, 5) = answer2.at(5, 4) = 1. / 2. * stressM.at(6) * ( m - 1 ) * pow(lambda1, m - 2);
224 answer2.at(4, 6) = answer2.at(6, 4) = 1. / 2. * stressM.at(5) * ( m - 1 ) * pow(lambda1, m - 2);
225 answer2.at(5, 6) = answer2.at(6, 5) = 1. / 2. * stressM.at(4) * ( m - 1 ) * pow(lambda1, m - 2);
226 } else if ( lambda1 == lambda2 ) { //two equal eigenvalues
227 gamma12 = 1. / 2. * lambda1P;
228 gamma13 = ( E1 - E3 ) / ( lambda1 - lambda3 );
229 gamma23 = ( E2 - E3 ) / ( lambda2 - lambda3 );
230 gamma113 = ( pow(lambda1, m - 1) * ( lambda1 - lambda3 ) - 2 * ( E1 - E3 ) ) / ( ( lambda1 - lambda3 ) * ( lambda1 - lambda3 ) );
231 gamma331 = ( pow(lambda3, m - 1) * ( lambda3 - lambda1 ) - 2 * ( E3 - E1 ) ) / ( ( lambda3 - lambda1 ) * ( lambda3 - lambda1 ) );
232
233
234 answer2.at(1, 1) = 2 * stressM.at(1) * ( m - 1 ) * pow(lambda1, m - 2);
235 answer2.at(2, 2) = 2 * stressM.at(2) * ( m - 1 ) * pow(lambda2, m - 2);
236 answer2.at(3, 3) = 2 * stressM.at(3) * ( m - 1 ) * pow(lambda3, m - 2);
237
238 answer2.at(4, 4) = stressM.at(2) * gamma113 + stressM.at(3) * gamma331;
239 answer2.at(5, 5) = stressM.at(1) * gamma113 + stressM.at(3) * gamma331;
240 answer2.at(6, 6) = 1. / 2. * ( stressM.at(1) + stressM.at(2) ) * ( m - 1 ) * pow(lambda1, m - 2);
241
242 answer2.at(1, 5) = answer2.at(5, 1) = 2. * stressM.at(5) * gamma113;
243 answer2.at(1, 6) = answer2.at(6, 1) = stressM.at(6) * ( m - 1 ) * pow(lambda1, m - 2);
244 answer2.at(2, 4) = answer2.at(4, 2) = 2. * stressM.at(4) * gamma113;
245 answer2.at(2, 6) = answer2.at(6, 2) = stressM.at(6) * ( m - 1 ) * pow(lambda1, m - 2);
246 answer2.at(3, 4) = answer2.at(4, 3) = 2. * stressM.at(4) * gamma331;
247 answer2.at(3, 5) = answer2.at(5, 3) = 2. * stressM.at(5) * gamma331;
248 answer2.at(4, 5) = answer2.at(5, 4) = stressM.at(6) * gamma113;
249 answer2.at(4, 6) = answer2.at(6, 4) = stressM.at(5) * gamma113;
250 answer2.at(5, 6) = answer2.at(6, 5) = stressM.at(4) * gamma113;
251 } else if ( lambda2 == lambda3 ) {
252 gamma23 = 1. / 2. * lambda2P;
253 gamma12 = ( E1 - E2 ) / ( lambda1 - lambda2 );
254 gamma13 = ( E1 - E3 ) / ( lambda1 - lambda3 );
255 gamma112 = ( pow(lambda1, m - 1) * ( lambda1 - lambda2 ) - 2 * ( E1 - E2 ) ) / ( ( lambda1 - lambda2 ) * ( lambda1 - lambda2 ) );
256 gamma221 = ( pow(lambda2, m - 1) * ( lambda2 - lambda1 ) - 2 * ( E2 - E1 ) ) / ( ( lambda2 - lambda1 ) * ( lambda2 - lambda1 ) );
257
258 answer2.at(1, 1) = 2 * stressM.at(1) * ( m - 1 ) * pow(lambda1, m - 2);
259 answer2.at(2, 2) = 2 * stressM.at(2) * ( m - 1 ) * pow(lambda2, m - 2);
260 answer2.at(3, 3) = 2 * stressM.at(3) * ( m - 1 ) * pow(lambda3, m - 2);
261
262 answer2.at(4, 4) = 1. / 2. * ( stressM.at(2) + stressM.at(3) ) * ( m - 1 ) * pow(lambda2, m - 2);
263 answer2.at(5, 5) = stressM.at(1) * gamma112 + stressM.at(3) * gamma221;
264 answer2.at(6, 6) = stressM.at(1) * gamma112 + stressM.at(2) * gamma221;
265
266 answer2.at(1, 5) = answer2.at(5, 1) = 2. * stressM.at(5) * gamma112;
267 answer2.at(1, 6) = answer2.at(6, 1) = 2. * stressM.at(6) * gamma112;
268 answer2.at(2, 4) = answer2.at(4, 2) = stressM.at(4) * ( m - 1 ) * pow(lambda2, m - 2);
269 answer2.at(2, 6) = answer2.at(6, 2) = 2. * stressM.at(6) * gamma221;
270 answer2.at(3, 4) = answer2.at(4, 3) = stressM.at(4) * ( m - 1 ) * pow(lambda2, m - 2);
271 answer2.at(3, 5) = answer2.at(5, 3) = 2. * stressM.at(5) * gamma221;
272 answer2.at(4, 5) = answer2.at(5, 4) = stressM.at(6) * gamma221;
273 answer2.at(4, 6) = answer2.at(6, 4) = stressM.at(5) * gamma221;
274 answer2.at(5, 6) = answer2.at(6, 5) = stressM.at(4) * gamma221;
275 } else if ( lambda1 == lambda3 ) {
276 gamma13 = 1. / 2. * lambda1P;
277 gamma12 = ( E1 - E2 ) / ( lambda1 - lambda2 );
278 gamma23 = ( E2 - E3 ) / ( lambda2 - lambda3 );
279 gamma223 = ( pow(lambda2, m - 1) * ( lambda2 - lambda3 ) - 2 * ( E2 - E3 ) ) / ( ( lambda2 - lambda3 ) * ( lambda2 - lambda3 ) );
280 gamma332 = ( pow(lambda3, m - 1) * ( lambda3 - lambda2 ) - 2 * ( E3 - E2 ) ) / ( ( lambda3 - lambda2 ) * ( lambda3 - lambda2 ) );
281
282 answer2.at(1, 1) = 2 * stressM.at(1) * ( m - 1 ) * pow(lambda1, m - 2);
283 answer2.at(2, 2) = 2 * stressM.at(2) * ( m - 1 ) * pow(lambda2, m - 2);
284 answer2.at(3, 3) = 2 * stressM.at(3) * ( m - 1 ) * pow(lambda3, m - 2);
285
286 answer2.at(4, 4) = stressM.at(2) * gamma223 + stressM.at(3) * gamma332;
287 answer2.at(5, 5) = 1. / 2. * ( stressM.at(1) + stressM.at(3) ) * ( m - 1 ) * pow(lambda3, m - 2);
288 answer2.at(6, 6) = stressM.at(1) * gamma332 + stressM.at(2) * gamma223;
289
290 answer2.at(1, 5) = answer2.at(5, 1) = stressM.at(5) * ( m - 1 ) * pow(lambda3, m - 2);
291 answer2.at(1, 6) = answer2.at(6, 1) = 2. * stressM.at(6) * gamma332;
292 answer2.at(2, 4) = answer2.at(4, 2) = 2. * stressM.at(4) * gamma223;
293 answer2.at(2, 6) = answer2.at(6, 2) = 2. * stressM.at(4) * gamma223;
294 answer2.at(3, 4) = answer2.at(4, 3) = 2. * stressM.at(4) * gamma332;
295 answer2.at(3, 5) = answer2.at(5, 3) = stressM.at(5) * ( m - 1 ) * pow(lambda3, m - 2);
296 answer2.at(4, 5) = answer2.at(5, 4) = stressM.at(6) * gamma332;
297 answer2.at(4, 6) = answer2.at(6, 4) = stressM.at(5) * gamma332;
298 answer2.at(5, 6) = answer2.at(6, 5) = stressM.at(4) * gamma332;
299 } else { //three different eigenvalues
300 gamma12 = ( E1 - E2 ) / ( lambda1 - lambda2 );
301 gamma13 = ( E1 - E3 ) / ( lambda1 - lambda3 );
302 gamma23 = ( E2 - E3 ) / ( lambda2 - lambda3 );
303
304 gamma112 = ( pow(lambda1, m - 1) * ( lambda1 - lambda2 ) - 2 * ( E1 - E2 ) ) / ( ( lambda1 - lambda2 ) * ( lambda1 - lambda2 ) );
305 gamma221 = ( pow(lambda2, m - 1) * ( lambda2 - lambda1 ) - 2 * ( E2 - E1 ) ) / ( ( lambda2 - lambda1 ) * ( lambda2 - lambda1 ) );
306 gamma113 = ( pow(lambda1, m - 1) * ( lambda1 - lambda3 ) - 2 * ( E1 - E3 ) ) / ( ( lambda1 - lambda3 ) * ( lambda1 - lambda3 ) );
307 gamma331 = ( pow(lambda3, m - 1) * ( lambda3 - lambda1 ) - 2 * ( E3 - E1 ) ) / ( ( lambda3 - lambda1 ) * ( lambda3 - lambda1 ) );
308 gamma223 = ( pow(lambda2, m - 1) * ( lambda2 - lambda3 ) - 2 * ( E2 - E3 ) ) / ( ( lambda2 - lambda3 ) * ( lambda2 - lambda3 ) );
309 gamma332 = ( pow(lambda3, m - 1) * ( lambda3 - lambda2 ) - 2 * ( E3 - E2 ) ) / ( ( lambda3 - lambda2 ) * ( lambda3 - lambda2 ) );
310
311 gamma = ( lambda1 * ( E2 - E3 ) + lambda2 * ( E3 - E1 ) + lambda3 * ( E1 - E2 ) ) / ( ( lambda1 - lambda2 ) * ( lambda2 - lambda3 ) * ( lambda3 - lambda1 ) );
312
313 answer2.at(1, 1) = 2 * stressM.at(1) * ( m - 1 ) * pow(lambda1, m - 2);
314 answer2.at(2, 2) = 2 * stressM.at(2) * ( m - 1 ) * pow(lambda2, m - 2);
315 answer2.at(3, 3) = 2 * stressM.at(3) * ( m - 1 ) * pow(lambda3, m - 2);
316
317 answer2.at(4, 4) = stressM.at(2) * gamma223 + stressM.at(3) * gamma332;
318 answer2.at(5, 5) = stressM.at(1) * gamma113 + stressM.at(3) * gamma331;
319 answer2.at(6, 6) = stressM.at(1) * gamma112 + stressM.at(2) * gamma221;
320
321 answer2.at(1, 5) = answer2.at(5, 1) = 2. * stressM.at(5) * gamma113;
322 answer2.at(1, 6) = answer2.at(6, 1) = 2. * stressM.at(6) * gamma112;
323 answer2.at(2, 4) = answer2.at(4, 2) = 2. * stressM.at(4) * gamma223;
324 answer2.at(2, 6) = answer2.at(6, 2) = 2. * stressM.at(6) * gamma221;
325 answer2.at(3, 4) = answer2.at(4, 3) = 2. * stressM.at(4) * gamma332;
326 answer2.at(3, 5) = answer2.at(5, 3) = 2. * stressM.at(5) * gamma331;
327 answer2.at(4, 5) = answer2.at(5, 4) = 2. * stressM.at(6) * gamma;
328 answer2.at(4, 6) = answer2.at(6, 4) = 2. * stressM.at(5) * gamma;
329 answer2.at(5, 6) = answer2.at(6, 5) = 2. * stressM.at(4) * gamma;
330 }
331
332
333 answer1.at(1, 1) = lambda1P;
334 answer1.at(2, 2) = lambda2P;
335 answer1.at(3, 3) = lambda3P;
336 answer1.at(4, 4) = gamma23;
337 answer1.at(5, 5) = gamma13;
338 answer1.at(6, 6) = gamma12;
339 return {answer1, answer2};
340}
341
342FloatMatrixF<9,9>
343LargeStrainMasterMaterial :: give3dMaterialStiffnessMatrix_dPdF(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const
344{
345 auto status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
346 auto sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
347
348 auto stiffness = sMat->give3dMaterialStiffnessMatrix(mode, gp, tStep);
350 stiffness.at(1, 4) = 2. * stiffness.at(1, 4);
351 stiffness.at(4, 1) = 2. * stiffness.at(4, 1);
352 stiffness.at(1, 5) = 2. * stiffness.at(1, 5);
353 stiffness.at(5, 1) = 2. * stiffness.at(5, 1);
354 stiffness.at(1, 6) = 2. * stiffness.at(1, 6);
355 stiffness.at(6, 1) = 2. * stiffness.at(6, 1);
356 stiffness.at(2, 4) = 2. * stiffness.at(2, 4);
357 stiffness.at(4, 2) = 2. * stiffness.at(4, 2);
358 stiffness.at(2, 5) = 2. * stiffness.at(2, 5);
359 stiffness.at(5, 2) = 2. * stiffness.at(5, 2);
360 stiffness.at(2, 6) = 2. * stiffness.at(2, 6);
361 stiffness.at(6, 2) = 2. * stiffness.at(6, 2);
362 stiffness.at(3, 4) = 2. * stiffness.at(3, 4);
363 stiffness.at(4, 3) = 2. * stiffness.at(4, 3);
364 stiffness.at(3, 5) = 2. * stiffness.at(3, 5);
365 stiffness.at(5, 3) = 2. * stiffness.at(5, 3);
366 stiffness.at(3, 6) = 2. * stiffness.at(3, 6);
367 stiffness.at(6, 3) = 2. * stiffness.at(6, 3);
368 stiffness.at(4, 4) = 4. * stiffness.at(4, 4);
369 stiffness.at(4, 5) = 4. * stiffness.at(4, 5);
370 stiffness.at(5, 4) = 4. * stiffness.at(5, 4);
371 stiffness.at(4, 6) = 4. * stiffness.at(4, 6);
372 stiffness.at(6, 4) = 4. * stiffness.at(6, 4);
373 stiffness.at(5, 5) = 4. * stiffness.at(5, 5);
374 stiffness.at(5, 6) = 4. * stiffness.at(5, 6);
375 stiffness.at(6, 5) = 4. * stiffness.at(6, 5);
376 stiffness.at(6, 6) = 4. * stiffness.at(6, 6);
378 auto junk = dot(FloatMatrixF<6,6>(stiffness), status->givePmatrix());
379 auto stiffness2 = dot(status->givePmatrix(), junk) + status->giveTLmatrix();
380
381 return convert_dSdE_2_dPdF_3D(stiffness2, status->giveTempStressVector(), status->giveTempFVector());
382}
383
384
385int
386LargeStrainMasterMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
387{
388 auto status = static_cast< LargeStrainMasterMaterialStatus * >( this->giveStatus(gp) );
389
390 if ( type == IST_StressTensor ) {
391 answer = status->giveStressVector();
392 return 1;
393 } else if ( type == IST_StrainTensor ) {
394 answer = status->giveStrainVector();
395 return 1;
396 } else {
397 StructuralMaterial *sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
398 return sMat->giveIPValue(answer, gp, type, tStep);
399 }
400}
401
402
403//=============================================================================
404
405LargeStrainMasterMaterialStatus :: LargeStrainMasterMaterialStatus(GaussPoint *g, Domain *d, int s) :
406 StructuralMaterialStatus(g),
407 domain(d),
408 slaveMat(s)
409{
410 Pmatrix = eye<6>();
411}
412
413
414void
415LargeStrainMasterMaterialStatus :: printOutputAt(FILE *file, TimeStep *tStep) const
416{
417 auto sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
418 auto mS = sMat->giveStatus(gp);
419
420 mS->printOutputAt(file, tStep);
421 // StructuralMaterialStatus :: printOutputAt(file, tStep);
422}
423
424
425// initializes temporary variables based on their values at the previous equlibrium state
426void LargeStrainMasterMaterialStatus :: initTempStatus()
427{
428 auto sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
429 auto mS = sMat->giveStatus(gp);
430 mS->initTempStatus();
431 //StructuralMaterialStatus :: initTempStatus();
432}
433
434
435// updates internal variables when equilibrium is reached
436void
437LargeStrainMasterMaterialStatus :: updateYourself(TimeStep *tStep)
438{
439 auto sMat = static_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
440 auto mS = sMat->giveStatus(gp);
441 mS->updateYourself(tStep);
442 // StructuralMaterialStatus :: updateYourself(tStep);
443}
444
445
446void
447LargeStrainMasterMaterialStatus :: saveContext(DataStream &stream, ContextMode mode)
448{
449 auto sMat = dynamic_cast< StructuralMaterial * >( domain->giveMaterial(slaveMat) );
450 auto mS = sMat->giveStatus(gp);
451 // save parent class status
452 mS->saveContext(stream, mode);
453}
454
455
456void
457LargeStrainMasterMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode)
458{
459 StructuralMaterialStatus :: restoreContext(stream, mode);
460}
461} // end namespace oofem
REGISTER_Material
#define REGISTER_Material(class)
Definition classfactory.h:135
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition femcmpnn.h:79
oofem::FloatArrayF::at
double & at(std::size_t i)
Definition floatarrayf.h:157
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatMatrixF::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrixf.h:224
oofem::FloatMatrix::jaco_
bool jaco_(FloatArray &eval, FloatMatrix &v, int nf)
Definition floatmatrix.C:1953
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::IntegrationPointStatus::gp
GaussPoint * gp
Associated integration point.
Definition integrationpointstatus.h:71
oofem::IntegrationPointStatus::saveContext
virtual void saveContext(DataStream &stream, ContextMode mode)
Definition integrationpointstatus.h:103
oofem::LargeStrainMasterMaterial::m
double m
Specifies the strain tensor.
Definition lsmastermat.h:72
oofem::LargeStrainMasterMaterial::constructTransformationMatrix
FloatMatrixF< 6, 6 > constructTransformationMatrix(const FloatMatrixF< 3, 3 > &eigenVectors) const
transformation matrices
Definition lsmastermat.C:147
oofem::LargeStrainMasterMaterial::slaveMat
int slaveMat
'slave' material model number.
Definition lsmastermat.h:70
oofem::LargeStrainMasterMaterial::constructL1L2TransformationMatrices
std::pair< FloatMatrixF< 6, 6 >, FloatMatrixF< 6, 6 > > constructL1L2TransformationMatrices(const FloatArrayF< 3 > &eigenValues, const FloatArrayF< 6 > &stress, double E1, double E2, double E3) const
Definition lsmastermat.C:196
oofem::MaterialStatus::updateYourself
void updateYourself(TimeStep *) override
Definition matstatus.h:104
oofem::MaterialStatus::initTempStatus
virtual void initTempStatus()
Definition matstatus.h:99
oofem::MaterialStatus::printOutputAt
void printOutputAt(FILE *file, TimeStep *tStep) const override
Print receiver's output to given stream.
Definition matstatus.h:93
oofem::Material::giveStatus
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Definition material.C:206
oofem::Material::initTempStatus
virtual void initTempStatus(GaussPoint *gp) const
Definition material.C:221
oofem::StructuralMaterialStatus::StructuralMaterialStatus
StructuralMaterialStatus(GaussPoint *g)
Constructor. Creates new StructuralMaterialStatus with IntegrationPoint g.
Definition structuralms.C:42
oofem::StructuralMaterial::StructuralMaterial
StructuralMaterial(int n, Domain *d)
Definition structuralmaterial.C:66
oofem::StructuralMaterial::convert_dSdE_2_dPdF_3D
static FloatMatrixF< 9, 9 > convert_dSdE_2_dPdF_3D(const FloatMatrixF< 6, 6 > &dSdE, const FloatArrayF< 6 > &S, const FloatArrayF< 9 > &F)
Definition structuralmaterial.C:431
oofem::StructuralMaterial::giveIPValue
int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override
Definition structuralmaterial.C:1996
IR_GIVE_OPTIONAL_FIELD
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Definition inputrecord.h:75
_IFT_LargeStrainMasterMaterial_m
#define _IFT_LargeStrainMasterMaterial_m
Definition lsmastermat.h:48
_IFT_LargeStrainMasterMaterial_slaveMat
#define _IFT_LargeStrainMasterMaterial_slaveMat
Definition lsmastermat.h:49
oofem::ContextMode
long ContextMode
Definition contextmode.h:43
oofem::to_voigt_strain_33
FloatArrayF< 6 > to_voigt_strain_33(const FloatMatrixF< 3, 3 > &t)
Definition floatmatrixf.h:848
oofem::Tdot
FloatMatrixF< N, P > Tdot(const FloatMatrixF< M, N > &a, const FloatMatrixF< M, P > &b)
Computes .
Definition floatmatrixf.h:677
oofem::TL
@ TL
Total Lagrange.
Definition fmode.h:44
oofem::from_voigt_form_9
FloatMatrixF< 3, 3 > from_voigt_form_9(const FloatArrayF< 9 > &v)
Definition floatmatrixf.h:774
oofem::from_voigt_stress_6
FloatMatrixF< 3, 3 > from_voigt_stress_6(const FloatArrayF< 6 > &v)
Definition floatmatrixf.h:806
oofem::to_voigt_form_33
FloatArrayF< 9 > to_voigt_form_33(const FloatMatrixF< 3, 3 > &t)
Definition floatmatrixf.h:787
oofem::dot
double dot(const FloatArray &x, const FloatArray &y)
Definition floatarray.C:1011
oofem::eye
FloatMatrixF< N, N > eye()
Constructs an identity matrix.
Definition floatmatrixf.h:944
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011