OOFEM 3.0
Loading...
Searching...
No Matches
iga.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 "inputrecord.h"
36#include "floatarray.h"
37#include "floatmatrix.h"
38#include "mathfem.h"
39#include "iga.h"
40#include "gausspoint.h"
41#include "feitspline.h"
42#include "parametermanager.h"
43#include "paramkey.h"
44
45
46#ifdef __OOFEG
47 #include "oofeggraphiccontext.h"
48 #include "oofegutils.h"
49 #include "sm/Elements/structuralelementevaluator.h"
50#endif
51
52namespace oofem {
53
54#ifdef __MPI_PARALLEL_MODE
55ParamKey IGAElement::IPK_IGAElement_knotSpanParallelMode(_IFT_IGAElement_KnotSpanParallelMode);
56#endif
57
58
59void IGAElement :: initializeFrom(InputRecord &ir, int priority)
60{
61 ParameterManager &ppm = domain->elementPPM;
62 Element :: initializeFrom(ir, priority); // read nodes , material, cross section
63 this->giveInterpolation()->initializeFrom(ir, ppm, this->number, priority); // read geometry
64#ifdef __MPI_PARALLEL_MODE
65 PM_UPDATE_PARAMETER(knotSpanParallelMode, ppm, ir, this->number, IPK_IGAElement_knotSpanParallelMode, priority);
66#endif
67}
68
69void
70IGAElement::initializeFinish ()
71{
72 int indx = 0;
73 ParameterManager &ppm = domain->elementPPM;
74
75#ifdef __MPI_PARALLEL_MODE
76 int numberOfKnotSpans = 0;
77#endif
78 Element::initializeFinish();
79 // set number of dofmanagers
80 this->numberOfDofMans = dofManArray.giveSize();
81 this->giveInterpolation()->postInitialize(ppm, this->number);
82
83
84 // generate individual IntegrationElements; one for each nonzero knot span
85 int nsd = this->giveNsd(this->giveGeometryType());
86 if ( nsd == 1 ) {
87 //HUHU
88 } else if ( nsd == 2 ) {
89 int numberOfKnotSpansU = this->giveInterpolation()->giveNumberOfKnotSpans(1);
90 int numberOfKnotSpansV = this->giveInterpolation()->giveNumberOfKnotSpans(2);
91#ifdef __MPI_PARALLEL_MODE
92 numberOfKnotSpans = numberOfKnotSpansU * numberOfKnotSpansV;
93#endif
94 const IntArray *knotMultiplicityU = this->giveInterpolation()->giveKnotMultiplicity(1);
95 const IntArray *knotMultiplicityV = this->giveInterpolation()->giveKnotMultiplicity(2);
96 const FloatArray *knotValuesU = this->giveInterpolation()->giveKnotValues(1);
97 const FloatArray *knotValuesV = this->giveInterpolation()->giveKnotValues(2);
98
99 FloatArray newgpcoords(2);
100 IntArray knotSpan(2);
101
102 int numberOfIntegrationRules = numberOfKnotSpansU * numberOfKnotSpansV;
103 integrationRulesArray.resize( numberOfIntegrationRules );
104
105 knotSpan.at(2) = -1;
106 for ( int vi = 1; vi <= numberOfKnotSpansV; vi++ ) {
107 double dv = knotValuesV->at(vi + 1) - knotValuesV->at(vi);
108 knotSpan.at(2) += knotMultiplicityV->at(vi);
109
110 knotSpan.at(1) = -1;
111 for ( int ui = 1; ui <= numberOfKnotSpansU; ui++ ) {
112 double du = knotValuesU->at(ui + 1) - knotValuesU->at(ui);
113 knotSpan.at(1) += knotMultiplicityU->at(ui);
114
115 integrationRulesArray [ indx ] = std::make_unique<IGAIntegrationElement>(indx, this, knotSpan);
116 integrationRulesArray [ indx ]->SetUpPointsOnSquare(numberOfGaussPoints, _PlaneStress); // HUHU _PlaneStress, rectangle
117
118 // remap local subelement gp coordinates into knot span coordinates and update integration weight
119 for ( auto &gp: *integrationRulesArray [ indx ] ) {
120 const FloatArray &gpcoords = gp->giveNaturalCoordinates();
121
122 newgpcoords.at(1) = knotValuesU->at(ui) + du * ( gpcoords.at(1) / 2.0 + 0.5 );
123 newgpcoords.at(2) = knotValuesV->at(vi) + dv * ( gpcoords.at(2) / 2.0 + 0.5 );
124 gp->setNaturalCoordinates(newgpcoords);
125 gp->setWeight(gp->giveWeight() / 4.0 * du * dv);
126 }
127
128 indx++;
129 }
130 }
131 } else if ( nsd == 3 ) {
132 int numberOfKnotSpansU = this->giveInterpolation()->giveNumberOfKnotSpans(1);
133 int numberOfKnotSpansV = this->giveInterpolation()->giveNumberOfKnotSpans(2);
134 int numberOfKnotSpansW = this->giveInterpolation()->giveNumberOfKnotSpans(3);
135#ifdef __MPI_PARALLEL_MODE
136 numberOfKnotSpans = numberOfKnotSpansU * numberOfKnotSpansV * numberOfKnotSpansW;
137#endif
138 const IntArray *knotMultiplicityU = this->giveInterpolation()->giveKnotMultiplicity(1);
139 const IntArray *knotMultiplicityV = this->giveInterpolation()->giveKnotMultiplicity(2);
140 const IntArray *knotMultiplicityW = this->giveInterpolation()->giveKnotMultiplicity(3);
141 const FloatArray *knotValuesU = this->giveInterpolation()->giveKnotValues(1);
142 const FloatArray *knotValuesV = this->giveInterpolation()->giveKnotValues(2);
143 const FloatArray *knotValuesW = this->giveInterpolation()->giveKnotValues(3);
144
145 FloatArray newgpcoords(3);
146 IntArray knotSpan(3);
147
148 int numberOfIntegrationRules = numberOfKnotSpansU * numberOfKnotSpansV * numberOfKnotSpansW;
149 integrationRulesArray.resize( numberOfIntegrationRules );
150
151 knotSpan.at(3) = -1;
152 for ( int wi = 1; wi <= numberOfKnotSpansW; wi++ ) {
153 double dw = knotValuesW->at(wi + 1) - knotValuesW->at(wi);
154 knotSpan.at(3) += knotMultiplicityW->at(wi);
155
156 knotSpan.at(2) = -1;
157 for ( int vi = 1; vi <= numberOfKnotSpansV; vi++ ) {
158 double dv = knotValuesV->at(vi + 1) - knotValuesV->at(vi);
159 knotSpan.at(2) += knotMultiplicityV->at(vi);
160
161 knotSpan.at(1) = -1;
162 for ( int ui = 1; ui <= numberOfKnotSpansU; ui++ ) {
163 double du = knotValuesU->at(ui + 1) - knotValuesU->at(ui);
164 knotSpan.at(1) += knotMultiplicityU->at(ui);
165
166 integrationRulesArray [ indx ] = std::make_unique<IGAIntegrationElement>(indx, this, knotSpan);
167 integrationRulesArray [ indx ]->SetUpPointsOnCube(numberOfGaussPoints, _3dMat);
168
169 // remap local subelement gp coordinates into knot span coordinates and update integration weight
170 for ( auto &gp: *integrationRulesArray [ indx ] ) {
171 const FloatArray &gpcoords = gp->giveNaturalCoordinates();
172
173 newgpcoords.at(1) = knotValuesU->at(ui) + du * ( gpcoords.at(1) / 2.0 + 0.5 );
174 newgpcoords.at(2) = knotValuesV->at(vi) + dv * ( gpcoords.at(2) / 2.0 + 0.5 );
175 newgpcoords.at(3) = knotValuesW->at(wi) + dw * ( gpcoords.at(3) / 2.0 + 0.5 );
176 gp->setNaturalCoordinates(newgpcoords);
177 gp->setWeight(gp->giveWeight() / 8.0 * du * dv * dw);
178 }
179
180 indx++;
181 }
182 }
183 }
184 } else {
185 throw ComponentInputException(ComponentInputException::ComponentType::ctElement, this->number, "Unsupported number of spatial dimensions");
186 }
187
188#ifdef __MPI_PARALLEL_MODE
189 if (ppm.checkIfSet(this->number, IPK_IGAElement_knotSpanParallelMode.getIndex()) == false) {
190 this->knotSpanParallelMode.resize(numberOfKnotSpans);
191 // set Element_local as default
192 for ( int i = 1; i <= numberOfKnotSpans; i++ ) {
193 knotSpanParallelMode.at(i) = Element_local;
194 }
195 }
196#endif
197}
198
199#ifdef __MPI_PARALLEL_MODE
200elementParallelMode
201IGAElement :: giveKnotSpanParallelMode(int knotSpanIndex) const
202{
203 elementParallelMode emode = this->giveParallelMode();
204 if ( emode == Element_remote ) {
205 return Element_remote;
206 } else if ( emode == Element_local ) {
207 return ( elementParallelMode ) this->knotSpanParallelMode.at(knotSpanIndex + 1);
208 } else {
209 OOFEM_ERROR("Cannot determine elementParallelMode");
210 }
211
212 return Element_local; //to make compiler happy
213}
214
215#endif // __MPI_PARALLEL_MODE
216
217
218// integration elements are setup in the same way as for IGAElement for now HUHU
219
220void IGATSplineElement :: initializeFrom(InputRecord &ir, int priority)
221{
222 ParameterManager &pm = domain->elementPPM;
223 TSplineInterpolation *interpol = static_cast< TSplineInterpolation * >( this->giveInterpolation() );
224 Element :: initializeFrom(ir, priority); // read nodes , material, cross section
225
226 interpol->initializeFrom(ir, pm, this->number, priority); // read geometry
227}
228
229void
230IGATSplineElement::initializeFinish()
231{
232 ParameterManager &pm = domain->elementPPM;
233 IGAElement::initializeFinish();
234 int indx = 0;
235 TSplineInterpolation *interpol = static_cast< TSplineInterpolation * >( this->giveInterpolation() );
236 // set number of dofmanagers
237 this->numberOfDofMans = dofManArray.giveSize();
238 // set number of control points before initialization HUHU HAHA
239 interpol->setNumberOfControlPoints(this->numberOfDofMans);
240 interpol->postInitialize(pm, this->number);
241
242 // generate individual IntegrationElements; one for each nonzero knot span
243 int nsd = giveNsd(this->giveGeometryType());
244 if ( nsd == 2 ) {
245 int numberOfKnotSpansU = this->giveInterpolation()->giveNumberOfKnotSpans(1);
246 int numberOfKnotSpansV = this->giveInterpolation()->giveNumberOfKnotSpans(2);
247 const IntArray *knotMultiplicityU = this->giveInterpolation()->giveKnotMultiplicity(1);
248 const IntArray *knotMultiplicityV = this->giveInterpolation()->giveKnotMultiplicity(2);
249 const FloatArray *knotValuesU = this->giveInterpolation()->giveKnotValues(1);
250 const FloatArray *knotValuesV = this->giveInterpolation()->giveKnotValues(2);
251
252 FloatArray newgpcoords(2);
253 IntArray knotSpan(2);
254
255 int numberOfIntegrationRules = numberOfKnotSpansU * numberOfKnotSpansV;
256 integrationRulesArray.resize( numberOfIntegrationRules );
257
258 knotSpan.at(2) = -1;
259 for ( int vi = 1; vi <= numberOfKnotSpansV; vi++ ) {
260 double dv = knotValuesV->at(vi + 1) - knotValuesV->at(vi);
261 knotSpan.at(2) += knotMultiplicityV->at(vi);
262
263 knotSpan.at(1) = -1;
264 for ( int ui = 1; ui <= numberOfKnotSpansU; ui++ ) {
265 double du = knotValuesU->at(ui + 1) - knotValuesU->at(ui);
266 knotSpan.at(1) += knotMultiplicityU->at(ui);
267
268 integrationRulesArray [ indx ] = std::make_unique<IGAIntegrationElement>(indx, this, knotSpan);
269 integrationRulesArray [ indx ]->SetUpPointsOnSquare(numberOfGaussPoints, _PlaneStress); // HUHU _PlaneStress, rectangle
270
271 // remap local subelement gp coordinates into knot span coordinates and update integration weight
272 for ( auto &gp: *integrationRulesArray [ indx ] ) {
273 const auto &gpcoords = gp->giveNaturalCoordinates();
274
275 newgpcoords.at(1) = knotValuesU->at(ui) + du * ( gpcoords.at(1) / 2.0 + 0.5 );
276 newgpcoords.at(2) = knotValuesV->at(vi) + dv * ( gpcoords.at(2) / 2.0 + 0.5 );
277 gp->setNaturalCoordinates(newgpcoords);
278 gp->setWeight(gp->giveWeight() / 4.0 * du * dv);
279 }
280
281 indx++;
282 }
283 }
284 } else {
285 throw ComponentInputException(ComponentInputException::ComponentType::ctElement, this->number, "Domain: unsupported number of spatial dimensions");
286 }
287}
288
289#ifdef __OOFEG
290
291 #define DRAW_MESH
292
293// if DRAW_MESH is defined only boundary of integration elements are drawn;
294// currently mesh is not properly drawn for tsplines
295// because integration elements (does not matter whether single span or multi span)
296// are generaly finer than T-mesh;
297
298void IGAElement :: drawRawGeometry(oofegGraphicContext &gc, TimeStep *tStep)
299{
300 WCRec p [ 8 ];
301 FEInterpolation *interp = this->giveInterpolation();
302 int nseg;
303
304 #ifdef DRAW_MESH
305 WCRec pp [ 2 ];
306 #endif
307
308 if ( !gc.testElementGraphicActivity(this) ) {
309 return;
310 }
311
312 EASValsSetLineWidth(OOFEG_RAW_GEOMETRY_WIDTH);
313 EASValsSetColor( gc.getElementColor() );
314 EASValsSetEdgeColor( gc.getElementEdgeColor() );
315 EASValsSetEdgeFlag(true);
316 EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
317 EASValsSetLineStyle(SOLID_STYLE);
318 EASValsSetFillStyle(FILL_SOLID);
319 //EASValsSetLineWidth(0);
320
321 #ifdef DRAW_MESH
322 nseg = 8;
323 #else
324 nseq = 4;
325 #endif
326
327 const FloatArray *knotVector = interp->giveKnotVector();
328 int nsd = this->giveNsd(this->giveGeometryType());
329
330 if ( nsd == 1 ) {
331 FloatArray c [ 2 ], cg [ 2 ];
332
333 for ( int j = 0; j < 2; j++ ) {
334 c [ j ].resize(1);
335 cg [ j ].resize(1);
336 }
337
338 // loop over individual integration rules (i.e., knot spans)
339 for ( auto &iRule: integrationRulesArray ) {
340 const IntArray *span = iRule->giveKnotSpan();
341 // divide span locally to get finer geometry rep.
342 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
343 for ( int i = 1; i <= nseg; i++ ) {
344 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
345 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
346
347 for ( int k = 0; k < 2; k++ ) {
348 interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
349 p [ k ].x = ( FPNum ) cg [ k ].at(1);
350 p [ k ].y = 0.;
351 p [ k ].z = 0.;
352 }
353
354 auto go = CreateLine3D(p);
355 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
356 EGAttachObject(go, ( EObjectP ) this);
357 EMAddGraphicsToModel(ESIModel(), go);
358 }
359 } // end loop over knot spans (irules)
360 } else if ( nsd == 2 ) {
361 FloatArray c [ 4 ], cg [ 4 ];
362
363 for ( int j = 0; j < 4; j++ ) {
364 c [ j ].resize(2);
365 cg [ j ].resize(2);
366 }
367
368 // loop over individual integration rules (i.e., knot spans)
369 for ( auto &iRule: integrationRulesArray ) {
370 const IntArray *span = iRule->giveKnotSpan();
371 // divide span locally to get finer geometry rep.
372 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
373 double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
374 for ( int i = 1; i <= nseg; i++ ) {
375 for ( int j = 1; j <= nseg; j++ ) {
376 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
377 c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
378 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
379 c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
380 c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
381 c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
382 c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
383 c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
384
385 for ( int k = 0; k < 4; k++ ) {
386 interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
387 p [ k ].x = ( FPNum ) cg [ k ].at(1);
388 p [ k ].y = ( FPNum ) cg [ k ].at(2);
389 p [ k ].z = 0.;
390 }
391
392 #ifdef DRAW_MESH
393 if ( i == 1 ) {
394 pp [ 0 ].x = p [ 0 ].x;
395 pp [ 0 ].y = p [ 0 ].y;
396 pp [ 0 ].z = p [ 0 ].z;
397
398 pp [ 1 ].x = p [ 3 ].x;
399 pp [ 1 ].y = p [ 3 ].y;
400 pp [ 1 ].z = p [ 3 ].z;
401
402 auto go = CreateLine3D(pp);
403 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
404 EGAttachObject(go, ( EObjectP ) this);
405 EMAddGraphicsToModel(ESIModel(), go);
406 }
407
408 if ( j == 1 ) {
409 pp [ 0 ].x = p [ 0 ].x;
410 pp [ 0 ].y = p [ 0 ].y;
411 pp [ 0 ].z = p [ 0 ].z;
412
413 pp [ 1 ].x = p [ 1 ].x;
414 pp [ 1 ].y = p [ 1 ].y;
415 pp [ 1 ].z = p [ 1 ].z;
416
417 auto go = CreateLine3D(pp);
418 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
419 EGAttachObject(go, ( EObjectP ) this);
420 EMAddGraphicsToModel(ESIModel(), go);
421 }
422
423 if ( i == nseg ) {
424 pp [ 0 ].x = p [ 1 ].x;
425 pp [ 0 ].y = p [ 1 ].y;
426 pp [ 0 ].z = p [ 1 ].z;
427
428 pp [ 1 ].x = p [ 2 ].x;
429 pp [ 1 ].y = p [ 2 ].y;
430 pp [ 1 ].z = p [ 2 ].z;
431
432 auto go = CreateLine3D(pp);
433 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
434 EGAttachObject(go, ( EObjectP ) this);
435 EMAddGraphicsToModel(ESIModel(), go);
436 }
437
438 if ( j == nseg ) {
439 pp [ 0 ].x = p [ 2 ].x;
440 pp [ 0 ].y = p [ 2 ].y;
441 pp [ 0 ].z = p [ 2 ].z;
442
443 pp [ 1 ].x = p [ 3 ].x;
444 pp [ 1 ].y = p [ 3 ].y;
445 pp [ 1 ].z = p [ 3 ].z;
446
447 auto go = CreateLine3D(pp);
448 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
449 EGAttachObject(go, ( EObjectP ) this);
450 EMAddGraphicsToModel(ESIModel(), go);
451 }
452
453 #else
454 auto go = CreateQuad3D(p);
455 EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
456 EGAttachObject(go, ( EObjectP ) this);
457 EMAddGraphicsToModel(ESIModel(), go);
458 #endif
459 }
460 }
461 } // end loop over knot spans (irules)
462 } else if ( nsd == 3 ) {
463 FloatArray c [ 8 ], cg [ 8 ];
464
465 for ( int j = 0; j < 8; j++ ) {
466 c [ j ].resize(3);
467 cg [ j ].resize(3);
468 }
469
470 // loop over individual integration rules (i.e., knot spans)
471 for ( auto &iRule: integrationRulesArray ) {
472 const IntArray *span = iRule->giveKnotSpan();
473 // divide span locally to get finer geometry rep.
474 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
475 double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
476 double dt = ( knotVector [ 2 ] [ span->at(3) + 1 ] - knotVector [ 2 ] [ span->at(3) ] ) / nseg;
477 for ( int i = 1; i <= nseg; i++ ) {
478 for ( int j = 1; j <= nseg; j++ ) {
479 for ( int k = 1; k <= nseg; k++ ) {
480 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
481 c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
482 c [ 0 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
483 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
484 c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
485 c [ 1 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
486 c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
487 c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
488 c [ 2 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
489 c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
490 c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
491 c [ 3 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
492 c [ 4 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
493 c [ 4 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
494 c [ 4 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
495 c [ 5 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
496 c [ 5 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
497 c [ 5 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
498 c [ 6 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
499 c [ 6 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
500 c [ 6 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
501 c [ 7 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
502 c [ 7 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
503 c [ 7 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
504
505 for ( int m = 0; m < 8; m++ ) {
506 interp->local2global( cg [ m ], c [ m ], FEIIGAElementGeometryWrapper( this, iRule->giveKnotSpan() ) );
507 p [ m ].x = ( FPNum ) cg [ m ].at(1);
508 p [ m ].y = ( FPNum ) cg [ m ].at(2);
509 p [ m ].z = ( FPNum ) cg [ m ].at(3);
510 }
511
512 #ifdef DRAW_MESH
513 if ( i == 1 && j == 1 ) {
514 pp [ 0 ].x = p [ 0 ].x;
515 pp [ 0 ].y = p [ 0 ].y;
516 pp [ 0 ].z = p [ 0 ].z;
517
518 pp [ 1 ].x = p [ 4 ].x;
519 pp [ 1 ].y = p [ 4 ].y;
520 pp [ 1 ].z = p [ 4 ].z;
521
522 #ifdef DRAW_VISIBLE_CONTOUR
523 #ifdef SPHERE_WITH_HOLE
524 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
525 for ( int ii = 0; ii < 2; ii++ ) {
526 xx += pp [ ii ].x;
527 yy += pp [ ii ].y;
528 zz += pp [ ii ].z;
529 }
530 xx /= 2.0;
531 yy /= 2.0;
532 zz /= 2.0;
533 rr = xx * xx + yy * yy;
534 r = rr + zz * zz;
535
536 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
537 if ( zz < 2.0001 || rr < 1.001 * 1.001 || yy < 0.0001 ) {
538 auto go = CreateLine3D(pp);
539 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
540 EGAttachObject(go, ( EObjectP ) this);
541 EMAddGraphicsToModel(ESIModel(), go);
542 }
543 }
544 #endif
545 #else
546 auto go = CreateLine3D(pp);
547 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
548 EGAttachObject(go, ( EObjectP ) this);
549 EMAddGraphicsToModel(ESIModel(), go);
550 #endif
551 }
552
553 if ( i == 1 && j == nseg ) {
554 pp [ 0 ].x = p [ 3 ].x;
555 pp [ 0 ].y = p [ 3 ].y;
556 pp [ 0 ].z = p [ 3 ].z;
557
558 pp [ 1 ].x = p [ 7 ].x;
559 pp [ 1 ].y = p [ 7 ].y;
560 pp [ 1 ].z = p [ 7 ].z;
561
562 #ifdef DRAW_VISIBLE_CONTOUR
563 #ifdef SPHERE_WITH_HOLE
564 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
565 for ( int ii = 0; ii < 2; ii++ ) {
566 xx += pp [ ii ].x;
567 yy += pp [ ii ].y;
568 zz += pp [ ii ].z;
569 }
570 xx /= 2.0;
571 yy /= 2.0;
572 zz /= 2.0;
573 rr = xx * xx + yy * yy;
574 r = rr + zz * zz;
575
576 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
577 if ( zz < 2.0001 || rr < 1.001 * 1.001 || yy < 0.0001 ) {
578 auto go = CreateLine3D(pp);
579 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
580 EGAttachObject(go, ( EObjectP ) this);
581 EMAddGraphicsToModel(ESIModel(), go);
582 }
583 }
584 #endif
585 #else
586 auto go = CreateLine3D(pp);
587 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
588 EGAttachObject(go, ( EObjectP ) this);
589 EMAddGraphicsToModel(ESIModel(), go);
590 #endif
591 }
592
593 if ( i == nseg && j == 1 ) {
594 pp [ 0 ].x = p [ 1 ].x;
595 pp [ 0 ].y = p [ 1 ].y;
596 pp [ 0 ].z = p [ 1 ].z;
597
598 pp [ 1 ].x = p [ 5 ].x;
599 pp [ 1 ].y = p [ 5 ].y;
600 pp [ 1 ].z = p [ 5 ].z;
601
602 #ifdef DRAW_VISIBLE_CONTOUR
603 #ifdef SPHERE_WITH_HOLE
604 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
605 for ( int ii = 0; ii < 2; ii++ ) {
606 xx += pp [ ii ].x;
607 yy += pp [ ii ].y;
608 zz += pp [ ii ].z;
609 }
610 xx /= 2.0;
611 yy /= 2.0;
612 zz /= 2.0;
613 rr = xx * xx + yy * yy;
614 r = rr + zz * zz;
615
616 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
617 if ( zz < 2.0001 || rr < 1.001 * 1.001 || yy < 0.0001 ) {
618 auto go = CreateLine3D(pp);
619 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
620 EGAttachObject(go, ( EObjectP ) this);
621 EMAddGraphicsToModel(ESIModel(), go);
622 }
623 }
624 #endif
625 #else
626 auto go = CreateLine3D(pp);
627 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
628 EGAttachObject(go, ( EObjectP ) this);
629 EMAddGraphicsToModel(ESIModel(), go);
630 #endif
631 }
632
633 if ( i == nseg && j == nseg ) {
634 pp [ 0 ].x = p [ 2 ].x;
635 pp [ 0 ].y = p [ 2 ].y;
636 pp [ 0 ].z = p [ 2 ].z;
637
638 pp [ 1 ].x = p [ 6 ].x;
639 pp [ 1 ].y = p [ 6 ].y;
640 pp [ 1 ].z = p [ 6 ].z;
641
642 #ifdef DRAW_VISIBLE_CONTOUR
643 #ifdef SPHERE_WITH_HOLE
644 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
645 for ( int ii = 0; ii < 2; ii++ ) {
646 xx += pp [ ii ].x;
647 yy += pp [ ii ].y;
648 zz += pp [ ii ].z;
649 }
650 xx /= 2.0;
651 yy /= 2.0;
652 zz /= 2.0;
653 rr = xx * xx + yy * yy;
654 r = rr + zz * zz;
655
656 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
657 if ( zz < 2.0001 || rr < 1.001 * 1.001 || yy < 0.0001 ) {
658 auto go = CreateLine3D(pp);
659 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
660 EGAttachObject(go, ( EObjectP ) this);
661 EMAddGraphicsToModel(ESIModel(), go);
662 }
663 }
664 #endif
665 #else
666 auto go = CreateLine3D(pp);
667 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
668 EGAttachObject(go, ( EObjectP ) this);
669 EMAddGraphicsToModel(ESIModel(), go);
670 #endif
671 }
672
673 if ( j == 1 && k == 1 ) {
674 pp [ 0 ].x = p [ 0 ].x;
675 pp [ 0 ].y = p [ 0 ].y;
676 pp [ 0 ].z = p [ 0 ].z;
677
678 pp [ 1 ].x = p [ 1 ].x;
679 pp [ 1 ].y = p [ 1 ].y;
680 pp [ 1 ].z = p [ 1 ].z;
681
682 #ifdef DRAW_VISIBLE_CONTOUR
683 #ifdef SPHERE_WITH_HOLE
684 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
685 for ( int ii = 0; ii < 2; ii++ ) {
686 xx += pp [ ii ].x;
687 yy += pp [ ii ].y;
688 zz += pp [ ii ].z;
689 }
690 xx /= 2.0;
691 yy /= 2.0;
692 zz /= 2.0;
693 rr = xx * xx + yy * yy;
694 r = rr + zz * zz;
695
696 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
697 auto go = CreateLine3D(pp);
698 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
699 EGAttachObject(go, ( EObjectP ) this);
700 EMAddGraphicsToModel(ESIModel(), go);
701 }
702 #endif
703 #else
704 auto go = CreateLine3D(pp);
705 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
706 EGAttachObject(go, ( EObjectP ) this);
707 EMAddGraphicsToModel(ESIModel(), go);
708 #endif
709 }
710
711 if ( j == 1 && k == nseg ) {
712 pp [ 0 ].x = p [ 4 ].x;
713 pp [ 0 ].y = p [ 4 ].y;
714 pp [ 0 ].z = p [ 4 ].z;
715
716 pp [ 1 ].x = p [ 5 ].x;
717 pp [ 1 ].y = p [ 5 ].y;
718 pp [ 1 ].z = p [ 5 ].z;
719
720 #ifdef DRAW_VISIBLE_CONTOUR
721 #ifdef SPHERE_WITH_HOLE
722 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
723 for ( int ii = 0; ii < 2; ii++ ) {
724 xx += pp [ ii ].x;
725 yy += pp [ ii ].y;
726 zz += pp [ ii ].z;
727 }
728 xx /= 2.0;
729 yy /= 2.0;
730 zz /= 2.0;
731 rr = xx * xx + yy * yy;
732 r = rr + zz * zz;
733
734 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
735 if ( yy < 1.5 || zz < 2.0001 ) {
736 auto go = CreateLine3D(pp);
737 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
738 EGAttachObject(go, ( EObjectP ) this);
739 EMAddGraphicsToModel(ESIModel(), go);
740 }
741 }
742 #endif
743 #else
744 auto go = CreateLine3D(pp);
745 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
746 EGAttachObject(go, ( EObjectP ) this);
747 EMAddGraphicsToModel(ESIModel(), go);
748 #endif
749 }
750
751 if ( j == nseg && k == 1 ) {
752 pp [ 0 ].x = p [ 3 ].x;
753 pp [ 0 ].y = p [ 3 ].y;
754 pp [ 0 ].z = p [ 3 ].z;
755
756 pp [ 1 ].x = p [ 2 ].x;
757 pp [ 1 ].y = p [ 2 ].y;
758 pp [ 1 ].z = p [ 2 ].z;
759
760 #ifdef DRAW_VISIBLE_CONTOUR
761 #ifdef SPHERE_WITH_HOLE
762 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
763 for ( int ii = 0; ii < 2; ii++ ) {
764 xx += pp [ ii ].x;
765 yy += pp [ ii ].y;
766 zz += pp [ ii ].z;
767 }
768 xx /= 2.0;
769 yy /= 2.0;
770 zz /= 2.0;
771 rr = xx * xx + yy * yy;
772 r = rr + zz * zz;
773
774 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
775 auto go = CreateLine3D(pp);
776 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
777 EGAttachObject(go, ( EObjectP ) this);
778 EMAddGraphicsToModel(ESIModel(), go);
779 }
780 #endif
781 #else
782 auto go = CreateLine3D(pp);
783 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
784 EGAttachObject(go, ( EObjectP ) this);
785 EMAddGraphicsToModel(ESIModel(), go);
786 #endif
787 }
788
789 if ( j == nseg && k == nseg ) {
790 pp [ 0 ].x = p [ 7 ].x;
791 pp [ 0 ].y = p [ 7 ].y;
792 pp [ 0 ].z = p [ 7 ].z;
793
794 pp [ 1 ].x = p [ 6 ].x;
795 pp [ 1 ].y = p [ 6 ].y;
796 pp [ 1 ].z = p [ 6 ].z;
797
798 #ifdef DRAW_VISIBLE_CONTOUR
799 #ifdef SPHERE_WITH_HOLE
800 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
801 for ( int ii = 0; ii < 2; ii++ ) {
802 xx += pp [ ii ].x;
803 yy += pp [ ii ].y;
804 zz += pp [ ii ].z;
805 }
806 xx /= 2.0;
807 yy /= 2.0;
808 zz /= 2.0;
809 rr = xx * xx + yy * yy;
810 r = rr + zz * zz;
811
812 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
813 if ( yy < 1.5 || zz < 2.0001 ) {
814 auto go = CreateLine3D(pp);
815 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
816 EGAttachObject(go, ( EObjectP ) this);
817 EMAddGraphicsToModel(ESIModel(), go);
818 }
819 }
820 #endif
821 #else
822 auto go = CreateLine3D(pp);
823 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
824 EGAttachObject(go, ( EObjectP ) this);
825 EMAddGraphicsToModel(ESIModel(), go);
826 #endif
827 }
828
829 if ( k == 1 && i == 1 ) {
830 pp [ 0 ].x = p [ 0 ].x;
831 pp [ 0 ].y = p [ 0 ].y;
832 pp [ 0 ].z = p [ 0 ].z;
833
834 pp [ 1 ].x = p [ 3 ].x;
835 pp [ 1 ].y = p [ 3 ].y;
836 pp [ 1 ].z = p [ 3 ].z;
837
838 #ifdef DRAW_VISIBLE_CONTOUR
839 #ifdef SPHERE_WITH_HOLE
840 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
841 for ( int ii = 0; ii < 2; ii++ ) {
842 xx += pp [ ii ].x;
843 yy += pp [ ii ].y;
844 zz += pp [ ii ].z;
845 }
846 xx /= 2.0;
847 yy /= 2.0;
848 zz /= 2.0;
849 rr = xx * xx + yy * yy;
850 r = rr + zz * zz;
851
852 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
853 auto go = CreateLine3D(pp);
854 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
855 EGAttachObject(go, ( EObjectP ) this);
856 EMAddGraphicsToModel(ESIModel(), go);
857 }
858 #endif
859 #else
860 auto go = CreateLine3D(pp);
861 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
862 EGAttachObject(go, ( EObjectP ) this);
863 EMAddGraphicsToModel(ESIModel(), go);
864 #endif
865 }
866
867 if ( k == 1 && i == nseg ) {
868 pp [ 0 ].x = p [ 1 ].x;
869 pp [ 0 ].y = p [ 1 ].y;
870 pp [ 0 ].z = p [ 1 ].z;
871
872 pp [ 1 ].x = p [ 2 ].x;
873 pp [ 1 ].y = p [ 2 ].y;
874 pp [ 1 ].z = p [ 2 ].z;
875
876 #ifdef DRAW_VISIBLE_CONTOUR
877 #ifdef SPHERE_WITH_HOLE
878 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
879 for ( int ii = 0; ii < 2; ii++ ) {
880 xx += pp [ ii ].x;
881 yy += pp [ ii ].y;
882 zz += pp [ ii ].z;
883 }
884 xx /= 2.0;
885 yy /= 2.0;
886 zz /= 2.0;
887 rr = xx * xx + yy * yy;
888 r = rr + zz * zz;
889
890 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
891 auto go = CreateLine3D(pp);
892 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
893 EGAttachObject(go, ( EObjectP ) this);
894 EMAddGraphicsToModel(ESIModel(), go);
895 }
896 #endif
897 #else
898 auto go = CreateLine3D(pp);
899 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
900 EGAttachObject(go, ( EObjectP ) this);
901 EMAddGraphicsToModel(ESIModel(), go);
902 #endif
903 }
904
905 if ( k == nseg && i == 1 ) {
906 pp [ 0 ].x = p [ 4 ].x;
907 pp [ 0 ].y = p [ 4 ].y;
908 pp [ 0 ].z = p [ 4 ].z;
909
910 pp [ 1 ].x = p [ 7 ].x;
911 pp [ 1 ].y = p [ 7 ].y;
912 pp [ 1 ].z = p [ 7 ].z;
913
914 #ifdef DRAW_VISIBLE_CONTOUR
915 #ifdef SPHERE_WITH_HOLE
916 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
917 for ( int ii = 0; ii < 2; ii++ ) {
918 xx += pp [ ii ].x;
919 yy += pp [ ii ].y;
920 zz += pp [ ii ].z;
921 }
922 xx /= 2.0;
923 yy /= 2.0;
924 zz /= 2.0;
925 rr = xx * xx + yy * yy;
926 r = rr + zz * zz;
927
928 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
929 if ( yy < 0.0001 ) {
930 auto go = CreateLine3D(pp);
931 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
932 EGAttachObject(go, ( EObjectP ) this);
933 EMAddGraphicsToModel(ESIModel(), go);
934 }
935 }
936 #endif
937 #else
938 auto go = CreateLine3D(pp);
939 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
940 EGAttachObject(go, ( EObjectP ) this);
941 EMAddGraphicsToModel(ESIModel(), go);
942 #endif
943 }
944
945 if ( k == nseg && i == nseg ) {
946 pp [ 0 ].x = p [ 5 ].x;
947 pp [ 0 ].y = p [ 5 ].y;
948 pp [ 0 ].z = p [ 5 ].z;
949
950 pp [ 1 ].x = p [ 6 ].x;
951 pp [ 1 ].y = p [ 6 ].y;
952 pp [ 1 ].z = p [ 6 ].z;
953
954 #ifdef DRAW_VISIBLE_CONTOUR
955 #ifdef SPHERE_WITH_HOLE
956 double xx = 0.0, yy = 0.0, zz = 0.0, rr, r;
957 for ( int ii = 0; ii < 2; ii++ ) {
958 xx += pp [ ii ].x;
959 yy += pp [ ii ].y;
960 zz += pp [ ii ].z;
961 }
962 xx /= 2.0;
963 yy /= 2.0;
964 zz /= 2.0;
965 rr = xx * xx + yy * yy;
966 r = rr + zz * zz;
967
968 if ( zz < 2.0001 /* || xx < 0.0001 */ || yy < 0.0001 || rr < 1.001 * 1.001 || r < 25.0 || r > 5.98 * 5.98 ) {
969 if ( yy < 2.0001 ) {
970 auto go = CreateLine3D(pp);
971 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
972 EGAttachObject(go, ( EObjectP ) this);
973 EMAddGraphicsToModel(ESIModel(), go);
974 }
975 }
976 #endif
977 #else
978 auto go = CreateLine3D(pp);
979 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
980 EGAttachObject(go, ( EObjectP ) this);
981 EMAddGraphicsToModel(ESIModel(), go);
982 #endif
983 }
984
985 #else
986 auto go = CreateHexahedron(p);
987 EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
988 EGAttachObject(go, ( EObjectP ) this);
989 EMAddGraphicsToModel(ESIModel(), go);
990 #endif
991 }
992 }
993 }
994 } // end loop over knot spans (irules)
995 } else {
996 OOFEM_ERROR("not implemented for nsd = %d", nsd);
997 }
998}
999
1000
1001void drawIGAPatchDeformedGeometry(Element *elem, StructuralElementEvaluator *se, oofegGraphicContext &gc, TimeStep *tStep, UnknownType)
1002{
1003 WCRec p [ 8 ];
1004 int nseg;
1005 FloatArray u;
1006 FloatMatrix N;
1007 FloatArray ur, d;
1008 IntArray lc;
1009 FEInterpolation *interp = elem->giveInterpolation();
1010 double defScale = gc.getDefScale();
1011
1012 if ( !gc.testElementGraphicActivity(elem) ) {
1013 return;
1014 }
1015
1016 EASValsSetLineWidth(OOFEG_DEFORMED_GEOMETRY_WIDTH);
1017 EASValsSetColor( gc.getDeformedElementColor() );
1018 EASValsSetEdgeColor( gc.getElementEdgeColor() );
1019 EASValsSetEdgeFlag(true);
1020 EASValsSetLayer(OOFEG_DEFORMED_GEOMETRY_LAYER);
1021 EASValsSetLineStyle(SOLID_STYLE);
1022 EASValsSetFillStyle(FILL_SOLID);
1023 EASValsSetLineWidth(0);
1024
1025 #ifdef DRAW_MESH
1026 //nseg = 8;
1027 nseg = 4;
1028 #else
1029 nseg = 4;
1030 #endif
1031
1032 const FloatArray *knotVector = interp->giveKnotVector();
1033 int nsd = interp->giveNsd(this->giveGeometryType());
1034
1035 se->computeVectorOf(VM_Total, tStep, u);
1036
1037 if ( nsd == 1 ) {
1038 FloatArray c [ 2 ], cg [ 2 ];
1039
1040
1041 for ( int j = 0; j < 2; j++ ) {
1042 c [ j ].resize(1);
1043 cg [ j ].resize(1);
1044 }
1045
1046 // loop over individual integration rules (i.e., knot spans)
1047 for ( int ir = 0; ir < elem->giveNumberOfIntegrationRules(); ir++ ) {
1048 IntegrationRule *iRule = elem->giveIntegrationRule(ir);
1049 const IntArray *span = iRule->giveKnotSpan();
1050 // divide span locally to get finer geometry rep.
1051 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
1052 for ( int i = 1; i <= nseg; i++ ) {
1053 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1054 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1055
1056 for ( int k = 0; k < 2; k++ ) {
1057 // create a dummy ip's
1058 GaussPoint gp(iRule, 999, c [ k ], 1.0, _PlaneStress);
1059
1060 // compute displacements at gp
1061 se->computeNMatrixAt(N, & gp);
1062
1063 // get local code numbers corresponding to ir
1064 se->giveIntegrationElementLocalCodeNumbers(lc, elem, gp.giveIntegrationRule());
1065 ur.resize( N.giveNumberOfColumns() );
1066 for ( int n = 1; n <= lc.giveSize(); n++ ) {
1067 ur.at(n) = u.at( lc.at(n) );
1068 }
1069
1070 // interpolate displacements
1071 d.beProductOf(N, ur);
1072
1073 interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( elem, iRule->giveKnotSpan() ) );
1074 p [ k ].x = ( FPNum ) ( cg [ k ].at(1) + d.at(1) * defScale );
1075 p [ k ].y = 0.;
1076 p [ k ].z = 0.;
1077 }
1078
1079 auto go = CreateLine3D(p);
1080 EGWithMaskChangeAttributes(WIDTH_MASK | STYLE_MASK | COLOR_MASK | LAYER_MASK, go);
1081 EGAttachObject(go, ( EObjectP ) elem);
1082 EMAddGraphicsToModel(ESIModel(), go);
1083 }
1084 } // end loop over knot spans (irules)
1085 } else if ( nsd == 2 ) {
1086 FloatArray c [ 4 ], cg [ 4 ];
1087
1088 for ( int j = 0; j < 4; j++ ) {
1089 c [ j ].resize(2);
1090 cg [ j ].resize(2);
1091 }
1092
1093 // loop over individual integration rules (i.e., knot spans)
1094 for ( int ir = 0; ir < elem->giveNumberOfIntegrationRules(); ir++ ) {
1095 IntegrationRule *iRule = elem->giveIntegrationRule(ir);
1096 const IntArray *span = iRule->giveKnotSpan();
1097 // divide span locally to get finer geometry rep.
1098 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
1099 double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
1100 for ( int i = 1; i <= nseg; i++ ) {
1101 for ( int j = 1; j <= nseg; j++ ) {
1102 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1103 c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1104 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1105 c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1106 c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1107 c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1108 c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1109 c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1110
1111 for ( int k = 0; k < 4; k++ ) {
1112 // create a dummy ip's
1113 GaussPoint gp(iRule, 999, c [ k ], 1.0, _PlaneStress);
1114
1115 // compute displacements at gp
1116 se->computeNMatrixAt(N, & gp);
1117
1118 // get local code numbers corresponding to ir
1119 se->giveIntegrationElementLocalCodeNumbers(lc, elem, gp.giveIntegrationRule());
1120 ur.resize( N.giveNumberOfColumns() );
1121 for ( int n = 1; n <= lc.giveSize(); n++ ) {
1122 ur.at(n) = u.at( lc.at(n) );
1123 }
1124
1125 // interpolate displacements
1126 d.beProductOf(N, ur);
1127
1128 interp->local2global( cg [ k ], c [ k ], FEIIGAElementGeometryWrapper( elem, iRule->giveKnotSpan() ) );
1129 p [ k ].x = ( FPNum ) ( cg [ k ].at(1) + d.at(1) * defScale );
1130 p [ k ].y = ( FPNum ) ( cg [ k ].at(2) + d.at(2) * defScale );
1131 p [ k ].z = 0.;
1132 }
1133
1134 auto go = CreateQuad3D(p);
1135 EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
1136 EGAttachObject(go, ( EObjectP ) elem);
1137 EMAddGraphicsToModel(ESIModel(), go);
1138 }
1139 }
1140 } // end loop over knot spans (irules)
1141 } else if ( nsd == 3 ) {
1142 FloatArray c [ 8 ], cg [ 8 ];
1143
1144 for ( int j = 0; j < 8; j++ ) {
1145 c [ j ].resize(3);
1146 cg [ j ].resize(3);
1147 }
1148
1149 // loop over individual integration rules (i.e., knot spans)
1150 for ( int ir = 0; ir < elem->giveNumberOfIntegrationRules(); ir++ ) {
1151 IntegrationRule *iRule = elem->giveIntegrationRule(ir);
1152 const IntArray *span = iRule->giveKnotSpan();
1153 // divide span locally to get finer geometry rep.
1154 double du = ( knotVector [ 0 ] [ span->at(1) + 1 ] - knotVector [ 0 ] [ span->at(1) ] ) / nseg;
1155 double dv = ( knotVector [ 1 ] [ span->at(2) + 1 ] - knotVector [ 1 ] [ span->at(2) ] ) / nseg;
1156 double dt = ( knotVector [ 2 ] [ span->at(3) + 1 ] - knotVector [ 2 ] [ span->at(3) ] ) / nseg;
1157 for ( int i = 1; i <= nseg; i++ ) {
1158 for ( int j = 1; j <= nseg; j++ ) {
1159 for ( int k = 1; k <= nseg; k++ ) {
1160 c [ 0 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1161 c [ 0 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1162 c [ 0 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
1163 c [ 1 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1164 c [ 1 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1165 c [ 1 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
1166 c [ 2 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1167 c [ 2 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1168 c [ 2 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
1169 c [ 3 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1170 c [ 3 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1171 c [ 3 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * ( k - 1 );
1172 c [ 4 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1173 c [ 4 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1174 c [ 4 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
1175 c [ 5 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1176 c [ 5 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * ( j - 1 );
1177 c [ 5 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
1178 c [ 6 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * i;
1179 c [ 6 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1180 c [ 6 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
1181 c [ 7 ].at(1) = knotVector [ 0 ] [ span->at(1) ] + du * ( i - 1 );
1182 c [ 7 ].at(2) = knotVector [ 1 ] [ span->at(2) ] + dv * j;
1183 c [ 7 ].at(3) = knotVector [ 2 ] [ span->at(3) ] + dt * k;
1184
1185 for ( int m = 0; m < 8; m++ ) {
1186 // create a dummy ip's
1187 GaussPoint gp(iRule, 999, c [ m ], 1.0, _3dMat);
1188
1189 // compute displacements at gp
1190 se->computeNMatrixAt(N, & gp);
1191
1192 // get local code numbers corresponding to ir
1193 se->giveIntegrationElementLocalCodeNumbers(lc, elem, gp.giveIntegrationRule());
1194 ur.resize( N.giveNumberOfColumns() );
1195 for ( int n = 1; n <= lc.giveSize(); n++ ) {
1196 ur.at(n) = u.at( lc.at(n) );
1197 }
1198
1199 // interpolate displacements
1200 d.beProductOf(N, ur);
1201
1202 interp->local2global( cg [ m ], c [ m ], FEIIGAElementGeometryWrapper( elem, iRule->giveKnotSpan() ) );
1203 p [ m ].x = ( FPNum ) ( cg [ m ].at(1) + d.at(1) * defScale );
1204 p [ m ].y = ( FPNum ) ( cg [ m ].at(2) + d.at(2) * defScale );
1205 p [ m ].z = ( FPNum ) ( cg [ m ].at(3) + d.at(3) * defScale );
1206 }
1207
1208 auto go = CreateHexahedron(p);
1209 EGWithMaskChangeAttributes(WIDTH_MASK | FILL_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
1210 EGAttachObject(go, ( EObjectP ) elem);
1211 EMAddGraphicsToModel(ESIModel(), go);
1212 }
1213 }
1214 }
1215 } // end loop over knot spans (irules)
1216 } else {
1217 OOFEM_ERROR("not implemented for nsd = %d", nsd);
1218 }
1219}
1220
1221
1222
1223
1224#endif
1225} // end namespace oofem
N
#define N(a, b)
oofem::BSplineInterpolation::postInitialize
void postInitialize(ParameterManager &pm, int elnum) override
Definition feibspline.C:86
oofem::Element::dofManArray
IntArray dofManArray
Array containing dofmanager numbers.
Definition element.h:138
oofem::Element::giveInterpolation
virtual FEInterpolation * giveInterpolation() const
Definition element.h:648
oofem::Element::giveNumberOfIntegrationRules
int giveNumberOfIntegrationRules()
Definition element.h:894
oofem::Element::numberOfDofMans
int numberOfDofMans
Number of dofmanagers.
Definition element.h:136
oofem::Element::giveIntegrationRule
virtual IntegrationRule * giveIntegrationRule(int i)
Definition element.h:899
oofem::Element::integrationRulesArray
std::vector< std ::unique_ptr< IntegrationRule > > integrationRulesArray
Definition element.h:157
oofem::Element::giveParallelMode
elementParallelMode giveParallelMode() const
Definition element.h:1139
oofem::Element::numberOfGaussPoints
int numberOfGaussPoints
Definition element.h:175
oofem::Element::initializeFinish
void initializeFinish() override
Definition element.C:786
oofem::Element::giveGeometryType
virtual Element_Geometry_Type giveGeometryType() const =0
oofem::FEInterpolation::giveNumberOfKnotSpans
virtual int giveNumberOfKnotSpans(int dim) const
Definition feinterpol.h:534
oofem::FEInterpolation::postInitialize
virtual void postInitialize(ParameterManager &pm, int elnum)
Definition feinterpol.h:184
oofem::FEInterpolation::giveKnotValues
virtual const FloatArray * giveKnotValues(int dim) const
Definition feinterpol.h:538
oofem::FEInterpolation::giveKnotVector
virtual const FloatArray * giveKnotVector() const
Definition feinterpol.h:530
oofem::FEInterpolation::giveKnotMultiplicity
virtual const IntArray * giveKnotMultiplicity(int dim) const
Definition feinterpol.h:542
oofem::FEInterpolation::giveNsd
virtual int giveNsd(const Element_Geometry_Type) const =0
oofem::FEInterpolation::local2global
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const =0
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition femcmpnn.h:79
oofem::FEMComponent::number
int number
Component number.
Definition femcmpnn.h:77
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Definition floatarray.C:689
oofem::GaussPoint::giveIntegrationRule
IntegrationRule * giveIntegrationRule()
Returns corresponding integration rule to receiver.
Definition gausspoint.h:185
oofem::IGAElement::giveNsd
virtual int giveNsd(const Element_Geometry_Type)=0
oofem::IGAElement::knotSpanParallelMode
IntArray knotSpanParallelMode
Definition iga.h:95
oofem::IGAElement::initializeFinish
void initializeFinish() override
Definition iga.C:70
oofem::IGAElement::IPK_IGAElement_knotSpanParallelMode
static ParamKey IPK_IGAElement_knotSpanParallelMode
Definition iga.h:94
oofem::IGATSplineElement::initializeFinish
void initializeFinish() override
Definition iga.C:230
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::IntArray::giveSize
int giveSize() const
Definition intarray.h:211
oofem::IntegrationRule::giveKnotSpan
virtual const IntArray * giveKnotSpan()
Returns receiver sub patch indices (if apply).
Definition integrationrule.h:245
oofem::ParameterManager::checkIfSet
bool checkIfSet(size_t componentIndex, size_t paramIndex)
Definition parametermanager.h:156
oofem::StructuralElementEvaluator::computeNMatrixAt
virtual void computeNMatrixAt(FloatMatrix &answer, GaussPoint *gp)=0
oofem::StructuralElementEvaluator::computeVectorOf
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Definition structuralelementevaluator.h:118
oofem::StructuralElementEvaluator::giveIntegrationElementLocalCodeNumbers
virtual int giveIntegrationElementLocalCodeNumbers(IntArray &answer, Element *elem, IntegrationRule *ie)
Definition structuralelementevaluator.C:81
oofem::TSplineInterpolation::setNumberOfControlPoints
void setNumberOfControlPoints(int num)
Definition feitspline.h:78
oofem::TSplineInterpolation::initializeFrom
void initializeFrom(InputRecord &ir, ParameterManager &pm, int elnum, int priority) override
Initializes receiver according to object description stored in input record.
Definition feitspline.C:46
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
_IFT_IGAElement_KnotSpanParallelMode
#define _IFT_IGAElement_KnotSpanParallelMode
Definition iga.h:45
oofem::elementParallelMode
elementParallelMode
Definition element.h:87
oofem::Element_remote
@ Element_remote
Element in active domain is only mirror of some remote element.
Definition element.h:89
oofem::Element_local
@ Element_local
Element is local, there are no contributions from other domains to this element.
Definition element.h:88
oofem::drawIGAPatchDeformedGeometry
void drawIGAPatchDeformedGeometry(Element *elem, StructuralElementEvaluator *se, oofegGraphicContext &gc, TimeStep *tStep, UnknownType)
Definition iga.C:1001
gc
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
OOFEG_DEFORMED_GEOMETRY_LAYER
#define OOFEG_DEFORMED_GEOMETRY_LAYER
Definition oofeggraphiccontext.h:102
OOFEG_DEFORMED_GEOMETRY_WIDTH
#define OOFEG_DEFORMED_GEOMETRY_WIDTH
Definition oofeggraphiccontext.h:95
OOFEG_RAW_GEOMETRY_WIDTH
#define OOFEG_RAW_GEOMETRY_WIDTH
Definition oofeggraphiccontext.h:94
OOFEG_RAW_GEOMETRY_LAYER
#define OOFEG_RAW_GEOMETRY_LAYER
Definition oofeggraphiccontext.h:101
PM_UPDATE_PARAMETER
#define PM_UPDATE_PARAMETER(_val, _pm, _ir, _componentnum, _paramkey, _prio)
Definition parametermanager.h:56
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