OOFEM 3.0
Loading...
Searching...
No Matches
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 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "sm/Elements/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
55namespace oofem {
56REGISTER_Element(BsplinePlaneStressElement);
57REGISTER_Element(NURBSPlaneStressElement);
58REGISTER_Element(TSplinePlaneStressElement);
59REGISTER_Element(NURBSSpace3dElement);
60
61
62BsplinePlaneStressElement :: BsplinePlaneStressElement(int n, Domain *aDomain) : IGAElement(n, aDomain), PlaneStressStructuralElementEvaluator(), interpolation(2) { }
63
64
65void BsplinePlaneStressElement :: initializeFrom(InputRecord &ir, int priority)
66{
67 //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
68 IGAElement :: initializeFrom(ir, priority);
69}
70
71
72int BsplinePlaneStressElement :: checkConsistency()
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
84NURBSPlaneStressElement :: NURBSPlaneStressElement(int n, Domain *aDomain) : IGAElement(n, aDomain), PlaneStressStructuralElementEvaluator(), interpolation(2) { }
85
86
87void NURBSPlaneStressElement :: initializeFrom(InputRecord &ir, int priority)
88{
89 //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
90 IGAElement :: initializeFrom(ir, priority);
91}
92
93
94int NURBSPlaneStressElement :: checkConsistency()
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
105TSplinePlaneStressElement :: TSplinePlaneStressElement(int n, Domain *aDomain) : IGATSplineElement(n, aDomain), PlaneStressStructuralElementEvaluator(), interpolation(2) { }
106
107
108
109
110NURBSSpace3dElement :: NURBSSpace3dElement(int n, Domain *aDomain) : IGAElement(n, aDomain), Space3dStructuralElementEvaluator(), interpolation(3) { }
111
112
113void NURBSSpace3dElement :: initializeFrom(InputRecord &ir, int priority)
114{
115 //PlaneStressStructuralElementEvaluator::initializeFrom(ir);
116 IGAElement :: initializeFrom(ir, priority);
117}
118
119
120int NURBSSpace3dElement :: checkConsistency()
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
141void BsplinePlaneStressElement :: drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
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 FloatArray *knotVector = interp->giveKnotVector();
157 const IntArray *span;
158 int j, nsd = this->giveNsd(this->giveGeometryType());
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
182 StructuralElementEvaluator :: computeVectorOf(VM_Total, tStep, u);
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
252void NURBSPlaneStressElement :: drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
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 FloatArray *knotVector = interp->giveKnotVector();
269 const IntArray *span;
270 int j, nsd = this->giveNsd(this->giveGeometryType());
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
294 StructuralElementEvaluator :: computeVectorOf(VM_Total, tStep, u);
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
454void TSplinePlaneStressElement :: drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
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 FloatArray *knotVector = interp->giveKnotVector();
471 const IntArray *span;
472 int j, nsd = this->giveNsd(this->giveGeometryType());
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
496 StructuralElementEvaluator :: computeVectorOf(VM_Total, tStep, u);
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
566void NURBSSpace3dElement :: drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
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 FloatArray *knotVector = interp->giveKnotVector();
584 const IntArray *span;
585 int j, nsd = this->giveNsd(this->giveGeometryType());
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
637 StructuralElementEvaluator :: computeVectorOf(VM_Total, tStep, u);
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
E
#define E(a, b)
REGISTER_Element
#define REGISTER_Element(class)
Definition classfactory.h:131
oofem::BSplineInterpolation::giveNumberOfControlPoints
virtual int giveNumberOfControlPoints(int dim) const
Definition feibspline.h:186
oofem::BsplinePlaneStressElement::giveNsd
int giveNsd(const Element_Geometry_Type) override
Definition igaelements.h:95
oofem::BsplinePlaneStressElement::giveGeometryType
Element_Geometry_Type giveGeometryType() const override
Definition igaelements.h:83
oofem::BsplinePlaneStressElement::interpolation
BSplineInterpolation interpolation
Definition igaelements.h:58
oofem::BsplinePlaneStressElement::giveInterpolation
FEInterpolation * giveInterpolation() const override
Definition igaelements.h:73
oofem::Element::integrationRulesArray
std::vector< std ::unique_ptr< IntegrationRule > > integrationRulesArray
Definition element.h:157
oofem::Element::giveNumberOfDofManagers
virtual int giveNumberOfDofManagers() const
Definition element.h:695
oofem::FEInterpolation::giveKnotVector
virtual const FloatArray * giveKnotVector() const
Definition feinterpol.h:530
oofem::FEInterpolation::local2global
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const =0
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatMatrix::beProductTOf
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:431
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::IGAElement::IGAElement
IGAElement(int n, Domain *aDomain)
Definition iga.h:98
oofem::IGATSplineElement::IGATSplineElement
IGATSplineElement(int n, Domain *aDomain)
Definition iga.h:121
oofem::IntArray::resize
void resize(int n)
Definition intarray.C:73
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::NURBSPlaneStressElement::giveInterpolation
FEInterpolation * giveInterpolation() const override
Definition igaelements.h:117
oofem::NURBSPlaneStressElement::interpolation
NURBSInterpolation interpolation
Definition igaelements.h:102
oofem::NURBSPlaneStressElement::giveNsd
int giveNsd(const Element_Geometry_Type) override
Definition igaelements.h:140
oofem::NURBSPlaneStressElement::giveGeometryType
Element_Geometry_Type giveGeometryType() const override
Definition igaelements.h:127
oofem::NURBSSpace3dElement::giveNsd
int giveNsd(const Element_Geometry_Type) override
Definition igaelements.h:227
oofem::NURBSSpace3dElement::giveInterpolation
FEInterpolation * giveInterpolation() const override
Definition igaelements.h:205
oofem::NURBSSpace3dElement::interpolation
NURBSInterpolation interpolation
Definition igaelements.h:190
oofem::NURBSSpace3dElement::giveGeometryType
Element_Geometry_Type giveGeometryType() const override
Definition igaelements.h:216
oofem::PlaneStressStructuralElementEvaluator::computeStressVector
void computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep) override
Definition planestresselementevaluator.C:95
oofem::Space3dStructuralElementEvaluator::computeStressVector
void computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep) override
Definition space3delementevaluator.C:100
oofem::StructuralElementEvaluator::computeStrainVector
void computeStrainVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, FloatArray &u)
Definition structuralelementevaluator.C:312
oofem::TSplinePlaneStressElement::giveInterpolation
FEInterpolation * giveInterpolation() const override
Definition igaelements.h:165
oofem::TSplinePlaneStressElement::giveGeometryType
Element_Geometry_Type giveGeometryType() const override
Definition igaelements.h:175
oofem::TSplinePlaneStressElement::giveNsd
int giveNsd(const Element_Geometry_Type) override
Definition igaelements.h:183
oofem::TSplinePlaneStressElement::interpolation
TSplineInterpolation interpolation
Definition igaelements.h:148
OOFEM_WARNING
#define OOFEM_WARNING(...)
Definition error.h:80
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
M_PI
#define M_PI
Definition mathfem.h:52
oofem::product
double product(const FloatArray &x)
Definition floatarray.C:1051
gc
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
OOFEG_VARPLOT_PATTERN_LAYER
#define OOFEG_VARPLOT_PATTERN_LAYER
Definition oofeggraphiccontext.h:105
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011