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