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