1298 {
1300
1303
1304
1305 auto dm =
simple->getDM();
1308
1309 auto set_section_monitor = [&](auto solver) {
1311 SNES snes;
1312 CHKERR TSGetSNES(solver, &snes);
1313 CHKERR SNESMonitorSet(snes,
1316 (void *)(snes_ctx_ptr.get()), nullptr);
1318 };
1319
1320 auto create_post_process_elements = [&]() {
1321 auto block_thermal_params = boost::make_shared<BlockedThermalParameters>();
1322
1323 auto block_thermoelastic_params =
1324 boost::make_shared<BlockedThermoElasticParameters>();
1325 auto coeff_expansion_ptr =
1326 block_thermoelastic_params->getCoeffExpansionPtr();
1327 auto ref_temp_ptr = block_thermoelastic_params->getRefTempPtr();
1328
1329 auto u_ptr = boost::make_shared<MatrixDouble>();
1330 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1331 auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
1332 auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
1333 auto vec_temp_ptr = boost::make_shared<VectorDouble>();
1334 auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
1335
1336 auto push_domain_ops = [&](auto &pp_fe) {
1338 auto &pip = pp_fe->getOpPtrVector();
1339
1341 Sev::verbose);
1343 pip, "MAT_THERMOELASTIC", block_thermoelastic_params, Sev::verbose);
1344
1346
1348 pip.push_back(
1350
1353 "U", mat_grad_ptr));
1354 auto elastic_common_ptr = commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
1355 mField, pip,
"U",
"MAT_ELASTIC", Sev::inform);
1357 pip.push_back(
1358 new
1359 typename H::template OpCalculateHenckyThermalStress<SPACE_DIM, IT, 0>(
1360 "U", vec_temp_ptr, elastic_common_ptr, coeff_expansion_ptr,
1361 ref_temp_ptr));
1364 elastic_common_ptr->getMatFirstPiolaStress(), mat_stress_ptr));
1365 pip.push_back(
1367 } else {
1368 mat_stress_ptr = elastic_common_ptr->getMatFirstPiolaStress();
1369 mat_strain_ptr = elastic_common_ptr->getMatLogC();
1370 }
1371
1373 };
1374
1375 auto push_post_proc_ops = [&](auto &pp_fe) {
1377 auto &pip = pp_fe->getOpPtrVector();
1379
1381 pip.push_back(
1382
1384
1385 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1386
1387 {{"T", vec_temp_ptr}},
1388
1389 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1390
1391 {},
1392
1393 {{"CAUCHY", mat_stress_ptr}, {"STRAIN", mat_strain_ptr}}
1394
1395 )
1396
1397 );
1398 } else {
1399 pip.push_back(
1400
1402
1403 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
1404
1405 {{"T", vec_temp_ptr}},
1406
1407 {{"U", u_ptr}, {"FLUX", mat_flux_ptr}},
1408
1409 {{"PIOLA", mat_stress_ptr}},
1410
1411 {{"HENCKY_STRAIN", mat_strain_ptr}}
1412
1413 )
1414
1415 );
1416 }
1417
1419 };
1420
1421 auto domain_post_proc = [&]() {
1423 return boost::shared_ptr<PostProcEle>();
1424 auto pp_fe = boost::make_shared<PostProcEle>(
mField);
1426 "push domain ops to domain element");
1428 "push post proc ops to domain element");
1429 return pp_fe;
1430 };
1431
1432 auto skin_post_proc = [&]() {
1434 return boost::shared_ptr<SkinPostProcEle>();
1435 auto pp_fe = boost::make_shared<SkinPostProcEle>(
mField);
1440 "push domain ops to side element");
1441 pp_fe->getOpPtrVector().push_back(op_side);
1443 "push post proc ops to skin element");
1444 return pp_fe;
1445 };
1446
1447 return std::make_pair(domain_post_proc(), skin_post_proc());
1448 };
1449
1450 auto monitor_ptr = boost::make_shared<FEMethod>();
1451
1452 auto set_time_monitor = [&](auto dm, auto solver, auto domain_post_proc_fe,
1453 auto skin_post_proc_fe) {
1455 monitor_ptr->preProcessHook = [&]() {
1457
1461 domain_post_proc_fe,
1462 monitor_ptr->getCacheWeakPtr());
1463 CHKERR domain_post_proc_fe->writeFile(
1464 "out_" + boost::lexical_cast<std::string>(monitor_ptr->ts_step) +
1465 ".h5m");
1466 }
1469 skin_post_proc_fe,
1470 monitor_ptr->getCacheWeakPtr());
1471 CHKERR skin_post_proc_fe->writeFile(
1472 "out_skin_" +
1473 boost::lexical_cast<std::string>(monitor_ptr->ts_step) + ".h5m");
1474 }
1475 }
1476
1477 struct AtomTestResult {
1478 bool fail = false;
1479 std::string msg = "";
1480 };
1481 AtomTestResult atom_test_result;
1482
1483 auto fail_atom_test = [&atom_test_result](const std::string &msg) {
1484 atom_test_result.fail = true;
1485 atom_test_result.msg = msg;
1486 };
1487
1488 struct AtomTestData{
1489 double expected = 0.0;
1491 };
1492
1494
1496 ->evalFEAtThePoint<SPACE_DIM>(
1501
1503 auto eval_num_vec =
1505 CHKERR VecZeroEntries(eval_num_vec);
1507 CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
1508 }
1509 CHKERR VecAssemblyBegin(eval_num_vec);
1510 CHKERR VecAssemblyEnd(eval_num_vec);
1511
1512 double eval_num;
1513 CHKERR VecSum(eval_num_vec, &eval_num);
1514 if (!(int)eval_num) {
1515 fail_atom_test(
1516 "did not find elements to evaluate the field, check the "
1517 "coordinates");
1518 }
1519 }
1520
1523 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1524 << "Eval point T: " << t_temp;
1525 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1526 AtomTestData atom_test_data;
1528 case 1:
1529 case 2:
1530 case 3:
1531 atom_test_data = {554.48, 1e-2};
1532 break;
1533 case 4:
1534 atom_test_data = {325.0, 2e-2};
1535 break;
1536 case 5:
1537 atom_test_data = {1.0, 1e-2};
1538 break;
1539 default:
1540 fail_atom_test("unknown atom test number");
1541 }
1542 if (fabs(t_temp - atom_test_data.expected) > atom_test_data.tol) {
1543 fail_atom_test("wrong temperature value");
1544 }
1545 }
1546 }
1549 auto t_flux = getFTensor1FromMat<SPACE_DIM>(*
fluxFieldPtr);
1550 auto flux_mag = sqrt(t_flux(
i) * t_flux(
i));
1551 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1552 << "Eval point FLUX magnitude: " << flux_mag;
1553 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1554 AtomTestData atom_test_data;
1556 case 1:
1557 case 2:
1558 case 3:
1559 atom_test_data = {27008.0, 2e1};
1560 break;
1561 case 4:
1562 atom_test_data = {150e3, 2.1e1};
1563 break;
1564 case 5:
1565 atom_test_data = {0.0, 1e-6};
1566 break;
1567 default:
1568 fail_atom_test("unknown atom test number");
1569 }
1570 if (fabs(flux_mag - atom_test_data.expected) > atom_test_data.tol) {
1571 fail_atom_test("wrong flux value");
1572 }
1573 }
1574 }
1577 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*
dispFieldPtr);
1578 auto disp_mag = sqrt(t_disp(
i) * t_disp(
i));
1579 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1580 << "Eval point U magnitude: " << disp_mag;
1581 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1582 AtomTestData atom_test_data;
1584 case 1:
1585 atom_test_data = {0.00345, 1e-5};
1586 break;
1587 case 2:
1588 case 3:
1589 atom_test_data = {0.00265, 1e-5};
1590 break;
1591 case 4:
1592 atom_test_data = {0.00075, 1e-5};
1593 break;
1594 case 5:
1595 atom_test_data = {std::sqrt(2) * (std::sqrt(std::exp(0.2)) - 1),
1596 1e-5};
1597 break;
1598 default:
1599 fail_atom_test("unknown atom test number");
1600 }
1601 if (fabs(disp_mag - atom_test_data.expected) > atom_test_data.tol) {
1602 fail_atom_test("wrong displacement value");
1603 }
1604 }
1605 }
1608 auto t_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(*
strainFieldPtr);
1609 auto t_strain_trace = t_strain(
i,
i);
1610 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1611 AtomTestData atom_test_data;
1613 case 1:
1614 atom_test_data = {0.00679, 1e-5};
1615 break;
1616 case 2:
1617 case 3:
1618 atom_test_data = {0.00522, 1e-5};
1619 break;
1620 case 4:
1621 atom_test_data = {0.00055, 1e-5};
1622 break;
1623 case 5:
1624 atom_test_data = {0.2, 1e-5};
1625 break;
1626 default:
1627 fail_atom_test("unknown atom test number");
1628 }
1629 if (fabs(t_strain_trace - atom_test_data.expected) >
1630 atom_test_data.tol) {
1631 fail_atom_test("wrong strain value");
1632 }
1633 }
1634 }
1636 double von_mises_stress;
1637 auto getVonMisesStress = [&](auto t_stress) {
1639 PetscBool plane_strain_flag = PETSC_FALSE;
1641 &plane_strain_flag, PETSC_NULLPTR);
1642
1643 const double s33 =
1645 ? t_stress(2, 2)
1647 (t_stress(0, 0) + t_stress(1, 1))
1648 : 0.0);
1649 von_mises_stress = std::sqrt(
1650 0.5 * ((t_stress(0, 0) - t_stress(1, 1)) *
1651 (t_stress(0, 0) - t_stress(1, 1)) +
1652 (t_stress(1, 1) - s33) * (t_stress(1, 1) - s33) +
1653 (s33 - t_stress(0, 0)) * (s33 - t_stress(0, 0)) +
1654 6.0 * (t_stress(0, 1) * t_stress(0, 1) +
1655 (
SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2)
1656 : 0.0) +
1657 (
SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0)
1658 : 0.0))));
1660 };
1662 auto t_stress = getFTensor2SymmetricFromMat<SPACE_DIM>(*
stressFieldPtr);
1663 CHKERR getVonMisesStress(t_stress);
1664 } else {
1665 auto t_stress =
1667 CHKERR getVonMisesStress(t_stress);
1668 }
1669 MOFEM_LOG(
"ThermoElasticSync", Sev::inform)
1670 << "Eval point von Mises Stress: " << von_mises_stress;
1671 if (
atom_test && fabs(monitor_ptr->ts_t - 10) < 1e-12) {
1672 AtomTestData atom_test_data;
1674 case 1:
1675 atom_test_data = {523.0, 5e-1};
1676 break;
1677 case 2:
1678 atom_test_data = {16.3, 6e-2};
1679 break;
1680 case 3:
1681 atom_test_data = {14.9, 5e-2};
1682 break;
1683 case 4:
1684 atom_test_data = {69.2, 1e-1};
1685 break;
1686 case 5:
1687 atom_test_data = {0.0, 5e-2};
1688 break;
1689 default:
1690 fail_atom_test("unknown atom test number");
1691 }
1692 if (fabs(von_mises_stress - atom_test_data.expected) >
1693 atom_test_data.tol) {
1694 fail_atom_test("wrong von Mises stress value");
1695 }
1696 }
1697 }
1699 }
1700
1701 if (atom_test_result.fail) {
1704 atom_test_result.msg.c_str());
1705 }
1706
1708 };
1709 auto null = boost::shared_ptr<FEMethod>();
1711 monitor_ptr, null);
1713 };
1714
1715 auto set_fieldsplit_preconditioner = [&](auto solver) {
1717
1718 SNES snes;
1719 CHKERR TSGetSNES(solver, &snes);
1720 KSP ksp;
1721 CHKERR SNESGetKSP(snes, &ksp);
1722 PC pc;
1723 CHKERR KSPGetPC(ksp, &pc);
1724 PetscBool is_pcfs = PETSC_FALSE;
1725 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1726
1727
1728 if (is_pcfs == PETSC_TRUE) {
1730 auto name_prb =
simple->getProblemName();
1731
1733 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"U", 0,
1736 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"FLUX", 0, 0,
1737 is_flux);
1739 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"T", 0, 0,
1740 is_T);
1742 CHKERR is_mng->isCreateProblemFieldAndRank(name_prb,
ROW,
"TBC", 0, 0,
1743 is_TBC);
1744 IS is_tmp, is_tmp2;
1745 CHKERR ISExpand(is_T, is_flux, &is_tmp);
1746 CHKERR ISExpand(is_TBC, is_tmp, &is_tmp2);
1747 CHKERR ISDestroy(&is_tmp);
1749
1752 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_TFlux);
1753 CHKERR PCFieldSplitSetIS(pc, PETSC_NULLPTR, is_u);
1754 }
1755
1757 };
1758
1760 CHKERR TSSetIJacobian(solver,
B,
B, PETSC_NULLPTR, PETSC_NULLPTR);
1763 CHKERR TSSetSolution(solver,
D);
1764 CHKERR TSSetFromOptions(solver);
1765
1766 CHKERR set_section_monitor(solver);
1767 CHKERR set_fieldsplit_preconditioner(solver);
1768
1769 auto [domain_post_proc_fe, skin_post_proc_fe] =
1770 create_post_process_elements();
1771 CHKERR set_time_monitor(dm, solver, domain_post_proc_fe, skin_post_proc_fe);
1772
1774 CHKERR TSSolve(solver, NULL);
1775
1777}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Section manager is used to create indexes and sections.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
intrusive_ptr for managing petsc objects
