OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
igaelements.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/igaelements.h"
36 #include "../sm/CrossSections/structuralcrosssection.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "domain.h"
40 #include "node.h"
41 #include "element.h"
42 #include "gausspoint.h"
43 #include "gaussintegrationrule.h"
44 #include "matresponsemode.h"
45 #include "crosssection.h"
46 #include "mathfem.h"
47 #include "iga/iga.h"
48 #include "classfactory.h"
49 
50 #ifdef __OOFEG
51  #include "oofeggraphiccontext.h"
52  #include "oofegutils.h"
53 #endif
54 
55 namespace oofem {
56 REGISTER_Element(BsplinePlaneStressElement);
57 REGISTER_Element(NURBSPlaneStressElement);
58 REGISTER_Element(TSplinePlaneStressElement);
59 REGISTER_Element(NURBSSpace3dElement);
60 
61 
63 
64 
66 {
67  //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
69 }
70 
71 
73 {
74  BSplineInterpolation *interpol = static_cast< BSplineInterpolation * >( this->giveInterpolation() );
75  if ( giveNumberOfDofManagers() != interpol->giveNumberOfControlPoints(1) * interpol->giveNumberOfControlPoints(2) ) {
76  OOFEM_WARNING("number of control points mismatch");
77  return 0;
78  }
79  return 1;
80 }
81 
82 
83 
85 
86 
88 {
89  //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
91 }
92 
93 
95 {
96  NURBSInterpolation *interpol = static_cast< NURBSInterpolation * >( this->giveInterpolation() );
97  if ( giveNumberOfDofManagers() != interpol->giveNumberOfControlPoints(1) * interpol->giveNumberOfControlPoints(2) ) {
98  OOFEM_WARNING("number of control points mismatch");
99  return 0;
100  }
101  return 1;
102 }
103 
104 
106 
107 
108 
109 
111 
112 
114 {
115  //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
116  return IGAElement :: initializeFrom(ir);
117 }
118 
119 
121 {
122  NURBSInterpolation *interpol = static_cast< NURBSInterpolation * >( this->giveInterpolation() );
123  if ( giveNumberOfDofManagers() != interpol->giveNumberOfControlPoints(1) * interpol->giveNumberOfControlPoints(2) * interpol->giveNumberOfControlPoints(3) ) {
124  OOFEM_WARNING("number of control points mismatch");
125  return 0;
126  }
127  return 1;
128 }
129 
130 // HUHU should be implemented by IGA element (it is the same for Bspline NURBS and TSpline)
131 // however in such a case it should be generalized in terms of appropriately multiplying
132 // nseq for those integration elements which span more tham just a single knot span
133 // the reason is to ensure compatible division to quads over which scalar quantity is interpolated
134 // bilinearly !!!
135 
136 #ifdef __OOFEG
137 
138 //#define COMPUTE_STRESS
139  #define COMPUTE_STRAIN
140 
142 {
143  int indx;
144  WCRec p [ 4 ];
145  GraphicObj *go;
146  double s [ 4 ];
147  FEInterpolation *interp = this->giveInterpolation();
148  FloatArray val;
149 
150  if ( !gc.testElementGraphicActivity(this) ) {
151  return;
152  }
153 
154  EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
155  EASValsSetFillStyle(FILL_SOLID);
156  const double *const *knotVector = interp->giveKnotVector();
157  const IntArray *span;
158  int j, nsd = this->giveNsd();
159  FloatArray c [ 4 ], cg [ 4 ], u;
160  IntArray sign [ 4 ];
161 
162  if ( nsd == 2 ) {
163  for ( j = 0; j < 4; j++ ) {
164  c [ j ].resize(2);
165  cg [ j ].resize(2);
166  sign [ j ].resize(2);
167  }
168  sign [ 0 ].at(1) = 1;
169  sign [ 0 ].at(2) = 1;
170  sign [ 1 ].at(1) = -1;
171  sign [ 1 ].at(2) = 1;
172  sign [ 2 ].at(1) = -1;
173  sign [ 2 ].at(2) = -1;
174  sign [ 3 ].at(1) = 1;
175  sign [ 3 ].at(2) = -1;
176  } else {
177  OOFEM_ERROR("not implemented for nsd = %d", nsd);
178  }
179 
180  indx = gc.giveIntVarIndx();
181 
183 
184  // loop over individual integration rules (i.e., knot spans)
185  for ( auto &iRule: integrationRulesArray ) {
186  span = iRule->giveKnotSpan();
187  if ( nsd == 2 ) {
188  // divide span locally to get finer geometry rep.
189  int i, j, k, nseg = 4;
190  double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
191  double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
192  for ( i = 1; i <= nseg; i++ ) {
193  for ( j = 1; j <= nseg; j++ ) {
194  c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
195  c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
196  c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
197  c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
198  c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
199  c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
200  c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
201  c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
202 
203  for ( k = 0; k < 4; k++ ) {
204  interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
205  p [ k ].x = ( FPNum ) cg [ k ].at(1);
206  p [ k ].y = ( FPNum ) cg [ k ].at(2);
207  p [ k ].z = 0.;
208 
209  #ifdef QUARTER_PLATE_WITH_HOLE_SINGLE_PATCH
210  //move sampling gp out of boundary to overcome degeneracy on quarter plate with hole modelled by single patch
211  if ( c [ k ].at(1) > 0.99999 && c [ k ].at(2) > 0.495 && c [ k ].at(2) < 0.505 ) {
212  c [ k ].at(1) += sign [ k ].at(1) * du / 10.0;
213  c [ k ].at(2) += sign [ k ].at(2) * dv / 10.0;
214  c [ k ].at(3) += sign [ k ].at(3) * dw / 10.0;
215  }
216  #endif
217 
218  // create a dummy ip's
219  GaussPoint gp(iRule.get(), 999, c [ k ], 1.0, _PlaneStress);
220  #ifdef COMPUTE_STRAIN
221  this->computeStrainVector(val, & gp, tStep, u);
222  #endif
223  #ifdef COMPUTE_STRESS
224  FloatArray strain;
225  this->computeStrainVector(strain, & gp, tStep, u);
226  this->computeStressVector(val, strain, & gp, tStep);
227  #endif
228  s [ k ] = val.at(indx);
229  }
230 
231  if ( ( std::isnan(s [ 0 ]) ) || ( std::isnan(s [ 1 ]) ) || ( std::isnan(s [ 2 ]) ) || ( std::isnan(s [ 3 ]) ) ) {
232  continue;
233  }
234 
235  if ( ( fabs(s [ 0 ]) > 1.e5 ) || ( fabs(s [ 1 ]) > 1.e5 ) || ( fabs(s [ 2 ]) > 1.e5 ) || ( fabs(s [ 3 ]) > 1.e5 ) ) {
236  continue;
237  }
238 
239  //printf ("QWD: %e %e %e %e\n", s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
240 
241  go = CreateQuadWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
242  EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
243  EGAttachObject(go, ( EObjectP ) this);
244  EMAddGraphicsToModel(ESIModel(), go);
245  }
246  }
247  }
248  } // end loop over knot spans (irules)
249 }
250 
251 
253 {
254  int indx;
255  WCRec p [ 4 ];
256  GraphicObj *go;
257  double s [ 4 ];
258  FEInterpolation *interp = this->giveInterpolation();
259  FloatArray val, u;
260 
261  if ( !gc.testElementGraphicActivity(this) ) {
262  return;
263  }
264 
265  EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
266  EASValsSetFillStyle(FILL_SOLID);
267  EASValsSetEdgeFlag(true);
268  const double *const *knotVector = interp->giveKnotVector();
269  const IntArray *span;
270  int j, nsd = this->giveNsd();
271  FloatArray c [ 4 ], cg [ 4 ];
272  IntArray sign [ 4 ];
273 
274  if ( nsd == 2 ) {
275  for ( j = 0; j < 4; j++ ) {
276  c [ j ].resize(2);
277  cg [ j ].resize(2);
278  sign [ j ].resize(2);
279  }
280  sign [ 0 ].at(1) = 1;
281  sign [ 0 ].at(2) = 1;
282  sign [ 1 ].at(1) = -1;
283  sign [ 1 ].at(2) = 1;
284  sign [ 2 ].at(1) = -1;
285  sign [ 2 ].at(2) = -1;
286  sign [ 3 ].at(1) = 1;
287  sign [ 3 ].at(2) = -1;
288  } else {
289  OOFEM_ERROR("not implemented for nsd = %d", nsd);
290  }
291 
292  indx = gc.giveIntVarIndx();
293 
295 
296  //double maxs=-1.0e10, mins=1.0e10;
297 
298  // loop over individual integration rules (i.e., knot spans)
299  for ( auto &iRule: integrationRulesArray ) {
300  span = iRule->giveKnotSpan();
301  if ( nsd == 2 ) {
302  // divide span locally to get finer geometry rep.
303  int i, j, k, nseg = 8;
304  double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
305  double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
306  for ( i = 1; i <= nseg; i++ ) {
307  for ( j = 1; j <= nseg; j++ ) {
308  c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
309  c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
310  c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
311  c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
312  c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
313  c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
314  c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
315  c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
316 
317  for ( k = 0; k < 4; k++ ) {
318  interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
319  p [ k ].x = ( FPNum ) cg [ k ].at(1);
320  p [ k ].y = ( FPNum ) cg [ k ].at(2);
321  p [ k ].z = 0.;
322 
323  //move sampling gp out of boundary to overcome degeneracy on quarter plate with hole modelled by single patch
324  if ( c [ k ].at(1) > 0.99999 && c [ k ].at(2) > 0.495 && c [ k ].at(2) < 0.505 ) {
325  c [ k ].at(1) += sign [ k ].at(1) * du / 10.0;
326  c [ k ].at(2) += sign [ k ].at(2) * dv / 10.0;
327  }
328 
329  // create a dummy ip's
330  GaussPoint gp(iRule.get(), 999, c [ k ], 1.0, _PlaneStress);
331  #ifdef COMPUTE_STRAIN
332  this->computeStrainVector(val, & gp, tStep, u);
333  #endif
334  #ifdef COMPUTE_STRESS
335  FloatArray strain;
336  this->computeStrainVector(strain, & gp, tStep, u);
337  this->computeStressVector(val, strain, & gp, tStep);
338  #endif
339  s [ k ] = val.at(indx);
340 
341  #ifdef QUARTER_PLATE_WITH_HOLE
342  if ( indx == 2 ) {
343  if ( cg [ k ].at(1) <= 1.0 + 1.0e-10 && cg [ k ].at(2) <= 1.0e-10 ) {
344  fprintf(stderr, "A: syy = %e\n", s [ k ]);
345  }
346  }
347  if ( indx == 1 ) {
348  if ( cg [ k ].at(1) <= 1.0e-10 && cg [ k ].at(2) <= 1.0 + 1.0e-10 ) {
349  fprintf(stderr, "B: sxx = %e\n", s [ k ]);
350  }
351  }
352 
353  double x, y, r, phi, rate, E, G, kap, ny;
354  x = cg [ k ].at(1);
355  y = cg [ k ].at(2);
356  if ( x < 1.0e-10 ) {
357  phi = M_PI / 2.0;
358  r = y;
359  } else {
360  phi = atan(y / x);
361  r = x / cos(phi);
362  }
363 
364  #ifdef STRESS
365  // exact stresses quarter plate with hole s0=1 a=1
366  rate = 1.0 / r;
367  rate *= rate;
368  if ( indx == 1 ) {
369  s [ k ] = 0.5 * ( 2.0 + 3.0 * rate * rate * cos(4.0 * phi) - rate * ( 3 * cos(2.0 * phi) + 2.0 * cos(4.0 * phi) ) );
370  }
371  if ( indx == 2 ) {
372  s [ k ] = 0.5 * ( -3.0 * rate * rate * cos(4.0 * phi) - rate * ( cos(2.0 * phi) - 2.0 * cos(4.0 * phi) ) );
373  }
374  if ( indx == 3 ) {
375  s [ k ] = 0.5 * ( 3.0 * rate * rate * sin(4.0 * phi) - rate * ( sin(2.0 * phi) + 2.0 * sin(4.0 * phi) ) );
376  }
377 
378  if ( indx == 2 ) {
379  if ( cg [ k ].at(1) <= 1.0 + 1.0e-10 && cg [ k ].at(2) <= 1.0e-10 ) {
380  fprintf(stderr, "A: syy = %e\n", s [ k ]);
381  }
382  }
383  if ( indx == 1 ) {
384  if ( cg [ k ].at(1) <= 1.0e-10 && cg [ k ].at(2) <= 1.0 + 1.0e-10 ) {
385  fprintf(stderr, "B: sxx = %e\n", s [ k ]);
386  }
387  }
388  #endif
389  #ifdef DISPL
390  // exact displ quarter plate with hole s0=1 a=1
391  E = 15000.0;
392  ny = 0.25;
393  G = E / ( 2.0 * ( 1.0 + ny ) );
394  kap = ( 3.0 - ny ) / ( 1.0 + ny );
395  rate = 1.0 / r;
396  if ( indx == 1 ) {
397  s [ k ] = ( r * ( kap + 1.0 ) * cos(phi) + 2.0 * rate * ( ( 1.0 + kap ) * cos(phi) + cos(3.0 * phi) ) - 2.0 * rate * rate * rate * cos(3.0 * phi) ) / ( 8.0 * G );
398  }
399  if ( indx == 2 ) {
400  s [ k ] = ( r * ( kap - 3.0 ) * sin(phi) + 2.0 * rate * ( ( 1.0 - kap ) * sin(phi) + sin(3.0 * phi) ) - 2.0 * rate * rate * rate * sin(3.0 * phi) ) / ( 8.0 * G );
401  }
402 
403  if ( indx == 1 ) {
404  if ( cg [ k ].at(1) <= 1.0 + 1.0e-10 && cg [ k ].at(2) <= 1.0e-10 ) {
405  fprintf(stderr, "A: ux = %e\n", s [ k ]);
406  }
407  }
408  if ( indx == 2 ) {
409  if ( cg [ k ].at(1) <= 1.0e-10 && cg [ k ].at(2) <= 1.0 + 1.0e-10 ) {
410  fprintf(stderr, "B: uy = %e\n", s [ k ]);
411  }
412  }
413  #endif
414 
415  if ( s [ k ] < mins ) {
416  mins = s [ k ];
417  }
418  if ( s [ k ] > maxs ) {
419  maxs = s [ k ];
420  }
421  #endif
422  }
423 
424  if ( ( std::isnan(s [ 0 ]) ) || ( std::isnan(s [ 1 ]) ) || ( std::isnan(s [ 2 ]) ) || ( std::isnan(s [ 3 ]) ) ) {
425  continue;
426  }
427 
428  if ( ( fabs(s [ 0 ]) > 1.e5 ) || ( fabs(s [ 1 ]) > 1.e5 ) || ( fabs(s [ 2 ]) > 1.e5 ) || ( fabs(s [ 3 ]) > 1.e5 ) ) {
429  continue;
430  }
431 
432  //printf ("QWD: %e %e %e %e\n", s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
433 
434  go = CreateQuadWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
435  EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
436  EGAttachObject(go, ( EObjectP ) this);
437  EMAddGraphicsToModel(ESIModel(), go);
438  }
439  }
440  }
441  } // end loop over knot spans (irules)
442 
443  #ifdef QUARTER_PLATE_WITH_HOLE
444  fprintf(stderr, "%d: %e %e %e %e\n", indx, mins, maxs, ( 10.0 * mins + maxs ) / 11.0, ( 10.0 * maxs + mins ) / 11.0);
445  #endif
446 }
447 
448 
449 // refinement of integration elements should be generalized in terms of appropriately multiplying
450 // nseq for those integration elements which span more tham just a single knot span
451 // the reason is to ensure compatible division to quads over which scalar quantity is interpolated
452 // bilinearly !!!
453 
455 {
456  int indx;
457  WCRec p [ 4 ];
458  GraphicObj *go;
459  double s [ 4 ];
460  FEInterpolation *interp = this->giveInterpolation();
461  FloatArray val;
462 
463  if ( !gc.testElementGraphicActivity(this) ) {
464  return;
465  }
466 
467  EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
468  EASValsSetFillStyle(FILL_SOLID);
469  EASValsSetEdgeFlag(true);
470  const double *const *knotVector = interp->giveKnotVector();
471  const IntArray *span;
472  int j, nsd = this->giveNsd();
473  FloatArray c [ 4 ], cg [ 4 ], u;
474  IntArray sign [ 4 ];
475 
476  if ( nsd == 2 ) {
477  for ( j = 0; j < 4; j++ ) {
478  c [ j ].resize(2);
479  cg [ j ].resize(2);
480  sign [ j ].resize(2);
481  }
482  sign [ 0 ].at(1) = 1;
483  sign [ 0 ].at(2) = 1;
484  sign [ 1 ].at(1) = -1;
485  sign [ 1 ].at(2) = 1;
486  sign [ 2 ].at(1) = -1;
487  sign [ 2 ].at(2) = -1;
488  sign [ 3 ].at(1) = 1;
489  sign [ 3 ].at(2) = -1;
490  } else {
491  OOFEM_ERROR("not implemented for nsd = %d", nsd);
492  }
493 
494  indx = gc.giveIntVarIndx();
495 
497 
498  // loop over individual integration rules (i.e., knot spans)
499  for ( auto &iRule: integrationRulesArray ) {
500  span = iRule->giveKnotSpan();
501  if ( nsd == 2 ) {
502  // divide span locally to get finer geometry rep.
503  int i, j, k, nseg = 4;
504  double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
505  double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
506  for ( i = 1; i <= nseg; i++ ) {
507  for ( j = 1; j <= nseg; j++ ) {
508  c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
509  c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
510  c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
511  c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
512  c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
513  c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
514  c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
515  c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
516 
517  for ( k = 0; k < 4; k++ ) {
518  interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
519  p [ k ].x = ( FPNum ) cg [ k ].at(1);
520  p [ k ].y = ( FPNum ) cg [ k ].at(2);
521  p [ k ].z = 0.;
522 
523  #ifdef QUARTER_PLATE_WITH_HOLE_SINGLE_PATCH
524  //move sampling gp out of boundary to overcome degeneracy on quarter plate with hole modelled by single patch
525  if ( c [ k ].at(1) > 0.99999 && c [ k ].at(2) > 0.495 && c [ k ].at(2) < 0.505 ) {
526  c [ k ].at(1) += sign [ k ].at(1) * du / 10.0;
527  c [ k ].at(2) += sign [ k ].at(2) * dv / 10.0;
528  }
529  #endif
530 
531  // create a dummy ip's
532  GaussPoint gp(iRule.get(), 999, c [ k ], 1.0, _PlaneStress);
533  #ifdef COMPUTE_STRAIN
534  this->computeStrainVector(val, & gp, tStep, u);
535  #endif
536  #ifdef COMPUTE_STRESS
537  FloatArray strain;
538  this->computeStrainVector(strain, & gp, tStep, u);
539  this->computeStressVector(val, strain, & gp, tStep);
540  #endif
541  s [ k ] = val.at(indx);
542  }
543 
544  if ( ( std::isnan(s [ 0 ]) ) || ( std::isnan(s [ 1 ]) ) || ( std::isnan(s [ 2 ]) ) || ( std::isnan(s [ 3 ]) ) ) {
545  continue;
546  }
547 
548  if ( ( fabs(s [ 0 ]) > 1.e5 ) || ( fabs(s [ 1 ]) > 1.e5 ) || ( fabs(s [ 2 ]) > 1.e5 ) || ( fabs(s [ 3 ]) > 1.e5 ) ) {
549  continue;
550  }
551 
552  //printf ("QWD: %e %e %e %e\n", s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
553 
554  go = CreateQuadWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ]);
555  EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
556  EGAttachObject(go, ( EObjectP ) this);
557  EMAddGraphicsToModel(ESIModel(), go);
558  }
559  }
560  }
561  } // end loop over knot spans (irules)
562 }
563 
564 
565 
567 {
568  int indx;
569  WCRec p [ 8 ];
570  GraphicObj *go;
571  double s [ 8 ];
572  FEInterpolation *interp = this->giveInterpolation();
573  FloatArray val, u;
574  //int huhu = 0;
575 
576  if ( !gc.testElementGraphicActivity(this) ) {
577  return;
578  }
579 
580  EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
581  EASValsSetFillStyle(FILL_SOLID);
582  EASValsSetEdgeFlag(true);
583  const double *const *knotVector = interp->giveKnotVector();
584  const IntArray *span;
585  int j, nsd = this->giveNsd();
586  FloatArray c [ 8 ], cg [ 8 ];
587  IntArray sign [ 8 ];
588 
589  if ( nsd == 3 ) {
590  for ( j = 0; j < 8; j++ ) {
591  c [ j ].resize(3);
592  cg [ j ].resize(3);
593  sign [ j ].resize(3);
594  }
595 
596  sign [ 0 ].at(1) = 1;
597  sign [ 0 ].at(2) = 1;
598  sign [ 0 ].at(3) = 1;
599  sign [ 1 ].at(1) = -1;
600  sign [ 1 ].at(2) = 1;
601  sign [ 1 ].at(3) = 1;
602  sign [ 2 ].at(1) = -1;
603  sign [ 2 ].at(2) = -1;
604  sign [ 2 ].at(3) = 1;
605  sign [ 3 ].at(1) = 1;
606  sign [ 3 ].at(2) = -1;
607  sign [ 3 ].at(3) = 1;
608  sign [ 4 ].at(1) = 1;
609  sign [ 4 ].at(2) = 1;
610  sign [ 4 ].at(3) = -1;
611  sign [ 5 ].at(1) = -1;
612  sign [ 5 ].at(2) = 1;
613  sign [ 5 ].at(3) = -1;
614  sign [ 6 ].at(1) = -1;
615  sign [ 6 ].at(2) = -1;
616  sign [ 6 ].at(3) = -1;
617  sign [ 7 ].at(1) = 1;
618  sign [ 7 ].at(2) = -1;
619  sign [ 7 ].at(3) = -1;
620  } else {
621  OOFEM_ERROR("not implemented for nsd = %d", nsd);
622  }
623 
624  indx = gc.giveIntVarIndx();
625 
626  #ifdef SPHERICAL_CS
627  if ( gc.giveIntVarType() == IST_StrainTensor ) {
628  huhu = 1;
629  }
630  #endif
631  #ifdef MISSES_STRESS
632  if ( gc.giveIntVarType() == IST_ErrorIndicatorLevel ) {
633  huhu = 2;
634  }
635  #endif
636 
638 
639  //double maxs=-1.0e10, mins=1.0e10;
640 
641  // loop over individual integration rules (i.e., knot spans)
642  for ( auto &iRule: integrationRulesArray ) {
643  span = iRule->giveKnotSpan();
644  if ( nsd == 3 ) {
645  // divide span locally to get finer geometry rep.
646  int i, j, k, m, nseg = 8;
647  double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
648  double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
649  double dw = ( knotVector [ 2 ] [ span->at(3) + 1 ] - knotVector [ 2 ] [ span->at(3) ] ) / nseg;
650 
651  #ifdef DRAW_VISIBLE_CONTOUR
652  WCRec pp [ 4 ];
653  double ss [ 4 ];
654  int kk, ii, nd [ 6 ] [ 4 ] = { { 0, 1, 2, 3 }, { 4, 5, 6, 7 }, { 0, 1, 5, 4 }, { 1, 2, 6, 5 }, { 2, 3, 7, 6 }, { 3, 0, 4, 7 } };
655  #endif
656 
657  for ( i = 1; i <= nseg; i++ ) {
658  for ( j = 1; j <= nseg; j++ ) {
659  for ( m = 1; m <= nseg; m++ ) {
660  c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
661  c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
662  c [ 0 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * ( m - 1 );
663  c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
664  c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
665  c [ 1 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * ( m - 1 );
666  c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
667  c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
668  c [ 2 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * ( m - 1 );
669  c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
670  c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
671  c [ 3 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * ( m - 1 );
672  c [ 4 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
673  c [ 4 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
674  c [ 4 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * m;
675  c [ 5 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
676  c [ 5 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
677  c [ 5 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * m;
678  c [ 6 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
679  c [ 6 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
680  c [ 6 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * m;
681  c [ 7 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
682  c [ 7 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
683  c [ 7 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dw * m;
684 
685  for ( k = 0; k < 8; k++ ) {
686  interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
687  p [ k ].x = ( FPNum ) cg [ k ].at(1);
688  p [ k ].y = ( FPNum ) cg [ k ].at(2);
689  p [ k ].z = ( FPNum ) cg [ k ].at(3);
690  }
691 
692  #ifdef DRAW_VISIBLE_CONTOUR
693  #ifdef SPHERE_WITH_HOLE
694 
695  // check whether drawing side visible in particular view !!!
696  int haha = 0;
697  for ( kk = 0; kk < 6; kk++ ) {
698  double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
699  for ( ii = 0; ii < 4; ii++ ) {
700  xx += ( pp [ ii ].x = p [ nd [ kk ] [ ii ] ].x );
701  yy += ( pp [ ii ].y = p [ nd [ kk ] [ ii ] ].y );
702  zz += ( pp [ ii ].z = p [ nd [ kk ] [ ii ] ].z );
703  ss [ ii ] = s [ nd [ kk ] [ ii ] ];
704  }
705  xx /= 4.0;
706  yy /= 4.0;
707  zz /= 4.0;
708  rr = xx * xx + yy * yy;
709  r = rr + zz * zz;
710 
711  if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.0 || r < 25.0 || r > 5.99 * 5.99 ) {
712  haha = 1;
713  }
714  }
715  if ( haha == 0 ) {
716  continue;
717  }
718  #endif
719  #endif
720 
721  for ( k = 0; k < 8; k++ ) {
722  // create a dummy ip's
723  GaussPoint gp(iRule.get(), 999, c [ k ], 1.0, _3dMat);
724  #ifdef COMPUTE_STRAIN
725  this->computeStrainVector(val, & gp, tStep, u);
726  #endif
727  #ifdef COMPUTE_STRESS
728  FloatArray strain;
729  this->computeStrainVector(strain, & gp, tStep, u);
730  this->computeStressVector(val, strain, & gp, tStep);
731  #endif
732  s [ k ] = val.at(indx);
733 
734  #ifdef MISSES_STRESS
735  if ( huhu == 2 ) {
736  double vonMisses;
737  vonMisses = sqrt( ( ( val.at(1) - val.at(2) ) * ( val.at(1) - val.at(2) ) + ( val.at(2) - val.at(3) ) * ( val.at(2) - val.at(3) ) + ( val.at(1) - val.at(3) ) * ( val.at(1) - val.at(3) ) + 6.0 * ( val.at(4) * val.at(4) + val.at(5) * val.at(5) + val.at(6) * val.at(6) ) ) / 2.0 );
738  s [ k ] = vonMisses;
739  }
740  #endif
741  #ifdef SPHERICAL_CS
742  if ( huhu == 1 ) {
743  double x, y, z, r, r2, rr;
744  FloatMatrix sigma(3, 3), T(3, 3), product(3, 3), result(3, 3);
745 
746  sigma.at(1, 1) = val.at(1);
747  sigma.at(1, 2) = val.at(6);
748  sigma.at(1, 3) = val.at(5);
749  sigma.at(2, 1) = val.at(6);
750  sigma.at(2, 2) = val.at(2);
751  sigma.at(2, 3) = val.at(4);
752  sigma.at(3, 1) = val.at(5);
753  sigma.at(3, 2) = val.at(4);
754  sigma.at(3, 3) = val.at(3);
755 
756  x = cg [ k ].at(1);
757  y = cg [ k ].at(2);
758  z = cg [ k ].at(3);
759  r2 = x * x + y * y + z * z;
760  r = sqrt(r2);
761  rr = sqrt(x * x + y * y);
762 
763  T.at(1, 1) = -z * z * y / rr / r2 - y * rr / r2;
764  T.at(1, 2) = z * z * x / rr / r2 + x * rr / r2;
765  T.at(1, 3) = 0.0;
766  T.at(2, 1) = -z * x / rr / r;
767  T.at(2, 2) = -z * y / rr / r;
768  T.at(2, 3) = rr / r;
769  T.at(3, 1) = x / r;
770  T.at(3, 2) = y / r;
771  T.at(3, 3) = z / r;
772 
773  product.beProductOf(T, sigma);
774  result.beProductTOf(product, T);
775 
776  val.at(1) = result.at(1, 1);
777  val.at(6) = result.at(1, 2);
778  val.at(5) = result.at(1, 3);
779  val.at(6) = result.at(2, 1);
780  val.at(2) = result.at(2, 2);
781  val.at(4) = result.at(2, 3);
782  val.at(5) = result.at(3, 1);
783  val.at(4) = result.at(3, 2);
784  val.at(3) = result.at(3, 3);
785 
786  s [ k ] = val.at(indx);
787  }
788  #endif
789  #ifdef SPHERE_WITH_HOLE
790  if ( s [ k ] < mins ) {
791  mins = s [ k ];
792  }
793  if ( s [ k ] > maxs ) {
794  maxs = s [ k ];
795  }
796  #endif
797  }
798 
799  if ( ( std::isnan(s [ 0 ]) ) || ( std::isnan(s [ 1 ]) ) || ( std::isnan(s [ 2 ]) ) || ( std::isnan(s [ 3 ]) ) ) {
800  continue;
801  }
802 
803  if ( ( std::isnan(s [ 4 ]) ) || ( std::isnan(s [ 5 ]) ) || ( std::isnan(s [ 6 ]) ) || ( std::isnan(s [ 7 ]) ) ) {
804  continue;
805  }
806 
807  if ( ( fabs(s [ 0 ]) > 1.e5 ) || ( fabs(s [ 1 ]) > 1.e5 ) || ( fabs(s [ 2 ]) > 1.e5 ) || ( fabs(s [ 3 ]) > 1.e5 ) ) {
808  continue;
809  }
810 
811  if ( ( fabs(s [ 4 ]) > 1.e5 ) || ( fabs(s [ 5 ]) > 1.e5 ) || ( fabs(s [ 6 ]) > 1.e5 ) || ( fabs(s [ 7 ]) > 1.e5 ) ) {
812  continue;
813  }
814 
815  //printf ("HWD: %e %e %e %e %e %e %e %e\n", s [ 0 ], s [ 1 ], s [ 2 ], s [ 3 ], s [ 4 ], s [ 5 ], s [ 6 ], s [ 7 ]);
816 
817  #ifdef DRAW_VISIBLE_CONTOUR
818  #ifdef SPHERE_WITH_HOLE
819  for ( kk = 0; kk < 6; kk++ ) {
820  double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
821  for ( ii = 0; ii < 4; ii++ ) {
822  xx += ( pp [ ii ].x = p [ nd [ kk ] [ ii ] ].x );
823  yy += ( pp [ ii ].y = p [ nd [ kk ] [ ii ] ].y );
824  zz += ( pp [ ii ].z = p [ nd [ kk ] [ ii ] ].z );
825  ss [ ii ] = s [ nd [ kk ] [ ii ] ];
826  }
827  xx /= 4.0;
828  yy /= 4.0;
829  zz /= 4.0;
830  rr = xx * xx + yy * yy;
831  r = rr + zz * zz;
832 
833  if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.0 || r < 25.0 || r > 5.99 * 5.99 ) {
834  go = CreateQuadWD3D(pp, ss [ 0 ], ss [ 1 ], ss [ 2 ], ss [ 3 ]);
835  EGWithMaskChangeAttributes(LAYER_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK, go);
836  EMAddGraphicsToModel(ESIModel(), go);
837  }
838  }
839  #endif
840  #else
841  go = CreateHexahedronWD(p, s);
842  EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
843  EGAttachObject(go, ( EObjectP ) this);
844  EMAddGraphicsToModel(ESIModel(), go);
845  #endif
846  }
847  }
848  }
849  }
850  } // end loop over knot spans (irules)
851 
852  #ifdef SPHERE_WITH_HOLE
853  fprintf(stderr, "%d %e %e %e %e\n", indx, mins, maxs, ( 10.0 * mins + maxs ) / 11.0, ( 10.0 * maxs + mins ) / 11.0);
854  #endif
855 }
856 
857 
858 #endif
859 } // end namespace oofem
BSplineInterpolation interpolation
Definition: igaelements.h:58
int testElementGraphicActivity(Element *)
Test if particular element passed fulfills various filtering criteria for its graphics output...
Interpolation for B-splines.
Definition: feibspline.h:60
virtual FEInterpolation * giveInterpolation() const
Definition: igaelements.h:73
Interpolation class for NURBS.
Definition: feinurbs.h:47
Class and object Domain.
Definition: domain.h:115
virtual int checkConsistency()
Performs consistency check.
Definition: igaelements.C:72
IGATSplineElement setups integration rules differently from IGAElement.
Definition: iga.h:123
General purpose 3d structural element evaluator.
virtual int giveNumberOfControlPoints(int dim)
Definition: feibspline.h:182
virtual FEInterpolation * giveInterpolation() const
Definition: igaelements.h:201
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual void drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
Definition: igaelements.C:454
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
bool isnan(double x)
Definition: mathfem.h:103
virtual void drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
Definition: igaelements.C:252
NURBSSpace3dElement(int n, Domain *aDomain)
Definition: igaelements.C:110
void computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, FloatArray &u)
Optimized version, allowing to pass element displacements as parameter.
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual int giveNumberOfDofManagers() const
Definition: element.h:656
virtual FEInterpolation * giveInterpolation() const
Definition: igaelements.h:115
Class representing a general abstraction for finite element interpolation class.
Definition: feinterpol.h:132
InternalStateType giveIntVarType()
Implements base IGAElement, supposed to be a parent class of all elements with B-spline or NURBS base...
Definition: iga.h:96
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: igaelements.C:87
NURBSPlaneStressElement(int n, Domain *aDomain)
Definition: igaelements.C:84
#define E(p)
Definition: mdm.C:368
#define M_PI
Definition: mathfem.h:52
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: igaelements.C:113
virtual void computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
Computes the stress vector.
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
TSplinePlaneStressElement(int n, Domain *aDomain)
Definition: igaelements.C:105
virtual int checkConsistency()
Performs consistency check.
Definition: igaelements.C:94
virtual void drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
Definition: igaelements.C:566
Class representing vector of real numbers.
Definition: floatarray.h:82
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
IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: iga.C:51
TSplineInterpolation interpolation
Definition: igaelements.h:145
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
virtual const double *const * giveKnotVector()
Returns the subdivision of patch parametric space.
Definition: feinterpol.h:457
Class representing the general Input Record.
Definition: inputrecord.h:101
BsplinePlaneStressElement(int n, Domain *aDomain)
Definition: igaelements.C:62
virtual int checkConsistency()
Performs consistency check.
Definition: igaelements.C:120
virtual void drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
Definition: igaelements.C:141
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
General purpose Plane stress structural element evaluator.
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
NURBSInterpolation interpolation
Definition: igaelements.h:100
Geometry wrapper for IGA elements.
Definition: iga.h:57
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: igaelements.C:65
the oofem namespace is to define a context or scope in which all oofem names are defined.
#define OOFEG_VARPLOT_PATTERN_LAYER
Class representing integration point in finite element program.
Definition: gausspoint.h:93
#define OOFEM_WARNING(...)
Definition: error.h:62
Class representing solution step.
Definition: timestep.h:80
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)=0
Evaluates global coordinates from given local ones.
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:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011