416 {
418
421 "Call stressToSpectralForm before evaluating the scalar "
422 "equation");
423 }
424 if (!std::isfinite(y)) {
426 "Non-finite Storakers scalar iterate");
427 }
428
429 state = {};
430 state.y = y;
433 const double log_minimum =
434 std::log(std::numeric_limits<double>::min());
435 const double log_maximum =
436 std::log(std::numeric_limits<double>::max());
437 const double negative_infinity =
438 -std::numeric_limits<double>::infinity();
439 const double log_linear_gamma =
440 mu1 > 0 ? std::log(2 *
mu1) + y : negative_infinity;
441 const double log_quadratic_gamma = std::log(
mu2) + 2 * y;
442 const double log_gamma =
443 logSumExp(log_linear_gamma, log_quadratic_gamma);
445 log_linear_gamma, std::log(2 *
mu2) + 2 * y);
446
447 int gamma_effective_sign = 1;
448 double log_abs_gamma_effective = log_gamma;
450 log_abs_gamma_effective =
454 if (log_gamma > log_minus_q) {
455 log_abs_gamma_effective =
456 log_gamma +
457 std::log(-std::expm1(log_minus_q - log_gamma));
458 } else if (log_gamma < log_minus_q) {
459 gamma_effective_sign = -1;
460 log_abs_gamma_effective =
461 log_minus_q +
462 std::log(-std::expm1(log_gamma - log_minus_q));
463 } else {
464 gamma_effective_sign = 0;
465 log_abs_gamma_effective = negative_infinity;
466 }
467 }
468
469 if (y > log_maximum || log_gamma > log_maximum ||
470 log_abs_gamma_effective > log_maximum ||
471 log_gamma_y > log_maximum) {
473 "Storakers scalar iterate is outside the floating-point "
474 "range");
475 }
476
477 state.z = y >= log_minimum ? std::exp(y) : 0;
478 state.gamma = log_gamma >= log_minimum ? std::exp(log_gamma) : 0;
479 state.gammaY =
480 log_gamma_y >= log_minimum ? std::exp(log_gamma_y) : 0;
482 state.gammaEffectiveSign = gamma_effective_sign;
483 state.logAbsGammaEffective = log_abs_gamma_effective;
484 if (gamma_effective_sign &&
485 log_abs_gamma_effective >= log_minimum) {
486 state.gammaEffective =
487 gamma_effective_sign * std::exp(log_abs_gamma_effective);
488 }
489
490 const double log_two_mu2 = std::log(2.) + std::log(
mu2);
491 double positive_gamma_term = 0;
492 if (gamma_effective_sign > 0) {
493 const double log_gamma_term =
494 std::log(2.) +
495 0.5 * (std::log(
mu2) + log_abs_gamma_effective);
496 if (log_gamma_term > log_maximum) {
498 "Storakers principal-stretch discriminant is outside the "
499 "floating-point range");
500 }
501 positive_gamma_term =
502 log_gamma_term >= log_minimum ? std::exp(log_gamma_term) : 0;
503 }
504
505 double log_derivative_sum = negative_infinity;
506 for (
int ii = 0; ii !=
SPACE_DIM; ++ii) {
507 double log_s = negative_infinity;
508 if (gamma_effective_sign > 0) {
509 state.tS(ii) = std::hypot(
tB(ii), positive_gamma_term);
510 log_s = std::log(state.tS(ii));
511 } else if (gamma_effective_sign == 0) {
514 "No positive stable Storakers stretch at gamma_eff=0");
515 }
516 state.tS(ii) =
tB(ii);
517 log_s = std::log(
tB(ii));
518 } else {
521 "Negative Storakers gamma_eff requires positive shifted "
522 "principal stresses");
523 }
524 const double log_b = std::log(
tB(ii));
525 const double log_negative_term =
526 std::log(4.) + std::log(
mu2) + log_abs_gamma_effective;
527 const double log_b_squared = 2 * log_b;
528 if (!(log_negative_term < log_b_squared)) {
530 "Storakers state is outside the stable principal-stretch "
531 "branch");
532 }
533 log_s =
534 log_b +
535 0.5 * std::log(-std::expm1(log_negative_term -
536 log_b_squared));
537 state.tS(ii) = log_s >= log_minimum ? std::exp(log_s) : 0;
538 }
539 if (!std::isfinite(state.tS(ii)) || state.tS(ii) <= 0) {
541 "Invalid Storakers principal-stretch discriminant");
542 }
543
545 const double log_b =
546 tB(ii) > 0 ? std::log(
tB(ii)) : negative_infinity;
547 state.tLogU(ii) =
549 } else {
550 const double log_denominator =
552 state.tLogU(ii) =
553 std::log(2.) + log_abs_gamma_effective - log_denominator;
554 }
555
556 if (!std::isfinite(state.tLogU(ii)) ||
557 state.tLogU(ii) < log_minimum ||
558 state.tLogU(ii) > log_maximum) {
560 "Invalid Storakers principal stretch");
561 }
562 state.tU(ii) = std::exp(state.tLogU(ii));
563
564 const double log_derivative_term =
565 -state.tLogU(ii) - log_s;
566 log_derivative_sum =
567 logSumExp(log_derivative_sum, log_derivative_term);
568 }
569
570 state.phi =
572 const double log_volumetric_derivative =
573 log_gamma_y + log_derivative_sum;
574 if (log_volumetric_derivative > log_maximum) {
576 "Invalid Storakers scalar derivative term");
577 }
578 const double volumetric_derivative =
579 log_volumetric_derivative >= log_minimum
580 ? std::exp(log_volumetric_derivative)
581 : 0;
583 if (!std::isfinite(state.phi) || !std::isfinite(state.dPhi) ||
584 state.dPhi <= 0) {
586 "Invalid Storakers scalar residual or derivative");
587 }
588
590 }
#define FTENSOR_INDEX(DIM, I)
FTensor::Tensor1< double, SPACE_DIM > tB
