OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
supgelement2.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 "supgelement2.h"
36 #include "domain.h"
37 #include "timestep.h"
38 #include "bodyload.h"
39 #include "intarray.h"
40 #include "floatarray.h"
41 #include "floatmatrix.h"
42 #include "fluidmodel.h"
43 #include "fluiddynamicmaterial.h"
44 #include "fluidcrosssection.h"
45 #include "integrationrule.h"
46 #include "node.h"
47 #include "dof.h"
48 
49 #ifdef __OOFEG
50  #include "oofeggraphiccontext.h"
51  #include "connectivitytable.h"
52 #endif
53 
54 namespace oofem {
56  SUPGElement(n, aDomain)
57  // Constructor. Creates an element with number n, belonging to aDomain.
58 { }
59 
60 
62 // Destructor.
63 { }
64 
67 {
68  //IRResultType result; // Required by IR_GIVE_FIELD macro
69 
71 }
72 
73 
74 void
76 {
78 }
79 
80 
81 void
83  CharType mtrx, TimeStep *tStep)
84 //
85 // returns characteristics matrix of receiver according to mtrx
86 //
87 {
88 #if 0
89  if ( mtrx == TangentStiffnessMatrix ) {
90  // support for stokes solver
91  IntArray vloc, ploc;
92  FloatMatrix h;
93  int size = this->computeNumberOfDofs();
94  this->giveLocalVelocityDofMap(vloc);
95  this->giveLocalPressureDofMap(ploc);
96  answer.resize(size, size);
97  answer.zero();
98  this->computeAdvectionDerivativeTerm_MB(h, tStep);
99  answer.assemble(h, vloc);
100  this->computeDiffusionDerivativeTerm_MB(h, TangentStiffness, tStep);
101  answer.assemble(h, vloc);
102  this->computePressureTerm_MB(h, tStep);
103  answer.assemble(h, vloc, ploc);
104  this->computeLinearAdvectionTerm_MC(h, tStep);
105  answer.assemble(h, ploc, vloc);
106  this->computeAdvectionDerivativeTerm_MC(h, tStep);
107  answer.assemble(h, ploc, vloc);
108  this->computeDiffusionDerivativeTerm_MC(h, tStep);
109  answer.assemble(h, ploc, vloc);
110  this->computePressureTerm_MC(h, tStep);
111  answer.assemble(h, ploc);
112  this->computeLSICStabilizationTerm_MB(h, tStep);
113  answer.assemble(h, vloc);
114  } else
115 #endif
116  {
117  OOFEM_ERROR("Unknown Type of characteristic mtrx.");
118  }
119 }
120 
121 
122 void
124  TimeStep *tStep)
125 //
126 // returns characteristics vector of receiver according to requested type
127 //
128 {
129  if ( mtrx == ExternalForcesVector ) {
130  // stokes flow
131  IntArray vloc, ploc;
132  FloatArray h;
133  int size = this->computeNumberOfDofs();
134  this->giveLocalVelocityDofMap(vloc);
135  this->giveLocalPressureDofMap(ploc);
136  answer.resize(size);
137  answer.zero();
138  this->computeBCRhsTerm_MB(h, tStep);
139  answer.assemble(h, vloc);
140  this->computeBCRhsTerm_MC(h, tStep);
141  answer.assemble(h, ploc);
142  } else
143 #if 0
144  if ( mtrx == InternalForcesVector ) {
145  // stokes flow
146  IntArray vloc, ploc;
147  FloatArray h;
148  int size = this->computeNumberOfDofs();
149  this->giveLocalVelocityDofMap(vloc);
150  this->giveLocalPressureDofMap(ploc);
151  answer.resize(size);
152  answer.zero();
153  this->computeAdvectionTerm_MB(h, tStep);
154  answer.assemble(h, vloc);
155  this->computeAdvectionTerm_MC(h, tStep);
156  answer.assemble(h, ploc);
157  this->computeDiffusionTerm_MB(h, tStep);
158  answer.assemble(h, vloc);
159  this->computeDiffusionTerm_MC(h, tStep);
160  answer.assemble(h, ploc);
161 
162  FloatMatrix m1;
163  FloatArray v, p;
164  // add lsic stabilization term
165  this->giveCharacteristicMatrix(m1, LSICStabilizationTerm_MB, tStep);
166  //m1.times( lscale / ( dscale * uscale * uscale ) );
167  this->computeVectorOfVelocities(VM_Total, tStep, v);
168  h.beProductOf(m1, v);
169  answer.assemble(h, vloc);
170  this->giveCharacteristicMatrix(m1, LinearAdvectionTerm_MC, tStep);
171  //m1.times( 1. / ( dscale * uscale ) );
172  h.beProductOf(m1, v);
173  answer.assemble(h, ploc);
174 
175 
176  // add pressure term
177  this->giveCharacteristicMatrix(m1, PressureTerm_MB, tStep);
178  this->computeVectorOfPressures(VM_Total, tStep, v);
179  h.beProductOf(m1, v);
180  answer.assemble(h, vloc);
181 
182  // pressure term
183  this->giveCharacteristicMatrix(m1, PressureTerm_MC, tStep);
184  this->computeVectorOfPressures(VM_Total, tStep, v);
185  h.beProductOf(m1, v);
186  answer.assemble(h, ploc);
187  } else
188 #endif
189  {
190  OOFEM_ERROR("Unknown Type of characteristic mtrx.");
191  }
192 }
193 
194 
195 int
197 //
198 // check internal consistency
199 // mainly tests, whether material and crossSection data
200 // are safe for conversion to "Structural" versions
201 //
202 {
203  int result = 1;
204  return result;
205 }
206 
207 
208 void
210 {
211  FloatArray stress;
212 
213  // force updating strains & stresses
214  for ( auto &iRule: integrationRulesArray ) {
215  for ( GaussPoint *gp: *iRule ) {
216  computeDeviatoricStress(stress, gp, tStep);
217  }
218  }
219 }
220 
221 
222 #ifdef __OOFEG
223 int
225  int node, TimeStep *tStep)
226 {
227  int indx = 1;
228  Node *n = this->giveNode(node);
229 
230  if ( type == IST_Velocity ) {
231  answer.resize( this->giveSpatialDimension() );
232  std::vector< Dof* >::const_iterator dofindx;
233  if ( ( dofindx = n->findDofWithDofId(V_u) ) != n->end() ) {
234  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
235  } else if ( ( dofindx = n->findDofWithDofId(V_v) ) != n->end() ) {
236  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
237  } else if ( ( dofindx = n->findDofWithDofId(V_w) ) != n->end() ) {
238  answer.at(indx++) = (*dofindx)->giveUnknown(VM_Total, tStep);
239  }
240 
241  return 1;
242  } else if ( type == IST_Pressure ) {
243  auto dofindx = n->findDofWithDofId(P_f);
244  if ( dofindx != n->end() ) {
245  answer.resize(1);
246  answer.at(1) = (*dofindx)->giveUnknown(VM_Total, tStep);
247  return 1;
248  } else {
249  return 0;
250  }
251  } else {
252  return Element :: giveInternalStateAtNode(answer, type, mode, node, tStep);
253  }
254 }
255 
256 #endif
257 
258 void
260 {
261  FloatMatrix n, b;
262 
263  answer.clear();
264 
265  int rule = 2;
266  /* consistent part + supg stabilization term */
267  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
268  this->computeNuMatrix(n, gp);
269  this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
270  double dV = this->computeVolumeAround(gp);
271  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
272  /* consistent part */
273  answer.plusProductUnsym(n, n, rho * dV);
274  /* supg stabilization */
275  answer.plusProductUnsym(b, n, rho * t_supg * dV);
276  }
277 }
278 
279 void
281 {
282  FloatMatrix n, b, bn;
283  FloatArray u, v;
284 
285  answer.clear();
286 
287  this->computeVectorOfVelocities(VM_Total, tStep, u);
288 
289  int rule = 2;
290  /* consistent part + supg stabilization term */
291  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
292  this->computeNuMatrix(n, gp);
293  this->computeUDotGradUMatrix( bn, gp, tStep->givePreviousStep() );
294  this->computeUDotGradUMatrix(b, gp, tStep);
295  v.beProductOf(b, u);
296  double dV = this->computeVolumeAround(gp);
297  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
298  /* consistent part */
299  answer.plusProduct(n, v, rho * dV);
300 
301  /* supg stabilization */
302  answer.plusProduct(bn, v, t_supg * rho * dV);
303  }
304 }
305 
306 void
308 {
309  FloatMatrix n, b, bn, grad_u, grad_uN;
310 
311  answer.clear();
312 
313  int rule = 2;
314  /* consistent part + supg stabilization term */
315  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
316  this->computeNuMatrix(n, gp);
317  this->computeUDotGradUMatrix( bn, gp, tStep->givePreviousStep() );
318  this->computeUDotGradUMatrix(b, gp, tStep);
319  double dV = this->computeVolumeAround(gp);
320  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
321 
322  this->computeGradUMatrix(grad_u, gp, tStep);
323 
324  /* consistent part */
325  answer.plusProductUnsym(n, b, rho * dV);
326 
327  grad_uN.beProductOf(grad_u, n);
328  answer.plusProductUnsym(n, grad_uN, rho * dV);
329  /* supg stabilization */
330  answer.plusProductUnsym(bn, b, t_supg * rho * dV);
331  answer.plusProductUnsym(bn, grad_uN, t_supg * rho * dV);
332  }
333 }
334 
335 void
337 {
338  FloatArray u, eps, stress, dDB_u;
339  FloatMatrix b, un_gu, dDB;
340  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
341 
342  answer.clear();
343 
344  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
345  this->computeVectorOfVelocities(VM_Total, tStep, u);
346 
347  int rule = 1;
348  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
349  double dV = this->computeVolumeAround(gp);
350  this->computeBMatrix(b, gp);
351  this->computeDivTauMatrix(dDB, gp, tStep);
352  this->computeUDotGradUMatrix( un_gu, gp, tStep->givePreviousStep() );
353  eps.beProductOf(b, u);
354  mat->computeDeviatoricStressVector(stress, gp, eps, tStep);
355  dDB_u.beProductOf(dDB, u);
356  /* consistent part */
357  answer.plusProduct(b, stress, dV / Re);
358 
359  /* SUPG term */
360  answer.plusProduct(un_gu, dDB_u, ( -1.0 ) * t_supg * dV);
361  }
362 }
363 
364 void
366 {
367  FloatMatrix _db, _d, _b, dDB, un_gu;
368  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
369 
370  answer.clear();
371 
372  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
373 
374  int rule = 1;
375  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
376  double dV = this->computeVolumeAround(gp);
377  this->computeBMatrix(_b, gp);
378  mat->giveDeviatoricStiffnessMatrix(_d, mode, gp, tStep);
379 
380  this->computeDivTauMatrix(dDB, gp, tStep);
381  this->computeUDotGradUMatrix( un_gu, gp, tStep->givePreviousStep() );
382  /* standard term */
383  _db.beProductOf(_d, _b);
384  answer.plusProductUnsym(_b, _db, dV / Re); //answer.plusProduct (_b,_db,area);
385 
386  /* SUPG term */
387  answer.plusProductUnsym( un_gu, dDB, t_supg * dV * ( -1.0 ) * ( 1. / Re ) );
388  }
389 }
390 
391 
392 void
394 {
395  FloatMatrix gu, np, b;
396 
397  answer.clear();
398 
399  int rule = 1;
400  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
401  double dV = this->computeVolumeAround(gp);
402  this->computeDivUMatrix(gu, gp);
403  this->computeNpMatrix(np, gp);
404 
405  /*alternative computing*/
406  //this->computeNuMatrix(gu, gp);
407  //this->computeGradPMatrix(np, gp);
408 
409  /* standard term */
410  answer.plusProductUnsym(gu, np, ( -1.0 ) * dV);
411  }
412 
413  for ( GaussPoint *gp: *this->integrationRulesArray [ 1 ] ) {
414  double dV = this->computeVolumeAround(gp);
415  this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
416  this->computeGradPMatrix(np, gp);
417 
418  // supg term
419  answer.plusProductUnsym(b, np, t_supg * dV);
420  }
421 }
422 void
424 {
425  FloatMatrix b;
426 
427  answer.clear();
428 
429  int rule = 0;
430  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
431  double dV = this->computeVolumeAround(gp);
432  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
433  this->computeDivUMatrix(b, gp);
434 
435  answer.plusProductSymmUpper(b, b, dV * rho * t_lsic);
436  }
437 
438  answer.symmetrized();
439 }
440 
441 void
443 {
444  FloatMatrix gu, np;
445 
446  answer.clear();
447 
448  int rule = 0;
449  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
450  double dV = this->computeVolumeAround(gp);
451  this->computeDivUMatrix(gu, gp);
452  this->computeNpMatrix(np, gp);
453 
454  /* standard term */
455  answer.plusProductUnsym(np, gu, dV);
456  }
457 }
458 
459 void
461 {
462  // N_epsilon (due to PSPG stabilization)
463  FloatMatrix g, b;
464  FloatArray u, v;
465 
466  answer.clear();
467 
468  int rule = 1;
469  this->computeVectorOfVelocities(VM_Total, tStep, u);
470 
471  /* pspg stabilization term */
472  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
473  this->computeGradPMatrix(g, gp);
474  this->computeUDotGradUMatrix(b, gp, tStep);
475  v.beProductOf(b, u);
476  double dV = this->computeVolumeAround(gp);
477  answer.plusProduct(g, v, t_pspg * dV);
478  }
479 }
480 
481 
482 void
484 {
485  FloatMatrix g, b;
486 
487  answer.clear();
488 
489  int rule = 1;
490  /* pspg stabilization term */
491  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
492  this->computeGradPMatrix(g, gp);
493  this->computeUDotGradUMatrix(b, gp, tStep);
494  double dV = this->computeVolumeAround(gp);
495  answer.plusProductUnsym(g, b, dV * t_pspg);
496  }
497 }
498 
499 void
501 {
502  FloatMatrix dDB, g;
503 
504  answer.clear();
505 
506  int rule = 1;
507  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
508  double dV = this->computeVolumeAround(gp);
509  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
510 
511  this->computeDivTauMatrix(dDB, gp, tStep);
512  this->computeGradPMatrix(g, gp);
513 
514  answer.plusProductUnsym(g, dDB, ( -1.0 ) * dV * t_pspg / rho);
515  }
516 
517  answer.zero();
518 }
519 
520 void
522 {
523  answer.clear();
524 
525  /*
526  * FloatMatrix dDB, _d, g;
527  * double sum, dV, coeff, rho, isd, nsd = this->giveNumberOfSpatialDimensions();
528  * GaussPoint *gp;
529  * int i,k;
530  * FloatArray u, dDB_u;
531  *
532  * this->computeVectorOfVelocities(VM_Total, tStep, u);
533  *
534  * for ( GaussPoint *gp: *this->integrationRulesArray [ 1 ] ) {
535  * dV = this->computeVolumeAround(gp);
536  * rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
537  *
538  * coeff = (-1.0) * dV * t_pspg / rho;
539  *
540  * this->computeGradPMatrix(g, gp);
541  * this->computeDivTauMatrix(dDB, gp, tStep);
542  *
543  * dDB_u.beProductOf(dDB, u);
544  *
545  * for ( i = 1; i <= pndofs; i++ ) {
546  * for ( sum = 0, isd = 1; isd <= nsd; isd++ ) {
547  * sum += g.at(isd, i) * dDB_u.at(isd);
548  * }
549  *
550  * answer.at(i) += coeff * sum;
551  * }
552  *
553  *
554  *
555  * }
556  */
557 }
558 
559 
560 void
562 {
563  FloatMatrix g, n;
564 
565  answer.clear();
566  // pspg stabilization term: M_\epsilon term
567 
568  int rule = 0;
569  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
570  this->computeGradPMatrix(g, gp);
571  this->computeNuMatrix(n, gp);
572  double dV = this->computeVolumeAround(gp);
573  answer.plusProductUnsym(g, n, dV * t_pspg);
574  }
575 }
576 
577 void
579 {
580  FloatMatrix g;
581 
582  answer.clear();
583 
584  int rule = 0;
585  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
586  this->computeGradPMatrix(g, gp);
587  double dV = this->computeVolumeAround(gp);
588  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
589  answer.plusProductSymmUpper(g, g, dV * t_pspg / rho);
590  }
591 
592  answer.symmetrized();
593 }
594 
595 
596 void
598 {
599  int nLoads;
600 
601  answer.clear();
602 
603  int rule = 0;
604  FloatArray gVector, helpLoadVector;
605  FloatMatrix b, nu;
606 
607  // add body load (gravity) termms
608  nLoads = this->giveBodyLoadArray()->giveSize();
609  for ( int i = 1; i <= nLoads; i++ ) {
610  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
611  bcGeomType ltype = load->giveBCGeoType();
612  if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
613  load->computeComponentArrayAt(gVector, tStep, VM_Total);
614  if ( gVector.giveSize() ) {
615  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
616  this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
617  this->computeNuMatrix(nu, gp);
618  double dV = this->computeVolumeAround(gp);
619  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
620  answer.plusProduct(b, gVector, t_supg * rho * dV);
621  answer.plusProduct(nu, gVector, rho * dV);
622  }
623  }
624  }
625  }
626 
627  // integrate tractions
628  // if no traction bc applied but side marked as with traction load
629  // then zero traction is assumed !!!
630 
631  // loop over boundary load array
632  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
633  for ( int i = 1; i <= nLoads; i++ ) {
634  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
635  int id = boundaryLoadArray.at(i * 2);
636  Load *load = domain->giveLoad(n);
637  bcGeomType ltype = load->giveBCGeoType();
638  if ( ltype == EdgeLoadBGT ) {
639  this->computeEdgeLoadVector_MB(helpLoadVector, load, id, tStep);
640  if ( helpLoadVector.giveSize() ) {
641  answer.add(helpLoadVector);
642  }
643  } else if ( ltype == SurfaceLoadBGT ) {
644  this->computeSurfaceLoadVector_MB(helpLoadVector, load, id, tStep);
645  if ( helpLoadVector.giveSize() ) {
646  answer.add(helpLoadVector);
647  }
648  } else {
649  OOFEM_ERROR("unsupported load type class");
650  }
651  }
652 }
653 
654 
655 void
657 {
658  int nLoads;
659  FloatArray gVector, helpLoadVector;
660  FloatMatrix g;
661 
662  int rule = 1;
663 
664  answer.clear();
665 
666  nLoads = this->giveBodyLoadArray()->giveSize();
667  for ( int i = 1; i <= nLoads; i++ ) {
668  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
669  bcGeomType ltype = load->giveBCGeoType();
670  if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
671  load->computeComponentArrayAt(gVector, tStep, VM_Total);
672  if ( gVector.giveSize() ) {
673  for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
674  this->computeGradPMatrix(g, gp);
675  double dV = this->computeVolumeAround(gp);
676  answer.plusProduct(g, gVector, t_pspg * dV);
677  }
678  }
679  }
680  }
681 
682  // integrate tractions
683  // if no traction bc applied but side marked as with traction load
684  // then zero traction is assumed !!!
685 
686  // loop over boundary load array
687  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
688  for ( int i = 1; i <= nLoads; i++ ) {
689  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
690  int id = boundaryLoadArray.at(i * 2);
691  Load *load = domain->giveLoad(n);
692  bcGeomType ltype = load->giveBCGeoType();
693  if ( ltype == EdgeLoadBGT ) {
694  this->computeEdgeLoadVector_MC(helpLoadVector, load, id, tStep);
695  if ( helpLoadVector.giveSize() ) {
696  answer.add(helpLoadVector);
697  }
698  } else if ( ltype == SurfaceLoadBGT ) {
699  this->computeSurfaceLoadVector_MC(helpLoadVector, load, id, tStep);
700  if ( helpLoadVector.giveSize() ) {
701  answer.add(helpLoadVector);
702  }
703  } else {
704  OOFEM_ERROR("unsupported load type class");
705  }
706  }
707 }
708 
709 
710 void
712 {
713  if ( type != ExternalForcesVector ) {
714  answer.clear();
715  return;
716  }
717 
718  int rule = 1;
719  FloatArray un, nV;
720  FloatArray mb, mc;
721 
722  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
723 
724  if ( load->giveBCValType() == ForceLoadBVT ) {
725  FloatMatrix b, g, nu;
726  FloatArray gVector;
727  load->computeComponentArrayAt(gVector, tStep, VM_Total);
728 
729  for ( auto &gp : *integrationRulesArray [ rule ] ) {
730  double dV = this->computeVolumeAround(gp);
731  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
732 
733  this->computeNuMatrix(nu, gp);
734  mb.plusProduct(nu, gVector, rho * dV);
735  if ( t_supg != 0 ) {
736  this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
737  mb.plusProduct(b, gVector, t_supg * rho * dV);
738  }
739 
740  if ( t_pspg != 0. ) {
741  this->computeGradPMatrix(g, gp);
742  mc.plusProduct(g, gVector, t_pspg * dV);
743  }
744  }
745  }
746 
747  IntArray vloc, ploc;
748  this->giveLocalVelocityDofMap(vloc);
749  this->giveLocalPressureDofMap(ploc);
750  answer.resize(this->computeNumberOfDofs());
751  answer.zero();
752  answer.assemble(mb, vloc);
753  answer.assemble(mc, ploc);
754 }
755 
756 void
758 {
759  OOFEM_ERROR("not implemented");
760 }
761 
762 void
764 {
765  OOFEM_ERROR("not implemented");
766 }
767 
768 void
770 {
771  OOFEM_ERROR("not implemented");
772 }
773 
774 void
776 {
777  OOFEM_ERROR("not implemented");
778 }
779 
780 
781 void
783 {
784  FloatArray u;
785  FloatMatrix b;
786  this->computeVectorOfVelocities(VM_Total, tStep, u);
787 
788  this->computeBMatrix(b, gp);
789  answer.beProductOf(b, u);
790 }
791 
792 void
794 {
795  FloatArray eps;
796 
797  // compute deviatoric strain
798  this->computeDeviatoricStrain(eps, gp, tStep);
799  // call material to compute stress
800  static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeDeviatoricStressVector(answer, gp, eps, tStep);
801 }
802 } // end namespace oofem
CrossSection * giveCrossSection()
Definition: element.C:495
InternalStateType
Type representing the physical meaning of element or constitutive model internal variable.
IntArray * giveBoundaryLoadArray()
Returns array containing load numbers of boundary loads acting on element.
Definition: element.C:381
Class and object Domain.
Definition: domain.h:115
virtual void computeEdgeLoadVector_MB(FloatArray &answer, Load *load, int id, TimeStep *tStep)
Definition: supgelement2.C:757
virtual int computeNumberOfDofs()
Computes or simply returns total number of element&#39;s local DOFs.
Definition: element.h:424
Fluid cross-section.
Abstract base class for all fluid materials.
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
virtual void computeLinearAdvectionTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes the linear advection term for mass conservation equation.
Definition: supgelement2.C:442
bcGeomType
Type representing the geometric character of loading.
Definition: bcgeomtype.h:40
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
Class implementing element body load, acting over whole element volume (e.g., the dead weight)...
Definition: bodyload.h:49
virtual void computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes Rhs terms due to boundary conditions.
Definition: supgelement2.C:656
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
Base class for fluid problems.
Definition: fluidmodel.h:46
virtual void computeNpMatrix(FloatMatrix &answer, GaussPoint *gp)=0
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
virtual void computeAdvectionTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes advection terms for mass conservation equation.
Definition: supgelement2.C:460
virtual void computeGradPMatrix(FloatMatrix &answer, GaussPoint *gp)=0
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: supgelement2.C:66
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: supgelement.C:80
virtual void computeDeviatoricStrain(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Definition: supgelement2.C:782
virtual void giveLocalPressureDofMap(IntArray &map)
Definition: supgelement.h:209
virtual void computeEdgeLoadVector_MC(FloatArray &answer, Load *load, int id, TimeStep *tStep)
Definition: supgelement2.C:769
virtual void updateInternalState(TimeStep *tStep)
Updates element state after equilibrium in time step has been reached.
Definition: supgelement2.C:209
virtual void computeComponentArrayAt(FloatArray &answer, TimeStep *tStep, ValueModeType mode)
Computes boundary condition value - its components values at given time.
Definition: load.C:82
General stabilized SUPG/PSPG element for CFD analysis.
Definition: supgelement.h:58
virtual void computeLSICStabilizationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes SLIC stabilization term for momentum balance equation(s).
Definition: supgelement2.C:423
SUPGElement2(int n, Domain *aDomain)
Definition: supgelement2.C:55
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
MatResponseMode
Describes the character of characteristic material matrix.
virtual void computeDiffusionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes diffusion derivative terms for mass conservation equation.
Definition: supgelement2.C:500
virtual void computePressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes pressure terms for mass conservation equation.
Definition: supgelement2.C:578
virtual void computeSurfaceLoadVector_MB(FloatArray &answer, Load *load, int id, TimeStep *tStep)
Definition: supgelement2.C:763
Distributed body load.
Definition: bcgeomtype.h:43
virtual void computeAccelerationTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes acceleration terms for mass conservation equation.
Definition: supgelement2.C:561
virtual void computeDivTauMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)=0
virtual void computeAdvectionDerivativeTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes the derivative of advection terms for momentum balance equations(s) with respect to nodal ve...
Definition: supgelement2.C:307
virtual void computePressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes pressure terms for momentum balance equations(s).
Definition: supgelement2.C:393
void plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:780
virtual int checkConsistency()
Performs consistency check.
Definition: supgelement2.C:196
virtual void computeVectorOfPressures(ValueModeType mode, TimeStep *tStep, FloatArray &pressures)
Definition: fmelement.C:55
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
Definition: element.h:518
virtual void computeDeviatoricStressVector(FloatArray &stress_dev, double &epsp_vol, GaussPoint *gp, const FloatArray &eps, double pressure, TimeStep *tStep)
Computes the deviatoric stress vector and volumetric strain rate from given deviatoric strain and pre...
Distributed edge load.
Definition: bcgeomtype.h:44
virtual void computeVectorOfVelocities(ValueModeType mode, TimeStep *tStep, FloatArray &velocities)
Definition: fmelement.C:48
virtual void computeSurfaceLoadVector_MC(FloatArray &answer, Load *load, int id, TimeStep *tStep)
Definition: supgelement2.C:775
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual void computeUDotGradUMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)=0
double t_supg
Stabilization coefficients, updated for each solution step in updateStabilizationCoeffs() ...
Definition: supgelement.h:70
IntArray * giveBodyLoadArray()
Returns array containing load numbers of loads acting on element.
Definition: element.C:372
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
IntArray bodyLoadArray
Array containing indexes of loads (body loads and boundary loads are kept separately), that apply on receiver.
Definition: element.h:160
virtual void computeAdvectionTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes nonlinear advection terms for momentum balance equations(s).
Definition: supgelement2.C:280
void plusProductSymmUpper(const FloatMatrix &a, const FloatMatrix &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:698
virtual void computeGradUMatrix(FloatMatrix &answer, GaussPoint *gp, TimeStep *tStep)=0
virtual int giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode, int node, TimeStep *tStep)
Returns internal state variable (like stress,strain) at node of element in Reduced form...
Definition: element.C:1661
virtual void computeDivUMatrix(FloatMatrix &answer, GaussPoint *gp)=0
TimeStep * givePreviousStep()
Returns pointer to previous solution step.
Definition: timestep.C:114
virtual void computeAdvectionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes the derivative of advection terms for mass conservation equation with respect to nodal veloc...
Definition: supgelement2.C:483
int giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode, int node, TimeStep *tStep)
Returns internal state variable (like stress,strain) at node of element in Reduced form...
Definition: supgelement2.C:224
virtual void computeAccelerationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes acceleration terms (generalized mass matrix with stabilization terms) for momentum balance e...
Definition: supgelement2.C:259
Distributed surface load.
Definition: bcgeomtype.h:45
virtual void computeDiffusionDerivativeTerm_MB(FloatMatrix &answer, MatResponseMode mode, TimeStep *tStep)
Computes the derivative of diffusion terms for momentum balance equations(s) with respect to nodal ve...
Definition: supgelement2.C:365
virtual void giveLocalVelocityDofMap(IntArray &map)
Definition: supgelement.h:208
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual void computeDeviatoricStress(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
Definition: supgelement2.C:793
virtual int giveSpatialDimension()
Returns the element spatial dimension (1, 2, or 3).
Definition: element.C:1347
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
virtual void computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes Rhs terms due to boundary conditions.
Definition: supgelement2.C:597
void assemble(const FloatArray &fe, const IntArray &loc)
Assembles the array fe (typically, the load vector of a finite element) into the receiver, using loc as location array.
Definition: floatarray.C:551
CharType
Definition: chartype.h:87
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
virtual void computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes diffusion terms for momentum balance equations(s).
Definition: supgelement2.C:336
std::vector< Dof * >::const_iterator findDofWithDofId(DofIDItem dofID) const
Finds index of DOF with required physical meaning of receiver.
Definition: dofmanager.C:266
virtual ~SUPGElement2()
Definition: supgelement2.C:61
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual bcGeomType giveBCGeoType() const
Returns geometry character of boundary condition.
virtual void computeNuMatrix(FloatMatrix &answer, GaussPoint *gp)=0
virtual void computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
Computes the contribution of the given body load (volumetric).
Definition: supgelement2.C:711
Class representing the a dynamic Input Record.
std::vector< std::unique_ptr< IntegrationRule > > integrationRulesArray
List of integration rules of receiver (each integration rule contains associated integration points a...
Definition: element.h:170
virtual void computeBMatrix(FloatMatrix &anwer, GaussPoint *gp)=0
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Load is base abstract class for all loads.
Definition: load.h:61
virtual void computeDiffusionTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes diffusion terms for mass conservation equation.
Definition: supgelement2.C:521
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: supgelement.C:65
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: supgelement2.C:75
std::vector< Dof * >::iterator end()
Definition: dofmanager.h:158
void assemble(const FloatMatrix &src, const IntArray &loc)
Assembles the contribution using localization array into receiver.
Definition: floatmatrix.C:251
void clear()
Sets size of receiver to be an empty matrix. It will have zero rows and zero columns size...
Definition: floatmatrix.h:516
virtual bcValType giveBCValType() const
Returns receiver load type.
Class implementing node in finite element mesh.
Definition: node.h:87
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
Load * giveLoad(int n)
Service for accessing particular domain load.
Definition: domain.C:222
void symmetrized()
Initializes the lower half of the receiver according to the upper half.
Definition: floatmatrix.C:1576
Class representing integration point in finite element program.
Definition: gausspoint.h:93
IntArray boundaryLoadArray
Definition: element.h:160
virtual void giveCharacteristicMatrix(FloatMatrix &answer, CharType, TimeStep *tStep)
Computes characteristic matrix of receiver of requested type in given time step.
Definition: supgelement2.C:82
Class representing solution step.
Definition: timestep.h:80
InternalStateMode
Determines the mode of internal variable.
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
virtual void giveCharacteristicVector(FloatArray &answer, CharType, ValueModeType, TimeStep *tStep)
Computes characteristic vector of receiver of requested type in given time step.
Definition: supgelement2.C:123
void plusProduct(const FloatMatrix &b, const FloatArray &s, double dV)
Adds the product .
Definition: floatarray.C:226
virtual void giveDeviatoricStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)=0
Computes the deviatoric stiffness; .

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:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011