OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
tr1_2d_supg_axi.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 "tr1_2d_supg_axi.h"
36 #include "fei2dtrlin.h"
37 #include "fluidmodel.h"
38 #include "node.h"
39 #include "material.h"
40 #include "gausspoint.h"
41 #include "gaussintegrationrule.h"
42 #include "floatmatrix.h"
43 #include "floatarray.h"
44 #include "domain.h"
45 #include "mathfem.h"
46 #include "engngm.h"
47 #include "timestep.h"
48 #include "bodyload.h"
49 #include "boundaryload.h"
50 #include "fluiddynamicmaterial.h"
51 #include "fluidcrosssection.h"
52 #include "crosssection.h"
53 #include "classfactory.h"
54 
55 #ifdef __OOFEG
56  #include "oofeggraphiccontext.h"
57  #include "connectivitytable.h"
58 #endif
59 
60 namespace oofem {
61 REGISTER_Element(TR1_2D_SUPG_AXI);
62 
64  // Constructor.
65 { }
66 
68 // Destructor
69 { }
70 
71 void
73 // Sets up the array containing the four Gauss points of the receiver.
74 {
75  if ( integrationRulesArray.size() == 0 ) {
76  integrationRulesArray.resize( 1 );
77  integrationRulesArray [ 0 ].reset( new GaussIntegrationRule(1, this, 1, 3) );
79  }
80 }
81 
82 
83 void
85 {
86  answer.resize(6, 6);
87  answer.zero();
88  FloatArray un, n;
89  double _val, u, v, dV, rho;
90  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
91 
92  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
93  dV = this->computeVolumeAround(gp);
94  rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
95  this->computeNVector(n, gp);
96 
97  for ( int i = 1; i <= 3; i++ ) {
98  for ( int j = 1; j <= 3; j++ ) {
99  /* consistent mass */
100  _val = n.at(i) * n.at(j) * rho * dV;
101  answer.at(2 * i - 1, 2 * j - 1) += _val;
102  answer.at(2 * i, 2 * j) += _val;
103 
104  /* SUPG stabilization term */
105  u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
106  v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
107  answer.at(2 * i - 1, 2 * j - 1) += ( u * b [ i - 1 ] + v * c [ j - 1 ] ) * n.at(j) * rho * t_supg * dV;
108  answer.at(2 * i, 2 * j) += ( u * b [ i - 1 ] + v * c [ j - 1 ] ) * n.at(j) * rho * t_supg * dV;
109  }
110  }
111  }
112 }
113 
114 void
116 {
117  answer.resize(6);
118  answer.zero();
119 
120  double dudx, dudy, dvdx, dvdy, _u, _v;
121  FloatArray u, un, n;
122 
123  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
124  this->computeVectorOfVelocities(VM_Total, tStep, u);
125 
126  dudx = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
127  dudy = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
128  dvdx = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
129  dvdy = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
130 
131  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
132  double dV = this->computeVolumeAround(gp);
133  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
134  this->computeNVector(n, gp);
135 
136  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
137  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
138 
139  // standard galerkin term
140  for ( int i = 1; i <= 3; i++ ) {
141  answer.at(2 * i - 1) += rho * n.at(i) * ( _u * dudx + _v * dudy ) * dV;
142  answer.at(2 * i) += rho * n.at(i) * ( _u * dvdx + _v * dvdy ) * dV;
143  }
144 
145  // supg stabilization term
146  for ( int i = 1; i <= 3; i++ ) {
147  answer.at(2 * i - 1) += ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( _u * dudx + _v * dudy ) * rho * t_supg * dV;
148  answer.at(2 * i) += ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( _u * dvdx + _v * dvdy ) * rho * t_supg * dV;
149  }
150  }
151 }
152 
153 void
155 {
156  answer.resize(6, 6);
157  answer.zero();
158 
159  FloatArray u, un, n;
160  this->computeVectorOfVelocities(VM_Total, tStep, u);
161  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
162  double dudx [ 2 ] [ 2 ];
163  int w_dof_addr, u_dof_addr, d1j, d2j, dij;
164  double _u, _v, dV, rho;
165 
166  dudx [ 0 ] [ 0 ] = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
167  dudx [ 0 ] [ 1 ] = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
168  dudx [ 1 ] [ 0 ] = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
169  dudx [ 1 ] [ 1 ] = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
170 
171  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
172  dV = this->computeVolumeAround(gp);
173  rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
174  this->computeNVector(n, gp);
175 
176  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
177  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
178 
179  // dN(v)/dv
180  for ( int i = 1; i <= 2; i++ ) { // test function index
181  for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
182  for ( int j = 1; j <= 2; j++ ) { // velocity vector component
183  for ( int m = 1; m <= 3; m++ ) { // nodal components
184  w_dof_addr = ( k - 1 ) * 2 + i;
185  u_dof_addr = ( m - 1 ) * 2 + j;
186  d1j = ( j == 1 );
187  d2j = ( j == 2 );
188  dij = ( i == j );
189  answer.at(w_dof_addr, u_dof_addr) += dV * rho * n.at(k) * ( 0.0 * d1j * n.at(m) * dudx [ i - 1 ] [ 0 ] + dij * _u * b [ m - 1 ] +
190  0.0 * d2j * n.at(m) * dudx [ i - 1 ] [ 0 ] + dij * _v * c [ m - 1 ] );
191  }
192  }
193  }
194  }
195 
196  // stabilization term dN_delta/du
197  for ( int i = 1; i <= 2; i++ ) { // test function index
198  for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
199  for ( int j = 1; j <= 2; j++ ) { // velocity vector component
200  for ( int m = 1; m <= 3; m++ ) { // nodal components
201  w_dof_addr = ( k - 1 ) * 2 + i;
202  u_dof_addr = ( m - 1 ) * 2 + j;
203  //d1j = ( j == 1 );
204  //d2j = ( j == 2 );
205  dij = ( i == j );
206  answer.at(w_dof_addr, u_dof_addr) += dV * t_supg * rho *
207  ( _u * b [ k - 1 ] + _v * c [ k - 1 ] ) * ( dij * _u * b [ m - 1 ] + dij * _v * c [ m - 1 ] );
208  }
209  }
210  }
211  }
212  }
213 }
214 
215 
216 void
218 {
219  answer.resize(6);
220  answer.zero();
221  FloatArray u, eps, stress;
222  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
223  this->computeVectorOfVelocities(VM_Total, tStep, u);
224  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
225  FloatMatrix _b;
226 
227  // stabilization term K_delta
228  FloatArray un, n;
229  double _u, _v, _r;
230  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
231  // end k_delta declaration
232 
233  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
234  double dV = this->computeVolumeAround(gp);
235  this->computeBMtrx(_b, gp);
236  eps.beProductOf(_b, u);
237  mat->computeDeviatoricStressVector(stress, gp, eps, tStep);
238  answer.plusProduct(_b, stress, dV / Re);
239 
240 #if 1
241  // stabilization term k_delta
242  _r = this->computeRadiusAt(gp);
243  this->computeNVector(n, gp);
244  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
245  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
246 
247  for ( int i = 1; i <= 3; i++ ) {
248  answer.at(2 * i - 1) -= t_supg * ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( stress.at(1) / _r ) * dV / Re;
249  answer.at(2 * i) -= t_supg * ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( stress.at(4) / _r ) * dV / Re;
250  }
251 
252 #endif
253  }
254 }
255 
256 
257 void
259 {
260  answer.resize(6, 6);
261  answer.zero();
262  FloatMatrix _db, _d, _b;
263  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
264  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
265 
266  // stabilization term K_delta
267  FloatArray un, u, n, eps, stress;
268  double _u, _v, _r;
269  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
270  this->computeVectorOfVelocities(VM_Total, tStep, u);
271 
272 
273  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
274  double dV = this->computeVolumeAround(gp);
275  this->computeBMtrx(_b, gp);
276  mat->giveDeviatoricStiffnessMatrix(_d, mode, gp, tStep);
277  _db.beProductOf(_d, _b);
278  answer.plusProductUnsym(_b, _db, dV);
279  //answer.plusProductSymmUpper (_bs,_db,dV*t_supg);
280  // }
281 
282 #if 1
283  _r = this->computeRadiusAt(gp);
284  this->computeNVector(n, gp);
285  eps.beProductOf(_b, u);
286  mat->computeDeviatoricStressVector(stress, gp, eps, tStep);
287  //_mu = mat->giveCharacteristicValue(MRM_Viscosity, gp, tStep);
288 
289  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
290  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
291 
292  for ( int i = 1; i <= 3; i++ ) {
293  for ( int j = 1; j <= 6; j++ ) {
294  //answer.at(2*i-1,j) -= t_supg*(_u*b[i-1]+_v*c[i-1])*(_d.at(1,1)*_b.at(1,j)/_r - stress.at(1)/_r/_r)*dV;
295  //answer.at(2*i,j) -= t_supg*(_u*b[i-1]+_v*c[i-1])*(_d.at(4,4)*_b.at(4,j)/_r - stress.at(4)/_r/_r)*dV;
296 
297  answer.at(2 * i - 1, j) -= t_supg * ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( _db.at(1, j) / _r ) * dV;
298  answer.at(2 * i, j) -= t_supg * ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * ( _db.at(4, j) / _r ) * dV;
299  }
300  }
301 
302 #endif
303  }
304 
305 
306 
307  answer.times(1. / Re);
308 }
309 
310 
311 void
313 {
314  answer.resize(6, 3);
315  answer.zero();
316  FloatArray un, n;
317  double _u, _v;
318 
319  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
320 
321  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
322  double dV = this->computeVolumeAround(gp);
323  double _r = this->computeRadiusAt(gp);
324  this->computeNVector(n, gp);
325  // G matrix
326  for ( int i = 1; i <= 3; i++ ) {
327  for ( int j = 1; j <= 3; j++ ) {
328  answer.at(2 * i - 1, j) -= b [ i - 1 ] * n.at(j) * dV;
329  answer.at(2 * i, j) -= c [ i - 1 ] * n.at(j) * dV;
330 
331  answer.at(2 * i - 1, j) -= n.at(i) * n.at(j) * dV / _r;
332  }
333  }
334 
335  // stabilization term (G_\delta mtrx)
336  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
337  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
338  for ( int i = 1; i <= 3; i++ ) {
339  for ( int j = 1; j <= 3; j++ ) {
340  answer.at(2 * i - 1, j) += ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * b [ j - 1 ] * dV * t_supg;
341  answer.at(2 * i, j) += ( _u * b [ i - 1 ] + _v * c [ i - 1 ] ) * c [ j - 1 ] * dV * t_supg;
342  }
343  }
344  }
345 }
346 
347 void
349 {
350  answer.resize(6, 6);
351  answer.zero();
352  double n[] = {
353  b [ 0 ], c [ 0 ], b [ 1 ], c [ 1 ], b [ 2 ], c [ 2 ]
354  };
355 
356  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
357  double dV = this->computeVolumeAround(gp);
358  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
359 
360  for ( int i = 1; i <= 6; i++ ) {
361  for ( int j = 1; j <= 6; j++ ) {
362  answer.at(i, j) += dV * t_lsic * rho * n [ i - 1 ] * n [ j - 1 ];
363  }
364  }
365  }
366 }
367 
368 
369 
370 void
372 {
373  answer.resize(3, 6);
374  answer.zero();
375 
376  FloatArray n;
377  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
378  double dV = this->computeVolumeAround(gp);
379  double _r = this->computeRadiusAt(gp);
380  this->computeNVector(n, gp);
381  for ( int i = 1; i <= 3; i++ ) {
382  for ( int j = 1; j <= 3; j++ ) {
383  answer.at(j, 2 * i - 1) += b [ i - 1 ] * n.at(j) * dV;
384  answer.at(j, 2 * i) += c [ i - 1 ] * n.at(j) * dV;
385 
386  answer.at(i, 1 + ( j - 1 ) * 2) += n.at(i) * n.at(j) * dV / _r;
387  }
388  }
389  }
390 }
391 
392 void
394 {
395  // N_epsilon (due to PSPG stabilization)
396  double dudx, dudy, dvdx, dvdy, _u, _v;
397  FloatArray u, un, n;
398 
399  answer.resize(3);
400  answer.zero();
401 
402  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
403  this->computeVectorOfVelocities(VM_Total, tStep, u);
404  dudx = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
405  dudy = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
406  dvdx = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
407  dvdy = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
408 
409  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
410  double dV = this->computeVolumeAround(gp);
411  this->computeNVector(n, gp);
412 
413  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
414  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
415 
416  for ( int i = 1; i <= 3; i++ ) {
417  answer.at(i) += t_pspg * dV * ( b [ i - 1 ] * ( _u * dudx + _v * dudy ) + c [ i - 1 ] * ( _u * dvdx + _v * dvdy ) );
418  }
419  }
420 }
421 
422 
423 void
425 {
426  answer.resize(3, 6);
427  answer.zero();
428  int w_dof_addr, u_dof_addr, d1j, d2j, km1, mm1;
429  FloatArray u, un, n;
430  double _u, _v;
431 
432  this->computeVectorOfVelocities(VM_Total, tStep, u);
433  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
434 
435  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
436  double dV = this->computeVolumeAround(gp);
437  this->computeNVector(n, gp);
438 
439  _u = n.at(1) * un.at(1) + n.at(2) * un.at(3) + n.at(3) * un.at(5);
440  _v = n.at(1) * un.at(2) + n.at(2) * un.at(4) + n.at(3) * un.at(6);
441 
442  for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
443  km1 = k - 1;
444  for ( int j = 1; j <= 2; j++ ) { // velocity vector component
445  for ( int m = 1; m <= 3; m++ ) { // nodal components
446  w_dof_addr = k;
447  u_dof_addr = ( m - 1 ) * 2 + j;
448  mm1 = m - 1;
449  d1j = ( j == 1 );
450  d2j = ( j == 2 );
451  answer.at(w_dof_addr, u_dof_addr) += t_pspg * dV * ( b [ km1 ] * ( _u * d1j * b [ mm1 ] + _v * d1j * c [ mm1 ] ) + c [ km1 ] * ( _u * d2j * b [ mm1 ] + _v * d2j * c [ mm1 ] ) );
452  }
453  }
454  }
455  }
456 }
457 
459 {
460  answer.resize(3);
461  answer.zero();
462 
463 #if 1
464  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
465  FloatArray eps, stress, u;
466  FloatMatrix _b;
467  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
468 
469  // stabilization term K_eps
470  this->computeVectorOfVelocities(VM_Total, tStep, u);
471 
472  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
473  double dV = this->computeVolumeAround(gp);
474  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
475  double _r = this->computeRadiusAt(gp);
476  this->computeBMtrx(_b, gp);
477  eps.beProductOf(_b, u);
478  mat->computeDeviatoricStressVector(stress, gp, eps, tStep);
479  stress.times(1. / Re);
480  for ( int i = 1; i <= 3; i++ ) {
481  answer.at(i) -= t_pspg * ( b [ i - 1 ] * stress.at(1) + c [ i - 1 ] * stress.at(4) ) * dV / rho / _r;
482  }
483  }
484 
485 #endif
486 }
487 
489 {
490  answer.resize(3, 6);
491  answer.zero();
492 
493 #if 1
494  double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
495  FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
496  FloatMatrix _d, _b, _db;
497  //FloatArray eps, stress,u;
498 
499  // stabilization term K_eps
500  //this -> computeVectorOfVelocities(VM_Total,tStep, u) ;
501  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
502  double dV = this->computeVolumeAround(gp);
503  double rho = mat->give('d', gp);
504  double _r = this->computeRadiusAt(gp);
505  this->computeBMtrx(_b, gp);
506  mat->giveDeviatoricStiffnessMatrix(_d, TangentStiffness, gp, tStep);
507  _db.beProductOf(_d, _b);
508  //eps.beProductOf (_b, u);
509  //mat->computeDeviatoricStressVector (stress,gp,eps,tStep);
510 
511  for ( int i = 1; i <= 3; i++ ) {
512  for ( int j = 1; j <= 6; j++ ) {
513  answer.at(i, j) -= t_pspg * ( b [ i - 1 ] * ( _db.at(1, j) / _r ) + c [ i - 1 ] * ( _db.at(4, j) / _r ) ) * dV / rho;
514  }
515  }
516  }
517 
518  answer.times(1. / Re);
519 #endif
520 }
521 
522 
523 void
525 {
526  answer.resize(3, 6);
527  answer.zero();
528  FloatArray n;
529  // M_\epsilon
530 
531  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
532  double dV = this->computeVolumeAround(gp);
533  this->computeNVector(n, gp);
534 
535  for ( int i = 1; i <= 3; i++ ) {
536  for ( int j = 1; j <= 3; j++ ) {
537  answer.at(i, 2 * j - 1) += t_pspg * dV * b [ i - 1 ] * n.at(j);
538  answer.at(i, 2 * j) += t_pspg * dV * c [ i - 1 ] * n.at(j);
539  }
540  }
541  }
542 }
543 
544 
545 void
547 {
548  answer.resize(3, 3);
549  answer.zero();
550 
551  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
552  double dV = this->computeVolumeAround(gp);
553  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
554  double coeff = t_pspg / rho;
555 
556  for ( int i = 1; i <= 3; i++ ) {
557  for ( int j = 1; j <= 3; j++ ) {
558  answer.at(i, j) += dV * coeff * ( b [ i - 1 ] * b [ j - 1 ] + c [ i - 1 ] * c [ j - 1 ] );
559  }
560  }
561  }
562 }
563 
564 void
566 {
567  /* one should call material driver instead */
568  FloatArray u;
569  FloatMatrix _b;
570 
571  this->computeVectorOfVelocities(VM_Total, tStep, u);
572  this->computeBMtrx(_b, gp);
573  answer.beProductOf(_b, u);
574 }
575 
576 void
578 {
579  /* one should call material driver instead */
580  FloatArray eps;
581 
582  this->computeDeviatoricStrain(eps, gp, tStep);
583  static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeDeviatoricStressVector(answer, gp, eps, tStep);
584 }
585 
586 
587 
588 void
590 {
591  answer.resize(6);
592  answer.zero();
593 
594  int nLoads;
595  FloatArray un, gVector;
596  FloatArray nV;
597 
598  // add body load (gravity) termms
599  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
600  nLoads = this->giveBodyLoadArray()->giveSize();
601  for ( int i = 1; i <= nLoads; i++ ) {
602  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
603  bcGeomType ltype = load->giveBCGeoType();
604  if ( ltype == BodyLoadBGT && load->giveBCValType() == ForceLoadBVT ) {
605  load->computeComponentArrayAt(gVector, tStep, VM_Total);
606  if ( gVector.giveSize() ) {
607  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
608  double dV = this->computeVolumeAround(gp);
609  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
610  double coeff = rho * dV;
611  this->computeNVector(nV, gp);
612  double u = nV.at(1) * un.at(1) + nV.at(2) * un.at(3) + nV.at(3) * un.at(5);
613  double v = nV.at(1) * un.at(2) + nV.at(2) * un.at(4) + nV.at(3) * un.at(6);
614 
615  answer.at(1) += coeff * ( gVector.at(1) * ( nV.at(1) + t_supg * ( b [ 0 ] * u + c [ 0 ] * v ) ) );
616  answer.at(2) += coeff * ( gVector.at(2) * ( nV.at(1) + t_supg * ( b [ 0 ] * u + c [ 0 ] * v ) ) );
617  answer.at(3) += coeff * ( gVector.at(1) * ( nV.at(2) + t_supg * ( b [ 1 ] * u + c [ 1 ] * v ) ) );
618  answer.at(4) += coeff * ( gVector.at(2) * ( nV.at(2) + t_supg * ( b [ 1 ] * u + c [ 1 ] * v ) ) );
619  answer.at(5) += coeff * ( gVector.at(1) * ( nV.at(3) + t_supg * ( b [ 2 ] * u + c [ 2 ] * v ) ) );
620  answer.at(6) += coeff * ( gVector.at(2) * ( nV.at(3) + t_supg * ( b [ 2 ] * u + c [ 2 ] * v ) ) );
621  }
622  }
623  }
624  }
625 
626  //answer.times(rc);
627 
628  // loop over sides
629  int n1, n2, code, sid;
630  double tx, ty, l, side_r;
631  //IntArray nodecounter (3);
632  for ( int j = 1; j <= boundarySides.giveSize(); j++ ) {
633  code = boundaryCodes.at(j);
634  sid = boundarySides.at(j);
635  if ( ( code & FMElement_PrescribedTractionBC ) ) {
636  FloatArray t, coords(1);
637  int n, id;
638  // integrate tractions
639  n1 = sid;
640  n2 = ( n1 == 3 ? 1 : n1 + 1 );
641 
642  tx = giveNode(n2)->giveCoordinate(1) - giveNode(n1)->giveCoordinate(1);
643  ty = giveNode(n2)->giveCoordinate(2) - giveNode(n1)->giveCoordinate(2);
644  l = sqrt(tx * tx + ty * ty);
645  // radius at side center
646  side_r = 0.5 * ( giveNode(n2)->giveCoordinate(1) + giveNode(n1)->giveCoordinate(1) );
647 
648  // if no traction bc applied but side marked as with traction load
649  // then zero traction is assumed !!!
650 
651  // loop over boundary load array
652  int numLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
653  for ( int i = 1; i <= numLoads; i++ ) {
654  n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
655  id = boundaryLoadArray.at(i * 2);
656  if ( id != sid ) {
657  continue;
658  }
659 
660  auto load = dynamic_cast< BoundaryLoad * >( domain->giveLoad(n) );
661  if ( load ) {
662  load->computeValueAt(t, tStep, coords, VM_Total);
663 
664  // here it is assumed constant traction, one point integration only
665  // n1 (u,v)
666  answer.at(n1 * 2 - 1) += t.at(1) * side_r * l / 2.;
667  answer.at(n1 * 2) += t.at(2) * side_r * l / 2.;
668  // n2 (u,v)
669  answer.at(n2 * 2 - 1) += t.at(1) * side_r * l / 2.;
670  answer.at(n2 * 2) += t.at(2) * side_r * l / 2.;
671 
672  //answer.at(n1)+= (t.at(1)*nx + t.at(2)*ny) * side_r * l/2.;
673  //answer.at(n2)+= (t.at(1)*nx + t.at(2)*ny) * side_r * l/2.;
674  }
675  }
676  }
677  }
678 }
679 
680 
681 void
683 {
684  int nLoads;
685  FloatArray gVector;
686 
687  answer.resize(3);
688  answer.zero();
689  nLoads = this->giveBodyLoadArray()->giveSize();
690  for ( int i = 1; i <= nLoads; i++ ) {
691  Load *load = domain->giveLoad( bodyLoadArray.at(i) );
692  bcGeomType ltype = load->giveBCGeoType();
693  if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
694  load->computeComponentArrayAt(gVector, tStep, VM_Total);
695  if ( gVector.giveSize() ) {
696  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
697  double dV = this->computeVolumeAround(gp);
698  double coeff = t_pspg * dV;
699 
700  answer.at(1) += coeff * ( b [ 0 ] * gVector.at(1) + c [ 0 ] * gVector.at(2) );
701  answer.at(2) += coeff * ( b [ 1 ] * gVector.at(1) + c [ 1 ] * gVector.at(2) );
702  answer.at(3) += coeff * ( b [ 2 ] * gVector.at(1) + c [ 2 ] * gVector.at(2) );
703  }
704  }
705  }
706  }
707 }
708 
709 
710 void
712 {
713  if ( type != ExternalForcesVector ) {
714  answer.clear();
715  return;
716  }
717 
718  FloatArray un, nV;
719 
720  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
721 
722  answer.resize(9);
723 
724  if ( load->giveBCValType() == ForceLoadBVT ) {
725  FloatArray gVector;
726  load->computeComponentArrayAt(gVector, tStep, VM_Total);
727 
728  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
729  double dV = this->computeVolumeAround(gp);
730  double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity(gp);
731  double coeff;
732  this->computeNVector(nV, gp);
733  double u = nV.at(1) * un.at(1) + nV.at(2) * un.at(3) + nV.at(3) * un.at(5);
734  double v = nV.at(1) * un.at(2) + nV.at(2) * un.at(4) + nV.at(3) * un.at(6);
735 
736  coeff = rho * dV;
737  answer.at(1) += coeff * gVector.at(1) * ( nV.at(1) + t_supg * ( b [ 0 ] * u + c [ 0 ] * v ) );
738  answer.at(2) += coeff * gVector.at(2) * ( nV.at(1) + t_supg * ( b [ 0 ] * u + c [ 0 ] * v ) );
739  answer.at(4) += coeff * gVector.at(1) * ( nV.at(2) + t_supg * ( b [ 1 ] * u + c [ 1 ] * v ) );
740  answer.at(5) += coeff * gVector.at(2) * ( nV.at(2) + t_supg * ( b [ 1 ] * u + c [ 1 ] * v ) );
741  answer.at(7) += coeff * gVector.at(1) * ( nV.at(3) + t_supg * ( b [ 2 ] * u + c [ 2 ] * v ) );
742  answer.at(8) += coeff * gVector.at(2) * ( nV.at(3) + t_supg * ( b [ 2 ] * u + c [ 2 ] * v ) );
743 
744  coeff = t_pspg * dV;
745  answer.at(3) += coeff * ( b [ 0 ] * gVector.at(1) + c [ 0 ] * gVector.at(2) );
746  answer.at(6) += coeff * ( b [ 1 ] * gVector.at(1) + c [ 1 ] * gVector.at(2) );
747  answer.at(9) += coeff * ( b [ 2 ] * gVector.at(1) + c [ 2 ] * gVector.at(2) );
748  }
749  }
750 }
751 
752 
753 void
755 {
756  int nLoads;
757  FloatMatrix helpMatrix;
758 
759  answer.clear();
760 
761  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
762 
763  for ( int i = 1; i <= nLoads; i++ ) {
764  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
765  int side = boundaryLoadArray.at(i * 2);
766  Load *load = domain->giveLoad(n);
767  bcType boundarytype = load->giveType();
768  if ( boundarytype == SlipWithFriction ) {
769  this->computeSlipWithFrictionBCTerm_MB(helpMatrix, load, side, tStep);
770  } else if ( boundarytype == PenetrationWithResistance ) {
771  this->computePenetrationWithResistanceBCTerm_MB(helpMatrix, load, side, tStep);
772  } else {
773  helpMatrix.clear();
774  // OOFEM_ERROR("unsupported load type class");
775  }
776 
777  answer.add(helpMatrix);
778  }
779 }
780 
781 
782 void
784 {
785  int node1, node2;
786  double l, t1, t2, _t1, _t2;
787  double beta;
788  FloatArray nt(6), n;
789  answer.clear();
790 
791  BoundaryLoad *edgeLoad = static_cast< BoundaryLoad * >(load);
792  beta = edgeLoad->giveProperty('a', tStep);
793  node1 = side;
794  node2 = ( node1 == 3 ? 1 : node1 + 1 );
795 
796  _t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
797  _t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
798  l = sqrt(_t1 * _t1 + _t2 * _t2);
799 
800  t1 = _t1 / l;
801  t2 = _t2 / l;
802 
803  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
804  double dV = this->computeVolumeAround(gp);
805  this->computeNVector(n, gp);
806 
807  nt.at(1) = n.at(1) * t1;
808  nt.at(2) = n.at(1) * t2;
809  nt.at(3) = n.at(2) * t1;
810  nt.at(4) = n.at(2) * t2;
811  nt.at(5) = n.at(3) * t1;
812  nt.at(6) = n.at(3) * t2;
813 
814  answer.plusDyadSymmUpper(nt, beta * dV);
815  }
816  answer.symmetrized();
817 }
818 
819 void
821 {
822  int node1, node2;
823  double l, n1, n2, _t1, _t2;
824  double alpha;
825  FloatArray nt(6), n;
826  answer.clear();
827 
828  BoundaryLoad *edgeLoad = static_cast< BoundaryLoad * >(load);
829  alpha = edgeLoad->giveProperty('a', tStep);
830  node1 = side;
831  node2 = ( node1 == 3 ? 1 : node1 + 1 );
832 
833 
834  _t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
835  _t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
836  l = sqrt(_t1 * _t1 + _t2 * _t2);
837 
838  n1 = _t2 / l;
839  n2 = -_t1 / l;
840 
841 
842  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
843  double dV = this->computeVolumeAround(gp);
844  this->computeNVector(n, gp);
845 
846  nt.at(1) = n.at(1) * n1;
847  nt.at(2) = n.at(1) * n2;
848  nt.at(3) = n.at(2) * n1;
849  nt.at(4) = n.at(2) * n2;
850  nt.at(5) = n.at(3) * n1;
851  nt.at(6) = n.at(3) * n2;
852 
853  answer.plusDyadSymmUpper(nt, dV / alpha);
854  }
855  answer.symmetrized();
856 }
857 
858 void
860 {
861  int node1, node2;
862  double l, n1, n2, t1, t2, dV;
863  FloatArray nt(6), n;
864  answer.clear();
865  //beta
866  //area
867  node1 = side;
868  node2 = ( node1 == 3 ? 1 : node1 + 1 );
869 
870  t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
871  t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
872  l = sqrt(t1 * t1 + t2 * t2);
873 
874  n1 = t2 / l;
875  n2 = -t1 / l;
876 
877 
878  for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
879  dV = this->computeVolumeAround(gp);
880  this->computeNVector(n, gp);
881 
882  nt.at(1) = n.at(1) * n1;
883  nt.at(2) = n.at(1) * n2;
884  nt.at(3) = n.at(2) * n1;
885  nt.at(4) = n.at(2) * n2;
886  nt.at(5) = n.at(3) * n1;
887  nt.at(6) = n.at(3) * n2;
888 
889  answer.plusDyadUnsym(nt, n, dV);
890  }
891 
892  answer.negated();
893 }
894 
895 
896 
897 void
899 {
900  int nLoads = 0;
901  FloatMatrix helpMatrix;
902  // loop over boundary load array
903 
904  answer.clear();
905 
906  nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
907 
908  for ( int i = 1; i <= nLoads; i++ ) {
909  int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
910  int side = boundaryLoadArray.at(i * 2);
911  Load *load = domain->giveLoad(n);
912  bcType boundarytype = load->giveType();
913  if ( boundarytype == OutFlowBC ) {
914  this->computeOutFlowBCTerm_MB(helpMatrix, side, tStep);
915  } else {
916  helpMatrix.clear();
917  // OOFEM_ERROR("unsupported load type class");
918  }
919 
920  answer.add(helpMatrix);
921  }
922 }
923 
924 
925 double
927 {
928  FloatArray n;
929  this->computeNVector(n, gp);
930 
931  return n.at(1) * this->giveNode(1)->giveCoordinate(1)
932  + n.at(2) * this->giveNode(2)->giveCoordinate(1)
933  + n.at(3) * this->giveNode(3)->giveCoordinate(1);
934 }
935 
937 {
938  _b.resize(4, 6);
939  double _r = this->computeRadiusAt(gp);
940 
941 
942  _b.at(1, 1) = b [ 0 ];
943  _b.at(1, 2) = 0.;
944  _b.at(1, 3) = b [ 1 ];
945  _b.at(1, 4) = 0.;
946  _b.at(1, 5) = b [ 2 ];
947  _b.at(1, 6) = 0.;
948  _b.at(2, 1) = 1. / _r;
949  _b.at(2, 2) = 0.;
950  _b.at(2, 3) = 1. / _r;
951  _b.at(2, 4) = 0.;
952  _b.at(2, 5) = 1. / _r;
953  _b.at(2, 6) = 0.;
954  _b.at(3, 1) = 0.0;
955  _b.at(3, 2) = c [ 0 ];
956  _b.at(3, 3) = 0.0;
957  _b.at(3, 4) = c [ 1 ];
958  _b.at(3, 5) = 0.0;
959  _b.at(3, 6) = c [ 2 ];
960  _b.at(4, 1) = c [ 0 ];
961  _b.at(4, 2) = b [ 0 ];
962  _b.at(4, 3) = c [ 1 ];
963  _b.at(4, 4) = b [ 1 ];
964  _b.at(4, 5) = c [ 2 ];
965  _b.at(4, 6) = b [ 2 ];
966 }
967 
968 void
970 {
972 }
973 
974 double
976 {
977  double _r, weight, detJ;
978 
980  weight = gp->giveWeight();
981  _r = computeRadiusAt(gp);
982 
983  return detJ * weight * _r;
984 }
985 
986 void
988 {
990 
991  this->rc = ( this->giveNode(1)->giveCoordinate(1) +
992  this->giveNode(2)->giveCoordinate(1) +
993  this->giveNode(3)->giveCoordinate(1) ) / 3.0;
994  //this->rc = 1.0;
995 }
996 
997 
998 
999 void
1001 {
1002  /* UGN-Based Stabilization */
1003  double h_ugn, sum = 0.0, vnorm, t_sugn1, t_sugn2, t_sugn3, u_1, u_2;
1004  double dscale, uscale, lscale, tscale, dt;
1005  //bool zeroFlag = false;
1006  int im1;
1007  FloatArray u;
1008 
1013 
1014  this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), u);
1015  u.times(uscale);
1016  double nu;
1017 
1018  // compute averaged viscosity based on rule of mixture
1019  GaussPoint *gp;
1020  if ( integrationRulesArray [ 0 ]->giveNumberOfIntegrationPoints() ) {
1021  gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
1022  } else {
1023  gp = integrationRulesArray [ 1 ]->getIntegrationPoint(0);
1024  }
1025 
1026  nu = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->giveEffectiveViscosity( gp, tStep->givePreviousStep() );
1028 
1029  dt = tStep->giveTimeIncrement() * tscale;
1030 
1031  for ( int i = 1; i <= 3; i++ ) {
1032  im1 = i - 1;
1033  sum += fabs(u.at( ( im1 ) * 2 + 1 ) * b [ im1 ] / lscale + u.at(im1 * 2 + 2) * c [ im1 ] / lscale);
1034  }
1035 
1036  /*
1037  * u_1=(u.at(1)+u.at(3)+u.at(5))/3.0;
1038  * u_2=(u.at(2)+u.at(4)+u.at(6))/3.0;
1039  * vnorm=sqrt(u_1*u_1+u_2*u_2);
1040  */
1041  vnorm = 0.;
1042  for ( int i = 1; i <= 3; i++ ) {
1043  im1 = i - 1;
1044  u_1 = u.at( ( im1 ) * 2 + 1 );
1045  u_2 = u.at( ( im1 ) * 2 + 2 );
1046  vnorm = max( vnorm, sqrt(u_1 * u_1 + u_2 * u_2) );
1047  }
1048 
1049  if ( ( vnorm == 0.0 ) || ( sum < vnorm * 1e-10 ) ) {
1050  //t_sugn1 = inf;
1051  t_sugn2 = dt / 2.0;
1052  //t_sugn3 = inf;
1053  this->t_supg = 1. / sqrt( 1. / ( t_sugn2 * t_sugn2 ) );
1054  this->t_pspg = this->t_supg;
1055  this->t_lsic = 0.0;
1056  } else {
1057  h_ugn = 2.0 * vnorm / sum;
1058  t_sugn1 = 1. / sum;
1059  t_sugn2 = dt / 2.0;
1060  t_sugn3 = h_ugn * h_ugn / 4.0 / nu;
1061 
1062  this->t_supg = 1. / sqrt( 1. / ( t_sugn1 * t_sugn1 ) + 1. / ( t_sugn2 * t_sugn2 ) + 1. / ( t_sugn3 * t_sugn3 ) );
1063  this->t_pspg = this->t_supg;
1064 
1065  //this->t_lsic = h_ugn * vnorm * z / 2.0;
1066  this->t_lsic = 0.0;
1067  }
1068 
1069  this->t_supg *= uscale / lscale;
1070  this->t_pspg *= 1. / ( lscale * dscale );
1071  this->t_lsic *= ( dscale * uscale ) / ( lscale * lscale );
1072 
1073  //this->t_lsic = 0.0;
1074 }
1075 
1076 void
1078 {
1079  int neg = 0, pos = 0, zero = 0, si = 0;
1080  double x1, x2, x3, y1, y2, y3;
1081 
1082  answer.resize(2);
1083  for ( double ifi: fi ) {
1084  if ( ifi >= 0. ) {
1085  pos++;
1086  } else if ( ifi < 0.0 ) {
1087  neg++;
1088  } else {
1089  zero++;
1090  }
1091  }
1092 
1093  if ( neg == 0 ) { // all level set values positive
1094  answer.at(1) = 1.0;
1095  answer.at(2) = 0.0;
1096  } else if ( pos == 0 ) { // all level set values negative
1097  answer.at(1) = 0.0;
1098  answer.at(2) = 1.0;
1099  } else if ( zero == 3 ) {
1100  // ???????
1101  answer.at(1) = 1.0;
1102  answer.at(2) = 0.0;
1103  } else {
1104  // zero level set inside
1105  // find the vertex with level set sign different from other two
1106  for ( int i = 1; i <= 3; i++ ) {
1107  if ( ( pos > neg ) && ( fi.at(i) < 0.0 ) ) {
1108  si = i;
1109  break;
1110  }
1111 
1112  if ( ( pos < neg ) && ( fi.at(i) >= 0.0 ) ) {
1113  si = i;
1114  break;
1115  }
1116  }
1117 
1118  if ( si && ( ( pos + neg ) == 3 ) ) {
1119  x1 = this->giveNode(si)->giveCoordinate(1);
1120  y1 = this->giveNode(si)->giveCoordinate(2);
1121 
1122  // compute intersections
1123  int prev_node = ( si > 1 ) ? si - 1 : 3;
1124  int next_node = ( si < 3 ) ? si + 1 : 1;
1125 
1126  double t = fi.at(si) / ( fi.at(si) - fi.at(next_node) );
1127  x2 = x1 + t * ( this->giveNode(next_node)->giveCoordinate(1) - this->giveNode(si)->giveCoordinate(1) );
1128  y2 = y1 + t * ( this->giveNode(next_node)->giveCoordinate(2) - this->giveNode(si)->giveCoordinate(2) );
1129 
1130  t = fi.at(si) / ( fi.at(si) - fi.at(prev_node) );
1131  x3 = x1 + t * ( this->giveNode(prev_node)->giveCoordinate(1) - this->giveNode(si)->giveCoordinate(1) );
1132  y3 = y1 + t * ( this->giveNode(prev_node)->giveCoordinate(2) - this->giveNode(si)->giveCoordinate(2) );
1133 
1134 
1135  // compute volume associated to triangle (x1,y1; x2,y2; x3,y3)
1136  double __volume = 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 ) * ( ( x1 + x2 + x3 ) / 3. );
1137  double volume = this->area * ( ( this->giveNode(1)->giveCoordinate(1) + this->giveNode(2)->giveCoordinate(1) + this->giveNode(3)->giveCoordinate(1) ) / 3. );
1138  if ( fabs(__volume) / volume > 1.00001 ) {
1139  OOFEM_ERROR("internal consistency error");
1140  }
1141 
1142  // prevent some roundoff errors
1143  if ( fabs(__volume) > volume ) {
1144  __volume = sgn(__volume) * volume;
1145  }
1146 
1147  if ( pos > neg ) {
1148  // negative area computed
1149  answer.at(2) = fabs(__volume) / volume;
1150  answer.at(1) = 1.0 - answer.at(2);
1151  } else {
1152  // postive area computed
1153  answer.at(1) = fabs(__volume) / volume;
1154  answer.at(2) = 1.0 - answer.at(1);
1155  }
1156  } else {
1157  OOFEM_ERROR("internal consistency error");
1158  }
1159  }
1160 }
1161 } // end namespace oofem
CrossSection * giveCrossSection()
Definition: element.C:495
virtual void computePressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes pressure terms for momentum balance equations(s).
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 evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: fei2dtrlin.C:43
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...
Fluid cross-section.
virtual void computeAccelerationTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes acceleration terms for mass conservation equation.
#define FMElement_PrescribedTractionBC
Definition: fmelement.h:44
virtual void computeAdvectionTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes nonlinear advection terms for momentum balance equations(s).
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 computeDiffusionTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes diffusion terms for mass conservation equation.
const FloatArray & giveSubPatchCoordinates()
Returns local sub-patch coordinates of the receiver.
Definition: gausspoint.h:144
bcGeomType
Type representing the geometric character of loading.
Definition: bcgeomtype.h:40
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1) ...
Definition: mathfem.h:91
Class implementing element body load, acting over whole element volume (e.g., the dead weight)...
Definition: bodyload.h:49
void clear()
Clears receiver (zero size).
Definition: floatarray.h:206
Base class for fluid problems.
Definition: fluidmodel.h:46
IntArray boundarySides
Array of boundary sides.
Definition: supgelement.h:62
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
void negated()
Changes sign of receiver values.
Definition: floatmatrix.C:1605
virtual double giveProperty(int aProperty, TimeStep *tStep, const std::map< std::string, FunctionArgument > &valDict)
Definition: boundaryload.C:137
virtual void computeComponentArrayAt(FloatArray &answer, TimeStep *tStep, ValueModeType mode)
Computes boundary condition value - its components values at given time.
Definition: load.C:82
void plusDyadUnsym(const FloatArray &a, const FloatArray &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:811
virtual void computeDeviatoricStress(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
virtual double giveCoordinate(int i)
Definition: node.C:82
virtual void updateStabilizationCoeffs(TimeStep *tStep)
Updates the stabilization coefficients used for CBS and SUPG algorithms.
virtual void computeAdvectionTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes advection terms for mass conservation equation.
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
MatResponseMode
Describes the character of characteristic material matrix.
virtual void LS_PCS_computeVOFFractions(FloatArray &answer, FloatArray &fi)
Returns VOF fractions for each material on element according to nodal values of level set function (p...
Distributed body load.
Definition: bcgeomtype.h:43
virtual double computeRadiusAt(GaussPoint *gp)
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
virtual void computeBCLhsTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - velocity.
void plusProductUnsym(const FloatMatrix &a, const FloatMatrix &b, double dV)
Adds to the receiver the product .
Definition: floatmatrix.C:780
virtual double give(int aProperty, GaussPoint *gp)
Returns the value of material property &#39;aProperty&#39;.
Definition: material.C:52
bcType
Type representing the type of bc.
Definition: bctype.h:40
virtual void computeSlipWithFrictionBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
Computes Lhs term due to applied slip with friction bc.
virtual void computeAdvectionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes the derivative of advection terms for mass conservation equation with respect to nodal veloc...
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...
Abstract base class representing a boundary load (force, momentum, ...) that acts directly on a bound...
Definition: boundaryload.h:110
virtual void computeVectorOfVelocities(ValueModeType mode, TimeStep *tStep, FloatArray &velocities)
Definition: fmelement.C:48
virtual void computeOutFlowBCTerm_MB(FloatMatrix &answer, int side, TimeStep *tStep)
Computes Lhs contribution due to outflow BC.
#define OOFEM_ERROR(...)
Definition: error.h:61
REGISTER_Element(LSpace)
virtual double giveVariableScale(VarScaleType varId)
Returns the scale factor for given variable type.
Definition: engngm.h:1087
virtual void initGeometry()
Definition: tr1_2d_supg.C:1504
TR1_2D_SUPG_AXI(int n, Domain *d)
virtual double giveWeight()
Returns integration weight of receiver.
Definition: gausspoint.h:181
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
double t_supg
Stabilization coefficients, updated for each solution step in updateStabilizationCoeffs() ...
Definition: supgelement.h:70
Wrapper around element definition to provide FEICellGeometry interface.
Definition: feinterpol.h:95
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
virtual void computeBMtrx(FloatMatrix &answer, GaussPoint *gp)
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 computeBCLhsPressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - pressure.
virtual void computeAccelerationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes acceleration terms (generalized mass matrix with stabilization terms) for momentum balance e...
TimeStep * givePreviousStep()
Returns pointer to previous solution step.
Definition: timestep.C:114
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
virtual void computeDiffusionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes diffusion derivative terms for mass conservation equation.
virtual void computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
Computes the contribution of the given body load (volumetric).
virtual void computeLSICStabilizationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes SLIC stabilization term for momentum balance equation(s).
virtual int setupIntegrationPoints(IntegrationRule &irule, int npoints, Element *element)
Sets up integration rule for the given element.
Definition: crosssection.C:54
Class representing vector of real numbers.
Definition: floatarray.h:82
void plusDyadSymmUpper(const FloatArray &a, double dV)
Adds to the receiver the dyadic product .
Definition: floatmatrix.C:756
virtual void computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes Rhs terms due to boundary conditions.
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual void computePressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes pressure terms for mass conservation equation.
static FEI2dTrLin interp
Definition: tr1_2d_supg.h:74
virtual double giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the determinant of the transformation.
Definition: fei2dtrlin.C:147
virtual void computeDeviatoricStrain(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
virtual void computeAdvectionDerivativeTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes the derivative of advection terms for momentum balance equations(s) with respect to nodal ve...
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
CharType
Definition: chartype.h:87
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
void computeNVector(FloatArray &answer, GaussPoint *gp)
virtual void computeLinearAdvectionTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes the linear advection term for mass conservation equation.
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
virtual void computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)
Computes diffusion terms for momentum balance equations(s).
double rc
Radius at element center.
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual bcGeomType giveBCGeoType() const
Returns geometry character of boundary condition.
IntArray boundaryCodes
Boundary sides codes.
Definition: supgelement.h:64
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
virtual void computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)
Computes Rhs terms due to boundary conditions.
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
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Class representing 2d linear triangular element for solving incompressible fluid with SUPG solver...
Definition: tr1_2d_supg.h:68
Load is base abstract class for all loads.
Definition: load.h:61
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
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 initGeometry()
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.
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
virtual void computeValueAt(FloatArray &answer, TimeStep *tStep, const FloatArray &coords, ValueModeType mode)
Computes components values of load at given point - global coordinates (coordinates given)...
Definition: boundaryload.C:58
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
Class representing solution step.
Definition: timestep.h:80
virtual void computePenetrationWithResistanceBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
Computes Lhs contribution due to applied Penetration bc.
Class representing Gaussian-quadrature integration rule.
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
void plusProduct(const FloatMatrix &b, const FloatArray &s, double dV)
Adds the product .
Definition: floatarray.C:226
virtual void computeGaussPoints()
Initializes the array of integration rules member variable.
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:32 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011