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