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

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011