OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mitc4.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 "../sm/Elements/Shells/mitc4.h"
36 #include "../sm/Materials/structuralms.h"
37 #include "../sm/Materials/structuralmaterial.h"
38 #include "../sm/CrossSections/structuralcrosssection.h"
39 #include "fei2dquadlin.h"
40 #include "node.h"
41 #include "material.h"
42 #include "crosssection.h"
43 #include "gausspoint.h"
44 #include "../sm/CrossSections/variablecrosssection.h"
45 #include "gaussintegrationrule.h"
46 #include "floatmatrix.h"
47 #include "floatarray.h"
48 #include "intarray.h"
49 #include "load.h"
50 #include "boundaryload.h"
51 #include "mathfem.h"
52 #include "classfactory.h"
53 #include "connectivitytable.h"
54 
55 namespace oofem {
56 REGISTER_Element(MITC4Shell);
57 
58 FEI2dQuadLin MITC4Shell :: interp_lin(1, 2);
59 
60 MITC4Shell :: MITC4Shell(int n, Domain *aDomain) :
63 {
64  numberOfDofMans = 4;
65  nPointsXY = 4;
66  nPointsZ = 2;
67 
69 }
70 
71 
74 
75 
78 {
79  return & interp_lin;
80 }
81 
82 
83 Interface *
85 {
86  if ( interface == ZZNodalRecoveryModelInterfaceType ) {
87  return static_cast< ZZNodalRecoveryModelInterface * >(this);
88  } else if ( interface == SPRNodalRecoveryModelInterfaceType ) {
89  return static_cast< SPRNodalRecoveryModelInterface * >(this);
90  } else if ( interface == NodalAveragingRecoveryModelInterfaceType ) {
91  return static_cast< NodalAveragingRecoveryModelInterface * >(this);
92  } else if ( interface == SpatialLocalizerInterfaceType ) {
93  return static_cast< SpatialLocalizerInterface * >(this);
94  }
95 
96  return NULL;
97 }
98 
99 
100 void
102 {
103  pap.resize(numberOfDofMans);
104  for ( int i = 1; i <= numberOfDofMans; i++ ) {
105  pap.at(i) = this->giveNode(i)->giveNumber();
106  }
107 }
108 
109 
110 void
112 {
113  int found = 0;
114  answer.resize(1);
115 
116  for ( int i = 1; i <= numberOfDofMans; i++ ) {
117  if ( this->giveNode(i)->giveNumber() == pap ) {
118  found = 1;
119  }
120  }
121 
122  if ( found ) {
123  answer.at(1) = pap;
124  } else {
125  OOFEM_ERROR("unknown node number %d", pap);
126  }
127 }
128 
129 
130 int
132 {
134 }
135 
136 
139 {
140  return SPRPatchType_3dBiLin;
141 }
142 
143 
144 void
146 // Sets up the array containing the eight Gauss points of the receiver.
147 {
148  if ( integrationRulesArray.size() == 0 ) {
149  integrationRulesArray.resize(1);
150  integrationRulesArray [ 0 ].reset( new GaussIntegrationRule(1, this, 1, 10) );
152  }
153 }
154 
155 
156 void
158 {
159  if ( directorType == 0 ) { // normal to the midplane
160  FloatArray e1, e2, e3;
161  this->computeLocalBaseVectors(e1, e2, e3);
162  V1 = e3;
163  V2 = e3;
164  V3 = e3;
165  V4 = e3;
166  } else if ( directorType == 1 ) { // nodal average
167  FloatArray e1, e2, e3;
168  FloatMatrix directors;
169  FloatArray nodeDir;
170  nodeDir.resize(3);
171  directors.resize(3, 4);
172  int csNum = this->giveCrossSection()->giveNumber();
173  // OOFEM_ERROR("DirectorType 1 not implemented yet.");
174  IntArray neighbours;
175  ConnectivityTable *conTable = this->giveDomain()->giveConnectivityTable();
176  IntArray node(1);
177  int count;
178 
179  for ( int i = 1; i <= 4; i++ ) {
180  count = 0;
181  nodeDir.zero();
182  node.at(1) = this->giveNode(i)->giveNumber();
183  conTable->giveNodeNeighbourList(neighbours, node);
184 
185  for ( int j = 1; j <= neighbours.giveSize(); j++ ) {
186  MITC4Shell *neighbour = dynamic_cast< MITC4Shell * >( this->giveDomain()->giveElement( neighbours.at(j) ) );
187  if ( neighbour ) {
188  if ( neighbour->giveCrossSection()->giveNumber() == csNum ) {
189  neighbour->computeLocalBaseVectors(e1, e2, e3);
190  nodeDir.add(e3);
191  count++;
192  }
193  }
194  }
195  directors.at(1, i) = nodeDir.at(1) / count;
196  directors.at(2, i) = nodeDir.at(2) / count;
197  directors.at(3, i) = nodeDir.at(3) / count;
198  }
199  V1.beColumnOf(directors, 1);
200  V2.beColumnOf(directors, 2);
201  V3.beColumnOf(directors, 3);
202  V4.beColumnOf(directors, 4);
203 
204  V1.normalize();
205  V2.normalize();
206  V3.normalize();
207  V4.normalize();
208  } else if ( directorType == 2 ) { // specified at crosssection
209  V1.resize(3);
210  V2.resize(3);
211  V3.resize(3);
212  V4.resize(3);
213  FloatArray *c1, *c2, *c3, *c4;
214  c1 = this->giveNode(1)->giveCoordinates();
215  c2 = this->giveNode(2)->giveCoordinates();
216  c3 = this->giveNode(3)->giveCoordinates();
217  c4 = this->giveNode(4)->giveCoordinates();
218  V1.at(1) = this->giveCrossSection()->give(CS_DirectorVectorX, * c1, this, false);
219  V1.at(2) = this->giveCrossSection()->give(CS_DirectorVectorY, * c1, this, false);
220  V1.at(3) = this->giveCrossSection()->give(CS_DirectorVectorZ, * c1, this, false);
221 
222  V2.at(1) = this->giveCrossSection()->give(CS_DirectorVectorX, * c2, this, false);
223  V2.at(2) = this->giveCrossSection()->give(CS_DirectorVectorY, * c2, this, false);
224  V2.at(3) = this->giveCrossSection()->give(CS_DirectorVectorZ, * c2, this, false);
225 
226  V3.at(1) = this->giveCrossSection()->give(CS_DirectorVectorX, * c3, this, false);
227  V3.at(2) = this->giveCrossSection()->give(CS_DirectorVectorY, * c3, this, false);
228  V3.at(3) = this->giveCrossSection()->give(CS_DirectorVectorZ, * c3, this, false);
229 
230  V4.at(1) = this->giveCrossSection()->give(CS_DirectorVectorX, * c4, this, false);
231  V4.at(2) = this->giveCrossSection()->give(CS_DirectorVectorY, * c4, this, false);
232  V4.at(3) = this->giveCrossSection()->give(CS_DirectorVectorZ, * c4, this, false);
233 
234  V1.normalize();
235  V2.normalize();
236  V3.normalize();
237  V4.normalize();
238  } else {
239  OOFEM_ERROR("Unsupported directorType");
240  }
241 }
242 
243 
244 void
246 {
247  FloatArray V1g, V2g, V3g, V4g;
248  this->giveDirectorVectors(V1g, V2g, V3g, V4g);
250 
255 }
256 
257 
258 void
260 // Returns the [6x24] displacement interpolation matrix {N} of the receiver,
261 // evaluated at gp.
262 // Zeroes in rows 4, 5, 6.
263 {
264  FloatArray h(4);
265 
266  interp_lin.evalN( h, iLocCoord, FEIElementGeometryWrapper(this) );
267 
268  double a1, a2, a3, a4;
269  this->giveThickness(a1, a2, a3, a4);
270 
271  FloatArray V1, V2, V3, V4;
272  this->giveLocalDirectorVectors(V1, V2, V3, V4);
273 
274  FloatArray e2 = {
275  0, 1, 0
276  };
277 
278  FloatArray V11(3), V12(3), V13(3), V14(3), V21(3), V22(3), V23(3), V24(3);
279  V11.beVectorProductOf(e2, V1);
280  V11.normalize();
281  V12.beVectorProductOf(e2, V2);
282  V12.normalize();
283  V13.beVectorProductOf(e2, V3);
284  V13.normalize();
285  V14.beVectorProductOf(e2, V4);
286  V14.normalize();
287 
288  V21.beVectorProductOf(V1, V11);
289  V22.beVectorProductOf(V2, V12);
290  V23.beVectorProductOf(V3, V13);
291  V24.beVectorProductOf(V4, V14);
292 
293  answer.resize(6, 24);
294  answer.zero();
295 
296  answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = h.at(1);
297  answer.at(1, 4) = -iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V21.at(1);
298  answer.at(1, 5) = iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V11.at(1);
299  answer.at(2, 4) = -iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V21.at(2);
300  answer.at(2, 5) = iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V11.at(2);
301  answer.at(3, 4) = -iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V21.at(3);
302  answer.at(3, 5) = iLocCoord.at(3) / 2.0 *a1 *h.at(1) * V11.at(3);
303 
304  answer.at(1, 7) = answer.at(2, 8) = answer.at(3, 9) = h.at(2);
305  answer.at(1, 10) = -iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V22.at(1);
306  answer.at(1, 11) = iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V12.at(1);
307  answer.at(2, 10) = -iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V22.at(2);
308  answer.at(2, 11) = iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V12.at(2);
309  answer.at(3, 10) = -iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V22.at(3);
310  answer.at(3, 11) = iLocCoord.at(3) / 2.0 *a2 *h.at(2) * V12.at(3);
311 
312  answer.at(1, 13) = answer.at(2, 14) = answer.at(3, 15) = h.at(3);
313  answer.at(1, 16) = -iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V23.at(1);
314  answer.at(1, 17) = iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V13.at(1);
315  answer.at(2, 16) = -iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V23.at(2);
316  answer.at(2, 17) = iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V13.at(2);
317  answer.at(3, 16) = -iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V23.at(3);
318  answer.at(3, 17) = iLocCoord.at(3) / 2.0 *a3 *h.at(3) * V13.at(3);
319 
320  answer.at(1, 19) = answer.at(2, 20) = answer.at(3, 21) = h.at(4);
321  answer.at(1, 23) = -iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V24.at(1);
322  answer.at(1, 24) = iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V14.at(1);
323  answer.at(2, 23) = -iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V24.at(2);
324  answer.at(2, 24) = iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V14.at(2);
325  answer.at(3, 23) = -iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V24.at(3);
326  answer.at(3, 24) = iLocCoord.at(3) / 2.0 *a4 *h.at(4) * V14.at(3);
327 }
328 
329 
330 void
332 {
333  StructuralCrossSection *scs = dynamic_cast< StructuralCrossSection * >( this->giveCrossSection() );
334  //SimpleCrossSection *cs = dynamic_cast< SimpleCrossSection * >( this->giveCrossSection() );
335  scs->give3dDegeneratedShellStiffMtrx(answer, rMode, gp, tStep);
336 }
337 
338 
339 void
340 MITC4Shell :: giveNodeCoordinates(double &x1, double &x2, double &x3, double &x4,
341  double &y1, double &y2, double &y3, double &y4,
342  double &z1, double &z2, double &z3, double &z4)
343 {
344  FloatArray nc1(3), nc2(3), nc3(3), nc4(3);
345 
346  this->giveLocalCoordinates( nc1, * ( this->giveNode(1)->giveCoordinates() ) );
347  this->giveLocalCoordinates( nc2, * ( this->giveNode(2)->giveCoordinates() ) );
348  this->giveLocalCoordinates( nc3, * ( this->giveNode(3)->giveCoordinates() ) );
349  this->giveLocalCoordinates( nc4, * ( this->giveNode(4)->giveCoordinates() ) );
350 
351  x1 = nc1.at(1);
352  x2 = nc2.at(1);
353  x3 = nc3.at(1);
354  x4 = nc4.at(1);
355 
356  y1 = nc1.at(2);
357  y2 = nc2.at(2);
358  y3 = nc3.at(2);
359  y4 = nc4.at(2);
360 
361  z1 = nc1.at(3);
362  z2 = nc2.at(3);
363  z3 = nc3.at(3);
364  z4 = nc4.at(3);
365 }
366 
367 void
369 // Returns global coordinates given in global vector
370 // transformed into local coordinate system of the
371 // receiver
372 {
373  FloatArray offset;
374  // test the parametr
375  if ( global.giveSize() != 3 ) {
376  OOFEM_ERROR("cannot transform coordinates - size mismatch");
377  exit(1);
378  }
379 
381 
382  offset = global;
383  offset.subtract( * this->giveNode(1)->giveCoordinates() );
384  answer.beProductOf(GtoLRotationMatrix, offset);
385 }
386 
389 {
390  IRResultType result; // Required by IR_GIVE_FIELD macro
391 
395 
396  directorType = 0; // default
398 
399  return this->NLStructuralElement :: initializeFrom(ir);
400 }
401 
402 
403 void
405 {
406  answer = {
407  D_u, D_v, D_w, R_u, R_v, R_w
408  };
409 }
410 
411 double
413 // Returns the portion of the receiver which is attached to gp.
414 {
415  double detJ, weight;
416  FloatMatrix jacobianMatrix(3, 3);
417 
418  FloatArray lcoords(3);
419  lcoords.at(1) = gp->giveNaturalCoordinate(1);
420  lcoords.at(2) = gp->giveNaturalCoordinate(2);
421  lcoords.at(3) = gp->giveNaturalCoordinate(3);
422 
423  weight = gp->giveWeight();
424 
425  this->giveJacobian(lcoords, jacobianMatrix);
426 
427  detJ = jacobianMatrix.giveDeterminant();
428  return detJ * weight;
429 }
430 
431 
432 void
434 // Returns the jacobianMatrix
435 {
436  FloatArray h(4);
437  FloatMatrix dn(4, 2);
438  FloatMatrix dndx(4, 2);
439 
440  jacobianMatrix.resize(3, 3);
441  // get node coordinates
442  double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
443  this->giveNodeCoordinates(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4);
444 
445  double r3 = lcoords.at(3);
446 
447  // get local director vector
448  FloatArray V1, V2, V3, V4;
449  this->giveLocalDirectorVectors(V1, V2, V3, V4);
450 
451  // get thickness
452  double a1, a2, a3, a4;
453  this->giveThickness(a1, a2, a3, a4);
454 
455  interp_lin.evalN( h, lcoords, FEIElementGeometryWrapper(this) );
456  interp_lin.evaldNdxi(dn, lcoords, FEIElementGeometryWrapper(this) );
457 
458  FloatArray hk1(4);
459  // derivatives of interpolation functions
460  // dh(r1,r2)/dr1
461  hk1.at(1) = dn.at(1, 1);
462  hk1.at(2) = dn.at(2, 1);
463  hk1.at(3) = dn.at(3, 1);
464  hk1.at(4) = dn.at(4, 1);
465 
466  FloatArray hk2(4);
467  // dh(r1,r2)/dr2
468  hk2.at(1) = dn.at(1, 2);
469  hk2.at(2) = dn.at(2, 2);
470  hk2.at(3) = dn.at(3, 2);
471  hk2.at(4) = dn.at(4, 2);
472 
473  double h1, h2, h3, h4;
474  // interpolation functions - h(r1,r2)
475  h1 = h.at(1);
476  h2 = h.at(2);
477  h3 = h.at(3);
478  h4 = h.at(4);
479 
480  // Jacobian Matrix
481  jacobianMatrix.at(1, 1) = hk1.at(1) * x1 + hk1.at(2) * x2 + hk1.at(3) * x3 + hk1.at(4) * x4 + r3 / 2. * ( a1 * hk1.at(1) * V1.at(1) + a2 * hk1.at(2) * V2.at(1) + a3 * hk1.at(3) * V3.at(1) + a4 * hk1.at(4) * V4.at(1) );
482  jacobianMatrix.at(1, 2) = hk1.at(1) * y1 + hk1.at(2) * y2 + hk1.at(3) * y3 + hk1.at(4) * y4 + r3 / 2. * ( a1 * hk1.at(1) * V1.at(2) + a2 * hk1.at(2) * V2.at(2) + a3 * hk1.at(3) * V3.at(2) + a4 * hk1.at(4) * V4.at(2) );
483  jacobianMatrix.at(1, 3) = r3 / 2. * ( a1 * hk1.at(1) * V1.at(3) + a2 * hk1.at(2) * V2.at(3) + a3 * hk1.at(3) * V3.at(3) + a4 * hk1.at(4) * V4.at(3) );
484  jacobianMatrix.at(2, 1) = hk2.at(1) * x1 + hk2.at(2) * x2 + hk2.at(3) * x3 + hk2.at(4) * x4 + r3 / 2. * ( a1 * hk2.at(1) * V1.at(1) + a2 * hk2.at(2) * V2.at(1) + a3 * hk2.at(3) * V3.at(1) + a4 * hk2.at(4) * V4.at(1) );
485  jacobianMatrix.at(2, 2) = hk2.at(1) * y1 + hk2.at(2) * y2 + hk2.at(3) * y3 + hk2.at(4) * y4 + r3 / 2. * ( a1 * hk2.at(1) * V1.at(2) + a2 * hk2.at(2) * V2.at(2) + a3 * hk2.at(3) * V3.at(2) + a4 * hk2.at(4) * V4.at(2) );
486  jacobianMatrix.at(2, 3) = r3 / 2. * ( a1 * hk2.at(1) * V1.at(3) + a2 * hk2.at(2) * V2.at(3) + a3 * hk2.at(3) * V3.at(3) + a4 * hk2.at(4) * V4.at(3) );
487  jacobianMatrix.at(3, 1) = 1. / 2. * ( a1 * h1 * V1.at(1) + a2 * h2 * V2.at(1) + a3 * h3 * V3.at(1) + a4 * h4 * V4.at(1) );
488  jacobianMatrix.at(3, 2) = 1. / 2. * ( a1 * h1 * V1.at(2) + a2 * h2 * V2.at(2) + a3 * h3 * V3.at(2) + a4 * h4 * V4.at(2) );
489  jacobianMatrix.at(3, 3) = 1. / 2. * ( a1 * h1 * V1.at(3) + a2 * h2 * V2.at(3) + a3 * h3 * V3.at(3) + a4 * h4 * V4.at(3) );
490 }
491 
492 void
494 {
495  // This element adds an additional stiffness for the so called drilling dofs.
496  NLStructuralElement :: computeStiffnessMatrix(answer, rMode, tStep);
497 
498  bool drillType = this->giveStructuralCrossSection()->give( CS_DrillingType, this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0) );
499  if ( drillType == 1 ) {
500  double relDrillCoeff = this->giveStructuralCrossSection()->give( CS_RelDrillingStiffness, this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0) );
501  FloatArray n;
502  IntArray drillDofs = {
503  6, 12, 18, 24
504  };
505  int j = 1;
506  while ( answer.at(j, j) == 0 ) {
507  j++;
508  }
509  drillCoeff = answer.at(j, j);
510  // find the smallest non-zero number on the diagonal
511  for ( int i = j; i <= 24; i++ ) {
512  if ( drillCoeff > answer.at(i, i) && answer.at(i, i) != 0 ) {
513  drillCoeff = answer.at(i, i);
514  }
515  }
516 
517  if ( relDrillCoeff == 0.0 ) {
518  relDrillCoeff = 0.001; // default
519  }
520  drillCoeff *= relDrillCoeff;
521 
522  FloatMatrix drillStiffness;
523  drillStiffness.resize(4, 4);
524  drillStiffness.zero();
525  for ( int i = 1; i <= 4; i++ ) {
526  drillStiffness.at(i, i) = drillCoeff;
527  }
528 
529  /*
530  * for ( GaussPoint *gp : *integrationRulesArray [ 0 ] ) {
531  * double dV = this->computeVolumeAround(gp);
532  * // double drillCoeff = this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp);
533  * double coeff = drillCoeff;
534  * if ( this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp)> 0)
535  * drillCoeff *= this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp);
536  * // Drilling stiffness is here for improved numerical properties
537  * this->interp_lin.evalN( n, gp->giveNaturalCoordinates(), FEIVoidCellGeometry() );
538  * for ( int j = 0; j < 4; j++ ) {
539  * n(j) -= 0.25;
540  * }
541  * drillStiffness.plusDyadSymmUpper(n, coeff * dV);
542  * }*/
543 
544  // drillStiffness.symmetrized();
545  answer.assemble(drillStiffness, drillDofs);
546  }
547 }
548 
549 void
550 MITC4Shell :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
551 {
552  // This element adds an additional stiffness for the so called drilling dofs.
553  NLStructuralElement :: giveInternalForcesVector(answer, tStep, useUpdatedGpRecord);
554 
555  bool drillType = this->giveStructuralCrossSection()->give( CS_DrillingType, this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0) );
556 
557  if ( drillType == 1 ) {
558  FloatArray n, tmp;
559  FloatArray drillUnknowns, drillMoment;
560  IntArray drillDofs = {
561  6, 12, 18, 24
562  };
563  this->computeVectorOf(VM_Total, tStep, tmp);
564  drillUnknowns.beSubArrayOf(tmp, drillDofs);
565 
566  /*
567  * for ( GaussPoint *gp : *integrationRulesArray [ 0 ] ) {
568  * double dV = this->computeVolumeAround(gp);
569  * // double drillCoeff = this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp);
570  * double coeff = drillCoeff;
571  * if ( this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp)> 0)
572  * drillCoeff *= this->giveStructuralCrossSection()->give(CS_DrillingStiffness, gp);
573  *
574  * this->interp_lin.evalN( n, gp->giveNaturalCoordinates(), FEIVoidCellGeometry() );
575  * for ( int j = 0; j < 4; j++ ) {
576  * n(j) -= 0.25;
577  * }
578  * double dtheta = n.dotProduct(drillUnknowns);
579  * drillMoment.add(coeff * dV * dtheta, n);
580  * }
581  */
582 
583  FloatMatrix drillStiffness;
584  drillStiffness.resize(4, 4);
585  drillStiffness.zero();
586  for ( int i = 1; i <= 4; i++ ) {
587  drillStiffness.at(i, i) = drillCoeff;
588  }
589 
590  drillMoment.beProductOf(drillStiffness, drillUnknowns);
591 
592  answer.assemble(drillMoment, drillDofs);
593  }
594 }
595 
596 void
598 // Returns the [6x20] strain-displacement matrix {B} of the receiver,
599 // evaluated at gp.
600 {
601  FloatArray h(4);
602  FloatMatrix jacobianMatrix(3, 3);
603  FloatMatrix dn(4, 2);
604  FloatMatrix dndx(4, 2);
605 
606  answer.resize(6, 24);
607  // get node coordinates
608  double x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4;
609  this->giveNodeCoordinates(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4);
610 
611  // get gp coordinates
612  double r1 = gp->giveNaturalCoordinate(1);
613  double r2 = gp->giveNaturalCoordinate(2);
614  double r3 = gp->giveNaturalCoordinate(3);
615 
616  FloatArray lcoords(3);
617  lcoords.at(1) = r1;
618  lcoords.at(2) = r2;
619  lcoords.at(3) = r3;
620 
621  // get director vector
622  FloatArray V1, V2, V3, V4;
623  this->giveLocalDirectorVectors(V1, V2, V3, V4);
624 
625  // get thickness
626  double a1, a2, a3, a4;
627  this->giveThickness(a1, a2, a3, a4);
628 
629  interp_lin.evalN( h, lcoords, FEIElementGeometryWrapper(this) );
630  interp_lin.evaldNdxi(dn, lcoords, FEIElementGeometryWrapper(this) );
631 
632  FloatArray hkx, hky;
633  this->givedNdx(hkx, hky, lcoords);
634 
635  // Jacobian Matrix
636  this->giveJacobian(lcoords, jacobianMatrix);
637  FloatMatrix inv(3, 3);
638  FloatMatrix inv2(2, 2);
639  inv.beInverseOf(jacobianMatrix);
640 
641  double sb = 2 * inv.at(1, 1) * inv.at(3, 3);
642  double sa = 2 * inv.at(1, 2) * inv.at(3, 3);
643  double cb = 2 * inv.at(2, 1) * inv.at(3, 3);
644  double ca = 2 * inv.at(2, 2) * inv.at(3, 3);
645 
646  FloatArray e2 = {
647  0, 1, 0
648  };
649 
650  FloatArray V11(3), V12(3), V13(3), V14(3), V21(3), V22(3), V23(3), V24(3);
651  V11.beVectorProductOf(e2, V1);
652  V11.normalize();
653  V12.beVectorProductOf(e2, V2);
654  V12.normalize();
655  V13.beVectorProductOf(e2, V3);
656  V13.normalize();
657  V14.beVectorProductOf(e2, V4);
658  V14.normalize();
659 
660  V21.beVectorProductOf(V1, V11);
661  V22.beVectorProductOf(V2, V12);
662  V23.beVectorProductOf(V3, V13);
663  V24.beVectorProductOf(V4, V14);
664 
665 
666  answer.zero();
667  answer.at(4, 1) = 1. / 32. * ( ( a1 * V1.at(1) + a2 * V2.at(1) ) * ( cb * ( 1. + r2 ) ) + ( a1 * V1.at(1) + a4 * V4.at(1) ) * ( ca * ( 1. + r1 ) ) );
668  answer.at(4, 2) = 1. / 32. * ( ( a1 * V1.at(2) + a2 * V2.at(2) ) * ( cb * ( 1. + r2 ) ) + ( a1 * V1.at(2) + a4 * V4.at(2) ) * ( ca * ( 1. + r1 ) ) );
669  answer.at(4, 3) = 1. / 32. * ( ( a1 * V1.at(3) + a2 * V2.at(3) ) * ( cb * ( 1. + r2 ) ) + ( a1 * V1.at(3) + a4 * V4.at(3) ) * ( ca * ( 1. + r1 ) ) );
670  answer.at(4, 4) = -a1 / 32. * ( ( V21.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( cb * ( 1. + r2 ) ) ) + ( V21.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( ca * ( 1. + r1 ) ) ) );
671  answer.at(4, 5) = a1 / 32. * ( ( V11.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( cb * ( 1. + r2 ) ) ) + ( V11.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( ca * ( 1. + r1 ) ) ) );
672 
673  answer.at(4, 7) = 1. / 32. * ( -( a1 * V1.at(1) + a2 * V2.at(1) ) * ( cb * ( 1. + r2 ) ) + ( a2 * V2.at(1) + a3 * V3.at(1) ) * ( ca * ( 1. - r1 ) ) );
674  answer.at(4, 8) = 1. / 32. * ( -( a1 * V1.at(2) + a2 * V2.at(2) ) * ( cb * ( 1. + r2 ) ) + ( a2 * V2.at(2) + a3 * V3.at(2) ) * ( ca * ( 1. - r1 ) ) );
675  answer.at(4, 9) = 1. / 32. * ( -( a1 * V1.at(3) + a2 * V2.at(3) ) * ( cb * ( 1. + r2 ) ) + ( a2 * V2.at(3) + a3 * V3.at(3) ) * ( ca * ( 1. - r1 ) ) );
676  answer.at(4, 10) = -a1 / 32. * ( ( V21.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( cb * ( 1. + r2 ) ) ) + ( V21.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( ca * ( 1. - r1 ) ) ) );
677  answer.at(4, 11) = a1 / 32. * ( ( V11.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( cb * ( 1. + r2 ) ) ) + ( V11.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( ca * ( 1. - r1 ) ) ) );
678 
679  answer.at(4, 13) = 1. / 32. * ( -( a3 * V3.at(1) + a4 * V4.at(1) ) * ( cb * ( 1. - r2 ) ) - ( a2 * V2.at(1) + a3 * V3.at(1) ) * ( ca * ( 1. - r1 ) ) );
680  answer.at(4, 14) = 1. / 32. * ( -( a3 * V3.at(2) + a4 * V4.at(2) ) * ( cb * ( 1. - r2 ) ) - ( a2 * V2.at(2) + a3 * V3.at(2) ) * ( ca * ( 1. - r1 ) ) );
681  answer.at(4, 15) = 1. / 32. * ( -( a3 * V3.at(3) + a4 * V4.at(3) ) * ( cb * ( 1. - r2 ) ) - ( a2 * V2.at(3) + a3 * V3.at(3) ) * ( ca * ( 1. - r1 ) ) );
682  answer.at(4, 16) = -a1 / 32. * ( ( V21.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( cb * ( 1. - r2 ) ) ) + ( V21.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( ca * ( 1. - r1 ) ) ) );
683  answer.at(4, 17) = a1 / 32. * ( ( V11.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( cb * ( 1. - r2 ) ) ) + ( V11.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( ca * ( 1. - r1 ) ) ) );
684 
685  answer.at(4, 19) = 1. / 32. * ( ( a3 * V3.at(1) + a4 * V4.at(1) ) * ( cb * ( 1. - r2 ) ) - ( a1 * V1.at(1) + a4 * V4.at(1) ) * ( ca * ( 1. + r1 ) ) );
686  answer.at(4, 20) = 1. / 32. * ( ( a3 * V3.at(2) + a4 * V4.at(2) ) * ( cb * ( 1. - r2 ) ) - ( a1 * V1.at(2) + a4 * V4.at(2) ) * ( ca * ( 1. + r1 ) ) );
687  answer.at(4, 21) = 1. / 32. * ( ( a3 * V3.at(3) + a4 * V4.at(3) ) * ( cb * ( 1. - r2 ) ) - ( a1 * V1.at(3) + a4 * V4.at(3) ) * ( ca * ( 1. + r1 ) ) );
688  answer.at(4, 22) = -a1 / 32. * ( ( V21.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( cb * ( 1. - r2 ) ) ) + ( V21.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( ca * ( 1. + r1 ) ) ) );
689  answer.at(4, 23) = a1 / 32. * ( ( V11.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( cb * ( 1. - r2 ) ) ) + ( V11.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( ca * ( 1. + r1 ) ) ) );
690 
691  answer.at(5, 1) = 1. / 32. * ( ( a1 * V1.at(1) + a2 * V2.at(1) ) * ( sb * ( 1. + r2 ) ) + ( a1 * V1.at(1) + a4 * V4.at(1) ) * ( sa * ( 1. + r1 ) ) );
692  answer.at(5, 2) = 1. / 32. * ( ( a1 * V1.at(2) + a2 * V2.at(2) ) * ( sb * ( 1. + r2 ) ) + ( a1 * V1.at(2) + a4 * V4.at(2) ) * ( sa * ( 1. + r1 ) ) );
693  answer.at(5, 3) = 1. / 32. * ( ( a1 * V1.at(3) + a2 * V2.at(3) ) * ( sb * ( 1. + r2 ) ) + ( a1 * V1.at(3) + a4 * V4.at(3) ) * ( sa * ( 1. + r1 ) ) );
694  answer.at(5, 4) = -a1 / 32. * ( ( V21.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( sb * ( 1. + r2 ) ) ) + ( V21.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( sa * ( 1. + r1 ) ) ) );
695  answer.at(5, 5) = a1 / 32. * ( ( V11.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( sb * ( 1. + r2 ) ) ) + ( V11.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( sa * ( 1. + r1 ) ) ) );
696 
697  answer.at(5, 7) = 1. / 32. * ( -( a1 * V1.at(1) + a2 * V2.at(1) ) * ( sb * ( 1. + r2 ) ) + ( a2 * V2.at(1) + a3 * V3.at(1) ) * ( sa * ( 1. - r1 ) ) );
698  answer.at(5, 8) = 1. / 32. * ( -( a1 * V1.at(2) + a2 * V2.at(2) ) * ( sb * ( 1. + r2 ) ) + ( a2 * V2.at(2) + a3 * V3.at(2) ) * ( sa * ( 1. - r1 ) ) );
699  answer.at(5, 9) = 1. / 32. * ( -( a1 * V1.at(3) + a2 * V2.at(3) ) * ( sb * ( 1. + r2 ) ) + ( a2 * V2.at(3) + a3 * V3.at(3) ) * ( sa * ( 1. - r1 ) ) );
700  answer.at(5, 10) = -a1 / 32. * ( ( V21.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( sb * ( 1. + r2 ) ) ) + ( V21.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( sa * ( 1. - r1 ) ) ) );
701  answer.at(5, 11) = a1 / 32. * ( ( V11.dotProduct({ x1 - x2, y1 - y2, z1 - z2 }) * ( sb * ( 1. + r2 ) ) ) + ( V11.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( sa * ( 1. - r1 ) ) ) );
702 
703  answer.at(5, 13) = 1. / 32. * ( -( a3 * V3.at(1) + a4 * V4.at(1) ) * ( sb * ( 1. - r2 ) ) - ( a2 * V2.at(1) + a3 * V3.at(1) ) * ( sa * ( 1. - r1 ) ) );
704  answer.at(5, 14) = 1. / 32. * ( -( a3 * V3.at(2) + a4 * V4.at(2) ) * ( sb * ( 1. - r2 ) ) - ( a2 * V2.at(2) + a3 * V3.at(2) ) * ( sa * ( 1. - r1 ) ) );
705  answer.at(5, 15) = 1. / 32. * ( -( a3 * V3.at(3) + a4 * V4.at(3) ) * ( sb * ( 1. - r2 ) ) - ( a2 * V2.at(3) + a3 * V3.at(3) ) * ( sa * ( 1. - r1 ) ) );
706  answer.at(5, 16) = -a1 / 32. * ( ( V21.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( sb * ( 1. - r2 ) ) ) + ( V21.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( sa * ( 1. - r1 ) ) ) );
707  answer.at(5, 17) = a1 / 32. * ( ( V11.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( sb * ( 1. - r2 ) ) ) + ( V11.dotProduct({ x2 - x3, y2 - y3, z2 - z3 }) * ( sa * ( 1. - r1 ) ) ) );
708 
709  answer.at(5, 19) = 1. / 32. * ( ( a3 * V3.at(1) + a4 * V4.at(1) ) * ( sb * ( 1. - r2 ) ) - ( a1 * V1.at(1) + a4 * V4.at(1) ) * ( sa * ( 1. + r1 ) ) );
710  answer.at(5, 20) = 1. / 32. * ( ( a3 * V3.at(2) + a4 * V4.at(2) ) * ( sb * ( 1. - r2 ) ) - ( a1 * V1.at(2) + a4 * V4.at(2) ) * ( sa * ( 1. + r1 ) ) );
711  answer.at(5, 21) = 1. / 32. * ( ( a3 * V3.at(3) + a4 * V4.at(3) ) * ( sb * ( 1. - r2 ) ) - ( a1 * V1.at(3) + a4 * V4.at(3) ) * ( sa * ( 1. + r1 ) ) );
712  answer.at(5, 22) = -a1 / 32. * ( ( V21.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( sb * ( 1. - r2 ) ) ) + ( V21.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( sa * ( 1. + r1 ) ) ) );
713  answer.at(5, 23) = a1 / 32. * ( ( V11.dotProduct({ x4 - x3, y4 - y3, z4 - z3 }) * ( sb * ( 1. - r2 ) ) ) + ( V11.dotProduct({ x1 - x4, y1 - y4, z1 - z4 }) * ( sa * ( 1. + r1 ) ) ) );
714 
715 
716 
717  answer.at(1, 1) = hkx.at(1);
718  answer.at(1, 4) = -r3 / 2. *a1 *hkx.at(1) * V21.at(1);
719  answer.at(1, 5) = r3 / 2. *a1 *hkx.at(1) * V11.at(1);
720  answer.at(1, 7) = hkx.at(2);
721  answer.at(1, 10) = -r3 / 2. *a2 *hkx.at(2) * V22.at(1);
722  answer.at(1, 11) = r3 / 2. *a2 *hkx.at(2) * V12.at(1);
723  answer.at(1, 13) = hkx.at(3);
724  answer.at(1, 16) = -r3 / 2. *a3 *hkx.at(3) * V23.at(1);
725  answer.at(1, 17) = r3 / 2. *a3 *hkx.at(3) * V13.at(1);
726  answer.at(1, 19) = hkx.at(4);
727  answer.at(1, 22) = -r3 / 2. *a4 *hkx.at(4) * V24.at(1);
728  answer.at(1, 23) = r3 / 2. *a4 *hkx.at(4) * V14.at(1);
729 
730  answer.at(2, 2) = hky.at(1);
731  answer.at(2, 4) = -r3 / 2. *a1 *hky.at(1) * V21.at(2);
732  answer.at(2, 5) = r3 / 2. *a1 *hky.at(1) * V11.at(2);
733 
734  answer.at(2, 8) = hky.at(2);
735  answer.at(2, 10) = -r3 / 2. *a2 *hky.at(2) * V22.at(2);
736  answer.at(2, 11) = r3 / 2. *a2 *hky.at(2) * V12.at(2);
737 
738  answer.at(2, 14) = hky.at(3);
739  answer.at(2, 16) = -r3 / 2. *a3 *hky.at(3) * V23.at(2);
740  answer.at(2, 17) = r3 / 2. *a3 *hky.at(3) * V13.at(2);
741 
742  answer.at(2, 20) = hky.at(4);
743  answer.at(2, 22) = -r3 / 2. *a4 *hky.at(4) * V24.at(2);
744  answer.at(2, 23) = r3 / 2. *a4 *hky.at(4) * V14.at(2);
745 
746  answer.at(6, 1) = hky.at(1);
747  answer.at(6, 2) = hkx.at(1);
748  answer.at(6, 4) = -r3 / 2. * a1 * ( hkx.at(1) * V21.at(2) + hky.at(1) * V21.at(1) );
749  answer.at(6, 5) = r3 / 2. * a1 * ( hky.at(1) * V11.at(1) + hky.at(1) * V11.at(2) );
750 
751  answer.at(6, 7) = hky.at(2);
752  answer.at(6, 8) = hkx.at(2);
753  answer.at(6, 10) = -r3 / 2. * a2 * ( hkx.at(2) * V22.at(2) + hky.at(2) * V22.at(1) );
754  answer.at(6, 11) = r3 / 2. * a2 * ( hky.at(2) * V12.at(1) + hky.at(2) * V12.at(2) );
755 
756  answer.at(6, 13) = hky.at(3);
757  answer.at(6, 14) = hkx.at(3);
758  answer.at(6, 16) = -r3 / 2. * a3 * ( hkx.at(3) * V23.at(2) + hky.at(3) * V23.at(1) );
759  answer.at(6, 17) = r3 / 2. * a3 * ( hky.at(3) * V13.at(1) + hky.at(3) * V13.at(2) );
760 
761  answer.at(6, 19) = hky.at(4);
762  answer.at(6, 20) = hkx.at(4);
763  answer.at(6, 22) = -r3 / 2. * a4 * ( hkx.at(4) * V24.at(2) + hky.at(4) * V24.at(1) );
764  answer.at(6, 23) = r3 / 2. * a4 * ( hky.at(4) * V14.at(1) + hky.at(4) * V14.at(2) );
765 }
766 
767 void
768 MITC4Shell :: giveThickness(double &a1, double &a2, double &a3, double &a4)
769 {
770  FloatArray *c1, *c2, *c3, *c4;
771  c1 = this->giveNode(1)->giveCoordinates();
772  c2 = this->giveNode(2)->giveCoordinates();
773  c3 = this->giveNode(3)->giveCoordinates();
774  c4 = this->giveNode(4)->giveCoordinates();
775 
776  a1 = this->giveCrossSection()->give(CS_Thickness, * c1, this, false);
777  a2 = this->giveCrossSection()->give(CS_Thickness, * c2, this, false);
778  a3 = this->giveCrossSection()->give(CS_Thickness, * c3, this, false);
779  a4 = this->giveCrossSection()->give(CS_Thickness, * c4, this, false);
780 }
781 
782 
783 const FloatMatrix *
785 // Returns the rotation matrix of the receiver of the size [3,3]
786 // coords(local) = T * coords(global)
787 //
788 // local coordinate (described by vector triplet e1',e2',e3') is defined as follows:
789 //
790 // e1' : [N4-N3] Ni - means i - th node
791 // help : [N2-N3]
792 // e3' : e1' x help
793 // e2' : e3' x e1'
794 {
795  if ( !GtoLRotationMatrix.isNotEmpty() ) {
796  FloatArray e1, e2, e3;
797 
798 
799  this->computeLocalBaseVectors(e1, e2, e3);
800 
802 
803  for ( int i = 1; i <= 3; i++ ) {
804  GtoLRotationMatrix.at(1, i) = e1.at(i);
805  GtoLRotationMatrix.at(2, i) = e2.at(i);
806  GtoLRotationMatrix.at(3, i) = e3.at(i);
807  }
808  }
809 
810  return & GtoLRotationMatrix;
811 }
812 
813 void
815 {
816  FloatArray help;
817  FloatArray coordA, coordB;
818 
819  // compute A - (node2+node3)/2
820  coordA.beDifferenceOf( * this->giveNode(2)->giveCoordinates(), * this->giveNode(3)->giveCoordinates() );
821  coordA.times(0.5);
822  coordA.add( * this->giveNode(3)->giveCoordinates() );
823 
824  // compute B - (node1+node4)/2
825  coordB.beDifferenceOf( * this->giveNode(1)->giveCoordinates(), * this->giveNode(4)->giveCoordinates() );
826  coordB.times(0.5);
827  coordB.add( * this->giveNode(4)->giveCoordinates() );
828 
829  // compute e1' = [B-A]
830  e1.beDifferenceOf(coordB, coordA);
831 
832 
833  // compute A - (node3+node4)/2
834  coordA.beDifferenceOf( * this->giveNode(4)->giveCoordinates(), * this->giveNode(3)->giveCoordinates() );
835  coordA.times(0.5);
836  coordA.add( * this->giveNode(3)->giveCoordinates() );
837 
838  // compute B - (node2+node1)/2
839  coordB.beDifferenceOf( * this->giveNode(1)->giveCoordinates(), * this->giveNode(2)->giveCoordinates() );
840  coordB.times(0.5);
841  coordB.add( * this->giveNode(2)->giveCoordinates() );
842 
843  // compute help = [B-A]
844  help.beDifferenceOf(coordB, coordA);
845 
846  // let us normalize e1'
847  e1.normalize();
848 
849  // compute e3' : vector product of e1' x help
850  e3.beVectorProductOf(e1, help);
851  e3.normalize();
852 
853  // now from e3' x e1' compute e2'
854  e2.beVectorProductOf(e3, e1);
855 }
856 
857 void
859 // Returns the rotation matrix of the reciever of the size [3,3]
860 // {alpha_i,beta_i} = Ti * {rotL_xi, rotL_yi, rotL_zi}
861 // alpha_i, beta_i - rotations about the components of director vector at node i
862 // r1_i, r2_i, r3_i, - rotations about local coordinates e1', e2', e3'
863 //
864 // local coordinate (described by vector triplet e1',e2',e3') is defined as follows:
865 //
866 // e1' : [N4-N3] Ni - means i - th node
867 // help : [N2-N3]
868 // e3' : e1' x help
869 // e2' : e3' x e1'
870 {
871  FloatArray e1, e2, e3;
872 
873  this->computeLocalBaseVectors(e1, e2, e3);
874 
875  // get director vector
876  FloatArray V1, V2, V3, V4;
877  this->giveDirectorVectors(V1, V2, V3, V4);
878 
879  FloatArray V11(3), V12(3), V13(3), V14(3), V21(3), V22(3), V23(3), V24(3);
880  V11.beVectorProductOf(e2, V1);
881  V11.normalize();
882  V12.beVectorProductOf(e2, V2);
883  V12.normalize();
884  V13.beVectorProductOf(e2, V3);
885  V13.normalize();
886  V14.beVectorProductOf(e2, V4);
887  V14.normalize();
888 
889  V21.beVectorProductOf(V1, V11);
890  V22.beVectorProductOf(V2, V12);
891  V23.beVectorProductOf(V3, V13);
892  V24.beVectorProductOf(V4, V14);
893 
894  answer1.resize(3, 3);
895  answer2.resize(3, 3);
896  answer3.resize(3, 3);
897  answer4.resize(3, 3);
898 
899  answer1.at(1, 1) = V11.dotProduct(e1);
900  answer1.at(1, 2) = V11.dotProduct(e2);
901  answer1.at(1, 3) = V11.dotProduct(e3);
902  answer1.at(2, 1) = V21.dotProduct(e1);
903  answer1.at(2, 2) = V21.dotProduct(e2);
904  answer1.at(2, 3) = V21.dotProduct(e3);
905  answer1.at(3, 1) = V1.dotProduct(e1);
906  answer1.at(3, 2) = V1.dotProduct(e2);
907  answer1.at(3, 3) = V1.dotProduct(e3);
908 
909  answer2.at(1, 1) = V12.dotProduct(e1);
910  answer2.at(1, 2) = V12.dotProduct(e2);
911  answer2.at(1, 3) = V12.dotProduct(e3);
912  answer2.at(2, 1) = V22.dotProduct(e1);
913  answer2.at(2, 2) = V22.dotProduct(e2);
914  answer2.at(2, 3) = V22.dotProduct(e3);
915  answer2.at(3, 1) = V2.dotProduct(e1);
916  answer2.at(3, 2) = V2.dotProduct(e2);
917  answer2.at(3, 3) = V2.dotProduct(e3);
918 
919  answer3.at(1, 1) = V13.dotProduct(e1);
920  answer3.at(1, 2) = V13.dotProduct(e2);
921  answer3.at(1, 3) = V13.dotProduct(e3);
922  answer3.at(2, 1) = V23.dotProduct(e1);
923  answer3.at(2, 2) = V23.dotProduct(e2);
924  answer3.at(2, 3) = V23.dotProduct(e3);
925  answer3.at(3, 1) = V3.dotProduct(e1);
926  answer3.at(3, 2) = V3.dotProduct(e2);
927  answer3.at(3, 3) = V3.dotProduct(e3);
928 
929  answer4.at(1, 1) = V14.dotProduct(e1);
930  answer4.at(1, 2) = V14.dotProduct(e2);
931  answer4.at(1, 3) = V14.dotProduct(e3);
932  answer4.at(2, 1) = V24.dotProduct(e1);
933  answer4.at(2, 2) = V24.dotProduct(e2);
934  answer4.at(2, 3) = V24.dotProduct(e3);
935  answer4.at(3, 1) = V4.dotProduct(e1);
936  answer4.at(3, 2) = V4.dotProduct(e2);
937  answer4.at(3, 3) = V4.dotProduct(e3);
938 }
939 
940 bool
942 // Returns the rotation matrix of the receiver of the size [24,24]
943 // r(local) = T * r(global)
944 // for one node (r written transposed): {u,v,w,alpha,beta} = T * {u,v,w,r1,r2,r3}
945 
946 {
948 
949  answer.resize(24, 24);
950  answer.zero();
951 
952  FloatMatrix LtoDir1, LtoDir2, LtoDir3, LtoDir4;
953  this->computeLToDirectorRotationMatrix(LtoDir1, LtoDir2, LtoDir3, LtoDir4);
954 
955  for ( int i = 0; i <= 3; i++ ) {
956  answer.setSubMatrix(GtoLRotationMatrix, i * 6 + 1, i * 6 + 1);
957  }
959  FloatMatrix help;
960 
961  help.beProductOf(LtoDir1, GtoLRotationMatrix);
962  answer.setSubMatrix(help, 4, 4);
963 
964  help.beProductOf(LtoDir2, GtoLRotationMatrix);
965  answer.setSubMatrix(help, 10, 10);
966 
967  help.beProductOf(LtoDir3, GtoLRotationMatrix);
968  answer.setSubMatrix(help, 16, 16);
969 
970  help.beProductOf(LtoDir4, GtoLRotationMatrix);
971  answer.setSubMatrix(help, 22, 22);
972 
973  return 1;
974 }
975 
976 
977 void
979 {
980  this->giveStructuralCrossSection()->giveRealStress_3dDegeneratedShell(answer, gp, strain, tStep);
981 }
982 
983 
984 void
986 // returns characteristic tensor of the receiver at given gp and tStep
987 {
988  answer.resize(3, 3);
989  answer.zero();
991  StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveStructuralCrossSection()->giveMaterial(gp) );
992 
993  if ( type == GlobalForceTensor ) {
994  FloatArray stress, localStress, localStrain;
995  this->computeStrainVector(localStrain, gp, tStep);
996  this->computeStressVector(localStress, localStrain, gp, tStep);
997  mat->transformStressVectorTo(stress, GtoLRotationMatrix, localStress, false);
998 
999  answer.at(1, 1) = stress.at(1);
1000  answer.at(2, 2) = stress.at(2);
1001  answer.at(3, 3) = stress.at(3);
1002  answer.at(1, 2) = stress.at(4);
1003  answer.at(2, 1) = stress.at(4);
1004  answer.at(2, 3) = stress.at(5);
1005  answer.at(3, 2) = stress.at(5);
1006  answer.at(1, 3) = stress.at(6);
1007  answer.at(3, 1) = stress.at(6);
1008  } else if ( type == GlobalStrainTensor ) {
1009  FloatArray strain, localStrain;
1010  this->computeStrainVector(localStrain, gp, tStep);
1011  mat->transformStrainVectorTo(strain, GtoLRotationMatrix, localStrain, false);
1012 
1013  answer.at(1, 1) = strain.at(1);
1014  answer.at(2, 2) = strain.at(2);
1015  answer.at(3, 3) = strain.at(3);
1016  answer.at(2, 3) = strain.at(4) / 2.;
1017  answer.at(3, 2) = strain.at(4) / 2.;
1018  answer.at(1, 3) = strain.at(5) / 2.;
1019  answer.at(3, 1) = strain.at(5) / 2.;
1020  answer.at(1, 2) = strain.at(6) / 2.;
1021  answer.at(2, 1) = strain.at(6) / 2.;
1022  } else {
1023  OOFEM_ERROR("unsupported tensor mode");
1024  exit(1);
1025  }
1026 }
1027 
1028 void
1030 // Performs end-of-step operations.
1031 {
1032  FloatArray v;
1033  GaussPoint *gp;
1034 
1035  fprintf(file, "element %d (%8d):\n", this->giveLabel(), number);
1036 
1037  for ( int i = 0; i < nPointsXY; i++ ) {
1038  fprintf(file, " GP %d :", i + 1);
1039 
1040  this->giveMidplaneIPValue(v, i, IST_ShellForceTensor, tStep);
1041  fprintf(file, " forces ");
1042  for ( auto &val : v ) {
1043  fprintf(file, " %.4e", val);
1044  }
1045 
1046  this->giveMidplaneIPValue(v, i, IST_ShellMomentTensor, tStep);
1047  fprintf(file, "\n moments ");
1048  for ( auto &val : v ) {
1049  fprintf(file, " %.4e", val);
1050  }
1051 
1052  this->giveMidplaneIPValue(v, i, IST_ShellStrainTensor, tStep);
1053  fprintf(file, "\n strains ");
1054  for ( auto &val : v ) {
1055  fprintf(file, " %.4e", val);
1056  }
1057 
1058  this->giveMidplaneIPValue(v, i, IST_CurvatureTensor, tStep);
1059  fprintf(file, "\n curvatures ");
1060  for ( auto &val : v ) {
1061  fprintf(file, " %.4e", val);
1062  }
1063 
1064  for ( int j = 0; j < nPointsZ; j++ ) {
1065  gp = integrationRulesArray [ 0 ]->getIntegrationPoint(nPointsZ * i + j);
1066 
1067  fprintf(file, "\n GP %d.%d :", i + 1, j + 1);
1068 
1069  this->giveIPValue(v, gp, IST_StrainTensor, tStep);
1070  fprintf(file, " strains ");
1071  for ( auto &val : v ) {
1072  fprintf(file, " %.4e", val);
1073  }
1074 
1075  this->giveIPValue(v, gp, IST_StressTensor, tStep);
1076  fprintf(file, "\n stresses ");
1077  for ( auto &val : v ) {
1078  fprintf(file, " %.4e", val);
1079  }
1080  }
1081  fprintf(file, "\n");
1082  }
1083 }
1084 
1085 void
1087 {
1088  GaussPoint *gp = NULL;
1089 
1090  if ( type == IST_ShellMomentTensor || type == IST_ShellForceTensor ) {
1091  double J, thickness, z, w;
1092  FloatArray mLocal;
1093  mLocal.resize(6);
1094  mLocal.zero();
1095 
1096  for ( int i = 0; i < nPointsZ; i++ ) {
1097  gp = integrationRulesArray [ 0 ]->getIntegrationPoint(nPointsZ * gpXY + i);
1098  thickness = this->giveCrossSection()->give(CS_Thickness, gp->giveGlobalCoordinates(), this, false);
1099  J = thickness / 2.0;
1100  if ( type == IST_ShellMomentTensor ) {
1101  z = gp->giveNaturalCoordinates().at(3) * ( thickness / 2 );
1102  } else if ( type == IST_ShellForceTensor ) {
1103  z = 1;
1104  }
1105  w = gp->giveWeight() * J * z;
1106 
1107  FloatArray localStress, localStrain;
1108  this->computeStrainVector(localStrain, gp, tStep);
1109  this->computeStressVector(localStress, localStrain, gp, tStep);
1110  mLocal.add(w, localStress);
1111  }
1112 
1113  // local to global
1114  StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveStructuralCrossSection()->giveMaterial(gp) );
1115  this->computeGtoLRotationMatrix();
1116  mat->transformStressVectorTo(answer, GtoLRotationMatrix, mLocal, false);
1117  } else if ( type == IST_CurvatureTensor ) {
1118  FloatArray h;
1119  FloatMatrix dn;
1120  FloatArray coords;
1121 
1122  gp = integrationRulesArray [ 0 ]->getIntegrationPoint(nPointsZ * gpXY);
1123  coords = gp->giveNaturalCoordinates();
1124 
1125  StructuralMaterial *mat = dynamic_cast< StructuralMaterial * >( this->giveStructuralCrossSection()->giveMaterial(gp) );
1126 
1127  FloatArray hkx, hky;
1128  this->givedNdx(hkx, hky, coords);
1129 
1130  FloatArray dofs(24);
1131 
1132 
1133  FloatArray rotX(4), rotY(4);
1134 
1135  this->computeVectorOf(VM_Total, tStep, dofs);
1136  for ( int i = 0; i < 4; i++ ) {
1137  rotX(i) = dofs.at(i * 6 + 4);
1138  rotY(i) = dofs.at(i * 6 + 5);
1139  }
1140  FloatArray cLocal(6);
1141  cLocal.zero();
1142  cLocal.at(1) = rotY.dotProduct(hkx);
1143  cLocal.at(2) = -rotX.dotProduct(hky);
1144  cLocal.at(6) = rotY.dotProduct(hky) - rotX.dotProduct(hkx);
1145 
1146  mat->transformStrainVectorTo(answer, GtoLRotationMatrix, cLocal, false);
1147  } else if ( type == IST_ShellStrainTensor ) {
1148  FloatArray coords;
1149  gp = integrationRulesArray [ 0 ]->getIntegrationPoint(nPointsZ * gpXY);
1150  coords = gp->giveNaturalCoordinates();
1151  coords.at(3) = 0; //set to midplane
1152  IntegrationRule *iRule = new GaussIntegrationRule(1, this, 1, 10);
1153  GaussPoint *midGP = new GaussPoint( iRule, 1, coords, 1, this->giveMaterialMode() );
1154 
1155  this->giveIPValue(answer, midGP, IST_StrainTensor, tStep);
1156  } else {
1157  OOFEM_ERROR("MITC4Shell :: giveMidplaneIPValue - unknown type");
1158  }
1159 
1160  return;
1161 }
1162 
1163 void
1165 {
1166  FloatArray h, hk1(4), hk2(4);
1167  FloatMatrix dn, dndx, jacobianMatrix, inv, inv2;
1168 
1169  interp_lin.evalN( h, coords, FEIElementGeometryWrapper(this) );
1170  interp_lin.evaldNdxi(dn, coords, FEIElementGeometryWrapper(this) );
1171 
1172  // derivatives of interpolation functions
1173  // dh(r1,r2)/dr1
1174  hk1.at(1) = dn.at(1, 1);
1175  hk1.at(2) = dn.at(2, 1);
1176  hk1.at(3) = dn.at(3, 1);
1177  hk1.at(4) = dn.at(4, 1);
1178 
1179  // dh(r1,r2)/dr2
1180  hk2.at(1) = dn.at(1, 2);
1181  hk2.at(2) = dn.at(2, 2);
1182  hk2.at(3) = dn.at(3, 2);
1183  hk2.at(4) = dn.at(4, 2);
1184 
1185  // Jacobian Matrix
1186  this->giveJacobian(coords, jacobianMatrix);
1187  inv.beInverseOf(jacobianMatrix);
1188 
1189  inv2.beSubMatrixOf(inv, 1, 2, 1, 2);
1190  dndx.beProductTOf(dn, inv2);
1191 
1192  hkx.resize(4);
1193  hkx.at(1) = dndx.at(1, 1);
1194  hkx.at(2) = dndx.at(2, 1);
1195  hkx.at(3) = dndx.at(3, 1);
1196  hkx.at(4) = dndx.at(4, 1);
1197 
1198  hky.resize(4);
1199  hky.at(1) = dndx.at(1, 2);
1200  hky.at(2) = dndx.at(2, 2);
1201  hky.at(3) = dndx.at(3, 2);
1202  hky.at(4) = dndx.at(4, 2);
1203 
1204  return;
1205 }
1206 
1207 void
1209 {
1211 }
1212 
1213 int
1215 {
1216  FloatMatrix globTensor;
1217  CharTensor cht;
1218 
1219  answer.resize(6);
1220 
1221  if ( type == IST_StrainTensor ) {
1222  cht = GlobalStrainTensor;
1223 
1224  this->giveCharacteristicTensor(globTensor, cht, gp, tStep);
1225 
1226  answer.at(1) = globTensor.at(1, 1); //xx
1227  answer.at(2) = globTensor.at(2, 2); //yy
1228  answer.at(3) = globTensor.at(3, 3); //zz
1229  answer.at(4) = 2 * globTensor.at(2, 3); //yz
1230  answer.at(5) = 2 * globTensor.at(1, 3); //xz
1231  answer.at(6) = 2 * globTensor.at(1, 2); //xy
1232 
1233  return 1;
1234  } else if ( type == IST_StressTensor ) {
1235  cht = GlobalForceTensor;
1236 
1237  this->giveCharacteristicTensor(globTensor, cht, gp, tStep);
1238 
1239  answer.at(1) = globTensor.at(1, 1); //xx
1240  answer.at(2) = globTensor.at(2, 2); //yy
1241  answer.at(3) = globTensor.at(3, 3); //zz
1242  answer.at(4) = globTensor.at(2, 3); //yz
1243  answer.at(5) = globTensor.at(1, 3); //xz
1244  answer.at(6) = globTensor.at(1, 2); //xy
1245 
1246  return 1;
1247  } else if ( type == IST_ShellMomentTensor || type == IST_ShellForceTensor || type == IST_CurvatureTensor || type == IST_ShellStrainTensor ) {
1248  int gpnXY = ( gp->giveNumber() - 1 ) / 2;
1249  this->giveMidplaneIPValue(answer, gpnXY, type, tStep);
1250 
1251  return 1;
1252  } else {
1253  return NLStructuralElement :: giveIPValue(answer, gp, type, tStep);
1254  }
1255 }
1256 
1257 
1258 bool
1260 //converts global coordinates to local planar area coordinates,
1261 //does not return a coordinate in the thickness direction, but
1262 //does check that the point is in the element thickness
1263 {
1264  // rotate the input point Coordinate System into the element CS
1265  FloatArray inputCoords_ElCS;
1266  std :: vector< FloatArray >lc(3);
1267  FloatArray llc;
1268  this->giveLocalCoordinates( inputCoords_ElCS, const_cast< FloatArray & >(coords) );
1269  for ( int _i = 0; _i < 4; _i++ ) {
1270  this->giveLocalCoordinates( lc [ _i ], * this->giveNode(_i + 1)->giveCoordinates() );
1271  }
1272  bool inplane = interp_lin.global2local( llc, inputCoords_ElCS, FEIVertexListGeometryWrapper(lc) ) > 0;
1273  answer.resize(2);
1274  answer.at(1) = inputCoords_ElCS.at(1);
1275  answer.at(2) = inputCoords_ElCS.at(2);
1276  GaussPoint _gp(NULL, 1, answer, 2.0, _2dPlate);
1277  // now check if the third local coordinate is within the thickness of element
1278  bool outofplane = ( fabs( inputCoords_ElCS.at(3) ) <= this->giveCrossSection()->give(CS_Thickness, & _gp) / 2. );
1279 
1280  return inplane && outofplane;
1281 }
1282 
1283 
1284 
1285 int
1287 {
1288  //double l1 = lcoords.at(1);
1289  //double l2 = lcoords.at(2);
1290  //double l3 = lcoords.at(3);
1291  FloatMatrix N;
1292 
1293  computeNmatrixAt(lcoords, N);
1294 
1295  answer.resize(3);
1296  for ( int _i = 1; _i <= 3; _i++ ) {
1297  answer.at(_i) = N.at(_i,_i) * this->giveNode(1)->giveCoordinate(_i) + N.at(_i,_i+6) *this->giveNode(2)->giveCoordinate(_i) + N.at(_i,_i+12)*this->giveNode(3)->giveCoordinate(_i) + N.at(_i,_i+18)*this->giveNode(4)->giveCoordinate(_i);
1298  }
1299  return true;
1300 }
1301 
1302 
1303 int
1305 // Returns the rotation matrix of the receiver of the size [5,6]
1306 // f(local) = T * f(global)
1307 {
1308  this->computeGtoLRotationMatrix();
1309 
1310  answer.resize(6, 6);
1311  answer.zero();
1312 
1313  for ( int i = 1; i <= 3; i++ ) {
1314  answer.at(1, i) = answer.at(4, i + 3) = GtoLRotationMatrix.at(1, i);
1315  answer.at(2, i) = answer.at(5, i + 3) = GtoLRotationMatrix.at(2, i);
1316  answer.at(3, i) = answer.at(6, i + 3) = GtoLRotationMatrix.at(3, i);
1317  }
1318 
1319  return 1;
1320 }
1321 
1322 
1323 void
1325  InternalStateType type, TimeStep *tStep)
1326 {
1327  double x1 = 0.0, x2 = 0.0, y = 0.0;
1328  FloatMatrix A(3, 3);
1329  FloatMatrix b, r;
1330  FloatArray val;
1331  double u, v;
1332 
1333 
1334  int size = 0;
1335 
1336  for ( GaussPoint *gp : *integrationRulesArray [ 0 ] ) {
1337  giveIPValue(val, gp, type, tStep);
1338  if ( size == 0 ) {
1339  size = val.giveSize();
1340  b.resize(3, size);
1341  r.resize(3, size);
1342  A.zero();
1343  r.zero();
1344  }
1345 
1346  const FloatArray &coord = gp->giveNaturalCoordinates();
1347  u = coord.at(1);
1348  v = coord.at(2);
1349 
1350  A.at(1, 1) += 1;
1351  A.at(1, 2) += u;
1352  A.at(1, 3) += v;
1353  A.at(2, 1) += u;
1354  A.at(2, 2) += u * u;
1355  A.at(2, 3) += u * v;
1356  A.at(3, 1) += v;
1357  A.at(3, 2) += v * u;
1358  A.at(3, 3) += v * v;
1359 
1360  for ( int j = 1; j <= size; j++ ) {
1361  y = val.at(j);
1362  r.at(1, j) += y;
1363  r.at(2, j) += y * u;
1364  r.at(3, j) += y * v;
1365  }
1366  }
1367 
1368  A.solveForRhs(r, b);
1369 
1370  switch ( node ) {
1371  case 1:
1372  x1 = 1.0;
1373  x2 = 1.0;
1374  break;
1375  case 2:
1376  x1 = -1.0;
1377  x2 = 1.0;
1378  break;
1379  case 3:
1380  x1 = -1.0;
1381  x2 = -1.0;
1382  break;
1383  case 4:
1384  x1 = 1.0;
1385  x2 = -1.0;
1386  break;
1387  default:
1388  OOFEM_ERROR("unsupported node");
1389  }
1390 
1391  answer.resize(size);
1392  for ( int j = 1; j <= size; j++ ) {
1393  answer.at(j) = b.at(1, j) + x1 *b.at(2, j) + x2 *b.at(3, j);
1394  }
1395 }
1396 
1397 
1398 void
1400 {
1401  if ( iEdge == 1 ) { // edge between nodes 1 2
1402  answer = {
1403  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12
1404  };
1405  } else if ( iEdge == 2 ) { // edge between nodes 2 3
1406  answer = {
1407  7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18
1408  };
1409  } else if ( iEdge == 3 ) { // edge between nodes 3 4
1410  answer = {
1411  13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24
1412  };
1413  } else if ( iEdge == 4 ) { // edge between nodes 4 1
1414  answer = {
1415  19, 20, 21, 22, 23, 24, 1, 2, 3, 4, 5, 6
1416  };
1417  } else {
1418  OOFEM_ERROR("wrong edge number");
1419  }
1420 }
1421 
1422 
1423 double
1425 {
1426  std :: vector< FloatArray >lc = {
1428  };
1429  this->giveNodeCoordinates( lc [ 0 ].at(1), lc [ 1 ].at(1), lc [ 2 ].at(1), lc [ 3 ].at(1),
1430  lc [ 0 ].at(2), lc [ 1 ].at(2), lc [ 2 ].at(2), lc [ 3 ].at(2),
1431  lc [ 0 ].at(3), lc [ 1 ].at(3), lc [ 2 ].at(3), lc [ 3 ].at(3) );
1432 
1433 
1435  return detJ * gp->giveWeight();
1436 }
1437 
1438 /*
1439  * void
1440  * MITC4Shell :: computeEdgeIpGlobalCoords(FloatArray &answer, GaussPoint *gp, int iEdge)
1441  * {
1442  * std :: vector< FloatArray > lc = {FloatArray(3), FloatArray(3), FloatArray(3), FloatArray(3)};
1443  * this->giveNodeCoordinates(lc[0].at(1), lc[1].at(1), lc[2].at(1), lc[3].at(1),
1444  * lc[0].at(2), lc[1].at(2), lc[2].at(2), lc[3].at(2),
1445  * lc[0].at(3), lc[1].at(3), lc[2].at(3), lc[3].at(3));
1446  *
1447  * FloatArray local;
1448  * this->interp_lin.edgeLocal2global( local, iEdge, gp->giveNaturalCoordinates(), FEIVertexListGeometryWrapper(lc) );
1449  * local.resize(3);
1450  * local.at(3) = 0.;
1451  *
1452  * answer = local; // MITC4 - todo
1453  * // answer.beProductOf(this->lcsMatrix, local);
1454  * }
1455  */
1456 
1457 int
1459 {
1460  FloatArray e1, e2, e3, xl, yl;
1461  this->computeLocalBaseVectors(e1, e2, e3);
1462 
1463  IntArray edgeNodes;
1464  this->interp_lin.computeLocalEdgeMapping(edgeNodes, iEdge);
1465 
1466  xl.beDifferenceOf( * this->giveNode( edgeNodes.at(2) )->giveCoordinates(), * this->giveNode( edgeNodes.at(1) )->giveCoordinates() );
1467 
1468  xl.normalize();
1469  yl.beVectorProductOf(e3, xl);
1470 
1471  answer.resize(6, 6);
1472  answer.zero();
1473 
1474  answer.at(1, 1) = answer.at(4, 4) = e1.dotProduct(xl);
1475  answer.at(1, 2) = answer.at(4, 5) = e1.dotProduct(yl);
1476  answer.at(1, 3) = answer.at(4, 6) = e1.dotProduct(e3);
1477  answer.at(2, 1) = answer.at(5, 4) = e2.dotProduct(xl);
1478  answer.at(2, 2) = answer.at(5, 5) = e2.dotProduct(yl);
1479  answer.at(2, 3) = answer.at(5, 6) = e2.dotProduct(e3);
1480  answer.at(3, 1) = answer.at(6, 4) = e3.dotProduct(xl);
1481  answer.at(3, 2) = answer.at(6, 5) = e3.dotProduct(yl);
1482  answer.at(3, 3) = answer.at(6, 6) = e3.dotProduct(e3);
1483 
1484 
1485  return 1;
1486 }
1487 
1488 
1489 
1490 void
1492 {
1493  const FloatArray &coords2 = sgp->giveNaturalCoordinates();
1494  FloatArray coords(3);
1495  coords.at(1) = coords2.at(1);
1496  coords.at(2) = coords2.at(2);
1497  coords.at(3) = 0.0;
1498  this->computeNmatrixAt(coords, answer);
1499 }
1500 
1501 void
1503 {
1504  int i;
1505  answer.resize(24);
1506  answer.zero();
1507  if ( iSurf == 1 ) {
1508  for ( i = 1; i <= 24; i++ ) {
1509  answer.at(i) = i;
1510  }
1511  } else {
1512  OOFEM_ERROR("wrong surface number");
1513  }
1514 }
1515 
1518 {
1519  IntegrationRule *iRule = new GaussIntegrationRule(1, this, 1, 1);
1520  int npoints = iRule->getRequiredNumberOfIntegrationPoints(_Square, approxOrder);
1521  iRule->SetUpPointsOnSquare(npoints, _Unknown);
1522  return iRule;
1523 }
1524 
1525 
1526 
1527 double
1529 {
1530  double detJ, weight;
1531  FloatMatrix jacobianMatrix(2, 2);
1532  FloatMatrix dn;
1533  FloatArray lcoords(2);
1534  lcoords.at(1) = gp->giveNaturalCoordinate(1);
1535  lcoords.at(2) = gp->giveNaturalCoordinate(2);
1536  FloatArray x(4), y(4), z(4);
1537  this->giveNodeCoordinates( x.at(1), x.at(2), x.at(3), x.at(4), y.at(1), y.at(2), y.at(3), y.at(4), z.at(1), z.at(2), z.at(3), z.at(4) );
1538 
1539 
1540  weight = gp->giveWeight();
1541 
1542  interp_lin.evaldNdxi(dn, lcoords, FEIElementGeometryWrapper(this) );
1543 
1544  for ( int i = 1; i <= dn.giveNumberOfRows(); i++ ) {
1545  double xl = x.at(i);
1546  double yl = y.at(i);
1547 
1548  jacobianMatrix.at(1, 1) += dn.at(i, 1) * xl;
1549  jacobianMatrix.at(1, 2) += dn.at(i, 1) * yl;
1550  jacobianMatrix.at(2, 1) += dn.at(i, 2) * xl;
1551  jacobianMatrix.at(2, 2) += dn.at(i, 2) * yl;
1552  }
1553 
1554  detJ = jacobianMatrix.giveDeterminant();
1555  return detJ * weight;
1556 }
1557 
1558 void
1559 MITC4Shell :: computeEdgeNMatrix(FloatMatrix &answer, int boundaryID, const FloatArray &lcoords)
1560 {
1561  FloatArray n_vec;
1562  this->giveInterpolation()->boundaryEdgeEvalN( n_vec, boundaryID, lcoords, FEIElementGeometryWrapper(this) );
1563  answer.beNMatrixOf(n_vec, 6);
1564 }
1565 
1566 
1567 void
1568 MITC4Shell :: computeSurfaceNMatrix(FloatMatrix &answer, int boundaryID, const FloatArray &lcoords)
1569 {
1570  FloatArray n_vec;
1571  this->giveInterpolation()->boundarySurfaceEvalN( n_vec, boundaryID, lcoords, FEIElementGeometryWrapper(this) );
1572  answer.beNMatrixOf(n_vec, 6);
1573 }
1574 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
void giveDirectorVectors(FloatArray &V1, FloatArray &V2, FloatArray &V3, FloatArray &V4)
Definition: mitc4.C:157
CrossSection * giveCrossSection()
Definition: element.C:495
Relative penalty stiffness for drilling DOFs.
Definition: crosssection.h:69
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
The element interface required by NodalAvergagingRecoveryModel.
virtual void computeConstitutiveMatrixAt(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
Computes constitutive matrix of receiver.
Definition: mitc4.C:331
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
void giveCharacteristicTensor(FloatMatrix &answer, CharTensor type, GaussPoint *gp, TimeStep *tStep)
Definition: mitc4.C:985
virtual void setupIRForMassMtrxIntegration(IntegrationRule &iRule)
Setup Integration Rule Gauss Points for Mass Matrix integration.
Definition: mitc4.C:1208
virtual void giveEdgeDofMapping(IntArray &answer, int iEdge) const
Assembles edge dof mapping mask, which provides mapping between edge local DOFs and "global" element ...
Definition: mitc4.C:1399
int number
Component number.
Definition: femcmpnn.h:80
The element interface required by ZZNodalRecoveryModel.
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Default implementation using Newton&#39;s method to find the local coordinates.
Definition: fei2dquadlin.C:120
virtual SPRPatchType SPRNodalRecoveryMI_givePatchType()
Definition: mitc4.C:138
virtual void giveRealStress_3dDegeneratedShell(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedStrain, TimeStep *tStep)
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
int nlGeometry
Flag indicating if geometrical nonlinearities apply.
void beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
Computes vector product (or cross product) of vectors given as parameters, , and stores the result in...
Definition: floatarray.C:415
Class and object Domain.
Definition: domain.h:115
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Returns local vector of unknowns.
Definition: element.C:86
virtual void giveSurfaceDofMapping(IntArray &answer, int iSurf) const
Assembles surface dof mapping mask, which provides mapping between surface local DOFs and "global" el...
Definition: mitc4.C:1502
virtual void printOutputAt(FILE *file, TimeStep *tStep)
Prints output of receiver to stream, for given time step.
Definition: mitc4.C:1029
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
Definition: element.h:822
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mitc4.C:388
double drillCoeff
Definition: mitc4.h:84
void giveLocalDirectorVectors(FloatArray &V1, FloatArray &V2, FloatArray &V3, FloatArray &V4)
Definition: mitc4.C:245
Director vector component in y-axis.
Definition: crosssection.h:75
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
Wrapper around cell with vertex coordinates stored in FloatArray**.
Definition: feinterpol.h:115
The element interface required by ZZNodalRecoveryModel.
Abstract base class for "structural" finite elements with geometrical nonlinearities.
void zero()
Sets all component to zero.
Definition: intarray.C:52
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
virtual void computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
Computes the stress vector of receiver at given integration point, at time step tStep.
Definition: mitc4.C:978
ConnectivityTable * giveConnectivityTable()
Returns receiver&#39;s associated connectivity table.
Definition: domain.C:1170
virtual Interface * giveInterface(InterfaceType interface)
Interface requesting service.
Definition: mitc4.C:84
const FloatMatrix * computeGtoLRotationMatrix()
Definition: mitc4.C:784
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge Jacobian of transformation between local and global coordinates.
Definition: feinterpol2d.C:175
void givedNdx(FloatArray &hkx, FloatArray &hky, FloatArray coords)
Definition: mitc4.C:1164
virtual void boundarySurfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the array of edge interpolation functions (shape functions) at given point.
virtual void SPRNodalRecoveryMI_giveDofMansDeterminedByPatch(IntArray &answer, int pap)
Definition: mitc4.C:111
virtual void giveDofManDofIDMask(int inode, IntArray &) const
Returns dofmanager dof mask for node.
Definition: mitc4.C:404
#define _IFT_NLStructuralElement_nlgeoflag
This class implements an quad element based on Mixed Interpolation of Tensorial Components (MITC)...
Definition: mitc4.h:71
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
virtual double giveCoordinate(int i)
Definition: node.C:82
void computeLocalBaseVectors(FloatArray &e1, FloatArray &e2, FloatArray &e3)
Definition: mitc4.C:814
static FEI2dQuadLin interp_lin
Element geometry approximation.
Definition: mitc4.h:77
virtual void giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
Evaluates nodal representation of real internal forces.
Definition: mitc4.C:550
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 giveNodeCoordinates(double &x1, double &x2, double &x3, double &x4, double &y1, double &y2, double &y3, double &y4, double &z1, double &z2, double &z3, double &z4)
Definition: mitc4.C:340
FloatMatrix GtoLRotationMatrix
Transformation Matrix form GtoL(3,3) is stored at the element level for computation efficiency...
Definition: mitc4.h:82
Abstract base class representing integration rule.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
virtual double computeEdgeVolumeAround(GaussPoint *gp, int iEdge)
Computes volume related to integration point on local edge.
Definition: mitc4.C:1424
void setSubMatrix(const FloatMatrix &src, int sr, int sc)
Adds the given matrix as sub-matrix to receiver.
Definition: floatmatrix.C:557
void beColumnOf(const FloatMatrix &mat, int col)
Reciever will be set to a given column in a matrix.
Definition: floatarray.C:1114
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
bool isNotEmpty() const
Tests for empty matrix.
Definition: floatmatrix.h:162
double giveNaturalCoordinate(int i) const
Returns i-th natural element coordinate of receiver.
Definition: gausspoint.h:136
void computeLToDirectorRotationMatrix(FloatMatrix &answer1, FloatMatrix &answer2, FloatMatrix &answer3, FloatMatrix &answer4)
Definition: mitc4.C:858
virtual void computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
Computes the stiffness matrix of receiver.
Definition: mitc4.C:493
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
Definition: mitc4.C:1214
Director vector component in z-axis.
Definition: crosssection.h:76
virtual void computeEdgeNMatrix(FloatMatrix &answer, int boundaryID, const FloatArray &lcoords)
computes edge interpolation matrix
Definition: mitc4.C:1559
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
void giveLocalCoordinates(FloatArray &answer, FloatArray &global)
Definition: mitc4.C:368
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
int computeLoadGToLRotationMtrx(FloatMatrix &answer)
Returns transformation matrix from global coordinate system to local element coordinate system for el...
Definition: mitc4.C:1304
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
Definition: fei2dquadlin.C:295
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
DofIDItem
Type representing particular dof type.
Definition: dofiditem.h:86
int giveNumber()
Returns number of receiver.
Definition: gausspoint.h:184
virtual void computeGaussPoints()
Initializes the array of integration rules member variable.
Definition: mitc4.C:145
virtual double giveWeight()
Returns integration weight of receiver.
Definition: gausspoint.h:181
virtual integrationDomain giveIntegrationDomain() const
Returns integration domain for receiver, used to initialize integration point over receiver volume...
Definition: mitc4.h:149
StructuralCrossSection * giveStructuralCrossSection()
Helper function which returns the structural cross-section for the element.
void giveThickness(double &a1, double &a2, double &a3, double &a4)
Definition: mitc4.C:768
virtual void NodalAveragingRecoveryMI_computeNodalValue(FloatArray &answer, int node, InternalStateType type, TimeStep *tStep)
Computes the element value in given node.
Definition: mitc4.C:1324
Class representing connectivity table.
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
virtual void boundaryEdgeEvalN(FloatArray &answer, int boundary, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates the basis functions on the requested boundary.
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
virtual void computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Compute strain vector of receiver evaluated at given integration point at given time step from elemen...
#define N(p, q)
Definition: mdm.C:367
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
virtual IntegrationRule * GetSurfaceIntegrationRule(int approxOrder)
Definition: mitc4.C:1517
void giveMidplaneIPValue(FloatArray &answer, int gpXY, InternalStateType type, TimeStep *tStep)
Definition: mitc4.C:1086
MITC4Shell(int n, Domain *d)
Definition: mitc4.C:60
const FloatArray & giveGlobalCoordinates()
Definition: gausspoint.h:160
virtual int setupIntegrationPoints(IntegrationRule &irule, int npoints, Element *element)
Sets up integration rule for the given element.
Definition: crosssection.C:54
int numberOfGaussPoints
Number of integration points as specified by nip.
Definition: element.h:188
void beNMatrixOf(const FloatArray &n, int nsd)
Assigns the receiver to be a repeated diagonal matrix.
Definition: floatmatrix.C:505
virtual int computeGlobalCoordinates(FloatArray &answer, const FloatArray &lcoords)
Computes the global coordinates from given element&#39;s local coordinates.
Definition: mitc4.C:1286
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
#define _IFT_MITC4Shell_directorType
Definition: mitc4.h:48
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
void giveNodeNeighbourList(IntArray &answer, IntArray &nodeList)
Returns list of elements sharing given nodes.
virtual double give(CrossSectionProperty a, GaussPoint *gp)
Returns the value of cross section property at given point.
Definition: crosssection.C:151
virtual void computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
Computes the stiffness matrix of receiver.
static void transformStrainVectorTo(FloatArray &answer, const FloatMatrix &base, const FloatArray &strainVector, bool transpose=false)
Transforms 3d strain vector into another coordinate system.
void beSubArrayOf(const FloatArray &src, const IntArray &indx)
Extract sub vector form src array and stores the result into receiver.
Definition: floatarray.C:379
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 resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual int computeLoadLEToLRotationMatrix(FloatMatrix &answer, int iEdge, GaussPoint *gp)
Returns transformation matrix from local edge c.s to element local coordinate system of load vector c...
Definition: mitc4.C:1458
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
Definition: mitc4.C:412
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
void giveJacobian(FloatArray lcoords, FloatMatrix &jacobianMatrix)
Definition: mitc4.C:433
Class Interface.
Definition: interface.h:82
int giveNumberOfIntegrationPoints() const
Returns number of integration points of receiver.
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void SPRNodalRecoveryMI_giveSPRAssemblyPoints(IntArray &pap)
Definition: mitc4.C:101
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
virtual int SPRNodalRecoveryMI_giveNumberOfIP()
Definition: mitc4.C:131
virtual void computeNmatrixAt(const FloatArray &iLocCoord, FloatMatrix &answer)
Computes interpolation matrix for element unknowns.
Definition: mitc4.C:259
The spatial localizer element interface associated to spatial localizer.
virtual void give3dDegeneratedShellStiffMtrx(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
Method for computing 3d shell stiffness matrix on degenerated shell elements.
Abstract base class for all "structural" constitutive models.
std::vector< std::unique_ptr< IntegrationRule > > integrationRulesArray
List of integration rules of receiver (each integration rule contains associated integration points a...
Definition: element.h:170
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
int setUpIntegrationPoints(integrationDomain intdomain, int nPoints, MaterialMode matMode)
Initializes the receiver.
virtual FloatArray * giveCoordinates()
Definition: node.h:114
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Domain * giveDomain() const
Definition: femcmpnn.h:100
Type of artificially added drilling stiffness for drilling DOFs.
Definition: crosssection.h:70
InterfaceType
Enumerative type, used to identify interface type.
Definition: interfacetype.h:43
virtual void giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord=0)
Evaluates nodal representation of real internal forces.
int directorType
Definition: mitc4.h:83
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
virtual void computeSurfaceNMatrixAt(FloatMatrix &answer, int iSurf, GaussPoint *sgp)
Definition: mitc4.C:1491
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
Abstract base class for all structural cross section models.
virtual void computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int=1, int=ALL_STRAINS)
Computes the geometrical matrix of receiver in given integration point.
Definition: mitc4.C:597
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
Definition: fei2dquadlin.C:350
int giveLabel() const
Definition: element.h:1055
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
virtual bool computeLocalCoordinates(FloatArray &answer, const FloatArray &coords)
Computes the element local coordinates from given global coordinates.
Definition: mitc4.C:1259
int giveNumber() const
Definition: femcmpnn.h:107
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
double normalize()
Normalizes receiver.
Definition: floatarray.C:828
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
static void transformStressVectorTo(FloatArray &answer, const FloatMatrix &base, const FloatArray &stressVector, bool transpose=false)
Transforms 3d stress vector into another coordinate system.
Director vector component in x-axis.
Definition: crosssection.h:74
Class representing integration point in finite element program.
Definition: gausspoint.h:93
Class representing solution step.
Definition: timestep.h:80
virtual Material * giveMaterial(IntegrationPoint *ip)
Returns the material associated with the GP.
#define _IFT_MITC4Shell_nipZ
Definition: mitc4.h:47
int numberOfDofMans
Number of dofmanagers.
Definition: element.h:149
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
virtual MaterialMode giveMaterialMode()
Returns material mode for receiver integration points.
Definition: mitc4.h:150
virtual double computeSurfaceVolumeAround(GaussPoint *gp, int iSurf)
Computes volume related to integration point on local surface.
Definition: mitc4.C:1528
virtual void computeSurfaceNMatrix(FloatMatrix &answer, int boundaryID, const FloatArray &lcoords)
Computes surface interpolation matrix.
Definition: mitc4.C:1568
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: fei2dquadlin.C:59
const FloatArray & giveNaturalCoordinates()
Returns coordinate array of receiver.
Definition: gausspoint.h:138
virtual FEInterpolation * giveInterpolation() const
Definition: mitc4.C:73
#define _IFT_Element_nip
Definition: element.h:72
Class representing Gaussian-quadrature integration rule.
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:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011