OOFEM 3.0
Loading...
Searching...
No Matches
tr1_2d_supg.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.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 "intarray.h"
45#include "domain.h"
46#include "mathfem.h"
47#include "engngm.h"
48#include "fm/Materials/fluiddynamicmaterial.h"
49#include "fluidcrosssection.h"
50#include "load.h"
51#include "timestep.h"
52#include "boundaryload.h"
53#include "geotoolbox.h"
54#include "materialinterface.h"
55#include "contextioerr.h"
56#include "reinforcement.h"
57#include "crosssection.h"
58#include "dynamicinputrecord.h"
59#include "classfactory.h"
60#include "paramkey.h"
61
62#ifdef __OOFEG
63 #include "oofeggraphiccontext.h"
64 #include "connectivitytable.h"
65#endif
66
67namespace oofem {
68#define TRSUPG_ZERO_VOF 1.e-8
69
70REGISTER_Element(TR1_2D_SUPG);
71
72ParamKey TR1_2D_SUPG::IPK_TR1_2D_SUPG_vof("vof");
73ParamKey TR1_2D_SUPG::IPK_TR1_2D_SUPG_pvof("pvof");
74ParamKey TR1_2D_SUPG::IPK_TR1_2D_SUPG_mat0("mat0");
75ParamKey TR1_2D_SUPG::IPK_TR1_2D_SUPG_mat1("mat1");
76
77FEI2dTrLin TR1_2D_SUPG :: interp(1, 2);
78
79TR1_2D_SUPG :: TR1_2D_SUPG(int n, Domain *aDomain) :
80 SUPGElement(n, aDomain), SpatialLocalizerInterface(this), ZZNodalRecoveryModelInterface(this), LEPlicElementInterface()
81{
82 numberOfDofMans = 3;
83}
84
85int
86TR1_2D_SUPG :: computeNumberOfDofs()
87{
88 return 9;
89}
90
91void
92TR1_2D_SUPG :: giveDofManDofIDMask(int inode, IntArray &answer) const
93{
94 answer = {V_u, V_v, P_f};
95}
96
97FEInterpolation *
98TR1_2D_SUPG :: giveInterpolation() const { return & interp; }
99
100void
101TR1_2D_SUPG :: initializeFrom(InputRecord &ir, int priority)
102{
103 SUPGElement :: initializeFrom(ir, priority);
104
105 ParameterManager &ppm = this->giveDomain()->elementPPM;
106 bool flag;
107
108 PM_UPDATE_PARAMETER_AND_REPORT(vof, ppm, ir, this->number, IPK_TR1_2D_SUPG_pvof, priority, flag) ;
109 if (flag && (vof > 0.0)) {
110 setPermanentVolumeFraction(vof);
111 this->temp_vof = vof;
112 } else {
113 this->vof = 0.0;
114 PM_UPDATE_PARAMETER_AND_REPORT(vof, ppm, ir, this->number, IPK_TR1_2D_SUPG_vof, priority, flag) ;
115 if (flag) {
116 this->temp_vof = this->vof;
117 }
118 }
119
120 this->initGeometry();
121}
122
123
124void
125TR1_2D_SUPG :: giveInputRecord(DynamicInputRecord &input)
126{
127 SUPGElement :: giveInputRecord(input);
128 if ( this->permanentVofFlag ) {
129 input.setField(this->vof, IPK_TR1_2D_SUPG_pvof.getNameCStr());
130 } else {
131 input.setField(this->vof, IPK_TR1_2D_SUPG_vof.getNameCStr());
132 }
133}
134
135
136void
137TR1_2D_SUPG :: computeGaussPoints()
138// Sets up the array containing the four Gauss points of the receiver.
139{
140 if ( integrationRulesArray.size() == 0 ) {
141 integrationRulesArray.resize(1);
142 integrationRulesArray [ 0 ] = std::make_unique<GaussIntegrationRule>(1, this, 1, 3);
143 this->giveCrossSection()->setupIntegrationPoints(* integrationRulesArray [ 0 ], 1, this);
144 }
145}
146
147
148void
149TR1_2D_SUPG :: computeAccelerationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
150{
151 answer.resize(6, 6);
152 answer.zero();
153 FloatArray un;
154 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
155
156 double ar6 = rho * area / 6.0;
157 double ar12 = rho * area / 12.0;
158 double usum, vsum, coeff;
159
160 /* consistent mass */
161
162 answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = ar6;
163 answer.at(4, 4) = answer.at(5, 5) = answer.at(6, 6) = ar6;
164
165 answer.at(1, 3) = answer.at(1, 5) = ar12;
166 answer.at(3, 1) = answer.at(3, 5) = ar12;
167 answer.at(5, 1) = answer.at(5, 3) = ar12;
168
169 answer.at(2, 4) = answer.at(2, 6) = ar12;
170 answer.at(4, 2) = answer.at(4, 6) = ar12;
171 answer.at(6, 2) = answer.at(6, 4) = ar12;
172
173 /* SUPG stabilization term */
174 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
175
176 usum = un.at(1) + un.at(3) + un.at(5);
177 vsum = un.at(2) + un.at(4) + un.at(6);
178 coeff = rho * t_supg * area / 12.;
179
180 answer.at(1, 1) += coeff * ( b [ 0 ] * ( usum + un.at(1) ) + c [ 0 ] * ( vsum + un.at(2) ) );
181 answer.at(1, 3) += coeff * ( b [ 0 ] * ( usum + un.at(3) ) + c [ 0 ] * ( vsum + un.at(4) ) );
182 answer.at(1, 5) += coeff * ( b [ 0 ] * ( usum + un.at(5) ) + c [ 0 ] * ( vsum + un.at(6) ) );
183
184 answer.at(2, 2) += coeff * ( b [ 0 ] * ( usum + un.at(1) ) + c [ 0 ] * ( vsum + un.at(2) ) );
185 answer.at(2, 4) += coeff * ( b [ 0 ] * ( usum + un.at(3) ) + c [ 0 ] * ( vsum + un.at(4) ) );
186 answer.at(2, 6) += coeff * ( b [ 0 ] * ( usum + un.at(5) ) + c [ 0 ] * ( vsum + un.at(6) ) );
187
188 answer.at(3, 1) += coeff * ( b [ 1 ] * ( usum + un.at(1) ) + c [ 1 ] * ( vsum + un.at(2) ) );
189 answer.at(3, 3) += coeff * ( b [ 1 ] * ( usum + un.at(3) ) + c [ 1 ] * ( vsum + un.at(4) ) );
190 answer.at(3, 5) += coeff * ( b [ 1 ] * ( usum + un.at(5) ) + c [ 1 ] * ( vsum + un.at(6) ) );
191
192 answer.at(4, 2) += coeff * ( b [ 1 ] * ( usum + un.at(1) ) + c [ 1 ] * ( vsum + un.at(2) ) );
193 answer.at(4, 4) += coeff * ( b [ 1 ] * ( usum + un.at(3) ) + c [ 1 ] * ( vsum + un.at(4) ) );
194 answer.at(4, 6) += coeff * ( b [ 1 ] * ( usum + un.at(5) ) + c [ 1 ] * ( vsum + un.at(6) ) );
195
196 answer.at(5, 1) += coeff * ( b [ 2 ] * ( usum + un.at(1) ) + c [ 2 ] * ( vsum + un.at(2) ) );
197 answer.at(5, 3) += coeff * ( b [ 2 ] * ( usum + un.at(3) ) + c [ 2 ] * ( vsum + un.at(4) ) );
198 answer.at(5, 5) += coeff * ( b [ 2 ] * ( usum + un.at(5) ) + c [ 2 ] * ( vsum + un.at(6) ) );
199
200 answer.at(6, 2) += coeff * ( b [ 2 ] * ( usum + un.at(1) ) + c [ 2 ] * ( vsum + un.at(2) ) );
201 answer.at(6, 4) += coeff * ( b [ 2 ] * ( usum + un.at(3) ) + c [ 2 ] * ( vsum + un.at(4) ) );
202 answer.at(6, 6) += coeff * ( b [ 2 ] * ( usum + un.at(5) ) + c [ 2 ] * ( vsum + un.at(6) ) );
203}
204
205
206void
207TR1_2D_SUPG :: computeAdvectionTerm_MB(FloatArray &answer, TimeStep *tStep)
208{
209 answer.resize(6);
210 answer.zero();
211
212 FloatArray u, un;
213 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
214 double dudx, dudy, dvdx, dvdy, usum, vsum, coeff;
215 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
216 this->computeVectorOfVelocities(VM_Total, tStep, u);
217
218 dudx = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
219 dudy = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
220 dvdx = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
221 dvdy = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
222
223 usum = un.at(1) + un.at(3) + un.at(5);
224 vsum = un.at(2) + un.at(4) + un.at(6);
225 coeff = rho * area / 12.0;
226
227 // standard galerkin term
228 answer.at(1) = coeff * ( dudx * ( usum + un.at(1) ) + dudy * ( vsum + un.at(2) ) );
229 answer.at(3) = coeff * ( dudx * ( usum + un.at(3) ) + dudy * ( vsum + un.at(4) ) );
230 answer.at(5) = coeff * ( dudx * ( usum + un.at(5) ) + dudy * ( vsum + un.at(6) ) );
231 answer.at(2) = coeff * ( dvdx * ( usum + un.at(1) ) + dvdy * ( vsum + un.at(2) ) );
232 answer.at(4) = coeff * ( dvdx * ( usum + un.at(3) ) + dvdy * ( vsum + un.at(4) ) );
233 answer.at(6) = coeff * ( dvdx * ( usum + un.at(5) ) + dvdy * ( vsum + un.at(6) ) );
234
235 // supg stabilization term
236 coeff = t_supg * rho * area / 12.0;
237 double u1u1 = usum * usum + un.at(1) * un.at(1) + un.at(3) * un.at(3) + un.at(5) * un.at(5);
238 double u1u2 = usum * vsum + un.at(1) * un.at(2) + un.at(3) * un.at(4) + un.at(5) * un.at(6);
239 double u2u2 = vsum * vsum + un.at(2) * un.at(2) + un.at(4) * un.at(4) + un.at(6) * un.at(6);
240 answer.at(1) += coeff * ( b [ 0 ] * ( dudx * u1u1 + dudy * u1u2 ) + c [ 0 ] * ( dudx * u1u2 + dudy * u2u2 ) );
241 answer.at(3) += coeff * ( b [ 1 ] * ( dudx * u1u1 + dudy * u1u2 ) + c [ 1 ] * ( dudx * u1u2 + dudy * u2u2 ) );
242 answer.at(5) += coeff * ( b [ 2 ] * ( dudx * u1u1 + dudy * u1u2 ) + c [ 2 ] * ( dudx * u1u2 + dudy * u2u2 ) );
243
244 answer.at(2) += coeff * ( b [ 0 ] * ( dvdx * u1u1 + dvdy * u1u2 ) + c [ 0 ] * ( dvdx * u1u2 + dvdy * u2u2 ) );
245 answer.at(4) += coeff * ( b [ 1 ] * ( dvdx * u1u1 + dvdy * u1u2 ) + c [ 1 ] * ( dvdx * u1u2 + dvdy * u2u2 ) );
246 answer.at(6) += coeff * ( b [ 2 ] * ( dvdx * u1u1 + dvdy * u1u2 ) + c [ 2 ] * ( dvdx * u1u2 + dvdy * u2u2 ) );
247}
248
249
250void
251TR1_2D_SUPG :: computeAdvectionDerivativeTerm_MB(FloatMatrix &answer, TimeStep *tStep)
252{
253 answer.resize(6, 6);
254 answer.zero();
255
256 FloatArray u, un;
257 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
258 this->computeVectorOfVelocities(VM_Total, tStep, u);
259 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
260
261 double dudx [ 2 ] [ 2 ], usum [ 2 ];
262 double coeff, ar12 = area / 12.;
263 int w_dof_addr, u_dof_addr, d1j, d2j, dij;
264
265 dudx [ 0 ] [ 0 ] = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
266 dudx [ 0 ] [ 1 ] = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
267 dudx [ 1 ] [ 0 ] = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
268 dudx [ 1 ] [ 1 ] = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
269 usum [ 0 ] = un.at(1) + un.at(3) + un.at(5);
270 usum [ 1 ] = un.at(2) + un.at(4) + un.at(6);
271
272 // dN(v)/dv
273 for ( int i = 1; i <= 2; i++ ) { // test function index
274 for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
275 for ( int j = 1; j <= 2; j++ ) { // velocity vector component
276 for ( int m = 1; m <= 3; m++ ) { // nodal components
277 w_dof_addr = ( k - 1 ) * 2 + i;
278 u_dof_addr = ( m - 1 ) * 2 + j;
279 d1j = ( j == 1 );
280 d2j = ( j == 2 );
281 dij = ( i == j );
282 coeff = ( m == k ) ? area / 6. : area / 12.;
283 answer.at(w_dof_addr, u_dof_addr) = rho * ( 0.0 * d1j * dudx [ i - 1 ] [ 0 ] * coeff + dij * b [ m - 1 ] * ar12 * ( usum [ 0 ] + un.at( ( k - 1 ) * 2 + 1 ) ) +
284 0.0 * d2j * dudx [ i - 1 ] [ 1 ] * coeff + dij * c [ m - 1 ] * ar12 * ( usum [ 1 ] + un.at( ( k - 1 ) * 2 + 2 ) ) );
285 }
286 }
287 }
288 }
289
290 // stabilization term dN_delta/du
291 for ( int i = 1; i <= 2; i++ ) { // test function index
292 for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
293 for ( int j = 1; j <= 2; j++ ) { // velocity vector component
294 for ( int m = 1; m <= 3; m++ ) { // nodal components
295 w_dof_addr = ( k - 1 ) * 2 + i;
296 u_dof_addr = ( m - 1 ) * 2 + j;
297 d1j = ( j == 1 );
298 d2j = ( j == 2 );
299 dij = ( i == j );
300 answer.at(w_dof_addr, u_dof_addr) += t_supg * rho *
301 (
302 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
303 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
304 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
305 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
306
307 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
308 dij * b [ k - 1 ] * b [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 0 ] + un.at(1) * un.at(1) + un.at(3) * un.at(3) + un.at(5) * un.at(5) ) +
309 0.0 * d2j * b [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
310 dij * b [ k - 1 ] * c [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 1 ] + un.at(1) * un.at(2) + un.at(3) * un.at(4) + un.at(5) * un.at(6) ) +
311
312 0.0 * d1j * c [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
313 dij * c [ k - 1 ] * b [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 1 ] + un.at(1) * un.at(2) + un.at(3) * un.at(4) + un.at(5) * un.at(6) ) +
314 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
315 dij * c [ k - 1 ] * c [ m - 1 ] * ar12 * ( usum [ 1 ] * usum [ 1 ] + un.at(2) * un.at(2) + un.at(4) * un.at(4) + un.at(6) * un.at(6) )
316 );
317 }
318 }
319 }
320 }
321}
322
323
324void
325TR1_2D_SUPG :: computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)
326{
327 double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
328
329 FloatArray u;
330 this->computeVectorOfVelocities(VM_Total, tStep, u);
331
332 FloatArrayF<3> eps = {
333 b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5),
334 c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6),
335 b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6) + c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5),
336 };
337 auto stress = (1. / Re) * static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeDeviatoricStress2D(
338 eps, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep);
339
340 // \int dNu/dxj \Tau_ij
341 answer.resize(6);
342 for ( int i = 0; i < 3; i++ ) {
343 //rh1p(1,lok) = -geome(7,ia)*( sigxx(ia)*geome(lok,ia) + sigxy(ia)*geome(lok1,ia) )*0.5d+00;
344 //rh1p(2,lok) = -geome(7,ia)*( sigxy(ia)*geome(lok,ia) + sigyy(ia)*geome(lok1,ia) )*0.5d+00;
345 answer.at( ( i ) * 2 + 1 ) = area * ( stress.at(1) * b [ i ] + stress.at(3) * c [ i ] );
346 answer.at( ( i + 1 ) * 2 ) = area * ( stress.at(3) * b [ i ] + stress.at(2) * c [ i ] );
347 }
348}
349
350
351void
352TR1_2D_SUPG :: computeDiffusionDerivativeTerm_MB(FloatMatrix &answer, MatResponseMode mode, TimeStep *tStep)
353{
354 FloatMatrix _db, _d, _b(3, 6);
355 double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
356
357 _b.at(1, 1) = b [ 0 ];
358 _b.at(1, 2) = 0.;
359 _b.at(1, 3) = b [ 1 ];
360 _b.at(1, 4) = 0.;
361 _b.at(1, 5) = b [ 2 ];
362 _b.at(1, 6) = 0.;
363 _b.at(2, 1) = 0.;
364 _b.at(2, 2) = c [ 0 ];
365 _b.at(2, 3) = 0.;
366 _b.at(2, 4) = c [ 1 ];
367 _b.at(2, 5) = 0.;
368 _b.at(2, 6) = c [ 2 ];
369 _b.at(3, 1) = c [ 0 ];
370 _b.at(3, 2) = b [ 0 ];
371 _b.at(3, 3) = c [ 1 ];
372 _b.at(3, 4) = b [ 1 ];
373 _b.at(3, 5) = c [ 2 ];
374 _b.at(3, 6) = b [ 2 ];
375
376 _d = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeTangent2D(
377 mode, integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep);
378 _db.beProductOf(_d, _b);
379 answer.resize(6, 6);
380 answer.zero();
381 answer.plusProductUnsym(_b, _db, area); //answer.plusProduct (_b,_db,area);
382 answer.times(1. / Re);
383}
384
385void
386TR1_2D_SUPG :: computePressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
387{
388 answer.resize(6, 3);
389 answer.zero();
390 FloatArray p, un;
391 double usum, vsum;
392 double ar3 = area / 3.0, coeff;
393
394 this->computeVectorOfPressures(VM_Total, tStep, p);
395
396 // G matrix
397 answer.at(1, 1) = answer.at(1, 2) = answer.at(1, 3) = -b [ 0 ] * ar3;
398 answer.at(3, 1) = answer.at(3, 2) = answer.at(3, 3) = -b [ 1 ] * ar3;
399 answer.at(5, 1) = answer.at(5, 2) = answer.at(5, 3) = -b [ 2 ] * ar3;
400
401 answer.at(2, 1) = answer.at(2, 2) = answer.at(2, 3) = -c [ 0 ] * ar3;
402 answer.at(4, 1) = answer.at(4, 2) = answer.at(4, 3) = -c [ 1 ] * ar3;
403 answer.at(6, 1) = answer.at(6, 2) = answer.at(6, 3) = -c [ 2 ] * ar3;
404
405 // stabilization term (G_\delta mtrx)
406 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
407 usum = un.at(1) + un.at(3) + un.at(5);
408 vsum = un.at(2) + un.at(4) + un.at(6);
409 coeff = ar3 * t_supg;
410
411 answer.at(1, 1) += coeff * ( usum * b [ 0 ] * b [ 0 ] + vsum * c [ 0 ] * b [ 0 ] );
412 answer.at(1, 2) += coeff * ( usum * b [ 0 ] * b [ 1 ] + vsum * c [ 0 ] * b [ 1 ] );
413 answer.at(1, 3) += coeff * ( usum * b [ 0 ] * b [ 2 ] + vsum * c [ 0 ] * b [ 2 ] );
414
415 answer.at(3, 1) += coeff * ( usum * b [ 1 ] * b [ 0 ] + vsum * c [ 1 ] * b [ 0 ] );
416 answer.at(3, 2) += coeff * ( usum * b [ 1 ] * b [ 1 ] + vsum * c [ 1 ] * b [ 1 ] );
417 answer.at(3, 3) += coeff * ( usum * b [ 1 ] * b [ 2 ] + vsum * c [ 1 ] * b [ 2 ] );
418
419 answer.at(5, 1) += coeff * ( usum * b [ 2 ] * b [ 0 ] + vsum * c [ 2 ] * b [ 0 ] );
420 answer.at(5, 2) += coeff * ( usum * b [ 2 ] * b [ 1 ] + vsum * c [ 2 ] * b [ 1 ] );
421 answer.at(5, 3) += coeff * ( usum * b [ 2 ] * b [ 2 ] + vsum * c [ 2 ] * b [ 2 ] );
422
423 answer.at(2, 1) += coeff * ( usum * b [ 0 ] * c [ 0 ] + vsum * c [ 0 ] * c [ 0 ] );
424 answer.at(2, 2) += coeff * ( usum * b [ 0 ] * c [ 1 ] + vsum * c [ 0 ] * c [ 1 ] );
425 answer.at(2, 3) += coeff * ( usum * b [ 0 ] * c [ 2 ] + vsum * c [ 0 ] * c [ 2 ] );
426
427 answer.at(4, 1) += coeff * ( usum * b [ 1 ] * c [ 0 ] + vsum * c [ 1 ] * c [ 0 ] );
428 answer.at(4, 2) += coeff * ( usum * b [ 1 ] * c [ 1 ] + vsum * c [ 1 ] * c [ 1 ] );
429 answer.at(4, 3) += coeff * ( usum * b [ 1 ] * c [ 2 ] + vsum * c [ 1 ] * c [ 2 ] );
430
431 answer.at(6, 1) += coeff * ( usum * b [ 2 ] * c [ 0 ] + vsum * c [ 2 ] * c [ 0 ] );
432 answer.at(6, 2) += coeff * ( usum * b [ 2 ] * c [ 1 ] + vsum * c [ 2 ] * c [ 1 ] );
433 answer.at(6, 3) += coeff * ( usum * b [ 2 ] * c [ 2 ] + vsum * c [ 2 ] * c [ 2 ] );
434}
435
436
437void
438TR1_2D_SUPG :: computeLSICStabilizationTerm_MB(FloatMatrix &answer, TimeStep *tStep)
439{
440 answer.resize(6, 6);
441 answer.zero();
442 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
443 double coeff = area * t_lsic * rho;
444 double n[] = {
445 b [ 0 ], c [ 0 ], b [ 1 ], c [ 1 ], b [ 2 ], c [ 2 ]
446 };
447
448 for ( int i = 1; i <= 6; i++ ) {
449 for ( int j = 1; j <= 6; j++ ) {
450 answer.at(i, j) = coeff * n [ i - 1 ] * n [ j - 1 ];
451 }
452 }
453}
454
455
456void
457TR1_2D_SUPG :: computeLinearAdvectionTerm_MC(FloatMatrix &answer, TimeStep *tStep)
458{
459 answer.resize(3, 6);
460 answer.zero();
461 double ar3 = area / 3.0;
462
463 // G^T matrix
464 answer.at(1, 1) = answer.at(2, 1) = answer.at(3, 1) = b [ 0 ] * ar3;
465 answer.at(1, 3) = answer.at(2, 3) = answer.at(3, 3) = b [ 1 ] * ar3;
466 answer.at(1, 5) = answer.at(2, 5) = answer.at(3, 5) = b [ 2 ] * ar3;
467
468 answer.at(1, 2) = answer.at(2, 2) = answer.at(3, 2) = c [ 0 ] * ar3;
469 answer.at(1, 4) = answer.at(2, 4) = answer.at(3, 4) = c [ 1 ] * ar3;
470 answer.at(1, 6) = answer.at(2, 6) = answer.at(3, 6) = c [ 2 ] * ar3;
471}
472
473void
474TR1_2D_SUPG :: computeAdvectionTerm_MC(FloatArray &answer, TimeStep *tStep)
475{
476 // N_epsilon (due to PSPG stabilization)
477 double coeff = t_pspg * area / 3.0;
478 double dudx, dudy, dvdx, dvdy, usum, vsum;
479 FloatArray u, un;
480
481 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
482 this->computeVectorOfVelocities(VM_Total, tStep, u);
483
484 dudx = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
485 dudy = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
486 dvdx = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
487 dvdy = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
488
489 usum = un.at(1) + un.at(3) + un.at(5);
490 vsum = un.at(2) + un.at(4) + un.at(6);
491
492 answer.resize(3);
493
494 answer.at(1) = coeff * ( b [ 0 ] * ( dudx * usum + dudy * vsum ) + c [ 0 ] * ( dvdx * usum + dvdy * vsum ) );
495 answer.at(2) = coeff * ( b [ 1 ] * ( dudx * usum + dudy * vsum ) + c [ 1 ] * ( dvdx * usum + dvdy * vsum ) );
496 answer.at(3) = coeff * ( b [ 2 ] * ( dudx * usum + dudy * vsum ) + c [ 2 ] * ( dvdx * usum + dvdy * vsum ) );
497}
498
499
500void
501TR1_2D_SUPG :: computeAdvectionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)
502{
503 answer.resize(3, 6);
504 answer.zero();
505 int w_dof_addr, u_dof_addr, d1j, d2j, km1, mm1;
506 FloatArray u, un;
507
508 this->computeVectorOfVelocities(VM_Total, tStep, u);
509 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
510
511 double dudx [ 2 ] [ 2 ], usum [ 2 ];
512 double coeff;
513
514 dudx [ 0 ] [ 0 ] = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
515 dudx [ 0 ] [ 1 ] = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
516 dudx [ 1 ] [ 0 ] = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
517 dudx [ 1 ] [ 1 ] = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
518 usum [ 0 ] = un.at(1) + un.at(3) + un.at(5);
519 usum [ 1 ] = un.at(2) + un.at(4) + un.at(6);
520
521 // dN_epsilon(v)/dv
522 coeff = t_pspg * area / 3.;
523 for ( int k = 1; k <= 3; k++ ) { // nodal val of function w
524 km1 = k - 1;
525 for ( int j = 1; j <= 2; j++ ) { // velocity vector component
526 for ( int m = 1; m <= 3; m++ ) { // nodal components
527 w_dof_addr = k;
528 u_dof_addr = ( m - 1 ) * 2 + j;
529 mm1 = m - 1;
530 d1j = ( j == 1 );
531 d2j = ( j == 2 );
532 answer.at(w_dof_addr, u_dof_addr) = coeff * ( 0.0 * d1j * b [ km1 ] * dudx [ 0 ] [ 0 ] + d1j * b [ km1 ] * b [ mm1 ] * usum [ 0 ] +
533 0.0 * d2j * b [ km1 ] * dudx [ 0 ] [ 1 ] + d1j * b [ km1 ] * c [ mm1 ] * usum [ 1 ] +
534 0.0 * d1j * c [ km1 ] * dudx [ 1 ] [ 0 ] + d2j * c [ km1 ] * b [ mm1 ] * usum [ 0 ] +
535 0.0 * d2j * c [ km1 ] * dudx [ 1 ] [ 1 ] + d2j * c [ km1 ] * c [ mm1 ] * usum [ 1 ] );
536 }
537 }
538 }
539}
540
541void
542TR1_2D_SUPG :: computeAccelerationTerm_MC(FloatMatrix &answer, TimeStep *tStep)
543{
544 answer.resize(3, 6);
545 answer.zero();
546 double coeff = t_pspg * area / 3.0;
547 // M_\epsilon
548
549 answer.at(1, 1) = answer.at(1, 3) = answer.at(1, 5) = coeff * b [ 0 ];
550 answer.at(1, 2) = answer.at(1, 4) = answer.at(1, 6) = coeff * c [ 0 ];
551 answer.at(2, 1) = answer.at(2, 3) = answer.at(2, 5) = coeff * b [ 1 ];
552 answer.at(2, 2) = answer.at(2, 4) = answer.at(2, 6) = coeff * c [ 1 ];
553 answer.at(3, 1) = answer.at(3, 3) = answer.at(3, 5) = coeff * b [ 2 ];
554 answer.at(3, 2) = answer.at(3, 4) = answer.at(3, 6) = coeff * c [ 2 ];
555}
556
557void
558TR1_2D_SUPG :: computePressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)
559{
560 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
561 double coeff = t_pspg * area / rho;
562 answer.resize(3, 3);
563
564
565 for ( int i = 1; i <= 3; i++ ) {
566 for ( int j = 1; j <= 3; j++ ) {
567 answer.at(i, j) = coeff * ( b [ i - 1 ] * b [ j - 1 ] + c [ i - 1 ] * c [ j - 1 ] );
568 }
569 }
570}
571
572
573void
574TR1_2D_SUPG :: computeSlipWithFrictionBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
575{
576 int node1, node2, node3;
577 int d1 = 1;
578 int d2 = 1;
579 int d3 = 1;
580 double l, t1, t2, _t1, _t2;
581 double beta;
582
583 answer.resize(6, 6);
584 answer.zero();
585 //beta
586 //area
587 BoundaryLoad *edgeLoad = static_cast< BoundaryLoad * >(load);
588 beta = edgeLoad->giveProperty('a', tStep);
589 node1 = side;
590 node2 = ( node1 == 3 ? 1 : node1 + 1 );
591
592 node3 = ( node2 == 3 ? 1 : node2 + 1 );
593
594 switch ( node3 ) {
595 case 1:
596 d1 = 0;
597 case 2:
598 d2 = 0;
599 case 3:
600 d3 = 0;
601 }
602
603 _t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
604 _t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
605 l = sqrt(_t1 * _t1 + _t2 * _t2);
606
607 t1 = _t1 / l;
608 t2 = _t2 / l;
609
610 double t11 = t1 * t1;
611 double t12 = t1 * t2;
612 double t22 = t2 * t2;
613
614 double d12 = d1 * d2;
615 double d13 = d1 * d3;
616 double d23 = d2 * d3;
617
618 answer.at(1, 1) = ( l / 3 ) * beta * d1 * t11;
619 answer.at(1, 2) = ( l / 3 ) * beta * d1 * t12;
620 answer.at(2, 1) = ( l / 3 ) * beta * d1 * t12;
621 answer.at(2, 2) = ( l / 3 ) * beta * d1 * t22;
622
623 answer.at(1, 3) = ( l / 6 ) * beta * d12 * t11;
624 answer.at(1, 4) = ( l / 6 ) * beta * d12 * t12;
625 answer.at(2, 3) = ( l / 6 ) * beta * d12 * t12;
626 answer.at(2, 4) = ( l / 6 ) * beta * d12 * t22;
627
628 answer.at(1, 5) = ( l / 6 ) * beta * d13 * t11;
629 answer.at(1, 6) = ( l / 6 ) * beta * d13 * t12;
630 answer.at(2, 5) = ( l / 6 ) * beta * d13 * t12;
631 answer.at(2, 6) = ( l / 6 ) * beta * d13 * t22;
632
633 answer.at(3, 1) = ( l / 6 ) * beta * d12 * t11;
634 answer.at(3, 2) = ( l / 6 ) * beta * d12 * t12;
635 answer.at(4, 1) = ( l / 6 ) * beta * d12 * t12;
636 answer.at(4, 2) = ( l / 6 ) * beta * d12 * t22;
637
638 answer.at(3, 3) = ( l / 3 ) * beta * d2 * t11;
639 answer.at(3, 4) = ( l / 3 ) * beta * d2 * t12;
640 answer.at(4, 3) = ( l / 3 ) * beta * d2 * t12;
641 answer.at(4, 4) = ( l / 3 ) * beta * d2 * t22;
642
643 answer.at(3, 5) = ( l / 6 ) * beta * d23 * t11;
644 answer.at(3, 6) = ( l / 6 ) * beta * d23 * t12;
645 answer.at(4, 5) = ( l / 6 ) * beta * d23 * t12;
646 answer.at(4, 6) = ( l / 6 ) * beta * d23 * t22;
647
648 answer.at(5, 1) = ( l / 6 ) * beta * d13 * t11;
649 answer.at(5, 2) = ( l / 6 ) * beta * d13 * t12;
650 answer.at(6, 1) = ( l / 6 ) * beta * d13 * t12;
651 answer.at(6, 2) = ( l / 6 ) * beta * d13 * t22;
652
653 answer.at(5, 3) = ( l / 6 ) * beta * d23 * t11;
654 answer.at(5, 4) = ( l / 6 ) * beta * d23 * t12;
655 answer.at(6, 3) = ( l / 6 ) * beta * d23 * t12;
656 answer.at(6, 4) = ( l / 6 ) * beta * d23 * t22;
657
658 answer.at(5, 5) = ( l / 3 ) * beta * d3 * t11;
659 answer.at(5, 6) = ( l / 3 ) * beta * d3 * t12;
660 answer.at(6, 5) = ( l / 3 ) * beta * d3 * t12;
661 answer.at(6, 6) = ( l / 3 ) * beta * d3 * t22;
662
663 //answer.negated();
664}
665
666
667void
668TR1_2D_SUPG :: computePenetrationWithResistanceBCTerm_MB(FloatMatrix &answer, Load *load, int side, TimeStep *tStep)
669{
670 int node1, node2, node3;
671 int d1 = 1;
672 int d2 = 1;
673 int d3 = 1;
674 double l, n1, n2, _t1, _t2;
675 double alpha;
676
677 answer.resize(6, 6);
678 answer.zero();
679
680 BoundaryLoad *edgeLoad = static_cast< BoundaryLoad * >(load);
681 alpha = edgeLoad->giveProperty('a', tStep);
682 node1 = side;
683 node2 = ( node1 == 3 ? 1 : node1 + 1 );
684
685 node3 = ( node2 == 3 ? 1 : node2 + 1 );
686
687 switch ( node3 ) {
688 case 1:
689 d1 = 0;
690 case 2:
691 d2 = 0;
692 case 3:
693 d3 = 0;
694 }
695
696 _t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
697 _t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
698 l = sqrt(_t1 * _t1 + _t2 * _t2);
699
700 n1 = _t2 / l;
701 n2 = -_t1 / l;
702
703 double n11 = n1 * n1;
704 double n12 = n1 * n2;
705 double n22 = n2 * n2;
706
707 double d12 = d1 * d2;
708 double d13 = d1 * d3;
709 double d23 = d2 * d3;
710
711 answer.at(1, 1) = ( l / 3 ) * ( 1 / alpha ) * d1 * n11;
712 answer.at(1, 2) = ( l / 3 ) * ( 1 / alpha ) * d1 * n12;
713 answer.at(2, 1) = ( l / 3 ) * ( 1 / alpha ) * d1 * n12;
714 answer.at(2, 2) = ( l / 3 ) * ( 1 / alpha ) * d1 * n22;
715
716 answer.at(1, 3) = ( l / 6 ) * ( 1 / alpha ) * d12 * n11;
717 answer.at(1, 4) = ( l / 6 ) * ( 1 / alpha ) * d12 * n12;
718 answer.at(2, 3) = ( l / 6 ) * ( 1 / alpha ) * d12 * n12;
719 answer.at(2, 4) = ( l / 6 ) * ( 1 / alpha ) * d12 * n22;
720
721 answer.at(1, 5) = ( l / 6 ) * ( 1 / alpha ) * d13 * n11;
722 answer.at(1, 6) = ( l / 6 ) * ( 1 / alpha ) * d13 * n12;
723 answer.at(2, 5) = ( l / 6 ) * ( 1 / alpha ) * d13 * n12;
724 answer.at(2, 6) = ( l / 6 ) * ( 1 / alpha ) * d13 * n22;
725
726 answer.at(3, 1) = ( l / 6 ) * ( 1 / alpha ) * d12 * n11;
727 answer.at(3, 2) = ( l / 6 ) * ( 1 / alpha ) * d12 * n12;
728 answer.at(4, 1) = ( l / 6 ) * ( 1 / alpha ) * d12 * n12;
729 answer.at(4, 2) = ( l / 6 ) * ( 1 / alpha ) * d12 * n22;
730
731 answer.at(3, 3) = ( l / 3 ) * ( 1 / alpha ) * d2 * n11;
732 answer.at(3, 4) = ( l / 3 ) * ( 1 / alpha ) * d2 * n12;
733 answer.at(4, 3) = ( l / 3 ) * ( 1 / alpha ) * d2 * n12;
734 answer.at(4, 4) = ( l / 3 ) * ( 1 / alpha ) * d2 * n22;
735
736 answer.at(3, 5) = ( l / 6 ) * ( 1 / alpha ) * d23 * n11;
737 answer.at(3, 6) = ( l / 6 ) * ( 1 / alpha ) * d23 * n12;
738 answer.at(4, 5) = ( l / 6 ) * ( 1 / alpha ) * d23 * n12;
739 answer.at(4, 6) = ( l / 6 ) * ( 1 / alpha ) * d23 * n22;
740
741 answer.at(5, 1) = ( l / 6 ) * ( 1 / alpha ) * d13 * n11;
742 answer.at(5, 2) = ( l / 6 ) * ( 1 / alpha ) * d13 * n12;
743 answer.at(6, 1) = ( l / 6 ) * ( 1 / alpha ) * d13 * n12;
744 answer.at(6, 2) = ( l / 6 ) * ( 1 / alpha ) * d13 * n22;
745
746 answer.at(5, 3) = ( l / 6 ) * ( 1 / alpha ) * d23 * n11;
747 answer.at(5, 4) = ( l / 6 ) * ( 1 / alpha ) * d23 * n12;
748 answer.at(6, 3) = ( l / 6 ) * ( 1 / alpha ) * d23 * n12;
749 answer.at(6, 4) = ( l / 6 ) * ( 1 / alpha ) * d23 * n22;
750
751 answer.at(5, 5) = ( l / 3 ) * ( 1 / alpha ) * d3 * n11;
752 answer.at(5, 6) = ( l / 3 ) * ( 1 / alpha ) * d3 * n12;
753 answer.at(6, 5) = ( l / 3 ) * ( 1 / alpha ) * d3 * n12;
754 answer.at(6, 6) = ( l / 3 ) * ( 1 / alpha ) * d3 * n22;
755
756 //answer.negated();
757}
758
759void
760TR1_2D_SUPG :: computeOutFlowBCTerm_MB(FloatMatrix &answer, int side, TimeStep *tStep)
761{
762 int node1, node2, node3;
763 int d1 = 1;
764 int d2 = 1;
765 int d3 = 1;
766 double l, n1, n2, t1, t2;
767
768 answer.resize(6, 3);
769 answer.zero();
770 //beta
771 //area
772 node1 = side;
773 node2 = ( node1 == 3 ? 1 : node1 + 1 );
774
775 node3 = ( node2 == 3 ? 1 : node2 + 1 );
776
777 switch ( node3 ) {
778 case 1:
779 d1 = 0;
780 case 2:
781 d2 = 0;
782 case 3:
783 d3 = 0;
784 }
785
786
787 t1 = giveNode(node2)->giveCoordinate(1) - giveNode(node1)->giveCoordinate(1);
788 t2 = giveNode(node2)->giveCoordinate(2) - giveNode(node1)->giveCoordinate(2);
789 l = sqrt(t1 * t1 + t2 * t2);
790
791 n1 = t2 / l;
792 n2 = -t1 / l;
793
794 answer.at(1, 1) = ( l / 3 ) * d1 * d1 * n1;
795 answer.at(1, 2) = ( l / 6 ) * d1 * d2 * n1;
796 answer.at(1, 3) = ( l / 6 ) * d1 * d3 * n1;
797 answer.at(2, 1) = ( l / 3 ) * d1 * d1 * n2;
798 answer.at(2, 2) = ( l / 6 ) * d1 * d2 * n2;
799 answer.at(2, 3) = ( l / 6 ) * d1 * d3 * n2;
800
801 answer.at(3, 1) = ( l / 6 ) * d2 * d1 * n1;
802 answer.at(3, 2) = ( l / 3 ) * d2 * d2 * n1;
803 answer.at(3, 3) = ( l / 6 ) * d2 * d3 * n1;
804 answer.at(4, 1) = ( l / 6 ) * d2 * d1 * n2;
805 answer.at(4, 2) = ( l / 3 ) * d2 * d2 * n2;
806 answer.at(4, 3) = ( l / 6 ) * d2 * d3 * n2;
807
808 answer.at(5, 1) = ( l / 6 ) * d3 * d1 * n1;
809 answer.at(5, 2) = ( l / 6 ) * d3 * d2 * n1;
810 answer.at(5, 3) = ( l / 3 ) * d3 * d3 * n1;
811 answer.at(6, 1) = ( l / 6 ) * d3 * d1 * n2;
812 answer.at(6, 2) = ( l / 6 ) * d3 * d2 * n2;
813 answer.at(6, 3) = ( l / 3 ) * d3 * d3 * n2;
814
815 answer.negated();
816}
817
818void
819TR1_2D_SUPG :: computeHomogenizedReinforceTerm_MB(FloatMatrix &answer, Load *load, TimeStep *tStep)
820{
821 double coeffx, coeffy, usum [ 2 ];
822 FloatArray un;
823
824 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
825
826
827 usum [ 0 ] = un.at(1) + un.at(3) + un.at(5);
828 usum [ 1 ] = un.at(2) + un.at(4) + un.at(6);
829 Reinforcement *reinfload = dynamic_cast< Reinforcement * >(load);
830 double kx = reinfload->givePermeability()->at(1);
831 double ky = reinfload->givePermeability()->at(2);
832 double mu_0 = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->
833 giveEffectiveViscosity( integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep );
834 coeffx = area * mu_0 / ( 12.0 * kx );
835 coeffy = area * mu_0 / ( 12.0 * ky );
836 for ( int i = 1; i <= 3; i++ ) {
837 answer.at(2 * i - 1, 2 * i - 1) -= coeffx;
838 answer.at(2 * i, 2 * i) -= coeffy;
839 for ( int j = 1; j <= 3; j++ ) {
840 answer.at(2 * i - 1, 2 * j - 1) -= coeffx * ( 1.0 + 4.0 * t_supg * ( b [ i - 1 ] * usum [ 0 ] + c [ i - 1 ] * usum [ 1 ] ) );
841 answer.at(2 * i, 2 * j) -= coeffy * ( 1.0 + 4.0 * t_supg * ( b [ i - 1 ] * usum [ 0 ] + c [ i - 1 ] * usum [ 1 ] ) ); //pozor na b[i], c[i]!!!
842 }
843 }
844}
845
846void
847TR1_2D_SUPG :: computeHomogenizedReinforceTerm_MC(FloatMatrix &answer, Load *load, TimeStep *tStep)
848{
849 double coeffx, coeffy;
850
851 Reinforcement *reinfload = dynamic_cast< Reinforcement * >(load);
852 double kx = reinfload->givePermeability()->at(1);
853 double ky = reinfload->givePermeability()->at(2);
854 double mu_0 = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->
855 giveEffectiveViscosity( integrationRulesArray [ 0 ]->getIntegrationPoint(0), tStep );
856 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
857 coeffx = area * mu_0 / ( 3.0 * kx * rho );
858 coeffy = area * mu_0 / ( 3.0 * ky * rho );
859 for ( int i = 1; i <= 3; i++ ) {
860 for ( int j = 1; j <= 3; j++ ) {
861 answer.at(i, 2 * j - 1) -= coeffx * t_pspg * b [ i - 1 ];
862 answer.at(i, 2 * j) -= coeffy * t_pspg * c [ i - 1 ];
863 }
864 }
865}
866
867void
868TR1_2D_SUPG :: computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
869{
870 answer.resize(6);
871 answer.zero();
872
873 double usum [ 2 ];
874 bcGeomType ltype;
875 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
876 FloatArray un, gVector;
877
878 // add body load (gravity) termms
879 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
880
881
882 usum [ 0 ] = un.at(1) + un.at(3) + un.at(5);
883 usum [ 1 ] = un.at(2) + un.at(4) + un.at(6);
884 double coeff = rho * area / 3.0;
885 int nLoads = this->giveBodyLoadArray()->giveSize();
886 for ( int i = 1; i <= nLoads; i++ ) {
887 Load *load = domain->giveLoad( bodyLoadArray.at(i) );
888 ltype = load->giveBCGeoType();
889 if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
890 load->computeComponentArrayAt(gVector, tStep, VM_Total);
891 if ( gVector.giveSize() ) {
892 answer.at(1) += coeff * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 0 ] * usum [ 0 ] + c [ 0 ] * usum [ 1 ] ) ) );
893 answer.at(2) += coeff * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 0 ] * usum [ 0 ] + c [ 0 ] * usum [ 1 ] ) ) );
894 answer.at(3) += coeff * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 1 ] * usum [ 0 ] + c [ 1 ] * usum [ 1 ] ) ) );
895 answer.at(4) += coeff * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 1 ] * usum [ 0 ] + c [ 1 ] * usum [ 1 ] ) ) );
896 answer.at(5) += coeff * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 2 ] * usum [ 0 ] + c [ 2 ] * usum [ 1 ] ) ) );
897 answer.at(6) += coeff * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 2 ] * usum [ 0 ] + c [ 2 ] * usum [ 1 ] ) ) );
898 }
899 }
900
901 if ( ltype == BodyLoadBGT && load->giveBCValType() == ReinforceBVT ) {
902 Reinforcement *rload = dynamic_cast< Reinforcement * >( domain->giveLoad( bodyLoadArray.at(i) ) );
903 double phi = rload->givePorosity();
904 double alpha = rload->giveshapefactor();
905 double kx = rload->givePermeability()->at(1);
906 double ky = rload->givePermeability()->at(2);
907 double tau_0 = static_cast< FluidCrossSection * >( this->giveCrossSection() )->
908 giveFluidMaterial()->give( YieldStress, integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
909 gVector.resize(2);
910 gVector.at(1) = tau_0 * sqrt(kx * phi) / ( kx * alpha );
911 gVector.at(2) = tau_0 * sqrt(ky * phi) / ( ky * alpha );
912
913 answer.at(1) -= area * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 0 ] * usum [ 0 ] + c [ 0 ] * usum [ 1 ] ) ) );
914 answer.at(2) -= area * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 0 ] * usum [ 0 ] + c [ 0 ] * usum [ 1 ] ) ) );
915 answer.at(3) -= area * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 1 ] * usum [ 0 ] + c [ 1 ] * usum [ 1 ] ) ) );
916 answer.at(4) -= area * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 1 ] * usum [ 0 ] + c [ 1 ] * usum [ 1 ] ) ) );
917 answer.at(5) -= area * ( gVector.at(1) * ( 1.0 + t_supg * ( b [ 2 ] * usum [ 0 ] + c [ 2 ] * usum [ 1 ] ) ) );
918 answer.at(6) -= area * ( gVector.at(2) * ( 1.0 + t_supg * ( b [ 2 ] * usum [ 0 ] + c [ 2 ] * usum [ 1 ] ) ) );
919 }
920 }
921
922 // loop over sides
923 int n1, n2;
924 double tx, ty, l;
925
926 if ( true ) {
927 FloatArray t, coords(1);
928 // integrate tractions
929
930 // if no traction bc applied but side marked as with traction load
931 // then zero traction is assumed !!!
932
933 // loop over boundary load array
934 int numLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
935 for ( int i = 1; i <= numLoads; i++ ) {
936 int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
937 int id = boundaryLoadArray.at(i * 2);
938
939 n1 = id;
940 n2 = ( n1 == 3 ? 1 : n1 + 1 );
941
942 tx = giveNode(n2)->giveCoordinate(1) - giveNode(n1)->giveCoordinate(1);
943 ty = giveNode(n2)->giveCoordinate(2) - giveNode(n1)->giveCoordinate(2);
944 l = sqrt(tx * tx + ty * ty);
945
946 auto load = dynamic_cast< BoundaryLoad * >( domain->giveLoad(n) );
947 auto loadtype = load->giveType();
948 if ( loadtype == TransmissionBC ) {
949 load->computeValueAt(t, tStep, coords, VM_Total);
950
951 // here it is assumed constant traction, one point integration only
952 // n1 (u,v)
953 answer.at( ( n1 - 1 ) * 2 + 1 ) += t.at(1) * l / 2.;
954 answer.at(n1 * 2) += t.at(2) * l / 2.;
955 // n2 (u,v)
956 answer.at( ( n2 - 1 ) * 2 + 1 ) += t.at(1) * l / 2.;
957 answer.at(n2 * 2) += t.at(2) * l / 2.;
958
959 //answer.at(n1)+= (t.at(1)*nx + t.at(2)*ny) * l/2.;
960 //answer.at(n2)+= (t.at(1)*nx + t.at(2)*ny) * l/2.;
961 }
962 }
963 }
964}
965
966void
967TR1_2D_SUPG :: computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)
968{
969 int nLoads;
970 double coeff;
971 FloatArray gVector;
972
973 answer.resize(3);
974 answer.zero();
975 coeff = t_pspg * area;
976 nLoads = this->giveBodyLoadArray()->giveSize();
977 for ( int i = 1; i <= nLoads; i++ ) {
978 Load *load = domain->giveLoad( bodyLoadArray.at(i) );
979 bcGeomType ltype = load->giveBCGeoType();
980 if ( ltype == BodyLoadBGT && load->giveBCValType() == ForceLoadBVT ) {
981 load->computeComponentArrayAt(gVector, tStep, VM_Total);
982 if ( gVector.giveSize() ) {
983 answer.at(1) += coeff * ( b [ 0 ] * gVector.at(1) + c [ 0 ] * gVector.at(2) );
984 answer.at(2) += coeff * ( b [ 1 ] * gVector.at(1) + c [ 1 ] * gVector.at(2) );
985 answer.at(3) += coeff * ( b [ 2 ] * gVector.at(1) + c [ 2 ] * gVector.at(2) );
986 }
987 }
988
989 if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ReinforceBVT ) ) {
990 Reinforcement *rload = dynamic_cast< Reinforcement * >( domain->giveLoad( bodyLoadArray.at(i) ) );
991 double phi = rload->givePorosity();
992 double alpha = rload->giveshapefactor();
993 double kx = rload->givePermeability()->at(1);
994 double ky = rload->givePermeability()->at(2);
995 double tau_0 = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->
996 give( YieldStress, integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
997 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->
998 giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
999
1000 gVector.resize(2);
1001 gVector.at(1) = tau_0 * sqrt(kx * phi) / ( kx * alpha * rho );
1002 gVector.at(2) = tau_0 * sqrt(ky * phi) / ( ky * alpha * rho );
1003
1004 answer.at(1) -= coeff * ( b [ 0 ] * gVector.at(1) + c [ 0 ] * gVector.at(2) );
1005 answer.at(2) -= coeff * ( b [ 1 ] * gVector.at(1) + c [ 1 ] * gVector.at(2) );
1006 answer.at(3) -= coeff * ( b [ 2 ] * gVector.at(1) + c [ 2 ] * gVector.at(2) );
1007 }
1008 }
1009}
1010
1011void
1012TR1_2D_SUPG :: computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
1013{
1014 if ( type != ExternalForcesVector ) {
1015 answer.clear();
1016 return;
1017 }
1018
1019 FloatArray un;
1020 double coeff;
1021
1022 // compute averaged viscosity based on rule of mixture
1023 GaussPoint *gp;
1024 if ( integrationRulesArray [ 0 ]->giveNumberOfIntegrationPoints() ) {
1025 gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
1026 } else {
1027 gp = integrationRulesArray [ 1 ]->getIntegrationPoint(0);
1028 }
1029
1030 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->
1031 giveDensity( gp );
1032
1033 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), un);
1034
1035 double u1sum = un.at(1) + un.at(3) + un.at(5);
1036 double u2sum = un.at(2) + un.at(4) + un.at(6);
1037
1038 answer.resize(9);
1039
1040 if ( load->giveBCValType() == ForceLoadBVT ) {
1041 FloatArray gVector;
1042 load->computeComponentArrayAt(gVector, tStep, VM_Total);
1043
1044 // MB
1045 coeff = rho * area / 3.0;
1046 answer.at(1) = coeff * gVector.at(1) * ( 1.0 + t_supg * ( b [ 0 ] * u1sum + c [ 0 ] * u2sum ) );
1047 answer.at(2) = coeff * gVector.at(2) * ( 1.0 + t_supg * ( b [ 0 ] * u1sum + c [ 0 ] * u2sum ) );
1048 answer.at(4) = coeff * gVector.at(1) * ( 1.0 + t_supg * ( b [ 1 ] * u1sum + c [ 1 ] * u2sum ) );
1049 answer.at(5) = coeff * gVector.at(2) * ( 1.0 + t_supg * ( b [ 1 ] * u1sum + c [ 1 ] * u2sum ) );
1050 answer.at(7) = coeff * gVector.at(1) * ( 1.0 + t_supg * ( b [ 2 ] * u1sum + c [ 2 ] * u2sum ) );
1051 answer.at(8) = coeff * gVector.at(2) * ( 1.0 + t_supg * ( b [ 2 ] * u1sum + c [ 2 ] * u2sum ) );
1052
1053 // MC
1054 coeff = t_pspg * area;
1055 answer.at(3) = coeff * ( b [ 0 ] * gVector.at(1) + c [ 0 ] * gVector.at(2) );
1056 answer.at(6) = coeff * ( b [ 1 ] * gVector.at(1) + c [ 1 ] * gVector.at(2) );
1057 answer.at(9) = coeff * ( b [ 2 ] * gVector.at(1) + c [ 2 ] * gVector.at(2) );
1058
1059 } else if ( load->giveBCValType() == ReinforceBVT ) {
1060 Reinforcement *rload = dynamic_cast< Reinforcement * >( load );
1061 FloatArray t;
1062 double phi = rload->givePorosity();
1063 double alpha = rload->giveshapefactor();
1064 double kx = rload->givePermeability()->at(1);
1065 double ky = rload->givePermeability()->at(2);
1066 double tau_0 = static_cast< FluidCrossSection * >( this->giveCrossSection() )->
1067 giveFluidMaterial()->give( YieldStress, gp );
1068
1069 t.resize(2);
1070 t.at(1) = tau_0 * sqrt(kx * phi) / ( kx * alpha );
1071 t.at(2) = tau_0 * sqrt(ky * phi) / ( ky * alpha );
1072
1073 // MB
1074 answer.at(1) = area * t.at(1) * ( 1.0 + t_supg * ( b [ 0 ] * u1sum + c [ 0 ] * u2sum ) );
1075 answer.at(2) = area * t.at(2) * ( 1.0 + t_supg * ( b [ 0 ] * u1sum + c [ 0 ] * u2sum ) );
1076 answer.at(4) = area * t.at(1) * ( 1.0 + t_supg * ( b [ 1 ] * u1sum + c [ 1 ] * u2sum ) );
1077 answer.at(5) = area * t.at(2) * ( 1.0 + t_supg * ( b [ 1 ] * u1sum + c [ 1 ] * u2sum ) );
1078 answer.at(7) = area * t.at(1) * ( 1.0 + t_supg * ( b [ 2 ] * u1sum + c [ 2 ] * u2sum ) );
1079 answer.at(8) = area * t.at(2) * ( 1.0 + t_supg * ( b [ 2 ] * u1sum + c [ 2 ] * u2sum ) );
1080
1081 // MC
1082 coeff = t_pspg * area / rho;
1083 answer.at(3) = coeff * ( b [ 0 ] * t.at(1) + c [ 0 ] * t.at(2) );
1084 answer.at(6) = coeff * ( b [ 1 ] * t.at(1) + c [ 1 ] * t.at(2) );
1085 answer.at(9) = coeff * ( b [ 2 ] * t.at(1) + c [ 2 ] * t.at(2) );
1086 }
1087}
1088
1089void
1090TR1_2D_SUPG :: updateStabilizationCoeffs(TimeStep *tStep)
1091{
1092 //TR1_2D_SUPG :: updateStabilizationCoeffs (tStep);
1093#if 0
1094 int i, j, k, m, km1, mm1, d1j, d2j, dij, w_dof_addr, u_dof_addr;
1095 double __g_norm, __gamma_norm, __gammav_norm, __beta_norm, __betav_norm, __c_norm, __ctilda_norm, __e_norm, __k_norm, __Re;
1096 double __t_p1, __t_p2, __t_p3, __t_s1, __t_s2, __t_s3;
1097 double nu, rho;
1098 double dudx [ 2 ] [ 2 ], usum [ 2 ];
1099 FloatArray u, un, a;
1100
1101 // compute averaged viscosity based on rule of mixture
1102 GaussPoint *gp;
1103 if ( integrationRulesArray [ 0 ]->giveNumberOfIntegrationPoints() ) {
1104 gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
1105 } else {
1106 gp = integrationRulesArray [ 1 ]->getIntegrationPoint(0);
1107 }
1108
1109 nu = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->
1110 giveCharacteristicValue( MRM_Viscosity, gp, tStep->givePreviousStep() );
1111 rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
1112
1113 //this -> computeVectorOfVelocities(VM_Total,tStep->givePreviousStep(),un) ;
1114 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), u);
1115 this->computeVectorOfVelocities(VM_Acceleration, tStep->givePreviousStep(), a);
1116
1117 un = u;
1118 usum [ 0 ] = un.at(1) + un.at(3) + un.at(5);
1119 usum [ 1 ] = un.at(2) + un.at(4) + un.at(6);
1120
1121 FloatMatrix __tmp;
1122 FloatArray __tmpvec, n;
1123 // assemble g matrix
1124 __tmp.resize(6, 3);
1125 double ar3 = area / 3.0;
1126
1127 __tmp.at(1, 1) = __tmp.at(1, 2) = __tmp.at(1, 3) = b [ 0 ] * ar3;
1128 __tmp.at(3, 1) = __tmp.at(3, 2) = __tmp.at(3, 3) = b [ 1 ] * ar3;
1129 __tmp.at(5, 1) = __tmp.at(5, 2) = __tmp.at(5, 3) = b [ 2 ] * ar3;
1130
1131 __tmp.at(2, 1) = __tmp.at(2, 2) = __tmp.at(2, 3) = c [ 0 ] * ar3;
1132 __tmp.at(4, 1) = __tmp.at(4, 2) = __tmp.at(4, 3) = c [ 1 ] * ar3;
1133 __tmp.at(6, 1) = __tmp.at(6, 2) = __tmp.at(6, 3) = c [ 2 ] * ar3;
1134
1135 __g_norm = __tmp.computeFrobeniusNorm();
1136
1137 // assemble \gamma matrix (advectionTerm of mass conservation eq)
1138 __tmp.resize(3, 6);
1139 dudx [ 0 ] [ 0 ] = b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5);
1140 dudx [ 0 ] [ 1 ] = c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5);
1141 dudx [ 1 ] [ 0 ] = b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6);
1142 dudx [ 1 ] [ 1 ] = c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6);
1143
1144
1145 // dN_epsilon(v)/dv
1146 double coeff = area / 3.;
1147 for ( k = 1; k <= 3; k++ ) { // nodal val of function w
1148 km1 = k - 1;
1149 for ( j = 1; j <= 2; j++ ) { // velocity vector component
1150 for ( m = 1; m <= 3; m++ ) { // nodal components
1151 w_dof_addr = k;
1152 u_dof_addr = ( m - 1 ) * 2 + j;
1153 mm1 = m - 1;
1154 d1j = ( j == 1 );
1155 d2j = ( j == 2 );
1156 __tmp.at(w_dof_addr, u_dof_addr) = coeff * ( 0.0 * d1j * b [ km1 ] * dudx [ 0 ] [ 0 ] + d1j * b [ km1 ] * b [ mm1 ] * usum [ 0 ] +
1157 0.0 * d2j * b [ km1 ] * dudx [ 0 ] [ 1 ] + d1j * b [ km1 ] * c [ mm1 ] * usum [ 1 ] +
1158 0.0 * d1j * c [ km1 ] * dudx [ 1 ] [ 0 ] + d2j * c [ km1 ] * b [ mm1 ] * usum [ 0 ] +
1159 0.0 * d2j * c [ km1 ] * dudx [ 1 ] [ 1 ] + d2j * c [ km1 ] * c [ mm1 ] * usum [ 1 ] );
1160 }
1161 }
1162 }
1163
1164 __gamma_norm = __tmp.computeFrobeniusNorm();
1165 __tmpvec.beProductOf(__tmp, u);
1166 __gammav_norm = __tmpvec.computeNorm();
1167 // compute beta mtrx (acceleration term of mass conservation eq)
1168 __tmp.resize(3, 6);
1169 __tmp.zero();
1170 __tmp.at(1, 1) = __tmp.at(1, 3) = __tmp.at(1, 5) = ar3 * b [ 0 ];
1171 __tmp.at(1, 2) = __tmp.at(1, 4) = __tmp.at(1, 6) = ar3 * c [ 0 ];
1172 __tmp.at(2, 1) = __tmp.at(2, 3) = __tmp.at(2, 5) = ar3 * b [ 1 ];
1173 __tmp.at(2, 2) = __tmp.at(2, 4) = __tmp.at(2, 6) = ar3 * c [ 1 ];
1174 __tmp.at(3, 1) = __tmp.at(3, 3) = __tmp.at(3, 5) = ar3 * b [ 2 ];
1175 __tmp.at(3, 2) = __tmp.at(3, 4) = __tmp.at(3, 6) = ar3 * c [ 2 ];
1176 __beta_norm = __tmp.computeFrobeniusNorm();
1177 __tmpvec.beProductOf(__tmp, a);
1178 __betav_norm = __tmpvec.computeNorm();
1179 // compute c mtrx (advection term of momentum balance)
1180 // standard galerkin term
1181 __tmp.resize(6, 6);
1182 __tmp.zero();
1183
1184 for ( i = 1; i <= 2; i++ ) { // test function index
1185 for ( k = 1; k <= 3; k++ ) { // nodal val of function w
1186 for ( j = 1; j <= 2; j++ ) { // velocity vector component
1187 for ( m = 1; m <= 3; m++ ) { // nodal components
1188 w_dof_addr = ( k - 1 ) * 2 + i;
1189 u_dof_addr = ( m - 1 ) * 2 + j;
1190 d1j = ( j == 1 );
1191 d2j = ( j == 2 );
1192 dij = ( i == j );
1193 coeff = ( m == k ) ? area / 6. : area / 12.;
1194 __tmp.at(w_dof_addr, u_dof_addr) = rho * ( 0.0 * d1j * dudx [ i - 1 ] [ 0 ] * coeff + dij * b [ m - 1 ] * ( area / 12. ) * ( usum [ 0 ] + un.at( ( k - 1 ) * 2 + 1 ) ) +
1195 0.0 * d2j * dudx [ i - 1 ] [ 1 ] * coeff + dij * c [ m - 1 ] * ( area / 12. ) * ( usum [ 1 ] + un.at( ( k - 1 ) * 2 + 2 ) ) );
1196 }
1197 }
1198 }
1199 }
1200
1201 __c_norm = __tmp.computeFrobeniusNorm();
1202 // compute ctilda mtrx (stabilization acceleration term of momentum balance)
1203 __tmp.resize(6, 6);
1204 __tmp.zero();
1205 coeff = rho * area / 12.;
1206
1207 __tmp.at(1, 1) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(1) ) + c [ 0 ] * ( usum [ 1 ] + un.at(2) ) );
1208 __tmp.at(1, 3) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(3) ) + c [ 0 ] * ( usum [ 1 ] + un.at(4) ) );
1209 __tmp.at(1, 5) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(5) ) + c [ 0 ] * ( usum [ 1 ] + un.at(6) ) );
1210
1211 __tmp.at(2, 2) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(1) ) + c [ 0 ] * ( usum [ 1 ] + un.at(2) ) );
1212 __tmp.at(2, 4) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(3) ) + c [ 0 ] * ( usum [ 1 ] + un.at(4) ) );
1213 __tmp.at(2, 6) = coeff * ( b [ 0 ] * ( usum [ 0 ] + un.at(5) ) + c [ 0 ] * ( usum [ 1 ] + un.at(6) ) );
1214
1215 __tmp.at(3, 1) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(1) ) + c [ 1 ] * ( usum [ 1 ] + un.at(2) ) );
1216 __tmp.at(3, 3) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(3) ) + c [ 1 ] * ( usum [ 1 ] + un.at(4) ) );
1217 __tmp.at(3, 5) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(5) ) + c [ 1 ] * ( usum [ 1 ] + un.at(6) ) );
1218
1219 __tmp.at(4, 2) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(1) ) + c [ 1 ] * ( usum [ 1 ] + un.at(2) ) );
1220 __tmp.at(4, 4) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(3) ) + c [ 1 ] * ( usum [ 1 ] + un.at(4) ) );
1221 __tmp.at(4, 6) = coeff * ( b [ 1 ] * ( usum [ 0 ] + un.at(5) ) + c [ 1 ] * ( usum [ 1 ] + un.at(6) ) );
1222
1223 __tmp.at(5, 1) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(1) ) + c [ 2 ] * ( usum [ 1 ] + un.at(2) ) );
1224 __tmp.at(5, 3) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(3) ) + c [ 2 ] * ( usum [ 1 ] + un.at(4) ) );
1225 __tmp.at(5, 5) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(5) ) + c [ 2 ] * ( usum [ 1 ] + un.at(6) ) );
1226
1227 __tmp.at(6, 2) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(1) ) + c [ 2 ] * ( usum [ 1 ] + un.at(2) ) );
1228 __tmp.at(6, 4) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(3) ) + c [ 2 ] * ( usum [ 1 ] + un.at(4) ) );
1229 __tmp.at(6, 6) = coeff * ( b [ 2 ] * ( usum [ 0 ] + un.at(5) ) + c [ 2 ] * ( usum [ 1 ] + un.at(6) ) );
1230
1231 __ctilda_norm = __tmp.computeFrobeniusNorm();
1232
1233 // compute e mtrx (advection term of momentum balance)
1234 __tmp.resize(6, 6);
1235 __tmp.zero();
1236 double __n[] = {
1237 b [ 0 ], c [ 0 ], b [ 1 ], c [ 1 ], b [ 2 ], c [ 2 ]
1238 };
1239
1240 for ( i = 1; i <= 6; i++ ) {
1241 for ( j = 1; j <= 6; j++ ) {
1242 __tmp.at(i, j) = coeff * __n [ i - 1 ] * __n [ j - 1 ];
1243 }
1244 }
1245
1246 __e_norm = __tmp.computeFrobeniusNorm();
1247 // compute element level Reynolds number
1248 // compute k~ norm first
1249 __tmp.resize(6, 6);
1250 __tmp.zero();
1251 FloatMatrix _b(3, 6), _d, _db;
1252 double ar12 = area / 12.;
1253 for ( i = 1; i <= 2; i++ ) { // test function index
1254 for ( k = 1; k <= 3; k++ ) { // nodal val of function w
1255 for ( j = 1; j <= 2; j++ ) { // velocity vector component
1256 for ( m = 1; m <= 3; m++ ) { // nodal components
1257 w_dof_addr = ( k - 1 ) * 2 + i;
1258 u_dof_addr = ( m - 1 ) * 2 + j;
1259 d1j = ( j == 1 );
1260 d2j = ( j == 2 );
1261 dij = ( i == j );
1262 __tmp.at(w_dof_addr, u_dof_addr) += rho *
1263 (
1264 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
1265 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
1266 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
1267 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
1268
1269 0.0 * d1j * b [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
1270 dij * b [ k - 1 ] * b [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 0 ] + un.at(1) * un.at(1) + un.at(3) * un.at(3) + un.at(5) * un.at(5) ) +
1271 0.0 * d2j * b [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 0 ] + un.at( ( m - 1 ) * 2 + 1 ) ) +
1272 dij * b [ k - 1 ] * c [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 1 ] + un.at(1) * un.at(2) + un.at(3) * un.at(4) + un.at(5) * un.at(6) ) +
1273
1274 0.0 * d1j * c [ k - 1 ] * dudx [ i - 1 ] [ 0 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
1275 dij * c [ k - 1 ] * b [ m - 1 ] * ar12 * ( usum [ 0 ] * usum [ 1 ] + un.at(1) * un.at(2) + un.at(3) * un.at(4) + un.at(5) * un.at(6) ) +
1276 0.0 * d2j * c [ k - 1 ] * dudx [ i - 1 ] [ 1 ] * ar12 * ( usum [ 1 ] + un.at( ( m - 1 ) * 2 + 2 ) ) +
1277 dij * c [ k - 1 ] * c [ m - 1 ] * ar12 * ( usum [ 1 ] * usum [ 1 ] + un.at(2) * un.at(2) + un.at(4) * un.at(4) + un.at(6) * un.at(6) )
1278 );
1279 }
1280 }
1281 }
1282 }
1283
1284 __k_norm = __tmp.computeFrobeniusNorm();
1285 double u_1, u_2, vnorm = 0.;
1286 int im1;
1287 for ( i = 1; i <= 3; i++ ) {
1288 im1 = i - 1;
1289 u_1 = u.at( ( im1 ) * 2 + 1 );
1290 u_2 = u.at( ( im1 ) * 2 + 2 );
1291 vnorm = max( vnorm, sqrt(u_1 * u_1 + u_2 * u_2) );
1292 }
1293
1294 if ( vnorm == 0.0 ) {
1295 //t_sugn1 = inf;
1296 double t_sugn2 = tStep->giveTimeIncrement() / 2.0;
1297 //t_sugn3 = inf;
1298 this->t_supg = 1. / sqrt( 1. / ( t_sugn2 * t_sugn2 ) );
1299 this->t_pspg = this->t_supg;
1300 this->t_lsic = 0.0;
1301 //printf ("%e %e\n",this->t_supg,this->t_pspg);
1302 } else {
1303 __Re = vnorm * vnorm * __c_norm / __k_norm / nu;
1304
1305 __t_p1 = __g_norm / __gamma_norm;
1306 __t_p2 = tStep->giveTimeIncrement() * __g_norm / 2.0 / __beta_norm;
1307 __t_p3 = __t_p1 * __Re;
1308 this->t_pspg = 1. / sqrt( 1. / ( __t_p1 * __t_p1 ) + 1. / ( __t_p2 * __t_p2 ) + 1. / ( __t_p3 * __t_p3 ) );
1309
1310
1311 __t_s1 = __c_norm / __k_norm;
1312 __t_s2 = tStep->giveTimeIncrement() * __c_norm / __ctilda_norm / 2.0;
1313 __t_s3 = __t_s1 * __Re;
1314 this->t_supg = 1. / sqrt( 1. / ( __t_s1 * __t_s1 ) + 1. / ( __t_s2 * __t_s2 ) + 1. / ( __t_s3 * __t_s3 ) );
1315
1316 //printf ("%e %e\n",this->t_supg,this->t_pspg);
1317 this->t_lsic = __c_norm / __e_norm;
1318 this->t_lsic = 0.0;
1319 }
1320
1321#else
1322 /* UGN-Based Stabilization */
1323 double h_ugn, sum = 0.0, vnorm, t_sugn1, t_sugn2, t_sugn3, u_1, u_2, z, Re_ugn;
1324 double dscale, uscale, lscale, tscale, dt;
1325 //bool zeroFlag = false;
1326 int im1;
1327 FloatArray u;
1328
1329 uscale = domain->giveEngngModel()->giveVariableScale(VST_Velocity);
1330 lscale = domain->giveEngngModel()->giveVariableScale(VST_Length);
1331 tscale = domain->giveEngngModel()->giveVariableScale(VST_Time);
1332 dscale = domain->giveEngngModel()->giveVariableScale(VST_Density);
1333
1334 this->computeVectorOfVelocities(VM_Total, tStep->givePreviousStep(), u);
1335 u.times(uscale);
1336 double nu;
1337
1338 // compute averaged viscosity based on rule of mixture
1339 GaussPoint *gp;
1340 if ( integrationRulesArray [ 0 ]->giveNumberOfIntegrationPoints() ) {
1341 gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
1342 } else {
1343 gp = integrationRulesArray [ 1 ]->getIntegrationPoint(0);
1344 }
1345
1346 nu = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->giveEffectiveViscosity( gp, tStep->givePreviousStep() );
1347 nu *= domain->giveEngngModel()->giveVariableScale(VST_Viscosity);
1348
1349 dt = tStep->giveTimeIncrement() * tscale;
1350
1351 for ( int i = 1; i <= 3; i++ ) {
1352 im1 = i - 1;
1353 sum += fabs(u.at( ( im1 ) * 2 + 1 ) * b [ im1 ] / lscale + u.at(im1 * 2 + 2) * c [ im1 ] / lscale);
1354 }
1355
1356 /*
1357 * u_1=(u.at(1)+u.at(3)+u.at(5))/3.0;
1358 * u_2=(u.at(2)+u.at(4)+u.at(6))/3.0;
1359 * vnorm=sqrt(u_1*u_1+u_2*u_2);
1360 */
1361 vnorm = 0.;
1362 for ( int i = 1; i <= 3; i++ ) {
1363 im1 = i - 1;
1364 u_1 = u.at( ( im1 ) * 2 + 1 );
1365 u_2 = u.at( ( im1 ) * 2 + 2 );
1366 vnorm = max( vnorm, sqrt(u_1 * u_1 + u_2 * u_2) );
1367 }
1368
1369 if ( ( vnorm == 0.0 ) || ( sum < vnorm * 1e-10 ) ) {
1370 //t_sugn1 = inf;
1371 t_sugn2 = dt / 2.0;
1372 //t_sugn3 = inf;
1373 this->t_supg = 1. / sqrt( 1. / ( t_sugn2 * t_sugn2 ) );
1374 this->t_pspg = this->t_supg;
1375 this->t_lsic = 0.0;
1376 } else {
1377 h_ugn = 2.0 * vnorm / sum;
1378 t_sugn1 = 1. / sum;
1379 t_sugn2 = dt / 2.0;
1380 t_sugn3 = h_ugn * h_ugn / 4.0 / nu;
1381
1382 this->t_supg = 1. / sqrt( 1. / ( t_sugn1 * t_sugn1 ) + 1. / ( t_sugn2 * t_sugn2 ) + 1. / ( t_sugn3 * t_sugn3 ) );
1383 this->t_pspg = this->t_supg;
1384
1385 Re_ugn = vnorm * h_ugn / ( 2. * nu );
1386 z = ( Re_ugn <= 3. ) ? Re_ugn / 3. : 1.0;
1387 this->t_lsic = h_ugn * vnorm * z / 2.0;
1388 }
1389
1390 // if (this->number == 1) {
1391 // printf ("t_supg %e t_pspg %e t_lsic %e\n", t_supg, t_pspg, t_lsic);
1392 // }
1393
1394
1395 this->t_supg *= uscale / lscale;
1396 this->t_pspg *= 1. / ( lscale * dscale );
1397 this->t_lsic *= ( dscale * uscale ) / ( lscale * lscale );
1398
1399 this->t_lsic = 0.0;
1400
1401#endif
1402
1403 //this->t_lsic=0.0;
1404 //this->t_pspg=0.0;
1405}
1406
1407/*
1408 * void
1409 * TR1_2D_SUPG :: updateStabilizationCoeffs (TimeStep* tStep)
1410 * {
1411 * double h_ugn, sum=0.0, snorm, vnorm, t_sugn1, t_sugn2, t_sugn3, u_1, u_2, z, Re_ugn;
1412 * bool zeroFlag = false;
1413 * int i,im1;
1414 * FloatArray u;
1415 * this -> computeVectorOfVelocities(VM_Total,tStep, u) ;
1416 * double nu = this->giveMaterial()->giveCharacteristicValue(MRM_Viscosity, integrationRulesArray[0]->getIntegrationPoint(0), tStep);
1417 *
1418 * // UGN-Based Stabilization
1419 *
1420 * for (i=1;i<=3;i++) {
1421 * im1=i-1;
1422 * snorm = sqrt(u.at(im1*2+1)*u.at(im1*2+1) + u.at(im1*2+2)*u.at(im1*2+2));
1423 * if (snorm == 0.0) zeroFlag = true;
1424 * sum+= fabs(u.at((im1)*2+1)*b[im1] + u.at(im1*2+2)*c[im1])/snorm;
1425 * }
1426 * u_1=(u.at(1)+u.at(3)+u.at(5))/3.0;
1427 * u_2=(u.at(2)+u.at(4)+u.at(6))/3.0;
1428 * vnorm=sqrt(u_1*u_1+u_2*u_2);
1429 *
1430 * if (zeroFlag) {
1431 *
1432 * //t_sugn1 = inf;
1433 * t_sugn2 = tStep->giveTimeIncrement()/2.0;
1434 * //t_sugn3 = inf;
1435 * this->t_supg=1./sqrt(1./(t_sugn2*t_sugn2));
1436 * this->t_pspg=this->t_supg;
1437 * this->t_lsic=0.0;
1438 *
1439 * } else {
1440 * h_ugn = 2.0*vnorm/sum;
1441 * //t_sugn1 = h_ugn/(2.*vnorm);
1442 * t_sugn1 = 1./sum;
1443 * t_sugn2 = tStep->giveTimeIncrement()/2.0;
1444 * t_sugn3 = h_ugn*h_ugn/4.0/nu;
1445 *
1446 * this->t_supg=1./sqrt(1./(t_sugn1*t_sugn1) + 1./(t_sugn2*t_sugn2) + 1./(t_sugn3*t_sugn3));
1447 * this->t_pspg=this->t_supg;
1448 *
1449 * Re_ugn = vnorm*h_ugn/(2.*nu);
1450 * z = (Re_ugn <= 3.)?Re_ugn/3. : 1.0 ;
1451 * this->t_lsic=h_ugn*vnorm*z/2.0;
1452 * this->t_lsic = 0.0;
1453 * }
1454 * //printf ("Element %d ( t_supg %e, t_pspg %e, t_lsic %e)\n", this->giveNumber(), t_supg, t_pspg, t_lsic);
1455 * }
1456 *
1457 */
1458
1459double
1460TR1_2D_SUPG :: computeCriticalTimeStep(TimeStep *tStep)
1461{
1462 return 1.e6;
1463}
1464
1465
1466Interface *
1467TR1_2D_SUPG :: giveInterface(InterfaceType interface)
1468{
1469 if ( interface == ZZNodalRecoveryModelInterfaceType ) {
1470 return static_cast< ZZNodalRecoveryModelInterface * >(this);
1471 } else if ( interface == NodalAveragingRecoveryModelInterfaceType ) {
1472 return static_cast< NodalAveragingRecoveryModelInterface * >(this);
1473 } else if ( interface == SPRNodalRecoveryModelInterfaceType ) {
1474 return static_cast< SPRNodalRecoveryModelInterface * >(this);
1475 } else if ( interface == SpatialLocalizerInterfaceType ) {
1476 return static_cast< SpatialLocalizerInterface * >(this);
1477 } else if ( interface == EIPrimaryFieldInterfaceType ) {
1478 return static_cast< EIPrimaryFieldInterface * >(this);
1479 } else if ( interface == LEPlicElementInterfaceType ) {
1480 return static_cast< LEPlicElementInterface * >(this);
1481 } else if ( interface == LevelSetPCSElementInterfaceType ) {
1482 return static_cast< LevelSetPCSElementInterface * >(this);
1483 }
1484
1485 return NULL;
1486}
1487
1488
1489void
1490TR1_2D_SUPG :: computeDeviatoricStrain(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
1491{
1492 /* one should call material driver instead */
1493 FloatArray u(6);
1494 answer.resize(3);
1495
1496
1497 this->computeVectorOfVelocities(VM_Total, tStep, u);
1498
1499 answer.at(1) = ( b [ 0 ] * u.at(1) + b [ 1 ] * u.at(3) + b [ 2 ] * u.at(5) );
1500 answer.at(2) = ( c [ 0 ] * u.at(2) + c [ 1 ] * u.at(4) + c [ 2 ] * u.at(6) );
1501 answer.at(3) = ( b [ 0 ] * u.at(2) + b [ 1 ] * u.at(4) + b [ 2 ] * u.at(6) + c [ 0 ] * u.at(1) + c [ 1 ] * u.at(3) + c [ 2 ] * u.at(5) );
1502}
1503
1504void
1505TR1_2D_SUPG :: initGeometry()
1506{
1507 Node *node1, *node2, *node3;
1508 double x1, x2, x3, y1, y2, y3;
1509
1510 node1 = giveNode(1);
1511 node2 = giveNode(2);
1512 node3 = giveNode(3);
1513
1514 // init geometry data
1515 x1 = node1->giveCoordinate(1);
1516 x2 = node2->giveCoordinate(1);
1517 x3 = node3->giveCoordinate(1);
1518
1519 y1 = node1->giveCoordinate(2);
1520 y2 = node2->giveCoordinate(2);
1521 y3 = node3->giveCoordinate(2);
1522
1523 this->area = 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 );
1524
1525 if ( area < 0.0 ) {
1526 OOFEM_ERROR("Area is negative, check element numbering orientation");
1527 }
1528
1529 b [ 0 ] = ( y2 - y3 ) / ( 2. * area );
1530 c [ 0 ] = ( x3 - x2 ) / ( 2. * area );
1531 b [ 1 ] = ( y3 - y1 ) / ( 2. * area );
1532 c [ 1 ] = ( x1 - x3 ) / ( 2. * area );
1533 b [ 2 ] = ( y1 - y2 ) / ( 2. * area );
1534 c [ 2 ] = ( x2 - x1 ) / ( 2. * area );
1535}
1536
1537
1538int
1539TR1_2D_SUPG :: checkConsistency()
1540{
1541 return SUPGElement :: checkConsistency();
1542}
1543
1544void
1545TR1_2D_SUPG :: computeNMtrx(FloatArray &answer, GaussPoint *gp)
1546{
1547 double l1, l2;
1548 answer.resize(3);
1549
1550 answer.at(1) = l1 = gp->giveNaturalCoordinate(1);
1551 answer.at(2) = l2 = gp->giveNaturalCoordinate(2);
1552 answer.at(3) = 1.0 - l1 - l2;
1553}
1554
1555/*
1556 * double
1557 * TR1_2D_SUPG :: computeCriticalTimeStep (TimeStep* tStep)
1558 * {
1559 * FloatArray u;
1560 * double dt1, dt2, dt;
1561 * double Re = static_cast<FluidModel*>(domain->giveEngngModel())->giveReynoldsNumber();
1562 *
1563 * this -> computeVectorOfVelocities(VM_Total,tStep, u) ;
1564 *
1565 * double vn1 = sqrt(u.at(1)*u.at(1)+u.at(2)*u.at(2));
1566 * double vn2 = sqrt(u.at(3)*u.at(3)+u.at(4)*u.at(4));
1567 * double vn3 = sqrt(u.at(5)*u.at(5)+u.at(6)*u.at(6));
1568 * double veln = max (vn1, max(vn2,vn3));
1569 *
1570 * double l1 = 1.0/(sqrt(b[0]*b[0]+c[0]*c[0]));
1571 * double l2 = 1.0/(sqrt(b[1]*b[1]+c[1]*c[1]));
1572 * double l3 = 1.0/(sqrt(b[2]*b[2]+c[2]*c[2]));
1573 *
1574 * double ln = min (l1, min (l2,l3));
1575 *
1576 * // viscous limit
1577 * dt2 = 0.5*ln*ln*Re;
1578 * if (veln != 0.0) {
1579 * dt1 = ln/veln;
1580 * dt = dt1*dt2/(dt1+dt2);
1581 * } else {
1582 * dt = dt2;
1583 * }
1584 * return dt;
1585 * }
1586 */
1587
1588
1589/* Note this LEPLIC implementation will not work for axi elements as the true volume
1590 * is not taken into account (the effect of radius) !
1591 */
1592
1593
1594double
1595TR1_2D_SUPG :: computeLEPLICVolumeFraction(const FloatArray &n, const double p, LEPlic *matInterface, bool updFlag)
1596{
1597 Polygon pg;
1598 double answer, volume = computeMyVolume(matInterface, updFlag);
1599 this->formVolumeInterfacePoly(pg, matInterface, n, p, updFlag);
1600 answer = fabs(pg.computeVolume() / volume);
1601 if ( answer > 1.0 ) {
1602 return 1.0;
1603 } else {
1604 return answer;
1605 }
1606}
1607
1608void
1609TR1_2D_SUPG :: formMaterialVolumePoly(Polygon &matvolpoly, LEPlic *matInterface,
1610 const FloatArray &normal, const double p, bool updFlag)
1611{
1612 double x, y;
1613 Vertex v;
1614
1615 matvolpoly.clear();
1616
1617 if ( this->vof <= TRSUPG_ZERO_VOF ) {
1618 return;
1619 } else if ( this->vof >= ( 1 - TRSUPG_ZERO_VOF ) ) {
1620 for ( int i = 1; i <= 3; i++ ) {
1621 if ( updFlag ) {
1622 x = matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() );
1623 y = matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() );
1624 } else {
1625 x = this->giveNode(i)->giveCoordinate(1);
1626 y = this->giveNode(i)->giveCoordinate(2);
1627 }
1628
1629 v.setCoords(x, y);
1630 matvolpoly.addVertex(v);
1631 }
1632
1633 return;
1634 }
1635
1636 this->formVolumeInterfacePoly(matvolpoly, matInterface, normal, p, updFlag);
1637}
1638
1639
1640void
1641TR1_2D_SUPG :: formVolumeInterfacePoly(Polygon &matvolpoly, LEPlic *matInterface,
1642 const FloatArray &normal, const double p, bool updFlag)
1643{
1644 int next;
1645 bool nodeIn [ 3 ];
1646 double nx = normal.at(1), ny = normal.at(2), x, y;
1647 double tx, ty;
1648 Vertex v;
1649
1650 matvolpoly.clear();
1651
1652 for ( int i = 1; i <= 3; i++ ) {
1653 if ( updFlag ) {
1654 x = matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() );
1655 y = matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() );
1656 } else {
1657 x = this->giveNode(i)->giveCoordinate(1);
1658 y = this->giveNode(i)->giveCoordinate(2);
1659 }
1660
1661 if ( ( nx * x + ny * y + p ) >= 0. ) {
1662 nodeIn [ i - 1 ] = true;
1663 } else {
1664 nodeIn [ i - 1 ] = false;
1665 }
1666 }
1667
1668 if ( nodeIn [ 0 ] && nodeIn [ 1 ] && nodeIn [ 2 ] ) { // all nodes inside
1669 for ( int i = 1; i <= 3; i++ ) {
1670 if ( updFlag ) {
1671 x = matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() );
1672 y = matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() );
1673 } else {
1674 x = this->giveNode(i)->giveCoordinate(1);
1675 y = this->giveNode(i)->giveCoordinate(2);
1676 }
1677
1678 v.setCoords(x, y);
1679 matvolpoly.addVertex(v);
1680 }
1681
1682 return;
1683 } else if ( !( nodeIn [ 0 ] || nodeIn [ 1 ] || nodeIn [ 2 ] ) ) { // all nodes outside
1684 return;
1685 } else {
1686 for ( int i = 1; i <= 3; i++ ) {
1687 next = i < 3 ? i + 1 : 1;
1688 if ( nodeIn [ i - 1 ] ) {
1689 if ( updFlag ) {
1690 v.setCoords( matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() ),
1691 matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() ) );
1692 } else {
1693 v.setCoords( this->giveNode(i)->giveCoordinate(1),
1694 this->giveNode(i)->giveCoordinate(2) );
1695 }
1696
1697 matvolpoly.addVertex(v);
1698 }
1699
1700 if ( nodeIn [ next - 1 ] ^ nodeIn [ i - 1 ] ) {
1701 // compute intersection with (i,next) edge
1702 if ( updFlag ) {
1703 x = matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() );
1704 y = matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() );
1705 tx = matInterface->giveUpdatedXCoordinate( this->giveNode(next)->giveNumber() ) - x;
1706 ty = matInterface->giveUpdatedYCoordinate( this->giveNode(next)->giveNumber() ) - y;
1707 } else {
1708 x = this->giveNode(i)->giveCoordinate(1);
1709 y = this->giveNode(i)->giveCoordinate(2);
1710 tx = this->giveNode(next)->giveCoordinate(1) - x;
1711 ty = this->giveNode(next)->giveCoordinate(2) - y;
1712 }
1713
1714 double s, sd = nx * tx + ny * ty;
1715 if ( fabs(sd) > 1.e-6 ) {
1716 s = ( -p - ( nx * x + ny * y ) ) / sd;
1717 v.setCoords(x + tx * s, y + ty * s);
1718 matvolpoly.addVertex(v);
1719 } else {
1720 // pathological case - lines are parallel
1721 if ( nodeIn [ i - 1 ] ) {
1722 if ( updFlag ) {
1723 v.setCoords( matInterface->giveUpdatedXCoordinate( this->giveNode(next)->giveNumber() ),
1724 matInterface->giveUpdatedYCoordinate( this->giveNode(next)->giveNumber() ) );
1725 } else {
1726 v.setCoords( this->giveNode(next)->giveCoordinate(1), this->giveNode(next)->giveCoordinate(2) );
1727 }
1728
1729 matvolpoly.addVertex(v);
1730 } else {
1731 v.setCoords(x, y);
1732 matvolpoly.addVertex(v);
1733 if ( updFlag ) {
1734 v.setCoords( matInterface->giveUpdatedXCoordinate( this->giveNode(next)->giveNumber() ),
1735 matInterface->giveUpdatedYCoordinate( this->giveNode(next)->giveNumber() ) );
1736 } else {
1737 v.setCoords( this->giveNode(next)->giveCoordinate(1), this->giveNode(next)->giveCoordinate(2) );
1738 }
1739
1740 matvolpoly.addVertex(v);
1741 }
1742 }
1743 }
1744 } // end loop over elem nodes
1745 }
1746}
1747
1748
1749double
1750TR1_2D_SUPG :: truncateMatVolume(const Polygon &matvolpoly, double &volume)
1751{
1752 Polygon me, clip;
1753 Graph g;
1754
1755 this->formMyVolumePoly(me, NULL, false);
1756 g.clip(clip, me, matvolpoly);
1757#ifdef __OOFEG
1758 EASValsSetColor( gc [ 0 ].getActiveCrackColor() );
1759 //GraphicObj *go = clip.draw(::gc[OOFEG_DEBUG_LAYER],true);
1760 clip.draw(gc [ OOFEG_DEBUG_LAYER ], true);
1761 //EVFastRedraw(myview);
1762#endif
1763 volume = clip.computeVolume();
1764 return volume / area;
1765}
1766
1767void
1768TR1_2D_SUPG :: formMyVolumePoly(Polygon &me, LEPlic *matInterface, bool updFlag)
1769{
1770 double x, y;
1771 Vertex v;
1772
1773 me.clear();
1774
1775 for ( int i = 1; i <= 3; i++ ) {
1776 if ( updFlag ) {
1777 x = matInterface->giveUpdatedXCoordinate( this->giveNode(i)->giveNumber() );
1778 y = matInterface->giveUpdatedYCoordinate( this->giveNode(i)->giveNumber() );
1779 } else {
1780 x = this->giveNode(i)->giveCoordinate(1);
1781 y = this->giveNode(i)->giveCoordinate(2);
1782 }
1783
1784 v.setCoords(x, y);
1785 me.addVertex(v);
1786 }
1787}
1788
1789
1790double
1791TR1_2D_SUPG :: computeMyVolume(LEPlic *matInterface, bool updFlag)
1792{
1793 double x1, x2, x3, y1, y2, y3;
1794 if ( updFlag ) {
1795 x1 = matInterface->giveUpdatedXCoordinate( this->giveNode(1)->giveNumber() );
1796 x2 = matInterface->giveUpdatedXCoordinate( this->giveNode(2)->giveNumber() );
1797 x3 = matInterface->giveUpdatedXCoordinate( this->giveNode(3)->giveNumber() );
1798
1799 y1 = matInterface->giveUpdatedYCoordinate( this->giveNode(1)->giveNumber() );
1800 y2 = matInterface->giveUpdatedYCoordinate( this->giveNode(2)->giveNumber() );
1801 y3 = matInterface->giveUpdatedYCoordinate( this->giveNode(3)->giveNumber() );
1802 return 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 );
1803 } else {
1804 return area;
1805 }
1806}
1807
1808double
1809TR1_2D_SUPG :: computeVolumeAround(GaussPoint *gp)
1810// Returns the portion of the receiver which is attached to gp.
1811{
1812 double determinant = fabs( 4 * area * area * ( c [ 1 ] * b [ 0 ] - c [ 0 ] * b [ 1 ] ) );
1813
1814 return gp->giveWeight() * determinant;
1815}
1816
1817
1818double
1819TR1_2D_SUPG :: computeCriticalLEPlicTimeStep(TimeStep *tStep)
1820{
1821 FloatArray u;
1822 double Re = static_cast< FluidModel * >( domain->giveEngngModel() )->giveReynoldsNumber();
1823
1824 this->computeVectorOfVelocities(VM_Total, tStep, u);
1825
1826 double vn1 = sqrt( u.at(1) * u.at(1) + u.at(2) * u.at(2) );
1827 double vn2 = sqrt( u.at(3) * u.at(3) + u.at(4) * u.at(4) );
1828 double vn3 = sqrt( u.at(5) * u.at(5) + u.at(6) * u.at(6) );
1829 double veln = max( vn1, max(vn2, vn3) );
1830
1831 double l1 = 1.0 / ( sqrt(b [ 0 ] * b [ 0 ] + c [ 0 ] * c [ 0 ]) );
1832 double l2 = 1.0 / ( sqrt(b [ 1 ] * b [ 1 ] + c [ 1 ] * c [ 1 ]) );
1833 double l3 = 1.0 / ( sqrt(b [ 2 ] * b [ 2 ] + c [ 2 ] * c [ 2 ]) );
1834
1835 double ln = min( l1, min(l2, l3) );
1836
1837 if ( veln != 0.0 ) {
1838 return ln / veln;
1839 } else {
1840 return 0.5 * ln * ln * Re;
1841 }
1842}
1843
1844
1845void
1846TR1_2D_SUPG :: giveElementCenter(LEPlic *mat_interface, FloatArray &center, bool upd)
1847{
1848 FloatArray coords;
1849
1850 center.resize(2);
1851 center.zero();
1852 if ( upd ) {
1853 for ( int i = 1; i <= 3; i++ ) {
1854 mat_interface->giveUpdatedCoordinate( coords, giveNode(i)->giveNumber() );
1855 center.add(coords);
1856 }
1857 } else {
1858 for ( int i = 1; i <= 3; i++ ) {
1859 center.at(1) += this->giveNode(i)->giveCoordinate(1);
1860 center.at(2) += this->giveNode(i)->giveCoordinate(2);
1861 }
1862 }
1863
1864 center.times(1. / 3.);
1865}
1866
1867int
1868TR1_2D_SUPG :: EIPrimaryFieldI_evaluateFieldVectorAt(FloatArray &answer, PrimaryField &pf,
1869 const FloatArray &coords, IntArray &dofId, ValueModeType mode,
1870 TimeStep *tStep)
1871{
1872 int indx, es;
1873 double sum;
1874 FloatArray elemvector, f, lc;
1875 //FloatMatrix n;
1876 IntArray elemdofs;
1877 // determine element dof ids
1878 this->giveElementDofIDMask(elemdofs);
1879 es = elemdofs.giveSize();
1880 // first evaluate element unknown vector
1881 this->computeVectorOf(pf, elemdofs, mode, tStep, elemvector);
1882 // determine corresponding local coordinates
1883 if ( this->computeLocalCoordinates(lc, coords) ) {
1884 // compute interpolation matrix
1885 // this->computeNmatrixAt(n, &lc);
1886 // compute answer
1887 answer.resize( dofId.giveSize() );
1888 for ( int i = 1; i <= dofId.giveSize(); i++ ) {
1889 if ( ( indx = elemdofs.findFirstIndexOf( dofId.at(i) ) ) ) {
1890 sum = 0.0;
1891 for ( int j = 1; j <= 3; j++ ) {
1892 sum += lc.at(j) * elemvector.at(es * ( j - 1 ) + indx);
1893 }
1894
1895 answer.at(i) = sum;
1896 } else {
1897 // OOFEM_ERROR("unknown dof id encountered");
1898 answer.at(i) = 0.0;
1899 }
1900 }
1901
1902 return 0; // ok
1903 } else {
1904 OOFEM_ERROR("target point not in receiver volume");
1905 return 1; // fail
1906 }
1907}
1908
1909
1910void
1911TR1_2D_SUPG :: updateYourself(TimeStep *tStep)
1912{
1913 SUPGElement :: updateYourself(tStep);
1914 LEPlicElementInterface :: updateYourself(tStep);
1915}
1916
1917int
1918TR1_2D_SUPG :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
1919{
1920 if ( type == IST_VOFFraction ) {
1921 MaterialInterface *mi = domain->giveEngngModel()->giveMaterialInterface( domain->giveNumber() );
1922 if ( mi ) {
1923 FloatArray val;
1924 mi->giveElementMaterialMixture( val, gp->giveElement()->giveNumber() );
1925 answer.resize(1);
1926 answer.at(1) = val.at(1);
1927 return 1;
1928 } else {
1929 answer.resize(1);
1930 answer.at(1) = 1.0;
1931 return 1;
1932 }
1933 } else {
1934 return SUPGElement :: giveIPValue(answer, gp, type, tStep);
1935 }
1936}
1937
1938void
1939TR1_2D_SUPG :: NodalAveragingRecoveryMI_computeNodalValue(FloatArray &answer, int node,
1940 InternalStateType type, TimeStep *tStep)
1941{
1942 GaussPoint *gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
1943 this->giveIPValue(answer, gp, type, tStep);
1944}
1945
1946void
1947TR1_2D_SUPG :: SPRNodalRecoveryMI_giveSPRAssemblyPoints(IntArray &pap)
1948{
1949 pap.resize(3);
1950 pap.at(1) = this->giveNode(1)->giveNumber();
1951 pap.at(2) = this->giveNode(2)->giveNumber();
1952 pap.at(3) = this->giveNode(3)->giveNumber();
1953}
1954
1955void
1956TR1_2D_SUPG :: SPRNodalRecoveryMI_giveDofMansDeterminedByPatch(IntArray &answer, int pap)
1957{
1958 answer.resize(1);
1959 if ( ( pap == this->giveNode(1)->giveNumber() ) ||
1960 ( pap == this->giveNode(2)->giveNumber() ) ||
1961 ( pap == this->giveNode(3)->giveNumber() ) ) {
1962 answer.at(1) = pap;
1963 } else {
1964 OOFEM_ERROR("node unknown");
1965 }
1966}
1967
1968int
1969TR1_2D_SUPG :: SPRNodalRecoveryMI_giveNumberOfIP() { return 1; }
1970
1971
1972SPRPatchType
1973TR1_2D_SUPG :: SPRNodalRecoveryMI_givePatchType()
1974{
1975 return SPRPatchType_2dxy;
1976}
1977
1978
1979
1980void
1981TR1_2D_SUPG :: printOutputAt(FILE *file, TimeStep *tStep)
1982// Performs end-of-step operations.
1983{
1984 SUPGElement :: printOutputAt(file, tStep);
1985 double rho = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveDensity( integrationRulesArray [ 0 ]->getIntegrationPoint(0) );
1986 fprintf(file, "\telement_status { VOF %e, density %e }\n\n", this->giveVolumeFraction(), rho);
1987}
1988
1989
1990void TR1_2D_SUPG :: saveContext(DataStream &stream, ContextMode mode)
1991{
1992 SUPGElement :: saveContext(stream, mode);
1993 LEPlicElementInterface :: saveContext(stream, mode);
1994}
1995
1996
1997void TR1_2D_SUPG :: restoreContext(DataStream &stream, ContextMode mode)
1998{
1999 SUPGElement :: restoreContext(stream, mode);
2000 LEPlicElementInterface :: restoreContext(stream, mode);
2001}
2002
2003
2004double
2005TR1_2D_SUPG :: LS_PCS_computeF(LevelSetPCS *ls, TimeStep *tStep)
2006{
2007 double answer;
2008 FloatArray fi(3), un;
2009
2010 this->computeVectorOfVelocities(VM_Total, tStep, un);
2011
2012 for ( int i = 1; i <= 3; i++ ) {
2013 fi.at(i) = ls->giveLevelSetDofManValue( dofManArray.at(i) );
2014 }
2015
2016 double fix = b [ 0 ] * fi.at(1) + b [ 1 ] * fi.at(2) + b [ 2 ] * fi.at(3);
2017 double fiy = c [ 0 ] * fi.at(1) + c [ 1 ] * fi.at(2) + c [ 2 ] * fi.at(3);
2018 double norm = sqrt(fix * fix + fiy * fiy);
2019
2020 answer = ( 1. / 3. ) * ( fix * ( un.at(1) + un.at(3) + un.at(5) ) + fiy * ( un.at(2) + un.at(4) + un.at(6) ) ) / norm;
2021 return answer;
2022}
2023
2024
2025double
2026TR1_2D_SUPG :: LS_PCS_computeS(LevelSetPCS *ls, TimeStep *tStep)
2027{
2028#if 0
2029 double fi, answer, eps = 0.0;
2030 FloatArray s(3);
2031
2032 for ( int i = 1; i <= 3; i++ ) {
2033 fi = ls->giveLevelSetDofManValue( dofManArray.at(i) );
2034 s.at(i) = fi / sqrt(fi * fi + eps * eps);
2035 }
2036
2037 answer = ( 1. / 3. ) * ( s.at(1) + s.at(2) + s.at(3) );
2038 return answer;
2039
2040#else
2041
2042 int neg = 0, pos = 0, zero = 0, si = 0;
2043 double x1, x2, x3, y1, y2, y3;
2044 FloatArray fi(3);
2045
2046 for ( int i = 1; i <= 3; i++ ) {
2047 fi.at(i) = ls->giveLevelSetDofManValue( dofManArray.at(i) );
2048 }
2049
2050
2051 for ( int i = 1; i <= 3; i++ ) {
2052 if ( fi.at(i) >= 0. ) {
2053 pos++;
2054 }
2055
2056 if ( fi.at(i) < 0.0 ) {
2057 neg++;
2058 }
2059
2060 if ( fi.at(i) == 0. ) {
2061 zero++;
2062 }
2063 }
2064
2065 if ( zero == 3 ) {
2066 return 0.0;
2067 } else if ( neg == 0 ) { // all level set values positive
2068 return 1.0; //return area;
2069 } else if ( pos == 0 ) { // all level set values negative
2070 return -1.0; //return -area;
2071 } else {
2072 // zero level set inside
2073 // find the vertex vith level set sign different from other two
2074 for ( int i = 1; i <= 3; i++ ) {
2075 if ( ( pos > neg ) && ( fi.at(i) < 0.0 ) ) {
2076 si = i;
2077 break;
2078 }
2079
2080 if ( ( pos < neg ) && ( fi.at(i) >= 0.0 ) ) {
2081 si = i;
2082 break;
2083 }
2084 }
2085
2086 if ( si ) {
2087 x1 = this->giveNode(si)->giveCoordinate(1);
2088 y1 = this->giveNode(si)->giveCoordinate(2);
2089
2090 // compute intersections
2091 int prev_node = ( si > 1 ) ? si - 1 : 3;
2092 int next_node = ( si < 3 ) ? si + 1 : 1;
2093
2094 //double l = distance( *this->giveNode(si)->giveCoordinates(), *this->giveNode(next_node)->giveCoordinates() );
2095 double t = fi.at(si) / ( fi.at(si) - fi.at(next_node) );
2096 x2 = x1 + t * ( this->giveNode(next_node)->giveCoordinate(1) - x1 );
2097 y2 = y1 + t * ( this->giveNode(next_node)->giveCoordinate(2) - y1 );
2098
2099 //l = distance( this->giveNode(si)->giveCoordinates(), *this->giveNode(prev_node)->giveCoordinates() );
2100 t = fi.at(si) / ( fi.at(si) - fi.at(prev_node) );
2101 x3 = x1 + t * ( this->giveNode(prev_node)->giveCoordinate(1) - x1 );
2102 y3 = y1 + t * ( this->giveNode(prev_node)->giveCoordinate(2) - y1 );
2103
2104 // compute area
2105 double __area = fabs( 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 ) );
2106
2107 if ( pos > neg ) {
2108 // negative area computed
2109 return ( ( area - __area ) - __area ) / area;
2110 } else {
2111 // postive area computed
2112 return ( __area - ( area - __area ) ) / area;
2113 }
2114 } else {
2115 OOFEM_ERROR("internal consistency error");
2116 return 0.0;
2117 }
2118 }
2119
2120#endif
2121}
2122
2123
2124void
2125TR1_2D_SUPG :: LS_PCS_computedN(FloatMatrix &answer)
2126{
2127 answer.resize(3, 2);
2128
2129 for ( int i = 1; i <= 3; i++ ) {
2130 answer.at(i, 1) = b [ i - 1 ];
2131 answer.at(i, 2) = c [ i - 1 ];
2132 }
2133}
2134
2135
2136void
2137TR1_2D_SUPG :: LS_PCS_computeVOFFractions(FloatArray &answer, FloatArray &fi)
2138{
2139 int neg = 0, pos = 0, zero = 0, si = 0;
2140 double x1, x2, x3, y1, y2, y3;
2141
2142 answer.resize(2);
2143 for ( int i = 1; i <= 3; i++ ) {
2144 if ( fi.at(i) >= 0. ) {
2145 pos++;
2146 }
2147
2148 if ( fi.at(i) < 0.0 ) {
2149 neg++;
2150 }
2151
2152 if ( fi.at(i) == 0. ) {
2153 zero++;
2154 }
2155 }
2156
2157 if ( neg == 0 ) { // all level set values positive
2158 answer.at(1) = 1.0;
2159 answer.at(2) = 0.0;
2160 } else if ( pos == 0 ) { // all level set values negative
2161 answer.at(1) = 0.0;
2162 answer.at(2) = 1.0;
2163 } else if ( zero == 3 ) {
2164 // ???????
2165 answer.at(1) = 1.0;
2166 answer.at(2) = 0.0;
2167 } else {
2168 // zero level set inside
2169 // find the vertex with level set sign different from other two
2170 for ( int i = 1; i <= 3; i++ ) {
2171 if ( ( pos > neg ) && ( fi.at(i) < 0.0 ) ) {
2172 si = i;
2173 break;
2174 }
2175
2176 if ( ( pos < neg ) && ( fi.at(i) >= 0.0 ) ) {
2177 si = i;
2178 break;
2179 }
2180 }
2181
2182 if ( si && ( ( pos + neg ) == 3 ) ) {
2183 x1 = this->giveNode(si)->giveCoordinate(1);
2184 y1 = this->giveNode(si)->giveCoordinate(2);
2185
2186 // compute intersections
2187 int prev_node = ( si > 1 ) ? si - 1 : 3;
2188 int next_node = ( si < 3 ) ? si + 1 : 1;
2189
2190 double t = fi.at(si) / ( fi.at(si) - fi.at(next_node) );
2191 x2 = x1 + t * ( this->giveNode(next_node)->giveCoordinate(1) - this->giveNode(si)->giveCoordinate(1) );
2192 y2 = y1 + t * ( this->giveNode(next_node)->giveCoordinate(2) - this->giveNode(si)->giveCoordinate(2) );
2193
2194 t = fi.at(si) / ( fi.at(si) - fi.at(prev_node) );
2195 x3 = x1 + t * ( this->giveNode(prev_node)->giveCoordinate(1) - this->giveNode(si)->giveCoordinate(1) );
2196 y3 = y1 + t * ( this->giveNode(prev_node)->giveCoordinate(2) - this->giveNode(si)->giveCoordinate(2) );
2197
2198 // compute area
2199 double __area = 0.5 * ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 );
2200 if ( fabs(__area) / area > 1.00001 ) {
2201 OOFEM_ERROR("internal consistency error");
2202 }
2203
2204 // prevent some roundoff errors
2205 if ( fabs(__area) > area ) {
2206 __area = sgn(__area) * area;
2207 }
2208
2209 if ( pos > neg ) {
2210 // negative area computed
2211 answer.at(2) = fabs(__area) / area;
2212 answer.at(1) = 1.0 - answer.at(2);
2213 } else {
2214 // postive area computed
2215 answer.at(1) = fabs(__area) / area;
2216 answer.at(2) = 1.0 - answer.at(1);
2217 }
2218 } else {
2219 OOFEM_ERROR("internal consistency error");
2220 }
2221 }
2222}
2223
2224
2225void
2226TR1_2D_SUPG :: giveLocalVelocityDofMap(IntArray &map)
2227{
2228 map = {1, 2, 4, 5, 7, 8};
2229}
2230
2231void
2232TR1_2D_SUPG :: giveLocalPressureDofMap(IntArray &map)
2233{
2234 map = {3, 6, 9};
2235}
2236
2237
2238void
2239TR1_2D_SUPG :: computeDeviatoricStress(FloatArray &answer, const FloatArray &eps, GaussPoint *gp, TimeStep *tStep)
2240{
2241 answer = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeDeviatoricStress2D(eps, gp, tStep);
2242}
2243
2244void
2245TR1_2D_SUPG :: computeTangent(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
2246{
2247 answer = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial()->computeTangent2D(mode, gp, tStep);
2248}
2249
2250
2251#ifdef __OOFEG
2252int
2253TR1_2D_SUPG :: giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode,
2254 int node, TimeStep *tStep)
2255{
2256 /*
2257 * if (type == IST_VOFFraction) {
2258 * answer.resize(1);
2259 * answer.at(1) = this->giveTempVolumeFraction();
2260 * return 1;
2261 * } else if (type == IST_Density) {
2262 * answer.resize(1);
2263 * answer.at(1) = this->giveMaterial()->giveCharacteristicValue(MRM_Density, integrationRulesArray[0]-> getIntegrationPoint(0), tStep);
2264 * return 1;
2265 *
2266 * } else
2267 */
2268 return SUPGElement :: giveInternalStateAtNode(answer, type, mode, node, tStep);
2269}
2270
2271void
2272TR1_2D_SUPG :: drawRawGeometry(oofegGraphicContext &gc, TimeStep *tStep)
2273{
2274 WCRec p [ 3 ];
2275 GraphicObj *go;
2276
2277 if ( !gc.testElementGraphicActivity(this) ) {
2278 return;
2279 }
2280
2281 EASValsSetLineWidth(OOFEG_RAW_GEOMETRY_WIDTH);
2282 EASValsSetColor( gc.getElementColor() );
2283 EASValsSetEdgeColor( gc.getElementEdgeColor() );
2284 EASValsSetEdgeFlag(true);
2285 EASValsSetLayer(OOFEG_RAW_GEOMETRY_LAYER);
2286 p [ 0 ].x = ( FPNum ) this->giveNode(1)->giveCoordinate(1);
2287 p [ 0 ].y = ( FPNum ) this->giveNode(1)->giveCoordinate(2);
2288 p [ 0 ].z = 0.;
2289 p [ 1 ].x = ( FPNum ) this->giveNode(2)->giveCoordinate(1);
2290 p [ 1 ].y = ( FPNum ) this->giveNode(2)->giveCoordinate(2);
2291 p [ 1 ].z = 0.;
2292 p [ 2 ].x = ( FPNum ) this->giveNode(3)->giveCoordinate(1);
2293 p [ 2 ].y = ( FPNum ) this->giveNode(3)->giveCoordinate(2);
2294 p [ 2 ].z = 0.;
2295
2296 go = CreateTriangle3D(p);
2297 EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | EDGE_COLOR_MASK | EDGE_FLAG_MASK | LAYER_MASK, go);
2298 EGAttachObject(go, ( EObjectP ) this);
2299 EMAddGraphicsToModel(ESIModel(), go);
2300}
2301
2302void TR1_2D_SUPG :: drawScalar(oofegGraphicContext &gc, TimeStep *tStep)
2303{
2304 int i, indx, result = 0;
2305 WCRec p [ 3 ];
2306 GraphicObj *tr;
2307 FloatArray v1, v2, v3;
2308 double s [ 3 ];
2309
2310 if ( !gc.testElementGraphicActivity(this) ) {
2311 return;
2312 }
2313
2314 EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
2315
2316 // if ((gc.giveIntVarMode() == ISM_local) && (gc.giveIntVarType() == IST_VOFFraction)) {
2317 if ( ( gc.giveIntVarType() == IST_VOFFraction ) && ( gc.giveIntVarMode() == ISM_local ) ) {
2318 Polygon matvolpoly;
2319 this->formMaterialVolumePoly(matvolpoly, NULL, temp_normal, temp_p, false);
2320 EASValsSetColor( gc.getStandardSparseProfileColor() );
2321 //GraphicObj *go = matvolpoly.draw(gc,true,OOFEG_VARPLOT_PATTERN_LAYER);
2322 matvolpoly.draw(gc, true, OOFEG_VARPLOT_PATTERN_LAYER);
2323 return;
2324 }
2325
2326 if ( gc.giveIntVarMode() == ISM_recovered ) {
2327 result += this->giveInternalStateAtNode(v1, gc.giveIntVarType(), gc.giveIntVarMode(), 1, tStep);
2328 result += this->giveInternalStateAtNode(v2, gc.giveIntVarType(), gc.giveIntVarMode(), 2, tStep);
2329 result += this->giveInternalStateAtNode(v3, gc.giveIntVarType(), gc.giveIntVarMode(), 3, tStep);
2330 } else if ( gc.giveIntVarMode() == ISM_local ) {
2331 GaussPoint *gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
2332 result += giveIPValue(v1, gp, gc.giveIntVarType(), tStep);
2333 v2 = v1;
2334 v3 = v1;
2335 result *= 3;
2336 }
2337
2338 if ( result != 3 ) {
2339 return;
2340 }
2341
2342 indx = gc.giveIntVarIndx();
2343
2344 s [ 0 ] = v1.at(indx);
2345 s [ 1 ] = v2.at(indx);
2346 s [ 2 ] = v3.at(indx);
2347
2348 EASValsSetLayer(OOFEG_VARPLOT_PATTERN_LAYER);
2349
2350 if ( gc.getScalarAlgo() == SA_ISO_SURF ) {
2351 for ( i = 0; i < 3; i++ ) {
2352 p [ i ].x = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(1);
2353 p [ i ].y = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(2);
2354 p [ i ].z = 0.;
2355 }
2356
2357 //EASValsSetColor(gc.getYieldPlotColor(ratio));
2358 gc.updateFringeTableMinMax(s, 3);
2359 tr = CreateTriangleWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ]);
2360 EGWithMaskChangeAttributes(LAYER_MASK, tr);
2361 EMAddGraphicsToModel(ESIModel(), tr);
2362 } else if ( ( gc.getScalarAlgo() == SA_ZPROFILE ) || ( gc.getScalarAlgo() == SA_COLORZPROFILE ) ) {
2363 double landScale = gc.getLandScale();
2364
2365 for ( i = 0; i < 3; i++ ) {
2366 p [ i ].x = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(1);
2367 p [ i ].y = ( FPNum ) this->giveNode(i + 1)->giveCoordinate(2);
2368 p [ i ].z = s [ i ] * landScale;
2369 }
2370
2371 if ( gc.getScalarAlgo() == SA_ZPROFILE ) {
2372 EASValsSetColor( gc.getDeformedElementColor() );
2373 EASValsSetLineWidth(OOFEG_DEFORMED_GEOMETRY_WIDTH);
2374 EASValsSetFillStyle(FILL_SOLID);
2375 tr = CreateTriangle3D(p);
2376 EGWithMaskChangeAttributes(WIDTH_MASK | COLOR_MASK | FILL_MASK | LAYER_MASK, tr);
2377 } else {
2378 gc.updateFringeTableMinMax(s, 3);
2379 EASValsSetFillStyle(FILL_SOLID);
2380 tr = CreateTriangleWD3D(p, s [ 0 ], s [ 1 ], s [ 2 ]);
2381 EGWithMaskChangeAttributes(FILL_MASK | LAYER_MASK, tr);
2382 }
2383
2384 EMAddGraphicsToModel(ESIModel(), tr);
2385 }
2386}
2387
2388#endif
2389} // end namespace oofem
REGISTER_Element
#define REGISTER_Element(class)
Definition classfactory.h:131
oofem::BoundaryLoad::giveType
bcType giveType() const override
Definition boundaryload.h:166
oofem::BoundaryLoad::giveProperty
virtual double giveProperty(int aProperty, TimeStep *tStep, const std ::map< std ::string, FunctionArgument > &valDict) const
Definition boundaryload.C:132
oofem::CrossSection::give
virtual double give(CrossSectionProperty a, GaussPoint *gp) const
Definition crosssection.C:125
oofem::DofManager::giveCoordinate
double giveCoordinate(int i) const
Definition dofmanager.h:383
oofem::DynamicInputRecord::setField
void setField(int item, InputFieldType id)
Definition dynamicinputrecord.C:302
oofem::Element::giveElementDofIDMask
virtual void giveElementDofIDMask(IntArray &answer) const
Definition element.h:510
oofem::Element::giveNode
Node * giveNode(int i) const
Definition element.h:629
oofem::Element::boundaryLoadArray
IntArray boundaryLoadArray
Definition element.h:147
oofem::Element::dofManArray
IntArray dofManArray
Array containing dofmanager numbers.
Definition element.h:138
oofem::Element::giveBodyLoadArray
IntArray * giveBodyLoadArray()
Returns array containing load numbers of loads acting on element.
Definition element.C:411
oofem::Element::numberOfDofMans
int numberOfDofMans
Number of dofmanagers.
Definition element.h:136
oofem::Element::computeVectorOf
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Definition element.C:103
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::giveCharacteristicValue
virtual double giveCharacteristicValue(CharType type, TimeStep *tStep)
Definition element.C:678
oofem::Element::computeLocalCoordinates
virtual bool computeLocalCoordinates(FloatArray &answer, const FloatArray &gcoords)
Definition element.C:1240
oofem::Element::bodyLoadArray
IntArray bodyLoadArray
Definition element.h:147
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition femcmpnn.h:97
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition femcmpnn.h:79
oofem::FEMComponent::number
int number
Component number.
Definition femcmpnn.h:77
oofem::FEMComponent::giveNumber
int giveNumber() const
Definition femcmpnn.h:104
oofem::FMElement::computeVectorOfPressures
void computeVectorOfPressures(ValueModeType mode, TimeStep *tStep, FloatArray &pressures)
Definition fmelement.C:51
oofem::FMElement::computeVectorOfVelocities
void computeVectorOfVelocities(ValueModeType mode, TimeStep *tStep, FloatArray &velocities)
Definition fmelement.C:44
oofem::FloatArray::computeNorm
double computeNorm() const
Definition floatarray.C:861
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
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::add
void add(const FloatArray &src)
Definition floatarray.C:218
oofem::FloatArray::times
void times(double s)
Definition floatarray.C:834
oofem::FloatMatrix::times
void times(double f)
Definition floatmatrix.C:1155
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
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::zero
void zero()
Zeroes all coefficient of receiver.
Definition floatmatrix.C:1064
oofem::FloatMatrix::computeFrobeniusNorm
double computeFrobeniusNorm() const
Definition floatmatrix.C:1170
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::FluidCrossSection::giveFluidMaterial
FluidDynamicMaterial * giveFluidMaterial()
Definition fluidcrosssection.C:94
oofem::FluidDynamicMaterial::computeDeviatoricStress2D
virtual std::pair< FloatArrayF< 3 >, double > computeDeviatoricStress2D(const FloatArrayF< 3 > &eps, double pressure, GaussPoint *gp, TimeStep *tStep) const
Definition fluiddynamicmaterial.C:52
oofem::FluidModel::giveReynoldsNumber
virtual double giveReynoldsNumber()=0
oofem::GaussPoint::giveNaturalCoordinate
double giveNaturalCoordinate(int i) const
Returns i-th natural element coordinate of receiver.
Definition gausspoint.h:136
oofem::GaussPoint::giveElement
Element * giveElement()
Returns corresponding element to receiver.
Definition gausspoint.h:187
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::Graph::clip
void clip(Polygon &result, const Polygon &a, const Polygon &b)
Definition geotoolbox.C:575
oofem::IntArray::resize
void resize(int n)
Definition intarray.C:73
oofem::IntArray::findFirstIndexOf
int findFirstIndexOf(int value) const
Definition intarray.C:280
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
oofem::IntArray::giveSize
int giveSize() const
Definition intarray.h:211
oofem::LEPlicElementInterface::normal
FloatArray normal
Interface segment normal.
Definition leplic.h:67
oofem::LEPlicElementInterface::setPermanentVolumeFraction
void setPermanentVolumeFraction(double v)
Definition leplic.h:106
oofem::LEPlicElementInterface::p
double p
Line constant of line segment representing interface.
Definition leplic.h:65
oofem::LEPlicElementInterface::vof
double vof
Volume fraction of reference fluid in element.
Definition leplic.h:63
oofem::LEPlic::giveUpdatedXCoordinate
double giveUpdatedXCoordinate(int num)
Definition leplic.h:185
oofem::LEPlic::giveUpdatedCoordinate
void giveUpdatedCoordinate(FloatArray &answer, int num)
Definition leplic.h:179
oofem::LEPlic::giveUpdatedYCoordinate
double giveUpdatedYCoordinate(int num)
Definition leplic.h:186
oofem::LevelSetPCS::giveLevelSetDofManValue
double giveLevelSetDofManValue(int i)
Returns level set value in specific node.
Definition levelsetpcs.h:168
oofem::Load::computeComponentArrayAt
virtual void computeComponentArrayAt(FloatArray &answer, TimeStep *tStep, ValueModeType mode)
Definition load.C:84
oofem::MaterialInterface::giveElementMaterialMixture
virtual void giveElementMaterialMixture(FloatArray &answer, int ielem)=0
oofem::Polygon::addVertex
void addVertex(Vertex v)
Definition geotoolbox.h:96
oofem::Polygon::draw
GraphicObj * draw(oofegGraphicContext &, bool filled, int layer=OOFEG_DEBUG_LAYER)
Definition geotoolbox.C:152
oofem::Polygon::computeVolume
double computeVolume() const
Definition geotoolbox.C:72
oofem::Reinforcement::givePermeability
FloatArray * givePermeability()
Definition reinforcement.h:89
oofem::Reinforcement::givePorosity
double givePorosity()
Accessor.
Definition reinforcement.h:87
oofem::SUPGElement::SUPGElement
SUPGElement(int n, Domain *aDomain)
Definition supgelement.C:61
oofem::SpatialLocalizerInterface::SpatialLocalizerInterface
SpatialLocalizerInterface(Element *element)
Definition spatiallocalizer.h:65
oofem::TR1_2D_SUPG
Definition tr1_2d_supg.h:69
oofem::TR1_2D_SUPG::IPK_TR1_2D_SUPG_vof
static ParamKey IPK_TR1_2D_SUPG_vof
Definition tr1_2d_supg.h:78
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::formMaterialVolumePoly
void formMaterialVolumePoly(Polygon &matvolpoly, LEPlic *matInterface, const FloatArray &normal, const double p, bool updFlag) override
Assembles the true element material polygon (takes receiver vof into accout).
Definition tr1_2d_supg.C:1609
oofem::TR1_2D_SUPG::area
double area
Definition tr1_2d_supg.h:76
oofem::TR1_2D_SUPG::IPK_TR1_2D_SUPG_pvof
static ParamKey IPK_TR1_2D_SUPG_pvof
Definition tr1_2d_supg.h:79
oofem::TR1_2D_SUPG::giveInternalStateAtNode
int giveInternalStateAtNode(FloatArray &answer, InternalStateType type, InternalStateMode mode, int node, TimeStep *tStep) override
Definition tr1_2d_supg.C:2253
oofem::TR1_2D_SUPG::initGeometry
virtual void initGeometry()
Definition tr1_2d_supg.C:1505
oofem::TR1_2D_SUPG::giveIPValue
int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep) override
Definition tr1_2d_supg.C:1918
oofem::TR1_2D_SUPG::IPK_TR1_2D_SUPG_mat1
static ParamKey IPK_TR1_2D_SUPG_mat1
Definition tr1_2d_supg.h:81
oofem::TR1_2D_SUPG::formMyVolumePoly
void formMyVolumePoly(Polygon &myPoly, LEPlic *mat_interface, bool updFlag) override
Assembles receiver volume.
Definition tr1_2d_supg.C:1768
oofem::TR1_2D_SUPG::c
double c[3]
Definition tr1_2d_supg.h:75
oofem::TR1_2D_SUPG::formVolumeInterfacePoly
void formVolumeInterfacePoly(Polygon &matvolpoly, LEPlic *matInterface, const FloatArray &normal, const double p, bool updFlag) override
Assembles receiver material polygon based solely on given interface line.
Definition tr1_2d_supg.C:1641
oofem::TR1_2D_SUPG::IPK_TR1_2D_SUPG_mat0
static ParamKey IPK_TR1_2D_SUPG_mat0
Definition tr1_2d_supg.h:80
oofem::TR1_2D_SUPG::computeMyVolume
double computeMyVolume(LEPlic *matInterface, bool updFlag) override
Computes the volume of receiver.
Definition tr1_2d_supg.C:1791
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::Vertex::setCoords
void setCoords(double x, double y)
Definition geotoolbox.h:77
oofem::ZZNodalRecoveryModelInterface::ZZNodalRecoveryModelInterface
ZZNodalRecoveryModelInterface(Element *element)
Constructor.
Definition zznodalrecoverymodel.h:106
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
YieldStress
#define YieldStress
Definition matconst.h:80
oofem::ContextMode
long ContextMode
Definition contextmode.h:43
oofem::norm
double norm(const FloatArray &x)
Definition floatarray.C:1026
oofem::min
FloatArrayF< N > min(const FloatArrayF< N > &a, const FloatArrayF< N > &b)
Definition floatarrayf.h:522
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::InternalStateMode
InternalStateMode
Determines the mode of internal variable.
Definition internalstatemode.h:40
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::ReinforceBVT
@ ReinforceBVT
Definition bcvaltype.h:49
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::SPRNodalRecoveryModelInterfaceType
@ SPRNodalRecoveryModelInterfaceType
Definition interfacetype.h:51
oofem::ZZNodalRecoveryModelInterfaceType
@ ZZNodalRecoveryModelInterfaceType
Definition interfacetype.h:49
oofem::EIPrimaryFieldInterfaceType
@ EIPrimaryFieldInterfaceType
Definition interfacetype.h:60
oofem::LevelSetPCSElementInterfaceType
@ LevelSetPCSElementInterfaceType
Definition interfacetype.h:75
oofem::LEPlicElementInterfaceType
@ LEPlicElementInterfaceType
Definition interfacetype.h:74
oofem::SpatialLocalizerInterfaceType
@ SpatialLocalizerInterfaceType
Definition interfacetype.h:57
oofem::NodalAveragingRecoveryModelInterfaceType
@ NodalAveragingRecoveryModelInterfaceType
Definition interfacetype.h:50
oofem::TransmissionBC
@ TransmissionBC
Neumann type (prescribed flux).
Definition bctype.h:43
gc
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
OOFEG_DEBUG_LAYER
#define OOFEG_DEBUG_LAYER
Definition oofeggraphiccontext.h:110
OOFEG_VARPLOT_PATTERN_LAYER
#define OOFEG_VARPLOT_PATTERN_LAYER
Definition oofeggraphiccontext.h:105
OOFEG_DEFORMED_GEOMETRY_WIDTH
#define OOFEG_DEFORMED_GEOMETRY_WIDTH
Definition oofeggraphiccontext.h:95
OOFEG_RAW_GEOMETRY_WIDTH
#define OOFEG_RAW_GEOMETRY_WIDTH
Definition oofeggraphiccontext.h:94
OOFEG_RAW_GEOMETRY_LAYER
#define OOFEG_RAW_GEOMETRY_LAYER
Definition oofeggraphiccontext.h:101
PM_UPDATE_PARAMETER_AND_REPORT
#define PM_UPDATE_PARAMETER_AND_REPORT(_val, _pm, _ir, _componentnum, _paramkey, _prio, _flag)
Definition parametermanager.h:66
TRSUPG_ZERO_VOF
#define TRSUPG_ZERO_VOF
Definition tr1_2d_cbs.C:66
tr1_2d_supg.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