OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
huertaerrorestimator.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 "../sm/ErrorEstimators/huertaerrorestimator.h"
36 #include "../sm/EngineeringModels/adaptnlinearstatic.h"
37 #include "../sm/EngineeringModels/adaptlinearstatic.h"
38 #include "domain.h"
39 #include "node.h"
40 #include "element.h"
41 #include "load.h"
42 #include "floatarray.h"
43 #include "floatmatrix.h"
44 #include "mathfem.h"
45 #include "function.h"
46 #include "timestep.h"
47 #include "metastep.h"
48 #include "integrationrule.h"
49 #include "connectivitytable.h"
50 #include "crosssection.h"
51 #include "dof.h"
52 #include "util.h"
53 #include "eleminterpunknownmapper.h"
54 #include "verbose.h"
55 #include "datastream.h"
56 #include "contextioerr.h"
57 #include "timer.h"
58 #include "calmls.h"
59 #include "nrsolver.h"
60 #include "errorestimatortype.h"
61 #include "classfactory.h"
62 #include "dynamicdatareader.h"
63 #include "dynamicinputrecord.h"
64 #include "heavisidetimefunction.h"
65 #include "outputmanager.h"
66 #include "boundarycondition.h"
67 #include "feinterpol.h"
68 #include "gausspoint.h"
69 #include "unknownnumberingscheme.h"
70 
71 #include <vector>
72 #include <string>
73 
74 
75 namespace oofem {
76 REGISTER_ErrorEstimator(HuertaErrorEstimator, EET_HEE);
77 
78 //#define STIFFNESS_TYPE TangentStiffnessMatrix
79 #define STIFFNESS_TYPE ElasticStiffnessMatrix
80 //#define STIFFNESS_TYPE SecantStiffnessMatrix
81 
82 #define ERROR_EXCESS 0.1 // 10 %
83 
84 
85 // debug defines
86 
87 //#define USE_FILE
88 
89 #define PRINT_ERROR
90 
91 #ifdef VERBOSE
92  #define INFO
93 //#define TIME_INFO
94  #define EXACT_ERROR
95  #define PRINT_ERROR
96 #endif
97 
98 #ifdef USE_FILE
99  #define USE_INPUT_FILE
100 //#define USE_OUTPUT_FILE
101 //#define USE_CONTEXT_FILE
102 #endif
103 
104 #ifdef PRINT_ERROR
105 //#define PRINT_FINE_ERROR
106  #define PRINT_COARSE_ERROR
107 #endif
108 
109 static bool masterRun = true;
110 
111 static bool exactFlag = false;
112 
113 #ifdef EXACT_ERROR
114 static bool wholeFlag = false, huertaFlag = false;
115 
116 double exactENorm, coarseUNorm, fineUNorm, mixedNorm;
117 
118  #ifdef PRINT_ERROR
119 static int finePos;
120 static FloatArray exactFineError;
121 static FloatArray exactCoarseError;
122  #endif
123 #endif
124 
125 //static FloatArray uNormArray;
126 
127 static DynamicDataReader refinedReader("huerta");
128 
129 static int impCSect, perCSect;
130 static FloatArray impPos;
131 
132 static int globalNelems;
133 
134 
135 int
136 HuertaErrorEstimator :: estimateError(EE_ErrorMode err_mode, TimeStep *tStep)
137 {
138  Domain *d = this->domain;
139  EngngModel *model = d->giveEngngModel();
140  int ielem, nelems = d->giveNumberOfElements();
141  int inode, nnodes = d->giveNumberOfDofManagers();
142  double pe;
143  IntArray localNodeIdArray, globalNodeIdArray; // these arrays are declared here to
144  // prevent their repeated creation for
145  // each element and patch problem
146 
147  if ( this->stateCounter == tStep->giveSolutionStateCounter() ) {
148  return 1;
149  }
150 
151  this->globalENorm = 0.0;
152  this->globalUNorm = 0.0;
153  this->globalWENorm = 0.0;
154 
155  // initiate to prevent oofeg failure (which loads eNorms from context)
156  this->eNorms.resize(nelems);
157  this->eNorms.zero();
158 
159  if ( initialSkipSteps != 0 ) {
160  initialSkipSteps--;
161 
162  skippedNelems = nelems;
163  this->stateCounter = tStep->giveSolutionStateCounter();
164  OOFEM_LOG_RELEVANT( "Relative error estimate [step number %5d]: skipped\n",
165  d->giveEngngModel()->giveCurrentStep()->giveNumber() );
166  return 1;
167  }
168 
169  if ( stepsToSkip != 0 ) {
170  stepsToSkip--;
171  skippedSteps++;
172 
173  skippedNelems = nelems;
174  this->stateCounter = tStep->giveSolutionStateCounter();
175  OOFEM_LOG_RELEVANT( "\nRelative error estimate [step number %5d]: skipped\n",
176  d->giveEngngModel()->giveCurrentStep()->giveNumber() );
177  return 1;
178  }
179 
180  OOFEM_LOG_INFO( "[%d] Estimating error [step number %5d]\n",
181  d->giveEngngModel()->giveRank(),
182  d->giveEngngModel()->giveCurrentStep()->giveNumber() );
183 
184  if ( dynamic_cast< AdaptiveLinearStatic * >( d->giveEngngModel() ) ) {
185  this->mode = HEE_linear;
186  } else if ( dynamic_cast< AdaptiveNonLinearStatic * >( d->giveEngngModel() ) ) {
187  this->mode = HEE_nlinear;
188  } else {
189  OOFEM_ERROR("Unsupported analysis type");
190  this->mode = HEE_linear;
191  }
192 
193  // check if each node has default number of dofs
194  // check if first load time function is constant
195  // check if only constant edge or surface load is used
196 
197  this->buildRefinedMesh();
198 
199  localNodeIdArray.resize(this->refinedMesh.nodes);
200  localNodeIdArray.zero();
201 
202  this->skippedNelems = 0;
203 
204 #ifdef EXACT_ERROR
205  if ( exactFlag == true ) {
206  #ifdef PRINT_ERROR
207  finePos = 0;
208  exactFineError.resize(this->refinedMesh.elems);
209  exactCoarseError.resize(nelems);
210  #endif
211 
212  wholeFlag = true;
213  this->solveRefinedWholeProblem(localNodeIdArray, globalNodeIdArray, tStep);
214  wholeFlag = false;
215 
216  if ( huertaFlag == false ) {
217  OOFEM_LOG_INFO("\n\n");
218  OOFEM_LOG_INFO("Exact eNorm2: %15.8e\n", exactENorm);
219  OOFEM_LOG_INFO("Exact coarse uNorm2: %15.8e\n", coarseUNorm);
220  // fprintf(stdout, "Exact mixed euNorm: %15.8e\n", mixedNorm);
221  OOFEM_LOG_INFO("Exact fine uNorm2: %15.8e\n", fineUNorm);
222 
223  pe = sqrt( exactENorm / ( exactENorm + coarseUNorm ) );
224  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
225  OOFEM_LOG_INFO("Exact relative error (coarse): %6.3f%% (L2 norm)\n", pe * 100.0);
226  }
227 
228  if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
229  OOFEM_LOG_INFO("Exact relative error (coarse): %6.3f%% (energy norm)\n", pe * 100.0);
230  }
231 
232  pe = sqrt(exactENorm / fineUNorm);
233  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
234  OOFEM_LOG_INFO("Exact relative error (fine): %6.3f%% (L2 norm)\n", pe * 100.0);
235  }
236 
237  if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
238  OOFEM_LOG_INFO("Exact relative error (fine): %6.3f%% (energy norm)\n", pe * 100.0);
239  }
240 
241  exit(1);
242  }
243  }
244 
245 #endif
246 
247  primaryUnknownError.resize( this->refinedMesh.nodes * d->giveDofManager(1)->giveNumberOfDofs() );
248  primaryUnknownError.zero();
249 
250  // uNormArray.resize(nelems);
251 
252  // first loop over patches and then over elements;
253  // this is advantageous, because when orthogonalizing element and patch error I have
254  // to assemble again the stifness matrix that can be then used (as whole) for evaluation
255  // of error in the energy norm;
256  // when using opposite way, I have to evaluate error on fine elements and than sum
257  // contributions to coarse elements;
258  // in this case it would be better to evaluate the error on the patch as a whole
259  // and assoociating it with the corresponding node;
260  // but this would complicate remeshing criterion
261 
262  // however there is problem that in nonlinear analysis stiffness matrix of the same fine
263  // element is different for element and for patch problem because each may be at slightly
264  // different point of the solution
265 
266  // freopen("/dev/null", "w", stdout);
267 
268  for ( inode = 1; inode <= nnodes; inode++ ) {
269  this->solveRefinedPatchProblem(inode, localNodeIdArray, globalNodeIdArray, tStep);
270  }
271 
272  for ( ielem = 1; ielem <= nelems; ielem++ ) {
273  this->solveRefinedElementProblem(ielem, localNodeIdArray, globalNodeIdArray, tStep);
274  }
275 
276 #ifdef __PARALLEL_MODE
277  #ifdef __USE_MPI
278  double buffer_out [ 4 ], buffer_in [ 4 ];
279 
280  buffer_out [ 0 ] = globalENorm;
281  buffer_out [ 1 ] = globalUNorm;
282  buffer_out [ 2 ] = ( double ) nelems;
283  buffer_out [ 3 ] = ( double ) skippedNelems;
284 
285  MPI_Allreduce(buffer_out, buffer_in, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
286 
287  globalENorm = buffer_in [ 0 ];
288  globalUNorm = buffer_in [ 1 ];
289  globalNelems = ( int ) ( buffer_in [ 2 ] + 0.1 );
290  skippedNelems = ( int ) ( buffer_in [ 3 ] + 0.1 );
291  #endif
292 #else
293  globalNelems = nelems;
294 #endif
295 
296 #ifdef PRINT_COARSE_ERROR
297  OOFEM_LOG_DEBUG("\n");
298  if ( exactFlag == true ) {
299  OOFEM_LOG_DEBUG(" elem a_Error2 x_Error2 a/x rate2 \n");
300  OOFEM_LOG_DEBUG("---------------------------------------------------------\n");
301  } else {
302  OOFEM_LOG_DEBUG(" elem a_Error2 \n");
303  OOFEM_LOG_DEBUG("-----------------------\n");
304  }
305 
306  for ( ielem = 1; ielem <= nelems; ielem++ ) {
307  if ( exactFlag == false ) {
308  OOFEM_LOG_DEBUG("[%d] %5d: %15.8e %s\n", d->giveEngngModel()->giveRank(), ielem,
309  this->eNorms.at(ielem) * this->eNorms.at(ielem),
310  ( this->skipRegion( d->giveElement(ielem)->giveRegionNumber() ) != 0 ) ? "(skipped)" :
311  ( d->giveElement(ielem)->giveParallelMode() == Element_remote ) ? "(remote)" : "");
312  }
313 
314  #ifdef EXACT_ERROR
315  else {
316  if ( fabs( exactCoarseError.at(ielem) ) > 1.0e-30 && this->eNorms.at(ielem) != 0.0 ) {
317  OOFEM_LOG_DEBUG("%5d: %15.8e %15.8e %15.8e %s\n", ielem,
318  this->eNorms.at(ielem) * this->eNorms.at(ielem), exactCoarseError.at(ielem),
319  this->eNorms.at(ielem) / sqrt( exactCoarseError.at(ielem) ),
320  ( this->skipRegion( d->giveElement(ielem)->giveRegionNumber() ) != 0 ) ? "(skipped)" : "");
321  } else {
322  OOFEM_LOG_DEBUG("%5d: %15.8e %15.8e N/A %s\n", ielem,
323  this->eNorms.at(ielem) * this->eNorms.at(ielem), exactCoarseError.at(ielem),
324  ( this->skipRegion( d->giveElement(ielem)->giveRegionNumber() ) != 0 ) ? "(skipped)" : "");
325  }
326  }
327  #endif
328  }
329 
330 #endif
331 
332  double pwe = 0.0;
333  if ( wError == true ) {
334  // calculate weighted error
335  double elemErrLimit;
336 
337  elemErrLimit = sqrt( ( this->globalENorm + this->globalUNorm ) / ( globalNelems - skippedNelems ) ) * requiredError;
338  if ( elemErrLimit != 0.0 ) {
339  double eerror, iratio;
340 
341  for ( ielem = 1; ielem <= nelems; ielem++ ) {
342  eerror = this->eNorms.at(ielem);
343  iratio = eerror / elemErrLimit;
344  globalWENorm += eerror * eerror * iratio;
345  }
346 
347 #ifdef __PARALLEL_MODE
348  #ifdef __USE_MPI
349  double myGlobalWENorm = globalWENorm;
350  MPI_Allreduce(& myGlobalWENorm, & globalWENorm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
351  #endif
352 #endif
353 
354  pwe = sqrt( globalWENorm / ( this->globalENorm + this->globalUNorm ) );
355  }
356  }
357 
358  OOFEM_LOG_INFO("\n");
359  OOFEM_LOG_INFO("Global eNorm2: %15.8e\n", this->globalENorm);
360  if ( wError == true ) {
361  OOFEM_LOG_INFO("Global wNorm2: %15.8e\n", globalWENorm);
362  }
363 
364  OOFEM_LOG_INFO("Global uNorm2: %15.8e\n", this->globalUNorm);
365 
366  // report the error estimate
367  pe = sqrt( this->globalENorm / ( this->globalENorm + this->globalUNorm ) );
368  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
369  OOFEM_LOG_INFO("Relative error estimate [step number %5d]: %6.3f%% (L2 norm)\n",
370  d->giveEngngModel()->giveCurrentStep()->giveNumber(), pe * 100.0);
371  }
372 
373  if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
374  OOFEM_LOG_INFO("Relative error estimate [step number %5d]: %6.3f%% (energy norm)\n",
375  d->giveEngngModel()->giveCurrentStep()->giveNumber(), pe * 100.0);
376  }
377 
378  if ( wError == true ) {
379  OOFEM_LOG_INFO("Relative error estimate [step number %5d]: %6.3f%% (weighted)\n",
380  d->giveEngngModel()->giveCurrentStep()->giveNumber(), pwe * 100.0);
381  }
382  this->globalErrorEstimate = pe;
383 
384  //fflush(stdout);
385 
386  if ( maxSkipSteps != 0 && this->mode == HEE_nlinear ) {
387  if ( lastError < 0.0 ) {
388  lastError = pe;
389  stepsToSkip = 0;
390  } else {
391  if ( requiredError > pe ) {
392  if ( pe <= lastError ) {
393  stepsToSkip = maxSkipSteps;
394  } else {
395  // estimate number of steps to skip using linear extrapolation
396  stepsToSkip = ( int ) ( ( requiredError - pe ) / ( pe - lastError ) * ( skippedSteps + 1 ) );
397  // make the number of steps to skip safe with respect to how many steps have been skipped last time
398  stepsToSkip = stepsToSkip / ( skippedSteps + 1 );
399  }
400 
401  // make the number of steps to skip safe with respect to the absolute value of the error
402  // (decrease by 0.5 for pe equal requiredError)
403  stepsToSkip = ( int ) ( stepsToSkip * ( ( 0.5 - 1.0 ) / requiredError * pe + 1 ) + 0.5 );
404  if ( skippedSteps == 0 ) {
405  stepsToSkip--;
406  }
407 
408  if ( stepsToSkip > maxSkipSteps ) {
409  stepsToSkip = maxSkipSteps;
410  }
411 
412  if ( stepsToSkip < 0 ) {
413  stepsToSkip = 0;
414  }
415  } else {
416  // allow error excess and prevent recursion
417  if ( requiredError * ( 1.00 + ERROR_EXCESS ) < pe && skippedSteps != 0 ) {
418  /* HUHU CHEATING do not return just by one step */
419  if ( maxSkipSteps == 1 ) {
420  return 1;
421  }
422 
423  int tStepNumber, curNumber = tStep->giveNumber();
424 
425  tStepNumber = curNumber - skippedSteps;
426 
427  OOFEM_LOG_INFO("Returning to step %5d\n", tStepNumber);
428 
429  skippedSteps = 0; // prevent recursion when returning
430  maxSkipSteps /= 2; // prevent oscillation
431 
432  // I trying to avoid repeated solution of (all, not only the one to which I am returning) previous steps
433  // IMPORTANT: I must restore context only from steps corresponding to this session !!!
434  // this means I cannot go in previous adaptive run, because there was a different domain
435  // it would be much cleaner to call restore from engng model
436  while ( tStepNumber < curNumber ) {
437  try {
438  FileDataStream stream(model->giveContextFileName(tStepNumber, 0), false);
439  model->restoreContext(stream, CM_State );
440  } catch(ContextIOERR & c) {
441  c.print();
442  exit(1);
443  }
444 
445  stepsToSkip = 0;
446  model->giveCurrentStep()->incrementStateCounter();
447  this->estimateError( temporaryEM, model->giveCurrentStep() );
448 
449  if ( lastError > requiredError ) {
450  return 1;
451  }
452 
453  tStepNumber += ( skippedSteps = stepsToSkip ) + 1;
454  }
455 
456  return 1;
457  }
458  }
459 
460  lastError = pe;
461  }
462 
463  skippedSteps = 0;
464  }
465 
466 #ifdef EXACT_ERROR
467  if ( exactFlag == true ) {
468  #ifdef __PARALLEL_MODE
469  #ifdef __USE_MPI
470  buffer_out [ 0 ] = exactENorm;
471  buffer_out [ 1 ] = coarseUNorm;
472  buffer_out [ 2 ] = mixedNorm;
473  buffer_out [ 3 ] = fineUNorm;
474 
475  MPI_Allreduce(buffer_out, buffer_in, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
476 
477  exactENorm = buffer_in [ 0 ];
478  coarseUNorm = buffer_in [ 1 ];
479  mixedNorm = buffer_in [ 2 ];
480  fineUNorm = buffer_in [ 3 ];
481  #endif
482  #endif
483 
484  OOFEM_LOG_INFO("\n");
485  OOFEM_LOG_INFO("Exact eNorm2: %15.8e\n", exactENorm);
486  OOFEM_LOG_INFO("Exact coarse uNorm2: %15.8e\n", coarseUNorm);
487  //fprintf(stdout, "Exact mixed euNorm: %15.8e\n", mixedNorm);
488  OOFEM_LOG_INFO("Exact fine uNorm2: %15.8e\n", fineUNorm);
489 
490  pe = sqrt( exactENorm / ( exactENorm + coarseUNorm ) );
491  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
492  OOFEM_LOG_INFO("Exact relative error (coarse): %6.3f%% (L2 norm)\n", pe * 100.0);
493  }
494 
495  if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
496  OOFEM_LOG_INFO("Exact relative error (coarse): %6.3f%% (energy norm)\n", pe * 100.0);
497  }
498 
499  pe = sqrt(exactENorm / fineUNorm);
500  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
501  OOFEM_LOG_INFO("Exact relative error (fine): %6.3f%% (L2 norm)\n", pe * 100.0);
502  }
503 
504  if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
505  OOFEM_LOG_INFO("Exact relative error (fine): %6.3f%% (energy norm)\n", pe * 100.0);
506  }
507  }
508 
509 #endif
510 
511  this->globalENorm = sqrt(this->globalENorm);
512  this->globalUNorm = sqrt(this->globalUNorm);
513  if ( wError == true ) {
514  this->globalWENorm = sqrt(this->globalWENorm);
515  }
516 
517  this->stateCounter = tStep->giveSolutionStateCounter();
518 
519  return 1;
520 }
521 
522 
523 
524 double
525 HuertaErrorEstimator :: giveElementError(EE_ErrorType type, Element *elem, TimeStep *tStep)
526 {
527  this->estimateError(equilibratedEM, tStep);
528  if ( type == primaryUnknownET ) {
529  return this->eNorms.at( elem->giveNumber() );
530  }
531 
532  return 0.0;
533 }
534 
535 
536 
537 double
538 HuertaErrorEstimator :: giveValue(EE_ValueType type, TimeStep *tStep)
539 {
540  this->estimateError(equilibratedEM, tStep);
541  if ( type == globalErrorEEV ) {
542  return this->globalENorm;
543  } else if ( type == globalNormEEV ) {
544  return this->globalUNorm;
545  } else if ( type == globalWeightedErrorEEV ) {
546  return this->globalWENorm;
547  } else if ( type == relativeErrorEstimateEEV ) {
548  return this->globalErrorEstimate;
549  }
550  return 0.0;
551 }
552 
553 
554 RemeshingCriteria *
555 HuertaErrorEstimator :: giveRemeshingCrit()
556 {
557  if ( !this->rc ) {
558  this->rc.reset( new HuertaRemeshingCriteria(1, this) );
559  }
560 
561  return this->rc.get();
562 }
563 
564 
565 IRResultType
566 HuertaErrorEstimator :: initializeFrom(InputRecord *ir)
567 {
568  IRResultType result; // Required by IR_GIVE_FIELD macro
569  int n, level, wErrorFlag = 0;
570 
571  ErrorEstimator :: initializeFrom(ir);
572  n = 1;
573  IR_GIVE_OPTIONAL_FIELD(ir, n, _IFT_HuertaErrorEstimator_normtype);
574  if ( n == 0 ) {
575  this->normType = L2Norm;
576  } else {
577  this->normType = EnergyNorm; // default
578  }
579 
580  level = this->refineLevel;
581  IR_GIVE_OPTIONAL_FIELD(ir, level, _IFT_HuertaErrorEstimator_refinelevel);
582  if ( level >= 0 ) {
583  this->refineLevel = level;
584  }
585 
586  IR_GIVE_FIELD(ir, this->requiredError, _IFT_HuertaErrorEstimator_requirederror);
587 
588  IR_GIVE_OPTIONAL_FIELD(ir, maxSkipSteps, _IFT_HuertaErrorEstimator_skipsteps);
589  if ( maxSkipSteps < 0 ) {
590  maxSkipSteps = 0;
591  }
592 
593  if ( maxSkipSteps > 5 ) {
594  maxSkipSteps = 5;
595  }
596 
597  IR_GIVE_OPTIONAL_FIELD(ir, initialSkipSteps, _IFT_HuertaErrorEstimator_initialskipsteps);
598  if ( initialSkipSteps < 0 ) {
599  initialSkipSteps = 0;
600  }
601 
602  IR_GIVE_OPTIONAL_FIELD(ir, wErrorFlag, _IFT_HuertaErrorEstimator_werror);
603  if ( wErrorFlag != 0 ) {
604  wError = true;
605  }
606 
607  if ( masterRun == true ) { // prevent overwriting of static variables
608  masterRun = false;
609 
610  perCSect = 0;
611  IR_GIVE_OPTIONAL_FIELD(ir, perCSect, _IFT_HuertaErrorEstimator_perfectCSect);
612 
613  impCSect = 0;
614  IR_GIVE_OPTIONAL_FIELD(ir, impCSect, _IFT_HuertaErrorEstimator_impCSect);
615  IR_GIVE_OPTIONAL_FIELD(ir, impPos, _IFT_HuertaErrorEstimator_impPos);
616 
617  if ( impCSect != 0 && perCSect == 0 ) {
618  OOFEM_ERROR("Missing perfect material specification (through cross-section)");
619  }
620 
621 #ifdef EXACT_ERROR
622  n = 0;
623  IR_GIVE_OPTIONAL_FIELD(ir, n, _IFT_HuertaErrorEstimator_exact);
624  if ( n > 0 ) {
625  exactFlag = true;
626  if ( n != 1 ) {
627  huertaFlag = true; // run also error estimate
628  }
629  } else {
630  exactFlag = false;
631  }
632 
633 #endif
634  }
635 
636  return this->giveRemeshingCrit()->initializeFrom(ir);
637 }
638 
639 
640 contextIOResultType
641 HuertaErrorEstimator :: saveContext(DataStream &stream, ContextMode mode, void *obj)
642 {
643  contextIOResultType iores;
644  TimeStep *tStep = this->domain->giveEngngModel()->giveCurrentStep();
645 
646  if ( this->stateCounter != tStep->giveSolutionStateCounter() ) {
647  this->estimateError(equilibratedEM, tStep);
648  }
649 
650  // save parent class status
651  if ( ( iores = ErrorEstimator :: saveContext(stream, mode, obj) ) != CIO_OK ) {
652  THROW_CIOERR(iores);
653  }
654 
655  if ( ( iores = this->eNorms.storeYourself(stream) ) != CIO_OK ) {
656  THROW_CIOERR(iores);
657  }
658 
659  // write a raw data
660  if ( !stream.write(stateCounter) ) {
661  THROW_CIOERR(CIO_IOERR);
662  }
663 
664  return CIO_OK;
665 }
666 
667 
668 contextIOResultType
669 HuertaErrorEstimator :: restoreContext(DataStream &stream, ContextMode mode, void *obj)
670 {
671  contextIOResultType iores;
672 
673  // read parent class status
674  if ( ( iores = ErrorEstimator :: restoreContext(stream, mode, obj) ) != CIO_OK ) {
675  THROW_CIOERR(iores);
676  }
677 
678  if ( ( iores = eNorms.restoreYourself(stream) ) != CIO_OK ) {
679  THROW_CIOERR(iores);
680  }
681 
682  // read raw data
683  if ( !stream.read(stateCounter) ) {
684  THROW_CIOERR(CIO_IOERR);
685  }
686 
687  return CIO_OK;
688 }
689 
690 
691 HuertaRemeshingCriteria :: HuertaRemeshingCriteria(int n, ErrorEstimator *e) : RemeshingCriteria(n, e)
692 {
693  this->mode = primaryUnknownBased;
694  this->stateCounter = 0;
695  this->refineCoeff = 1.0;
696  this->noRemesh = false;
697  this->wError = false;
698 }
699 
700 
701 double
702 HuertaRemeshingCriteria :: giveRequiredDofManDensity(int num, TimeStep *tStep, int relative)
703 {
704  double size;
705 
706  this->estimateMeshDensities(tStep);
707  size = this->nodalDensities.at(num);
708  size = max(minElemSize, size);
709  if ( relative ) {
710  return size / this->giveDofManDensity(num);
711  }
712 
713  return size;
714 }
715 
716 
717 RemeshingStrategy
718 HuertaRemeshingCriteria :: giveRemeshingStrategy(TimeStep *tStep)
719 {
720  this->estimateMeshDensities(tStep);
721  return this->remeshingStrategy;
722 }
723 
724 
725 #define MAX_COARSE_RATE 2.0
726 #define MAX_REFINE_RATE 5.0
727 
728 int
729 HuertaRemeshingCriteria :: estimateMeshDensities(TimeStep *tStep)
730 {
731  int nelem, nnode, elemPolyOrder, skipped;
732  double globValNorm = 0.0, globValErrorNorm = 0.0, globValWErrorNorm = 0.0;
733  double globValNorm2, globValErrorNorm2, globValWErrorNorm2;
734  double eerror, iratio, currDensity, elemSize, elemErrLimit, percentError;
735  Element *ielem;
736  EE_ErrorType errorType = unknownET;
737  bool refine = false;
738  int result;
739  double sval, maxVal;
740  FloatArray val;
741 
742  static int run = 0;
743 
744  if ( stateCounter == tStep->giveSolutionStateCounter() ) {
745  return 1;
746  }
747 
748  this->remeshingStrategy = NoRemeshing_RS;
749  if ( this->noRemesh == true ) {
750  // remember time stamp
751  stateCounter = tStep->giveSolutionStateCounter();
752  return 1;
753  }
754 
755  nelem = this->domain->giveNumberOfElements();
756  nnode = this->domain->giveNumberOfDofManagers();
757 
758  skipped = this->ee->giveNumberOfSkippedElements();
759 
760 #ifndef __PARALLEL_MODE
761  globalNelems = nelem;
762 #endif
763 
764  if ( skipped == globalNelems ) {
765  // remember time stamp
766  stateCounter = tStep->giveSolutionStateCounter();
767  return 1;
768  }
769 
770  // compute element error limit based on equally distribution error principle
771  if ( mode == primaryUnknownBased ) {
772  globValNorm = this->ee->giveValue(globalNormEEV, tStep);
773  globValErrorNorm = this->ee->giveValue(globalErrorEEV, tStep);
774  globValWErrorNorm = this->ee->giveValue(globalWeightedErrorEEV, tStep);
775  errorType = primaryUnknownET;
776  } else {
777  OOFEM_ERROR("unsupported mode");
778  }
779 
780  globValNorm2 = globValNorm * globValNorm;
781  globValErrorNorm2 = globValErrorNorm * globValErrorNorm;
782  globValWErrorNorm2 = globValWErrorNorm * globValWErrorNorm;
783 
784  elemErrLimit = sqrt( ( globValNorm2 + globValErrorNorm2 ) / ( globalNelems - skipped ) ) * this->requiredError;
785  if ( elemErrLimit == 0.0 ) {
786  return 0;
787  }
788 
789  run++;
790 
791  this->nodalDensities.resize(nnode);
792 
793  std :: vector< int > connectedElems(nnode);
794  for ( int i = 0; i < nnode; i++ ) {
795  connectedElems [ i ] = 0;
796  }
797 
798  percentError = sqrt( globValErrorNorm2 / ( globValErrorNorm2 + globValNorm2 ) );
799  // test if solution is allowable
800  if ( percentError > this->requiredError ) {
801  this->remeshingStrategy = RemeshingFromPreviousState_RS;
802  } else {
803  if ( run > 100 ) {
804  OOFEM_LOG_INFO("hehe\n");
805  for ( int i = 1; i <= nelem; i++ ) {
806  ielem = domain->giveElement(i);
807 
808  // skipping these elements will result in reusing their current size;
809  // in adaptive analysis this leads to monotonical decrease of their size !!!
810  // it is better to either try to enlarge these elements (as they have zero error)
811  // or to presribe zero nodal mesh density (but this can be handled only by t3d)
812  /*
813  * if(skipped != 0){
814  * if(this -> ee -> skipRegion(ielem -> giveRegionNumber()) != 0)continue;
815  * }
816  */
817 
818  currDensity = ielem->computeMeanSize();
819 
820  //#ifdef HUHU
821  // toto je treba udelat obecne
822  // zjistit, zda material na prvku je nelokalni, po te z prvku vytahnout danou velicinu
823  // a pokud presahne danou mez pouzit mensi z elemSize a dane size
824  // chtelo by tez zaridit spusteni remeshingu, kdyz tato situace nastane -> combined
825  // huerta + error indicator
826  maxVal = 0.0;
827  for ( GaussPoint *gp: *ielem->giveDefaultIntegrationRulePtr() ) {
828  result = ielem->giveIPValue(val, gp, IST_PrincipalDamageTensor, tStep);
829  if ( result ) {
830  sval = val.computeNorm();
831  maxVal = max(maxVal, sval);
832  }
833  }
834 
835  if ( maxVal > 0.75 ) {
836  if ( currDensity > 1.0 ) {
837  OOFEM_LOG_INFO("step %d elem %d dam %e size %e\n", tStep->giveNumber(), i, maxVal, currDensity);
838  this->remeshingStrategy = RemeshingFromPreviousState_RS;
839  }
840  }
841  }
842  }
843  }
844 
845  if ( this->remeshingStrategy == NoRemeshing_RS ) {
846  // check whether to remesh because of low level of error
847  if ( this->ee->giveErrorEstimatorType() == EET_HEE ) {
848  /* HUHU toto je divne !!! pak se neuplatni refine viz dale */
849  if ( static_cast< HuertaErrorEstimator * >(this->ee)->giveAnalysisMode() == HuertaErrorEstimator :: HEE_linear ) {
850  stateCounter = tStep->giveSolutionStateCounter();
851  return 1;
852  }
853 
854  if ( wError == false ) {
855  if ( percentError <= this->requiredError ) {
856  if ( percentError >= this->requiredError * 0.5 * this->refineCoeff || run <= 5 ) {
857  for ( int i = 1; i <= nnode; i++ ) {
858  nodalDensities.at(i) = this->giveDofManDensity(i);
859  }
860 
861  stateCounter = tStep->giveSolutionStateCounter();
862  OOFEM_LOG_INFO("huhu\n");
863  return 1;
864  }
865  }
866  } else {
867  double pwe = sqrt( globValWErrorNorm2 / ( globValErrorNorm2 + globValNorm2 ) );
868  if ( pwe <= this->requiredError * 1.1 ) {
869  if ( pwe >= this->requiredError * 0.5 * this->refineCoeff || run <= 5 ) {
870  for ( int i = 1; i <= nnode; i++ ) {
871  nodalDensities.at(i) = this->giveDofManDensity(i);
872  }
873 
874  stateCounter = tStep->giveSolutionStateCounter();
875  return 1;
876  }
877  }
878  }
879 
880  OOFEM_LOG_INFO("haha\n");
881  this->remeshingStrategy = RemeshingFromPreviousState_RS;
882  }
883  }
884 
885  /* HUHU ma zamezit aby se opakovane restartovalo ze stejneho kroku, na coz davky nejsou pripravene */
886  if ( run == 1 && tStep->giveNumber() != 1 ) {
887  this->remeshingStrategy = NoRemeshing_RS;
888  for ( int i = 1; i <= nnode; i++ ) {
889  nodalDensities.at(i) = this->giveDofManDensity(i);
890  }
891 
892  stateCounter = tStep->giveSolutionStateCounter();
893  return 1;
894  }
895 
896  for ( int i = 1; i <= nelem; i++ ) {
897  ielem = domain->giveElement(i);
898 
899  // skipping these elements will result in reusing their current size;
900  // in adaptive analysis this leads to monotonical decrease of their size !!!
901  // it is better to either try to enlarge these elements (as they have zero error)
902  // or to presribe zero nodal mesh density (but this can be handled only by t3d)
903  /*
904  * if(skipped != 0){
905  * if(this -> ee -> skipRegion(ielem -> giveRegionNumber()) != 0)continue;
906  * }
907  */
908 
909  eerror = this->ee->giveElementError(errorType, ielem, tStep);
910  iratio = eerror / elemErrLimit;
911  if ( iratio > 1.0 ) {
912  refine = true;
913 
914  // checking the step number avoids application of the refine coefficient for linear problems
915  // or linear stage in nonlinear problems
916  // if(tStep -> giveNumber() != 1)iratio /= this -> refineCoeff;
917  iratio /= this->refineCoeff;
918  // limit the refinement (note there is also minElemSize)
919  // if (iratio > MAX_REFINE_RATE)iratio = MAX_REFINE_RATE;
920  } else {
921  // limit the coarsening
922  if ( iratio < 1.0 / MAX_COARSE_RATE ) {
923  iratio = 1.0 / MAX_COARSE_RATE;
924  }
925  }
926 
927  currDensity = ielem->computeMeanSize();
928  elemPolyOrder = ielem->giveInterpolation()->giveInterpolationOrder();
929  elemSize = currDensity / pow(iratio, 1.0 / elemPolyOrder);
930  //#ifdef HUHU
931  // toto je treba udelat obecne
932  // zjistit, zda material na prvku je nelokalni, po te z prvku vytahnout danou velicinu
933  // a pokud presahne danou mez pouzit mensi z elemSize a dane size
934  // chtelo by tez zaridit spusteni remeshingu, kdyz tato situace nastane -> combined
935  // huerta + error indicator
936  maxVal = 0.0;
937  for ( GaussPoint *gp: *ielem->giveDefaultIntegrationRulePtr() ) {
938  result = ielem->giveIPValue(val, gp, IST_PrincipalDamageTensor, tStep);
939  if ( result ) {
940  sval = val.computeNorm();
941  maxVal = max(maxVal, sval);
942  }
943  }
944 
945  if ( maxVal > 0.75 ) {
946  elemSize = min(elemSize, 1.0);
947  }
948 
949  //#endif
950  for ( int j = 1; j <= ielem->giveNumberOfDofManagers(); j++ ) {
951  int jnode = ielem->giveDofManager(j)->giveNumber();
952  if ( connectedElems [ jnode - 1 ] != 0 ) {
953  //this -> nodalDensities.at(jnode) = min(this -> nodalDensities.at(jnode), elemSize);
954  this->nodalDensities.at(jnode) += elemSize;
955  } else {
956  this->nodalDensities.at(jnode) = elemSize;
957  }
958 
959  connectedElems [ jnode - 1 ]++;
960  }
961  }
962 
963  // init non-initialized nodes -> those in skip regions
964  for ( int i = 0; i < nnode; i++ ) {
965  if ( connectedElems [ i ] == 0 ) {
966  this->nodalDensities.at(i + 1) = this->giveDofManDensity(i + 1);
967  } else {
968  this->nodalDensities.at(i + 1) /= connectedElems [ i ];
969  }
970  }
971 
972  if ( refine == true ) {
973  if ( this->ee->giveErrorEstimatorType() == EET_HEE ) {
974  if ( static_cast< HuertaErrorEstimator * >(this->ee)->giveAnalysisMode() == HuertaErrorEstimator :: HEE_linear ) {
975  this->remeshingStrategy = RemeshingFromPreviousState_RS;
976  }
977  }
978  }
979 
980  // remember time stamp
981  stateCounter = tStep->giveSolutionStateCounter();
982  return 1;
983 }
984 
985 
986 IRResultType
987 HuertaRemeshingCriteria :: initializeFrom(InputRecord *ir)
988 {
989  IRResultType result; // Required by IR_GIVE_FIELD macro
990  double coeff;
991  int noRemeshFlag = 0, wErrorFlag = 0;
992 
993  IR_GIVE_FIELD(ir, this->requiredError, _IFT_HuertaRemeshingCriteria_requirederror);
994  IR_GIVE_FIELD(ir, this->minElemSize, _IFT_HuertaRemeshingCriteria_minelemsize);
995 
996  IR_GIVE_OPTIONAL_FIELD(ir, noRemeshFlag, _IFT_HuertaRemeshingCriteria_noremesh);
997  if ( noRemeshFlag != 0 ) {
998  this->noRemesh = true;
999  }
1000 
1001 
1002  IR_GIVE_OPTIONAL_FIELD(ir, wErrorFlag, _IFT_HuertaRemeshingCriteria_werror);
1003  if ( wErrorFlag != 0 ) {
1004  this->wError = true;
1005  }
1006 
1007  coeff = this->refineCoeff;
1008  IR_GIVE_OPTIONAL_FIELD(ir, coeff, _IFT_HuertaRemeshingCriteria_refinecoeff);
1009  if ( coeff > 0.0 && coeff <= 1.0 ) {
1010  this->refineCoeff = coeff;
1011  }
1012 
1013  return IRRT_OK;
1014 }
1015 
1016 
1017 double
1018 HuertaRemeshingCriteria :: giveDofManDensity(int num)
1019 {
1020  int isize;
1021  ConnectivityTable *ct = domain->giveConnectivityTable();
1022  const IntArray *con;
1023  double density;
1024 
1025  con = ct->giveDofManConnectivityArray(num);
1026  isize = con->giveSize();
1027 
1028 #if 0
1029  // Minimum density
1030  for ( i = 1; i <= isize; i++ ) {
1031  Element *ielem = domain->giveElement( con->at(i) );
1032  if ( i == 1 ) {
1033  density = ielem->computeMeanSize();
1034  } else {
1035  density = min( density, ielem->computeMeanSize() );
1036  }
1037  }
1038 #endif
1039 
1040  // Average density
1041 
1042  density = 0.0;
1043  for ( int i = 1; i <= isize; i++ ) {
1044  Element *ielem = domain->giveElement( con->at(i) );
1045 
1046  density += ielem->computeMeanSize();
1047  }
1048 
1049  density /= isize;
1050 
1051  return density;
1052 }
1053 
1054 
1055 void
1056 HuertaErrorEstimator :: buildRefinedMesh()
1057 {
1058  int nelems;
1059  Domain *d = this->domain;
1060 
1061  if ( this->refinedMesh.completed == 1 ) {
1062  return;
1063  }
1064 
1065  nelems = d->giveNumberOfElements();
1066  this->refinedElementList.reserve(nelems);
1067 
1068  for ( int ielem = 1; ielem <= nelems; ielem++ ) {
1069  this->refinedElementList.emplace_back(d, ielem, this->refineLevel);
1070  }
1071 
1072  if ( refinedMesh.refineMeshGlobally(d, this->refineLevel, this->refinedElementList) != 0 ) {
1073  OOFEM_ERROR("refineMeshGlobally failed");
1074  }
1075 
1076  this->refinedMesh.completed = 1;
1077 }
1078 
1079 
1080 
1081 void
1082 HuertaErrorEstimatorInterface :: setupRefinedElementProblem1D(Element *element, RefinedElement *refinedElement,
1083  int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
1084  HuertaErrorEstimatorInterface :: SetupMode mode, TimeStep *tStep, int nodes,
1085  FloatArray **corner, FloatArray &midNode,
1086  int &localNodeId, int &localElemId, int &localBcId,
1087  IntArray &controlNode, IntArray &controlDof,
1088  HuertaErrorEstimator :: AnalysisMode aMode, const char *edgetype)
1089 {
1090  std :: list< FloatArray >newNodes;
1091  IntArray *connectivity, boundary(1);
1092  int startNode, endNode, pos, nd, bc, dofs;
1093 
1094  if ( nodeId != 0 ) {
1095  startNode = endNode = nodeId;
1096  } else {
1097  startNode = 1;
1098  endNode = nodes;
1099  }
1100 
1101  dofs = element->giveDomain()->giveDofManager(1)->giveNumberOfDofs();
1102 
1103  if ( mode == HuertaErrorEstimatorInterface :: CountMode ) {
1104  for ( int inode = startNode; inode <= endNode; inode++ ) {
1105  localElemId += ( level + 1 );
1106 
1107  connectivity = refinedElement->giveFineNodeArray(inode);
1108  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1109 
1110  pos = 1;
1111  for ( int m = 0; m < level + 2; m++, pos++ ) {
1112  nd = connectivity->at(pos);
1113  if ( localNodeIdArray.at(nd) == 0 ) {
1114  localNodeIdArray.at(nd) = ++localNodeId;
1115 
1116  // count boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1117 
1118  bc = 0;
1119  if ( nodeId == 0 ) {
1120  if ( m == 0 && boundary.at(1) == 0 ) {
1121  bc = 1;
1122  }
1123  } else {
1124  if ( m == level + 1 ) {
1125  bc = 1;
1126  }
1127  }
1128 
1129 #ifdef EXACT_ERROR
1130  if ( wholeFlag == true ) {
1131  bc = 0;
1132  }
1133 
1134 #endif
1135 
1136  if ( bc == 1 ) {
1137  localBcId += dofs;
1138  }
1139  }
1140  }
1141  }
1142 
1143  return;
1144  }
1145 
1146  if ( mode == HuertaErrorEstimatorInterface :: NodeMode ) {
1147  double x, y, z, u, du = 1.0 / ( level + 1 );
1148  double xc, yc, zc, xm, ym, zm;
1149  int idof, sideNumBc, bcId, bcDofId;
1150  IntArray sideBcDofId, dofIdArray, *loadArray;
1151  FloatMatrix *lcs;
1152  bool hasBc;
1153 
1154  element->giveElementDofIDMask(dofIdArray);
1155 
1156  for ( int inode = startNode; inode <= endNode; inode++ ) {
1157  xc = corner [ inode - 1 ]->at(1);
1158  yc = corner [ inode - 1 ]->at(2);
1159  zc = corner [ inode - 1 ]->at(3);
1160 
1161  xm = midNode.at(1);
1162  ym = midNode.at(2);
1163  zm = midNode.at(3);
1164 
1165  Node *node = element->giveNode(inode);
1166 
1167  if ( node->giveNumberOfDofs() != dofs ) {
1168  OOFEM_ERROR("Dof mismatch");
1169  }
1170 
1171  connectivity = refinedElement->giveFineNodeArray(inode);
1172  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1173  hasBc = refinedElement->giveBcDofArray1D(inode, element, sideBcDofId, sideNumBc, tStep);
1174 
1175  pos = 1;
1176  u = 0.0;
1177  for ( int m = 0; m < level + 2; m++, u += du, pos++ ) {
1178  nd = connectivity->at(pos);
1179  if ( localNodeIdArray.at(nd) == 0 ) {
1180  DynamicInputRecord *ir = new DynamicInputRecord();
1181 
1182  localNodeIdArray.at(nd) = ++localNodeId;
1183  globalNodeIdArray.at(localNodeId) = nd;
1184 
1185  x = xc * ( 1.0 - u ) + xm * u;
1186  y = yc * ( 1.0 - u ) + ym * u;
1187  z = zc * ( 1.0 - u ) + zm * u;
1188 
1189  FloatArray coord = {x, y, z};
1190  newNodes.push_back(coord);
1191 
1192  ir->setRecordKeywordField(_IFT_Node_Name, localNodeId);
1193  ir->setField(coord, _IFT_Node_coords);
1194 
1195  if ( ( lcs = node->giveLocalCoordinateTriplet() ) != NULL ) {
1196  FloatArray lcs_vec = {lcs->at(1, 1), lcs->at(1, 2), lcs->at(1, 3),
1197  lcs->at(2, 1), lcs->at(2, 2), lcs->at(2, 3)};
1198  ir->setField(lcs_vec, _IFT_Node_lcs);
1199  }
1200 
1201  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1202  bc = 0;
1203  if ( nodeId == 0 ) {
1204  if ( m == 0 && boundary.at(1) == 0 ) {
1205  bc = 1;
1206  }
1207  } else {
1208  if ( m == level + 1 ) {
1209  bc = 1;
1210  }
1211  }
1212 
1213 #ifdef EXACT_ERROR
1214  if ( wholeFlag == true ) {
1215  bc = 0;
1216  }
1217 
1218 #endif
1219 
1220  // jak jsou razeny stupne volnosti v bc a ic (je-li jich jiny pocet, nez default dofs)
1221  // razeni je zrejme shodne s razenim dof v dofmanageru
1222 
1223  if ( bc == 1 ) {
1224  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
1225  IntArray bcs, dofids;
1226  for ( Dof *dof: *node ) {
1227  bcs.followedBy(++localBcId);
1228  dofids.followedBy(dof->giveDofID());
1229  }
1230  ir->setField(bcs, _IFT_DofManager_bc);
1231  ir->setField(dofids, _IFT_DofManager_dofidmask);
1232  }
1233  } else {
1234  if ( hasBc == true ) {
1235  // it is necessary to reproduce bc from coarse mesh
1236 
1237  if ( m == 0 ) { // at node
1238  IntArray bcs, dofids;
1239  for ( Dof *nodeDof: *node ) {
1240  bcDofId = 0;
1241  if ( nodeDof->hasBc(tStep) != 0 ) {
1242  bcDofId = nodeDof->giveBcId();
1243  }
1244  bcs.followedBy(bcDofId);
1245  dofids.followedBy(nodeDof->giveDofID());
1246  }
1247  ir->setField(bcs, _IFT_DofManager_bc);
1248  ir->setField(dofids, _IFT_DofManager_dofidmask);
1249 
1250  // copy node load
1251  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
1252  ir->setField(* loadArray, _IFT_DofManager_load);
1253  }
1254  } else {
1255  if ( sideNumBc != 0 ) {
1256  IntArray bcs;
1257 
1258  // I rely on the fact that bc dofs to be reproduced are ordered with respect to the dof ordering of the corner node
1259 
1260  bcId = 1;
1261  for ( idof = 1; idof <= dofs; idof++ ) {
1262  bcDofId = 0;
1263  if ( bcId <= sideNumBc ) {
1264  auto nodeDof = node->giveDofWithID( sideBcDofId.at(bcId) );
1265  if ( nodeDof->giveDofID() == dofIdArray.at(idof) ) {
1266  bcDofId = nodeDof->giveBcId();
1267  bcId++;
1268  }
1269  }
1270 
1271  bcs.at(idof) = bcDofId;
1272  }
1273  ir->setField(bcs, _IFT_DofManager_bc);
1274  //ir->setField(dofids, _IFT_DofManager_dofidmask);
1275  }
1276  }
1277  } else {
1278  // copy node load
1279 
1280  if ( m == 0 ) {
1281  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
1282  ir->setField(* loadArray, _IFT_DofManager_load);
1283  }
1284  }
1285  }
1286  }
1287 
1288  refinedReader.insertInputRecord(DataReader :: IR_dofmanRec, ir);
1289  }
1290  }
1291  }
1292 
1293  return;
1294  }
1295 
1296  std :: vector< FloatArray > newNodesVec( newNodes.begin(), newNodes.end() );
1297  if ( mode == HuertaErrorEstimatorInterface :: ElemMode ) {
1298  int csect, csect2;
1299  int nd1, nd2;
1300  IntArray boundaryLoadArray;
1301  bool hasLoad;
1302 
1303  csect = element->giveCrossSection()->giveNumber();
1304 
1305  for ( int inode = startNode; inode <= endNode; inode++ ) {
1306  connectivity = refinedElement->giveFineNodeArray(inode);
1307  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1308  hasLoad = refinedElement->giveBoundaryLoadArray1D(inode, element, boundaryLoadArray);
1309 
1310  for ( int m = 0; m < level + 1; m++ ) {
1311  DynamicInputRecord *ir = new DynamicInputRecord();
1312  localElemId++;
1313  ir->setRecordKeywordField(edgetype, localElemId);
1314 
1315  nd = m + 1;
1316 
1317  nd1 = localNodeIdArray.at( connectivity->at(nd) );
1318  nd2 = localNodeIdArray.at( connectivity->at(nd + 1) );
1319 
1321  csect2 = csect;
1322  if ( impCSect != 0 && impCSect == csect ) {
1323  FloatArray coordinates1, coordinates2;
1324  coordinates1 = newNodesVec [ nd1 - 1 ];
1325  coordinates2 = newNodesVec [ nd2 - 1 ];
1326 
1327  csect2 = impCSect;
1328  for ( int i = 1; i <= impPos.giveSize(); i++ ) {
1329  if ( impPos.at(i) < min( coordinates1.at(i), coordinates2.at(i) ) ) {
1330  csect2 = perCSect;
1331  break;
1332  }
1333 
1334  if ( impPos.at(i) > max( coordinates1.at(i), coordinates2.at(i) ) ) {
1335  csect2 = perCSect;
1336  break;
1337  }
1338  }
1339  }
1340 
1341  ir->setField(IntArray{nd1, nd2}, "nodes");
1342  ir->setField(csect2, "crosssect");
1343 
1344  // copy body and boundary loads
1345 
1346  if ( element->giveBodyLoadArray()->giveSize() != 0 ) {
1347  ir->setField(* element->giveBodyLoadArray(), "bodyloads");
1348  }
1349 
1350  if ( hasLoad == true ) {
1351  if ( boundaryLoadArray.giveSize() != 0 ) {
1352  ir->setField(boundaryLoadArray, "boundaryloads");
1353  }
1354  }
1355 
1356  refinedReader.insertInputRecord(DataReader :: IR_elemRec, ir);
1357  }
1358  }
1359  }
1360 
1361  if ( mode == HuertaErrorEstimatorInterface :: BCMode ) {
1362  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
1363  double u, du = 1.0 / ( level + 1 );
1364  double xc, yc, zc, xm, ym, zm;
1365  FloatArray globCoord(3);
1366  FloatMatrix Nmatrix;
1367  FloatArray uCoarse, uFine;
1368 
1369  MaterialMode mmode = element->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0)->giveMaterialMode();
1370 
1371  // create a fictitious integration point
1372  IntegrationRule ir(1, element);
1373  GaussPoint gp(&ir, 1, {}, 1.0, mmode);
1374 
1375  for ( int inode = startNode; inode <= endNode; inode++ ) {
1376  xc = corner [ inode - 1 ]->at(1);
1377  yc = corner [ inode - 1 ]->at(2);
1378  zc = corner [ inode - 1 ]->at(3);
1379 
1380  xm = midNode.at(1);
1381  ym = midNode.at(2);
1382  zm = midNode.at(3);
1383 
1384  connectivity = refinedElement->giveFineNodeArray(inode);
1385  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1386 
1387  // get corner displacements
1388  element->computeVectorOf(VM_Total, tStep, uCoarse);
1389 
1390  pos = 1;
1391  u = 0.0;
1392  for ( int m = 0; m < level + 2; m++, u += du, pos++ ) {
1393  nd = connectivity->at(pos);
1394  if ( localNodeIdArray.at(nd) > 0 ) {
1395  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
1396 
1397  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1398  bc = 0;
1399  if ( nodeId == 0 ) {
1400  if ( m == 0 && boundary.at(1) == 0 ) {
1401  bc = 1;
1402  }
1403  } else {
1404  if ( m == level + 1 ) {
1405  bc = 1;
1406  }
1407  }
1408 
1409 #ifdef EXACT_ERROR
1410  if ( wholeFlag == true ) {
1411  bc = 0;
1412  }
1413 
1414 #endif
1415 
1416  if ( bc == 1 ) {
1417  globCoord.at(1) = xc * ( 1.0 - u ) + xm * u;
1418  globCoord.at(2) = yc * ( 1.0 - u ) + ym * u;
1419  globCoord.at(3) = zc * ( 1.0 - u ) + zm * u;
1420 
1421  // this effectively rewrites the local coordinates of the fictitious integration point
1422  FloatArray lcoord;
1423  element->computeLocalCoordinates(lcoord, globCoord);
1424  gp.setNaturalCoordinates(lcoord);
1425 
1426  // get N matrix at the fictitious integration point
1427  this->HuertaErrorEstimatorI_computeNmatrixAt(&gp, Nmatrix);
1428  // get displacement at the fictitious integration point
1429  uFine.beProductOf(Nmatrix, uCoarse);
1430 
1431  // first loadtime function must be constant 1.0
1432  for ( int idof = 1; idof <= dofs; idof++ ) {
1433  DynamicInputRecord *ir = new DynamicInputRecord();
1434  ir->setRecordKeywordField(_IFT_BoundaryCondition_Name, ++localBcId);
1435  ir->setField(1, _IFT_GeneralBoundaryCondition_timeFunct);
1436  ir->setField(uFine.at(idof), _IFT_BoundaryCondition_PrescribedValue);
1437  refinedReader.insertInputRecord(DataReader :: IR_bcRec, ir);
1438  }
1439  }
1440  }
1441  }
1442  }
1443 
1444  } else {
1445  for ( int inode = startNode; inode <= endNode; inode++ ) {
1446  connectivity = refinedElement->giveFineNodeArray(inode);
1447  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1448 
1449  pos = 1;
1450  for ( int m = 0; m < level + 2; m++, pos++ ) {
1451  nd = connectivity->at(pos);
1452  if ( localNodeIdArray.at(nd) > 0 ) {
1453  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
1454 
1455  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1456  bc = 0;
1457  if ( nodeId == 0 ) {
1458  if ( m == 0 && boundary.at(1) == 0 ) {
1459  bc = 1;
1460  }
1461  } else {
1462  if ( m == level + 1 ) {
1463  bc = 1;
1464  }
1465  }
1466 
1467 #ifdef EXACT_ERROR
1468  if ( wholeFlag == true ) {
1469  bc = 0;
1470  }
1471 
1472 #endif
1473 
1474  if ( bc == 1 ) {
1475  for ( int idof = 1; idof <= dofs; idof++ ) {
1476  localBcId++;
1477  controlNode.at(localBcId) = -localNodeIdArray.at(nd);
1478  controlDof.at(localBcId) = idof;
1479  }
1480  }
1481  }
1482  }
1483  }
1484  }
1485  return;
1486  }
1487 }
1488 
1489 
1490 
1491 void
1492 HuertaErrorEstimatorInterface :: setupRefinedElementProblem2D(Element *element, RefinedElement *refinedElement,
1493  int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
1494  HuertaErrorEstimatorInterface :: SetupMode mode, TimeStep *tStep, int nodes,
1495  FloatArray **corner, FloatArray *midSide, FloatArray &midNode,
1496  int &localNodeId, int &localElemId, int &localBcId,
1497  IntArray &controlNode, IntArray &controlDof,
1498  HuertaErrorEstimator :: AnalysisMode aMode, const char *quadtype)
1499 {
1500  IntArray *connectivity, boundary(2);
1501  int startNode, endNode, inode, n, m, pos, nd, bc, dofs;
1502  Node *node;
1503 
1504  if ( nodeId != 0 ) {
1505  startNode = endNode = nodeId;
1506  } else {
1507  startNode = 1;
1508  endNode = nodes;
1509  }
1510 
1511  dofs = element->giveDomain()->giveDofManager(1)->giveNumberOfDofs();
1512 
1513  if ( mode == HuertaErrorEstimatorInterface :: CountMode ) {
1514  for ( inode = startNode; inode <= endNode; inode++ ) {
1515  localElemId += ( level + 1 ) * ( level + 1 );
1516 
1517  connectivity = refinedElement->giveFineNodeArray(inode);
1518  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1519 
1520  pos = 1;
1521  for ( n = 0; n < level + 2; n++ ) {
1522  for ( m = 0; m < level + 2; m++, pos++ ) {
1523  nd = connectivity->at(pos);
1524  if ( localNodeIdArray.at(nd) == 0 ) {
1525  localNodeIdArray.at(nd) = ++localNodeId;
1526 
1527  // count boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1528  bc = 0;
1529  if ( nodeId == 0 ) {
1530  if ( m == 0 && boundary.at(1) == 0 ) {
1531  bc = 1;
1532  }
1533 
1534  if ( n == 0 && boundary.at(2) == 0 ) {
1535  bc = 1;
1536  }
1537  } else {
1538  if ( n == level + 1 || m == level + 1 ) {
1539  bc = 1;
1540  }
1541 
1542  /* HUHU CHEATING */
1543  if ( bc == 0 ) {
1544  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
1545  if ( m == 0 && boundary.at(1) == 0 ) {
1546  bc = 1;
1547  }
1548 
1549  if ( n == 0 && boundary.at(2) == 0 ) {
1550  bc = 1;
1551  }
1552  }
1553  }
1554  }
1555 
1556 #ifdef EXACT_ERROR
1557  if ( wholeFlag == true ) {
1558  bc = 0;
1559  }
1560 
1561 #endif
1562 
1563  if ( bc == 1 ) {
1564  localBcId += dofs;
1565  }
1566  }
1567  }
1568  }
1569  }
1570 
1571  return;
1572  }
1573 
1574  if ( mode == HuertaErrorEstimatorInterface :: NodeMode ) {
1575  int s1, s2, idof, index;
1576  double x, y, z, u, v, du = 1.0 / ( level + 1 ), dv = 1.0 / ( level + 1 );
1577  double xc, yc, zc, xs1, ys1, zs1, xs2, ys2, zs2, xm, ym, zm;
1578  int bcId, bcDofId;
1579  std::vector< IntArray >sideBcDofIdList;
1580  IntArray sideNumBc(2), dofIdArray, * loadArray;
1581  bool hasBc;
1582  FloatMatrix *lcs;
1583 
1584  element->giveElementDofIDMask(dofIdArray);
1585 
1586  sideBcDofIdList.resize(2);
1587 
1588  for ( inode = startNode; inode <= endNode; inode++ ) {
1589  s1 = inode;
1590  if ( ( s2 = inode - 1 ) == 0 ) {
1591  s2 = nodes;
1592  }
1593 
1594  xc = corner [ inode - 1 ]->at(1);
1595  yc = corner [ inode - 1 ]->at(2);
1596  zc = corner [ inode - 1 ]->at(3);
1597 
1598  xs1 = midSide [ s1 - 1 ].at(1);
1599  ys1 = midSide [ s1 - 1 ].at(2);
1600  zs1 = midSide [ s1 - 1 ].at(3);
1601 
1602  xs2 = midSide [ s2 - 1 ].at(1);
1603  ys2 = midSide [ s2 - 1 ].at(2);
1604  zs2 = midSide [ s2 - 1 ].at(3);
1605 
1606  xm = midNode.at(1);
1607  ym = midNode.at(2);
1608  zm = midNode.at(3);
1609 
1610  node = element->giveNode(inode);
1611 
1612  if ( node->giveNumberOfDofs() != dofs ) {
1613  OOFEM_ERROR("Dof mismatch");
1614  }
1615 
1616  connectivity = refinedElement->giveFineNodeArray(inode);
1617  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1618  hasBc = refinedElement->giveBcDofArray2D(inode, element, sideBcDofIdList, sideNumBc, tStep);
1619 
1620  pos = 1;
1621  v = 0.0;
1622  for ( n = 0; n < level + 2; n++, v += dv ) {
1623  u = 0.0;
1624  for ( m = 0; m < level + 2; m++, u += du, pos++ ) {
1625  nd = connectivity->at(pos);
1626  if ( localNodeIdArray.at(nd) == 0 ) {
1627  DynamicInputRecord *ir = new DynamicInputRecord();
1628  ir->setRecordKeywordField("node", localNodeId);
1629 
1630  localNodeIdArray.at(nd) = ++localNodeId;
1631  globalNodeIdArray.at(localNodeId) = nd;
1632 
1633  x = ( xc * ( 1.0 - u ) + xs1 * u ) * ( 1.0 - v ) + ( xs2 * ( 1.0 - u ) + xm * u ) * v;
1634  y = ( yc * ( 1.0 - u ) + ys1 * u ) * ( 1.0 - v ) + ( ys2 * ( 1.0 - u ) + ym * u ) * v;
1635  z = ( zc * ( 1.0 - u ) + zs1 * u ) * ( 1.0 - v ) + ( zs2 * ( 1.0 - u ) + zm * u ) * v;
1636 
1637  ir->setField(FloatArray{x, y, z}, "coords");
1638 
1639  if ( ( lcs = node->giveLocalCoordinateTriplet() ) != NULL ) {
1640  FloatArray lcs_vec = {lcs->at(1, 1), lcs->at(1, 2), lcs->at(1, 3),
1641  lcs->at(2, 1), lcs->at(2, 2), lcs->at(2, 3)};
1642  ir->setField(lcs_vec, _IFT_Node_lcs);
1643  }
1644 
1645  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1646  bc = 0;
1647  if ( nodeId == 0 ) {
1648  if ( m == 0 && boundary.at(1) == 0 ) {
1649  bc = 1;
1650  }
1651 
1652  if ( n == 0 && boundary.at(2) == 0 ) {
1653  bc = 1;
1654  }
1655  } else {
1656  if ( n == level + 1 || m == level + 1 ) {
1657  bc = 1;
1658  }
1659 
1660  /* HUHU CHEATING */
1661  if ( bc == 0 ) {
1662  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
1663  if ( m == 0 && boundary.at(1) == 0 ) {
1664  bc = 1;
1665  }
1666 
1667  if ( n == 0 && boundary.at(2) == 0 ) {
1668  bc = 1;
1669  }
1670  }
1671  }
1672  }
1673 
1674 #ifdef EXACT_ERROR
1675  if ( wholeFlag == true ) {
1676  bc = 0;
1677  }
1678 
1679 #endif
1680 
1681  // jak jsou razeny stupne volnosti v bc a ic (je-li jich jiny pocet, nez default dofs)
1682  // razeni je zrejme shodne s razenim dof v dofmanageru
1683 
1684  if ( bc == 1 ) {
1685  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
1686  FloatArray bcs(dofs);
1687  for ( idof = 0; idof < dofs; idof++ ) {
1688  bcs.at(idof) = ++localBcId;
1689  }
1690  ir->setField(bcs, "bc");
1691  }
1692  } else {
1693  if ( hasBc == true && ( m == 0 || n == 0 ) ) {
1694  // it is necessary to reproduce bc from coarse mesh
1695 
1696  if ( m == 0 && n == 0 ) { // at node
1697  IntArray bcs, dofids;
1698 
1699  for ( Dof *nodeDof: *node ) {
1700  bcDofId = 0;
1701  if ( nodeDof->hasBc(tStep) != 0 ) {
1702  bcDofId = nodeDof->giveBcId();
1703  }
1704 
1705  bcs.followedBy(bcDofId);
1706  dofids.followedBy(nodeDof->giveDofID());
1707  }
1708  ir->setField(bcs, _IFT_DofManager_bc);
1709  ir->setField(dofids, _IFT_DofManager_dofidmask);
1710 
1711  // copy node load
1712 
1713  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
1714  ir->setField(* loadArray, "load");
1715  }
1716  } else {
1717  index = 0;
1718  if ( m == 0 && sideNumBc.at(1) != 0 ) {
1719  index = 1; // along next side
1720  }
1721 
1722  if ( n == 0 && sideNumBc.at(2) != 0 ) {
1723  index = 2; // along prev side
1724  }
1725 
1726  if ( index != 0 ) {
1727  FloatArray bcs(dofs);
1728 
1729  // I rely on the fact that bc dofs to be reproduced are ordered with respect to the dof ordering of the corner node
1730 
1731  const IntArray &sideBcDofId = sideBcDofIdList[index-1];
1732  bcId = 1;
1733  for ( idof = 1; idof <= dofs; idof++ ) {
1734  bcDofId = 0;
1735  if ( bcId <= sideNumBc.at(index) ) {
1736  auto nodeDof = node->giveDofWithID( sideBcDofId.at(bcId) );
1737  if ( nodeDof->giveDofID() == dofIdArray.at(idof) ) {
1738  bcDofId = nodeDof->giveBcId();
1739  bcId++;
1740  }
1741  }
1742 
1743  bcs.at(idof) = bcDofId;
1744  }
1745  ir->setField(bcs, "bc");
1746  }
1747  }
1748  } else {
1749  // copy node load
1750 
1751  if ( m == 0 && n == 0 ) {
1752  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
1753  ir->setField(* loadArray, "load");
1754  }
1755  }
1756  }
1757  }
1758 
1759  refinedReader.insertInputRecord(DataReader :: IR_dofmanRec, ir);
1760  }
1761  }
1762  }
1763  }
1764 
1765  return;
1766  }
1767 
1768  if ( mode == HuertaErrorEstimatorInterface :: ElemMode ) {
1769  int csect, iside, index;
1770  int nd1, nd2, nd3, nd4;
1771  IntArray *loadArray;
1772  std::vector< IntArray >boundaryLoadList;
1773  bool hasLoad;
1774 
1775  csect = element->giveCrossSection()->giveNumber();
1776 
1777  boundaryLoadList.resize(2);
1778 
1779  for ( inode = startNode; inode <= endNode; inode++ ) {
1780  connectivity = refinedElement->giveFineNodeArray(inode);
1781  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1782  hasLoad = refinedElement->giveBoundaryLoadArray2D(inode, element, boundaryLoadList);
1783 
1784  for ( n = 0; n < level + 1; n++ ) {
1785  for ( m = 0; m < level + 1; m++ ) {
1786  DynamicInputRecord *ir = new DynamicInputRecord();
1787  ir->setRecordKeywordField(quadtype, localElemId);
1788 
1789  localElemId++;
1790 
1791  nd = n * ( level + 2 ) + m + 1;
1792 
1793  nd1 = localNodeIdArray.at( connectivity->at(nd) );
1794  nd2 = localNodeIdArray.at( connectivity->at(nd + 1) );
1795  nd3 = localNodeIdArray.at( connectivity->at(nd + level + 3) );
1796  nd4 = localNodeIdArray.at( connectivity->at(nd + level + 2) );
1797 
1798  ir->setField(IntArray{nd1, nd2, nd3, nd4}, "nodes");
1799  ir->setField(csect, "crosssect");
1800 
1801  // copy body and boundary loads
1802 
1803  if ( ( loadArray = element->giveBodyLoadArray() )->giveSize() != 0 ) {
1804  ir->setField(* loadArray, "bodyloads");
1805  }
1806 
1807  // boundary load is not copied on non-boundary sides
1808 
1809  if ( hasLoad == true && ( m == 0 || n == 0 ) ) {
1810  if ( m == 0 && n == 0 ) {
1811  int loads = 0;
1812  for ( iside = 1; iside <= 2; iside++ ) {
1813  if ( boundary.at(iside) == 0 ) {
1814  continue;
1815  }
1816 
1817  loads += boundaryLoadList[iside-1].giveSize();
1818  }
1819 
1820  if ( loads != 0 ) {
1821  for ( iside = 1; iside <= 2; iside++ ) {
1822  if ( boundary.at(iside) == 0 ) {
1823  continue;
1824  }
1825 
1826  if ( boundaryLoadList[iside-1].giveSize() == 0 ) {
1827  continue;
1828  }
1829 
1830  ir->setField(boundaryLoadList[iside-1], "boundaryloads");
1831  }
1832  }
1833  } else {
1834  index = 0;
1835  if ( m == 0 && boundary.at(1) != 0 && boundaryLoadList[0].giveSize() != 0 ) {
1836  index = 1;
1837  }
1838 
1839  if ( n == 0 && boundary.at(2) != 0 && boundaryLoadList[1].giveSize() != 0 ) {
1840  index = 2;
1841  }
1842 
1843  if ( index != 0 ) {
1844  ir->setField(boundaryLoadList[index-1], "boundaryloads");
1845  }
1846  }
1847  }
1848 
1849  refinedReader.insertInputRecord(DataReader :: IR_elemRec, ir);
1850  }
1851  }
1852  }
1853 
1854  return;
1855  }
1856 
1857  if ( mode == HuertaErrorEstimatorInterface :: BCMode ) {
1858  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
1859  int s1, s2, idof;
1860  double u, v, du = 1.0 / ( level + 1 ), dv = 1.0 / ( level + 1 );
1861  double xc, yc, zc, xs1, ys1, zs1, xs2, ys2, zs2, xm, ym, zm;
1862  GaussPoint *gp;
1863  FloatArray globCoord(3);
1864  FloatMatrix Nmatrix;
1865  FloatArray uCoarse, uFine;
1866 
1867  MaterialMode mmode = element->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0)->giveMaterialMode();
1868 
1869  // create a fictitious integration point
1870  IntegrationRule ir(0, element);
1871  gp = new GaussPoint( &ir, 1, {}, 1.0, mmode);
1872 
1873  for ( inode = startNode; inode <= endNode; inode++ ) {
1874  s1 = inode;
1875  if ( ( s2 = inode - 1 ) == 0 ) {
1876  s2 = nodes;
1877  }
1878 
1879  xc = corner [ inode - 1 ]->at(1);
1880  yc = corner [ inode - 1 ]->at(2);
1881  zc = corner [ inode - 1 ]->at(3);
1882 
1883  xs1 = midSide [ s1 - 1 ].at(1);
1884  ys1 = midSide [ s1 - 1 ].at(2);
1885  zs1 = midSide [ s1 - 1 ].at(3);
1886 
1887  xs2 = midSide [ s2 - 1 ].at(1);
1888  ys2 = midSide [ s2 - 1 ].at(2);
1889  zs2 = midSide [ s2 - 1 ].at(3);
1890 
1891  xm = midNode.at(1);
1892  ym = midNode.at(2);
1893  zm = midNode.at(3);
1894 
1895  connectivity = refinedElement->giveFineNodeArray(inode);
1896  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1897 
1898  // get corner displacements
1899  element->computeVectorOf(VM_Total, tStep, uCoarse);
1900 
1901  pos = 1;
1902  v = 0.0;
1903  for ( n = 0; n < level + 2; n++, v += dv ) {
1904  u = 0.0;
1905  for ( m = 0; m < level + 2; m++, u += du, pos++ ) {
1906  nd = connectivity->at(pos);
1907  if ( localNodeIdArray.at(nd) > 0 ) {
1908  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
1909 
1910  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1911  bc = 0;
1912  if ( nodeId == 0 ) {
1913  if ( m == 0 && boundary.at(1) == 0 ) {
1914  bc = 1;
1915  }
1916 
1917  if ( n == 0 && boundary.at(2) == 0 ) {
1918  bc = 1;
1919  }
1920  } else {
1921  if ( n == level + 1 || m == level + 1 ) {
1922  bc = 1;
1923  }
1924 
1925  /* HUHU CHEATING */
1926  if ( bc == 0 ) {
1927  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
1928  if ( m == 0 && boundary.at(1) == 0 ) {
1929  bc = 1;
1930  }
1931 
1932  if ( n == 0 && boundary.at(2) == 0 ) {
1933  bc = 1;
1934  }
1935  }
1936  }
1937  }
1938 
1939 #ifdef EXACT_ERROR
1940  if ( wholeFlag == true ) {
1941  bc = 0;
1942  }
1943 
1944 #endif
1945 
1946  if ( bc == 1 ) {
1947  globCoord.at(1) = ( xc * ( 1.0 - u ) + xs1 * u ) * ( 1.0 - v ) + ( xs2 * ( 1.0 - u ) + xm * u ) * v;
1948  globCoord.at(2) = ( yc * ( 1.0 - u ) + ys1 * u ) * ( 1.0 - v ) + ( ys2 * ( 1.0 - u ) + ym * u ) * v;
1949  globCoord.at(3) = ( zc * ( 1.0 - u ) + zs1 * u ) * ( 1.0 - v ) + ( zs2 * ( 1.0 - u ) + zm * u ) * v;
1950 
1951  // this effectively rewrites the local coordinates of the fictitious integration point
1952  FloatArray lcoord;
1953  element->computeLocalCoordinates(lcoord, globCoord);
1954  gp->setNaturalCoordinates(lcoord);
1955  // get N matrix at the fictitious integration point
1956  this->HuertaErrorEstimatorI_computeNmatrixAt(gp, Nmatrix);
1957  // get displacement at the fictitious integration point
1958  uFine.beProductOf(Nmatrix, uCoarse);
1959 
1960  // first loadtime function must be constant 1.0
1961  for ( idof = 1; idof <= dofs; idof++ ) {
1962  DynamicInputRecord *ir = new DynamicInputRecord();
1963  ir->setRecordKeywordField("BoundaryCondition", ++localBcId);
1964  ir->setField(1, "Function");
1965  ir->setField(uFine.at(idof), "prescribedvalue");
1966  refinedReader.insertInputRecord(DataReader :: IR_bcRec, ir);
1967  }
1968  }
1969  }
1970  }
1971  }
1972  }
1973 
1974  delete gp;
1975  } else {
1976  int idof;
1977 
1978  for ( inode = startNode; inode <= endNode; inode++ ) {
1979  connectivity = refinedElement->giveFineNodeArray(inode);
1980  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
1981 
1982  pos = 1;
1983  for ( n = 0; n < level + 2; n++ ) {
1984  for ( m = 0; m < level + 2; m++, pos++ ) {
1985  nd = connectivity->at(pos);
1986  if ( localNodeIdArray.at(nd) > 0 ) {
1987  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
1988 
1989  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
1990  bc = 0;
1991  if ( nodeId == 0 ) {
1992  if ( m == 0 && boundary.at(1) == 0 ) {
1993  bc = 1;
1994  }
1995 
1996  if ( n == 0 && boundary.at(2) == 0 ) {
1997  bc = 1;
1998  }
1999  } else {
2000  if ( n == level + 1 || m == level + 1 ) {
2001  bc = 1;
2002  }
2003 
2004  /* HUHU CHEATING */
2005  if ( bc == 0 ) {
2006  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
2007  if ( m == 0 && boundary.at(1) == 0 ) {
2008  bc = 1;
2009  }
2010 
2011  if ( n == 0 && boundary.at(2) == 0 ) {
2012  bc = 1;
2013  }
2014  }
2015  }
2016  }
2017 
2018 #ifdef EXACT_ERROR
2019  if ( wholeFlag == true ) {
2020  bc = 0;
2021  }
2022 
2023 #endif
2024 
2025  if ( bc == 1 ) {
2026  for ( idof = 1; idof <= dofs; idof++ ) {
2027  localBcId++;
2028  controlNode.at(localBcId) = -localNodeIdArray.at(nd);
2029  controlDof.at(localBcId) = idof;
2030  }
2031  }
2032  }
2033  }
2034  }
2035  }
2036  }
2037 
2038  return;
2039  }
2040 }
2041 
2042 
2043 
2044 void
2045 HuertaErrorEstimatorInterface :: setupRefinedElementProblem3D(Element *element, RefinedElement *refinedElement,
2046  int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
2047  HuertaErrorEstimatorInterface :: SetupMode mode, TimeStep *tStep, int nodes,
2048  FloatArray **corner, FloatArray *midSide, FloatArray *midFace, FloatArray &midNode,
2049  int &localNodeId, int &localElemId, int &localBcId,
2050  int hexaSideNode [ 1 ] [ 3 ], int hexaFaceNode [ 1 ] [ 3 ],
2051  IntArray &controlNode, IntArray &controlDof,
2052  HuertaErrorEstimator :: AnalysisMode aMode, const char *hexatype)
2053 {
2054  IntArray *connectivity, boundary(3);
2055  int startNode, endNode, inode, k, n, m, pos, nd, bc, dofs;
2056  Node *node;
2057 
2058  if ( nodeId != 0 ) {
2059  startNode = endNode = nodeId;
2060  } else {
2061  startNode = 1;
2062  endNode = nodes;
2063  }
2064 
2065  dofs = element->giveDomain()->giveDofManager(1)->giveNumberOfDofs();
2066 
2067  if ( mode == HuertaErrorEstimatorInterface :: CountMode ) {
2068 #ifdef DEBUG
2069  if ( nodes / 2 * 2 != nodes ) {
2070  abort();
2071  }
2072 
2073  // OOFEM_ERROR("unexpected situation");
2074 
2075  /* number of internal quad faces = nodes * 3 / 2;
2076  *
2077  * at each node there are 3 internal quad faces and each internal quad face is shared by two internal hexas;
2078  * it is clear that this does not work for element with odd number of nodes such as pyramid;
2079  * however, pyramid is not decomposable into hexas and therefore should not provide this interface */
2080 #endif
2081 
2082  for ( inode = startNode; inode <= endNode; inode++ ) {
2083  localElemId += ( level + 1 ) * ( level + 1 ) * ( level + 1 );
2084 
2085  connectivity = refinedElement->giveFineNodeArray(inode);
2086  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
2087 
2088  pos = 1;
2089  for ( k = 0; k < level + 2; k++ ) {
2090  for ( n = 0; n < level + 2; n++ ) {
2091  for ( m = 0; m < level + 2; m++, pos++ ) {
2092  nd = connectivity->at(pos);
2093  if ( localNodeIdArray.at(nd) == 0 ) {
2094  localNodeIdArray.at(nd) = ++localNodeId;
2095 
2096  // count boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
2097  bc = 0;
2098  if ( nodeId == 0 ) {
2099  if ( m == 0 && boundary.at(1) == 0 ) {
2100  bc = 1;
2101  }
2102 
2103  if ( n == 0 && boundary.at(2) == 0 ) {
2104  bc = 1;
2105  }
2106 
2107  if ( k == 0 && boundary.at(3) == 0 ) {
2108  bc = 1;
2109  }
2110  } else {
2111  if ( k == level + 1 || n == level + 1 || m == level + 1 ) {
2112  bc = 1;
2113  }
2114 
2115  /* HUHU CHEATING */
2116  if ( bc == 0 ) {
2117  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
2118  if ( m == 0 && boundary.at(1) == 0 ) {
2119  bc = 1;
2120  }
2121 
2122  if ( n == 0 && boundary.at(2) == 0 ) {
2123  bc = 1;
2124  }
2125 
2126  if ( k == 0 && boundary.at(3) == 0 ) {
2127  bc = 1;
2128  }
2129  }
2130  }
2131  }
2132 
2133 #ifdef EXACT_ERROR
2134  if ( wholeFlag == true ) {
2135  bc = 0;
2136  }
2137 
2138 #endif
2139 
2140  if ( bc == 1 ) {
2141  localBcId += dofs;
2142  }
2143  }
2144  }
2145  }
2146  }
2147  }
2148 
2149  return;
2150  }
2151 
2152  if ( mode == HuertaErrorEstimatorInterface :: NodeMode ) {
2153  int s1, s2, s3, f1, f2, f3, idof, index;
2154  double x, y, z, u, v, w, du = 1.0 / ( level + 1 ), dv = 1.0 / ( level + 1 ), dw = 1.0 / ( level + 1 );
2155  double xc, yc, zc, xm, ym, zm;
2156  double xs1, ys1, zs1, xs2, ys2, zs2, xs3, ys3, zs3;
2157  double xf1, yf1, zf1, xf2, yf2, zf2, xf3, yf3, zf3;
2158  int bcId, bcDofId;
2159  std::vector < IntArray >sideBcDofIdList, faceBcDofIdList;
2160  IntArray sideNumBc(3), faceNumBc(3), dofIdArray, * loadArray;
2161  bool hasBc;
2162  FloatMatrix *lcs;
2163 
2164  element->giveElementDofIDMask(dofIdArray);
2165 
2166  sideBcDofIdList.resize(3);
2167  faceBcDofIdList.resize(3);
2168 
2169  for ( inode = startNode; inode <= endNode; inode++ ) {
2170  s1 = hexaSideNode [ inode - 1 ] [ 0 ];
2171  s2 = hexaSideNode [ inode - 1 ] [ 1 ];
2172  s3 = hexaSideNode [ inode - 1 ] [ 2 ];
2173  f1 = hexaFaceNode [ inode - 1 ] [ 0 ];
2174  f2 = hexaFaceNode [ inode - 1 ] [ 1 ];
2175  f3 = hexaFaceNode [ inode - 1 ] [ 2 ];
2176 
2177  xc = corner [ inode - 1 ]->at(1);
2178  yc = corner [ inode - 1 ]->at(2);
2179  zc = corner [ inode - 1 ]->at(3);
2180 
2181  xs1 = midSide [ s1 - 1 ].at(1);
2182  ys1 = midSide [ s1 - 1 ].at(2);
2183  zs1 = midSide [ s1 - 1 ].at(3);
2184 
2185  xs2 = midSide [ s2 - 1 ].at(1);
2186  ys2 = midSide [ s2 - 1 ].at(2);
2187  zs2 = midSide [ s2 - 1 ].at(3);
2188 
2189  xs3 = midSide [ s3 - 1 ].at(1);
2190  ys3 = midSide [ s3 - 1 ].at(2);
2191  zs3 = midSide [ s3 - 1 ].at(3);
2192 
2193  xf1 = midFace [ f1 - 1 ].at(1);
2194  yf1 = midFace [ f1 - 1 ].at(2);
2195  zf1 = midFace [ f1 - 1 ].at(3);
2196 
2197  xf2 = midFace [ f2 - 1 ].at(1);
2198  yf2 = midFace [ f2 - 1 ].at(2);
2199  zf2 = midFace [ f2 - 1 ].at(3);
2200 
2201  xf3 = midFace [ f3 - 1 ].at(1);
2202  yf3 = midFace [ f3 - 1 ].at(2);
2203  zf3 = midFace [ f3 - 1 ].at(3);
2204 
2205  xm = midNode.at(1);
2206  ym = midNode.at(2);
2207  zm = midNode.at(3);
2208 
2209  node = element->giveNode(inode);
2210 
2211  if ( node->giveNumberOfDofs() != dofs ) {
2212  OOFEM_ERROR("Dof mismatch");
2213  }
2214 
2215  connectivity = refinedElement->giveFineNodeArray(inode);
2216  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
2217  hasBc = refinedElement->giveBcDofArray3D(inode, element, sideBcDofIdList, sideNumBc,
2218  faceBcDofIdList, faceNumBc, tStep);
2219  pos = 1;
2220  w = 0.0;
2221  for ( k = 0; k < level + 2; k++, w += dw ) {
2222  v = 0.0;
2223  for ( n = 0; n < level + 2; n++, v += dv ) {
2224  u = 0.0;
2225  for ( m = 0; m < level + 2; m++, u += du, pos++ ) {
2226  nd = connectivity->at(pos);
2227  if ( localNodeIdArray.at(nd) == 0 ) {
2228  DynamicInputRecord *ir = new DynamicInputRecord();
2229  ir->setRecordKeywordField("node", localNodeId);
2230 
2231  localNodeIdArray.at(nd) = ++localNodeId;
2232  globalNodeIdArray.at(localNodeId) = nd;
2233 
2234  x = ( ( xc * ( 1.0 - u ) + xs1 * u ) * ( 1.0 - v ) + ( xs2 * ( 1.0 - u ) + xf1 * u ) * v ) * ( 1.0 - w )
2235  + ( ( xs3 * ( 1.0 - u ) + xf2 * u ) * ( 1.0 - v ) + ( xf3 * ( 1.0 - u ) + xm * u ) * v ) * w;
2236  y = ( ( yc * ( 1.0 - u ) + ys1 * u ) * ( 1.0 - v ) + ( ys2 * ( 1.0 - u ) + yf1 * u ) * v ) * ( 1.0 - w )
2237  + ( ( ys3 * ( 1.0 - u ) + yf2 * u ) * ( 1.0 - v ) + ( yf3 * ( 1.0 - u ) + ym * u ) * v ) * w;
2238  z = ( ( zc * ( 1.0 - u ) + zs1 * u ) * ( 1.0 - v ) + ( zs2 * ( 1.0 - u ) + zf1 * u ) * v ) * ( 1.0 - w )
2239  + ( ( zs3 * ( 1.0 - u ) + zf2 * u ) * ( 1.0 - v ) + ( zf3 * ( 1.0 - u ) + zm * u ) * v ) * w;
2240 
2241  ir->setField(FloatArray{x, y, z}, "coords");
2242 
2243  if ( ( lcs = node->giveLocalCoordinateTriplet() ) != NULL ) {
2244  FloatArray lcs_vec = {lcs->at(1, 1), lcs->at(1, 2), lcs->at(1, 3),
2245  lcs->at(2, 1), lcs->at(2, 2), lcs->at(2, 3)};
2246  ir->setField(lcs_vec, _IFT_Node_lcs);
2247  }
2248 
2249  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
2250  bc = 0;
2251  if ( nodeId == 0 ) {
2252  if ( m == 0 && boundary.at(1) == 0 ) {
2253  bc = 1;
2254  }
2255 
2256  if ( n == 0 && boundary.at(2) == 0 ) {
2257  bc = 1;
2258  }
2259 
2260  if ( k == 0 && boundary.at(3) == 0 ) {
2261  bc = 1;
2262  }
2263  } else {
2264  if ( k == level + 1 || n == level + 1 || m == level + 1 ) {
2265  bc = 1;
2266  }
2267 
2268  /* HUHU CHEATING */
2269  if ( bc == 0 ) {
2270  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
2271  if ( m == 0 && boundary.at(1) == 0 ) {
2272  bc = 1;
2273  }
2274 
2275  if ( n == 0 && boundary.at(2) == 0 ) {
2276  bc = 1;
2277  }
2278 
2279  if ( k == 0 && boundary.at(3) == 0 ) {
2280  bc = 1;
2281  }
2282  }
2283  }
2284  }
2285 
2286 #ifdef EXACT_ERROR
2287  if ( wholeFlag == true ) {
2288  bc = 0;
2289  }
2290 
2291 #endif
2292 
2293  // jak jsou razeny stupne volnosti v bc a ic (je-li jich jiny pocet, nez default dofs)
2294  // razeni je zrejme shodne s razenim dof v dofmanageru
2295 
2296  if ( bc == 1 ) {
2297  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
2298  FloatArray bcs(dofs);
2299  for ( idof = 0; idof < dofs; idof++ ) {
2300  bcs.at(idof) = ++localBcId;
2301  }
2302  ir->setField(bcs, "bc");
2303  }
2304  } else {
2305  if ( hasBc == true && ( m == 0 || n == 0 || k == 0 ) ) {
2306  // it is necessary to reproduce bc from coarse mesh
2307 
2308  if ( m == 0 && n == 0 && k == 0 ) { // at node
2309  IntArray bcs, dofids;
2310 
2311  for ( Dof *nodeDof: *node ) {
2312  bcDofId = 0;
2313  if ( nodeDof->hasBc(tStep) != 0 ) {
2314  bcDofId = nodeDof->giveBcId();
2315  }
2316 
2317  bcs.followedBy(bcDofId);
2318  dofids.followedBy(nodeDof->giveDofID());
2319  }
2320  ir->setField(bcs, _IFT_DofManager_bc);
2321  ir->setField(dofids, _IFT_DofManager_dofidmask);
2322 
2323  // copy node load
2324 
2325  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
2326  ir->setField(* loadArray, "load");
2327  }
2328  } else {
2329  if ( ( m == 0 && n == 0 ) || ( m == 0 && k == 0 ) || ( n == 0 && k == 0 ) ) {
2330  index = 0;
2331  if ( n == 0 && k == 0 && sideNumBc.at(1) != 0 ) {
2332  index = 1;
2333  }
2334 
2335  if ( m == 0 && k == 0 && sideNumBc.at(2) != 0 ) {
2336  index = 2;
2337  }
2338 
2339  if ( m == 0 && n == 0 && sideNumBc.at(3) != 0 ) {
2340  index = 3;
2341  }
2342 
2343  if ( index != 0 ) {
2344  FloatArray bcs(dofs);
2345 
2346  // I rely on the fact that bc dofs to be reproduced are ordered with respect to the dof ordering of the corner node
2347 
2348  const IntArray &sideBcDofId = sideBcDofIdList[index-1];
2349  bcId = 1;
2350  for ( idof = 1; idof <= dofs; idof++ ) {
2351  bcDofId = 0;
2352  if ( bcId <= sideNumBc.at(index) ) {
2353  Dof *nodeDof = node->giveDofWithID( sideBcDofId.at(bcId) );
2354  if ( nodeDof->giveDofID() == dofIdArray.at(idof) ) {
2355  bcDofId = nodeDof->giveBcId();
2356  bcId++;
2357  }
2358  }
2359 
2360  bcs.at(idof) = bcDofId;
2361  }
2362  ir->setField(bcs, "bc");
2363  }
2364  } else {
2365  index = 0;
2366  if ( m == 0 && faceNumBc.at(1) != 0 ) {
2367  index = 1;
2368  }
2369 
2370  if ( n == 0 && faceNumBc.at(2) != 0 ) {
2371  index = 2;
2372  }
2373 
2374  if ( k == 0 && faceNumBc.at(3) != 0 ) {
2375  index = 3;
2376  }
2377 
2378  if ( index != 0 ) {
2379  FloatArray bcs(dofs);
2380 
2381  // I rely on the fact that bc dofs to be reproduced are ordered with respect to the dof ordering of the corner node
2382 
2383  const IntArray &faceBcDofId = faceBcDofIdList[index-1];
2384  bcId = 1;
2385  for ( idof = 1; idof <= dofs; idof++ ) {
2386  bcDofId = 0;
2387  if ( bcId <= faceNumBc.at(index) ) {
2388  Dof *nodeDof = node->giveDofWithID( faceBcDofId.at(bcId) );
2389  if ( nodeDof->giveDofID() == dofIdArray.at(idof) ) {
2390  bcDofId = nodeDof->giveBcId();
2391  bcId++;
2392  }
2393  }
2394 
2395  bcs.at(idof) = bcDofId;
2396  }
2397  ir->setField(bcs, "bc");
2398  }
2399  }
2400  }
2401  } else {
2402  // copy node load
2403 
2404  if ( m == 0 && n == 0 ) {
2405  if ( ( loadArray = node->giveLoadArray() )->giveSize() != 0 ) {
2406  ir->setField(* loadArray, "loads");
2407  }
2408  }
2409  }
2410  }
2411 
2412  refinedReader.insertInputRecord(DataReader :: IR_dofmanRec, ir);
2413  }
2414  }
2415  }
2416  }
2417  }
2418 
2419  return;
2420  }
2421 
2422  if ( mode == HuertaErrorEstimatorInterface :: ElemMode ) {
2423  int csect, iside;
2424  int nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd8;
2425  IntArray *loadArray;
2426  std::vector< IntArray >boundaryLoadList;
2427  bool hasLoad;
2428 
2429  csect = element->giveCrossSection()->giveNumber();
2430 
2431  boundaryLoadList.resize(3);
2432 
2433  for ( inode = startNode; inode <= endNode; inode++ ) {
2434  connectivity = refinedElement->giveFineNodeArray(inode);
2435  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
2436  hasLoad = refinedElement->giveBoundaryLoadArray3D(inode, element, boundaryLoadList);
2437 
2438  for ( k = 0; k < level + 1; k++ ) {
2439  for ( n = 0; n < level + 1; n++ ) {
2440  for ( m = 0; m < level + 1; m++ ) {
2441  DynamicInputRecord *ir = new DynamicInputRecord();
2442 
2443  localElemId++;
2444 
2445  nd = k * ( level + 2 ) * ( level + 2 ) + n * ( level + 2 ) + m + 1;
2446 
2447  nd1 = localNodeIdArray.at( connectivity->at(nd) );
2448  nd2 = localNodeIdArray.at( connectivity->at(nd + 1) );
2449  nd3 = localNodeIdArray.at( connectivity->at(nd + level + 3) );
2450  nd4 = localNodeIdArray.at( connectivity->at(nd + level + 2) );
2451 
2452  nd += ( level + 2 ) * ( level + 2 );
2453 
2454  nd5 = localNodeIdArray.at( connectivity->at(nd) );
2455  nd6 = localNodeIdArray.at( connectivity->at(nd + 1) );
2456  nd7 = localNodeIdArray.at( connectivity->at(nd + level + 3) );
2457  nd8 = localNodeIdArray.at( connectivity->at(nd + level + 2) );
2458 
2459  ir->setRecordKeywordField(hexatype, localElemId);
2460  ir->setField(IntArray{nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd8}, "nodes");
2461  ir->setField(csect, "crosssect");
2462 
2463  // copy body and boundary loads
2464 
2465  if ( ( loadArray = element->giveBodyLoadArray() )->giveSize() != 0 ) {
2466  ir->setField(* loadArray, "bodyloads");
2467  }
2468 
2469  // boundary load is not copied on non-boundary sides
2470 
2471  if ( hasLoad == true && ( m == 0 || n == 0 || k == 0 ) ) {
2472  int loads = 0;
2473  for ( iside = 1; iside <= 3; iside++ ) {
2474  if ( boundary.at(iside) == 0 ) {
2475  continue;
2476  }
2477 
2478  if ( m != 0 && iside == 1 ) {
2479  continue;
2480  }
2481 
2482  if ( n != 0 && iside == 2 ) {
2483  continue;
2484  }
2485 
2486  if ( k != 0 && iside == 3 ) {
2487  continue;
2488  }
2489 
2490  loads += boundaryLoadList[iside-1].giveSize();
2491  }
2492 
2493  if ( loads != 0 ) {
2494  for ( iside = 1; iside <= 3; iside++ ) {
2495  if ( boundary.at(iside) == 0 ) {
2496  continue;
2497  }
2498 
2499  if ( m != 0 && iside == 1 ) {
2500  continue;
2501  }
2502 
2503  if ( n != 0 && iside == 2 ) {
2504  continue;
2505  }
2506 
2507  if ( k != 0 && iside == 3 ) {
2508  continue;
2509  }
2510 
2511  if ( boundaryLoadList[iside-1].giveSize() == 0 ) {
2512  continue;
2513  }
2514 
2515  ir->setField(boundaryLoadList[iside-1], "boundaryloads");
2516  }
2517  }
2518  }
2519 
2520  refinedReader.insertInputRecord(DataReader :: IR_elemRec, ir);
2521  }
2522  }
2523  }
2524  }
2525 
2526  return;
2527  }
2528 
2529  if ( mode == HuertaErrorEstimatorInterface :: BCMode ) {
2530  if ( aMode == HuertaErrorEstimator :: HEE_linear ) {
2531  int s1, s2, s3, f1, f2, f3, idof;
2532  double u, v, w, du = 1.0 / ( level + 1 ), dv = 1.0 / ( level + 1 ), dw = 1.0 / ( level + 1 );
2533  double xc, yc, zc, xm, ym, zm;
2534  double xs1, ys1, zs1, xs2, ys2, zs2, xs3, ys3, zs3;
2535  double xf1, yf1, zf1, xf2, yf2, zf2, xf3, yf3, zf3;
2536  FloatArray globCoord(3);
2537  FloatMatrix Nmatrix;
2538  FloatArray uCoarse, uFine;
2539 
2540  MaterialMode mmode = element->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0)->giveMaterialMode();
2541 
2542  // create a fictitious integration point
2543  IntegrationRule irule(0, element);
2544  auto gp = new GaussPoint(&irule, 1, {}, 1.0, mmode);
2545  for ( inode = startNode; inode <= endNode; inode++ ) {
2546  s1 = hexaSideNode [ inode - 1 ] [ 0 ];
2547  s2 = hexaSideNode [ inode - 1 ] [ 1 ];
2548  s3 = hexaSideNode [ inode - 1 ] [ 2 ];
2549  f1 = hexaFaceNode [ inode - 1 ] [ 0 ];
2550  f2 = hexaFaceNode [ inode - 1 ] [ 1 ];
2551  f3 = hexaFaceNode [ inode - 1 ] [ 2 ];
2552 
2553  xc = corner [ inode - 1 ]->at(1);
2554  yc = corner [ inode - 1 ]->at(2);
2555  zc = corner [ inode - 1 ]->at(3);
2556 
2557  xs1 = midSide [ s1 - 1 ].at(1);
2558  ys1 = midSide [ s1 - 1 ].at(2);
2559  zs1 = midSide [ s1 - 1 ].at(3);
2560 
2561  xs2 = midSide [ s2 - 1 ].at(1);
2562  ys2 = midSide [ s2 - 1 ].at(2);
2563  zs2 = midSide [ s2 - 1 ].at(3);
2564 
2565  xs3 = midSide [ s3 - 1 ].at(1);
2566  ys3 = midSide [ s3 - 1 ].at(2);
2567  zs3 = midSide [ s3 - 1 ].at(3);
2568 
2569  xf1 = midFace [ f1 - 1 ].at(1);
2570  yf1 = midFace [ f1 - 1 ].at(2);
2571  zf1 = midFace [ f1 - 1 ].at(3);
2572 
2573  xf2 = midFace [ f2 - 1 ].at(1);
2574  yf2 = midFace [ f2 - 1 ].at(2);
2575  zf2 = midFace [ f2 - 1 ].at(3);
2576 
2577  xf3 = midFace [ f3 - 1 ].at(1);
2578  yf3 = midFace [ f3 - 1 ].at(2);
2579  zf3 = midFace [ f3 - 1 ].at(3);
2580 
2581  xm = midNode.at(1);
2582  ym = midNode.at(2);
2583  zm = midNode.at(3);
2584 
2585  connectivity = refinedElement->giveFineNodeArray(inode);
2586  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
2587 
2588  // get corner displacements
2589  element->computeVectorOf(VM_Total, tStep, uCoarse);
2590 
2591  pos = 1;
2592  w = 0.0;
2593  for ( k = 0; k < level + 2; k++, w += dw ) {
2594  v = 0.0;
2595  for ( n = 0; n < level + 2; n++, v += dv ) {
2596  u = 0.0;
2597  for ( m = 0; m < level + 2; m++, u += du, pos++ ) {
2598  nd = connectivity->at(pos);
2599  if ( localNodeIdArray.at(nd) > 0 ) {
2600  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
2601 
2602  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
2603  bc = 0;
2604  if ( nodeId == 0 ) {
2605  if ( m == 0 && boundary.at(1) == 0 ) {
2606  bc = 1;
2607  }
2608 
2609  if ( n == 0 && boundary.at(2) == 0 ) {
2610  bc = 1;
2611  }
2612 
2613  if ( k == 0 && boundary.at(3) == 0 ) {
2614  bc = 1;
2615  }
2616  } else {
2617  if ( k == level + 1 || n == level + 1 || m == level + 1 ) {
2618  bc = 1;
2619  }
2620 
2621  /* HUHU CHEATING */
2622  if ( bc == 0 ) {
2623  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
2624  if ( m == 0 && boundary.at(1) == 0 ) {
2625  bc = 1;
2626  }
2627 
2628  if ( n == 0 && boundary.at(2) == 0 ) {
2629  bc = 1;
2630  }
2631 
2632  if ( k == 0 && boundary.at(3) == 0 ) {
2633  bc = 1;
2634  }
2635  }
2636  }
2637  }
2638 
2639 #ifdef EXACT_ERROR
2640  if ( wholeFlag == true ) {
2641  bc = 0;
2642  }
2643 
2644 #endif
2645 
2646  if ( bc == 1 ) {
2647  globCoord.at(1) = ( ( xc * ( 1.0 - u ) + xs1 * u ) * ( 1.0 - v ) + ( xs2 * ( 1.0 - u ) + xf1 * u ) * v ) * ( 1.0 - w )
2648  + ( ( xs3 * ( 1.0 - u ) + xf2 * u ) * ( 1.0 - v ) + ( xf3 * ( 1.0 - u ) + xm * u ) * v ) * w;
2649  globCoord.at(2) = ( ( yc * ( 1.0 - u ) + ys1 * u ) * ( 1.0 - v ) + ( ys2 * ( 1.0 - u ) + yf1 * u ) * v ) * ( 1.0 - w )
2650  + ( ( ys3 * ( 1.0 - u ) + yf2 * u ) * ( 1.0 - v ) + ( yf3 * ( 1.0 - u ) + ym * u ) * v ) * w;
2651  globCoord.at(3) = ( ( zc * ( 1.0 - u ) + zs1 * u ) * ( 1.0 - v ) + ( zs2 * ( 1.0 - u ) + zf1 * u ) * v ) * ( 1.0 - w )
2652  + ( ( zs3 * ( 1.0 - u ) + zf2 * u ) * ( 1.0 - v ) + ( zf3 * ( 1.0 - u ) + zm * u ) * v ) * w;
2653 
2654  // this effectively rewrites the local coordinates of the fictitious integration point
2655  FloatArray lcoord;
2656  element->computeLocalCoordinates(lcoord, globCoord);
2657  gp->setNaturalCoordinates(lcoord);
2658  // get N matrix at the fictitious integration point
2659  this->HuertaErrorEstimatorI_computeNmatrixAt(gp, Nmatrix);
2660  // get displacement at the fictitious integration point
2661  uFine.beProductOf(Nmatrix, uCoarse);
2662 
2663  // first loadtime function must be constant 1.0
2664  for ( idof = 1; idof <= dofs; idof++ ) {
2665  DynamicInputRecord *ir = new DynamicInputRecord();
2666  ir->setRecordKeywordField("boundarycondition", ++localBcId);
2667  ir->setField(1, "Function");
2668  ir->setField(uFine.at(idof), "prescribedvalue");
2669  refinedReader.insertInputRecord(DataReader :: IR_bcRec, ir);
2670  }
2671  }
2672  }
2673  }
2674  }
2675  }
2676  }
2677 
2678  } else {
2679  int idof;
2680 
2681  for ( inode = startNode; inode <= endNode; inode++ ) {
2682  connectivity = refinedElement->giveFineNodeArray(inode);
2683  refinedElement->giveBoundaryFlagArray(inode, element, boundary);
2684 
2685  pos = 1;
2686  for ( k = 0; k < level + 2; k++ ) {
2687  for ( n = 0; n < level + 2; n++ ) {
2688  for ( m = 0; m < level + 2; m++, pos++ ) {
2689  nd = connectivity->at(pos);
2690  if ( localNodeIdArray.at(nd) > 0 ) {
2691  localNodeIdArray.at(nd) = -localNodeIdArray.at(nd); // prevent repeated bc specification
2692 
2693  // setup boundary conditions (for element (nodeId == 0), for patch (nodeId != 0))
2694  bc = 0;
2695  if ( nodeId == 0 ) {
2696  if ( m == 0 && boundary.at(1) == 0 ) {
2697  bc = 1;
2698  }
2699 
2700  if ( n == 0 && boundary.at(2) == 0 ) {
2701  bc = 1;
2702  }
2703 
2704  if ( k == 0 && boundary.at(3) == 0 ) {
2705  bc = 1;
2706  }
2707  } else {
2708  if ( k == level + 1 || n == level + 1 || m == level + 1 ) {
2709  bc = 1;
2710  }
2711 
2712  /* HUHU CHEATING */
2713  if ( bc == 0 ) {
2714  if ( element->giveNode(nodeId)->giveParallelMode() == DofManager_shared ) {
2715  if ( m == 0 && boundary.at(1) == 0 ) {
2716  bc = 1;
2717  }
2718 
2719  if ( n == 0 && boundary.at(2) == 0 ) {
2720  bc = 1;
2721  }
2722 
2723  if ( k == 0 && boundary.at(3) == 0 ) {
2724  bc = 1;
2725  }
2726  }
2727  }
2728  }
2729 
2730 #ifdef EXACT_ERROR
2731  if ( wholeFlag == true ) {
2732  bc = 0;
2733  }
2734 
2735 #endif
2736 
2737  if ( bc == 1 ) {
2738  for ( idof = 1; idof <= dofs; idof++ ) {
2739  localBcId++;
2740  controlNode.at(localBcId) = -localNodeIdArray.at(nd);
2741  controlDof.at(localBcId) = idof;
2742  }
2743  }
2744  }
2745  }
2746  }
2747  }
2748  }
2749  }
2750 
2751  return;
2752  }
2753 }
2754 
2755 
2756 
2757 
2758 
2759 void
2760 HuertaErrorEstimator :: solveRefinedElementProblem(int elemId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
2761  TimeStep *tStep)
2762 {
2763  int contextFlag = 0;
2764  Element *element;
2765  RefinedElement *refinedElement;
2766  HuertaErrorEstimatorInterface *interface;
2767  EngngModel *problem, *refinedProblem;
2768  int localNodeId, localElemId, localBcId, localf;
2769  int mats, csects, loads, funcs, nlbarriers;
2770  int inode, idof, dofs, pos, ielem, size;
2771  IntArray dofIdArray;
2772  FloatArray nodeSolution, uCoarse, elementVector, patchVector, coarseVector;
2773  FloatArray elementError, patchError, coarseSolution;
2774  Domain *domain = this->domain, *refinedDomain;
2775  Node *node;
2776  TimeStep *refinedTStep;
2777  double coarseSol;
2778  FloatMatrix mat;
2779  EIPrimaryUnknownMapper mapper;
2780  Dof *nodeDof;
2781  double coeff, elementNorm, patchNorm, mixedNorm, eNorm = 0.0, uNorm = 0.0;
2782  IntArray controlNode, controlDof;
2783 
2784 #ifdef TIME_INFO
2785  Timer timer;
2786  double et_setup, et_init, et_solve, et_error;
2787  timer.startTimer();
2788 #endif
2789 
2790  element = domain->giveElement(elemId);
2791 
2792  if ( element->giveParallelMode() == Element_remote ) {
2793  this->skippedNelems++;
2794  this->eNorms.at(elemId) = 0.0;
2795  // uNormArray.at(elemId) = 0.0;
2796  return;
2797  }
2798 
2799  if ( this->skipRegion( element->giveRegionNumber() ) != 0 ) {
2800  this->skippedNelems++;
2801  this->eNorms.at(elemId) = 0.0;
2802  // uNormArray.at(elemId) = 0.0;
2803 
2804 #ifdef INFO
2805  // printf("\nElement no %d: skipped [step number %5d]\n", elemId, tStep -> giveNumber());
2806 #endif
2807 
2808  return;
2809  }
2810 
2811 #ifdef INFO
2812  OOFEM_LOG_DEBUG( "[%d] Element no %d: estimating error [step number %5d]\n",
2813  domain->giveEngngModel()->giveRank(), elemId, tStep->giveNumber() );
2814 #endif
2815 
2816  refinedElement = &this->refinedElementList[elemId-1];
2817  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
2818  if ( interface == NULL ) {
2819  OOFEM_ERROR("Element has no Huerta error estimator interface defined");
2820  }
2821 
2822  problem = domain->giveEngngModel();
2823 
2824  mats = domain->giveNumberOfMaterialModels();
2825  csects = domain->giveNumberOfCrossSectionModels();
2826  loads = domain->giveNumberOfBoundaryConditions();
2827  funcs = domain->giveNumberOfFunctions();
2828  nlbarriers = domain->giveNumberOfNonlocalBarriers();
2829 
2830  localNodeIdArray.zero();
2831  localNodeId = 0;
2832  localElemId = 0;
2833  localBcId = 0;
2834  localf = 0;
2835 
2836  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
2837  localNodeIdArray, globalNodeIdArray,
2838  HuertaErrorEstimatorInterface :: CountMode, tStep,
2839  localNodeId, localElemId, localBcId,
2840  controlNode, controlDof,
2841  this->mode);
2842 
2843  if ( this->mode == HEE_nlinear ) {
2844  controlDof.resize(localBcId);
2845  controlNode.resize(localBcId);
2846 
2847  localBcId = 0;
2848  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
2849  localNodeIdArray, globalNodeIdArray,
2850  HuertaErrorEstimatorInterface :: BCMode, tStep,
2851  localNodeId, localElemId, localBcId,
2852  controlNode, controlDof,
2853  this->mode);
2854 
2855  localBcId = 0;
2856  localf = 1;
2857  }
2858 
2859  setupRefinedProblemProlog("element", elemId, localNodeIdArray, localNodeId, localElemId,
2860  mats, csects, loads + localBcId, funcs + localf,
2861  controlNode, controlDof, tStep);
2862 
2863  globalNodeIdArray.resize(localNodeId);
2864 
2865  localNodeIdArray.zero();
2866  localNodeId = 0;
2867  localElemId = 0;
2868  localBcId = loads;
2869 
2870  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
2871  localNodeIdArray, globalNodeIdArray,
2872  HuertaErrorEstimatorInterface :: NodeMode, tStep,
2873  localNodeId, localElemId, localBcId,
2874  controlNode, controlDof,
2875  this->mode);
2876  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
2877  localNodeIdArray, globalNodeIdArray,
2878  HuertaErrorEstimatorInterface :: ElemMode, tStep,
2879  localNodeId, localElemId, localBcId,
2880  controlNode, controlDof,
2881  this->mode);
2882 
2883 
2884  setupRefinedProblemEpilog1(csects, mats, loads, nlbarriers);
2885 
2886  if ( this->mode == HEE_linear ) {
2887  localBcId = loads;
2888  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
2889  localNodeIdArray, globalNodeIdArray,
2890  HuertaErrorEstimatorInterface :: BCMode, tStep,
2891  localNodeId, localElemId, localBcId,
2892  controlNode, controlDof,
2893  this->mode);
2894  }
2895 
2896  setupRefinedProblemEpilog2(funcs);
2897 
2898 #ifdef TIME_INFO
2899  timer.stopTimer();
2900  et_setup = timer.getUtime();
2901 #endif
2902 
2903  dofs = domain->giveDofManager(1)->giveNumberOfDofs();
2904  domain->giveElement(1)->giveElementDofIDMask(dofIdArray);
2905 
2906 #ifdef USE_INPUT_FILE
2907  std :: ostringstream fileName;
2908  fileName << "/home/dr/Huerta/element_" << elemId << ".in";
2909  refinedReader.writeToFile( fileName.str().c_str() );
2910 #endif
2911 
2912 #ifdef TIME_INFO
2913  timer.startTimer();
2914 #endif
2915  refinedProblem = InstanciateProblem(refinedReader, _processor, contextFlag);
2916  refinedReader.finish();
2917 #ifdef TIME_INFO
2918  timer.stopTimer();
2919  et_init = timer.getUtime();
2920 #endif
2921 
2922 #ifdef DEBUG
2923  refinedProblem->checkConsistency();
2924 #endif
2925 
2926  refinedDomain = refinedProblem->giveDomain(1);
2927 
2928  // solve the problem first and then map the coarse solution;
2929  // when mapping the coarse solution at first, tstep is NULL and it cannot be accessed;
2930  // this makes some overhead for nonlinear problems because the coarse solution is mapped twice
2931  // when initiating the solution and after the solution
2932 
2933 #ifdef TIME_INFO
2934  timer.startTimer();
2935 #endif
2936  if ( this->mode == HEE_linear ) {
2937  refinedProblem->solveYourself();
2938  refinedProblem->terminateAnalysis();
2939  } else {
2940  AdaptiveNonLinearStatic *prob = dynamic_cast< AdaptiveNonLinearStatic * >(refinedProblem);
2941  if ( prob ) {
2942  static_cast< AdaptiveNonLinearStatic * >(refinedProblem)->initializeAdaptiveFrom(problem);
2943  } else {
2944  OOFEM_ERROR("Refined problem must be of the type AdaptiveNonLinearStatic");
2945  }
2946  }
2947 
2948 #ifdef TIME_INFO
2949  timer.stopTimer();
2950  et_solve = timer.getUtime();
2951 #endif
2952 
2953 #ifdef TIME_INFO
2954  timer.startTimer();
2955 #endif
2956  refinedTStep = refinedProblem->giveCurrentStep();
2957 
2958  size = refinedDomain->giveNumberOfDofManagers() * dofs;
2959  coarseSolution.resize(size);
2960  elementError.resize(size);
2961  patchError.resize(size);
2962 
2963  // map coarse solution
2964  uCoarse.resize( refinedProblem->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );
2965  uCoarse.zero();
2966  mapper.mapAndUpdate(uCoarse, VM_Total, domain, refinedDomain, tStep);
2967 
2968  // get coarse solution and element and patch error (including BC !!!)
2969  pos = 1;
2970  for ( inode = 1; inode <= localNodeId; inode++ ) {
2971  node = refinedDomain->giveNode(inode);
2972  node->giveUnknownVector(nodeSolution, dofIdArray, VM_Total, refinedTStep);
2973  for ( idof = 1; idof <= dofs; idof++, pos++ ) {
2974  double sol = nodeSolution.at(idof);
2975  nodeDof = node->giveDofWithID(dofIdArray.at(idof));
2976  if ( nodeDof->hasBc(refinedTStep) == 0 ) {
2977  coarseSol = uCoarse.at( nodeDof->__giveEquationNumber() );
2978  } else {
2979  // coarse solution is identical with fine solution at BC
2980  coarseSol = sol;
2981  }
2982 
2983  coarseSolution.at(pos) = coarseSol;
2984  elementError.at(pos) = sol - coarseSol;
2985  patchError.at(pos) = primaryUnknownError.at( ( globalNodeIdArray.at(inode) - 1 ) * dofs + idof ) - coarseSol;
2986  }
2987  }
2988 
2989  /* enforce zero error on element and patch boundary (this is just for 1D !!!)
2990  * (for nonlinear problem there may be a nonzero error due to the tolerance) */
2991  /*
2992  * elementError.at(1) = 0.0;
2993  * elementError.at(this -> refineLevel + 3) = 0.0;
2994  *
2995  * patchError.at(this -> refineLevel + 2) = 0.0;
2996  */
2997  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
2998  FloatMatrix Nmatrix;
2999  FloatArray elementVectorGp, patchVectorGp, coarseVectorGp;
3000 
3001  eNorm = uNorm = 0.0;
3002  for ( ielem = 1; ielem <= localElemId; ielem++ ) {
3003  element = refinedDomain->giveElement(ielem);
3004  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3005 
3006  elementNorm = patchNorm = mixedNorm = 0.0;
3007  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
3008  double dV = element->computeVolumeAround(gp);
3009 
3010  interface->HuertaErrorEstimatorI_computeNmatrixAt(gp, Nmatrix);
3011 
3012  this->extractVectorFrom(element, elementError, elementVector, dofs, refinedTStep);
3013  elementVectorGp.beProductOf(Nmatrix, elementVector);
3014  this->extractVectorFrom(element, patchError, patchVector, dofs, refinedTStep);
3015  patchVectorGp.beProductOf(Nmatrix, patchVector);
3016  this->extractVectorFrom(element, coarseSolution, coarseVector, dofs, refinedTStep);
3017  coarseVectorGp.beProductOf(Nmatrix, coarseVector);
3018 
3019  elementNorm += elementVectorGp.computeSquaredNorm() * dV;
3020  mixedNorm += elementVectorGp.dotProduct(patchVectorGp) * dV;
3021  patchNorm += patchVectorGp.computeSquaredNorm() * dV;
3022  uNorm += coarseVectorGp.computeSquaredNorm() * dV;
3023  }
3024 
3025  if ( fabs(elementNorm) < 1.0e-30 ) {
3026  if ( elementNorm == 0.0 ) {
3027  coeff = 0.0;
3028  } else {
3029  if ( fabs(mixedNorm) > 1.0e6 * fabs(elementNorm) ) {
3030  OOFEM_ERROR("division by zero");
3031  }
3032 
3033  coeff = mixedNorm / elementNorm;
3034  }
3035  } else {
3036  coeff = mixedNorm / elementNorm;
3037  }
3038 
3039  eNorm += ( 1.0 + coeff * coeff ) * elementNorm + patchNorm - 2.0 * coeff * mixedNorm;
3040  }
3041  } else if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
3042  FloatArray tmpVector;
3043 
3044 #ifdef PRINT_FINE_ERROR
3045  OOFEM_LOG_DEBUG("\n");
3046  if ( exactFlag == true ) {
3047  OOFEM_LOG_DEBUG(" elem sub e_Error2 p_Error2 a_Error2 x_Error2 \n");
3048  OOFEM_LOG_DEBUG("------------------------------------------------------------------------------\n");
3049  } else {
3050  OOFEM_LOG_DEBUG(" elem sub e_Error2 p_Error2 a_Error2 \n");
3051  OOFEM_LOG_DEBUG("-------------------------------------------------------------\n");
3052  }
3053 
3054 #endif
3055 
3056  eNorm = uNorm = 0.0;
3057  for ( ielem = 1; ielem <= localElemId; ielem++ ) {
3058  element = refinedDomain->giveElement(ielem);
3059  element->giveCharacteristicMatrix(mat, STIFFNESS_TYPE, refinedTStep);
3060 
3061  this->extractVectorFrom(element, elementError, elementVector, dofs, refinedTStep);
3062  this->extractVectorFrom(element, patchError, patchVector, dofs, refinedTStep);
3063  this->extractVectorFrom(element, coarseSolution, coarseVector, dofs, refinedTStep);
3064 
3065  tmpVector.beProductOf(mat, elementVector);
3066  elementNorm = tmpVector.dotProduct(elementVector);
3067  mixedNorm = tmpVector.dotProduct(patchVector);
3068 
3069  tmpVector.beProductOf(mat, patchVector);
3070  patchNorm = tmpVector.dotProduct(patchVector);
3071 
3072  if ( fabs(elementNorm) < 1.0e-30 ) {
3073  if ( elementNorm == 0.0 ) {
3074  coeff = 0.0;
3075  } else {
3076  if ( fabs(mixedNorm) > 1.0e6 * fabs(elementNorm) ) {
3077  OOFEM_ERROR("division by zero");
3078  }
3079 
3080  coeff = mixedNorm / elementNorm;
3081  }
3082  } else {
3083  coeff = mixedNorm / elementNorm;
3084  }
3085 
3086  eNorm += ( 1.0 + coeff * coeff ) * elementNorm + patchNorm - 2.0 * coeff * mixedNorm;
3087 
3088  tmpVector.beProductOf(mat, coarseVector);
3089  uNorm += tmpVector.dotProduct(coarseVector);
3090 
3091 #ifdef PRINT_FINE_ERROR
3092  double pEnorm = coeff * coeff * elementNorm + patchNorm - 2.0 * coeff * mixedNorm;
3093  if ( exactFlag == false ) {
3094  OOFEM_LOG_DEBUG("%5d: %3d %15.8e %15.8e %15.8e\n",
3095  elemId, ielem, elementNorm, pEnorm, elementNorm + pEnorm);
3096  }
3097 
3098  #ifdef EXACT_ERROR
3099  else {
3100  OOFEM_LOG_DEBUG( "%5d: %3d %15.8e %15.8e %15.8e %15.8e\n",
3101  elemId, ielem, elementNorm, pEnorm, elementNorm + pEnorm, exactFineError.at(++finePos) );
3102  }
3103  #endif
3104 #endif
3105  }
3106 
3107 #ifdef PRINT_FINE_ERROR
3108  OOFEM_LOG_DEBUG("\n");
3109 #endif
3110  } else {
3111  OOFEM_ERROR("Unsupported norm type");
3112  }
3113 
3114  // update primaryUnknownError
3115  // (this is usefull for postprocessing only, but how to draw fine elements ?)
3116  // be carefull to not overwrite data needed in further calculations
3117  // pos = 1;
3118  // for(inode = 1; inode <= localNodeId; inode++){
3119  // for(idof = 1; idof <= dofs; idof++, pos++)primaryUnknownError.at((globalNodeIdArray.at(inode) - 1) * dofs + idof) =
3120  // elementError.at(pos) * (1.0 - coeff) + patchError.at(pos);
3121  // }
3122 
3123 
3124  double eeeNorm = 0.0;
3125 
3126 #ifdef HUHU
3127  int result;
3128  double sval, maxVal;
3129  int k;
3130 
3131  maxVal = 0.0;
3132  element = domain->giveElement(elemId);
3133  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
3134  result = element->giveIPValue(val, gp, IST_PrincipalDamageTensor, tStep);
3135  if ( result ) {
3136  sval = sqrt( dotProduct( val, val, val.giveSize() ) );
3137  maxVal = max(maxVal, sval);
3138  }
3139  }
3140 
3141  if ( maxVal > 0.25 ) {
3142  double rate, size = 1.0, currDensity;
3143 
3144  HuertaRemeshingCriteriaInterface *remeshInterface;
3145  remeshInterface = static_cast< HuertaRemeshingCriteriaInterface * >( element->giveInterface(HuertaRemeshingCriteriaInterfaceType) );
3146  if ( !remeshInterface ) {
3147  OOFEM_ERROR("element does not support HuertaRemeshingCriteriaInterface");
3148  }
3149 
3150  currDensity = remeshInterface->HuertaRemeshingCriteriaI_giveCharacteristicSize();
3151 
3152  rate = currDensity / size;
3153  if ( rate < 1.0 ) {
3154  rate = 1.0;
3155  }
3156 
3157  OOFEM_LOG_DEBUG("koko %d dam %e curr %e rate %e\n", elemId, maxVal, currDensity, rate);
3158 
3159  // eNorm *= rate * rate;
3160  eeeNorm = eNorm * rate;
3161  }
3162 
3163 #endif
3164 
3165  this->eNorms.at(elemId) = sqrt(eNorm);
3166  // uNormArray.at(elemId) = sqrt(uNorm);
3167 
3168  this->globalENorm += eNorm + eeeNorm;
3169  this->globalUNorm += uNorm;
3170 
3171 #ifdef TIME_INFO
3172  timer.stopTimer();
3173  et_error = timer.getUtime();
3174 
3175  OOFEM_LOG_DEBUG("HEE info: element %d: user time total %.2f s (setup %.2f s, init %.2f s, solve %.2f s, error %.2f s)\n",
3176  elemId,
3177  et_setup + et_init + et_solve + et_error,
3178  et_setup,
3179  et_init,
3180  et_solve,
3181  et_error);
3182 #endif
3183 
3184  delete refinedProblem;
3185 }
3186 
3187 
3188 
3189 void
3190 HuertaErrorEstimator :: solveRefinedPatchProblem(int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
3191  TimeStep *tStep)
3192 {
3193  int contextFlag = 0;
3194  Element *element;
3195  RefinedElement *refinedElement;
3196  HuertaErrorEstimatorInterface *interface;
3197  EngngModel *problem, *refinedProblem;
3198  int localNodeId, localElemId, localBcId, localf;
3199  int mats, csects, loads, funcs, nlbarriers;
3200  int inode, elemId, ielem, elems, skipped = 0;
3201  const IntArray *con;
3202  int dofs, pos;
3203  IntArray dofIdArray;
3204  FloatArray nodeSolution;
3205  Domain *domain = this->domain, *refinedDomain;
3206  TimeStep *refinedTStep;
3207  ConnectivityTable *ct = domain->giveConnectivityTable();
3208  IntArray controlNode, controlDof;
3209 
3210 #ifdef TIME_INFO
3211  Timer timer;
3212  double et_setup, et_init, et_solve, et_error;
3213 
3214  timer.startTimer();
3215 #endif
3216 
3217  dofManagerParallelMode parMode = domain->giveDofManager(nodeId)->giveParallelMode();
3218  if ( parMode == DofManager_remote || parMode == DofManager_null ) {
3219  return;
3220  }
3221 
3222  con = ct->giveDofManConnectivityArray(nodeId);
3223  elems = con->giveSize();
3224 
3225  for ( ielem = 1; ielem <= elems; ielem++ ) {
3226  elemId = con->at(ielem);
3227  element = domain->giveElement(elemId);
3228 
3229  /* HUHU CHEATING */
3230  if ( element->giveParallelMode() == Element_remote ) {
3231  skipped++;
3232  continue;
3233  }
3234 
3235  if ( this->skipRegion( element->giveRegionNumber() ) != 0 ) {
3236  skipped++;
3237  }
3238  }
3239 
3240  if ( skipped == elems ) {
3241 #ifdef INFO
3242  // printf("\nPatch no %d: skipped [step number %5d]\n", nodeId, tStep -> giveNumber());
3243 #endif
3244 
3245  return;
3246  }
3247 
3248 #ifdef INFO
3249  OOFEM_LOG_INFO( "[%d] Patch no %d: estimating error [step number %5d]\n",
3250  domain->giveEngngModel()->giveRank(), nodeId, tStep->giveNumber() );
3251 #endif
3252 
3253  problem = domain->giveEngngModel();
3254 
3255  mats = domain->giveNumberOfMaterialModels();
3256  csects = domain->giveNumberOfCrossSectionModels();
3257  loads = domain->giveNumberOfBoundaryConditions();
3258  funcs = domain->giveNumberOfFunctions();
3259  nlbarriers = domain->giveNumberOfNonlocalBarriers();
3260 
3261  localNodeIdArray.zero();
3262  localNodeId = 0;
3263  localElemId = 0;
3264  localBcId = 0;
3265  localf = 0;
3266 
3267  for ( ielem = 1; ielem <= elems; ielem++ ) {
3268  elemId = con->at(ielem);
3269  element = domain->giveElement(elemId);
3270 
3271  /* HUHU CHEATING */
3272  if ( element->giveParallelMode() == Element_remote ) {
3273  continue;
3274  }
3275 
3276  refinedElement = &this->refinedElementList.at(elemId);
3277  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3278  if ( interface == NULL ) {
3279  OOFEM_ERROR("Element has no Huerta error estimator interface defined");
3280  }
3281 
3282  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3283  if ( element->giveNode(inode)->giveNumber() != nodeId ) {
3284  continue;
3285  }
3286 
3287  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, inode,
3288  localNodeIdArray, globalNodeIdArray,
3289  HuertaErrorEstimatorInterface :: CountMode, tStep,
3290  localNodeId, localElemId, localBcId,
3291  controlNode, controlDof,
3292  this->mode);
3293  break;
3294  }
3295  }
3296 
3297  if ( this->mode == HEE_nlinear ) {
3298  controlDof.resize(localBcId);
3299  controlNode.resize(localBcId);
3300 
3301  localBcId = 0;
3302  for ( ielem = 1; ielem <= elems; ielem++ ) {
3303  elemId = con->at(ielem);
3304  element = domain->giveElement(elemId);
3305 
3306  /* HUHU CHEATING */
3307  if ( element->giveParallelMode() == Element_remote ) {
3308  continue;
3309  }
3310 
3311  refinedElement = &this->refinedElementList.at(elemId);
3312  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3313 
3314  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3315  if ( element->giveNode(inode)->giveNumber() != nodeId ) {
3316  continue;
3317  }
3318 
3319  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, inode,
3320  localNodeIdArray, globalNodeIdArray,
3321  HuertaErrorEstimatorInterface :: BCMode, tStep,
3322  localNodeId, localElemId, localBcId,
3323  controlNode, controlDof,
3324  this->mode);
3325  break;
3326  }
3327  }
3328 
3329  localBcId = 0;
3330  localf = 1;
3331  }
3332 
3333  setupRefinedProblemProlog("patch", nodeId, localNodeIdArray, localNodeId, localElemId,
3334  mats, csects, loads + localBcId, funcs + localf,
3335  controlNode, controlDof, tStep);
3336 
3337  globalNodeIdArray.resize(localNodeId);
3338 
3339  localNodeIdArray.zero();
3340  localNodeId = 0;
3341  localElemId = 0;
3342  localBcId = loads;
3343 
3344  for ( ielem = 1; ielem <= elems; ielem++ ) {
3345  elemId = con->at(ielem);
3346  element = domain->giveElement(elemId);
3347 
3348  /* HUHU CHEATING */
3349  if ( element->giveParallelMode() == Element_remote ) {
3350  continue;
3351  }
3352 
3353  refinedElement = &this->refinedElementList.at(elemId);
3354  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3355 
3356  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3357  if ( element->giveNode(inode)->giveNumber() != nodeId ) {
3358  continue;
3359  }
3360 
3361  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, inode,
3362  localNodeIdArray, globalNodeIdArray,
3363  HuertaErrorEstimatorInterface :: NodeMode, tStep,
3364  localNodeId, localElemId, localBcId,
3365  controlNode, controlDof,
3366  this->mode);
3367  break;
3368  }
3369  }
3370 
3371  for ( ielem = 1; ielem <= elems; ielem++ ) {
3372  elemId = con->at(ielem);
3373  element = domain->giveElement(elemId);
3374 
3375  /* HUHU CHEATING */
3376  if ( element->giveParallelMode() == Element_remote ) {
3377  continue;
3378  }
3379 
3380  refinedElement = &this->refinedElementList.at(elemId);
3381  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3382 
3383  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3384  if ( element->giveNode(inode)->giveNumber() != nodeId ) {
3385  continue;
3386  }
3387 
3388  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, inode,
3389  localNodeIdArray, globalNodeIdArray,
3390  HuertaErrorEstimatorInterface :: ElemMode, tStep,
3391  localNodeId, localElemId, localBcId,
3392  controlNode, controlDof,
3393  this->mode);
3394  break;
3395  }
3396  }
3397 
3398  setupRefinedProblemEpilog1(csects, mats, loads, nlbarriers);
3399 
3400  if ( this->mode == HEE_linear ) {
3401  localBcId = loads;
3402  for ( ielem = 1; ielem <= elems; ielem++ ) {
3403  elemId = con->at(ielem);
3404  element = domain->giveElement(elemId);
3405 
3406  /* HUHU CHEATING */
3407  if ( element->giveParallelMode() == Element_remote ) {
3408  continue;
3409  }
3410 
3411  refinedElement = &this->refinedElementList.at(elemId);
3412  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3413 
3414  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3415  if ( element->giveNode(inode)->giveNumber() != nodeId ) {
3416  continue;
3417  }
3418 
3419  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, inode,
3420  localNodeIdArray, globalNodeIdArray,
3421  HuertaErrorEstimatorInterface :: BCMode, tStep,
3422  localNodeId, localElemId, localBcId,
3423  controlNode, controlDof,
3424  this->mode);
3425  break;
3426  }
3427  }
3428  }
3429 
3430  setupRefinedProblemEpilog2(funcs);
3431 
3432 #ifdef TIME_INFO
3433  timer.stopTimer();
3434  et_setup = timer.getUtime();
3435 #endif
3436 
3437  dofs = domain->giveDofManager(1)->giveNumberOfDofs();
3438  domain->giveElement(1)->giveElementDofIDMask(dofIdArray);
3439 
3440 #ifdef USE_INPUT_FILE
3441  std :: ostringstream fileName;
3442  fileName << "/home/dr/Huerta/patch_" << nodeId << ".in";
3443  refinedReader.writeToFile( fileName.str().c_str() );
3444 #endif
3445 
3446 #ifdef TIME_INFO
3447  timer.startTimer();
3448 #endif
3449  refinedProblem = InstanciateProblem(refinedReader, _processor, contextFlag);
3450  refinedReader.finish();
3451 #ifdef TIME_INFO
3452  timer.stopTimer();
3453  et_init = timer.getUtime();
3454 #endif
3455 
3456 #ifdef DEBUG
3457  refinedProblem->checkConsistency();
3458 #endif
3459 
3460  refinedDomain = refinedProblem->giveDomain(1);
3461 
3462 #ifdef TIME_INFO
3463  timer.startTimer();
3464 #endif
3465  if ( this->mode == HEE_linear ) {
3466  refinedProblem->solveYourself();
3467  refinedProblem->terminateAnalysis();
3468  } else {
3469  AdaptiveNonLinearStatic *prob = dynamic_cast< AdaptiveNonLinearStatic * >(refinedProblem);
3470  if ( prob ) {
3471  static_cast< AdaptiveNonLinearStatic * >(refinedProblem)->initializeAdaptiveFrom(problem);
3472  } else {
3473  OOFEM_ERROR("Refined problem must be of the type AdaptiveNonLinearStatic");
3474  }
3475  }
3476 
3477 #ifdef TIME_INFO
3478  timer.stopTimer();
3479  et_solve = timer.getUtime();
3480 #endif
3481 
3482  //fprintf(stdout, "\n");
3483 
3484 #ifdef TIME_INFO
3485  timer.startTimer();
3486 #endif
3487  refinedTStep = refinedProblem->giveCurrentStep();
3488 
3489  // store fine solution in primaryUnknownError
3490  for ( int inode = 1; inode <= localNodeId; inode++ ) {
3491  refinedDomain->giveNode(inode)->giveUnknownVector(nodeSolution, dofIdArray, VM_Total, refinedTStep);
3492  pos = globalNodeIdArray.at(inode);
3493  for ( int idof = 1; idof <= dofs; idof++ ) {
3494  primaryUnknownError.at( ( pos - 1 ) * dofs + idof ) = nodeSolution.at(idof);
3495  }
3496  }
3497 
3498 #ifdef TIME_INFO
3499  timer.stopTimer();
3500  et_error = timer.getUtime();
3501 
3502  OOFEM_LOG_DEBUG("HEE info: patch %d: user time total %.2f s (setup %.2f s, init %.2f s, solve %.2f s, error %.2f s)\n",
3503  nodeId,
3504  et_setup + et_init + et_solve + et_error,
3505  et_setup,
3506  et_init,
3507  et_solve,
3508  et_error);
3509 #endif
3510 
3511  delete refinedProblem;
3512 }
3513 
3514 
3515 #ifndef EXACT_ERROR
3516 void
3517 HuertaErrorEstimator :: solveRefinedWholeProblem(IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
3518  TimeStep *tStep) { }
3519 #else
3520 void
3521 HuertaErrorEstimator :: solveRefinedWholeProblem(IntArray &localNodeIdArray, IntArray &globalNodeIdArray,
3522  TimeStep *tStep)
3523 {
3524  int contextFlag = 0;
3525  Element *element;
3526  RefinedElement *refinedElement;
3527  HuertaErrorEstimatorInterface *interface;
3528  EngngModel *refinedProblem;
3529  int localNodeId, localElemId, localBcId, localf;
3530  int mats, csects, loads, funcs, nlbarriers;
3531  int inode, idof, dofs, elemId, ielem, elems, size;
3532  IntArray dofIdArray;
3533  FloatArray nodeSolution, uCoarse, errorVector, coarseVector, fineVector;
3534  FloatArray fineSolution, coarseSolution, errorSolution;
3535  Domain *domain = this->domain, *refinedDomain;
3536  Node *node;
3537  TimeStep *refinedTStep;
3538  FloatMatrix mat;
3539  EIPrimaryUnknownMapper mapper;
3540  Dof *nodeDof;
3541  IntArray controlNode, controlDof;
3542 
3543  #ifdef TIME_INFO
3544  Timer timer;
3545  double et_setup, et_init, et_solve, et_error;
3546 
3547  timer.startTimer();
3548  #endif
3549  #ifdef INFO
3550  OOFEM_LOG_INFO( "Whole 0: estimating error [step number %5d]\n", tStep->giveNumber() );
3551  #endif
3552 
3553  elems = domain->giveNumberOfElements();
3554 
3555  mats = domain->giveNumberOfMaterialModels();
3556  csects = domain->giveNumberOfCrossSectionModels();
3557  loads = domain->giveNumberOfBoundaryConditions();
3558  funcs = domain->giveNumberOfFunctions();
3559  nlbarriers = domain->giveNumberOfNonlocalBarriers();
3560 
3561  localNodeIdArray.zero();
3562  localNodeId = 0;
3563  localElemId = 0;
3564  localBcId = 0;
3565  localf = 0;
3566 
3567  for ( elemId = 1; elemId <= elems; elemId++ ) {
3568  element = domain->giveElement(elemId);
3569  refinedElement = &this->refinedElementList.at(elemId);
3570  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3571  if ( interface == NULL ) {
3572  OOFEM_ERROR("Element has no Huerta error estimator interface defined");
3573  }
3574 
3575  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
3576  localNodeIdArray, globalNodeIdArray,
3577  HuertaErrorEstimatorInterface :: CountMode, tStep,
3578  localNodeId, localElemId, localBcId,
3579  controlNode, controlDof,
3580  this->mode);
3581  }
3582 
3583  if ( this->mode == HEE_nlinear ) {
3584  controlDof.resize(localBcId);
3585  controlNode.resize(localBcId);
3586 
3587  localBcId = 0;
3588  for ( elemId = 1; elemId <= elems; elemId++ ) {
3589  element = domain->giveElement(elemId);
3590  refinedElement = &this->refinedElementList.at(elemId);
3591  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3592  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
3593  localNodeIdArray, globalNodeIdArray,
3594  HuertaErrorEstimatorInterface :: BCMode, tStep,
3595  localNodeId, localElemId, localBcId,
3596  controlNode, controlDof,
3597  this->mode);
3598  }
3599 
3600  localBcId = 0;
3601  localf = 1;
3602  }
3603 
3604  setupRefinedProblemProlog("whole", 0, localNodeIdArray, localNodeId, localElemId,
3605  mats, csects, loads + localBcId, funcs + localf,
3606  controlNode, controlDof, tStep);
3607 
3608  globalNodeIdArray.resize(localNodeId);
3609 
3610  localNodeIdArray.zero();
3611  localNodeId = 0;
3612  localElemId = 0;
3613  localBcId = loads;
3614 
3615  for ( elemId = 1; elemId <= elems; elemId++ ) {
3616  element = domain->giveElement(elemId);
3617  refinedElement = &this->refinedElementList.at(elemId);
3618  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3619  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
3620  localNodeIdArray, globalNodeIdArray,
3621  HuertaErrorEstimatorInterface :: NodeMode, tStep,
3622  localNodeId, localElemId, localBcId,
3623  controlNode, controlDof,
3624  this->mode);
3625  }
3626 
3627  for ( elemId = 1; elemId <= elems; elemId++ ) {
3628  element = domain->giveElement(elemId);
3629  refinedElement = &this->refinedElementList.at(elemId);
3630  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3631  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
3632  localNodeIdArray, globalNodeIdArray,
3633  HuertaErrorEstimatorInterface :: ElemMode, tStep,
3634  localNodeId, localElemId, localBcId,
3635  controlNode, controlDof,
3636  this->mode);
3637  }
3638 
3639  setupRefinedProblemEpilog1(csects, mats, loads, nlbarriers);
3640 
3641  if ( this->mode == HEE_linear ) {
3642  localBcId = loads;
3643  for ( elemId = 1; elemId <= elems; elemId++ ) {
3644  element = domain->giveElement(elemId);
3645  refinedElement = &this->refinedElementList.at(elemId);
3646  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3647  interface->HuertaErrorEstimatorI_setupRefinedElementProblem(refinedElement, this->refineLevel, 0,
3648  localNodeIdArray, globalNodeIdArray,
3649  HuertaErrorEstimatorInterface :: BCMode, tStep,
3650  localNodeId, localElemId, localBcId,
3651  controlNode, controlDof,
3652  this->mode);
3653  }
3654  }
3655 
3656  setupRefinedProblemEpilog2(funcs);
3657 
3658  #ifdef TIME_INFO
3659  timer.stopTimer();
3660  et_setup = timer.getUtime();
3661  #endif
3662 
3663  dofs = domain->giveDofManager(1)->giveNumberOfDofs();
3664  domain->giveElement(1)->giveElementDofIDMask(dofIdArray);
3665 
3666  #ifdef USE_INPUT_FILE
3667  std :: ostringstream fileName;
3668  fileName << "/home/dr/Huerta/whole_" << 0 << ".in";
3669  refinedReader.writeToFile( fileName.str().c_str() );
3670  #endif
3671 
3672  #ifdef TIME_INFO
3673  timer.startTimer();
3674  #endif
3675  refinedProblem = InstanciateProblem(refinedReader, _processor, contextFlag);
3676  #ifdef TIME_INFO
3677  timer.stopTimer();
3678  et_init = timer.getUtime();
3679  #endif
3680 
3681  #ifdef DEBUG
3682  refinedProblem->checkConsistency();
3683  #endif
3684 
3685  refinedDomain = refinedProblem->giveDomain(1);
3686 
3687  // solve the problem first and then map the coarse solution;
3688  // when mapping the coarse solution at first, tstep is NULL and it cannot be accessed;
3689 
3690  #ifdef TIME_INFO
3691  timer.startTimer();
3692  #endif
3693  refinedProblem->solveYourself();
3694  refinedProblem->terminateAnalysis();
3695  #ifdef TIME_INFO
3696  timer.stopTimer();
3697  et_solve = timer.getUtime();
3698  #endif
3699 
3700  //fprintf(stdout, "\n");
3701 
3702  #ifdef TIME_INFO
3703  timer.startTimer();
3704  #endif
3705  refinedTStep = refinedProblem->giveCurrentStep();
3706 
3707  size = refinedDomain->giveNumberOfDofManagers() * dofs;
3708  fineSolution.resize(size);
3709  coarseSolution.resize(size);
3710  errorSolution.resize(size);
3711 
3712  // map coarse solution
3713  uCoarse.resize( refinedProblem->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() ) );
3714  uCoarse.zero();
3715  mapper.mapAndUpdate(uCoarse, VM_Total, domain, refinedDomain, tStep);
3716 
3717  // get exact and coarse solution (including BC !!!)
3718  int pos = 1;
3719  for ( inode = 1; inode <= localNodeId; inode++ ) {
3720  node = refinedDomain->giveNode(inode);
3721  node->giveUnknownVector(nodeSolution, dofIdArray, VM_Total, refinedTStep);
3722  for ( idof = 1; idof <= dofs; idof++, pos++ ) {
3723  fineSolution.at(pos) = nodeSolution.at(idof);
3724  nodeDof = node->giveDofWithID(dofIdArray.at(idof));
3725  if ( nodeDof->hasBc(refinedTStep) == 0 ) {
3726  coarseSolution.at(pos) = uCoarse.at( nodeDof->__giveEquationNumber() );
3727  } else {
3728  // coarse solution is identical with fine solution at BC
3729  coarseSolution.at(pos) = fineSolution.at(pos);
3730  }
3731  }
3732  }
3733 
3734  errorSolution = fineSolution;
3735  errorSolution.subtract(coarseSolution);
3736 
3737  if ( this->normType == HuertaErrorEstimator :: L2Norm ) {
3738  FloatMatrix Nmatrix;
3739  FloatArray errorVectorGp, coarseVectorGp, fineVectorGp;
3740  /*
3741  * exactENorm = dotProduct(errorSolution.givePointer(), errorSolution.givePointer(), size);
3742  * coarseUNorm = dotProduct(coarseSolution.givePointer(), coarseSolution.givePointer(), size);
3743  * fineUNorm = dotProduct(fineSolution.givePointer(), fineSolution.givePointer(), size);
3744  */
3745  exactENorm = coarseUNorm = fineUNorm = mixedNorm = 0.0;
3746  for ( ielem = 1; ielem <= localElemId; ielem++ ) {
3747  element = refinedDomain->giveElement(ielem);
3748  interface = static_cast< HuertaErrorEstimatorInterface * >( element->giveInterface(HuertaErrorEstimatorInterfaceType) );
3749 
3750  for ( GaussPoint *gp: *element->giveDefaultIntegrationRulePtr() ) {
3751  double dV = element->computeVolumeAround(gp);
3752 
3753  interface->HuertaErrorEstimatorI_computeNmatrixAt(gp, Nmatrix);
3754 
3755  this->extractVectorFrom(element, errorSolution, errorVector, dofs, refinedTStep);
3756  errorVectorGp.beProductOf(Nmatrix, errorVector);
3757  exactENorm += errorVectorGp.computeSquaredNorm() * dV;
3758 
3759  this->extractVectorFrom(element, coarseSolution, coarseVector, dofs, refinedTStep);
3760  coarseVectorGp.beProductOf(Nmatrix, coarseVector);
3761  coarseUNorm += coarseVectorGp.computeSquaredNorm() * dV;
3762 
3763  mixedNorm += coarseVectorGp.dotProduct(errorVectorGp) * dV;
3764 
3765  this->extractVectorFrom(element, fineSolution, fineVector, dofs, refinedTStep);
3766  fineVectorGp.beProductOf(Nmatrix, fineVector);
3767  fineUNorm += fineVectorGp.computeSquaredNorm() * dV;
3768  }
3769  }
3770  } else if ( this->normType == HuertaErrorEstimator :: EnergyNorm ) {
3771  FloatArray tmpVector;
3772 
3773  #ifdef PRINT_ERROR
3774  double fineENorm, coarseENorm;
3775  int dim, pos, nelems;
3776 
3777  elemId = 1;
3778  element = domain->giveElement(elemId);
3779  dim = element->giveSpatialDimension();
3780  nelems = this->refineLevel + 1;
3781  while ( --dim ) {
3782  nelems *= this->refineLevel + 1;
3783  }
3784 
3785  nelems *= element->giveNumberOfNodes();
3786  coarseENorm = 0.0;
3787  pos = 0;
3788  #endif
3789 
3790  exactENorm = coarseUNorm = fineUNorm = mixedNorm = 0.0;
3791  for ( ielem = 1; ielem <= localElemId; ielem++ ) {
3792  if ( this->skipRegion( element->giveRegionNumber() ) != 0 ) {
3793  #ifdef PRINT_ERROR
3794  exactFineError.at(ielem) = 0.0;
3795  if ( ++pos == nelems ) {
3796  exactCoarseError.at(elemId) = coarseENorm;
3797  if ( ielem < localElemId ) {
3798  elemId++;
3799  element = domain->giveElement(elemId);
3800  dim = element->giveSpatialDimension();
3801  nelems = this->refineLevel + 1;
3802  while ( --dim ) {
3803  nelems *= this->refineLevel + 1;
3804  }
3805 
3806  nelems *= element->giveNumberOfNodes();
3807  coarseENorm = 0.0;
3808  pos = 0;
3809  }
3810  }
3811 
3812  #endif
3813 
3814  continue;
3815  }
3816 
3817  element = refinedDomain->giveElement(ielem);
3818  element->giveCharacteristicMatrix(mat, STIFFNESS_TYPE, refinedTStep);
3819 
3820  this->extractVectorFrom(element, errorSolution, errorVector, dofs, refinedTStep);
3821  tmpVector.beProductOf(mat, errorVector);
3822  exactENorm += tmpVector.dotProduct(errorVector); // Coarse and fine are identical? Also, unused.
3823  fineENorm = tmpVector.dotProduct(errorVector);
3824  coarseENorm += fineENorm;
3825 
3826  this->extractVectorFrom(element, coarseSolution, coarseVector, dofs, refinedTStep);
3827  tmpVector.beProductOf(mat, coarseVector);
3828  coarseUNorm += tmpVector.dotProduct(coarseVector);
3829 
3830  mixedNorm += tmpVector.dotProduct(errorVector);
3831 
3832  this->extractVectorFrom(element, fineSolution, fineVector, dofs, refinedTStep);
3833  tmpVector.beProductOf(mat, fineVector);
3834  fineUNorm += tmpVector.dotProduct(fineVector);
3835 
3836  #ifdef PRINT_ERROR
3837  exactFineError.at(ielem) = fineENorm;
3838  if ( ++pos == nelems ) {
3839  exactCoarseError.at(elemId) = coarseENorm;
3840  if ( ielem < localElemId ) {
3841  elemId++;
3842  element = domain->giveElement(elemId);
3843  dim = element->giveSpatialDimension();
3844  nelems = this->refineLevel + 1;
3845  while ( --dim ) {
3846  nelems *= this->refineLevel + 1;
3847  }
3848 
3849  nelems *= element->giveNumberOfNodes();
3850  coarseENorm = 0.0;
3851  pos = 0;
3852  }
3853  }
3854 
3855  #endif
3856  }
3857  } else {
3858  OOFEM_ERROR("Unsupported norm type");
3859  }
3860 
3861  #ifdef TIME_INFO
3862  timer.stopTimer();
3863  et_error = timer.getUtime();
3864 
3865  OOFEM_LOG_DEBUG("HEE info: whole 0: user time total %.2f s (setup %.2f s, init %.2f s, solve %.2f s, error %.2f s)\n",
3866  et_setup + et_init + et_solve + et_error,
3867  et_setup,
3868  et_init,
3869  et_solve,
3870  et_error);
3871  #endif
3872 
3873  delete refinedProblem;
3874 }
3875 #endif
3876 
3877 
3878 // pokud toto neni dostatecne obecne, da se funkce extractVectorFrom
3879 // do HuertaErrorEstimatorInterface a kazdy prvek si ji muze predefinovat
3880 
3881 void
3882 HuertaErrorEstimator :: extractVectorFrom(Element *element, FloatArray &vector, FloatArray &answer,
3883  int dofs, TimeStep *tStep)
3884 {
3885  int inode, idof, pos, p = 0;
3886 
3887  answer.resize(element->giveNumberOfDofManagers() * dofs);
3888  for ( inode = 1; inode <= element->giveNumberOfNodes(); inode++ ) {
3889  pos = ( element->giveDofManager(inode)->giveNumber() - 1 ) * dofs;
3890  for ( idof = 1; idof <= dofs; idof++ ) {
3891  answer.at(++p) = vector.at(pos + idof);
3892  }
3893  }
3894 }
3895 
3896 
3897 
3898 void
3899 HuertaErrorEstimator :: setupRefinedProblemProlog(const char *problemName, int problemId, IntArray &localNodeIdArray,
3900  int nodes, int elems, int csects, int mats, int loads, int funcs,
3901  IntArray &controlNode, IntArray &controlDof, TimeStep *tStep)
3902 {
3903  char line [ 1024 ];
3904  EngngModel *problem = this->domain->giveEngngModel();
3905  int i, nmstep, nsteps = 0;
3906  int ddfunc = 0, ddmSize = 0, ddvSize = 0, hpcSize = 0, hpcwSize = 0, renumber = 1;
3907  int controlMode = 0, hpcMode = 0, stiffMode = 0, maxIter = 30, reqIter = 3, manrmsteps = 0;
3908  double rtolv = -1.0, minStepLength = 0.0, initialStepLength = -1.0, stepLength = -1.0, psi = 1.0;
3909  IntArray ddm, hpc;
3910  FloatArray ddv, hpcw;
3911  IRResultType result; // Required by IR_GIVE_FIELD macro
3912 
3913 #if defined ( USE_OUTPUT_FILE ) || defined ( USE_CONTEXT_FILE )
3914  sprintf(line, "/home/dr/Huerta/%s_%d.out", problemName, problemId);
3915  skipUpdate = 0;
3916 #else
3917  sprintf(line, "/dev/null");
3918 #endif
3919  refinedReader.setOutputFileName(line);
3920 
3921  /* sprintf(skipUpdateString, "skipUpdate %d ", skipUpdate); */
3922 
3923 #ifdef USE_CONTEXT_FILE
3924  int contextOutputStep = 1;
3925 #endif
3926 
3927  sprintf(line, "Refined problem on %s %d", problemName, problemId);
3928  refinedReader.setDescription(line);
3929 
3930  if ( dynamic_cast< AdaptiveLinearStatic * >(problem) ) {
3931  DynamicInputRecord *ir = new DynamicInputRecord();
3932  ir->setRecordKeywordField(_IFT_LinearStatic_Name, 0);
3933  ir->setField(1, _IFT_EngngModel_nsteps);
3934  ir->setField(renumber, _IFT_EngngModel_renumberFlag);
3935 #ifdef USE_CONTEXT_FILE
3936  ir->setField(contextOutputStep, _IFT_EngngModel_contextoutputstep);
3937 #endif
3938  refinedReader.insertInputRecord(DataReader :: IR_emodelRec, ir);
3939  } else if ( dynamic_cast< AdaptiveNonLinearStatic * >(problem) ) {
3940  InputRecord *ir;
3941  nmstep = tStep->giveMetaStepNumber();
3942  ir = problem->giveMetaStep(nmstep)->giveAttributesRecord();
3943 
3944  IR_GIVE_OPTIONAL_FIELD(ir, stiffMode, _IFT_NonLinearStatic_stiffmode);
3945  IR_GIVE_OPTIONAL_FIELD(ir, controlMode, _IFT_NonLinearStatic_controlmode);
3946 
3947  switch ( controlMode ) {
3948  case 0:
3949  IR_GIVE_OPTIONAL_FIELD(ir, maxIter, _IFT_CylindricalALM_maxiter);
3950  IR_GIVE_OPTIONAL_FIELD(ir, reqIter, _IFT_CylindricalALM_reqiterations);
3951  IR_GIVE_OPTIONAL_FIELD(ir, minStepLength, _IFT_CylindricalALM_minsteplength);
3952  IR_GIVE_OPTIONAL_FIELD(ir, manrmsteps, _IFT_CylindricalALM_manrmsteps);
3953  IR_GIVE_FIELD(ir, stepLength, _IFT_CylindricalALM_steplength);
3954 
3955  initialStepLength = stepLength;
3956  IR_GIVE_OPTIONAL_FIELD(ir, initialStepLength, _IFT_CylindricalALM_initialsteplength);
3957 
3958  IR_GIVE_OPTIONAL_FIELD(ir, psi, _IFT_CylindricalALM_psi);
3959  IR_GIVE_OPTIONAL_FIELD(ir, hpcMode, _IFT_CylindricalALM_hpcmode);
3960 
3961  IR_GIVE_OPTIONAL_FIELD(ir, hpc, _IFT_CylindricalALM_hpc);
3962  IR_GIVE_OPTIONAL_FIELD(ir, hpcw, _IFT_CylindricalALM_hpcw);
3963  IR_GIVE_FIELD(ir, rtolv, _IFT_CylindricalALM_rtolv);
3964 
3965  hpcSize = hpc.giveSize();
3966  hpcwSize = hpcw.giveSize();
3967  break;
3968  case 1:
3969  IR_GIVE_OPTIONAL_FIELD(ir, maxIter, _IFT_NRSolver_maxiter);
3970  IR_GIVE_OPTIONAL_FIELD(ir, minStepLength, _IFT_NRSolver_minsteplength);
3971  IR_GIVE_OPTIONAL_FIELD(ir, manrmsteps, _IFT_NRSolver_manrmsteps);
3972 
3973  IR_GIVE_OPTIONAL_FIELD(ir, ddm, _IFT_NRSolver_ddm);
3974  IR_GIVE_OPTIONAL_FIELD(ir, ddv, _IFT_NRSolver_ddv);
3975  IR_GIVE_OPTIONAL_FIELD(ir, ddfunc, _IFT_NRSolver_ddfunc);
3976  IR_GIVE_FIELD(ir, rtolv, _IFT_NRSolver_rtolv);
3977 
3978  ddmSize = ddm.giveSize();
3979  ddvSize = ddv.giveSize();
3980  break;
3981  default:
3982  OOFEM_ERROR("Unsupported control mode");
3983  }
3984 
3985  if ( problemId != 0 ) {
3986  DynamicInputRecord *ir = new DynamicInputRecord();
3987 
3988  int bcSize = controlNode.giveSize(), ddActiveSize = 0;
3989 
3990  if ( controlMode == 1 ) {
3991  // count original dd active in refined problem
3992  for ( i = 1; i <= ddmSize; i += 2 ) {
3993  if ( localNodeIdArray.at( ddm.at(i) ) == 0 ) {
3994  continue;
3995  }
3996 
3997  ddActiveSize++;
3998  }
3999  }
4000 
4001  // note: it is impossible to check the solution using the external file;
4002  // first of all adaptnlinearstatic must be changed to nonlinearstatic to prevent adaptivity;
4003  // however the most severe problem is that in ddv only zeros are specified !!!
4004  // it would be desirable to print there coarse mesh displacement
4005 
4006  if ( nmstep == 1 ) {
4007  // do not specify context because that would result in recursive call to HEE
4008  // because adaptnlinearstatic analysis type is used
4009 
4010  ir->setRecordKeywordField(_IFT_AdaptiveNonLinearStatic_Name, 0);
4011  ir->setField(tStep->giveNumber(), _IFT_EngngModel_nsteps);
4012  ir->setField(renumber, _IFT_EngngModel_renumberFlag);
4013  ir->setField(rtolv, "rtolv");
4014  ir->setField(maxIter, "maxiter");
4015  ir->setField(reqIter, "reqiterations");
4016  ir->setField(minStepLength, "minsteplength");
4017  ir->setField(stiffMode, "stiffmode");
4018  ir->setField(manrmsteps, "manrmsteps");
4019  ir->setField(manrmsteps, "manrmsteps");
4020  //ir->setField(skipUpdate, "skipupdate");
4021 #ifdef USE_CONTEXT_FILE
4022  ir->setField(contextOutputStep, _IFT_EngngModel_contextoutputstep);
4023 #endif
4024 
4025  // this is not relevant but it is required
4026  // the refined problem is made adaptive only to enable call to initializeAdaptiveFrom
4027  ir->setField(3, "eetype");
4028  ir->setField(0, "meshpackage");
4029  ir->setField(0.1, "requiredError");
4030  ir->setField(0.01, "minelemsize");
4031  } else {
4032  // if there are more than just a single metastep, it is necessary to produce input with at least the
4033  // same number of metasteps as is the number of active metastep,
4034  // because initializeAdaptiveFrom uses a copy constructor for timestep;
4035  // only the active metastep should be filled by real values
4036 
4037  // do not specify context because that would result in recursive call to HEE
4038  // because adaptnlinearstatic analysis type is used
4039 
4040  ir->setRecordKeywordField(_IFT_AdaptiveNonLinearStatic_Name, 0);
4041  ir->setField(1, _IFT_EngngModel_nsteps);
4042  ir->setField(nmstep, _IFT_EngngModel_nmsteps);
4043  ir->setField(renumber, _IFT_EngngModel_renumberFlag);
4044  ir->setField(1, "equilmc");
4045  ir->setField(1, "controlmode");
4046  //ir->setField(skipUpdate, "skipupdate");
4047 #ifdef USE_CONTEXT_FILE
4048  ir->setField(contextOutputStep, _IFT_EngngModel_contextoutputstep);
4049 #endif
4050 
4051  // this is not relevant but it is required
4052  // the refined problem is made adaptive only to enable call to initializeAdaptiveFrom
4053  ir->setField(3, "eetype");
4054  ir->setField(0, "meshpackage");
4055  ir->setField(0.1, "requiredError");
4056  ir->setField(0.01, "minelemsize");
4057 
4058  refinedReader.insertInputRecord(DataReader :: IR_emodelRec, ir);
4059 
4060  // IMPORTANT: the total number of steps should be equal to the current step number
4061  // (to enable skipUpdate)
4062  // therefore the number of steps in the current metastep is modified
4063 
4064  // create fictitious metasteps
4065  for ( i = 1; i < nmstep; i++ ) {
4066  DynamicInputRecord *ir_meta = new DynamicInputRecord();
4067  ir_meta->setRecordKeywordField(_IFT_MetaStep_Name, i);
4068  ir_meta->setField(problem->giveMetaStep(i)->giveNumberOfSteps(), _IFT_MetaStep_nsteps);
4069  ir_meta->setField(rtolv, "rtolv");
4070  refinedReader.insertInputRecord(DataReader :: IR_mstepRec, ir_meta);
4071  nsteps += problem->giveMetaStep(i)->giveNumberOfSteps();
4072  }
4073 
4074  // create active metastep
4075  ir = new DynamicInputRecord();
4076  ir->setRecordKeywordField(_IFT_MetaStep_Name, nmstep);
4077  ir->setField(tStep->giveNumber() - nsteps, _IFT_MetaStep_nsteps);
4078  ir->setField(rtolv, "rtolv");
4079  ir->setField(maxIter, "maxiter");
4080  ir->setField(reqIter, "reqiterations");
4081  ir->setField(minStepLength, "minsteplength");
4082  ir->setField(stiffMode, "stiffmode");
4083  ir->setField(manrmsteps, "manrmsteps");
4084  ir->setField(1, "equilmc");
4085  ir->setField(1, "controlmode");
4086  }
4087 
4088  // apply refine problem bc as dd
4089  IntArray ddm_vals( 2 * ( bcSize + ddActiveSize ) );
4090 
4091  for ( i = 1; i <= bcSize; i++ ) {
4092  ddm_vals.at(i * 2 - 1) = controlNode.at(i);
4093  ddm_vals.at(i * 2) = controlDof.at(i);
4094  }
4095 
4096  if ( controlMode == 1 ) {
4097  // apply active original dd
4098  for ( i = 1; i <= ddmSize; i += 2 ) {
4099  if ( localNodeIdArray.at( ddm.at(i) ) == 0 ) {
4100  continue;
4101  }
4102 
4103  ddm_vals.at(bcSize * 2 + i * 2 - 1) = -localNodeIdArray.at( ddm.at(i) );
4104  ddm_vals.at(bcSize * 2 + i * 2) = ddm.at(i + 1);
4105  }
4106  }
4107  ir->setField(ddm_vals, "ddm");
4108 
4109  FloatArray ddv_vals(bcSize + ddActiveSize);
4110 
4111  // apply refined problem bc values
4112  for ( i = 1; i <= bcSize; i++ ) {
4113  ddv_vals.at(i) = 0.; // no further increment required, just to recover equilibrium of mapped state
4114  }
4115 
4116  if ( controlMode == 1 ) {
4117  // apply active original dd values
4118  for ( i = 1; i <= ddvSize; i++ ) {
4119  if ( localNodeIdArray.at( ddm.at(2 * i - 1) ) == 0 ) {
4120  continue;
4121  }
4122 
4123  ddv_vals.at(bcSize + i) = ddv.at(i);
4124  }
4125  }
4126  ir->setField(ddv_vals, "ddv");
4127 
4128  // use last funcs to control dd
4129  ir->setField(funcs, _IFT_NRSolver_ddfunc);
4130 
4131  refinedReader.insertInputRecord(DataReader :: IR_mstepRec, ir);
4132  } else {
4133  // it makes not too much sense to solve exact problem from beginning if adaptive restart is used
4134  // because the mesh may be already derefined in regions of no interest (but anyway ...)
4135  // however if adaptive restart is applied, number of current step does not correspond to the time
4136  // (step number = time + 1) because step number was increased when recovering equilibrium at the last time;
4137  // therefore problem -> giveCurrentStep() -> giveNumber() is not used and the number of steps is
4138  // recovered from the current time
4139 
4140  // do not prescribe skipUpdate here !!!
4141 
4142  // HUHU dodelat metasteps, dodelat paralelni zpracovani
4143  DynamicInputRecord *ir = new DynamicInputRecord();
4144  ir->setRecordKeywordField("nonlinearstatic", 0);
4145  ir->setField( ( int ) ( problem->giveCurrentStep()->giveTargetTime() + 1.5 ), "nsteps" );
4146  ir->setField(renumber, "renumber");
4147  ir->setField(rtolv, "rtolv");
4148  ir->setField(maxIter, "maxiter");
4149  ir->setField(reqIter, "reqiterations");
4150  ir->setField(minStepLength, "minsteplength");
4151  ir->setField(stiffMode, "stiffmode");
4152  ir->setField(manrmsteps, "manrmsteps");
4153  ir->setField(1, "equilmc");
4154  ir->setField(controlMode, "controlmode");
4155 #ifdef USE_CONTEXT_FILE
4156  ir->setField(contextOutputStep, _IFT_EngngModel_contextoutputstep);
4157 #endif
4158 
4159  switch ( controlMode ) {
4160  case 0:
4161  ir->setField(stepLength, "stepLength");
4162  ir->setField(initialStepLength, "initialStepLength");
4163  ir->setField(psi, "psi");
4164  ir->setField(hpcMode, "hpcmode");
4165 
4166  if ( hpcSize != 0 ) {
4167  // apply all original hpc
4168  IntArray hpc_vals(hpcSize);
4169  for ( i = 1; i <= hpcSize; i += 2 ) {
4170  hpc_vals.at(i * 2 - 1) = -localNodeIdArray.at( hpc.at(i) );
4171  hpc_vals.at(i * 2) = hpc.at(i + 1);
4172  }
4173  ir->setField(hpc_vals, "hpc");
4174  }
4175 
4176  if ( hpcwSize != 0 ) {
4177  ir->setField(hpcw, "hpcw");
4178  }
4179 
4180  break;
4181  case 1:
4182  if ( ddmSize != 0 ) {
4183  // apply all original dd
4184  IntArray ddm_vals(ddmSize);
4185  for ( i = 1; i <= ddmSize; i += 2 ) {
4186  ddm_vals.at(i * 2 - 1) = -localNodeIdArray.at( ddm.at(i) );
4187  ddm_vals.at(i * 2) = ddm.at(i + 1);
4188  }
4189  ir->setField("ddm");
4190  }
4191 
4192  if ( ddvSize != 0 ) {
4193  ir->setField(ddv, "ddv");
4194  }
4195 
4196  // use the original funcs to control dd
4197  ir->setField(ddfunc, _IFT_NRSolver_ddfunc);
4198  break;
4199  }
4200 
4201  refinedReader.insertInputRecord(DataReader :: IR_emodelRec, ir);
4202  }
4203  } else {
4204  OOFEM_ERROR("Unsupported analysis type");
4205  }
4206 
4207  DynamicInputRecord *ir = new DynamicInputRecord();
4208  ir->setField(this->domain->giveNumberOfSpatialDimensions(), _IFT_Domain_numberOfSpatialDimensions);
4209  if ( this->domain->isAxisymmetric() ) {
4210  ir->setField(_IFT_Domain_axisymmetric);
4211  }
4212 
4213 
4214  refinedReader.insertInputRecord(DataReader :: IR_domainCompRec, ir);
4215 
4216  ir = new DynamicInputRecord();
4217  ir->setRecordKeywordField(_IFT_OutputManager_Name, 0);
4218 #ifdef USE_OUTPUT_FILE
4219  ir->setField(_IFT_OutputManager_tstepall);
4220  ir->setField(_IFT_OutputManager_dofmanall);
4221  ir->setField(_IFT_OutputManager_elementall);
4222 #endif
4223 
4224  refinedReader.insertInputRecord(DataReader :: IR_outManRec, ir);
4225 
4226  ir = new DynamicInputRecord();
4227  ir->setRecordKeywordField("", 0);
4228  ir->setField(nodes, _IFT_Domain_ndofman);
4229  ir->setField(elems, _IFT_Domain_nelem);
4230  ir->setField(mats, _IFT_Domain_ncrosssect);
4231  ir->setField(csects, _IFT_Domain_nmat);
4232  ir->setField(loads, _IFT_Domain_nbc);
4233  ir->setField(funcs, _IFT_Domain_nfunct);
4234  refinedReader.insertInputRecord(DataReader :: IR_domainCompRec, ir);
4235 }
4236 
4237 
4238 
4239 void
4240 HuertaErrorEstimator :: setupRefinedProblemEpilog1(int csects, int mats, int loads, int nlbarriers)
4241 {
4242  Domain *domain = this->domain;
4243 
4244  /* copy csects, mats, loads */
4245 
4246  for ( int i = 1; i <= csects; i++ ) {
4247  DynamicInputRecord *ir = new DynamicInputRecord();
4248  domain->giveCrossSection(i)->giveInputRecord(* ir);
4249  refinedReader.insertInputRecord(DataReader :: IR_crosssectRec, ir);
4250  }
4251 
4252  for ( int i = 1; i <= mats; i++ ) {
4253  DynamicInputRecord *ir = new DynamicInputRecord();
4254  domain->giveMaterial(i)->giveInputRecord(* ir);
4255  refinedReader.insertInputRecord(DataReader :: IR_matRec, ir);
4256  }
4257 
4258  for ( int i = 1; i <= loads; i++ ) {
4259  DynamicInputRecord *ir = new DynamicInputRecord();
4260  domain->giveLoad(i)->giveInputRecord(* ir);
4261  refinedReader.insertInputRecord(DataReader :: IR_bcRec, ir);
4262  }
4263 }
4264 
4265 
4266 
4267 void
4268 HuertaErrorEstimator :: setupRefinedProblemEpilog2(int funcs)
4269 {
4270  Domain *domain = this->domain;
4271 
4272  /* copy tfuncs */
4273 
4274  for ( int i = 1; i <= funcs; i++ ) {
4275  DynamicInputRecord *ir = new DynamicInputRecord();
4276  domain->giveFunction(i)->giveInputRecord(* ir);
4277  refinedReader.insertInputRecord(DataReader :: IR_funcRec, ir);
4278  }
4279 
4280  if ( this->mode == HEE_nlinear ) {
4281  DynamicInputRecord *ir = new DynamicInputRecord();
4282  ir->setRecordKeywordField(_IFT_HeavisideTimeFunction_Name, funcs + 1);
4283  ir->setField(this->domain->giveEngngModel()->giveCurrentStep()->giveTargetTime() - 0.1, _IFT_HeavisideTimeFunction_origin);
4284  ir->setField(1., _IFT_HeavisideTimeFunction_value);
4285  refinedReader.insertInputRecord(DataReader :: IR_funcRec, ir);
4286  }
4287 }
4288 } // end namespace oofem
oofem::Element::giveCrossSection
CrossSection * giveCrossSection()
Definition: element.C:495
_IFT_HuertaErrorEstimator_normtype
#define _IFT_HuertaErrorEstimator_normtype
Definition: huertaerrorestimator.h:51
oofem::HuertaErrorEstimator::setupRefinedProblemEpilog1
void setupRefinedProblemEpilog1(int csects, int mats, int loads, int nlbarriers)
Definition: huertaerrorestimator.C:4240
_IFT_HuertaErrorEstimator_initialskipsteps
#define _IFT_HuertaErrorEstimator_initialskipsteps
Definition: huertaerrorestimator.h:55
_IFT_Domain_nfunct
#define _IFT_Domain_nfunct
Definition: domain.h:63
oofem::DynamicInputRecord::setField
void setField(int item, InputFieldType id)
Definition: dynamicinputrecord.C:299
oofem::RemeshingCriteria
The base class for all remeshing criteria.
Definition: remeshingcrit.h:61
oofem::GaussPoint::giveMaterialMode
MaterialMode giveMaterialMode()
Returns corresponding material mode of receiver.
Definition: gausspoint.h:191
oofem::EModelDefaultEquationNumbering
The representation of EngngModel default unknown numbering.
Definition: unknownnumberingscheme.h:80
oofem::HuertaRemeshingCriteria
The class representing Huerta remeshing criteria.
Definition: huertaerrorestimator.h:285
oofem::FloatArray::subtract
void subtract(const FloatArray &src)
Subtracts array src to receiver.
Definition: floatarray.C:258
oofem::DofManager_null
DofManager in active domain is shared only by remote elements (these are only introduced for nonlocal...
Definition: dofmanager.h:88
oofem::EngngModel::giveContextFileName
std::string giveContextFileName(int tStepNumber, int stepVersion) const
Returns the filename for the context file for the given step and version.
Definition: engngm.C:1691
_IFT_CylindricalALM_minsteplength
#define _IFT_CylindricalALM_minsteplength
Definition: calmls.h:54
oofem::Element::giveRegionNumber
int giveRegionNumber()
Definition: element.C:507
oofem::Element::giveIPValue
virtual int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
Returns the integration point corresponding value in full form.
Definition: element.C:1257
oofem::HuertaErrorEstimator::extractVectorFrom
void extractVectorFrom(Element *element, FloatArray &vector, FloatArray &answer, int dofs, TimeStep *tStep)
Extracts nodal vector from global vector for each dof of all element nodes.
Definition: huertaerrorestimator.C:3882
_IFT_LinearStatic_Name
#define _IFT_LinearStatic_Name
Definition: linearstatic.h:43
oofem::MetaStep::giveAttributesRecord
InputRecord * giveAttributesRecord()
Returns e-model attributes.
Definition: metastep.h:95
oofem::HuertaErrorEstimator::refinedMesh
RefinedMesh refinedMesh
Mesh refinement.
Definition: huertaerrorestimator.h:125
oofem::MetaStep::giveNumberOfSteps
int giveNumberOfSteps()
Returns number of Steps it represent.
Definition: metastep.h:91
oofem::FileDataStream
Implementation of FileDataStream representing DataStream interface to file i/o.
Definition: datastream.h:136
oofem::Domain
Class and object Domain.
Definition: domain.h:115
_IFT_OutputManager_tstepall
#define _IFT_OutputManager_tstepall
Definition: outputmanager.h:48
oofem::Element::computeVectorOf
void computeVectorOf(ValueModeType u, TimeStep *tStep, FloatArray &answer)
Returns local vector of unknowns.
Definition: element.C:86
oofem::Element::giveDefaultIntegrationRulePtr
virtual IntegrationRule * giveDefaultIntegrationRulePtr()
Access method for default integration rule.
Definition: element.h:822
oofem::EngngModel::giveNumberOfDomainEquations
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
_IFT_HuertaRemeshingCriteria_noremesh
#define _IFT_HuertaRemeshingCriteria_noremesh
Definition: huertaerrorestimator.h:67
oofem::DynamicDataReader
Class representing the implementation of a dynamic data reader for in-code use.
Definition: dynamicdatareader.h:53
_IFT_AdaptiveNonLinearStatic_Name
#define _IFT_AdaptiveNonLinearStatic_Name
Definition: adaptnlinearstatic.h:45
oofem::refinedReader
static DynamicDataReader refinedReader("huerta")
oofem::FEMComponent::domain
Domain * domain
Link to domain object, useful for communicating with other FEM components.
Definition: femcmpnn.h:82
_IFT_HuertaRemeshingCriteria_refinecoeff
#define _IFT_HuertaRemeshingCriteria_refinecoeff
Definition: huertaerrorestimator.h:69
oofem::Element::computeMeanSize
double computeMeanSize()
Computes the size of the element defined as its length.
Definition: element.C:1078
oofem::Domain::giveNumberOfDofManagers
int giveNumberOfDofManagers() const
Returns number of dof managers in domain.
Definition: domain.h:432
_IFT_NRSolver_manrmsteps
#define _IFT_NRSolver_manrmsteps
Definition: nrsolver.h:57
oofem::HuertaErrorEstimator::giveElementError
virtual double giveElementError(EE_ErrorType type, Element *elem, TimeStep *tStep)
Returns the element error.
Definition: huertaerrorestimator.C:525
oofem::HuertaErrorEstimator::solveRefinedWholeProblem
void solveRefinedWholeProblem(IntArray &localNodeIdArray, IntArray &globalNodeIdArray, TimeStep *tStep)
Solves the refined whole problem.
Definition: huertaerrorestimator.C:3517
_IFT_CylindricalALM_reqiterations
#define _IFT_CylindricalALM_reqiterations
Definition: calmls.h:58
oofem::HuertaRemeshingCriteria::giveRequiredDofManDensity
virtual double giveRequiredDofManDensity(int num, TimeStep *tStep, int relative=0)
Returns the required mesh size n given dof manager.
Definition: huertaerrorestimator.C:702
oofem::EIPrimaryUnknownMapper
The class implementing the primary unknown mapper using element interpolation functions.
Definition: eleminterpunknownmapper.h:50
oofem::REGISTER_ErrorEstimator
REGISTER_ErrorEstimator(CombinedZZSIErrorEstimator, EET_CZZSI)
oofem::HuertaErrorEstimator::estimateError
virtual int estimateError(EE_ErrorMode err_mode, TimeStep *tStep)
Estimates the error on associated domain at given time step.
Definition: huertaerrorestimator.C:136
oofem::Domain::giveNumberOfBoundaryConditions
int giveNumberOfBoundaryConditions() const
Returns number of boundary conditions in domain.
Definition: domain.h:440
oofem::Node::giveLocalCoordinateTriplet
FloatMatrix * giveLocalCoordinateTriplet()
Returns pointer to local coordinate triplet in node.
Definition: node.h:158
oofem::ErrorEstimator::rc
std::unique_ptr< RemeshingCriteria > rc
Definition: errorestimator.h:84
oofem::HuertaRemeshingCriteria::mode
HuertaRemeshingCriteriaModeType mode
Mode of receiver.
Definition: huertaerrorestimator.h:299
_IFT_HuertaErrorEstimator_refinelevel
#define _IFT_HuertaErrorEstimator_refinelevel
Definition: huertaerrorestimator.h:52
_IFT_CylindricalALM_initialsteplength
#define _IFT_CylindricalALM_initialsteplength
Definition: calmls.h:56
oofem::DataStream
The purpose of DataStream abstract class is to allow to store/restore context to different streams...
Definition: datastream.h:54
oofem::FloatArray::storeYourself
contextIOResultType storeYourself(DataStream &stream) const
Definition: floatarray.C:872
oofem::IntArray::zero
void zero()
Sets all component to zero.
Definition: intarray.C:52
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::max
int max(int i, int j)
Returns bigger value form two given decimals.
Definition: mathfem.h:71
oofem::HuertaErrorEstimator::globalWENorm
double globalWENorm
Global weighted error norm.
Definition: huertaerrorestimator.h:107
oofem::EngngModel::giveMetaStep
MetaStep * giveMetaStep(int i)
Returns the i-th meta step.
Definition: engngm.C:1765
CM_State
#define CM_State
Definition: contextmode.h:46
_IFT_NRSolver_ddm
#define _IFT_NRSolver_ddm
Definition: nrsolver.h:59
oofem::FEInterpolation::giveInterpolationOrder
int giveInterpolationOrder()
Returns the interpolation order.
Definition: feinterpol.h:159
_IFT_CylindricalALM_psi
#define _IFT_CylindricalALM_psi
Definition: calmls.h:51
oofem::HuertaErrorEstimator::globalErrorEstimate
double globalErrorEstimate
Global error estimate (relative)
Definition: huertaerrorestimator.h:111
_IFT_HuertaErrorEstimator_exact
#define _IFT_HuertaErrorEstimator_exact
Definition: huertaerrorestimator.h:60
_IFT_CylindricalALM_hpc
#define _IFT_CylindricalALM_hpc
Definition: calmls.h:62
oofem::Domain::giveConnectivityTable
ConnectivityTable * giveConnectivityTable()
Returns receiver&#39;s associated connectivity table.
Definition: domain.C:1170
oofem::Domain::isAxisymmetric
bool isAxisymmetric()
Returns true of axisymmetry is in effect.
Definition: domain.C:1074
oofem::HuertaErrorEstimatorInterface::setupRefinedElementProblem3D
void setupRefinedElementProblem3D(Element *element, RefinedElement *refinedElement, int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, HuertaErrorEstimatorInterface::SetupMode mode, TimeStep *tStep, int nodes, FloatArray **corner, FloatArray *midSide, FloatArray *midFace, FloatArray &midNode, int &localNodeId, int &localElemId, int &localBcId, int hexaSideNode[1][3], int hexaFaceNode[1][3], IntArray &controlNode, IntArray &controlDof, HuertaErrorEstimator::AnalysisMode aMode, const char *hexatype)
Definition: huertaerrorestimator.C:2045
oofem::HuertaRemeshingCriteria::stateCounter
StateCounterType stateCounter
Actual values (densities) state counter.
Definition: huertaerrorestimator.h:297
oofem::RefinedElement::giveFineNodeArray
IntArray * giveFineNodeArray(int node)
Definition: refinedelement.C:81
MAX_COARSE_RATE
#define MAX_COARSE_RATE
Definition: huertaerrorestimator.C:725
oofem::Domain::giveEngngModel
EngngModel * giveEngngModel()
Returns engineering model to which receiver is associated.
Definition: domain.C:433
_IFT_CylindricalALM_hpcmode
#define _IFT_CylindricalALM_hpcmode
Definition: calmls.h:61
_IFT_NRSolver_rtolv
#define _IFT_NRSolver_rtolv
Definition: nrsolver.h:63
_IFT_DofManager_dofidmask
#define _IFT_DofManager_dofidmask
Definition: dofmanager.h:52
oofem::TimeStep::giveTargetTime
double giveTargetTime()
Returns target time.
Definition: timestep.h:146
oofem::TimeStep::incrementStateCounter
void incrementStateCounter()
Updates solution state counter.
Definition: timestep.h:190
oofem::Element
Abstract base class for all finite elements.
Definition: element.h:145
oofem::FEMComponent::saveContext
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: femcmpnn.C:51
oofem::RefinedElement::giveBoundaryLoadArray3D
bool giveBoundaryLoadArray3D(int inode, Element *element, std::vector< IntArray > &boundaryLoadList)
Definition: refinedelement.C:469
oofem::EE_ErrorType
EE_ErrorType
Type characterizing different type of element errors.
Definition: errorestimator.h:67
oofem::CIO_IOERR
General IO error.
Definition: contextioresulttype.h:43
oofem::EngngModel::solveYourself
virtual void solveYourself()
Starts solution process.
Definition: engngm.C:501
_IFT_HeavisideTimeFunction_origin
#define _IFT_HeavisideTimeFunction_origin
Definition: heavisidetimefunction.h:44
oofem::Dof::__giveEquationNumber
virtual int __giveEquationNumber() const =0
Returns equation number of receiver, usually assigned by emodel.
oofem::Domain::giveNumberOfElements
int giveNumberOfElements() const
Returns number of elements in domain.
Definition: domain.h:434
oofem::HuertaRemeshingCriteria::giveRemeshingStrategy
virtual RemeshingStrategy giveRemeshingStrategy(TimeStep *tStep)
Determines, if the remeshing is needed, and if needed, the type of strategy used. ...
Definition: huertaerrorestimator.C:718
oofem::Material::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: material.C:110
oofem::CrossSection::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: crosssection.C:82
oofem::Element::giveCharacteristicMatrix
virtual void giveCharacteristicMatrix(FloatMatrix &answer, CharType type, TimeStep *tStep)
Computes characteristic matrix of receiver of requested type in given time step.
Definition: element.C:569
oofem::DynamicDataReader::writeToFile
void writeToFile(const char *fileName)
Writes all input records to file.
Definition: dynamicdatareader.C:90
oofem::Domain::giveNumberOfSpatialDimensions
int giveNumberOfSpatialDimensions()
Returns number of spatial dimensions.
Definition: domain.C:1067
oofem::MaterialMode
MaterialMode
Type representing material mode of integration point.
Definition: materialmode.h:89
_IFT_NRSolver_maxiter
#define _IFT_NRSolver_maxiter
Definition: nrsolver.h:54
oofem::masterRun
static bool masterRun
Definition: huertaerrorestimator.C:109
OOFEM_LOG_DEBUG
#define OOFEM_LOG_DEBUG(...)
Definition: logger.h:128
oofem::FEMComponent::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: femcmpnn.C:77
oofem::IntArray
Class implementing an array of integers.
Definition: intarray.h:61
oofem::IntArray::at
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
oofem::DataStream::read
virtual int read(int *data, int count)=0
Reads count integer values into array pointed by data.
oofem::HuertaErrorEstimator::solveRefinedPatchProblem
void solveRefinedPatchProblem(int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, TimeStep *tStep)
Solves the refined patch problem.
Definition: huertaerrorestimator.C:3190
_IFT_OutputManager_dofmanall
#define _IFT_OutputManager_dofmanall
Definition: outputmanager.h:50
oofem::HuertaErrorEstimator::requiredError
double requiredError
Required error to obtain.
Definition: huertaerrorestimator.h:129
_IFT_OutputManager_elementall
#define _IFT_OutputManager_elementall
Definition: outputmanager.h:51
oofem::Element::giveNumberOfDofManagers
virtual int giveNumberOfDofManagers() const
Definition: element.h:656
THROW_CIOERR
#define THROW_CIOERR(e)
Definition: contextioerr.h:61
oofem::Element::giveInterpolation
virtual FEInterpolation * giveInterpolation() const
Definition: element.h:629
oofem::dofManagerParallelMode
dofManagerParallelMode
In parallel mode, this type indicates the mode of DofManager.
Definition: dofmanager.h:80
oofem::HuertaErrorEstimator::stateCounter
StateCounterType stateCounter
Actual state counter.
Definition: huertaerrorestimator.h:117
oofem::IntegrationRule
Abstract base class representing integration rule.
Definition: integrationrule.h:94
oofem::HuertaErrorEstimator::normType
NormType normType
Type of norm used.
Definition: huertaerrorestimator.h:115
_IFT_BoundaryCondition_Name
#define _IFT_BoundaryCondition_Name
Definition: boundarycondition.h:48
OOFEM_LOG_RELEVANT
#define OOFEM_LOG_RELEVANT(...)
Definition: logger.h:126
_IFT_DofManager_bc
#define _IFT_DofManager_bc
Definition: dofmanager.h:54
oofem::HuertaErrorEstimatorInterface::HuertaErrorEstimatorI_computeNmatrixAt
virtual void HuertaErrorEstimatorI_computeNmatrixAt(GaussPoint *gp, FloatMatrix &answer)=0
oofem::HuertaRemeshingCriteria::estimateMeshDensities
virtual int estimateMeshDensities(TimeStep *tStep)
Estimates the nodal densities.
Definition: huertaerrorestimator.C:729
oofem::ErrorEstimator::giveErrorEstimatorType
ErrorEstimatorType giveErrorEstimatorType() const
Returns error estimation type of receiver.
Definition: errorestimator.h:142
_IFT_GeneralBoundaryCondition_timeFunct
#define _IFT_GeneralBoundaryCondition_timeFunct
Definition: generalboundarycondition.h:46
oofem::HuertaErrorEstimator::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: huertaerrorestimator.C:566
oofem::ErrorEstimator::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: errorestimator.C:60
oofem::Element::giveNumberOfNodes
virtual int giveNumberOfNodes() const
Returns number of nodes of receiver.
Definition: element.h:662
oofem::RemeshingCriteria::ee
ErrorEstimator * ee
Definition: remeshingcrit.h:64
_IFT_MetaStep_nsteps
#define _IFT_MetaStep_nsteps
Definition: metastep.h:44
_IFT_NonLinearStatic_stiffmode
#define _IFT_NonLinearStatic_stiffmode
Definition: nlinearstatic.h:46
oofem::DataStream::write
virtual int write(const int *data, int count)=0
Writes count integer values from array pointed by data.
oofem::HuertaErrorEstimator::restoreContext
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: huertaerrorestimator.C:669
oofem::HuertaRemeshingCriteria::wError
bool wError
Weighted error flag.
Definition: huertaerrorestimator.h:309
oofem::InstanciateProblem
EngngModel * InstanciateProblem(DataReader &dr, problemMode mode, int contextFlag, EngngModel *_master, bool parallelFlag)
Instanciates the new problem.
Definition: util.C:45
oofem::Domain::giveCrossSection
CrossSection * giveCrossSection(int n)
Service for accessing particular domain cross section model.
Definition: domain.C:339
_IFT_CylindricalALM_maxiter
#define _IFT_CylindricalALM_maxiter
Definition: calmls.h:52
oofem::Element::computeVolumeAround
virtual double computeVolumeAround(GaussPoint *gp)
Returns volume related to given integration point.
Definition: element.h:518
oofem::DofManager::giveUnknownVector
void giveUnknownVector(FloatArray &answer, const IntArray &dofMask, ValueModeType mode, TimeStep *tStep, bool padding=false)
Assembles the vector of unknowns in global c.s for given dofs of receiver.
Definition: dofmanager.C:685
oofem::Dof::giveBcId
virtual int giveBcId()=0
Returns the id of associated boundary condition, if there is any.
_IFT_Domain_nmat
#define _IFT_Domain_nmat
Definition: domain.h:59
_IFT_Domain_nelem
#define _IFT_Domain_nelem
Definition: domain.h:58
OOFEM_LOG_INFO
#define OOFEM_LOG_INFO(...)
Definition: logger.h:127
oofem::TimeStep::giveNumber
int giveNumber()
Returns receiver&#39;s number.
Definition: timestep.h:129
oofem::DofManager::giveNumberOfDofs
int giveNumberOfDofs() const
Definition: dofmanager.C:279
oofem::Domain::giveMaterial
Material * giveMaterial(int n)
Service for accessing particular domain material model.
Definition: domain.C:281
oofem::Domain::giveElement
Element * giveElement(int n)
Service for accessing particular domain fe element.
Definition: domain.C:160
oofem::FloatArray::dotProduct
double dotProduct(const FloatArray &x) const
Computes the dot product (or inner product) of receiver and argument.
Definition: floatarray.C:463
oofem::HuertaErrorEstimatorInterface
The element interface corresponding to HuertaErrorEstimator.
Definition: huertaerrorestimator.h:230
oofem::TimeStep::giveSolutionStateCounter
StateCounterType giveSolutionStateCounter()
Returns current solution state counter.
Definition: timestep.h:188
oofem::HuertaErrorEstimator::wError
bool wError
Weighted error flag.
Definition: huertaerrorestimator.h:131
_IFT_CylindricalALM_rtolv
#define _IFT_CylindricalALM_rtolv
Definition: calmls.h:71
oofem::ErrorEstimator::skippedNelems
int skippedNelems
Number of skipped elements.
Definition: errorestimator.h:93
oofem::HuertaErrorEstimator::giveRemeshingCrit
virtual RemeshingCriteria * giveRemeshingCrit()
Returns reference to associated remeshing criteria.
Definition: huertaerrorestimator.C:555
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition: error.h:61
oofem::HuertaErrorEstimator::mode
AnalysisMode mode
Linear analysis flag.
Definition: huertaerrorestimator.h:127
oofem::DynamicDataReader::finish
virtual void finish()
Allows to detach all data connections.
Definition: dynamicdatareader.C:84
oofem::HuertaErrorEstimator::solveRefinedElementProblem
void solveRefinedElementProblem(int elemId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, TimeStep *tStep)
Solves the refined element problem.
Definition: huertaerrorestimator.C:2760
oofem::HuertaErrorEstimator::setupRefinedProblemEpilog2
void setupRefinedProblemEpilog2(int tfuncs)
Definition: huertaerrorestimator.C:4268
oofem::HuertaErrorEstimator::setupRefinedProblemProlog
void setupRefinedProblemProlog(const char *problemName, int problemId, IntArray &localNodeIdArray, int nodes, int elems, int csects, int mats, int loads, int funcs, IntArray &controlNode, IntArray &controlDof, TimeStep *tStep)
Definition: huertaerrorestimator.C:3899
oofem::FloatArray::computeSquaredNorm
double computeSquaredNorm() const
Computes the square of the norm.
Definition: floatarray.C:846
oofem::ErrorEstimator::giveNumberOfSkippedElements
int giveNumberOfSkippedElements()
Returns number of elements skipped in error estimation.
Definition: errorestimator.h:125
oofem::HuertaRemeshingCriteria::nodalDensities
FloatArray nodalDensities
Array of nodal mesh densities.
Definition: huertaerrorestimator.h:293
_IFT_Domain_ndofman
#define _IFT_Domain_ndofman
Definition: domain.h:57
oofem::RefinedElement::giveBcDofArray1D
bool giveBcDofArray1D(int inode, Element *element, IntArray &sideBcDofId, int &sideNumBc, TimeStep *tStep)
Definition: refinedelement.C:150
oofem::HuertaErrorEstimatorInterface::setupRefinedElementProblem1D
void setupRefinedElementProblem1D(Element *element, RefinedElement *refinedElement, int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, HuertaErrorEstimatorInterface::SetupMode mode, TimeStep *tStep, int nodes, FloatArray **corner, FloatArray &midNode, int &localNodeId, int &localElemId, int &localBcId, IntArray &controlNode, IntArray &controlDof, HuertaErrorEstimator::AnalysisMode aMode, const char *edgetype)
Definition: huertaerrorestimator.C:1082
_IFT_Domain_axisymmetric
#define _IFT_Domain_axisymmetric
Definition: domain.h:71
oofem::FloatArray::restoreYourself
contextIOResultType restoreYourself(DataStream &stream)
Definition: floatarray.C:895
oofem::Dof::giveDofID
DofIDItem giveDofID() const
Returns DofID value of receiver, which determines type of of unknown connected to receiver (e...
Definition: dof.h:276
oofem::Domain::giveNumberOfCrossSectionModels
int giveNumberOfCrossSectionModels() const
Returns number of cross section models in domain.
Definition: domain.h:438
oofem::ConnectivityTable
Class representing connectivity table.
Definition: connectivitytable.h:55
oofem::Domain::giveNumberOfMaterialModels
int giveNumberOfMaterialModels() const
Returns number of material models in domain.
Definition: domain.h:436
oofem::Element::giveBodyLoadArray
IntArray * giveBodyLoadArray()
Returns array containing load numbers of loads acting on element.
Definition: element.C:372
oofem::FloatArray::beProductOf
void beProductOf(const FloatMatrix &aMatrix, const FloatArray &anArray)
Receiver becomes the result of the product of aMatrix and anArray.
Definition: floatarray.C:676
oofem::FEMComponent::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: femcmpnn.C:89
oofem::Load::giveInputRecord
virtual void giveInputRecord(DynamicInputRecord &input)
Setups the input record string of receiver.
Definition: load.C:129
_IFT_EngngModel_contextoutputstep
#define _IFT_EngngModel_contextoutputstep
Definition: engngm.h:69
oofem::HuertaErrorEstimator::globalENorm
double globalENorm
Global error norm.
Definition: huertaerrorestimator.h:105
_IFT_NRSolver_ddv
#define _IFT_NRSolver_ddv
Definition: nrsolver.h:60
oofem::RefinedElement::giveBoundaryLoadArray2D
bool giveBoundaryLoadArray2D(int inode, Element *element, std::vector< IntArray > &boundaryLoadList)
Definition: refinedelement.C:404
_IFT_EngngModel_nsteps
#define _IFT_EngngModel_nsteps
Definition: engngm.h:68
oofem::Element::giveElementDofIDMask
virtual void giveElementDofIDMask(IntArray &answer) const
Returns element dof mask for node.
Definition: element.h:498
oofem::FloatMatrix::at
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
_IFT_EngngModel_renumberFlag
#define _IFT_EngngModel_renumberFlag
Definition: engngm.h:70
_IFT_Domain_ncrosssect
#define _IFT_Domain_ncrosssect
Definition: domain.h:60
oofem::IntArray::resize
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
oofem::HuertaErrorEstimator::primaryUnknownError
FloatArray primaryUnknownError
Primary unknown nodal error.
Definition: huertaerrorestimator.h:119
oofem::RefinedElement::giveBcDofArray3D
bool giveBcDofArray3D(int inode, Element *element, std::vector< IntArray > &sideBcDofIdList, IntArray &sideNumBc, std::vector< IntArray > &faceBcDofIdList, IntArray &faceNumBc, TimeStep *tStep)
Definition: refinedelement.C:232
_IFT_Node_Name
#define _IFT_Node_Name
Definition: node.h:49
_IFT_HuertaErrorEstimator_impPos
#define _IFT_HuertaErrorEstimator_impPos
Definition: huertaerrorestimator.h:59
_IFT_HuertaRemeshingCriteria_werror
#define _IFT_HuertaRemeshingCriteria_werror
Definition: huertaerrorestimator.h:68
oofem::Domain::giveFunction
Function * giveFunction(int n)
Service for accessing particular domain load time function.
Definition: domain.C:268
_IFT_HuertaErrorEstimator_skipsteps
#define _IFT_HuertaErrorEstimator_skipsteps
Definition: huertaerrorestimator.h:54
oofem::ErrorEstimator
The base class for all error estimation or error indicator algorithms.
Definition: errorestimator.h:80
oofem::EngngModel::terminateAnalysis
void terminateAnalysis()
Performs analysis termination after finishing analysis.
Definition: engngm.C:1819
_IFT_HuertaErrorEstimator_requirederror
#define _IFT_HuertaErrorEstimator_requirederror
Definition: huertaerrorestimator.h:53
verbose.h
Initializes the variable VERBOSE, in order to get a few intermediate messages on screen: beginning an...
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::HuertaErrorEstimator::giveValue
virtual double giveValue(EE_ValueType type, TimeStep *tStep)
Returns the characteristic value of given type.
Definition: huertaerrorestimator.C:538
oofem::RemeshingStrategy
RemeshingStrategy
Type representing the remeshing strategy.
Definition: remeshingcrit.h:50
oofem::Element::giveParallelMode
elementParallelMode giveParallelMode() const
Return elementParallelMode of receiver.
Definition: element.h:1069
oofem::HuertaErrorEstimator::saveContext
virtual contextIOResultType saveContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Stores receiver state to output stream.
Definition: huertaerrorestimator.C:641
oofem::EngngModel::checkConsistency
virtual int checkConsistency()
Allows programmer to test some receiver&#39;s internal data, before computation begins.
Definition: engngm.h:995
oofem::Element::giveSpatialDimension
virtual int giveSpatialDimension()
Returns the element spatial dimension (1, 2, or 3).
Definition: element.C:1347
oofem::FloatMatrix
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
oofem::HuertaErrorEstimator::buildRefinedMesh
void buildRefinedMesh()
Builds refined mesh.
Definition: huertaerrorestimator.C:1056
ERROR_EXCESS
#define ERROR_EXCESS
Definition: huertaerrorestimator.C:82
oofem::ConnectivityTable::giveDofManConnectivityArray
const IntArray * giveDofManConnectivityArray(int dofman)
Definition: connectivitytable.C:95
oofem::globalNelems
static int globalNelems
Definition: huertaerrorestimator.C:132
oofem::EE_ErrorMode
EE_ErrorMode
Type determining whether temporary or equilibrated variables are used for error evaluation.
Definition: errorestimator.h:69
oofem::IRResultType
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
oofem::IntegrationRule::getIntegrationPoint
GaussPoint * getIntegrationPoint(int n)
Access particular integration point of receiver.
Definition: integrationrule.C:82
oofem::HuertaErrorEstimatorInterface::setupRefinedElementProblem2D
void setupRefinedElementProblem2D(Element *element, RefinedElement *refinedElement, int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, HuertaErrorEstimatorInterface::SetupMode mode, TimeStep *tStep, int nodes, FloatArray **corner, FloatArray *midSide, FloatArray &midNode, int &localNodeId, int &localElemId, int &localBcId, IntArray &controlNode, IntArray &controlDof, HuertaErrorEstimator::AnalysisMode aMode, const char *quadtype)
Definition: huertaerrorestimator.C:1492
_IFT_HuertaErrorEstimator_impCSect
#define _IFT_HuertaErrorEstimator_impCSect
Definition: huertaerrorestimator.h:58
oofem::HuertaErrorEstimator::eNorms
FloatArray eNorms
Cache storing element norms.
Definition: huertaerrorestimator.h:113
oofem::HuertaRemeshingCriteria::minElemSize
double minElemSize
Minimum element size alloved.
Definition: huertaerrorestimator.h:303
oofem::TimeStep::giveMetaStepNumber
int giveMetaStepNumber()
Returns receiver&#39;s meta step number.
Definition: timestep.h:135
oofem::GaussPoint::setNaturalCoordinates
void setNaturalCoordinates(const FloatArray &c)
Definition: gausspoint.h:139
_IFT_NonLinearStatic_controlmode
#define _IFT_NonLinearStatic_controlmode
Definition: nlinearstatic.h:44
oofem::FloatArray::computeNorm
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
oofem::ErrorEstimator::skipRegion
bool skipRegion(int reg)
Returns nonzero if region has been skipped in error estimation (user option).
Definition: errorestimator.C:79
oofem::InputRecord
Class representing the general Input Record.
Definition: inputrecord.h:101
oofem::Domain::giveNumberOfFunctions
int giveNumberOfFunctions() const
Returns number of load time functions in domain.
Definition: domain.h:444
oofem::ErrorEstimator::giveValue
virtual double giveValue(EE_ValueType type, TimeStep *tStep)=0
Returns the characteristic value of given type.
oofem::EIPrimaryUnknownMapper::mapAndUpdate
virtual int mapAndUpdate(FloatArray &answer, ValueModeType mode, Domain *oldd, Domain *newd, TimeStep *tStep)
Maps and updates the vector(s) of primary unknowns from old mesh oldd to new mesh newd...
Definition: eleminterpunknownmapper.C:52
oofem::HuertaRemeshingCriteria::refineCoeff
double refineCoeff
Refinement coefficient.
Definition: huertaerrorestimator.h:305
oofem::DofManager::giveDofWithID
Dof * giveDofWithID(int dofID) const
Returns DOF with given dofID; issues error if not present.
Definition: dofmanager.C:119
oofem::IntArray::followedBy
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
oofem::Timer
Class implementing single timer, providing wall clock and user time capabilities. ...
Definition: timer.h:46
_IFT_Domain_nbc
#define _IFT_Domain_nbc
Definition: domain.h:61
oofem::EngngModel::giveRank
int giveRank() const
Returns domain rank in a group of collaborating processes (0..groupSize-1)
Definition: engngm.h:1058
oofem::RefinedElement::giveBcDofArray2D
bool giveBcDofArray2D(int inode, Element *element, std::vector< IntArray > &sideBcDofIdList, IntArray &sideNumBc, TimeStep *tStep)
Definition: refinedelement.C:181
oofem::Domain::giveNumberOfNonlocalBarriers
int giveNumberOfNonlocalBarriers() const
Returns number of nonlocal integration barriers.
Definition: domain.h:448
oofem::FEMComponent::giveInterface
virtual Interface * giveInterface(InterfaceType t)
Interface requesting service.
Definition: femcmpnn.h:179
oofem::DynamicInputRecord
Class representing the a dynamic Input Record.
Definition: dynamicinputrecord.h:59
_IFT_Node_coords
#define _IFT_Node_coords
Definition: node.h:50
_IFT_OutputManager_Name
#define _IFT_OutputManager_Name
Definition: outputmanager.h:47
_IFT_HuertaRemeshingCriteria_requirederror
#define _IFT_HuertaRemeshingCriteria_requirederror
Definition: huertaerrorestimator.h:65
oofem::Element::computeLocalCoordinates
virtual bool computeLocalCoordinates(FloatArray &answer, const FloatArray &gcoords)
Computes the element local coordinates from given global coordinates.
Definition: element.C:1222
oofem::ContextMode
long ContextMode
Context mode (mask), defining the type of information written/read to/from context.
Definition: contextmode.h:43
oofem::Element_remote
Element in active domain is only mirror of some remote element.
Definition: element.h:102
_IFT_BoundaryCondition_PrescribedValue
#define _IFT_BoundaryCondition_PrescribedValue
[rn,optional] Prescribed value of all DOFs
Definition: boundarycondition.h:49
oofem::DofManager::giveLoadArray
IntArray * giveLoadArray()
Returns the array containing applied loadings of the receiver.
Definition: dofmanager.C:82
oofem::DynamicInputRecord::setRecordKeywordField
void setRecordKeywordField(std::string keyword, int number)
Definition: dynamicinputrecord.C:288
_IFT_HeavisideTimeFunction_value
#define _IFT_HeavisideTimeFunction_value
Definition: heavisidetimefunction.h:45
oofem::EE_ValueType
EE_ValueType
Type characterizing different type of errors.
Definition: errorestimator.h:65
oofem::Dof::hasBc
virtual bool hasBc(TimeStep *tStep)=0
Test if Dof has active boundary condition.
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition: femcmpnn.h:100
oofem::HuertaErrorEstimatorInterface::HuertaErrorEstimatorI_setupRefinedElementProblem
virtual void HuertaErrorEstimatorI_setupRefinedElementProblem(RefinedElement *refinedElement, int level, int nodeId, IntArray &localNodeIdArray, IntArray &globalNodeIdArray, HuertaErrorEstimatorInterface::SetupMode mode, TimeStep *tStep, int &localNodeId, int &localElemId, int &localBcId, IntArray &controlNode, IntArray &controlDof, HuertaErrorEstimator::AnalysisMode aMode)=0
oofem::min
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
oofem::FloatArray::push_back
void push_back(const double &iVal)
Add one element.
Definition: floatarray.h:118
oofem::FEMComponent::restoreContext
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode, void *obj=NULL)
Restores the receiver state previously written in stream.
Definition: femcmpnn.C:64
oofem::HuertaRemeshingCriteria::HuertaRemeshingCriteria
HuertaRemeshingCriteria(int n, ErrorEstimator *e)
Constructor.
Definition: huertaerrorestimator.C:691
oofem::HuertaRemeshingCriteria::requiredError
double requiredError
Required error to obtain.
Definition: huertaerrorestimator.h:301
oofem::exactFlag
static bool exactFlag
Definition: huertaerrorestimator.C:111
oofem::impPos
static FloatArray impPos
Definition: huertaerrorestimator.C:130
oofem::EngngModel
Abstract base class representing the "problem" under consideration.
Definition: engngm.h:181
oofem::DofManager_remote
DofManager in active domain is only mirror of some remote DofManager.
Definition: dofmanager.h:85
IR_GIVE_OPTIONAL_FIELD
#define IR_GIVE_OPTIONAL_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:78
oofem::EngngModel::giveCurrentStep
virtual TimeStep * giveCurrentStep(bool force=false)
Returns current time step.
Definition: engngm.h:683
oofem::IntArray::giveSize
int giveSize() const
Definition: intarray.h:203
oofem::FloatArray::giveSize
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
oofem::EngngModel::restoreContext
virtual contextIOResultType restoreContext(DataStream &stream, ContextMode mode)
Restores the state of model from output stream.
Definition: engngm.C:1592
_IFT_EngngModel_nmsteps
#define _IFT_EngngModel_nmsteps
Definition: engngm.h:72
_IFT_HeavisideTimeFunction_Name
#define _IFT_HeavisideTimeFunction_Name
Definition: heavisidetimefunction.h:43
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::ErrorEstimator::giveElementError
virtual double giveElementError(EE_ErrorType type, Element *elem, TimeStep *tStep)=0
Returns the element error.
oofem::Element::giveDofManager
DofManager * giveDofManager(int i) const
Definition: element.C:514
oofem::RefinedElement::giveBoundaryLoadArray1D
bool giveBoundaryLoadArray1D(int inode, Element *element, IntArray &boundaryLoadArray)
Definition: refinedelement.C:367
STIFFNESS_TYPE
#define STIFFNESS_TYPE
Definition: huertaerrorestimator.C:79
_IFT_HuertaRemeshingCriteria_minelemsize
#define _IFT_HuertaRemeshingCriteria_minelemsize
Definition: huertaerrorestimator.h:66
oofem::EngngModel::giveDomain
Domain * giveDomain(int n)
Service for accessing particular problem domain.
Definition: engngm.C:1720
oofem::Node
Class implementing node in finite element mesh.
Definition: node.h:87
_IFT_NRSolver_ddfunc
#define _IFT_NRSolver_ddfunc
Definition: nrsolver.h:61
IR_GIVE_FIELD
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
_IFT_DofManager_load
#define _IFT_DofManager_load
Definition: dofmanager.h:53
oofem::Dof
Abstract class Dof represents Degree Of Freedom in finite element mesh.
Definition: dof.h:93
_IFT_MetaStep_Name
#define _IFT_MetaStep_Name
Definition: metastep.h:43
oofem::FEMComponent::giveNumber
int giveNumber() const
Definition: femcmpnn.h:107
_IFT_NRSolver_minsteplength
#define _IFT_NRSolver_minsteplength
Definition: nrsolver.h:56
oofem::Domain::giveDofManager
DofManager * giveDofManager(int n)
Service for accessing particular domain dof manager.
Definition: domain.C:314
_IFT_HuertaErrorEstimator_perfectCSect
#define _IFT_HuertaErrorEstimator_perfectCSect
Definition: huertaerrorestimator.h:57
_IFT_CylindricalALM_hpcw
#define _IFT_CylindricalALM_hpcw
Definition: calmls.h:63
oofem::Element::giveNode
Node * giveNode(int i) const
Returns reference to the i-th node of element.
Definition: element.h:610
oofem::HuertaRemeshingCriteria::remeshingStrategy
RemeshingStrategy remeshingStrategy
Remeshing strategy proposed.
Definition: huertaerrorestimator.h:295
oofem::Domain::giveLoad
Load * giveLoad(int n)
Service for accessing particular domain load.
Definition: domain.C:222
_IFT_HuertaErrorEstimator_werror
#define _IFT_HuertaErrorEstimator_werror
Definition: huertaerrorestimator.h:56
oofem::impCSect
static int impCSect
Definition: huertaerrorestimator.C:129
oofem::Timer::startTimer
void startTimer()
Definition: timer.C:69
oofem::ContextIOERR
Context IO exception class.
Definition: contextioerr.h:46
oofem::Timer::getUtime
double getUtime()
Returns total user time elapsed in seconds.
Definition: timer.C:105
oofem::GaussPoint
Class representing integration point in finite element program.
Definition: gausspoint.h:93
oofem::DynamicDataReader::insertInputRecord
void insertInputRecord(InputRecordType type, InputRecord *record)
Main purpose of this class it the possibility to add new input records in code.
Definition: dynamicdatareader.C:53
oofem::TimeStep
Class representing solution step.
Definition: timestep.h:80
oofem::HuertaErrorEstimator::globalUNorm
double globalUNorm
Global norm of primary unknown.
Definition: huertaerrorestimator.h:109
oofem::DofManager_shared
DofManager is shared by neighboring partitions, it is necessary to sum contributions from all contrib...
Definition: dofmanager.h:82
oofem::DofManager::giveParallelMode
dofManagerParallelMode giveParallelMode() const
Return dofManagerParallelMode of receiver.
Definition: dofmanager.h:512
oofem::AdaptiveNonLinearStatic
This class implements Adaptive Non-LinearStatic Engineering problem.
Definition: adaptnlinearstatic.h:74
oofem::Timer::stopTimer
void stopTimer()
Definition: timer.C:77
oofem::HuertaRemeshingCriteria::giveDofManDensity
virtual double giveDofManDensity(int num)
Returns existing mesh size for given dof manager.
Definition: huertaerrorestimator.C:1018
_IFT_Node_lcs
#define _IFT_Node_lcs
Definition: node.h:51
_IFT_CylindricalALM_manrmsteps
#define _IFT_CylindricalALM_manrmsteps
Definition: calmls.h:60
oofem::HuertaRemeshingCriteria::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: huertaerrorestimator.C:987
oofem::DynamicDataReader::setOutputFileName
void setOutputFileName(const std::string &outputFileName)
Sets the output file name.
Definition: dynamicdatareader.h:81
oofem::HuertaErrorEstimator::refineLevel
int refineLevel
Refinement level.
Definition: huertaerrorestimator.h:121
oofem::DynamicDataReader::setDescription
void setDescription(const std::string &description)
Sets the description line.
Definition: dynamicdatareader.h:85
_IFT_Domain_numberOfSpatialDimensions
#define _IFT_Domain_numberOfSpatialDimensions
[in,optional] Specifies how many spatial dimensions the domain has.
Definition: domain.h:69
oofem::RefinedElement::giveBoundaryFlagArray
IntArray * giveBoundaryFlagArray(void)
Definition: refinedelement.h:61
_IFT_CylindricalALM_steplength
#define _IFT_CylindricalALM_steplength
Definition: calmls.h:55
oofem::FloatArray::resize
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
oofem::perCSect
static int perCSect
Definition: huertaerrorestimator.C:129
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:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011