OOFEM 3.0
Loading...
Searching...
No Matches
shell7basePhFi.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 "Shell7BasePhFi.h"
36#include "node.h"
37#include "load.h"
38#include "structuralms.h"
39#include "mathfem.h"
40#include "domain.h"
41#include "equationid.h"
42#include "gaussintegrationrule.h"
43#include "gausspoint.h"
44#include "feinterpol3d.h"
45#include "fei3dtrquad.h"
46#include "boundaryload.h"
47#include "constantpressureload.h"
48#include "constantsurfaceload.h"
49#include "vtkxmlexportmodule.h"
50#include "fracturemanager.h"
51
52#include "masterdof.h"
53
54#include <fstream>
55
56namespace oofem {
57
58const int nLayers = 5; //@todo: Generalize!
59const double disturB = 1e-8; //@todo: Generalize!
60
61
62Shell7BasePhFi :: Shell7BasePhFi(int n, Domain *aDomain) : Shell7Base(n, aDomain), PhaseFieldElement(n, aDomain){
63 this->numberOfLayers = nLayers;
64}
65
66void Shell7BasePhFi :: initializeFrom(InputRecord &ir, int priority)
67{
68 Shell7Base :: initializeFrom(ir, priority);
69}
70
71
72
73void
74Shell7BasePhFi :: postInitialize()
75{
76 Shell7Base :: postInitialize();
77}
78
79
80
81void
82Shell7BasePhFi :: giveDofManDofIDMask(int inode, EquationID ut, IntArray &answer) const
83{
84 // returns the total DofIDMask: [shell7Base IDMask, damage var. IdMask]
85 IntArray answer_d;
86 this->giveDofManDofIDMask_u(answer);
87 this->giveDofManDofIDMask_d(answer_d);
88 answer.followedBy(answer_d);
89
90}
91
92void
93Shell7BasePhFi :: giveDofManDofIDMask_u(IntArray &answer) const
94{
95 answer.setValues(7, D_u, D_v, D_w, W_u, W_v, W_w, Gamma);
96}
97
98void
99Shell7BasePhFi :: giveDofManDofIDMask_d(IntArray &answer) const
100{
101 // Returns [nextFreeDofID, nextFreeDofID+1, ..., nextFreeDofID + numberOfLayers]
102 answer.resize(0);
103 int sID = this->domain->giveNextFreeDofID(0);
104 for (int i = 0; i < numberOfLayers; i++) {
105 answer.followedBy(sID + i);
106 }
107}
108
109
110double
111Shell7BasePhFi :: computeDamageInLayerAt(int layer, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
112{
113 // d = N_d * a_d
114 FloatArray dVec, dVecRed;
115
116 this->computeDamageUnknowns(dVec, valueMode, stepN); // should be a vector with all damage nodal values for this element
117 // ordered such that all damage dofs associated with node 1 comes first
118 FloatArray Nvec, lcoords;
119 IntArray indx = computeDamageIndexArray(layer); // Since dVec only contains damage vars this index vector should be 1 3 5 7 9 11 for the first layer (with 2 layers in the cs)
120
121 dVecRed.beSubArrayOf(dVec, indx);
122 this->giveInterpolation()->evalN(Nvec, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
123
124 // Make sure returned damage is always between [0,1]
125 double d_temp = Nvec.dotProduct(dVecRed);
126 //if ( d_temp < 0.0 ) {
127 // return 0.0;
128 //} else if ( d_temp > 1.0 ) {
129 // return 1.0;
130 //} else {
131 return d_temp;
132 //}
133}
134
135double
136Shell7BasePhFi :: computeOldDamageInLayerAt(int layer, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
137{
138 // d = N_d * a_d
139 FloatArray dVec, dVecRed, ddVec;
140
141 this->computeDamageUnknowns(dVec, valueMode, stepN); // should be a vector with all damage nodal values for this element
142 // ordered such that all damage dofs associated with node 1 comes first
143 this->computeDamageUnknowns(ddVec, VM_Incremental, stepN);
144
145 dVec.subtract(ddVec);
146
147 FloatArray Nvec, lcoords;
148 IntArray indx = computeDamageIndexArray(layer); // Since dVec only contains damage vars this index vector should be 1 3 5 7 9 11 for the first layer (with 2 layers in the cs)
149
150 dVecRed.beSubArrayOf(dVec, indx);
151 this->giveInterpolation()->evalN(Nvec, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
152
153 // Make sure returned damage is always between [0,1]
154 double d_temp = Nvec.dotProduct(dVecRed);
155 //if ( d_temp < 0.0 ) {
156 // return 0.0;
157 //} else if ( d_temp > 1.0 ) {
158 // return 1.0;
159 //} else {
160 return d_temp;
161 //}
162}
163double
164Shell7BasePhFi :: computeDamageInLayerAt_dist(int layer, int index, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
165{
166 // d = N_d * a_d
167 FloatArray dVec, dVecRed;
168
169 this->computeDamageUnknowns(dVec, valueMode, stepN); // should be a vector with all damage nodal values for this element
170 // ordered such that all damage dofs associated with node 1 comes first
171 if (valueMode == VM_Total)
172 {
173 dVec.at(index) = dVec.at(index) + disturB;
174 } else if (valueMode == VM_Incremental)
175 {
176 dVec.at(index) = dVec.at(index) + disturB/(stepN->giveTimeIncrement());
177 }
178
179
180 FloatArray Nvec, lcoords;
181 IntArray indx = computeDamageIndexArray(layer); // Since dVec only contains damage vars this index vector should be 1 3 5 7 9 11 for the first layer (with 2 layers in the cs)
182
183 dVecRed.beSubArrayOf(dVec, indx);
184 this->giveInterpolation()->evalN(Nvec, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
185
186 // Make sure returned damage is always between [0,1]
187 double d_temp = Nvec.dotProduct(dVecRed);
188 // if ( d_temp < 0.0 ) {
189 // return 0.0;
190 // } else if ( d_temp > 1.0 ) {
191 // return 1.0;
192 // } else {
193 return d_temp;
194 // }
195}
196
197#if 0
198double
199Shell7BasePhFi :: computeDamageInLayerAt_dist(int layer, int index, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
200{
201 // d = N_d * a_d
202 FloatArray dVec, dVecRed;
203
204 this->computeDamageUnknowns(dVec, valueMode, stepN); // should be a vector with all damage nodal values for this element
205 // ordered such that all damage dofs associated with node 1 comes first
206 dVec.at(index) = dVec.at(index) + disturB;
207 FloatArray Nvec, lcoords;
208 IntArray indx = computeDamageIndexArray(layer); // Since dVec only contains damage vars this index vector should be 1 3 5 7 9 11 for the first layer (with 2 layers in the cs)
209
210 dVecRed.beSubArrayOf(dVec, indx);
211 this->giveInterpolation()->evalN(Nvec, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
212
213 // Make sure returned damage is always between [0,1]
214 double d_temp = Nvec.dotProduct(dVecRed);
215 if ( d_temp < 0.0 ) {
216 return 0.0;
217 } else if ( d_temp > 1.0 ) {
218 return 1.0;
219 } else {
220 return d_temp;
221 }
222}
223#endif
224
225IntArray
226Shell7BasePhFi :: computeDamageIndexArray(int layer)
227{
228 int numberOfNodes = this->giveNumberOfDofManagers();
229
230 IntArray indx(numberOfNodes);
231
232 for (int i = 1; i <= numberOfNodes; i++)
233 {
234 indx.at(i) = layer + (i-1)*this->layeredCS->giveNumberOfLayers(); // NEW JB
235 }
236 return indx;
237}
238
239double
240Shell7BasePhFi :: computeGInLayer(int layer, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
241{
242 // computes g = (1-d)^2 + r0
243 //double d = this->computeDamageInLayerAt(layer, gp, valueMode, stepN);
244 double d = this->computeOldDamageInLayerAt(layer, gp, valueMode, stepN);
245 //double d = this->computeDamageInLayerAt(layer, gp, valueMode, stepN);
246 double r0 = 1.0e-10;
247
248 return (1.0 - d) * (1.0 - d) + r0;
249
250}
251
252double
253Shell7BasePhFi :: computeGInLayer_dist(int layer, int indx, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
254{
255 // computes g = (1-d)^2 + r0
256 double d = this->computeDamageInLayerAt_dist(layer, indx, gp, valueMode, stepN);
257 double r0 = 1.0e-10;
258
259 return (1.0 - d) * (1.0 - d) + r0;
260
261}
262
263double
264Shell7BasePhFi :: computeGprimInLayer(int layer, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
265{
266 // compute g' =-2*(1-d)
267 double d = this->computeDamageInLayerAt(layer, gp, valueMode, stepN);
268 return -2.0 * (1.0 - d);
269
270}
271
272double
273Shell7BasePhFi :: computeGprimInLayer_dist(int layer, int indx, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
274{
275 // compute g' =-2*(1-d)
276 double d = this->computeDamageInLayerAt_dist(layer, indx, gp, valueMode, stepN);
277 return -2.0 * (1.0 - d);
278
279}
280
281double
282Shell7BasePhFi :: computeGbisInLayer(int layer, GaussPoint *gp, ValueModeType valueMode, TimeStep *stepN)
283{
284 // compute g'' = D(-2*(1-d))/Dd = 2
285 //@todo this method is trivial but here if one wants to change the expression for g
286 //return 2.0;
287 double d = this->computeDamageInLayerAt(layer, gp, valueMode, stepN);
288 //if ( d < 0.0 ) {
289 // return 0.0;
290 //} else if ( d > 1.0 ) {
291 // return 0.0;
292 //} else {
293 return 2.0;
294 //}
295}
296
297int
298Shell7BasePhFi :: giveNumberOfDofs()
299{
300 return this->giveNumberOfuDofs() + this->giveNumberOfdDofs();
301}
302
303int
304Shell7BasePhFi :: giveNumberOfuDofs()
305{
306 return 7 * this->giveNumberOfDofManagers();
307}
308
309int
310Shell7BasePhFi :: giveNumberOfdDofs()
311{
312 return this->giveNumberOfDofManagers() * this->numberOfLayers;
313}
314
315
316// Tangent matrices
317
318void
319Shell7BasePhFi :: computeStiffnessMatrix_uu(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
320{
321 Shell7Base :: computeStiffnessMatrix(answer, rMode, tStep);
322}
323
324
325void
326Shell7BasePhFi :: computeStiffnessMatrix_ud(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep) {
327
328 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
329 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
330 int ndofs_d = ndofs - ndofs_u;
331
332 answer.resize( ndofs_u, ndofs_d);
333 answer.zero();
334
335 int numberOfLayers = this->layeredCS->giveNumberOfLayers(); // conversion of types
336
337 FloatArray fu(ndofs_u), fd(ndofs_d), ftemp_u, sectionalForces, sectionalForces_dv;
338 FloatArray genEps, genEpsD, solVec, lCoords, Nd;
339 FloatMatrix B, Bd, K_ud(ndofs_u, this->numberOfDofMans), K_temp;
340 this->giveUpdatedSolutionVector(solVec, tStep); // => x, m, gam
341
342 IntArray ordering_temp(ndofs_u);
343 for ( int i = 1; i <= ndofs_u; i++ ) {
344 ordering_temp.at(i) = i;
345 }
346
347 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
348 fu.zero();
349 K_ud.zero();
350 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
351 Material *mat = domain->giveMaterial( this->layeredCS->giveLayerMaterial(layer) );
352 IntArray indx_d = computeDamageIndexArray(layer);
353 FloatArray dVecTotal, dVecLayer, Ddam_Dxi;
354
355 this->computeDamageUnknowns(dVecTotal, VM_Total, tStep); // vector with all damage nodal values for this element
356 dVecLayer.beSubArrayOf(dVecTotal, indx_d); // vector of damage variables for a given layer
357
358 for ( int j = 1; j <= iRuleL->giveNumberOfIntegrationPoints(); j++ ) {
359 GaussPoint *gp = iRuleL->getIntegrationPoint(j - 1);
360 lCoords = *gp->giveCoordinates();
361 this->computeBmatrixAt(lCoords, B);
362 this->computeNdvectorAt(lCoords, Nd); // (1x6)
363
364 this->computeGeneralizedStrainVectorNew(genEpsD, solVec, B);
365 this->computeGeneralizedStrainVectorNew(genEps, solVec, B); // used for computing the stress
366
367 double zeta = giveGlobalZcoord( *gp->giveCoordinates() );
368 this->computeSectionalForcesAt(sectionalForces, gp, mat, tStep, genEps, genEpsD, zeta); // these are per unit volume
369
370 // Computation of tangent: K_ud = \int Bu^t* (sec. forces) * g' * Nd
371 ftemp_u.beTProductOf(B,sectionalForces);
372 double dV = this->computeVolumeAroundLayer(gp, layer);
373 double Gprim = computeGprimInLayer(layer, gp, VM_Total, tStep);
374
375 fu.add(dV*Gprim, ftemp_u);
376 K_temp.beDyadicProductOf(fu, Nd);
377 K_ud.add(K_temp);
378 }
379
380 // Just place the columns in the right order
381 answer.assemble(K_ud, ordering_temp, indx_d);
382 //ordering_temp.printYourself();
383 //indx_d.printYourself();
384 }
385
386}
387
388void
389Shell7BasePhFi :: computeStiffnessMatrix_du(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep) {
390
391 answer.resize( this->numberOfDofMans*this->layeredCS->giveNumberOfLayers(), 42 );
392 answer.zero();
393
394}
395
396void
397Shell7BasePhFi :: computeStiffnessMatrix_dd(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
398{
399 // Computation of tangent: K_dd = \int Nd^t * ( -kp*neg_Mac'(alpha_dot)/delta_t + g_c/l + G''*Psibar) * Nd +
400 // \int Bd^t * ( g_c * l * [G^1 G^2]^t * [G^1 G^2] ) * Bd
401 // = K_dd1 + K_dd2
402 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
403 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
404 int ndofs_d = ndofs - ndofs_u;
405
406 answer.resize( ndofs_d, ndofs_d );
407 answer.zero();
408
409 int numberOfLayers = this->layeredCS->giveNumberOfLayers();
410
411 FloatArray lCoords;
412 FloatMatrix Bd, Nd, temp(this->giveNumberOfDofManagers(),this->giveNumberOfDofManagers());
413
414 IntArray ordering_temp(ndofs_u); // [1, 2, ..., 42]
415 for ( int i = 1; i <= ndofs_u; i++ ) {
416 ordering_temp.at(i) = i;
417 }
418
419 double l = this->giveInternalLength();
420 double g_c = this->giveCriticalEnergy();
421 double kp = this->givePenaltyParameter();
422 double Delta_t = tStep->giveTimeIncrement();
423
424 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
425 temp.zero();
426
427 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
428 IntArray indx_d = computeDamageIndexArray(layer);
429
430 for ( int j = 1; j <= iRuleL->giveNumberOfIntegrationPoints(); j++ ) {
431 GaussPoint *gp = iRuleL->getIntegrationPoint(j - 1);
432 lCoords = *gp->giveCoordinates();
433 this->computeBdmatrixAt(lCoords, Bd); // (2x6)
434 this->computeNdMatrixAt(lCoords, Nd); // (1x6)
435
436 double Gbis = computeGbisInLayer(layer, gp, VM_Total, tStep);
437 double Psibar = this->computeFreeEnergy( gp, tStep );
438 double dV = this->computeVolumeAroundLayer(gp, layer);
439
440
441 // K_dd1 = ( -kp*neg_Mac'(d_dot) / Delta_t + g_c/ l + G'' * Psibar ) * N^t*N;
442 double Delta_d = computeDamageInLayerAt(layer, gp, VM_Incremental, tStep);
443 double factorN = -kp * neg_MaCauleyPrime(Delta_d/Delta_t)/Delta_t + g_c / l + Gbis * Psibar;
444 temp.plusProductSymmUpper(Nd, Nd, factorN * dV);
445
446
447 // K_dd2 = g_c * l * Bd^t * Bd;
448 double factorB = g_c * l;
449 temp.plusProductSymmUpper(Bd, Bd, factorB * dV);
450 }
451 temp.symmetrized();
452 answer.assemble(temp, indx_d, indx_d);
453 }
454
455}
456
457
458void
459Shell7BasePhFi :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
460{
461
462 IntArray loc_u, loc_d;
463 FloatMatrix ansLocal;
464
465 loc_u = this->giveOrdering_Disp();
466 loc_d = this->giveOrdering_Damage();
467
468 int nDofs = this->computeNumberOfDofs();
469 answer.resize( nDofs, nDofs );
470 answer.zero();
471
472 FloatMatrix answer1, answer2, answer3, answer4, answer5, answer6;
473 this->computeStiffnessMatrix_uu(answer1, rMode, tStep);
474 //this->computeStiffnessMatrix_ud(answer2, rMode, tStep); //@todo: check why this gives worse convergence than the numerical
475 this->computeStiffnessMatrix_dd(answer4, rMode, tStep);
476 //this->computeStiffnessMatrixNum_ud(answer5, rMode, tStep);
477 //this->computeStiffnessMatrixNum_dd(answer6, rMode, tStep); // Yields same matrix as the analytical (answer 4)
478
479 //answer2.printYourself();
480 //answer4.printYourself();
481
482
483 //this->computeStiffnessMatrix_du(answer3, rMode, tStep); //symmetric
484 //answer3.beTranspositionOf(answer2);
485 answer3.beTranspositionOf(answer5);
486 //loc_u.printYourself();
487 //loc_d.printYourself();
488
489
490 answer.assemble( answer1, loc_u, loc_u );
491 //answer.assemble( answer5, loc_u, loc_d );
492 //answer.assemble( answer3, loc_d, loc_u );
493 answer.assemble( answer4, loc_d, loc_d );
494
495 //answer4.printYourself();
496
497}
498
499void
500Shell7BasePhFi :: computeStiffnessMatrixNum_ud(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
501{
502
503 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
504 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
505 int ndofs_d = ndofs - ndofs_u;
506 int upD = 1;
507
508 answer.resize( ndofs_u, ndofs_d);
509 answer.zero();
510
511 FloatArray force, force_dist, forceTemp, solVec(42), force_small;
512 IntArray indxDisp;
513
514
515
516 computeSectionalForces(force, tStep, solVec, upD);
517 //force.printYourself();
518
519 int numberOfLayers = this->layeredCS->giveNumberOfLayers();
520 int numdofMans = this->giveNumberOfDofManagers();
521
522 for ( int indx = 1; indx <= numberOfLayers*numdofMans; indx++ ) {
523 computeSectionalForces_dist(force_dist, indx, tStep, solVec, upD);
524 //force.printYourself();
525 //force_dist.printYourself();
526 //solVec.printYourself();
527 //force_dist.printYourself();
528 force_dist.subtract(force);
529
530 //force_dist.printYourself();
531
532 const IntArray &indxDisp = this->giveOrderingPhFi(Displacement);
533
534 force_small.beSubArrayOf(force_dist,indxDisp);
535 //force_small.printYourself();
536 force_small.times(1/disturB);
537 answer.addSubVectorCol(force_small,1,indx);
538 }
539}
540
541void
542Shell7BasePhFi :: computeStiffnessMatrixNum_dd(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
543{
544
545 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
546 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
547 int ndofs_d = ndofs - ndofs_u;
548 int upD = 1;
549
550 answer.resize( ndofs_d, ndofs_d);
551 answer.zero();
552
553 FloatArray force, force_dist, forceTemp, solVec(42), force_small;
554 IntArray indxDam;
555
556
557
558 computeSectionalForces(force, tStep, solVec, upD);
559 //force.printYourself();
560
561 int numberOfLayers = this->layeredCS->giveNumberOfLayers();
562 int numdofMans = this->giveNumberOfDofManagers();
563
564 for ( int indx = 1; indx <= numberOfLayers*numdofMans; indx++ ) {
565 computeSectionalForces_dist(force_dist, indx, tStep, solVec, upD);
566 //force_dist.printYourself();
567 force_dist.subtract(force);
568
569 //force_dist.printYourself();
570
571 const IntArray &indxDam = this->giveOrderingPhFi(Damage);
572
573 force_small.beSubArrayOf(force_dist,indxDam);
574 //force_small.printYourself();
575 force_small.times(1/disturB);
576 answer.addSubVectorCol(force_small,1,indx);
577 }
578}
579
580void
581Shell7BasePhFi :: new_computeBulkTangentMatrix(FloatMatrix &answer, FloatArray &solVec, FloatArray &solVecI, FloatArray &solVecJ, MatResponseMode rMode, TimeStep *tStep)
582{
583 //Shell7Base :: new_computeBulkTangentMatrix(answer, solVec, solVecI, solVecJ,rMode, tStep);
584 //return;
585 //@todo optimize method - remove I,J-stuff since XFEM does not need this anymore
586 FloatMatrix A [ 3 ] [ 3 ], lambdaI [ 3 ], lambdaJ [ 3 ];
587 FloatMatrix L(18,18), B;
588 FloatMatrix K;
589
590 int ndofs = Shell7BasePhFi :: giveNumberOfuDofs();
591 answer.resize(ndofs, ndofs);
592 answer.zero();
593
594 int numberOfLayers = this->layeredCS->giveNumberOfLayers();
595 FloatMatrix temp;
596 FloatArray genEpsI, genEpsJ, genEps, lCoords;
597
598 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
599 IntegrationRule *iRule = integrationRulesArray [ layer - 1 ];
600 StructuralMaterial *mat = static_cast< StructuralMaterial* >( domain->giveMaterial( this->layeredCS->giveLayerMaterial(layer) ) );
601
602 for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
603 GaussPoint *gp = iRule->getIntegrationPoint(i);
604 lCoords = *gp->giveCoordinates();
605
606 this->computeBmatrixAt(lCoords, B, 0, 0);
607 this->computeGeneralizedStrainVectorNew(genEpsI, solVecI, B);
608 this->computeGeneralizedStrainVectorNew(genEpsJ, solVecJ, B);
609 this->computeGeneralizedStrainVectorNew(genEps , solVec , B);
610 // Material stiffness
611 Shell7Base :: computeLinearizedStiffness(gp, mat, tStep, A, genEps);
612
613 double zeta = giveGlobalZcoord( *gp->giveCoordinates() );
614 this->computeLambdaGMatrices(lambdaI, genEpsI, zeta);
615 this->computeLambdaGMatrices(lambdaJ, genEpsJ, zeta);
616
617 // L = sum_{i,j} (lambdaI_i)^T * A^ij * lambdaJ_j
618 // @todo Naive implementation - should be optimized
619 // note: L will only be symmetric if lambdaI = lambdaJ
620 L.zero();
621 for ( int j = 0; j < 3; j++ ) {
622 for ( int k = 0; k < 3; k++ ) {
623 this->computeTripleProduct(temp, lambdaI [ j ], A [ j ][ k ], lambdaJ [ k ]);
624 L.add(temp);
625 }
626 }
627
628 this->computeTripleProduct(K, B, L, B);
629 double dV = this->computeVolumeAroundLayer(gp, layer);
630 double g = this->computeGInLayer(layer, gp, VM_Total, tStep); // Check so that the correct computeDamageAt method is used (in Shell7BasePhFi)
631 answer.add( g*dV, K );
632
633
634 }
635 }
636
637 //@todo Note! Since this block matrix is assembled outside this method with ordering_disp no assembly should be made inside here (=> double assembly) JB
638 //const IntArray &ordering = this->giveOrdering_All();
639 //answer.assemble(tempAnswer, ordering, ordering);
640}
641
642#if 0
643int
644 Shell7BasePhFi :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
645{
646 // Compute special IST quantities this element should support
647 switch (type) {
648 case IST_CauchyStressTensor:
649 this->computeCauchyStressVector(answer, gp, tStep); //@todo: add damage!!!
650 return 1;
651 default:
652 return Element :: giveIPValue(answer, gp, type, tStep);
653 }
654}
655#endif // 0
656
657
658// Internal forces
659
660#if 0
661void
662Shell7BasePhFi :: giveUpdatedSolutionVector_d(FloatArray &answer, TimeStep *tStep)
663{
664
665 IntArray dofIdArray;
666
667 //Shell7Base :: giveDofManDofIDMask(dummy, EID_MomentumBalance, dofIdArray);
668 this->giveDofManDofIDMask_d(dofIdArray);
669
670 this->computeVectorOfDofIDs(dofIdArray, VM_Total, tStep, temp);
671 answer.assemble( temp, this->giveOrdering(AllInv) );
672}
673#endif
674
675//
676
677
678void
679Shell7BasePhFi :: computeSectionalForces(FloatArray &answer, TimeStep *tStep, FloatArray &solVec, int useUpdatedGpRecord)
680{
681
682 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
683 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
684 int ndofs_d = ndofs - ndofs_u;
685
686 int numberOfLayers = this->layeredCS->giveNumberOfLayers(); // conversion of types
687 double sectionalForces_ds = 0;
688 FloatArray fu(ndofs_u), fd(ndofs_d), ftemp_u, ftemp_d(this->giveNumberOfDofManagers()), sectionalForces, sectionalForces_dv;
689 FloatArray genEps, genEpsD, totalSolVec, lCoords, Nd, f_debug(this->giveNumberOfDofManagers()), F_debug(ndofs_d);
690 FloatMatrix B, Bd;
691
692 totalSolVec = solVec;
693 sectionalForces.resize(2);
694 FloatArray dVecTotal, dVecLayer, Ddam_Dxi, gradd;
695 this->computeDamageUnknowns(dVecTotal, VM_Total, tStep); // vector with all damage nodal values for this element
696
697 fu.zero();
698 fd.zero();
699
700 FloatMatrix K_dd;
701 this->computeStiffnessMatrix_dd(K_dd, TangentStiffness, tStep);
702 F_debug.beProductOf(K_dd, dVecTotal);
703
704 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
705 ftemp_d.zero();
706 f_debug.zero();
707 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
708 Material *mat = domain->giveMaterial( this->layeredCS->giveLayerMaterial(layer) );
709 IntArray indx_d = computeDamageIndexArray(layer);
710 dVecLayer.beSubArrayOf(dVecTotal, indx_d); // vector of damage variables for a given layer
711
712 for ( int j = 1; j <= iRuleL->giveNumberOfIntegrationPoints(); j++ ) {
713 GaussPoint *gp = iRuleL->getIntegrationPoint(j - 1);
714 lCoords = *gp->giveCoordinates();
715 this->computeBmatrixAt(lCoords, B);
716 this->computeBdmatrixAt(lCoords, Bd); // (2x6)
717 this->computeNdvectorAt(lCoords, Nd); // (1x6)
718 gradd.beProductOf(Bd,dVecLayer); // [dalpha/dX1, dalpha/dX2]
719
720 this->computeGeneralizedStrainVectorNew(genEps, totalSolVec, B); // used for computing the stress
721 double zeta = giveGlobalZcoord( *gp->giveCoordinates() );
722
723 // Computation of sectional forces: fu = Bu^t*[N M T Ms Ts]^t
724 this->computeSectionalForcesAt(sectionalForces, gp, mat, tStep, genEps, genEps, zeta); // these are per unit volume
725 ftemp_u.beTProductOf(B,sectionalForces);
726 double dV = this->computeVolumeAroundLayer(gp, layer);
727 double g = this->computeGInLayer(layer, gp, VM_Total, tStep);
728 fu.add(dV*g, ftemp_u);
729
730 // Computation of sectional forces: fd = Nd^t * sectionalForces_ds + Bd^t * sectionalForces_dv
731 this->computeSectionalForcesAt_d(sectionalForces_ds, sectionalForces_dv, gp, mat, tStep, zeta, layer, gradd); // these are per unit volume
732
733 // "external force" for linear problem -(-2)*psibar
734 double Psibar = this->computeFreeEnergy( gp, tStep );
735 f_debug.add( -2.0*Psibar*dV, Nd );
736
737 ftemp_d.add( sectionalForces_ds*dV, Nd );
738 ftemp_d.plusProduct(Bd, sectionalForces_dv, dV);
739
740 }
741 // Assemble layer contribution into the correct place
742 F_debug.assemble(f_debug, indx_d);
743 fd.assemble(ftemp_d, indx_d);
744 }
745 //F_debug.printYourself();
746 //fd.printYourself();
747 answer.resize( ndofs );
748 answer.zero();
749 const IntArray &ordering_disp = this->giveOrderingPhFi(Displacement);
750 const IntArray &ordering_damage = this->giveOrderingPhFi(Damage);
751 answer.assemble(fu, ordering_disp); //ordering_disp contains only displacement related dofs, not damage
752 answer.assemble(fd, ordering_damage);
753 //answer.assemble(F_debug, ordering_damage);
754
755}
756
757void
758Shell7BasePhFi :: computeSectionalForces_dist(FloatArray &answer, int indx, TimeStep *tStep, FloatArray &solVec, int useUpdatedGpRecord)
759{
760
761 int ndofs = Shell7BasePhFi :: giveNumberOfDofs();
762 int ndofs_u = Shell7BasePhFi :: giveNumberOfuDofs();
763 int ndofs_d = ndofs - ndofs_u;
764
765 int numberOfLayers = this->layeredCS->giveNumberOfLayers(); // conversion of types
766 double sectionalForces_ds = 0;
767 FloatArray fu(ndofs_u), fd(ndofs_d), ftemp_u, ftemp_d(this->giveNumberOfDofManagers()), sectionalForces, sectionalForces_dv;
768 FloatArray genEps, genEpsD, totalSolVec, lCoords, Nd;
769 FloatMatrix B, Bd;
770 this->giveUpdatedSolutionVector(totalSolVec, tStep); // => x, m, gam
771
772 fu.zero();
773 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
774 ftemp_d.zero();
775 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
776 Material *mat = domain->giveMaterial( this->layeredCS->giveLayerMaterial(layer) );
777 IntArray indx_d = computeDamageIndexArray(layer);
778 FloatArray dVecTotal, dVecLayer, Ddam_Dxi, gradd;
779
780 this->computeDamageUnknowns(dVecTotal, VM_Total, tStep); // vector with all damage nodal values for this element
781 dVecTotal.at(indx) = dVecTotal.at(indx) + disturB;
782 dVecLayer.beSubArrayOf(dVecTotal, indx_d); // vector of damage variables for a given layer
783
784 //dVecLayer.printYourself();
785 for ( int j = 1; j <= iRuleL->giveNumberOfIntegrationPoints(); j++ ) {
786 GaussPoint *gp = iRuleL->getIntegrationPoint(j - 1);
787 lCoords = *gp->giveCoordinates();
788 this->computeBmatrixAt(lCoords, B);
789 this->computeBdmatrixAt(lCoords, Bd); // (2x6)
790 this->computeNdvectorAt(lCoords, Nd); // (1x6)
791 //Ddam_Dxi.beProductOf(Bd,dVecLayer); // [dalpha/dxi1, dalpha/dxi2] OLD
792 gradd.beProductOf(Bd,dVecLayer); // [dalpha/dX1, dalpha/dX2]
793
794 this->computeGeneralizedStrainVectorNew(genEpsD, solVec, B);
795 this->computeGeneralizedStrainVectorNew(genEps, totalSolVec, B); // used for computing the stress
796
797 double zeta = giveGlobalZcoord( *gp->giveCoordinates() );
798
799 // Computation of sectional forces: fu = Bu^t*[N M T Ms Ts]^t
800// this->computeSectionalForcesAt(sectionalForces, gp, mat, tStep, genEps, genEpsD, zeta); // these are per unit volume
801 this->computeSectionalForcesAt(sectionalForces, gp, mat, tStep, genEps, genEps, zeta); // these are per unit volume
802 ftemp_u.beTProductOf(B,sectionalForces);
803 double dV = this->computeVolumeAroundLayer(gp, layer);
804 double g = this->computeGInLayer_dist(layer, indx, gp, VM_Total, tStep);
805 fu.add(dV*g, ftemp_u);
806
807 // Computation of sectional forces: fd = Nd^t * sectionalForces_ds + Bd^t * sectionalForces_dv
808 //this->computeSectionalForcesAt_d(sectionalForces_ds, sectionalForces_dv, gp, mat, tStep, zeta, layer, Ddam_Dxi); // these are per unit volume // OLD
809 this->computeSectionalForcesAt_d_dist(sectionalForces_ds, sectionalForces_dv, indx, gp, mat, tStep, zeta, layer, gradd); // these are per unit volume
810
811 ftemp_d.add( sectionalForces_ds*dV, Nd );
812 //ftemp_d.plusProduct(BdX, sectionalForces_dv, dV); // OLD
813 ftemp_d.plusProduct(Bd, sectionalForces_dv, dV);
814
815 }
816 // Assemble layer contribution into the correct place
817 fd.assemble(ftemp_d, indx_d);
818
819 }
820
821 answer.resize( ndofs );
822 answer.zero();
823 const IntArray &ordering_disp = this->giveOrderingPhFi(Displacement);
824 const IntArray &ordering_damage = this->giveOrderingPhFi(Damage);
825 answer.assemble(fu, ordering_disp); //ordering_disp contains only displacement related dofs, not damage
826 answer.assemble(fd, ordering_damage);
827
828}
829
830double
831Shell7BasePhFi :: neg_MaCauley(double par)
832{
833 return 0.5 * ( abs(par) - par );
834}
835
836double
837Shell7BasePhFi :: neg_MaCauleyPrime(double par)
838{
839 return 0.5 * ( abs(par)/(par + 1.0e-12) - 1.0 ); // 0.5*(sign - 1) taken from Ragnars code
840}
841
842void
843Shell7BasePhFi :: computeSectionalForcesAt_d(double &sectionalForcesScal, FloatArray &sectionalForcesVec, IntegrationPoint *gp,
844 Material *mat, TimeStep *tStep, double zeta, int layer, FloatArray &gradd)
845{
846 //PhaseFieldCrossSection *cs = static_cast... // in the future...
847 sectionalForcesVec.zero();
848 //sectionalForcesScal = -kp*neg_Mac(alpha_dot) + g_c/l*d + G'*Psibar
849 double kp = this->givePenaltyParameter();
850 double Delta_t = tStep->giveTimeIncrement();
851 double d = computeDamageInLayerAt(layer, gp, VM_Total, tStep);
852 double Delta_d = computeDamageInLayerAt(layer, gp, VM_Incremental, tStep);
853 double l = this->giveInternalLength();
854 double g_c = this->giveCriticalEnergy();
855 double Gprim = computeGprimInLayer(layer, gp, VM_Total, tStep);
856 double Psibar = this->computeFreeEnergy( gp, tStep );
857
858 if (Psibar < 0) {
859 OOFEM_ERROR1("Shell7BasePhFi :: computeSectionalForcesAt_d - negative strain energy predicted")
860 }
861
862 sectionalForcesScal = -kp * neg_MaCauley(Delta_d/Delta_t) + g_c / l * d + Gprim * Psibar;
863
864 //sectionalForcesVec = grad(alpha) * g_c * l
865 sectionalForcesVec.beScaled(g_c*l,gradd);
866 //sectionalForcesVec = gradd * g_c * l;
867 //printf("scal %e vec %e %e\n", sectionalForcesScal, sectionalForcesVec.at(1), sectionalForcesVec.at(2));
868}
869
870void
871Shell7BasePhFi :: computeSectionalForcesAt_d_dist(double &sectionalForcesScal, FloatArray &sectionalForcesVec, int indx, IntegrationPoint *gp,
872 Material *mat, TimeStep *tStep, double zeta, int layer, FloatArray &gradd)
873{
874 //PhaseFieldCrossSection *cs = static_cast... // in the future...
875
876 //sectionalForcesScal = -kp*neg_Mac(alpha_dot) + g_c/l*d + G'*Psibar
877 double kp = this->givePenaltyParameter();
878 double Delta_t = tStep->giveTimeIncrement();
879 double d = computeDamageInLayerAt_dist(layer, indx, gp, VM_Total, tStep);
880 double Delta_d = computeDamageInLayerAt_dist(layer, indx, gp, VM_Incremental, tStep);
881 double l = this->giveInternalLength();
882 double g_c = this->giveCriticalEnergy();
883 double Gprim = computeGprimInLayer_dist(layer, indx, gp, VM_Total, tStep);
884 double Psibar = this->computeFreeEnergy( gp, tStep );
885
886 sectionalForcesScal = -kp * neg_MaCauley(Delta_d/Delta_t) + g_c / l * d + Gprim * Psibar;
887 sectionalForcesVec.zero();
888 //sectionalForcesVec = grad(alpha) * g_c * l
889 sectionalForcesVec.beScaled(g_c*l,gradd);
890 //sectionalForcesVec = gradd * g_c * l;
891
892}
893
894#if 0
895void
896Shell7BasePhFi :: computeBdmatrixAt(FloatArray &lCoords, FloatMatrix &answer)
897{
898 // Returns the matrix {B} of the receiver, evaluated at aGaussPoint. Such that
899 // B*a = [Dalpha/Dxi1, Dalpha/Dxi2]
900
901 FloatMatrix dNdxi;
902 this->fei->evaldNdxi( dNdxi, lCoords, FEIElementGeometryWrapper(this) );
903 answer.beTranspositionOf(dNdxi);
904
905
906}
907
908#else
909
910void
911Shell7BasePhFi :: computeBdmatrixAt(FloatArray &lCoords, FloatMatrix &answer)
912{
913 // Returns the matrix {B} of the receiver, evaluated at aGaussPoint. Such that
914 // B*a = [Dalpha/DX1, Dalpha/DX2]
915
916 FloatMatrix dNdxi;
917 this->fei->evaldNdxi( dNdxi, lCoords, FEIElementGeometryWrapper(this) );
918
919 FloatMatrix Gcon, Gcon_red, temp;
920 Shell7Base :: evalInitialContravarBaseVectorsAt(lCoords, Gcon); //[G^1 G^2 G^3]
921 Gcon_red.beSubMatrixOf(Gcon,1,2,1,2);
922
923 temp.beTranspositionOf(dNdxi);
924 answer.beProductOf(Gcon_red, temp); // dN/dX = [G^1 G^2] * dN/dxi
925
926}
927
928#endif
929
930void
931Shell7BasePhFi :: computeNdvectorAt(const FloatArray &iLocCoords, FloatArray &answer)
932{
933 // Returns the matrix {N} of the receiver, evaluated at aGaussPoint. Such that
934 // N*a = alpha
935 this->fei->evalN( answer, iLocCoords, FEIElementGeometryWrapper(this) );
936}
937
938
939void
940Shell7BasePhFi :: computeNdMatrixAt(const FloatArray &iLocCoords, FloatMatrix &answer)
941{
942 FloatArray Nvec;
943 this->computeNdvectorAt(iLocCoords, Nvec);
944 answer.resize(1,Nvec.giveSize());
945 for ( int i = 1; i<= Nvec.giveSize(); i++ ){
946 answer.at(1,i) = Nvec.at(i);
947 }
948
949}
950
951
952// Computation of solution vectors
953
954
955void
956Shell7BasePhFi :: computeVectorOfDofIDs(const IntArray &dofIdArray, ValueModeType u, TimeStep *stepN, FloatArray &answer)
957{
958 // Routine to extract the solution vector for an element given an dofid array.
959 // Size will be numberOfDofs and if a certain dofId does not exist a zero is used as value.
960 // This method may e.g. be used to obtain the enriched part of the solution vector
962 answer.resize( Shell7BasePhFi ::giveNumberOfDofs());
963 answer.zero();
964 int k = 0;
965 for ( int i = 1; i <= numberOfDofMans; i++ ) {
966 DofManager *dMan = this->giveDofManager(i);
967 for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {
968
969 if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
970 Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
973 answer.at(k+j) = d->giveUnknown(u, stepN); // Martin: modification to be used also for rates
974 }
975 }
976 k += 7;
977 }
978}
979
980
981
982
983
984// VTK export
985
986void
987Shell7BasePhFi :: giveShellExportData(ExportRegion &vtkPiece, IntArray &primaryVarsToExport, IntArray &internalVarsToExport, IntArray cellVarsToExport, TimeStep *tStep )
988{
989
990 int numCells = this->layeredCS->giveNumberOfLayers();
991 const int numCellNodes = 15; // quadratic wedge
992 int numTotalNodes = numCellNodes*numCells;
993
994 vtkPiece.setNumberOfCells(numCells);
995 vtkPiece.setNumberOfNodes(numTotalNodes);
996
997 std::vector <FloatArray> nodeCoords;
998 int val = 1;
999 int offset = 0;
1000 IntArray nodes(numCellNodes);
1001
1002 // Compute fictious node coords
1003 int nodeNum = 1;
1004 for ( int layer = 1; layer <= numCells; layer++ ) {
1005
1006
1007 // Node coordinates
1008 this->giveFictiousNodeCoordsForExport(nodeCoords, layer);
1009
1010 for ( int node = 1; node <= numCellNodes; node++ ) {
1011 vtkPiece.setNodeCoords(nodeNum, nodeCoords[node-1] );
1012 nodeNum += 1;
1013 }
1014
1015 // Connectivity
1016 for ( int i = 1; i <= numCellNodes; i++ ) {
1017 nodes.at(i) = val++;
1018 }
1019 vtkPiece.setConnectivity(layer, nodes);
1020
1021 // Offset
1022 offset += numCellNodes;
1023 vtkPiece.setOffset(layer, offset);
1024
1025 // Cell types
1026 vtkPiece.setCellType(layer, 26); // Quadratic wedge
1027 }
1028
1029
1030 // Export nodal variables from primary fields
1031 vtkPiece.setNumberOfPrimaryVarsToExport(primaryVarsToExport, numTotalNodes);
1032
1033 std::vector<FloatArray> updatedNodeCoords;
1034 FloatArray u(3), damage(1);
1035 std::vector<FloatArray> values;
1036 FloatArray dVecTotal, dVecLayer, damageArray;
1037 this->computeDamageUnknowns(dVecTotal, VM_Total, tStep); // vector with all damage nodal values for this element
1038
1039 for ( int fieldNum = 1; fieldNum <= primaryVarsToExport.giveSize(); fieldNum++ ) {
1040 UnknownType type = ( UnknownType ) primaryVarsToExport.at(fieldNum);
1041 nodeNum = 1;
1042 for ( int layer = 1; layer <= numCells; layer++ ) {
1043
1044 if ( type == DisplacementVector ) { // compute displacement as u = x - X
1045 this->giveFictiousNodeCoordsForExport(nodeCoords, layer);
1046 this->giveFictiousUpdatedNodeCoordsForExport(updatedNodeCoords, layer, tStep);
1047 for ( int j = 1; j <= numCellNodes; j++ ) {
1048 u = updatedNodeCoords[j-1];
1049 u.subtract(nodeCoords[j-1]);
1050 vtkPiece.setPrimaryVarInNode(fieldNum, nodeNum, u);
1051 nodeNum += 1;
1052 }
1053 } else if ( type == ScalarDamage ) {
1054 // compute damage in layer
1055 // set same damage in the nodes above and below the nodes on the midsurface
1056 IntArray indx_d = computeDamageIndexArray(layer);
1057
1058 dVecLayer.beSubArrayOf(dVecTotal, indx_d); // vector of damage variables for a given layer
1059 damageArray.setValues(15, dVecLayer.at(1), dVecLayer.at(2), dVecLayer.at(3), dVecLayer.at(1), dVecLayer.at(2), dVecLayer.at(3),
1060 dVecLayer.at(4), dVecLayer.at(5), dVecLayer.at(6), dVecLayer.at(4), dVecLayer.at(5), dVecLayer.at(6),
1061 dVecLayer.at(1), dVecLayer.at(2), dVecLayer.at(3) );
1062 for ( int j = 1; j <= numCellNodes; j++ ) {
1063 damage.at(1) = damageArray.at(j);
1064 vtkPiece.setPrimaryVarInNode(fieldNum, nodeNum, damage);
1065 nodeNum += 1;
1066 }
1067 } else {
1068 NodalRecoveryMI_recoverValues(values, layer, ( InternalStateType ) 1, tStep); // does not work well - fix
1069 for ( int j = 1; j <= numCellNodes; j++ ) {
1070 vtkPiece.setPrimaryVarInNode(fieldNum, nodeNum, values[j-1]);
1071 nodeNum += 1;
1072 }
1073 }
1074 }
1075 }
1076
1077 // Export nodal variables from internal fields
1078
1079 vtkPiece.setNumberOfInternalVarsToExport( internalVarsToExport, numTotalNodes );
1080 for ( int fieldNum = 1; fieldNum <= internalVarsToExport.giveSize(); fieldNum++ ) {
1081 InternalStateType type = ( InternalStateType ) internalVarsToExport.at(fieldNum);
1082 nodeNum = 1;
1083 //this->recoverShearStress(tStep);
1084 for ( int layer = 1; layer <= numCells; layer++ ) {
1085 recoverValuesFromIP(values, layer, type, tStep);
1086 for ( int j = 1; j <= numCellNodes; j++ ) {
1087 vtkPiece.setInternalVarInNode( fieldNum, nodeNum, values[j-1] );
1088 //ZZNodalRecoveryMI_recoverValues(el.nodeVars[fieldNum], layer, type, tStep);
1089 nodeNum += 1;
1090 }
1091 }
1092 }
1093
1094
1095 // Export cell variables
1096 FloatArray average;
1097 vtkPiece.setNumberOfCellVarsToExport(cellVarsToExport, numCells);
1098 for ( int i = 1; i <= cellVarsToExport.giveSize(); i++ ) {
1099 InternalStateType type = ( InternalStateType ) cellVarsToExport.at(i);;
1100
1101 for ( int layer = 1; layer <= numCells; layer++ ) {
1102 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
1103 VTKXMLExportModule::computeIPAverage(average, iRuleL, this, type, tStep);
1104
1105 if ( average.giveSize() == 6 ) {
1106 vtkPiece.setCellVar(i, layer, convV6ToV9Stress(average) );
1107 } else {
1108 vtkPiece.setCellVar(i, layer, average );
1109 }
1110
1111 }
1112
1113 }
1114
1115
1116
1117
1118}
1119
1120#if 0
1121void
1122Shell7BasePhFi :: recoverValuesFromIP(std::vector<FloatArray> &recoveredValues, int layer, InternalStateType type, TimeStep *tStep)
1123{
1124 // recover nodal values by coosing the ip closest to the node
1125
1126 //FEInterpolation *interpol = static_cast< FEInterpolation * >( &this->interpolationForExport );
1127
1128 // composite element interpolator
1129 FloatMatrix localNodeCoords;
1130 this->interpolationForExport.giveLocalNodeCoords(localNodeCoords);
1131
1132 int numNodes = localNodeCoords.giveNumberOfColumns();
1133 recoveredValues.resize(numNodes);
1134
1135 IntegrationRule *iRule = integrationRulesArray [ layer - 1 ];
1136 IntegrationPoint *ip;
1137
1138 // Find closest ip to the nodes
1139 IntArray closestIPArray(numNodes);
1140 FloatArray nodeCoords, ipCoords, ipValues;
1141
1142 for ( int i = 1; i <= numNodes; i++ ) {
1143 nodeCoords.beColumnOf(localNodeCoords, i);
1144 double distOld = 3.0; // should not be larger
1145 for ( int j = 0; j < iRule->giveNumberOfIntegrationPoints(); j++ ) {
1146 ip = iRule->getIntegrationPoint(j);
1147 ipCoords = *ip->giveCoordinates();
1148 double dist = distance(nodeCoords, ipCoords);
1149 if ( dist < distOld ) {
1150 closestIPArray.at(i) = j;
1151 distOld = dist;
1152 }
1153 }
1154 }
1155
1156 InternalStateValueType valueType = giveInternalStateValueType(type);
1157
1158 // recover ip values
1159 for ( int i = 1; i <= numNodes; i++ ) {
1160 ip = iRule->getIntegrationPoint( closestIPArray.at(i) );
1161 this->giveIPValue(ipValues, ip, type, tStep);
1162 if ( valueType == ISVT_TENSOR_S3 ) {
1163 recoveredValues[i-1].resize(9);
1164 recoveredValues[i-1] = convV6ToV9Stress(ipValues);
1165 } else if ( ipValues.giveSize() == 0 && type == IST_AbaqusStateVector) {
1166 recoveredValues[i-1].resize(23);
1167 recoveredValues[i-1].zero();
1168 } else if ( ipValues.giveSize() == 0 ) {
1169 recoveredValues[i-1].resize(giveInternalStateTypeSize(valueType));
1170 recoveredValues[i-1].zero();
1171
1172 } else {
1173 recoveredValues[i-1] = ipValues;
1174 }
1175 }
1176
1177}
1178
1179
1180void
1181Shell7BasePhFi :: recoverShearStress(TimeStep *tStep)
1182{
1183 // Recover shear stresses at ip by numerical integration of the momentum balance through the thickness
1184 std::vector<FloatArray> recoveredValues;
1185 int numberOfLayers = this->layeredCS->giveNumberOfLayers(); // conversion of types
1186 IntegrationRule *iRuleThickness = specialIntegrationRulesArray[ 0 ];
1187 FloatArray dS, Sold;
1188 FloatMatrix B, Smat(2,6); // 2 stress components * num of in plane ip ///@todo generalize
1189 Smat.zero();
1190 FloatArray Tcon(6), Trec(6); Tcon.zero(); Trec.zero();
1191
1192 for ( int layer = 1; layer <= numberOfLayers; layer++ ) {
1193 IntegrationRule *iRuleL = integrationRulesArray [ layer - 1 ];
1194 this->recoverValuesFromIP(recoveredValues, layer, IST_StressTensor, tStep);
1195 //this->ZZNodalRecoveryMI_recoverValues(recoveredValues, layer, IST_StressTensor, tStep);
1196 double thickness = this->layeredCS->giveLayerThickness(layer);
1197
1198 //set up vector of stresses in the ip's = [S1_xx, S1_yy, S1_xy, ..., Sn_xx, Sn_yy, Sn_xy]
1199 int numNodes = 15;
1200 FloatArray aS(numNodes*3);
1201 for ( int j = 1, pos = 0; j <= numNodes; j++, pos+=3 ) {
1202 aS.at(pos + 1) = recoveredValues[j-1].at(1); // S_xx
1203 aS.at(pos + 2) = recoveredValues[j-1].at(2); // S_yy
1204 aS.at(pos + 3) = recoveredValues[j-1].at(6); // S_xy
1205 }
1206 int numInPlaneIP = 6;
1207
1208 for ( int i = 0; i < iRuleThickness->giveNumberOfIntegrationPoints(); i++ ) {
1209 double dz = thickness * iRuleThickness->getIntegrationPoint(i)->giveWeight();
1210
1211 for ( int j = 0; j < numInPlaneIP; j++ ) {
1212
1213 int point = i*numInPlaneIP + j; // integration point number
1214 GaussPoint *gp = iRuleL->getIntegrationPoint(point);
1215
1216 this->computeBmatrixForStressRecAt(*gp->giveCoordinates(), B, layer);
1217 dS.beProductOf(B,aS*(-dz)); // stress increment
1218
1219 StructuralMaterialStatus* status = dynamic_cast< StructuralMaterialStatus* > ( gp->giveMaterialStatus() );
1220 Sold = status->giveStressVector();
1221
1222 Smat.at(1,j+1) += dS.at(1); // add increment from each level
1223 Smat.at(2,j+1) += dS.at(2);
1224
1225 //Tcon.at(j+1) += Sold.at(5)*dz;
1226
1227 // Replace old stresses with - this should probably not be done as it may affect the convergence in a nonlinear case
1228 Sold.at(5) = Smat.at(1,j+1); // S_xz
1229 Sold.at(4) = Smat.at(2,j+1); // S_yz
1230
1231
1232 status->letStressVectorBe(Sold);
1233 //Trec.at(j+1) += Sold.at(5)*dz;
1234 }
1235 }
1236
1237
1238 }
1239
1240}
1241
1242
1243void
1244Shell7BasePhFi :: computeBmatrixForStressRecAt(FloatArray &lcoords, FloatMatrix &answer, int layer)
1245{
1246 // Returns the special matrix {B} of the receiver, evaluated at aGaussPoint. Such that
1247 // B*a = [dS_xx/dx + dS_xy/dy, dS_yx/dx + dS_yy/dy ]^T, where a is the vector of in plane
1248 // stresses [S_xx, S_yy, S_xy]
1249
1250 // set up virtual cell geometry for an qwedge
1251 const int numNodes = 15;
1252 std::vector<FloatArray> nodes;
1253 giveFictiousNodeCoordsForExport(nodes, layer);
1254
1255 int VTKWedge2EL [] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };
1256 FloatArray *coords[numNodes];
1257
1258
1259 for ( int i = 1; i <= numNodes; i++ ) {
1260 int pos = VTKWedge2EL[ i-1 ];
1261 coords[ i - 1 ] = &nodes[ pos - 1];
1262
1263 }
1264
1265 FEInterpolation *interpol = static_cast< FEInterpolation * >( &this->interpolationForExport );
1266 FloatMatrix dNdx;
1267 interpol->evaldNdx( dNdx, lcoords, FEIVertexListGeometryWrapper(numNodes, (const FloatArray **)coords ) );
1268
1269 /*
1270 * 1 [d/dx 0 d/dy
1271 * 1 0 d/dy d/dx]
1272 */
1273 int ndofs = numNodes*3;
1274 answer.resize(2, ndofs);
1275 for ( int i = 1, j = 0; i <= numNodes; i++, j += 3 ) {
1276 answer.at(1, j + 1) = dNdx.at(i, 1);
1277 answer.at(1, j + 3) = dNdx.at(i, 2);
1278 answer.at(2, j + 2) = dNdx.at(i, 2);
1279 answer.at(2, j + 3) = dNdx.at(i, 1);
1280 }
1281
1282}
1283
1284
1285std::vector<FloatArray>
1286Shell7BasePhFi :: giveFictiousNodeCoordsForExport(int layer)
1287{
1288 // compute fictious node coords
1289 FloatMatrix localNodeCoords;
1290 this->interpolationForExport.giveLocalNodeCoords(localNodeCoords);
1291
1292 nodes.resize(localNodeCoords.giveNumberOfColumns());
1293 for ( int i = 1; i <= localNodeCoords.giveNumberOfColumns(); i++ ){
1294 FloatArray localCoords(3);
1295 localCoords.beColumnOf(localNodeCoords,i);
1296
1297 nodes[i-1] = this->vtkEvalInitialGlobalCoordinateAt(localCoords, layer);
1298 }
1299 return nodes;
1300}
1301
1302
1303std::vector<FloatArray>
1304Shell7BasePhFi :: giveFictiousCZNodeCoordsForExport(int interface)
1305{
1306 // compute fictious node coords
1307 FloatMatrix localNodeCoords;
1308 this->interpolationForCZExport.giveLocalNodeCoords(localNodeCoords);
1309
1310 std::vector<FloatArray> nodes(localNodeCoords.giveNumberOfColumns());
1311 for ( int i = 1; i <= localNodeCoords.giveNumberOfColumns(); i++ ){
1312 FloatArray localCoords(3);
1313 localCoords.beColumnOf(localNodeCoords,i);
1314
1315 localCoords.at(3) = 1.0;
1316 nodes[i-1] = this->vtkEvalInitialGlobalCoordinateAt(localCoords, interface);
1317 }
1318 return nodes;
1319}
1320
1321void
1322Shell7BasePhFi :: giveFictiousUpdatedNodeCoordsForExport(int layer, TimeStep *tStep)
1323{
1324 // compute fictious node coords
1325
1326 FloatMatrix localNodeCoords;
1327 this->interpolationForExport.giveLocalNodeCoords(localNodeCoords);
1328 std::vector<FloatArray> nodes(localNodeCoords.giveNumberOfColumns());
1329 for ( int i = 1; i <= localNodeCoords.giveNumberOfColumns(); i++ ){
1330 FloatArray localCoords;
1331 localCoords.beColumnOf(localNodeCoords,i);
1332
1333 nodes[i-1] = this->vtkEvalUpdatedGlobalCoordinateAt(localCoords, layer, tStep);
1334 }
1335 return nodes;
1336}
1337
1338
1339//void
1340//Shell7BasePhFi :: giveLocalNodeCoordsForExport(FloatArray &nodeLocalXi1Coords, FloatArray &nodeLocalXi2Coords, FloatArray &nodeLocalXi3Coords) {
1341// // Local coords for a quadratic wedge element (VTK cell type 26)
1342// double z = 0.999;
1343// nodeLocalXi1Coords.setValues(15, 1., 0., 0., 1., 0., 0., .5, 0., .5, .5, 0., .5, 1., 0., 0.);
1344// nodeLocalXi2Coords.setValues(15, 0., 1., 0., 0., 1., 0., .5, .5, 0., .5, .5, 0., 0., 1., 0.);
1345// nodeLocalXi3Coords.setValues(15, -z, -z, -z, z, z, z, -z, -z, -z, z, z, z, 0., 0., 0.);
1346//}
1347
1348#endif
1349
1350
1351
1352// Misc functions
1353
1354
1355
1356} // end namespace oofem
oofem::IntegrationPoint
GaussPoint IntegrationPoint
Definition gausspoint.h:272
oofem::giveInternalStateTypeSize
int giveInternalStateTypeSize(InternalStateValueType valType)
Definition cltypes.C:229
oofem::disturB
const double disturB
Definition shell7basePhFi.C:59
oofem::InternalStateValueType
InternalStateValueType
Determines the type of internal variable.
Definition internalstatevaluetype.h:40
oofem::giveInternalStateValueType
InternalStateValueType giveInternalStateValueType(InternalStateType type)
Definition cltypes.C:75
oofem::distance
double distance(const FloatArray &x, const FloatArray &y)
Definition floatarray.C:1016
oofem::nLayers
const int nLayers
Definition shell7basePhFi.C:58
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