OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "supg.h"
36 #include "sparsemtrx.h"
37 #include "nrsolver.h"
38 #include "timestep.h"
39 #include "element.h"
40 #include "dofmanager.h"
41 #include "dof.h"
42 #include "initialcondition.h"
43 #include "maskedprimaryfield.h"
44 #include "verbose.h"
45 #include "supgelement.h"
46 #include "classfactory.h"
47 #include "mathfem.h"
49 #include "leplic.h"
50 #include "levelsetpcs.h"
51 #include "datastream.h"
52 #include "function.h"
53 #include "contextioerr.h"
54 #include "timer.h"
55 #include "unknownnumberingscheme.h"
56 
57 namespace oofem {
58 /* define if implicit interface update required */
59 //#define SUPG_IMPLICIT_INTERFACE
60 
62 
63 
65  VectorAssembler(), lscale(l), dscale(d), uscale(u)
66 {}
67 
69 {
70  SUPGElement *eptr = static_cast< SUPGElement * >( &element );
71  FloatMatrix m1;
72  FloatArray vacc, vtot, ptot, h;
73  IntArray vloc, ploc;
74 
75  int size = eptr->computeNumberOfDofs();
76  eptr->giveLocalVelocityDofMap(vloc);
77  eptr->giveLocalPressureDofMap(ploc);
78  vec.resize(size);
79  vec.zero();
80 
81  eptr->computeVectorOfVelocities(VM_Acceleration, tStep, vacc);
82  eptr->computeVectorOfVelocities(VM_Total, tStep, vtot);
83  eptr->computeVectorOfPressures(VM_Total, tStep, ptot);
84 
85  // MB contributions:
86  // add (M+M_delta)*a
87  eptr->computeAccelerationTerm_MB(m1, tStep);
88  h.beProductOf(m1, vacc);
89  vec.assemble(h, vloc);
90  // add advection terms N+N_delta
91  eptr->computeAdvectionTerm_MB(h, tStep);
92  vec.assemble(h, vloc);
93  // add diffusion terms (K+K_delta)h
94  eptr->computeDiffusionTerm_MB(h, tStep);
95  vec.assemble(h, vloc);
96  // add lsic stabilization term
97  eptr->computeLSICStabilizationTerm_MB(m1, tStep);
98  m1.times( lscale / ( dscale * uscale * uscale ) );
99  h.beProductOf(m1, vtot);
100  vec.assemble(h, vloc);
101  eptr->computeBCLhsTerm_MB(m1, tStep);
102  if ( m1.isNotEmpty() ) {
103  h.beProductOf(m1, vtot);
104  vec.assemble(h, vloc);
105  }
106 
107  // add pressure term
108  eptr->computePressureTerm_MB(m1, tStep);
109  h.beProductOf(m1, ptot); // term due to prescribed pressure
110  vec.assemble(h, vloc);
111  eptr->computeBCLhsPressureTerm_MB(m1, tStep);
112  if ( m1.isNotEmpty() ) {
113  h.beProductOf(m1, ptot);
114  vec.assemble(h, vloc);
115  }
116 
117  // MC contributions:
118  // G^T term - linear advection term
119  eptr->computeLinearAdvectionTerm_MC(m1, tStep);
120  //m1.times(uscale/lscale);
121  m1.times( 1. / ( dscale * uscale ) );
122  eptr->computeVectorOfVelocities(VM_Total, tStep, vtot);
123  h.beProductOf(m1, vtot); // term due to prescribed velocity
124  vec.assemble(h, ploc);
125  // Diffusion term
126  eptr->computeDiffusionTerm_MC(h, tStep);
127  vec.assemble(h, ploc);
128  // AccelerationTerm
129  eptr->computeAccelerationTerm_MC(m1, tStep);
130  h.beProductOf(m1, vacc);
131  vec.assemble(h, ploc);
132  eptr->computeBCLhsPressureTerm_MC(m1, tStep);
133  if ( m1.isNotEmpty() ) {
134  h.beProductOf(m1, vacc);
135  vec.assemble(h, ploc);
136  }
137 
138  // advection N term (nonlinear)
139  eptr->computeAdvectionTerm_MC(h, tStep);
140  vec.assemble(h, ploc);
141  // pressure term
142  eptr->computePressureTerm_MC(m1, tStep);
143  h.beProductOf(m1, ptot);// term due to prescribed pressure
144  vec.assemble(h, ploc);
145 };
146 
147 SUPGTangentAssembler :: SUPGTangentAssembler(MatResponseMode m, double l, double d, double u, double a) :
148  MatrixAssembler(), rmode(m), lscale(l), dscale(d), uscale(u), alpha(a)
149 {}
150 
151 
153 {
154  SUPGElement *element = static_cast< SUPGElement * >( &el );
155  IntArray vloc, ploc;
156  FloatMatrix h;
157  int size = element->computeNumberOfDofs();
158  element->giveLocalVelocityDofMap(vloc);
159  element->giveLocalPressureDofMap(ploc);
160  answer.resize(size, size);
161  answer.zero();
162 
163  element->computeAccelerationTerm_MB(h, tStep);
164  answer.assemble(h, vloc);
165  element->computeAdvectionDerivativeTerm_MB(h, tStep);
166  h.times( alpha * tStep->giveTimeIncrement() );
167  answer.assemble(h, vloc);
168  element->computeDiffusionDerivativeTerm_MB(h, rmode, tStep);
169 
170  h.times( alpha * tStep->giveTimeIncrement() );
171  answer.assemble(h, vloc);
172  element->computePressureTerm_MB(h, tStep);
173  answer.assemble(h, vloc, ploc);
174  element->computeLSICStabilizationTerm_MB(h, tStep);
175  h.times( alpha * tStep->giveTimeIncrement() * lscale / ( dscale * uscale * uscale ) );
176  answer.assemble(h, vloc);
177  element->computeBCLhsTerm_MB(h, tStep);
178  if ( h.isNotEmpty() ) {
179  h.times( alpha * tStep->giveTimeIncrement() );
180  answer.assemble(h, vloc);
181  }
182 
183  element->computeBCLhsPressureTerm_MB(h, tStep);
184  if ( h.isNotEmpty() ) {
185  answer.assemble(h, vloc, ploc);
186  }
187 
188  // conservation eq part
189  element->computeLinearAdvectionTerm_MC(h, tStep);
190  h.times( alpha * tStep->giveTimeIncrement() * 1.0 / ( dscale * uscale ) );
191  answer.assemble(h, ploc, vloc);
192  element->computeAdvectionDerivativeTerm_MC(h, tStep);
193  h.times( alpha * tStep->giveTimeIncrement() );
194  answer.assemble(h, ploc, vloc);
195  element->computeAccelerationTerm_MC(h, tStep);
196  answer.assemble(h, ploc, vloc);
197  element->computeBCLhsPressureTerm_MC(h, tStep);
198  if ( h.isNotEmpty() ) {
199  answer.assemble(h, ploc, vloc);
200  }
201 
202  element->computeDiffusionDerivativeTerm_MC(h, tStep);
203  h.times( alpha * tStep->giveTimeIncrement() );
204  answer.assemble(h, ploc, vloc);
205  element->computePressureTerm_MC(h, tStep);
206  answer.assemble(h, ploc);
207 }
208 
209 
210 SUPG :: SUPG(int i, EngngModel * _master) : FluidModel(i, _master), accelerationVector()
211 {
212  initFlag = 1;
213  ndomains = 1;
214  consistentMassFlag = 0;
215  equationScalingFlag = false;
216  lscale = uscale = dscale = 1.0;
217 }
218 
219 
221 {
222 }
223 
224 
226 {
227  if ( !this->nMethod ) {
228  this->nMethod.reset( classFactory.createSparseLinSolver(solverType, this->giveDomain(1), this) );
229  if ( !this->nMethod ) {
230  OOFEM_ERROR("linear solver creation failed for lstype %d", solverType);
231  }
232  }
233  return this->nMethod.get();
234 }
235 
238 {
239  IRResultType result; // Required by IR_GIVE_FIELD macro
240 
241  result = FluidModel :: initializeFrom(ir);
242  if ( result != IRRT_OK ) {
243  return result;
244  }
245 
247  atolv = 1.e-15;
249 
250 
252 
253  maxiter = 200;
255 
256  int val = 0;
258  solverType = ( LinSystSolverType ) val;
259 
260  val = 0;
262  sparseMtrxType = ( SparseMtrxType ) val;
263 
265  deltaTF = 0;
267 
269 
270  alpha = 0.5;
272 
273  val = 0;
275  equationScalingFlag = val > 0;
276  if ( equationScalingFlag ) {
280  double vref = 1.0; // reference viscosity
281  Re = dscale * uscale * lscale / vref;
282  } else {
283  lscale = uscale = dscale = 1.0;
284  Re = 1.0;
285  }
286 
288  VelocityPressureField.reset( new DofDistributedPrimaryField(this, 1, FT_VelocityPressure, 1) );
289  } else {
290  VelocityPressureField.reset( new PrimaryField(this, 1, FT_VelocityPressure, 1) );
291  }
292 
293  val = 0;
295  if ( val == 1 ) {
296  this->materialInterface.reset( new LEPlic( 1, this->giveDomain(1) ) );
297  this->materialInterface->initializeFrom(ir);
298  // export velocity field
299  FieldManager *fm = this->giveContext()->giveFieldManager();
300  IntArray mask;
301  mask = {V_u, V_v, V_w};
302 
303  std :: shared_ptr< Field > _velocityField( new MaskedPrimaryField ( FT_Velocity, this->VelocityPressureField.get(), mask ) );
304  fm->registerField(_velocityField, FT_Velocity);
305 
306  //fsflag = 0;
307  //IR_GIVE_OPTIONAL_FIELD (ir, fsflag, _IFT_SUPG_fsflag, "fsflag");
308  } else if ( val == 2 ) {
309  // positive coefficient scheme level set alg
310  this->materialInterface.reset( new LevelSetPCS( 1, this->giveDomain(1) ) );
311  this->materialInterface->initializeFrom(ir);
312  }
313 
314  return IRRT_OK;
315 }
316 
317 
318 double
320 // returns unknown quantity like displaacement, velocity of equation eq
321 // This function translates this request to numerical method language
322 {
323  if ( this->requiresUnknownsDictionaryUpdate() ) {
324  int hash = this->giveUnknownDictHashIndx(mode, tStep);
325  if ( dof->giveUnknowns()->includes(hash) ) {
326  return dof->giveUnknowns()->at(hash);
327  } else {
328  OOFEM_ERROR("Dof unknowns dictionary does not contain unknown of value mode (%s)", __ValueModeTypeToString(mode));
329  }
330  } else {
331  int eq = dof->__giveEquationNumber();
332  if ( eq == 0 ) {
333  OOFEM_ERROR("invalid equation number");
334  }
335 
336  if ( mode == VM_Acceleration ) {
337  return accelerationVector.at(eq);
338  /*
339  * if (tStep->isTheCurrentTimeStep()) {
340  * return accelerationVector.at(eq);
341  * } else if (tStep == givePreviousStep()) {
342  * return previousAccelerationVector.at(eq);
343  * } else OOFEM_ERROR("VM_Acceleration history past previous step not supported");
344  */
345  } else if ( mode == VM_Velocity ) {
346  return VelocityPressureField->giveUnknownValue(dof, VM_Total, tStep);
347  } else {
348  return VelocityPressureField->giveUnknownValue(dof, mode, tStep);
349  }
350  }
351 
352  return 0;
353 }
354 
355 
356 void
358 {
359  // update element stabilization
360  for ( auto &elem : d->giveElements() ) {
361  static_cast< FMElement * >( elem.get() )->updateStabilizationCoeffs(tStep);
362  }
363 
364  if ( cmpn == InternalRhs ) {
365  this->internalForces.zero();
367  EModelDefaultEquationNumbering(), d, & this->eNorm);
369  return;
370  } else if ( cmpn == NonLinearLhs ) {
371  this->lhs->zero();
372  //if ( 1 ) { //if ((nite > 5)) // && (rnorm < 1.e4))
373  this->assemble( *lhs, tStep, SUPGTangentAssembler(TangentStiffness, lscale, dscale, uscale, alpha),
375  // } else {
376  // this->assemble( lhs, tStep, SUPGTangentAssembler(SecantStiffness),
377  // EModelDefaultEquationNumbering(), d );
378  // }
379  return;
380  } else {
381  OOFEM_ERROR("Unknown component");
382  }
383 }
384 
385 
386 double
388 {
389  if ( equationScalingFlag ) {
390  return this->Re;
391  } else {
392  return 1.0;
393  }
394 }
395 
396 
397 TimeStep *
399 {
400  if ( master && (!force)) {
402  } else {
403  if ( !stepWhenIcApply ) {
404  double dt = deltaT / this->giveVariableScale(VST_Time);
405 
407  0.0, dt, 0) );
408  }
409 
410  return stepWhenIcApply.get();
411  }
412 }
413 
414 TimeStep *
416 {
417  double dt = deltaT;
418 
419  Domain *domain = this->giveDomain(1);
420 
421  if ( !currentStep ) {
422  // first step -> generate initial step
424  }
425 
426  previousStep = std :: move(currentStep);
427 
428  if ( deltaTF ) {
429  dt *= domain->giveFunction(deltaTF)->evaluateAtTime(previousStep->giveNumber() + 1);
430  }
431 
432  // check for critical time step
433  for ( auto &elem : domain->giveElements() ) {
434  dt = min( dt, static_cast< SUPGElement & >( *elem ).computeCriticalTimeStep(previousStep.get()) );
435  }
436 
437  if ( materialInterface ) {
438  dt = min( dt, materialInterface->computeCriticalTimeStep(previousStep.get()) );
439  }
440 
441  // dt *= 0.6;
442  dt /= this->giveVariableScale(VST_Time);
443 
444  currentStep.reset( new TimeStep(*previousStep, dt) );
445 
446  OOFEM_LOG_INFO( "SolutionStep %d : t = %e, dt = %e\n", currentStep->giveNumber(),
447  currentStep->giveTargetTime() * this->giveVariableScale(VST_Time), dt * this->giveVariableScale(VST_Time) );
448  return currentStep.get();
449 }
450 
451 
452 void
454 {
456  FloatArray *solutionVector = NULL, *prevSolutionVector = NULL;
457  FloatArray externalForces(neq);
458  this->internalForces.resize(neq);
459 
460  if ( tStep->isTheFirstStep() ) {
462  if ( materialInterface ) {
463  materialInterface->initialize();
464  }
465 
466  this->applyIC(stepWhenIcApply);
467  //if (this->fsflag) this->updateDofManActivityMap(tStep);
468  }
469 
471  VelocityPressureField->advanceSolution(tStep);
472  solutionVector = VelocityPressureField->giveSolutionVector(tStep);
473  prevSolutionVector = VelocityPressureField->giveSolutionVector( tStep->givePreviousStep() );
474  }
475 
476  if ( initFlag ) {
478  //previousAccelerationVector.resize(neq);
480  solutionVector->resize(neq);
481  prevSolutionVector->resize(neq);
482  }
483 
485 
487  if ( !lhs ) {
488  OOFEM_ERROR("sparse matrix creation failed");
489  }
490 
491  lhs->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );
492 
493  if ( materialInterface ) {
495  }
496 
497  initFlag = 0;
498  } else if ( requiresUnknownsDictionaryUpdate() ) {
499  // rebuild lhs structure and resize solution vector
501  lhs->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );
502  }
503 
504 
506  *solutionVector = *prevSolutionVector;
507  }
508 
509  //previousAccelerationVector=accelerationVector;
510 
511  // evaluate element supg and sppg stabilization coeffs
513 
514  //
515  // predictor
516  //
517 
520  } else {
521  this->updateSolutionVectors_predictor(* solutionVector, accelerationVector, tStep);
522  }
523 
524  if ( tStep->giveNumber() != 1 ) {
525  if ( materialInterface ) {
526  //if (this->fsflag) updateDofManVals(tStep);
527 #ifdef SUPG_IMPLICIT_INTERFACE
528  #ifdef TIME_REPORT
529  Timer timer;
530  timer.startTimer();
531  #endif
532  materialInterface->updatePosition( this->giveCurrentStep() );
534  #ifdef TIME_REPORT
535  timer.stopTimer();
536  OOFEM_LOG_INFO("SUPG info: user time consumed by updating interfaces: %.2fs\n", timer.getUtime();
537  );
538  #endif
539 #else
540  //updateElementsForNewInterfacePosition (tStep);
541 #endif
542  //if (this->fsflag) this->updateDofManActivityMap(tStep);
543  }
544  }
545 
546  // evaluate element supg and sppg stabilization coeffs
547  // this -> evaluateElementStabilizationCoeffs (tStep);
548 
549  //
550  // assemble rhs (residual)
551  //
552  externalForces.zero();
553  this->assembleVector( externalForces, tStep, ExternalForceAssembler(), VM_Total,
556 
557 #if 0
558  this->updateSolutionVectors(* solutionVector, accelerationVector, incrementalSolutionVector, tStep);
559  this->giveNumericalMethod( this->giveCurrentMetaStep() );
561  double loadLevel;
562  int currentIterations;
563  this->updateComponent( tStep, InternalRhs, this->giveDomain(1) );
564  NM_Status status = this->nMethod->solve(*this->lhs,
565  externalForces,
566  NULL,
567  solutionVector,
569  this->internalForces,
570  this->eNorm,
571  loadLevel, // Only relevant for incrementalBCLoadVector?
573  currentIterations,
574  tStep);
575 
576  if ( !( status & NM_Success ) ) {
577  OOFEM_ERROR("No success in solving problem at time step", tStep->giveNumber());
578  }
579 
580 #else
581  int nite = 0;
582  double _absErrResid, err, avn = 0.0, aivn = 0.0, rnorm;
583  int jj, nnodes = this->giveDomain(1)->giveNumberOfDofManagers();
584  FloatArray rhs, internalForces(neq);
585 
586 
587  // algoritmic rhs part (assembled by e-model (in giveCharComponent service) from various element contribs)
588  internalForces.zero();
589  this->assembleVector( internalForces, tStep, SUPGInternalForceAssembler(lscale, dscale, uscale), VM_Total,
592 
593  rhs.beDifferenceOf(externalForces, internalForces);
594 
595  //
596  // corrector
597  //
598  //OOFEM_LOG_INFO("Iteration IncrSolErr RelResidErr AbsResidErr\n_______________________________________________________\n");
599  OOFEM_LOG_INFO("Iteration IncrSolErr RelResidErr (Rel_MB ,Rel_MC ) AbsResidErr\n__________________________________________________________________________________\n");
600  do {
601  nite++;
602  //
603  // Assemble lhs
604  //
605  //if (nite == 1)
606  {
607  // momentum balance part
608  lhs->zero();
609  if ( 1 ) { //if ((nite > 5)) // && (rnorm < 1.e4))
610  this->assemble( *lhs, tStep, SUPGTangentAssembler(TangentStiffness, lscale, dscale, uscale, alpha),
612  } else {
613  this->assemble( *lhs, tStep, SUPGTangentAssembler(SecantStiffness, lscale, dscale, uscale, alpha),
615  }
616  }
617  //if (this->fsflag) this->imposeAmbientPressureInOuterNodes(lhs,&rhs,tStep);
618 
619 #if 1
620  nMethod->solve(*lhs, rhs, incrementalSolutionVector);
621 #else
622 
623  SparseMtrx *__lhs = lhs->GiveCopy();
624 
625  nMethod->solve(*lhs, rhs, incrementalSolutionVector);
626 
627  // check solver
628  FloatArray __rhs(neq);
629  double __absErr = 0., __relErr = 0.;
630  __lhs->times(incrementalSolutionVector, __rhs);
631  for ( int i = 1; i <= neq; i++ ) {
632  if ( fabs( __rhs.at(i) - rhs.at(i) ) > __absErr ) {
633  __absErr = fabs( __rhs.at(i) - rhs.at(i) );
634  }
635 
636  if ( fabs( ( __rhs.at(i) - rhs.at(i) ) / rhs.at(i) ) > __relErr ) {
637  __relErr = fabs( ( __rhs.at(i) - rhs.at(i) ) / rhs.at(i) );
638  }
639  }
640 
641  OOFEM_LOG_INFO("SUPG: solver check: absoluteError %e, relativeError %e\n", __absErr, __relErr);
642  delete(__lhs);
643 #endif
644 
645 
646 
649  } else {
650  //update
651  this->updateSolutionVectors(* solutionVector, accelerationVector, incrementalSolutionVector, tStep);
654  } // end update
655 
656 #if 0
657  #ifdef SUPG_IMPLICIT_INTERFACE
658  if ( materialInterface ) {
659  #ifdef TIME_REPORT
660  Timer timer;
661  timer.startTimer();
662  #endif
663  //if (this->fsflag) updateDofManVals(tStep);
664  materialInterface->updatePosition( this->giveCurrentStep() );
666  #ifdef TIME_REPORT
667  timer.stopTimer();
668  OOFEM_LOG_INFO( "SUPG info: user time consumed by updating interfaces: %.2fs\n", timer.getUtime() );
669  #endif
670  //if (this->fsflag) this->updateDofManActivityMap(tStep);
671  }
672 
673  #endif
674 #else
675  if ( materialInterface ) {
676  //if (this->fsflag) updateDofManVals(tStep);
677  //if (this->fsflag) this->updateDofManActivityMap(tStep);
678  }
679 
680 #endif
681 
682  // check convergence and repeat iteration if desired
683  if ( avn < 1.e-12 ) {
684  err = aivn;
685  } else {
686  err = sqrt(aivn / avn);
687  }
688 
689  //
690  // assemble rhs (residual)
691  //
692  internalForces.zero();
693  this->assembleVector( internalForces, tStep, SUPGInternalForceAssembler(lscale, dscale, uscale), VM_Total,
696  rhs.beDifferenceOf(externalForces, internalForces);
697 
698  // check convergence and repeat iteration if desired
699  double rnorm_mb = 0.0, rnorm_mc = 0.0;
700  for ( int i = 1; i <= nnodes; i++ ) {
701  DofManager *inode = this->giveDomain(1)->giveDofManager(i);
702  for ( Dof *dof: *inode ) {
703  if ( ( jj = dof->__giveEquationNumber() ) ) {
704  DofIDItem type = dof->giveDofID();
705  double val = rhs.at(jj);
706  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
707  rnorm_mb += val * val;
708  } else {
709  rnorm_mc += val * val;
710  }
711  }
712  }
713  }
714 
715  rnorm_mb = sqrt(rnorm_mb);
716  rnorm_mc = sqrt(rnorm_mc);
717 
718  rnorm = rhs.computeNorm();
719 
720  _absErrResid = 0.0;
721  for ( int i = 1; i <= neq; i++ ) {
722  _absErrResid = max( _absErrResid, fabs( rhs.at(i) ) );
723  }
724 
725  //if ( requiresUnknownsDictionaryUpdate() ) {
726  // OOFEM_LOG_INFO("%-10d n/a %-15e %-15e\n", nite, rnorm, _absErrResid);
727  //} else {
728  OOFEM_LOG_INFO("%-10d %-15e %-15e (%-10e,%-10e) %-15e\n", nite, err, rnorm, rnorm_mb, rnorm_mc, _absErrResid);
729  //}
730 
731  if ( 0 ) {
732  // evaluate element supg and sppg stabilization coeffs
734  //
735  // assemble rhs (residual)
736  //
737  internalForces.zero();
738  this->assembleVector( internalForces, tStep, SUPGInternalForceAssembler(lscale, dscale, uscale), VM_Total,
741  rhs.beDifferenceOf(externalForces, internalForces);
742  }
743  } while ( ( rnorm > rtolv ) && ( _absErrResid > atolv ) && ( nite <= maxiter ) );
744 
745  if ( nite <= maxiter ) {
746  OOFEM_LOG_INFO("SUPG info: number of iterations: %d\n", nite);
747  } else {
748  OOFEM_WARNING("Convergence not reached, number of iterations: %d\n", nite);
749  if ( stopmaxiter ) {
750  exit(1);
751  }
752  }
753 #endif
754 
755 #ifndef SUPG_IMPLICIT_INTERFACE
756  if ( materialInterface ) {
757  #ifdef TIME_REPORT
758  Timer timer;
759  timer.startTimer();
760  #endif
761  materialInterface->updatePosition( this->giveCurrentStep() );
763  #ifdef TIME_REPORT
764  timer.stopTimer();
765  OOFEM_LOG_INFO( "SUPG info: user time consumed by updating interfaces: %.2fs\n", timer.getUtime() );
766  #endif
767  //if (this->fsflag) this->updateDofManActivityMap(tStep);
768  }
769 
770 #endif
771 
772  // update solution state counter
773  tStep->incrementStateCounter();
774 }
775 
776 
777 void
779 {
780  this->updateInternalState(tStep);
782  if ( materialInterface ) {
783  materialInterface->updateYourself(tStep);
784  }
785 
786  //previousSolutionVector = solutionVector;
787 }
788 
789 
790 void
792 {
793  for ( auto &domain: domainList ) {
795  for ( auto &dman : domain->giveDofManagers() ) {
796  this->updateDofUnknownsDictionary(dman.get(), tStep);
797  }
798  }
799 
800  for ( auto &elem : domain->giveElements() ) {
801  elem->updateInternalState(tStep);
802  }
803  }
804 }
805 
806 
809 {
810  contextIOResultType iores;
811 
812  if ( ( iores = EngngModel :: saveContext(stream, mode) ) != CIO_OK ) {
813  THROW_CIOERR(iores);
814  }
815 
816  if ( ( iores = VelocityPressureField->saveContext(stream, mode) ) != CIO_OK ) {
817  THROW_CIOERR(iores);
818  }
819 
820  if ( ( iores = accelerationVector.storeYourself(stream) ) != CIO_OK ) {
821  THROW_CIOERR(iores);
822  }
823 
824  if ( materialInterface ) {
825  if ( ( iores = materialInterface->saveContext(stream, mode) ) != CIO_OK ) {
826  THROW_CIOERR(iores);
827  }
828  }
829 
830  return CIO_OK;
831 }
832 
833 
836 {
837  contextIOResultType iores;
838 
839  if ( ( iores = EngngModel :: restoreContext(stream, mode) ) != CIO_OK ) {
840  THROW_CIOERR(iores);
841  }
842 
843  if ( ( iores = VelocityPressureField->restoreContext(stream, mode) ) != CIO_OK ) {
844  THROW_CIOERR(iores);
845  }
846 
847  if ( ( iores = accelerationVector.restoreYourself(stream) ) != CIO_OK ) {
848  THROW_CIOERR(iores);
849  }
850 
851  if ( materialInterface ) {
852  if ( ( iores = materialInterface->restoreContext(stream, mode) ) != CIO_OK ) {
853  THROW_CIOERR(iores);
854  }
855  }
856 
857  return CIO_OK;
858 }
859 
860 
861 int
863 {
864  // check internal consistency
865  // if success returns nonzero
866  Domain *domain = this->giveDomain(1);
867 
868  // check for proper element type
869  for ( auto &elem : domain->giveElements() ) {
870  if ( dynamic_cast< SUPGElement * >( elem.get() ) == NULL ) {
871  OOFEM_WARNING("Element %d has no SUPG base", elem->giveLabel());
872  return 0;
873  }
874  }
875 
876  int ret = EngngModel :: checkConsistency();
877  if ( ret == 0 ) {
878  return 0;
879  }
880 
881 
882  // scale boundary and initial conditions
883  if ( equationScalingFlag ) {
884  for ( auto &bc : domain->giveBcs() ) {
885  if ( bc->giveBCValType() == VelocityBVT ) {
886  bc->scale(1. / uscale);
887  } else if ( bc->giveBCValType() == PressureBVT ) {
888  bc->scale( 1. / this->giveVariableScale(VST_Pressure) );
889  } else if ( bc->giveBCValType() == ForceLoadBVT ) {
890  bc->scale( 1. / this->giveVariableScale(VST_Force) );
891  } else {
892  OOFEM_ERROR("unknown bc/ic type");
893  }
894  }
895 
896  for ( auto &ic : domain->giveIcs() ) {
897  if ( ic->giveICValType() == VelocityBVT ) {
898  ic->scale(VM_Total, 1. / uscale);
899  } else if ( ic->giveICValType() == PressureBVT ) {
900  ic->scale( VM_Total, 1. / this->giveVariableScale(VST_Pressure) );
901  } else {
902  OOFEM_ERROR("unknown bc/ic type");
903  }
904  }
905  }
906 
907  return 1;
908 }
909 
910 
911 void
913 {
916  //this->giveNumericalMethod( this->giveCurrentMetaStep() )->setDomain( this->giveDomain(1) );
917 }
918 
919 void
920 SUPG :: printDofOutputAt(FILE *stream, Dof *iDof, TimeStep *tStep)
921 {
922  double pscale = ( dscale * uscale * uscale );
923 
924  DofIDItem type = iDof->giveDofID();
925  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
926  iDof->printSingleOutputAt(stream, tStep, 'd', VM_Total, uscale);
927  } else if ( type == P_f ) {
928  iDof->printSingleOutputAt(stream, tStep, 'd', VM_Total, pscale);
929  } else {
930  OOFEM_ERROR("unsupported dof type");
931  }
932 }
933 
934 void
936 {
937  Domain *domain = this->giveDomain(1);
939  FloatArray *vp_vector = NULL;
940 
941 #ifdef VERBOSE
942  OOFEM_LOG_INFO("Applying initial conditions\n");
943 #endif
944  int nman = domain->giveNumberOfDofManagers();
945 
948 
950  VelocityPressureField->advanceSolution(stepWhenIcApply);
951  vp_vector = VelocityPressureField->giveSolutionVector(stepWhenIcApply);
952  vp_vector->resize(neq);
953  vp_vector->zero();
954 
955  for ( int j = 1; j <= nman; j++ ) {
956  DofManager *node = domain->giveDofManager(j);
957 
958  for ( Dof *dof: *node ) {
959  // ask for initial values obtained from
960  // bc (boundary conditions) and ic (initial conditions)
961  if ( !dof->isPrimaryDof() ) {
962  continue;
963  }
964 
965  int jj = dof->__giveEquationNumber();
966  DofIDItem type = dof->giveDofID();
967  if ( jj ) {
968  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
969  vp_vector->at(jj) = dof->giveUnknown(VM_Total, stepWhenIcApply);
970  } else {
971  vp_vector->at(jj) = dof->giveUnknown(VM_Total, stepWhenIcApply);
972  }
973  }
974  }
975  }
976  }
977 
978  //ICs and BCs are scaled already in CheckConsistency
979  //vp_vector->times(1./this->giveVariableScale(VST_Velocity));
980 
981  // update element state according to given ic
982 
983  //this->initElementsForNewStep (stepWhenIcApply);
984  if ( materialInterface ) {
985  this->updateElementsForNewInterfacePosition(stepWhenIcApply);
986  }
987 
988  for ( auto &elem : domain->giveElements() ) {
989  SUPGElement *element = static_cast< SUPGElement * >( elem.get() );
990  element->updateInternalState(stepWhenIcApply);
991  element->updateYourself(stepWhenIcApply);
992  }
993 }
994 
995 double
997 {
998  if ( varID == VST_Length ) {
999  return this->lscale;
1000  } else if ( varID == VST_Velocity ) {
1001  return this->uscale;
1002  } else if ( varID == VST_Density ) {
1003  return this->dscale;
1004  } else if ( varID == VST_Time ) {
1005  return ( lscale / uscale );
1006  } else if ( varID == VST_Pressure ) {
1007  return ( dscale * uscale * uscale );
1008  } else if ( varID == VST_Force ) {
1009  return ( uscale * uscale / lscale );
1010  } else if ( varID == VST_Viscosity ) {
1011  return 1.0;
1012  } else {
1013  OOFEM_ERROR("unknown variable type");
1014  }
1015 
1016  return 0.0;
1017 }
1018 
1019 void
1021 {
1022  Domain *domain = this->giveDomain(1);
1023 
1024  for ( auto &elem : domain->giveElements() ) {
1025  FMElement *ePtr = static_cast< FMElement * >( elem.get() );
1026  ePtr->updateStabilizationCoeffs(tStep);
1027  }
1028 }
1029 
1030 void
1032 {
1033  Domain *domain = this->giveDomain(1);
1034 
1035  OOFEM_LOG_DEBUG("SUPG :: updateElements - updating elements for interface position");
1036 
1037 
1038  for ( auto &elem : domain->giveElements() ) {
1039  SUPGElement *ePtr = static_cast< SUPGElement * >( elem.get() );
1041  }
1042 }
1043 
1044 
1045 void
1047 {
1048  // update DOF unknowns dictionary, where
1049  // unknowns are hold instead of keeping them in global unknowns
1050  // vectors in engng instancies
1051  // this is necessary, because during solution equation numbers for
1052  // particular DOFs may changed, and it is necesary to keep them
1053  // in DOF level.
1054 
1055  // empty -> all done in updateDofUnknownsDictionary_corrector (TimeStep* tStep);
1056 }
1057 
1058 
1059 void
1061 {
1062  double deltaT = tStep->giveTimeIncrement();
1063  Domain *domain = this->giveDomain(1);
1064 
1065  int nnodes = domain->giveNumberOfDofManagers();
1066 
1068  for ( int j = 1; j <= nnodes; j++ ) {
1069  DofManager *inode = domain->giveDofManager(j);
1070  for ( Dof *dof: *inode ) {
1071  double val, prev_val, accel;
1072  DofIDItem type = dof->giveDofID();
1073  if ( !dof->hasBc(tStep) ) {
1074  prev_val = dof->giveUnknown( VM_Total, tStep->givePreviousStep() );
1075  accel = dof->giveUnknown( VM_Acceleration, tStep->givePreviousStep() );
1076  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
1077  val = prev_val + deltaT * accel;
1078  } else {
1079  val = prev_val;
1080  }
1081 
1082  dof->updateUnknownsDictionary(tStep, VM_Total, val); // velocity
1083  dof->updateUnknownsDictionary(tStep, VM_Acceleration, accel); // acceleration
1084  } else {
1085  val = dof->giveBcValue(VM_Total, tStep);
1086  dof->updateUnknownsDictionary(tStep, VM_Total, val); // velocity
1087  dof->updateUnknownsDictionary(tStep, VM_Acceleration, 0.0); // acceleration
1088  }
1089  }
1090  }
1091  }
1092 }
1093 
1094 
1095 void
1097 {
1098  double deltaT = tStep->giveTimeIncrement();
1099  Domain *domain = this->giveDomain(1);
1100  int nnodes = domain->giveNumberOfDofManagers();
1101 
1103  for ( int j = 1; j <= nnodes; j++ ) {
1104  DofManager *inode = domain->giveDofManager(j);
1105  for ( Dof *iDof: *inode ) {
1106  DofIDItem type = iDof->giveDofID();
1107  if ( !iDof->hasBc(tStep) ) {
1108  double val, prev_val;
1109  prev_val = iDof->giveUnknown(VM_Total, tStep);
1110  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) {
1111  val = prev_val + deltaT *alpha *incrementalSolutionVector.at( iDof->__giveEquationNumber() );
1112  } else {
1113  val = prev_val + incrementalSolutionVector.at( iDof->__giveEquationNumber() );
1114  }
1115 
1116  iDof->updateUnknownsDictionary(tStep, VM_Total, val); // velocity
1117 
1118  prev_val = iDof->giveUnknown(VM_Acceleration, tStep);
1119  val = prev_val + incrementalSolutionVector.at( iDof->__giveEquationNumber() );
1120  iDof->updateUnknownsDictionary(tStep, VM_Acceleration, val); // acceleration
1121  }
1122  }
1123  }
1124  }
1125 }
1126 
1127 int
1129 {
1130  if ( ( tStep == this->giveCurrentStep() ) || ( tStep == this->givePreviousStep() ) ) {
1131  return ( tStep->giveNumber() % 2 ) * 100 + mode;
1132  } else {
1133  OOFEM_ERROR("unsupported solution step");
1134  }
1135 
1136  return 0;
1137 }
1138 
1139 void
1141 {
1142  double deltaT = tStep->giveTimeIncrement();
1143  Domain *domain = this->giveDomain(1);
1144 
1145  for ( auto &node : domain->giveDofManagers() ) {
1146  for ( Dof *iDof: *node ) {
1147  if ( !iDof->isPrimaryDof() ) {
1148  continue;
1149  }
1150 
1151  int jj = iDof->__giveEquationNumber();
1152  DofIDItem type = iDof->giveDofID();
1153 
1154  if ( jj ) {
1155  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) { // v = v + (deltaT)*a
1156  solutionVector.at(jj) += deltaT * accelerationVector.at(jj);
1157  }
1158  }
1159  }
1160  }
1161 
1162  for ( auto &elem : domain->giveElements() ) {
1163  int ndofman = elem->giveNumberOfInternalDofManagers();
1164  for ( int ji = 1; ji <= ndofman; ji++ ) {
1165  DofManager *dofman = elem->giveInternalDofManager(ji);
1166  for ( Dof *iDof: *dofman ) {
1167  if ( !iDof->isPrimaryDof() ) {
1168  continue;
1169  }
1170 
1171  int jj = iDof->__giveEquationNumber();
1172  DofIDItem type = iDof->giveDofID();
1173 
1174  if ( jj ) {
1175  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) { // v = v + (deltaT*alpha)*da
1176  solutionVector.at(jj) += deltaT * accelerationVector.at(jj);
1177  }
1178  }
1179  }
1180  } // end loop over elem internal dofmans
1181  } // end loop over elems
1182 }
1183 
1184 
1185 void
1187 {
1188  double deltaT = tStep->giveTimeIncrement();
1189  Domain *domain = this->giveDomain(1);
1190 
1191  accelerationVector.add(incrementalSolutionVector);
1192 
1193  for ( auto &node : domain->giveDofManagers() ) {
1194 
1195  for ( Dof *iDof: *node ) {
1196  if ( !iDof->isPrimaryDof() ) {
1197  continue;
1198  }
1199 
1200  int jj = iDof->__giveEquationNumber();
1201  DofIDItem type = iDof->giveDofID();
1202 
1203  if ( jj ) {
1204  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) { // v = v + (deltaT*alpha)*da
1205  solutionVector.at(jj) += deltaT * alpha * incrementalSolutionVector.at(jj);
1206  } else {
1207  solutionVector.at(jj) += incrementalSolutionVector.at(jj); // p = p + dp
1208  }
1209  }
1210  }
1211  } // end loop over dnam
1212 
1213  for ( auto &elem : domain->giveElements() ) {
1214  int ndofman = elem->giveNumberOfInternalDofManagers();
1215  for ( int ji = 1; ji <= ndofman; ji++ ) {
1216  DofManager *dofman = elem->giveInternalDofManager(ji);
1217  for ( Dof *iDof: *dofman ) {
1218  if ( !iDof->isPrimaryDof() ) {
1219  continue;
1220  }
1221 
1222  int jj = iDof->__giveEquationNumber();
1223  DofIDItem type = iDof->giveDofID();
1224 
1225  if ( jj ) {
1226  if ( ( type == V_u ) || ( type == V_v ) || ( type == V_w ) ) { // v = v + (deltaT*alpha)*da
1227  solutionVector.at(jj) += deltaT * alpha * incrementalSolutionVector.at(jj);
1228  } else {
1229  solutionVector.at(jj) += incrementalSolutionVector.at(jj); // p = p + dp
1230  }
1231  }
1232  }
1233  } // end loop over elem internal dofmans
1234  } // end loop over elems
1235 }
1236 
1237 
1238 #define __VOF_TRESHOLD 1.e-5
1239 
1240 #define __NODE_OUT 0
1241 #define __NODE_IN 1
1242 #define __NODE_OUT_ACTIVE 2
1243 #define __NODE_INTERPOL_CANDIDATE 3
1244 
1245 /*
1246  * void
1247  * SUPG::updateDofManActivityMap (TimeStep* tStep)
1248  * {
1249  * // loop over inactive nodes and test if some of its shared element has become active (temp_vof>0)
1250  * // then initialize its velocity and pressure
1251  *
1252  * Domain* domain = this->giveDomain(1);
1253  * Element* ielem;
1254  * int nnodes = domain->giveNumberOfDofManagers();
1255  * int nelems = domain->giveNumberOfElements();
1256  * int i, inode, ie;
1257  * //int inodes, _code;
1258  * LEPlicElementInterface *interface;
1259  * IntArray dofMask(2);
1260  * FloatArray vv;
1261  *
1262  * dofMask.at(1) = V_u;
1263  * dofMask.at(2) = V_v;
1264  *
1265  * printf ("SUPG::updateDofManActivityMap called\n");
1266  *
1267  * IntArray __old_DofManActivityMask (__DofManActivityMask);
1268  * __DofManActivityMask.resize(nnodes);
1269  * for (i=1; i<=nnodes;i++) __DofManActivityMask.at(i) = __NODE_OUT;
1270  *
1271  *
1272  * for (ie=1; ie<= nelems; ie++) {
1273  * ielem = domain->giveElement(ie);
1274  * if ((interface = (LEPlicElementInterface*) (ielem->giveInterface(LEPlicElementInterfaceType)))) {
1275  * if (interface->giveVolumeFraction() > 0.9999) {
1276  * // mark all its nodes as active
1277  * nnodes=ielem->giveNumberOfNodes();
1278  * for (i=1; i<=nnodes;i++) this->__DofManActivityMask.at(ielem->giveNode(i)->giveNumber())=__NODE_IN;
1279  * } else if (interface->giveVolumeFraction() > 0.0) {
1280  * nnodes=ielem->giveNumberOfNodes();
1281  * for (i=1; i<=nnodes;i++) {
1282  * inode = ielem->giveNode(i)->giveNumber();
1283  * //if (__DofManActivityMask.at(inode) == __NODE_OUT) {
1284  * FloatArray normal; double p,x,y;
1285  * interface->giveTempInterfaceNormal(normal);
1286  * p = interface->giveTempLineConstant();
1287  *
1288  * x=domain->giveNode(inode)->giveCoordinate(1);
1289  * y=domain->giveNode(inode)->giveCoordinate(2);
1290  *
1291  * if ((normal.at(1)*x+normal.at(2)*y+p) >= 0.) {
1292  * __DofManActivityMask.at(inode) = __NODE_IN;
1293  * } else {
1294  * if (__DofManActivityMask.at(inode) == __NODE_OUT) {
1295  * if (interface->giveVolumeFraction() < 0.5) __DofManActivityMask.at(inode) = __NODE_INTERPOL_CANDIDATE;
1296  * else __DofManActivityMask.at(inode) = __NODE_OUT_ACTIVE;
1297  * }
1298  * }
1299  * //}
1300  * }
1301  * }
1302  * }
1303  * }
1304  *
1305  * }
1306  *
1307  *
1308  * void
1309  * SUPG:: updateDofManVals (TimeStep* tStep)
1310  * {
1311  * Domain* domain = this->giveDomain(1);
1312  * int nnodes = domain->giveNumberOfDofManagers();
1313  * int inode;
1314  * const IntArray* ct;
1315  * IntArray dofMask(2);
1316  * FloatArray vv;
1317  * Dof* iDof;
1318  *
1319  * dofMask.at(1) = V_u;
1320  * dofMask.at(2) = V_v;
1321  *
1322  * printf ("SUPG::updateDofManVals called\n");
1323  * for (inode=1; inode <= nnodes; inode++) {
1324  * if (__DofManActivityMask.at(inode) == __NODE_INTERPOL_CANDIDATE) {
1325  * // loop over shared elements
1326  * ct = domain->giveConnectivityTable()->giveDofManConnectivityArray(inode);
1327  * // now we have to find active dofmans in the neighborhood to interpolate velocity
1328  * double vx=0.,vy=0.;
1329  * int _j, _k, _d, _cnt=0;
1330  * Element* _je;
1331  * for (_j=1; _j<=ct->giveSize(); _j++) {
1332  * _je=domain->giveElement(ct->at(_j));
1333  * for (_k=1; _k<=_je->giveNumberOfNodes(); _k++) {
1334  * int _code = __DofManActivityMask.at(_je->giveNode(_k)->giveNumber());
1335  * if ((_code == __NODE_IN) || (_code == __NODE_OUT_ACTIVE)) { // active node (not candidate)
1336  * _je->giveNode(_k)->giveUnknownVector(vv, dofMask, VM_Total, tStep);
1337  * vx += vv.at(1); vy += vv.at(2);
1338  * _cnt++;
1339  * }
1340  * }
1341  * }
1342  * if (_cnt) printf ("%d ", inode);
1343  * vx /= _cnt; vy /= _cnt; // average
1344  * for (Dof *iDof: domain->giveDofManager(inode)) {
1345  * DofIDItem type = iDof->giveDofID();
1346  * if (!iDof->hasBc(tStep)) {
1347  * if (type == V_u)
1348  * iDof->updateUnknownsDictionary (tStep, VM_Total, vx); // velocity
1349  * else if (type == V_v)
1350  * iDof->updateUnknownsDictionary (tStep, VM_Total, vy); // velocity
1351  * else if (type == V_w) {
1352  * OOFEM_ERROR("3d not yet supported");
1353  * } else if (type == P_f) {
1354  * //iDof->updateUnknownsDictionary (tStep, VM_Total, 0.0);
1355  * } else {_error2 ("unknown DOF type encountered (%s)", __DofIDItemToString (type));}
1356  * }
1357  * }
1358  * }
1359  * }
1360  * }
1361  *
1362  * void
1363  * SUPG::imposeAmbientPressureInOuterNodes(SparseMtrx* lhs, FloatArray* rhs, TimeStep* tStep)
1364  * {
1365  * Domain* domain = this->giveDomain(1);
1366  * int nnodes = domain->giveNumberOfDofManagers();
1367  * int inode, _code;
1368  * auto pdof;
1369  *
1370  * IntArray inMask(nnodes); inMask.zero();
1371  * for (inode=1;inode<=nnodes;inode++) {
1372  * _code = __DofManActivityMask.at(inode);
1373  * if ((_code == __NODE_OUT_ACTIVE) || (_code == __NODE_INTERPOL_CANDIDATE) || (_code == __NODE_OUT)) {
1374  * if ((pdof=domain->giveDofManager(inode)->findDofWithDofId(P_f)) != domain->giveDofManager(inode)->end() ) {
1375  * int _eq = (*pdof)->giveEquationNumber();
1376  * if (_eq) {
1377  * if (fabs(lhs->at(_eq,_eq)) > 0.0) {
1378  * lhs->at(_eq,_eq) *= 1.e6;
1379  * rhs->at(_eq) = 0.0;
1380  * } else {
1381  * lhs->at(_eq,_eq) = 1.0;
1382  * rhs->at(_eq) = 0.0;
1383  * }
1384  * }
1385  * }
1386  * }
1387  * }
1388  * }
1389  *
1390  * void
1391  * SUPG::__debug (TimeStep* tStep)
1392  * {
1393  * int i, in, id, nincr = 1000;
1394  * int neq = this -> giveNumberOfDomainEquations (1, EModelDefaultEquationNumbering());
1395  * IntArray loc;
1396  * SUPGElement* element = (SUPGElement*) giveDomain(1)->giveElement(1);
1397  * FloatArray vincr(6), F(6), fi(6), fprev(6), fapprox(6), *solutionVector;
1398  * FloatMatrix d(6,6);
1399  *
1400  * vincr.at(1) = 1.e-2;
1401  * vincr.at(2) = 1.e-7;
1402  * vincr.at(3) = 2.e-2;
1403  * element -> giveLocationArray (loc);
1404  * VelocityPressureField->advanceSolution(tStep);
1405  * solutionVector = VelocityPressureField->giveSolutionVector(tStep);
1406  * solutionVector->resize(neq);
1407  * // restrict vincr
1408  * for (i=1; i<=6; i++)
1409  * if (loc.at(i) <= 0) vincr.at(i) = 0.0;
1410  *
1411  *
1412  * for (in=1; in<=nincr; in++) {
1413  * for (i=1; i<=6; i++)
1414  * if ((id = loc.at(i))) solutionVector->at(id) += vincr.at(i);
1415  * // compute F from derivative
1416  * EngngModel::giveElementCharacteristicMatrix(d, 1, TangentDiffusionDerivativeTerm_MB, tStep, giveDomain(1));
1417  * fi.beProductOf (d, vincr);
1418  * fapprox = fprev;
1419  * fapprox.add(fi);
1420  * // compute true F
1421  * this->giveElementCharacteristicVector(F, 1, DiffusionTerm_MB, VM_Total, tStep, giveDomain(1));
1422  * // print
1423  * printf ("%d %e %e %e %e %e %e %e %e %e %e %e %e\n", in, F.at(1), F.at(2), F.at(3), F.at(4), F.at(5), F.at(6),
1424  * fapprox.at(1), fapprox.at(2), fapprox.at(3), fapprox.at(4), fapprox.at(5), fapprox.at(6));
1425  * // if (in==1) fapprox = F;
1426  * fprev = fapprox;
1427  * }
1428  *
1429  * }
1430  *
1431  */
1432 } // end namespace oofem
EngngModelTimer timer
E-model timer.
Definition: engngm.h:267
virtual void updateDofUnknownsDictionary(DofManager *dman, TimeStep *tStep)
Updates necessary values in Dofs unknown dictionaries.
Definition: supg.C:1046
LinSystSolverType
The values of this type should be related not to specific solvers, but more to specific packages that...
The representation of EngngModel default unknown numbering.
#define NM_Success
Numerical method exited with success.
Definition: nmstatus.h:47
virtual double giveUnknownComponent(ValueModeType mode, TimeStep *tStep, Domain *d, Dof *dof)
Returns requested unknown.
Definition: supg.C:319
virtual void printDofOutputAt(FILE *stream, Dof *iDof, TimeStep *tStep)
DOF printing routine.
Definition: supg.C:920
std::unique_ptr< TimeStep > currentStep
Current time step.
Definition: engngm.h:231
virtual ~SUPG()
Definition: supg.C:220
virtual void computePressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)=0
Computes pressure terms for momentum balance equations(s).
double lscale
Definition: supg.h:142
virtual void computeAdvectionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes the derivative of advection terms for mass conservation equation with respect to nodal veloc...
void initMetaStepAttributes(MetaStep *mStep)
Update e-model attributes attributes according to step metaStep attributes.
Definition: engngm.C:580
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode)
Stores the state of model to output stream.
Definition: supg.C:808
Class and object Domain.
Definition: domain.h:115
void updateElementsForNewInterfacePosition(TimeStep *tStep)
Definition: supg.C:1031
virtual int computeNumberOfDofs()
Computes or simply returns total number of element&#39;s local DOFs.
Definition: element.h:424
virtual int giveNumberOfDomainEquations(int di, const UnknownNumberingScheme &num)
Returns number of equations for given domain in active (current time step) time step.
Definition: engngm.C:391
Abstract class representing subset of DOFs (identified by DofId mask) of primary field.
virtual TimeStep * giveSolutionStepWhenIcApply(bool force=false)
Returns the solution step when Initial Conditions (IC) apply.
Definition: supg.C:398
virtual TimeStep * giveNextStep()
Returns next time step (next to current step) of receiver.
Definition: supg.C:415
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description in input reader.
Definition: engngm.C:262
virtual void updateDomainLinks()
Updates domain links after the domains of receiver have changed.
Definition: supg.C:912
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition: domain.h:432
Base class for all matrices stored in sparse format.
Definition: sparsemtrx.h:60
SparseMtrxType sparseMtrxType
Definition: supg.h:111
virtual void computeLSICStabilizationTerm_MB(FloatMatrix &answer, TimeStep *tStep)=0
Computes SLIC stabilization term for momentum balance equation(s).
LinSystSolverType solverType
Definition: supg.h:110
Class representing meta step.
Definition: metastep.h:62
SparseMtrx * createSparseMtrx(SparseMtrxType type)
Creates new instance of sparse matrix corresponding to given keyword.
Definition: classfactory.C:105
Abstract class representing field of primary variables (those, which are unknown and are typically as...
Definition: primaryfield.h:104
virtual void computeLinearAdvectionTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes the linear advection term for mass conservation equation.
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatarray.C:872
FloatArray internalForces
Definition: supg.h:119
std::vector< std::unique_ptr< Domain > > domainList
List of problem domains.
Definition: engngm.h:207
std::unique_ptr< TimeStep > previousStep
Previous time step.
Definition: engngm.h:233
virtual void updateInternalState(TimeStep *tStep)
Updates element state after equilibrium in time step has been reached.
Definition: supgelement.C:351
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
std::unique_ptr< SparseMtrx > lhs
Definition: supg.h:113
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
virtual void updateYourself(TimeStep *tStep)
Updates internal state after finishing time step.
Definition: supg.C:778
virtual int requiresUnknownsDictionaryUpdate()
Indicates if EngngModel requires Dofs dictionaries to be updated.
Definition: supg.h:184
ValueModeType
Type representing the mode of UnknownType or CharType, or similar types.
Definition: valuemodetype.h:78
VarScaleType
Type determining the scale corresponding to particular variable.
Definition: varscaletype.h:40
int consistentMassFlag
Definition: supg.h:138
Base class for fluid problems.
Definition: fluidmodel.h:46
std::unique_ptr< SparseLinearSystemNM > nMethod
Numerical method used to solve the problem.
Definition: supg.h:108
This base class is an abstraction for numerical algorithm.
Definition: nummet.h:80
virtual void updateElementForNewInterfacePosition(TimeStep *tStep)
Definition: supgelement.h:82
std::unique_ptr< TimeStep > stepWhenIcApply
Solution step when IC (initial conditions) apply.
Definition: engngm.h:229
bool equationScalingFlag
Definition: supg.h:140
virtual void updateYourself(TimeStep *tStep)
Updates element state after equilibrium in time step has been reached.
Definition: element.C:779
virtual bool hasField(InputFieldType id)=0
Returns true if record contains field identified by idString keyword.
#define _IFT_SUPG_deltatFunction
Definition: supg.h:51
void incrementStateCounter()
Updates solution state counter.
Definition: timestep.h:190
Abstract base class for all finite elements.
Definition: element.h:145
double dscale
Definition: supg.h:142
virtual int checkConsistency()
Allows programmer to test some receiver&#39;s internal data, before computation begins.
Definition: supg.C:862
Base class for dof managers.
Definition: dofmanager.h:113
std::vector< std::unique_ptr< DofManager > > & giveDofManagers()
Definition: domain.h:400
virtual void giveLocalPressureDofMap(IntArray &map)
Definition: supgelement.h:209
virtual int __giveEquationNumber() const =0
Returns equation number of receiver, usually assigned by emodel.
virtual void updateComponent(TimeStep *tStep, NumericalCmpn cmpn, Domain *d)
Updates components mapped to numerical method if necessary during solution process.
Definition: supg.C:357
General stabilized SUPG/PSPG element for CFD analysis.
Definition: supgelement.h:58
FieldManager * giveFieldManager()
Definition: engngm.h:131
SUPGTangentAssembler(MatResponseMode m, double l, double d, double u, double a)
Definition: supg.C:147
REGISTER_EngngModel(ProblemSequence)
virtual void solveYourselfAt(TimeStep *tStep)
Solves problem for given time step.
Definition: supg.C:453
unsigned long NM_Status
Mask defining NumMetod Status; which can be asked after finishing computation by Numerical Method...
Definition: nmstatus.h:44
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
Class implementing an array of integers.
Definition: intarray.h:61
MatResponseMode
Describes the character of characteristic material matrix.
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
virtual void giveUnknowns(FloatArray &masterUnknowns, ValueModeType mode, TimeStep *tStep)
The key method of class Dof.
Definition: dof.C:162
virtual void updateStabilizationCoeffs(TimeStep *tStep)
Updates the stabilization coefficients used for CBS and SUPG algorithms.
Definition: fmelement.h:64
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description in input reader.
Definition: supg.C:237
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
double giveTimeIncrement()
Returns solution step associated time increment.
Definition: timestep.h:150
int maxiter
Max number of iterations.
Definition: supg.h:128
bool isTheFirstStep()
Check if receiver is first step.
Definition: timestep.C:134
FloatArray eNorm
Definition: supg.h:120
Class representing field of primary variables, which are typically allocated on nodes.
bool isNotEmpty() const
Tests for empty matrix.
Definition: floatmatrix.h:162
virtual void assemble(SparseMtrx &answer, TimeStep *tStep, const MatrixAssembler &ma, const UnknownNumberingScheme &s, Domain *domain)
Assembles characteristic matrix of required type into given sparse matrix.
Definition: engngm.C:803
virtual void computeVectorOfPressures(ValueModeType mode, TimeStep *tStep, FloatArray &pressures)
Definition: fmelement.C:55
void updateSolutionVectors(FloatArray &solutionVector, FloatArray &accelerationVector, FloatArray &incrementalSolutionVector, TimeStep *tStep)
Definition: supg.C:1186
virtual void computeDiffusionTerm_MB(FloatArray &answer, TimeStep *tStep)=0
Computes diffusion terms for momentum balance equations(s).
MetaStep * giveCurrentMetaStep()
Returns current meta step.
Definition: engngm.C:1684
std::vector< std::unique_ptr< GeneralBoundaryCondition > > & giveBcs()
Definition: domain.h:322
#define _IFT_SUPG_atolv
Definition: supg.h:60
int giveNumber()
Returns receiver&#39;s number.
Definition: timestep.h:129
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
virtual void computeDiffusionTerm_MC(FloatArray &answer, TimeStep *tStep)=0
Computes diffusion terms for mass conservation equation.
Callback class for assembling specific types of vectors.
double uscale
Definition: supg.h:142
virtual void computeVectorOfVelocities(ValueModeType mode, TimeStep *tStep, FloatArray &velocities)
Definition: fmelement.C:48
double Re
Reynolds number.
Definition: supg.h:144
void updateDofUnknownsDictionary_corrector(TimeStep *tStep)
Definition: supg.C:1096
NumericalCmpn
Type representing numerical component.
Definition: numericalcmpn.h:46
#define _IFT_EngngModel_smtype
Definition: engngm.h:82
std::unique_ptr< PrimaryField > VelocityPressureField
Definition: supg.h:114
#define OOFEM_ERROR(...)
Definition: error.h:61
EngngModelContext * giveContext()
Context requesting service.
Definition: engngm.h:1078
Callback class for assembling specific types of matrices.
double computeSquaredNorm() const
Computes the square of the norm.
Definition: floatarray.C:846
int ndomains
Number of receiver domains.
Definition: engngm.h:205
DofIDItem
Type representing particular dof type.
Definition: dofiditem.h:86
void updateSolutionVectors_predictor(FloatArray &solutionVector, FloatArray &accelerationVector, TimeStep *tStep)
Definition: supg.C:1140
int updateSharedDofManagers(FloatArray &answer, const UnknownNumberingScheme &s, int ExchangeTag)
Exchanges necessary remote DofManagers data.
Definition: engngm.C:1957
SparseMtrxType
Enumerative type used to identify the sparse matrix type.
The reference incremental load vector is defined as totalLoadVector assembled at given time...
Implementation for assembling external forces vectors in standard monolithic FE-problems.
double atolv
Convergence tolerance.
Definition: supg.h:126
void times(double f)
Multiplies receiver by factor f.
Definition: floatmatrix.C:1594
double rtolv
Definition: supg.h:126
bool stopmaxiter
Flag if set to true (default), then when max number of iteration reached, computation stops otherwise...
Definition: supg.h:133
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
virtual void updateYourself(TimeStep *tStep)
Updates internal state after finishing time step.
Definition: engngm.C:612
DofIDItem giveDofID() const
Returns DofID value of receiver, which determines type of of unknown connected to receiver (e...
Definition: dof.h:276
#define _IFT_SUPG_uscale
Definition: supg.h:56
#define _IFT_SUPG_lscale
Definition: supg.h:55
#define _IFT_SUPG_alpha
Definition: supg.h:53
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
std::vector< std::unique_ptr< InitialCondition > > & giveIcs()
Definition: domain.h:329
MatResponseMode rmode
Definition: supg.h:91
EngngModel * master
Master e-model; if defined receiver is in maintained (slave) mode.
Definition: engngm.h:262
FloatArray accelerationVector
Definition: supg.h:116
#define _IFT_SUPG_rtolv
Definition: supg.h:59
Callback class for assembling SUPG tangent matrices.
Definition: supg.h:88
virtual void matrixFromElement(FloatMatrix &mat, Element &element, TimeStep *tStep) const
Definition: supg.C:152
TimeStep * givePreviousStep()
Returns pointer to previous solution step.
Definition: timestep.C:114
void updateInternalState(TimeStep *tStep)
Definition: supg.C:791
virtual void computeDiffusionDerivativeTerm_MB(FloatMatrix &answer, MatResponseMode mode, TimeStep *tStep)=0
Computes the derivative of diffusion terms for momentum balance equations(s) with respect to nodal ve...
Callback class for assembling SUPG internal forces.
Definition: supg.h:74
int deltaTF
Definition: supg.h:124
SUPG(int i, EngngModel *_master=NULL)
Definition: supg.C:210
#define _IFT_SUPG_miflag
Definition: supg.h:58
virtual void computeAccelerationTerm_MB(FloatMatrix &answer, TimeStep *tStep)=0
Computes acceleration terms (generalized mass matrix with stabilization terms) for momentum balance e...
virtual void giveLocalVelocityDofMap(IntArray &map)
Definition: supgelement.h:208
Function * giveFunction(int n)
Service for accessing particular domain load time function.
Definition: domain.C:268
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode)
Stores the state of model to output stream.
Definition: engngm.C:1527
Abstract base class representing Lagrangian-Eulerian (moving) material interfaces.
Definition: leplic.h:153
virtual void computeAccelerationTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes acceleration terms for mass conservation equation.
#define _IFT_SUPG_dscale
Definition: supg.h:57
Initializes the variable VERBOSE, in order to get a few intermediate messages on screen: beginning an...
void registerField(FieldPtr eField, FieldType key)
Registers the given field (the receiver is not assumed to own given field).
Definition: fieldmanager.C:42
Class representing vector of real numbers.
Definition: floatarray.h:82
SparseLinearSystemNM * createSparseLinSolver(LinSystSolverType st, Domain *d, EngngModel *m)
Creates new instance of SparseLinearSystemNM corresponding to given type.
Definition: classfactory.C:120
virtual int checkConsistency()
Allows programmer to test some receiver&#39;s internal data, before computation begins.
Definition: engngm.h:995
This abstract class represent a general base element class for fluid dynamic problems.
Definition: fmelement.h:54
std::unique_ptr< MaterialInterface > materialInterface
Definition: supg.h:147
virtual void computeBCLhsPressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - pressure.
Definition: supgelement.C:316
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
#define _IFT_EngngModel_lstype
Definition: engngm.h:81
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
void evaluateElementStabilizationCoeffs(TimeStep *tStep)
Definition: supg.C:1020
virtual void computeDiffusionDerivativeTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes diffusion derivative terms for mass conservation equation.
virtual void computeBCLhsPressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - pressure.
Definition: supgelement.C:288
Abstract base class representing Level Set representation of material interfaces. ...
Definition: levelsetpcs.h:114
void assemble(const FloatArray &fe, const IntArray &loc)
Assembles the array fe (typically, the load vector of a finite element) into the receiver, using loc as location array.
Definition: floatarray.C:551
virtual void vectorFromElement(FloatArray &vec, Element &element, TimeStep *tStep, ValueModeType mode) const
Definition: supg.C:68
int giveNumberOfTimeStepWhenIcApply()
Returns the time step number, when initial conditions should apply.
Definition: engngm.h:754
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual void reinitialize()
Reinitializes the receiver.
Definition: nummet.h:112
Class representing the general Input Record.
Definition: inputrecord.h:101
void updateDofUnknownsDictionary_predictor(TimeStep *tStep)
Definition: supg.C:1060
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
virtual void computeAdvectionTerm_MB(FloatArray &answer, TimeStep *tStep)=0
Computes nonlinear advection terms for momentum balance equations(s).
Class implementing single timer, providing wall clock and user time capabilities. ...
Definition: timer.h:46
#define _IFT_SUPG_cmflag
Definition: supg.h:52
#define _IFT_SUPG_scaleflag
Definition: supg.h:54
ClassFactory & classFactory
Definition: classfactory.C:59
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
#define _IFT_SUPG_stopmaxiter
Definition: supg.h:62
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode)
Restores the state of model from output stream.
Definition: supg.C:835
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
virtual TimeStep * givePreviousStep(bool force=false)
Returns previous time step.
Definition: engngm.h:693
virtual int giveUnknownDictHashIndx(ValueModeType mode, TimeStep *tStep)
This method is responsible for computing unique dictionary id (ie hash value) from given valueModeTyp...
Definition: supg.C:1128
virtual TimeStep * giveSolutionStepWhenIcApply(bool force=false)
Returns the solution step when Initial Conditions (IC) apply.
Definition: engngm.h:720
virtual double giveVariableScale(VarScaleType varId)
Returns the scale factor for given variable type.
Definition: supg.C:996
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
std::vector< std::unique_ptr< Element > > & giveElements()
Definition: domain.h:279
int initFlag
Definition: supg.h:137
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
virtual TimeStep * giveCurrentStep(bool force=false)
Returns current time step.
Definition: engngm.h:683
void applyIC(TimeStep *tStep)
Definition: supg.C:935
virtual void times(const FloatArray &x, FloatArray &answer) const
Evaluates .
Definition: sparsemtrx.h:137
virtual void updateDomainLinks()
Updates domain links after the domains of receiver have changed.
Definition: engngm.C:1521
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode)
Restores the state of model from output stream.
Definition: engngm.C:1592
the oofem namespace is to define a context or scope in which all oofem names are defined.
void assemble(const FloatMatrix &src, const IntArray &loc)
Assembles the contribution using localization array into receiver.
Definition: floatmatrix.C:251
Domain * giveDomain(int n)
Service for accessing particular problem domain.
Definition: engngm.C:1720
virtual double giveReynoldsNumber()
Definition: supg.C:387
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
Abstract class Dof represents Degree Of Freedom in finite element mesh.
Definition: dof.h:93
DofManager * giveDofManager(int n)
Service for accessing particular domain dof manager.
Definition: domain.C:314
double alpha
Integration constant.
Definition: supg.h:135
#define _IFT_SUPG_deltat
Definition: supg.h:50
void startTimer()
Definition: timer.C:69
virtual double evaluateAtTime(double t)
Returns the value of the function at given time.
Definition: function.C:76
double getUtime()
Returns total user time elapsed in seconds.
Definition: timer.C:105
#define OOFEM_WARNING(...)
Definition: error.h:62
void assembleVector(FloatArray &answer, TimeStep *tStep, const VectorAssembler &va, ValueModeType mode, const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms=NULL)
Assembles characteristic vector of required type from dofManagers, element, and active boundary condi...
Definition: engngm.C:986
virtual void computeAdvectionTerm_MC(FloatArray &answer, TimeStep *tStep)=0
Computes advection terms for mass conservation equation.
Class representing solution step.
Definition: timestep.h:80
virtual void computePressureTerm_MC(FloatMatrix &answer, TimeStep *tStep)=0
Computes pressure terms for mass conservation equation.
#define _IFT_SUPG_maxiter
Definition: supg.h:61
void stopTimer()
Definition: timer.C:77
double deltaT
Definition: supg.h:123
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
virtual void printSingleOutputAt(FILE *file, TimeStep *tStep, char ch, ValueModeType mode, double scale=1.0)
Prints Dof output (it prints value of unknown related to dof at given timeStep).
Definition: dof.C:76
SUPGInternalForceAssembler(double l, double d, double u)
Definition: supg.C:64
virtual void computeBCLhsTerm_MB(FloatMatrix &answer, TimeStep *tStep)
Computes Lhs terms due to boundary conditions - velocity.
Definition: supgelement.C:243
FloatArray incrementalSolutionVector
Definition: supg.h:117
virtual NumericalMethod * giveNumericalMethod(MetaStep *mStep)
Returns reference to receiver&#39;s numerical method.
Definition: supg.C:225
const char * __ValueModeTypeToString(ValueModeType _value)
Definition: cltypes.C:322
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
virtual void computeAdvectionDerivativeTerm_MB(FloatMatrix &answer, TimeStep *tStep)=0
Computes the derivative of advection terms for momentum balance equations(s) with respect to nodal ve...

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:31 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011