OOFEM 3.0
Loading...
Searching...
No Matches
solidshell.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 "sm/Elements/Shells/solidshell.h"
36#include <Materials/structuralms.h>
37#include "classfactory.h"
38#include "fei3dhexalin.h"
39#include "gausspoint.h"
40#include "gaussintegrationrule.h"
41#include "sm/CrossSections/structuralcrosssection.h"
42#include "parametermanager.h"
43#include "paramkey.h"
44
45namespace oofem {
46REGISTER_Element(SolidShell);
47
48FEI3dHexaLin SolidShell :: interpolation;
49ParamKey SolidShell :: IPK_SolidShell_EAS_type("eas_type");
50
51SolidShell :: SolidShell(int n, Domain *aDomain) : LSpace(n, aDomain)
52{
53 numberOfDofMans = 8;
54 numberOfGaussPoints = 8;
55 this->EAS_type = 0;
56}
57
58
59void
60SolidShell :: postInitialize()
61{
62 LSpace :: postInitialize();
63
64
65 int numEASparam = 0;
66
67 switch ( this->EAS_type) {
68 case 1:
69 numEASparam = 1;
70 break;
71
72 case 2:
73 numEASparam = 3;
74 break;
75 }
76
77 this->u_k.resize(numberOfDofMans*3);
78 this->u_k.zero();
79 this->fE.resize(numEASparam);
80 this->fE.zero();
81 this->KEC.resize(numEASparam, numberOfDofMans*3);
82 this->KEC.zero();
83 this->invKEE.resize(numEASparam, numEASparam);
84 this->invKEE.zero();
85 this->alpha.resize(numEASparam);
86 this->alpha.zero();
87}
88
89void
90SolidShell :: computeGaussPoints()
91// Sets up the array containing the four Gauss points of the receiver.
92{
93 if ( integrationRulesArray.size() == 0 ) {
94 integrationRulesArray.resize( 1 );
95 integrationRulesArray [ 0 ] = std::make_unique<GaussIntegrationRule>(1, this, 1, 6);
96 this->giveCrossSection()->setupIntegrationPoints(* integrationRulesArray [ 0 ], numberOfGaussPoints, this);
97 }
98}
99
100FEInterpolation *SolidShell :: giveInterpolation() const { return & interpolation; }
101
102
103void
104SolidShell :: initializeFrom(InputRecord &ir, int priority)
105{
106 ParameterManager &ppm = this->giveDomain()->elementPPM;
107 NLStructuralElement :: initializeFrom(ir, priority);
108 PM_UPDATE_PARAMETER(EAS_type, ppm, ir, this->number, IPK_SolidShell_EAS_type, priority) ;
109}
110
111
112
113void
114SolidShell :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
115// Returns the [ 6 x (8*3) ] strain-displacement matrix {B} of the receiver, eva-
116// luated at gp.
117// B matrix - 6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY :
118{
119 FEInterpolation *interp = this->giveInterpolation();
120 FloatMatrix dNdx;
121 const FloatArray &lCoords = gp->giveNaturalCoordinates();
122 interp->evaldNdx( dNdx, lCoords, FEIElementGeometryWrapper(this) );
123
124 answer.resize(6, dNdx.giveNumberOfRows() * 3);
125 answer.zero();
126
127 // Do ANS
128 // Need to evaluate strain at four points
129 FloatArray A = { 0.0, -1.0, 0.0};
130 FloatArray B = { 1.0, 0.0, 0.0};
131 FloatArray C = { 0.0, 1.0, 0.0};
132 FloatArray D = {-1.0, 0.0, 0.0};
133
134 FloatMatrix dNdxA, dNdxB, dNdxC, dNdxD;
135 interp->evaldNdx( dNdxA, A, FEIElementGeometryWrapper(this) );
136 interp->evaldNdx( dNdxB, B, FEIElementGeometryWrapper(this) );
137 interp->evaldNdx( dNdxC, C, FEIElementGeometryWrapper(this) );
138 interp->evaldNdx( dNdxD, D, FEIElementGeometryWrapper(this) );
139
140 FloatMatrix dNdx0;
141 interp->evaldNdx( dNdx0, {lCoords.at(1), lCoords.at(2), 0.0}, FEIElementGeometryWrapper(this) );
142
143 double NA = 0.5 * ( 1.0 - lCoords.at(2) );
144 double NC = 0.5 * ( 1.0 + lCoords.at(2) );
145 double NB = 0.5 * ( 1.0 - lCoords.at(1) );
146 double ND = 0.5 * ( 1.0 + lCoords.at(1) );
147
148 // E_xz = E(5) = N_A(gp) * B(A)(5,:) + N_C(gp) * B(C)(5,:)
149 // E_yz = E(4) = N_B(gp) * B(B)(4,:) + N_D(gp) * B(D)(4,:)
150
151 for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
152 answer.at(1, 3 * i - 2) = dNdx.at(i, 1);
153 answer.at(2, 3 * i - 1) = dNdx.at(i, 2);
154 answer.at(6, 3 * i - 2) = dNdx.at(i, 2);
155 answer.at(6, 3 * i - 1) = dNdx.at(i, 1);
156
157 // Evaluate thickness strain at the mid-surface
158 //answer.at(3, 3 * i - 0) = dNdx.at(i, 3);
159 answer.at(3, 3 * i - 0) = dNdx0.at(i, 3);
160
161 answer.at(4, 3 * i - 1) = NB * dNdxB.at(i, 3) + ND * dNdxD.at(i, 3);
162 answer.at(4, 3 * i - 0) = NB * dNdxB.at(i, 2) + ND * dNdxD.at(i, 2);
163
164 answer.at(5, 3 * i - 2) = NA * dNdxA.at(i, 3) + NC * dNdxC.at(i, 3);
165 answer.at(5, 3 * i - 0) = NA * dNdxA.at(i, 1) + NC * dNdxC.at(i, 1);
166 }
167}
168
169
170void
171SolidShell :: computeBHmatrixAt(FloatArray &lCoords, FloatMatrix &answer)
172// Returns the [ 6 x (8*3) ] strain-displacement matrix {B} of the receiver, eva-
173// luated at gp.
174// B matrix - 6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY :
175{
176#if 1
177 FEInterpolation *interp = this->giveInterpolation();
178 FloatMatrix dNdx;
179 interp->evaldNdx( dNdx, lCoords, FEIElementGeometryWrapper(this) );
180
181 answer.resize(9, dNdx.giveNumberOfRows() * 3);
182 answer.zero();
183 for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
184 answer.at(1, 3 * i - 2) = dNdx.at(i, 1); // du/dx
185 answer.at(2, 3 * i - 1) = dNdx.at(i, 2); // dv/dy
186 answer.at(3, 3 * i - 0) = dNdx.at(i, 3); // dw/dz
187 answer.at(4, 3 * i - 1) = dNdx.at(i, 3); // dv/dz
188 answer.at(7, 3 * i - 0) = dNdx.at(i, 2); // dw/dy
189 answer.at(5, 3 * i - 2) = dNdx.at(i, 3); // du/dz
190 answer.at(8, 3 * i - 0) = dNdx.at(i, 1); // dw/dx
191 answer.at(6, 3 * i - 2) = dNdx.at(i, 2); // du/dy
192 answer.at(9, 3 * i - 1) = dNdx.at(i, 1); // dv/dx
193 }
194
195#else
196 FEInterpolation *interp = this->giveInterpolation();
197 FloatMatrix dNdx;
198 interp->evaldNdx( dNdx, lCoords, FEIElementGeometryWrapper(this) );
199
200 answer.resize(9, dNdx.giveNumberOfRows() * 3);
201 answer.zero();
202
203 // Need to evaluate strain at four points
204 FloatArray A = { 0.0, -1.0, 0.0};
205 FloatArray B = { 1.0, 0.0, 0.0};
206 FloatArray C = { 0.0, 1.0, 0.0};
207 FloatArray D = {-1.0, 0.0, 0.0};
208
209 FloatMatrix dNdxA, dNdxB, dNdxC, dNdxD;
210 interp->evaldNdx( dNdxA, A, FEIElementGeometryWrapper(this) );
211 interp->evaldNdx( dNdxB, B, FEIElementGeometryWrapper(this) );
212 interp->evaldNdx( dNdxC, C, FEIElementGeometryWrapper(this) );
213 interp->evaldNdx( dNdxD, D, FEIElementGeometryWrapper(this) );
214
215 FloatMatrix dNdx0;
216 interp->evaldNdx( dNdx0, {lCoords.at(1), lCoords.at(2), 0.0}, FEIElementGeometryWrapper(this) );
217
218 double NA = 0.5 * ( 1.0 - lCoords.at(2) );
219 double NC = 0.5 * ( 1.0 + lCoords.at(2) );
220 double NB = 0.5 * ( 1.0 - lCoords.at(1) );
221 double ND = 0.5 * ( 1.0 + lCoords.at(1) );
222
223 // E_xz = E(5) = N_A(gp) * B(A)(5,:) + N_C(gp) * B(C)(5,:)
224 // E_yz = E(4) = N_B(gp) * B(B)(4,:) + N_D(gp) * B(D)(4,:)
225
226 for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
227 answer.at(1, 3 * i - 2) = dNdx.at(i, 1);
228 answer.at(2, 3 * i - 1) = dNdx.at(i, 2);
229 answer.at(6, 3 * i - 2) = dNdx.at(i, 2); // du/dy
230 answer.at(9, 3 * i - 1) = dNdx.at(i, 1); // dv/dx
231
232 // Evaluate thickness strain at the mid-surface
233 //answer.at(3, 3 * i - 0) = dNdx.at(i, 3);
234 answer.at(3, 3 * i - 0) = dNdx0.at(i, 3);
235
236 answer.at(4, 3 * i - 1) = NB * dNdxB.at(i, 3) + ND * dNdxD.at(i, 3);
237 answer.at(7, 3 * i - 0) = NB * dNdxB.at(i, 2) + ND * dNdxD.at(i, 2);
238
239 answer.at(5, 3 * i - 2) = NA * dNdxA.at(i, 3) + NC * dNdxC.at(i, 3);
240 answer.at(8, 3 * i - 0) = NA * dNdxA.at(i, 1) + NC * dNdxC.at(i, 1);
241 }
242#endif
243}
244
245
246void
247SolidShell :: computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
248// Returns the [ 6 x (8*3) ] strain-displacement matrix {B} of the receiver, eva-
249// luated at gp.
250// B matrix - 6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY :
251{
252 this->computeBHmatrixAt(gp->giveNaturalCoordinates(), answer);
253}
254
255
256void
257SolidShell :: computeEASBmatrixAt(GaussPoint *gp, FloatMatrix &answer)
258{
259 const FloatArray &lCoords = gp->giveNaturalCoordinates();
260 double xi = lCoords.at(1);
261 double eta = lCoords.at(2);
262 double zeta = lCoords.at(3);
263 FloatMatrix M;
264
265 switch ( this->EAS_type) {
266 case 1:
267 // alt 1
268 M.resize(6,1);
269 M.zero();
270 M.at(3,1) = zeta;
271 break;
272 case 2:
273 // alt 2
274 M.resize(6,3);
275 M.zero();
276 M.at(1,1) = xi;
277 M.at(2,2) = eta;
278 M.at(3,3) = zeta;
279 break;
280 }
281
282 FEInterpolation *interp = this->giveInterpolation();
283 FloatMatrix J0mat;
284 FloatArray center = {0.0, 0.0, 0.0};
285
286 interp->giveJacobianMatrixAt(J0mat, center, FEIElementGeometryWrapper(this) );
287// double J0 = interp->giveTransformationJacobian(center, FEIElementGeometryWrapper(this) );
288 double J0 = J0mat.giveDeterminant();
289 double J = interp->giveTransformationJacobian(lCoords, FEIElementGeometryWrapper(this) );
290
291 FloatMatrix T(6,6);
292 this->computeBondTransformationMatrix(T, J0mat);
293 //T.printYourself("T");
294 //T.beUnitMatrix();
295 answer.clear();
296 answer.plusProductUnsym(T, M, J/J0);
297}
298
299
300void
301SolidShell :: computeBondTransformationMatrix(FloatMatrix &answer, FloatMatrix &base)
302{
303 // given a bases {g1, g2, g3} compute the Bond (Voigt form) transformation matrix.
304 // TODO is this with OOFEM ordering of components?
305 FloatArray x, y, z;
306 x.beColumnOf(base,1);
307 y.beColumnOf(base,2);
308 z.beColumnOf(base,3);
309 answer = {
310 { x(0) * x(0), x(1) * x(1), x(2) * x(2), x(1) * x(2), x(0) * x(2), x(0) * x(1) },
311 { y(0) * y(0), y(1) * y(1), y(2) * y(2), y(1) * y(2), y(0) * y(2), y(0) * y(1) },
312 { z(0) * z(0), z(1) * z(1), z(2) * z(2), z(1) * z(2), z(0) * z(2), z(0) * z(1) },
313 { 2 * y(0) * z(0), 2 * y(1) * z(1), 2 * y(2) * z(2), y(2) * z(1) + y(1) * z(2), y(2) * z(0) + y(0) * z(2), y(1) * z(0) + y(0) * z(1) },
314 { 2 * x(0) * z(0), 2 * x(1) * z(1), 2 * x(2) * z(2), x(2) * z(1) + x(1) * z(2), x(2) * z(0) + x(0) * z(2), x(1) * z(0) + x(0) * z(1) },
315 { 2 * x(0) * y(0), 2 * x(1) * y(1), 2 * x(2) * y(2), x(2) * y(1) + x(1) * y(2), x(2) * y(0) + x(0) * y(2), x(1) * y(0) + x(0) * y(1) }
316 };
317}
318
319void
320SolidShell :: computeAlpha(FloatArray &answer, FloatArray &u)
321{
322 // compute alpha based on displacement update
323 FloatArray deltaU;
324 deltaU.beDifferenceOf(u,this->u_k);
325
326 FloatMatrix KEE_inv = this->invKEE;
327 FloatArray fE = this->fE;
328 FloatMatrix KEC = this->KEC;
329 FloatArray temp, deltaAlpha;
330 temp.beProductOf(KEC, deltaU);
331 temp.add(fE);
332 deltaAlpha.beProductOf(KEE_inv,temp);
333
334 FloatArray oldAlpha = this->alpha; // last converged values
335 answer = this->alpha - deltaAlpha;
336
337 // set current u-displcement
338 this->u_k = u;
339}
340
341void
342SolidShell :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
343{
344 if ( this->EAS_type ) {
345 FloatArray fC, fE, strain, u, vStrainC, vStrainE, vStress;
346 FloatMatrix KEE, KBE, BC, BE, D;
347 fE.clear();
348 fC.clear();
349
350 KEE.clear();
351
352 this->computeVectorOf(VM_Total, tStep, u);
353 FloatArray alpha;
354 this->computeAlpha(alpha, u );
355 StructuralCrossSection *cs = this->giveStructuralCrossSection();
356 for ( auto &gp: *this->giveDefaultIntegrationRulePtr() ) {
357 if ( nlGeometry == 0 ) {
358 computeBmatrixAt(gp, BC);
359 this->computeEASBmatrixAt(gp, BE);
360
361 vStrainC.beProductOf(BC, u);
362 vStrainE.beProductOf(BE, alpha);
363 vStrainC.add(vStrainE);
364 //this->computeStressVector(vStress, vStrainC, gp, tStep);
365 cs->giveRealStresses(vStress, gp, vStrainC, tStep);
366
367 } else if ( nlGeometry == 1 ) { // First Piola-Kirchhoff stress
368
369 this->computeBHmatrixAt(gp, BC);
370 FloatArray vFC, VFE;
371 vFC.beProductOf(BC, u);
372 //vFE.beProductOf(BE, alpha);
373 //vFC.add(VFE);
374 cs->giveFirstPKStresses(vStress, gp, vFC, tStep);
375
376 } else if ( nlGeometry == 2 ) { // Second Piola-Kirchhoff stress
377
378 this->computeBEmatrixAt(gp, BC, tStep);
379 this->computeEASBmatrixAt(gp, BE);
380
381 this->computeEVector(vStrainC, gp->giveNaturalCoordinates(), u);
382 vStrainE.beProductOf(BE, alpha);
383 vStrainC.add(vStrainE);
384 cs->giveRealStresses(vStress, gp, vStrainC, tStep);
385 }
386
387 double dV = this->computeVolumeAround(gp);
388
389 // Compute nodal internal forces at nodes as f = B^T*Stress dV
390 fC.plusProduct(BC, vStress, dV);
391 fE.plusProduct(BE, vStress, dV);
392 }
393 this->fE = fE;
394
395 FloatMatrix KEE_inv, KEC;
396 KEE_inv = this->invKEE;
397 KEC = this->KEC;
398
399 FloatArray temp, f;
400 temp.beProductOf(KEE_inv,fE);
401 answer = fC;
402 answer.plusProduct(KEC, temp, -1.0);
403 } else {
404 LSpace :: giveInternalForcesVector(answer, tStep, useUpdatedGpRecord);
405 }
406}
407
408
409void
410SolidShell :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
411{
412 if ( this->EAS_type ) {
413 FloatArray fC, fE, strain, u, vStrainC, vStrainE, vStress;
414 FloatMatrix KEE, KEC, KCC, KCC_geo, BC, BE, D, DBC, DBE;
415
416 KEC.clear();
417 KEE.clear();
418 KCC.clear();
419
420 StructuralCrossSection *cs = this->giveStructuralCrossSection();
421 for ( GaussPoint *gp: *this->giveDefaultIntegrationRulePtr() ) {
422 if ( nlGeometry == 0 ) {
423 this->computeBmatrixAt(gp, BC);
424 this->computeEASBmatrixAt(gp, BE);
425 double dV = this->computeVolumeAround(gp);
426
427 this->computeConstitutiveMatrixAt(D, rMode, gp, tStep);
428 DBC.beProductOf(D, BC);
429 DBE.beProductOf(D, BE);
430
431 KCC.plusProductUnsym(BC, DBC, dV);
432 KEC.plusProductUnsym(BE, DBC, dV);
433 KEE.plusProductUnsym(BE, DBE, dV);
434
435 } else if ( nlGeometry == 1 ) {
436 this->computeBHmatrixAt(gp, BC);
437 cs->giveStiffnessMatrix_dPdF(D, rMode, gp, tStep);
438
439 } else if ( nlGeometry == 2 ) {
440 this->computeBEmatrixAt(gp, BC, tStep);
441 this->computeEASBmatrixAt(gp, BE);
442 double dV = this->computeVolumeAround(gp);
443
444 this->computeConstitutiveMatrixAt(D, rMode, gp, tStep);
445 DBC.beProductOf(D, BC);
446 DBE.beProductOf(D, BE);
447
448 KCC.plusProductUnsym(BC, DBC, dV);
449 KEC.plusProductUnsym(BE, DBC, dV);
450 KEE.plusProductUnsym(BE, DBE, dV);
451
452 //TODO add geometric stiffness
453 this->computeGeometricStiffness(KCC_geo, gp, tStep);
454 KCC.add(dV, KCC_geo);
455 }
456 }
457
458 FloatMatrix KEE_inv;
459 KEE_inv.beInverseOf(KEE);
460 this->invKEE = KEE_inv;
461 this->KEC = KEC;
462
463 answer = KCC;
464
465 FloatMatrix K, tempmat;
466 tempmat.beProductOf(KEE_inv, KEC);
467 answer.plusProductUnsym(KEC, tempmat, -1.0);
468 } else {
469 LSpace :: computeStiffnessMatrix(answer, rMode, tStep);
470 }
471}
472
473
474void
475SolidShell :: computeBEmatrixAt(GaussPoint *gp, FloatMatrix &answer, TimeStep *tStep)
476{
477 FloatMatrix dN, F;
478 FloatArray vF;
479 // compute the derivatives of shape functions
480 this->interpolation.evaldNdx( dN, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
481
482 FloatArray u;
483 this->computeVectorOf(VM_Total, tStep, u); // solution vector
484
485 this->computeFVector(vF, gp->giveNaturalCoordinates(), u);
486 F.beMatrixForm(vF);
487
488 answer.resize(6, 24);
489 answer.zero();
490#if 1
491 for ( int i = 1, col = 0; i <= dN.giveNumberOfRows(); i++, col += 3 ) {
492 for ( int j = 1; j <= 3; j++ ) {
493 answer.at(1,col + j) = F.at(j,1) * dN.at(i, 1);
494 answer.at(2,col + j) = F.at(j,2) * dN.at(i, 2);
495 answer.at(3,col + j) = F.at(j,3) * dN.at(i, 3);
496 answer.at(4,col + j) = F.at(j,2) * dN.at(i, 3) + F.at(j,3) * dN.at(i, 2);
497 answer.at(5,col + j) = F.at(j,1) * dN.at(i, 3) + F.at(j,3) * dN.at(i, 1);
498 answer.at(6,col + j) = F.at(j,1) * dN.at(i, 2) + F.at(j,2) * dN.at(i, 1);
499 }
500 }
501#else
502
503 FEInterpolation *interp = this->giveInterpolation();
504 FloatMatrix dNdx;
505 const FloatArray &lCoords = gp->giveNaturalCoordinates();
506 interp->evaldNdx( dNdx, lCoords, FEIElementGeometryWrapper(this) );
507// Do ANS
508 // Need to evaluate strain at four points
509 FloatArray A = { 0.0, -1.0, 0.0};
510 FloatArray B = { 1.0, 0.0, 0.0};
511 FloatArray C = { 0.0, 1.0, 0.0};
512 FloatArray D = {-1.0, 0.0, 0.0};
513
514 FloatMatrix dNdxA, dNdxB, dNdxC, dNdxD;
515 interp->evaldNdx( dNdxA, A, FEIElementGeometryWrapper(this) );
516 interp->evaldNdx( dNdxB, B, FEIElementGeometryWrapper(this) );
517 interp->evaldNdx( dNdxC, C, FEIElementGeometryWrapper(this) );
518 interp->evaldNdx( dNdxD, D, FEIElementGeometryWrapper(this) );
519
520 // Eval deformation gradients
521 FloatArray vFA, vFB, vFC, vFD;
522 this->computeFVector(vFA, A, u);
523 this->computeFVector(vFB, B, u);
524 this->computeFVector(vFC, C, u);
525 this->computeFVector(vFD, D, u);
526 FloatMatrix FA, FB, FC, FD;
527 FA.beMatrixForm(vFA);
528 FB.beMatrixForm(vFB);
529 FC.beMatrixForm(vFC);
530 FD.beMatrixForm(vFD);
531
532 FloatMatrix dNdx0;
533 interp->evaldNdx( dNdx0, {lCoords.at(1), lCoords.at(2), 0.0}, FEIElementGeometryWrapper(this) );
534
535 double NA = 0.5 * ( 1.0 - lCoords.at(2) );
536 double NC = 0.5 * ( 1.0 + lCoords.at(2) );
537 double NB = 0.5 * ( 1.0 - lCoords.at(1) );
538 double ND = 0.5 * ( 1.0 + lCoords.at(1) );
539
540 // E_xz = E(5) = N_A(gp) * B(A)(5,:) + N_C(gp) * B(C)(5,:)
541 // E_yz = E(4) = N_B(gp) * B(B)(4,:) + N_D(gp) * B(D)(4,:)
542
543 for ( int i = 1, col = 0; i <= dN.giveNumberOfRows(); i++, col += 3 ) {
544 for ( int j = 1; j <= 3; j++ ) {
545 answer.at(1,col + j) = F.at(j,1) * dN.at(i, 1);
546 answer.at(2,col + j) = F.at(j,2) * dN.at(i, 2);
547 answer.at(3,col + j) = F.at(j,3) * dN.at(i, 3);
548 answer.at(6,col + j) = F.at(j,1) * dN.at(i, 2) + F.at(j,2) * dN.at(i, 1);
549
550 answer.at(4,col + j) = NB * ( FB.at(j,2) * dNdxB.at(i, 3) + FB.at(j,3) * dNdxB.at(i, 2) ) + ND * ( FD.at(j,2) * dNdxD.at(i, 3) + FD.at(j,3) * dNdxD.at(i, 2) ) ;
551 answer.at(5,col + j) = NA * ( FA.at(j,1) * dNdxA.at(i, 3) + FA.at(j,3) * dNdxA.at(i, 1) ) + NC * ( FC.at(j,1) * dNdxC.at(i, 3) + FC.at(j,3) * dNdxC.at(i, 1) );
552 }
553 }
554#endif
555}
556
557
558void
559SolidShell :: computeFVector(FloatArray &answer, FloatArray &lCoords, FloatArray &ae)
560{
561 FloatMatrix B;
562 this->computeBHmatrixAt(lCoords, B);
563 answer.beProductOf(B, ae);
564
565 answer.at(1) += 1.0;
566 answer.at(2) += 1.0;
567 answer.at(3) += 1.0;
568}
569
570
571void
572SolidShell :: computeEVector(FloatArray &answer, FloatArray &lCoords, FloatArray &u)
573{
574 FloatArray vF;
575 this->computeFVector(vF, lCoords, u);
576#if 0
577 FloatMatrix F, E;
578 F.beMatrixForm(vF);
579 E.beTProductOf(F, F);
580 E.at(1, 1) -= 1.0;
581 E.at(2, 2) -= 1.0;
582 E.at(3, 3) -= 1.0;
583 E.times(0.5);
584 answer.beSymVectorFormOfStrain(E);
585#else
586 FloatMatrix F, E;
587 F.beMatrixForm(vF);
588 E.beTProductOf(F, F);
589 E.at(1, 1) -= 1.0;
590 E.at(2, 2) -= 1.0;
591 E.at(3, 3) -= 1.0;
592 E.times(0.5);
593 answer.beSymVectorFormOfStrain(E);
594
595 // Need to evaluate strain at four points
596 FloatArray A = { 0.0, -1.0, 0.0};
597 FloatArray B = { 1.0, 0.0, 0.0};
598 FloatArray C = { 0.0, 1.0, 0.0};
599 FloatArray D = {-1.0, 0.0, 0.0};
600
601 // Eval deformation gradients
602 FloatArray vFA, vFB, vFC, vFD;
603 this->computeFVector(vFA, A, u);
604 this->computeFVector(vFB, B, u);
605 this->computeFVector(vFC, C, u);
606 this->computeFVector(vFD, D, u);
607 FloatMatrix FA, FB, FC, FD, EA, EB, EC, ED;
608 FA.beMatrixForm(vFA);
609 FB.beMatrixForm(vFB);
610 FC.beMatrixForm(vFC);
611 FD.beMatrixForm(vFD);
612
613 EA.beTProductOf(FA,FA); EA.times(0.5);
614 EB.beTProductOf(FB,FB); EB.times(0.5);
615 EC.beTProductOf(FC,FC); EC.times(0.5);
616 ED.beTProductOf(FD,FD); ED.times(0.5);
617
618// FloatMatrix dNdx0;
619// interp->evaldNdx( dNdx0, {lCoords.at(1), lCoords.at(2), 0.0}, FEIElementGeometryWrapper(this) );
620
621 double NA = 0.5 * ( 1.0 - lCoords.at(2) );
622 double NC = 0.5 * ( 1.0 + lCoords.at(2) );
623 double NB = 0.5 * ( 1.0 - lCoords.at(1) );
624 double ND = 0.5 * ( 1.0 + lCoords.at(1) );
625
626 answer.at(4) = NB * EB.at(2,3) + ND * ED.at(2,3);
627 answer.at(5) = NA * EA.at(2,3) + NC * EC.at(2,3);
628#endif
629}
630
631
632void
633SolidShell :: computeGeometricStiffness(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
634{
635 answer.resize(24,24);
636 answer.zero();
637
638 StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );
639 if ( !matStat ) { // no matstat created yet
640 return;
641 }
642 FloatArray vS;
643 vS = matStat->giveStressVector();
644
645 FloatMatrix S, dNdx, temp, G_ij;
646 S.beMatrixFormOfStress(vS);
647 this->interpolation.evaldNdx( dNdx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
648
649 temp.beProductTOf(S,dNdx);
650 G_ij.beProductOf(dNdx, temp);
651 FloatMatrix K_ij(3,3);
652 K_ij.zero();
653
654 for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
655 for ( int j = i; j <= dNdx.giveNumberOfRows(); j++ ) {
656 K_ij.at(1,1) = K_ij.at(2,2) = K_ij.at(3,3) = G_ij.at(i,j);
657 answer.assemble(K_ij, { (i-1)*3 + 1, (i-1)*3 + 2, (i-1)*3 + 3 }, { (j-1)*3 + 1, (j-1)*3 + 2, (j-1)*3 + 3 } );
658 }
659 }
660 answer.symmetrized();
661}
662
663} // end namespace oofem
E
#define E(a, b)
REGISTER_Element
#define REGISTER_Element(class)
Definition classfactory.h:131
oofem::Element::numberOfDofMans
int numberOfDofMans
Number of dofmanagers.
Definition element.h:136
oofem::Element::computeVectorOf
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Definition element.C:103
oofem::Element::integrationRulesArray
std::vector< std ::unique_ptr< IntegrationRule > > integrationRulesArray
Definition element.h:157
oofem::Element::numberOfGaussPoints
int numberOfGaussPoints
Definition element.h:175
oofem::Element::giveCrossSection
CrossSection * giveCrossSection()
Definition element.C:534
oofem::Element::giveDefaultIntegrationRulePtr
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Definition element.h:886
oofem::FEInterpolation::giveTransformationJacobian
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
Definition feinterpol.C:81
oofem::FEInterpolation::giveJacobianMatrixAt
virtual void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
Definition feinterpol.h:284
oofem::FEInterpolation::evaldNdx
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const =0
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition femcmpnn.h:97
oofem::FEMComponent::number
int number
Component number.
Definition femcmpnn.h:77
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::plusProduct
void plusProduct(const FloatMatrix &b, const FloatArray &s, double dV)
Definition floatarray.C:288
oofem::FloatArray::beSymVectorFormOfStrain
void beSymVectorFormOfStrain(const FloatMatrix &aMatrix)
Definition floatarray.C:1081
oofem::FloatArray::beColumnOf
void beColumnOf(const FloatMatrix &mat, int col)
Definition floatarray.C:1182
oofem::FloatArray::beDifferenceOf
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Definition floatarray.C:403
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:689
oofem::FloatArray::add
void add(const FloatArray &src)
Definition floatarray.C:218
oofem::FloatMatrix::times
void times(double f)
Definition floatmatrix.C:1155
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Definition floatmatrix.C:962
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
oofem::FloatMatrix::clear
*Sets size of receiver to be an empty matrix It will have zero rows and zero columns size void clear()
Definition floatmatrix.h:540
oofem::FloatMatrix::beProductTOf
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:431
oofem::FloatMatrix::plusProductUnsym
void plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)
Definition floatmatrix.C:691
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:365
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition floatmatrix.C:1064
oofem::FloatMatrix::beInverseOf
bool beInverseOf(const FloatMatrix &src)
Definition floatmatrix.C:768
oofem::FloatMatrix::giveDeterminant
double giveDeterminant() const
Definition floatmatrix.C:1085
oofem::FloatMatrix::beMatrixForm
void beMatrixForm(const FloatArray &aArray)
Definition floatmatrix.C:1689
oofem::FloatMatrix::giveNumberOfRows
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition floatmatrix.h:200
oofem::FloatMatrix::assemble
void assemble(const FloatMatrix &src, const IntArray &loc)
Definition floatmatrix.C:181
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::FloatMatrix::beTProductOf
void beTProductOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:398
oofem::GaussPoint::giveNaturalCoordinates
const FloatArray & giveNaturalCoordinates() const
Returns coordinate array of receiver.
Definition gausspoint.h:138
oofem::GaussPoint::giveMaterialStatus
IntegrationPointStatus * giveMaterialStatus(IntegrationPointStatusIDType key=IPSID_Default)
Definition gausspoint.h:204
oofem::LSpace::LSpace
LSpace(int n, Domain *d)
Definition lspace.C:62
oofem::NLStructuralElement::nlGeometry
int nlGeometry
Flag indicating if geometrical nonlinearities apply.
Definition nlstructuralelement.h:91
oofem::SolidShell::IPK_SolidShell_EAS_type
static ParamKey IPK_SolidShell_EAS_type
Definition solidshell.h:61
oofem::SolidShell::interpolation
static FEI3dHexaLin interpolation
Definition solidshell.h:60
oofem::SolidShell::computeAlpha
void computeAlpha(FloatArray &answer, FloatArray &u)
Definition solidshell.C:320
oofem::SolidShell::invKEE
FloatMatrix invKEE
Definition solidshell.h:92
oofem::SolidShell::KEC
FloatMatrix KEC
Definition solidshell.h:93
oofem::SolidShell::computeBEmatrixAt
virtual void computeBEmatrixAt(GaussPoint *gp, FloatMatrix &answer, TimeStep *tStep)
Definition solidshell.C:475
oofem::SolidShell::u_k
FloatArray u_k
Definition solidshell.h:91
oofem::SolidShell::computeEASBmatrixAt
virtual void computeEASBmatrixAt(GaussPoint *gp, FloatMatrix &answer)
Definition solidshell.C:257
oofem::SolidShell::computeBondTransformationMatrix
void computeBondTransformationMatrix(FloatMatrix &answer, FloatMatrix &base)
Definition solidshell.C:301
oofem::SolidShell::computeBmatrixAt
void computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int lowerIndx=1, int upperIndx=ALL_STRAINS) override
Definition solidshell.C:114
oofem::SolidShell::computeEVector
void computeEVector(FloatArray &answer, FloatArray &lCoords, FloatArray &ae)
Definition solidshell.C:572
oofem::SolidShell::computeFVector
void computeFVector(FloatArray &answer, FloatArray &lCoords, FloatArray &ae)
Definition solidshell.C:559
oofem::SolidShell::alpha
FloatArray alpha
Definition solidshell.h:90
oofem::SolidShell::computeBHmatrixAt
void computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer) override
Definition solidshell.C:247
oofem::SolidShell::computeGeometricStiffness
void computeGeometricStiffness(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)
Definition solidshell.C:633
oofem::SolidShell::giveInterpolation
FEInterpolation * giveInterpolation() const override
Definition solidshell.C:100
oofem::SolidShell::x
void x(int arg1)
oofem::SolidShell::fE
FloatArray fE
Definition solidshell.h:94
oofem::Structural3DElement::computeConstitutiveMatrixAt
void computeConstitutiveMatrixAt(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep) override
Definition structural3delement.C:252
oofem::Structural3DElement::computeVolumeAround
double computeVolumeAround(GaussPoint *gp) override
Definition structural3delement.C:331
oofem::StructuralCrossSection::giveRealStresses
FloatArray giveRealStresses(const FloatArray &reducedStrain, GaussPoint *gp, TimeStep *tStep) const
Definition structuralcrosssection.C:44
oofem::StructuralCrossSection::giveFirstPKStresses
virtual FloatArray giveFirstPKStresses(const FloatArray &reducedF, GaussPoint *gp, TimeStep *tStep) const
Definition structuralcrosssection.C:83
oofem::StructuralElement::giveStructuralCrossSection
StructuralCrossSection * giveStructuralCrossSection()
Helper function which returns the structural cross-section for the element.
Definition structuralelement.C:1236
oofem::StructuralMaterialStatus::giveStressVector
const FloatArray & giveStressVector() const
Returns the const pointer to receiver's stress vector.
Definition structuralms.h:105
S
#define S(p)
Definition mdm.C:469
PM_UPDATE_PARAMETER
#define PM_UPDATE_PARAMETER(_val, _pm, _ir, _componentnum, _paramkey, _prio)
Definition parametermanager.h:56
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