#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#include <cstdlib>
extern "C" {
}
};
};
IntegrationType::GAUSS;
boost::function<
double(
const double,
const double,
const double)>;
}
}
}
}
}
}
return std::exp(-
b_iso * tau);
}
double JC_ref_temp = 0.;
double JC_melt_temp = 1650.;
return (
temp - JC_ref_temp) / (JC_melt_temp - JC_ref_temp);
}
}
}
}
double JC_ref_temp = 0.;
double JC_melt_temp = 1650.;
(JC_melt_temp - JC_ref_temp);
}
template <typename T, int DIM>
inline auto
if (
C1_k < std::numeric_limits<double>::epsilon()) {
return t_alpha;
}
t_alpha(
i,
j) =
C1_k * t_plastic_strain(
i,
j);
return t_alpha;
}
template <int DIM>
return t_diff;
}
const char *
const ICTypes[] = {
"uniform",
"gaussian",
"linear",
"ICType", "IC_", nullptr};
11.4;
0.9;
double y, double z) {
double s =
};
double y, double z) {
};
double z) {
default:
}
};
#include <HenckyOps.hpp>
#include <PlasticOps.hpp>
#include <PlasticNaturalBCs.hpp>
#include <ThermoElasticOps.hpp>
#include <FiniteThermalOps.hpp>
#include <ThermalOps.hpp>
#include <ThermalConvection.hpp>
#include <ThermalRadiation.hpp>
#ifdef ADD_CONTACT
#ifdef ENABLE_PYTHON_BINDING
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
};
#include <ContactOps.hpp>
#endif
std::string output_name, int &counter = *(new int(0))) {
auto create_post_process_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(m_field);
auto push_vol_post_proc_ops = [&](auto &pp_fe) {
auto &pip = pp_fe->getOpPtrVector();
auto TAU_ptr = boost::make_shared<VectorDouble>();
auto T_ptr = boost::make_shared<VectorDouble>();
auto T_grad_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
auto U_ptr = boost::make_shared<MatrixDouble>();
auto FLUX_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
auto EP_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"TAU", TAU_ptr}, {"T", T_ptr}},
{{"U", U_ptr}, {"FLUX", FLUX_ptr}, {"T_GRAD", T_grad_ptr}},
{},
{{"EP", EP_ptr}}
)
);
};
CHKERR push_vol_post_proc_ops(pp_fe);
CHKERR pp_fe->writeFile(
"out_" + std::to_string(counter) +
"_" +
output_name + ".h5m");
} else {
auto &post_proc_moab = pp_fe->getPostProcMesh();
auto file_name = "out_" + std::to_string(counter) + "_" + output_name +
<< "Writing file " << file_name << std::endl;
CHKERR post_proc_moab.write_file(file_name.c_str(),
"VTK");
}
counter++;
};
CHKERR create_post_process_elements();
}
Vec sol, std::string output_name) {
auto create_post_process_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(m_field);
auto &pip = pp_fe->getOpPtrVector();
auto push_vol_post_proc_ops = [&](auto &pp_fe) {
auto &pip = pp_fe->getOpPtrVector();
auto TAU_ptr = boost::make_shared<VectorDouble>();
pip.push_back(
auto T_ptr = boost::make_shared<VectorDouble>();
auto U_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
auto FLUX_ptr = boost::make_shared<MatrixDouble>();
"FLUX", FLUX_ptr, smart_sol));
auto EP_ptr = boost::make_shared<MatrixDouble>();
"EP", EP_ptr, smart_sol));
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"TAU", TAU_ptr}, {"T", T_ptr}},
{{"U", U_ptr}, {"FLUX", FLUX_ptr}},
{},
{{"EP", EP_ptr}}
)
);
};
CHKERR push_vol_post_proc_ops(pp_fe);
CHKERR pp_fe->writeFile(
"out_" + output_name +
".h5m");
};
CHKERR create_post_process_elements();
}
const std::string &out_put_string,
const auto &out_put_quantity) {
for (
int i = 0;
i < size; ++
i) {
std::cout << "Proc " << rank << " " + out_put_string + " "
<< out_put_quantity << std::endl;
}
}
};
DomainRhsBCs::OpFlux<PlasticOps::DomainBCs, 1, SPACE_DIM>;
BoundaryRhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
BoundaryLhsBCs::OpFlux<PlasticOps::BoundaryBCs, 1, SPACE_DIM>;
DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1,
SPACE_DIM>;
DomainNaturalBCRhs::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 1>;
using OpForce = BoundaryNaturalBC::OpFlux<NaturalForceMeshsets, 1, SPACE_DIM>;
BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
IT>::OpEssentialRhs<HeatFluxCubitBcData, 3, SPACE_DIM>;
DomainNaturalBCRhs::OpFlux<SetTargetTemperature, 1, 1>;
DomainNaturalBCLhs::OpFlux<SetTargetTemperature, 1, 1>;
auto save_range = [](moab::Interface &moab,
const std::string name,
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1);
};
private:
};
};
static std::unordered_map<TS, ResizeCtx *>
ts_to_rctx;
: newEls(new_els) {}
int order_row, int order_col,
int order_data) {
constexpr int numNodes = 3;
constexpr int numEdges = 3;
auto &m_field = fe_raw_ptr->mField;
auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
auto fe_handle = fe_ptr->getFEEntityHandle();
auto set_base_quadrature = [&]() {
int rule = 2 * (order_data + 1);
}
}
fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
&fe_ptr->gaussPts(0, 0), 1);
&fe_ptr->gaussPts(1, 0), 1);
&fe_ptr->gaussPts(2, 0), 1);
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
false);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
1);
data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
std::copy(
Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
} else {
}
};
auto set_gauss_pts = [&](auto &ref_gauss_pts) {
fe_ptr->gaussPts.swap(ref_gauss_pts);
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
auto &data = fe_ptr->dataOnElement[
H1];
data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3);
double *shape_ptr =
&*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
&fe_ptr->gaussPts(1, 0), nb_gauss_pts);
};
auto refine_quadrature = [&]() {
moab::Core moab_ref;
double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
for (int nn = 0; nn != numNodes; nn++)
CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
CHKERR m_field_ref.
get_moab().get_adjacencies(ref_tri, 1,
true, edges,
moab::Interface::UNION);
for (int nn = 0; nn != numNodes; nn++) {
CHKERR moab_ref.side_element(tri, 0, nn, ent);
nodes_at_front.insert(ent);
}
CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
for (int ee = 0; ee != numEdges; ee++) {
CHKERR moab_ref.side_element(tri, 1, ee, ent);
CHKERR moab_ref.add_entities(meshset, &ent, 1);
}
CHKERR m_ref->addVerticesInTheMiddleOfEdges(ref_edges,
->updateMeshsetByEntitiesChildren(meshset,
BitRefLevel().set(1),
meshset, MBEDGE, true);
->updateMeshsetByEntitiesChildren(meshset,
BitRefLevel().set(1),
meshset, MBTRI, true);
->getEntitiesByTypeAndRefLevel(
#ifndef NDEBUG
#endif
int tt = 0;
for (Range::iterator tit = output_ref_tris.begin();
tit != output_ref_tris.end(); tit++, tt++) {
int num_nodes;
CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
}
auto &data = fe_ptr->dataOnElement[
H1];
const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
int gg = 0;
for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
double *tri_coords = &ref_coords(tt, 0);
CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
auto det = t_normal.
l2();
for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
for (
int dd = 0;
dd != 2;
dd++) {
ref_gauss_pts(dd, gg) = shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
}
<< ref_gauss_pts(0, gg) << ", " << ref_gauss_pts(1, gg);
ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
}
}
int num_nodes;
CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
true);
std::bitset<numNodes> all_nodes;
for (auto nn = 0; nn != numNodes; ++nn) {
all_nodes.set(nn);
}
};
}
private:
boost::make_shared<VectorDouble>();
boost::make_shared<MatrixDouble>();
boost::make_shared<MatrixDouble>();
boost::make_shared<MatrixDouble>();
boost::make_shared<MatrixDouble>();
boost::make_shared<MatrixDouble>();
boost::make_shared<VectorDouble>();
boost::make_shared<MatrixDouble>();
};
};
#ifdef ADD_CONTACT
#ifdef ENABLE_PYTHON_BINDING
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string thermal_block_name, std::string thermoelastic_block_name,
std::string thermoplastic_block_name,
Pip &pip, std::string u,
std::string ep, std::string tau, std::string temperature) {
typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
createCommonThermoPlasticOps<DIM, I, DomainEleOp>(
m_field, block_name, thermal_block_name, thermoelastic_block_name,
thermoplastic_block_name, pip, u, ep, tau, temperature,
scale,
auto m_D_ptr = common_hencky_ptr->matDPtr;
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
temperature, common_thermoplastic_ptr->getTempPtr()));
"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
pip.push_back(
new
typename H::template OpCalculateHenckyThermoPlasticStress<DIM, I, 0>(
u, common_thermoplastic_ptr->getTempPtr(), common_hencky_ptr,
common_thermoelastic_ptr->getCoeffExpansionPtr(),
common_thermoelastic_ptr->getRefTempPtr()));
if (common_hencky_ptr) {
u, common_hencky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(
}
auto inelastic_heat_frac_ptr =
common_thermoplastic_ptr->getInelasticHeatFractionPtr();
return (*inelastic_heat_frac_ptr);
};
temperature, common_hencky_ptr->getMatHenckyStress(),
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
new
typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
DIM,
I>(ep, common_plastic_ptr, m_D_ptr));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
std::string thermal_block_name, std::string thermoelastic_block_name,
std::string thermoplastic_block_name,
Pip &pip, std::string u,
std::string ep, std::string tau, std::string temperature) {
using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
using OpKCauchy =
typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
createCommonThermoPlasticOps<DIM, I, DomainEleOp>(
m_field, block_name, thermal_block_name, thermoelastic_block_name,
thermoplastic_block_name, pip, u, ep, tau, temperature,
scale,
auto m_D_ptr = common_hencky_ptr->matDPtr;
ep, common_plastic_ptr->getPlasticStrainDotPtr()));
tau, common_plastic_ptr->getPlasticTauDotPtr()));
pip.push_back(new typename P::template OpCalculatePlasticity<DIM, I>(
u, common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
if (common_hencky_ptr) {
pip.push_back(new typename H::template OpHenckyTangent<DIM, I, 0>(
u, common_hencky_ptr, m_D_ptr));
pip.push_back(
new OpKPiola(u, u, common_hencky_ptr->getMatTangent()));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM,
I>(
u, ep, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
} else {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticInternalForceLhs_dEP<DIM,
I>(
u, ep, common_plastic_ptr, m_D_ptr));
}
if (common_hencky_ptr) {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_LogStrain_dU<DIM,
I>(
tau, u, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM,
I>(
ep, u, common_plastic_ptr, common_hencky_ptr, m_D_ptr));
} else {
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_dU<DIM,
I>(
tau, u, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_dU<DIM,
I>(
ep, u, common_plastic_ptr, m_D_ptr));
}
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_dEP<DIM,
I>(
ep, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculatePlasticFlowLhs_dTAU<DIM,
I>(
ep, tau, common_plastic_ptr, m_D_ptr));
pip.push_back(
ep, temperature, common_thermoplastic_ptr));
pip.push_back(
new typename P::template Assembly<A>::
template OpCalculateConstraintsLhs_dEP<DIM,
I>(
tau, ep, common_plastic_ptr, m_D_ptr));
pip.push_back(
new typename P::template Assembly<
A>::template OpCalculateConstraintsLhs_dTAU<I>(tau, tau,
common_plastic_ptr));
pip.push_back(new typename H::template OpCalculateHenckyThermalStressdT<
u, temperature, common_hencky_ptr,
common_thermoelastic_ptr->getCoeffExpansionPtr()));
auto inelastic_heat_frac_ptr =
common_thermoplastic_ptr->getInelasticHeatFractionPtr();
return (*inelastic_heat_frac_ptr);
};
temperature, temperature, common_hencky_ptr, common_plastic_ptr,
common_thermoelastic_ptr->getCoeffExpansionPtr()));
temperature, ep, common_hencky_ptr, common_plastic_ptr,
temperature, u, common_hencky_ptr, common_plastic_ptr,
tau, temperature, common_thermoplastic_ptr));
}
template <int DIM, IntegrationType I, typename DomainEleOp>
std::string thermal_block_name, std::string thermoelastic_block_name,
std::string thermoplastic_block_name,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string u, std::string ep, std::string tau, std::string temperature,
bool with_rates = true) {
auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
auto common_thermal_ptr =
boost::make_shared<ThermoElasticOps::BlockedThermalParameters>();
auto common_thermoelastic_ptr =
boost::make_shared<ThermoElasticOps::BlockedThermoElasticParameters>();
auto common_thermoplastic_ptr =
boost::make_shared<ThermoPlasticOps::ThermoPlasticBlockedParameters>();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
auto make_d_mat = []() {
};
common_plastic_ptr->mDPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
auto m_D_ptr = common_plastic_ptr->mDPtr;
m_field, plastic_block_name, pip, m_D_ptr,
common_plastic_ptr->getParamsPtr(),
scale, sev),
"add mat block plastic ops");
m_field, pip, thermal_block_name, common_thermal_ptr,
"add mat block thermal ops");
m_field, pip, thermoelastic_block_name,
"add mat block thermal ops");
m_field, pip, thermoplastic_block_name,
common_thermoplastic_ptr, common_thermal_ptr, sev,
"add mat block thermoplastic ops");
auto common_hencky_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
common_plastic_ptr->mDPtr = common_hencky_ptr->matDPtr;
tau, common_plastic_ptr->getPlasticTauPtr()));
ep, common_plastic_ptr->getPlasticStrainPtr()));
u, common_plastic_ptr->mGradPtr));
temperature, common_thermoplastic_ptr->getTempPtr()));
"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()));
common_plastic_ptr->mGradPtr = common_hencky_ptr->matGradPtr;
common_plastic_ptr->mDPtr = common_hencky_ptr->matDPtr;
common_hencky_ptr->matLogCPlastic =
common_plastic_ptr->getPlasticStrainPtr();
common_plastic_ptr->mStrainPtr = common_hencky_ptr->getMatLogC();
common_plastic_ptr->mStressPtr = common_hencky_ptr->getMatHenckyStress();
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
u, common_hencky_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(u, common_hencky_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
u, common_hencky_ptr));
pip.push_back(
new
typename H::template OpCalculateHenckyThermoPlasticStress<DIM, I, 0>(
u, common_thermoplastic_ptr->getTempPtr(), common_hencky_ptr,
common_thermoelastic_ptr->getCoeffExpansionPtr(),
common_thermoelastic_ptr->getRefTempPtr()));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
u, common_hencky_ptr));
pip.push_back(new typename P::template OpCalculatePlasticSurface<DIM, I>(
u, common_plastic_ptr));
pip.push_back(new typename P::template OpCalculatePlasticHardening<DIM, I>(
u, common_plastic_ptr, common_thermoplastic_ptr));
return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
common_thermal_ptr, common_thermoelastic_ptr,
common_thermoplastic_ptr);
}
};
MOFEM_LOG(
"THERMOPLASTICITY", Sev::inform) <<
"Run step pre proc";
auto &m_field = ptr->ExRawPtr->mField;
auto get_norm = [&](auto x) {
double nrm;
CHKERR VecNorm(x, NORM_2, &nrm);
return nrm;
};
auto save_range = [](moab::Interface &moab,
const std::string name,
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(),
1);
};
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(d_vector, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d_vector, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
Range new_elements, all_old_els, flipped_els, adj_edges;
auto improve_element_quality = [&]() {
double spatial_coords[9], material_coords[9];
std::multimap<double, EntityHandle> el_q_map, el_J_increased_map,
el_J_decreased_map;
auto get_element_quality_measures = [&]() {
CHKERR moab.get_entities_by_type(0, MBVERTEX, verts);
std::vector<double> coords(verts.size() * 3);
CHKERR moab.get_coords(verts, coords.data());
auto t_x = getFTensor1FromPtr<3>(coords.data());
auto save_tag = [&](boost::shared_ptr<FieldEntity> ent_ptr) {
auto field_data = ent_ptr->getEntFieldData();
t_u = {field_data[0], field_data[1], 0.0};
} else {
t_u = {field_data[0], field_data[1], field_data[2]};
}
++t_x;
};
&verts);
double def_coord[3] = {0.0, 0.0, 0.0};
CHKERR moab.tag_get_handle(
"SpatialCoord", 3, MB_TYPE_DOUBLE,
th_spatial_coords,
MB_TAG_DENSE | MB_TAG_CREAT, def_coord);
CHKERR moab.tag_set_data(th_spatial_coords, verts, coords.data());
int num_nodes;
Range::iterator nit = all_els.begin();
for (int gg = 0; nit != all_els.end(); nit++, gg++) {
true);
CHKERR moab.get_coords(conn, num_nodes, material_coords);
CHKERR moab.tag_get_data(th_spatial_coords, conn, num_nodes,
spatial_coords);
Tools::triArea(spatial_coords) / Tools::triArea(material_coords);
<<
"Element: " << *nit <<
" Jacobian: " <<
J;
double q = Tools::areaLengthQuality(spatial_coords);
if (!std::isnormal(q))
q = -2;
if (
J < 0.1 && q < qual_tol && q > 0) {
el_J_decreased_map.insert(std::pair<double, EntityHandle>(
J, *nit));
}
else if (
J > 3.0 && q < qual_tol && q > 0) {
el_J_increased_map.insert(std::pair<double, EntityHandle>(
J, *nit));
} else if (q < qual_tol && q > 0) {
el_q_map.insert(std::pair<double, EntityHandle>(q, *nit));
}
}
double min_q = 1;
all_els, min_q, th_spatial_coords);
<< "Old minimum element quality: " << min_q;
auto pair = el_q_map.begin();
<< "New minimum element quality: " << pair->first;
};
CHKERR get_element_quality_measures();
std::vector<EntityHandle> new_connectivity;
auto get_string_from_vector = [](auto vec) {
std::string output_string = "";
for (auto el : vec) {
output_string += std::to_string(el) + " ";
}
return output_string;
};
std::string improved_reasons = "";
auto perform_edge_flip = [&]() {
for (auto pair = el_q_map.begin(); pair != el_q_map.end(); ++pair) {
double quality = pair->first;
element.insert(pair->second);
if (!flipped_els.contains(element)) {
double edge_length;
double longest_edge_length = 0;
for (auto edge : edges) {
if (edge_length > longest_edge_length) {
longest_edge_length = edge_length;
longest_edge = edge;
}
}
<< "Edge flip longest edge length: " << longest_edge_length
<< " for edge: " << longest_edge;
auto get_skin = [&]() {
body_ents);
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
Range boundary_ents = get_skin();
<< "Boundary entities: " << boundary_ents;
if (boundary_ents.contains(
Range(longest_edge, longest_edge)))
continue;
Range flip_candidate_els;
&longest_edge, 1,
SPACE_DIM,
false, flip_candidate_els);
flip_candidate_els);
Range neighbouring_el = subtract(flip_candidate_els, element);
->filterEntitiesByRefLevel(
#ifndef NDEBUG
if (neighbouring_el.size() != 1) {
SETERRQ(
"Should be 1 neighbouring element to bad element for edge "
"flip. Instead, there are %d",
neighbouring_el.size());
}
#endif
<< "flip_candidate_els: " << flip_candidate_els;
<< "Neighbouring element: " << neighbouring_el;
if (flipped_els.contains(neighbouring_el))
continue;
std::vector<EntityHandle> neighbouring_nodes;
&neighbouring_el.front(), 1, neighbouring_nodes, true);
std::vector<EntityHandle> element_nodes;
element_nodes, true);
CHKERR moab.tag_get_data(th_spatial_coords,
&neighbouring_nodes.front(), 3,
spatial_coords);
double neighbour_qual = Tools::areaLengthQuality(spatial_coords);
if (!std::isnormal(neighbour_qual))
neighbour_qual = -2;
<< "Element nodes before swap: "
<< get_string_from_vector(element_nodes);
<< "Neighbouring nodes before swap: "
<< get_string_from_vector(neighbouring_nodes);
std::vector<EntityHandle> reversed_neighbouring_nodes =
neighbouring_nodes;
std::reverse(reversed_neighbouring_nodes.begin(),
reversed_neighbouring_nodes.end());
int num_matches = 0;
std::vector<bool> mismatch_mask(element_nodes.size());
int loop_counter = 0;
while (num_matches != 2) {
std::rotate(reversed_neighbouring_nodes.begin(),
reversed_neighbouring_nodes.begin() + 1,
reversed_neighbouring_nodes.end());
std::transform(element_nodes.begin(), element_nodes.end(),
reversed_neighbouring_nodes.begin(),
mismatch_mask.begin(),
std::equal_to<EntityHandle>());
num_matches =
std::count(mismatch_mask.begin(), mismatch_mask.end(), true);
++loop_counter;
if (loop_counter > 3) {
"Not found matching nodes for edge flipping");
}
}
std::vector<EntityHandle> matched_elements(element_nodes.size());
std::transform(element_nodes.begin(), element_nodes.end(),
mismatch_mask.begin(), matched_elements.begin(),
return match ? el
: -1;
});
matched_elements.erase(std::remove(matched_elements.begin(),
matched_elements.end(), -1),
matched_elements.end());
std::vector<EntityHandle> mismatched_elements(element_nodes.size()),
neighbouring_mismatched_elements(neighbouring_nodes.size());
std::transform(element_nodes.begin(), element_nodes.end(),
mismatch_mask.begin(), mismatched_elements.begin(),
return match ? -1
: el;
});
std::transform(
reversed_neighbouring_nodes.begin(),
reversed_neighbouring_nodes.end(), mismatch_mask.begin(),
neighbouring_mismatched_elements.begin(),
return match ? -1 : el;
});
mismatched_elements.erase(std::remove(mismatched_elements.begin(),
mismatched_elements.end(),
-1),
mismatched_elements.end());
neighbouring_mismatched_elements.erase(
std::remove(neighbouring_mismatched_elements.begin(),
neighbouring_mismatched_elements.end(), -1),
neighbouring_mismatched_elements.end());
mismatched_elements.insert(mismatched_elements.end(),
neighbouring_mismatched_elements.begin(),
neighbouring_mismatched_elements.end());
<< "Reversed neighbouring nodes: "
<< get_string_from_vector(reversed_neighbouring_nodes);
<< "mismatch mask: " << get_string_from_vector(mismatch_mask);
<< "Old nodes are: "
<< get_string_from_vector(matched_elements);
<< "New nodes are: "
<< get_string_from_vector(mismatched_elements);
auto replace_correct_nodes =
[](std::vector<EntityHandle> &ABC,
std::vector<EntityHandle> &DBA,
const std::vector<EntityHandle> &AB,
const std::vector<EntityHandle> &CD) {
std::vector<std::vector<EntityHandle>> results;
for (
int i = 0;
i < 2; ++
i) {
for (
int j = 0;
j < 2; ++
j) {
if (std::find(ABC.begin(), ABC.end(), CD[
j]) ==
ABC.end()) {
std::vector<EntityHandle> tmp = ABC;
auto it = std::find(tmp.begin(), tmp.end(), AB[
i]);
if (it != tmp.end()) {
results.push_back(tmp);
}
}
}
}
if (results.size() != 2) {
SETERRQ(
"Failed to find two valid vertex replacements for edge "
"flip");
}
ABC = results[0];
DBA = results[1];
};
CHKERR replace_correct_nodes(element_nodes, neighbouring_nodes,
matched_elements, mismatched_elements);
<< "Element nodes after swap: "
<< get_string_from_vector(element_nodes);
<< "Neighbouring nodes after swap: "
<< get_string_from_vector(neighbouring_nodes);
CHKERR moab.tag_get_data(th_spatial_coords, &element_nodes.front(),
3, spatial_coords);
double new_qual = Tools::areaLengthQuality(spatial_coords);
if (!std::isnormal(new_qual))
new_qual = -2;
#ifndef NDEBUG
auto check_normal_direction = [&]() {
auto t_correct_normal =
auto [new_area, t_new_normal] =
Tools::triAreaAndNormal(spatial_coords);
t_diff(
i) = t_new_normal(
i) - t_correct_normal(
i);
if (t_diff(
i) * t_diff(
i) > 1e-6) {
SETERRQ(
"Direction of element to be created is wrong orientation");
}
};
CHKERR check_normal_direction();
#endif
CHKERR moab.tag_get_data(th_spatial_coords,
&neighbouring_nodes.front(), 3,
spatial_coords);
double new_neighbour_qual =
Tools::areaLengthQuality(spatial_coords);
if (!std::isnormal(new_neighbour_qual))
new_neighbour_qual = -2;
#ifndef NDEBUG
CHKERR check_normal_direction();
#endif
if (std::min(new_qual, new_neighbour_qual) >
std::min(quality, neighbour_qual)) {
<< "Element quality improved from " << quality << " and "
<< neighbour_qual << " to " << new_qual << " and "
<< new_neighbour_qual << " for elements" << element << " and "
<< neighbouring_el;
flipped_els.merge(flip_candidate_els);
new_connectivity.insert(new_connectivity.end(),
element_nodes.begin(),
element_nodes.end());
new_connectivity.insert(new_connectivity.end(),
neighbouring_nodes.begin(),
neighbouring_nodes.end());
}
}
}
improved_reasons += ", Edge flipping";
};
<< "Flipped elements: " << flipped_els;
<< "New connectivity: " << get_string_from_vector(new_connectivity);
CHKERR moab.get_entities_by_type(0, MBTRI, all_old_els,
true);
<< "Elements added to old_mesh: " << all_old_els;
old_mesh_bitreflevel,
BitRefLevel().set(), old_mesh_range, 1);
<< "Old mesh range: " << old_mesh_range;
#ifndef NDEBUG
(boost::format("old_mesh.vtk")).str(), old_mesh_range);
#endif
Range remaining_els_after_flip = subtract(all_old_els, flipped_els);
<< "Non-flipped elements: " << remaining_els_after_flip;
ReadUtilIface *read_util;
const int num_ele = flipped_els.size();
int num_nod_per_ele;
EntityType ent_type;
num_nod_per_ele = 3;
ent_type = MBTRI;
} else {
num_nod_per_ele = 4;
ent_type = MBTET;
}
if (flipped_els.size() > 0) {
auto new_conn = new_connectivity.begin();
for (auto e = 0; e != num_ele; ++e) {
for (
int n = 0;
n < num_nod_per_ele; ++
n) {
++new_conn;
}
if (adj_ele.size()) {
if (adj_ele.size() != 1) {
"Element duplication");
} else {
new_ele = adj_ele.front();
}
} else {
num_nod_per_ele, new_ele);
}
new_elements.insert(new_ele);
}
<< "New elements: " << new_elements;
<< "New elements from before: " << remaining_els_after_flip;
Range new_elements_nodes_vertices;
new_elements_nodes_vertices.merge(new_elements);
new_elements_nodes_vertices.merge(remaining_els_after_flip);
CHKERR moab.get_adjacencies(new_elements_nodes_vertices, 1,
true,
new_edges_range, moab::Interface::UNION);
Range new_vertices_range;
CHKERR moab.get_adjacencies(new_elements_nodes_vertices, 0,
true,
new_vertices_range, moab::Interface::UNION);
new_elements_nodes_vertices.merge(new_edges_range);
new_elements_nodes_vertices.merge(new_vertices_range);
<< "New mesh after flip: " << new_elements_nodes_vertices;
new_elements_nodes_vertices, new_mesh_bitreflevel, false);
#ifndef NDEBUG
new_mesh_bitreflevel,
BitRefLevel().set(), new_mesh_range, 1);
<< "New mesh range: " << new_mesh_range;
(boost::format("new_mesh.vtk")).str(),
new_mesh_range);
#endif
improved_reasons.erase(0, 2);
<< "Improved mesh quality by: " << improved_reasons;
};
};
auto solve_equil_sub_problem = [&]() {
if (flipped_els.size() == 0) {
}
auto U_ptr = boost::make_shared<MatrixDouble>();
auto EP_ptr = boost::make_shared<MatrixDouble>();
auto TAU_ptr = boost::make_shared<VectorDouble>();
auto T_ptr = boost::make_shared<VectorDouble>();
auto post_proc = [&](auto dm) {
auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
post_proc_fe->getOpPtrVector(), {H1, L2, HDIV}, "GEOMETRY");
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
EP_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{
{"TAU", TAU_ptr},
{"T", T_ptr},
},
{{"U", U_ptr}},
{},
{{"EP", EP_ptr}}
)
);
CHKERR post_proc_fe->writeFile(
"out_equilibrate.h5m");
};
auto solve_equilibrate_solution = [&]() {
auto set_domain_rhs = [&](auto &pip) {
"GEOMETRY");
typename B::template OpGradTimesSymTensor<1,
SPACE_DIM,
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
ptr->ExRawPtr
m_field, "MAT_PLASTIC", "MAT_THERMAL",
"MAT_THERMOELASTIC", "MAT_THERMOPLASTIC", pip, "U", "EP",
Sev::inform);
auto m_D_ptr = common_plastic_ptr->mDPtr;
"T", common_thermoplastic_ptr->getTempPtr()));
pip.push_back(
new
typename P::template OpCalculatePlasticityWithoutRates<
SPACE_DIM,
"U", common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
if (common_hencky_ptr) {
"U", common_hencky_ptr->getMatFirstPiolaStress()));
} else {
pip.push_back(
}
};
auto set_domain_lhs = [&](auto &pip) {
typename B::template OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
0>;
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
ptr->ExRawPtr
m_field, "MAT_PLASTIC", "MAT_THERMAL",
"MAT_THERMOELASTIC", "MAT_THERMOPLASTIC", pip, "U", "EP",
Sev::inform);
auto m_D_ptr = common_plastic_ptr->mDPtr;
pip.push_back(
new
typename P::template OpCalculatePlasticityWithoutRates<
SPACE_DIM,
"U", common_plastic_ptr, m_D_ptr, common_thermoplastic_ptr));
if (common_hencky_ptr) {
pip.push_back(
new typename H::template OpHenckyTangent<SPACE_DIM, IT, 0>(
"U", common_hencky_ptr, m_D_ptr));
pip.push_back(
new OpKPiola(
"U",
"U", common_hencky_ptr->getMatTangent()));
} else {
pip.push_back(
new OpKCauchy(
"U",
"U", m_D_ptr));
}
};
auto dm_sub =
createDM(m_field.get_comm(),
"DMMOFEM");
for (
auto f : {
"U",
"EP",
"TAU",
"T"}) {
}
return dm_sub;
};
auto sub_dm = create_sub_dm(
simple->getDM());
auto fe_rhs = boost::make_shared<DomainEle>(m_field);
auto fe_lhs = boost::make_shared<DomainEle>(m_field);
fe_rhs->getRuleHook = vol_rule;
fe_lhs->getRuleHook = vol_rule;
CHKERR set_domain_rhs(fe_rhs->getOpPtrVector());
CHKERR set_domain_lhs(fe_lhs->getOpPtrVector());
auto bc_mng = m_field.getInterface<
BcManager>();
"EQUIL_DM", "U");
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map)
MOFEM_LOG(
"PLASTICITY",
Sev::verbose) <<
"Marker " << bc.first;
auto time_scale = boost::make_shared<TimeScale>(
"", false, [](const double) { return 1; });
auto def_time_scale = [time_scale](const double time) {
return (!time_scale->argFileScale) ? time : 1;
};
auto def_file_name = [time_scale](const std::string &&name) {
return (!time_scale->argFileScale) ? name : "";
};
auto time_displacement_scaling = boost::make_shared<TimeScale>(
def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
PetscReal current_time;
TSGetTime(ts, ¤t_time);
pre_proc_ptr->ts_t = current_time;
auto get_bc_hook_rhs = [&]() {
ptr->ExRawPtr->mField, pre_proc_ptr,
{time_scale, time_displacement_scaling});
return hook;
};
auto get_bc_hook_lhs = [&]() {
ptr->ExRawPtr->mField, pre_proc_ptr,
{time_scale, time_displacement_scaling});
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
pre_proc_ptr->preProcessHook = get_bc_hook_lhs();
auto get_post_proc_hook_rhs = [&]() {
ptr->ExRawPtr->mField, post_proc_rhs_ptr, nullptr,
Sev::verbose)();
ptr->ExRawPtr->mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [&]() {
ptr->ExRawPtr->mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
snes_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_ptr);
snes_ctx_ptr->getPreProcSetOperators().push_front(pre_proc_ptr);
snes_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs_ptr);
snes_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs_ptr);
auto null_fe = boost::shared_ptr<FEMethod>();
null_fe, null_fe);
null_fe, null_fe);
CHKERR SNESSetDM(snes, sub_dm);
CHKERR SNESSetFromOptions(snes);
SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR SNESSolve(snes, PETSC_NULLPTR,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
SCATTER_REVERSE);
};
CHKERR solve_equilibrate_solution();
<< "Equilibrated solution after mapping";
};
CHKERR improve_element_quality();
if (flipped_els.size() > 0) {
CHKERR TSGetApplicationContext(ts, (
void **)&ctx);
if (!ctx)
}
CHKERR PetscBarrier((PetscObject)ts);
}
}
auto &m_field = ptr->ExRawPtr->mField;
<< "PostProc done";
}
return 0;
}
}
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
std::string thermoplastic_block_name,
boost::shared_ptr<ThermoPlasticOps::ThermoPlasticBlockedParameters>
blockedParamsPtr,
boost::shared_ptr<ThermoElasticOps::BlockedThermalParameters>
&blockedThermalParamsPtr,
OpMatThermoPlasticBlocks(
boost::shared_ptr<double> omega_0_ptr,
boost::shared_ptr<double> omega_h_ptr,
boost::shared_ptr<double> inelastic_heat_fraction_ptr,
boost::shared_ptr<double> temp_0_ptr,
boost::shared_ptr<double> heat_capacity,
boost::shared_ptr<double> heat_conductivity_scaling,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
omegaHPtr(omega_h_ptr),
inelasticHeatFractionPtr(inelastic_heat_fraction_ptr),
heatCapacityPtr(heat_capacity),
heatConductivityScalingPtr(heat_conductivity_scaling),
extractThermoPlasticBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
auto scale_param = [this](auto parameter, double scaling) {
return parameter * scaling;
};
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*omega0Ptr = b.omega_0;
*omegaHPtr = b.omega_h;
*inelasticHeatFractionPtr = scale_param(
b.inelastic_heat_fraction,
inelasticHeatFractionScaling(
*heatConductivityPtr / (*heatConductivityScalingPtr),
*heatCapacityPtr / (*heatCapacityScalingPtr)));
*temp0Ptr = b.temp_0;
}
}
*inelasticHeatFractionPtr =
inelasticHeatFractionScaling(
*heatConductivityPtr / (*heatConductivityScalingPtr),
*heatCapacityPtr / (*heatCapacityScalingPtr)));
}
private:
boost::shared_ptr<double> heatConductivityPtr;
boost::shared_ptr<double> heatCapacityPtr;
boost::shared_ptr<double> heatConductivityScalingPtr;
boost::shared_ptr<double> heatCapacityScalingPtr;
};
std::vector<BlockData> blockData;
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (
auto m : meshset_vec_ptr) {
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 3) {
"Expected that block has four attributes");
}
auto get_block_ents = [&]() {
m_field.
get_moab().get_entities_by_handle(
m->meshset, ents,
true);
return ents;
};
<<
m->getName() <<
": omega_0 = " << block_data[0]
<< " omega_h = " << block_data[1]
<< " inelastic_heat_fraction = " << block_data[2] << " temp_0 "
<< block_data[3];
blockData.push_back({block_data[0], block_data[1], block_data[2],
block_data[3], get_block_ents()});
}
}
boost::shared_ptr<double> omega0Ptr;
boost::shared_ptr<double> omegaHPtr;
boost::shared_ptr<double> inelasticHeatFractionPtr;
boost::shared_ptr<double> temp0Ptr;
};
pipeline.push_back(new OpMatThermoPlasticBlocks(
blockedParamsPtr->getOmega0Ptr(), blockedParamsPtr->getOmegaHPtr(),
blockedParamsPtr->getInelasticHeatFractionPtr(),
blockedParamsPtr->getTemp0Ptr(),
blockedThermalParamsPtr->getHeatConductivityPtr(),
blockedThermalParamsPtr->getHeatCapacityPtr(),
blockedThermalParamsPtr->getHeatConductivityScalingPtr(),
blockedThermalParamsPtr->getHeatCapacityScalingPtr(),
inelastic_heat_fraction_scaling_func, m_field, sev,
(boost::format("%s(.*)") % thermoplastic_block_name).str()
))
));
}
}
true);
auto get_ents_by_dim = [&](const auto dim) {
return domain_ents;
} else {
if (dim == 0)
else
return ents;
}
};
auto get_base = [&]() {
auto domain_ents = get_ents_by_dim(
SPACE_DIM);
if (domain_ents.empty())
switch (type) {
case MBQUAD:
case MBHEX:
case MBTRI:
case MBTET:
default:
}
};
const auto base = get_base();
#ifdef ADD_CONTACT
auto get_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_blocks = [&](auto skin) {
bool is_contact_block = false;
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true;
<<
"Find contact block set: " <<
m->getName();
auto meshset =
m->getMeshset();
Range contact_meshset_range;
meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
contact_meshset_range);
contact_range.merge(contact_meshset_range);
}
if (is_contact_block) {
<< "Nb entities in contact surface: " << contact_range.size();
skin = intersect(skin, contact_range);
}
return skin;
};
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
#endif
auto project_ho_geometry = [&]() {
};
PetscBool project_geometry = PETSC_TRUE;
&project_geometry, PETSC_NULLPTR);
if (project_geometry) {
}
}
auto get_command_line_parameters = [&]() {
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
<<
"Minumum quality for edge flipping: " <<
qual_tol;
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
(PetscInt *)&
ic_type, PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
PetscBool tau_order_is_set;
&tau_order_is_set);
PetscBool ep_order_is_set;
&ep_order_is_set);
PetscBool flux_order_is_set;
&flux_order_is_set);
PetscBool T_order_is_set;
&T_order_is_set);
PETSC_NULLPTR);
PETSC_NULLPTR);
PETSC_NULLPTR);
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Hardening " <<
H;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Viscous hardening " <<
visH;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation yield stress " <<
Qinf;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Saturation exponent " <<
b_iso;
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Kinematic hardening " <<
C1_k;
if (tau_order_is_set == PETSC_FALSE)
if (ep_order_is_set == PETSC_FALSE)
if (flux_order_is_set == PETSC_FALSE)
if (T_order_is_set == PETSC_FALSE)
MOFEM_LOG(
"PLASTICITY", Sev::inform) <<
"Approximation order " <<
order;
<<
"Ep approximation order " <<
ep_order;
PetscBool is_scale = PETSC_TRUE;
PETSC_NULLPTR);
MOFEM_LOG(
"THERMOPLASTICITY", Sev::inform) <<
"Scaling used? " << is_scale;
switch (is_scale) {
case PETSC_TRUE:
<< "Scaling does not yet work with radiation and convection BCs";
if (heat_cap != 0.0 && thermal_cond != 0.0)
return 1. / (thermal_cond * heat_cap);
else
return 1. / thermal_cond;
};
if (heat_cap != 0.0)
return 1. / (thermal_cond * heat_cap);
else
return 1.0;
};
[&](
double young_modulus,
double thermal_cond,
double heat_cap) {
if (heat_cap != 0.0)
else
};
break;
default:
return 1.0;
};
return 1.0;
};
double thermal_cond,
double heat_cap) { return 1.0; };
break;
}
<< "Thermal conductivity scale "
<< "Thermal capacity scale "
<< "Inelastic heat fraction scale "
#ifdef ADD_CONTACT
#endif
};
CHKERR get_command_line_parameters();
#ifdef ADD_CONTACT
#ifdef ENABLE_PYTHON_BINDING
sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdfPythonPtr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif
}
auto T_ptr = boost::make_shared<VectorDouble>();
auto post_proc = [&](auto dm) {
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
post_proc_fe->getOpPtrVector().push_back(
auto TAU_ptr = boost::make_shared<VectorDouble>();
auto EP_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"T", T_ptr}, {"TAU", TAU_ptr}},
{},
{},
{{"EP", EP_ptr}}
)
);
} else {
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{{"T", T_ptr}},
{},
{},
{}
)
);
}
CHKERR post_proc_fe->writeFile(
"out_init.h5m");
};
auto solve_init = [&]() {
auto set_domain_rhs = [&](auto &pip) {
"GEOMETRY");
"T",
nullptr, T_ptr,
nullptr, boost::make_shared<double>(
init_temp),
auto TAU_ptr = boost::make_shared<VectorDouble>();
auto EP_ptr = boost::make_shared<MatrixDouble>();
auto min_tau = boost::make_shared<double>(1.0);
auto max_tau = boost::make_shared<double>(2.0);
nullptr, min_tau, max_tau));
"EP", EP_ptr));
auto min_EP = boost::make_shared<double>(0.0);
auto max_EP = boost::make_shared<double>(0.01);
"EP", nullptr, EP_ptr, nullptr, min_EP, max_EP));
}
};
auto set_domain_lhs = [&](auto &pip) {
"GEOMETRY");
pip.push_back(new OpLhsScalarLeastSquaresProj("T", "T"));
pip.push_back(new OpLhsScalarLeastSquaresProj("TAU", "TAU"));
pip.push_back(
"EP", "EP"));
}
};
for (auto f : {"T"}) {
}
for (auto f : {"TAU", "EP"}) {
}
}
return dm_sub;
};
auto fe_rhs = boost::make_shared<DomainEle>(mField);
auto fe_lhs = boost::make_shared<DomainEle>(mField);
fe_rhs->getRuleHook = vol_rule;
fe_lhs->getRuleHook = vol_rule;
CHKERR set_domain_rhs(fe_rhs->getOpPtrVector());
CHKERR set_domain_lhs(fe_lhs->getOpPtrVector());
auto sub_dm = create_sub_dm(
simple->getDM());
auto null_fe = boost::shared_ptr<FEMethod>();
fe_lhs, null_fe, null_fe);
null_fe, null_fe);
CHKERR KSPSetFromOptions(ksp);
CHKERR KSPSolve(ksp, PETSC_NULLPTR,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
MOFEM_LOG(
"THERMAL", Sev::inform) <<
"Set thermoelastic initial conditions";
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0, true,
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1, true,
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2, true,
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3, true,
#ifdef ADD_CONTACT
for (auto b : {"FIX_X", "REMOVE_X"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
for (auto b : {"FIX_Y", "REMOVE_Y"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
for (auto b : {"FIX_Z", "REMOVE_Z"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
for (auto b : {"FIX_ALL", "REMOVE_ALL"})
CHKERR bc_mng->removeBlockDOFsOnEntities(
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3, false,
#endif
simple->getProblemName(),
"U");
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map)
MOFEM_LOG(
"PLASTICITY",
Sev::verbose) <<
"Marker " << bc.first;
}
auto get_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_flux_blocks = [&](auto skin, bool temp_bc = false) {
auto remove_cubit_blocks = [&](
auto c) {
) {
skin = subtract(skin, ents);
}
};
auto remove_named_blocks = [&](
auto n) {
(boost::format(
"%s(.*)") %
n).str()
))
) {
skin = subtract(skin, ents);
}
};
if (!temp_bc) {
"remove_cubit_blocks");
"remove_named_blocks");
}
"remove_cubit_blocks");
return skin;
};
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto remove_temp_bc_ents =
remove_flux_ents);
remove_temp_bc_ents);
MOFEM_LOG(
"SYNC", Sev::noisy) << remove_flux_ents << endl;
MOFEM_LOG(
"SYNC", Sev::noisy) << remove_temp_bc_ents << endl;
#ifndef NDEBUG
remove_flux_ents);
#endif
simple->getProblemName(),
"FLUX", remove_flux_ents);
simple->getProblemName(),
"TBC", remove_temp_bc_ents);
} else {
->removeDofsOnEntitiesNotDistributed(
simple->getProblemName(),
"FLUX",
remove_flux_ents);
->removeDofsOnEntitiesNotDistributed(
simple->getProblemName(),
"TBC",
remove_temp_bc_ents);
}
simple->getProblemName(),
"FLUX",
false);
}
auto boundary_marker =
bc_mng->getMergedBlocksMarker(vector<string>{"HEATFLUX"});
auto integration_rule_bc = [](
int,
int,
int ao) {
return 2 * ao; };
CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto thermal_block_params =
boost::make_shared<ThermoElasticOps::BlockedThermalParameters>();
auto heat_conductivity_ptr = thermal_block_params->getHeatConductivityPtr();
auto heat_capacity_ptr = thermal_block_params->getHeatCapacityPtr();
auto thermoelastic_block_params =
boost::make_shared<ThermoElasticOps::BlockedThermoElasticParameters>();
auto coeff_expansion_ptr = thermoelastic_block_params->getCoeffExpansionPtr();
auto ref_temp_ptr = thermoelastic_block_params->getRefTempPtr();
auto time_scale =
boost::make_shared<TimeScale>("", false, [](const double) { return 1; });
auto def_time_scale = [time_scale](const double time) {
return (!time_scale->argFileScale) ? time : 1;
};
auto def_file_name = [time_scale](const std::string &&name) {
return (!time_scale->argFileScale) ? name : "";
};
auto time_bodyforce_scaling = boost::make_shared<TimeScale>(
def_file_name("bodyforce_scale.txt"), false, def_time_scale);
auto time_heatsource_scaling = boost::make_shared<TimeScale>(
def_file_name("heatsource_scale.txt"), false, def_time_scale);
auto time_temperature_scaling = boost::make_shared<TimeScale>(
def_file_name("temperature_bc_scale.txt"), false, def_time_scale);
auto time_displacement_scaling = boost::make_shared<TimeScale>(
def_file_name("displacement_bc_scale.txt"), false, def_time_scale);
auto time_flux_scaling = boost::make_shared<TimeScale>(
def_file_name("flux_bc_scale.txt"), false, def_time_scale);
auto time_force_scaling = boost::make_shared<TimeScale>(
def_file_name("force_bc_scale.txt"), false, def_time_scale);
auto add_boundary_ops_lhs = [&](auto &pip) {
pip.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
"GEOMETRY");
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip,
mField,
"U", Sev::inform);
#ifdef ADD_CONTACT
ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
pip, "SIGMA", "U");
ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
mField, pip,
simple->getDomainFEName(),
"SIGMA",
"U",
"GEOMETRY",
vol_rule);
#endif
using T =
using OpConvectionLhsBC =
T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
using OpRadiationLhsBC =
T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
auto temp_bc_ptr = boost::make_shared<VectorDouble>();
T::AddFluxToPipeline<OpConvectionLhsBC>::add(pip,
mField,
"TBC",
"TBC",
"CONVECTION", Sev::verbose);
T::AddFluxToPipeline<OpRadiationLhsBC>::add(
pip,
mField,
"TBC",
"TBC", temp_bc_ptr,
"RADIATION", Sev::verbose);
struct OpTBCTimesNormalFlux
OpTBCTimesNormalFlux(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<Range> ents_ptr = nullptr)
:
OP(row_field_name, col_field_name,
OP::OPROWCOL, ents_ptr) {
this->sYmm = false;
this->assembleTranspose = true;
this->onlyTranspose = false;
}
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
auto a_mat_ptr = &*OP::locMat.data().begin();
int rr = 0;
for (; rr != OP::nbRows; rr++) {
const double alpha = t_w * t_row_base;
for (int cc = 0; cc != OP::nbCols; cc++) {
*a_mat_ptr += alpha * (t_col_base(
i) * t_normal(
i));
++t_col_base;
++a_mat_ptr;
}
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_normal;
++t_w;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locMat /= 2;
}
}
};
pip.push_back(new OpTBCTimesNormalFlux("TBC", "FLUX"));
};
auto add_boundary_ops_rhs = [&](auto &pip) {
"GEOMETRY");
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
pip,
mField,
"U", {time_scale, time_force_scaling},
"FORCE",
Sev::inform);
BoundaryNaturalBC::OpFlux<NaturalTemperatureMeshsets, 3, SPACE_DIM>;
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpTemperatureBC>::add(
pip,
mField,
"FLUX", {time_scale, time_temperature_scaling},
"TEMPERATURE", Sev::inform);
using OpFluxBC =
BoundaryNaturalBC::OpFlux<NaturalMeshsetType<HEATFLUXSET>, 1, 1>;
CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpFluxBC>::add(
pip,
mField,
"TBC", {time_scale, time_flux_scaling},
"FLUX",
Sev::inform);
using OpConvectionRhsBC =
T::OpFlux<ThermoElasticOps::ConvectionBcType<BLOCKSET>, 1, 1>;
using OpRadiationRhsBC =
T::OpFlux<ThermoElasticOps::RadiationBcType<BLOCKSET>, 1, 1>;
auto temp_bc_ptr = boost::make_shared<VectorDouble>();
T::AddFluxToPipeline<OpConvectionRhsBC>::add(
pip,
mField,
"TBC", temp_bc_ptr,
"CONVECTION", Sev::inform);
T::AddFluxToPipeline<OpRadiationRhsBC>::add(pip,
mField,
"TBC", temp_bc_ptr,
"RADIATION", Sev::inform);
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
struct OpTBCTimesNormalFlux
OpTBCTimesNormalFlux(
const std::string
field_name,
boost::shared_ptr<MatrixDouble> vec,
boost::shared_ptr<Range> ents_ptr = nullptr)
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_vec = getFTensor1FromMat<SPACE_DIM, 1>(*sourceVec);
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * (t_vec(
i) * t_normal(
i));
int rr = 0;
for (; rr != OP::nbRows; ++rr) {
OP::locF[rr] += alpha * t_row_base;
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_vec;
++t_normal;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locF /= 2;
}
}
protected:
boost::shared_ptr<MatrixDouble> sourceVec;
};
pip.push_back(new OpTBCTimesNormalFlux("TBC", mat_flux_ptr));
struct OpBaseNormalTimesTBC
OpBaseNormalTimesTBC(
const std::string
field_name,
boost::shared_ptr<VectorDouble> vec,
boost::shared_ptr<Range> ents_ptr = nullptr)
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_normal = OP::getFTensor1NormalsAtGaussPts();
auto t_vec = getFTensor0FromVec(*sourceVec);
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
const double alpha = t_w * t_vec;
int rr = 0;
for (; rr != OP::nbRows; ++rr) {
OP::locF[rr] += alpha * (t_row_base(
i) * t_normal(
i));
++t_row_base;
}
for (; rr < OP::nbRowBaseFunctions; ++rr)
++t_row_base;
++t_w;
++t_vec;
++t_normal;
}
EntityType fe_type = OP::getNumeredEntFiniteElementPtr()->getEntType();
if (fe_type == MBTRI) {
OP::locF /= 2;
}
}
protected:
boost::shared_ptr<VectorDouble> sourceVec;
};
pip.push_back(new OpBaseNormalTimesTBC("FLUX", temp_bc_ptr));
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
pip, "SIGMA", "U");
#endif
};
auto add_domain_ops_lhs = [&](auto &pip) {
pip.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
mField, pip,
"MAT_THERMAL", thermal_block_params,
"GEOMETRY");
auto get_inertia_and_mass_damping = [
this](
const double,
const double,
return (
rho /
scale) * fe_domain_lhs->ts_aa +
};
pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
}
CHKERR opThermoPlasticFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC",
"MAT_THERMAL",
"MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T");
auto hencky_common_data_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
pip.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
pip.push_back(
new OpHdivHdiv(
"FLUX",
"FLUX", resistance, mat_grad_ptr));
pip.push_back(
new OpHdivT(
"FLUX",
"T", []()
constexpr {
return -1; },
true));
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
new OpHdivU(
"FLUX",
"U", mat_flux_ptr, resistance, mat_grad_ptr));
auto op_capacity =
new OpCapacity(
"T",
"T", capacity);
op_capacity->feScalingFun = [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a;
};
pip.push_back(op_capacity);
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
CHKERR DomainNaturalBCLhs::AddFluxToPipeline<OpSetTemperatureLhs>::add(
pip,
mField,
"T", vec_temp_ptr,
"SETTEMP", Sev::verbose);
};
auto add_domain_ops_rhs = [&](auto &pip) {
pip.push_back(
new OpSetBc(
"FLUX",
true, boundary_marker));
mField, pip,
"MAT_THERMAL", thermal_block_params,
"GEOMETRY");
auto hencky_common_data_ptr =
HenckyOps::commonDataFactory<SPACE_DIM, IT, DomainEleOp>(
auto mat_D_ptr = hencky_common_data_ptr->matDPtr;
auto mat_grad_ptr = hencky_common_data_ptr->matGradPtr;
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto vec_temp_dot_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_div_ptr = boost::make_shared<VectorDouble>();
"FLUX", vec_temp_div_ptr));
pip.push_back(
auto resistance = [heat_conductivity_ptr](
const double,
const double,
return (1. / (*heat_conductivity_ptr));
};
auto capacity = [heat_capacity_ptr](
const double,
const double,
return -(*heat_capacity_ptr);
};
return -1.;
};
pip.push_back(
new OpHdivFlux(
"FLUX", mat_flux_ptr, resistance, mat_grad_ptr));
pip.push_back(
new OpHDivTemp(
"FLUX", vec_temp_ptr, unity));
pip.push_back(
new OpBaseDotT(
"T", vec_temp_dot_ptr, capacity));
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
{boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
Sev::inform);
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(
new OpInertiaForce(
"U", mat_acceleration, [](
double,
double,
double) {
}));
auto mat_velocity = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
}));
}
}
CHKERR opThermoPlasticFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC",
"MAT_THERMAL",
"MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T");
#ifdef ADD_CONTACT
CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
pip, "SIGMA", "U");
#endif
};
auto create_reaction_pipeline = [&](auto &pip) {
"GEOMETRY");
CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
mField,
"MAT_PLASTIC", pip,
"U",
"EP",
"TAU");
};
auto get_bc_hook_rhs = [&]() {
mField, pip_mng->getDomainRhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
auto get_bc_hook_lhs = [&]() {
mField, pip_mng->getDomainLhsFE(),
{time_scale, time_displacement_scaling});
return hook;
};
pip_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
pip_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
}
);
protected:
};
boost::shared_ptr<MoFEM::TsCtx>
dmts_ctx;
};
static std::unordered_map<TS, MyTsCtx *>
ts_to_ctx;
}
void *ctx) {
}
PetscErrorCode
MyTSResizeSetup(TS ts, PetscInt nstep, PetscReal time, Vec sol,
PetscBool *resize, void *ctx) {
PetscFunctionBegin;
PetscFunctionReturn(0);
}
Vec ts_vecsout[], void *ctx) {
PetscFunctionBegin;
MOFEM_LOG(
"REMESHING", Sev::verbose) <<
"number of vectors to map = " << nv;
for (PetscInt
i = 0;
i < nv; ++
i) {
double nrm;
CHKERR VecNorm(ts_vecsin[
i], NORM_2, &nrm);
<<
"Before remeshing: ts_vecsin[" <<
i <<
"] norm = " << nrm;
}
CHKERR DMSetFromOptions(intermediate_dm);
CHKERR DMSetUp(intermediate_dm);
Vec vec_in[nv], vec_out[nv];
for (PetscInt
i = 0;
i < nv; ++
i) {
CHKERR DMCreateGlobalVector(intermediate_dm, &vec_in[
i]);
CHKERR VecDuplicate(vec_in[
i], &vec_out[
i]);
}
VecScatter scatter_to_intermediate;
for (PetscInt
i = 0;
i < nv; ++
i) {
CHKERR scatter_mng->vecScatterCreate(
CHKERR VecScatterBegin(scatter_to_intermediate, ts_vecsin[
i], vec_in[
i],
INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterEnd(scatter_to_intermediate, ts_vecsin[
i], vec_in[
i],
INSERT_VALUES, SCATTER_FORWARD);
}
CHKERR VecScatterDestroy(&scatter_to_intermediate);
constexpr bool do_fake_mapping = false;
if (do_fake_mapping) {
SCATTER_FORWARD);
CHKERR VecCopy(vec_in[0], vec_out[0]);
CHKERR DMSetType(sub_dm,
"DMMOFEM");
MOFEM_LOG(
"REMESHING", Sev::inform) <<
"Created sub DpoM";
for (auto f : {"T", "TAU", "EP"}) {
}
auto fe_method = boost::make_shared<DomainEle>(m_field);
fe_method);
SCATTER_REVERSE);
} else {
}
auto ts_problem_name =
simple->getProblemName();
#ifndef NDEBUG
<< "Writing debug VTK files for remeshing verification";
"VTK", "");
"VTK", "");
"VTK", "");
"VTK", "");
}
#endif
VecScatter scatter_to_final;
for (PetscInt
i = 0;
i < nv; ++
i) {
CHKERR scatter_mng->vecScatterCreate(
CHKERR VecScatterBegin(scatter_to_final, vec_out[
i], ts_vecsout[
i],
INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecScatterEnd(scatter_to_final, vec_out[
i], ts_vecsout[
i],
INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecScatterDestroy(&scatter_to_final);
}
for (PetscInt
i = 0;
i < nv; ++
i) {
}
CHKERR bit_mng->addBitRefLevel(flipped_ents,
shift_mask.flip();
CHKERR bit_mng->shiftRightBitRef(1, shift_mask);
CHKERR TSSetIJacobian(ts,
B,
B,
nullptr,
nullptr);
}
for (PetscInt
i = 0;
i < nv; ++
i) {
double nrm;
CHKERR VecNorm(ts_vecsout[
i], NORM_2, &nrm);
<<
"After remeshing: ts_vecsout[" <<
i <<
"] norm = " << nrm;
}
PetscFunctionReturn(0);
}
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
(void *)(snes_ctx_ptr.get()), nullptr);
};
auto create_post_process_elements = [&]() {
auto pp_fe = boost::make_shared<PostProcEle>(
mField);
auto push_vol_ops = [this](auto &pip) {
"GEOMETRY");
return 1.;
};
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField,
"MAT_PLASTIC",
"MAT_THERMAL",
"MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T", 1., unity_2_args,
unity_2_args, unity_3_args, Sev::inform);
if (common_hencky_ptr) {
if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
}
return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
common_thermoplastic_ptr);
};
auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
auto &pip = pp_fe->getOpPtrVector();
auto [common_plastic_ptr, common_hencky_ptr, common_thermoplastic_ptr] =
p;
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"ISOTROPIC_HARDENING",
common_plastic_ptr->getPlasticIsoHardeningPtr()},
{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()},
{"T", common_thermoplastic_ptr->getTempPtr()}},
{{"U", u_ptr},
{"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()},
{"GEOMETRY", x_ptr}},
{{"GRAD", common_hencky_ptr->matGradPtr},
{"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
{{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()},
{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
{"HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()}}
)
);
} else {
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()},
{"T", common_thermoplastic_ptr->getTempPtr()}},
{{"U", u_ptr},
{"GEOMETRY", x_ptr},
{"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
{},
{{"STRAIN", common_plastic_ptr->mStrainPtr},
{"STRESS", common_plastic_ptr->mStressPtr},
{"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
}
};
PetscBool post_proc_vol = PETSC_FALSE;
PetscBool post_proc_skin = PETSC_TRUE;
case 2:
post_proc_vol = PETSC_TRUE;
post_proc_skin = PETSC_FALSE;
break;
case 3:
post_proc_vol = PETSC_FALSE;
post_proc_skin = PETSC_TRUE;
break;
default:
break;
}
&post_proc_vol, PETSC_NULLPTR);
&post_proc_skin, PETSC_NULLPTR);
auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
&post_proc_vol]() {
if (post_proc_vol == PETSC_FALSE)
return boost::shared_ptr<PostProcEle>();
auto pp_fe = boost::make_shared<PostProcEle>(mField);
push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
&post_proc_skin]() {
if (post_proc_skin == PETSC_FALSE)
return boost::shared_ptr<SkinPostProcEle>();
auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
pp_fe->getOpPtrVector().push_back(op_side);
pp_fe, push_vol_ops(op_side->getOpPtrVector())),
"push_vol_post_proc_ops");
return pp_fe;
};
return std::make_pair(vol_post_proc(), skin_post_proc());
};
auto scatter_create = [&](
auto D,
auto coeff) {
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULLPTR, &scatter);
};
boost::shared_ptr<SetPtsData> field_eval_data;
boost::shared_ptr<MatrixDouble> u_field_ptr;
std::array<double, 3> field_eval_coords = {0., 0., 0.};
int coords_dim = 3;
field_eval_coords.data(), &coords_dim,
boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
scalar_field_ptrs = boost::make_shared<
std::map<std::string, boost::shared_ptr<VectorDouble>>>();
boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
vector_field_ptrs = boost::make_shared<
std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
sym_tensor_field_ptrs = boost::make_shared<
std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
tensor_field_ptrs = boost::make_shared<
std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
auto test_monitor_ptr = boost::make_shared<FEMethod>();
field_eval_data =
field_eval_data,
simple->getDomainFEName());
field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
auto no_rule = [](
int,
int,
int) {
return -1; };
auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
return 1.;
};
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U", "EP",
"TAU", "T", 1., unity_2_args, unity_2_args, unity_3_args,
Sev::inform);
auto u_field_ptr = boost::make_shared<MatrixDouble>();
field_eval_fe_ptr->getOpPtrVector().push_back(
auto T_field_ptr = boost::make_shared<VectorDouble>();
field_eval_fe_ptr->getOpPtrVector().push_back(
auto FLUX_field_ptr = boost::make_shared<MatrixDouble>();
field_eval_fe_ptr->getOpPtrVector().push_back(
if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
scalar_field_ptrs->insert(std::make_pair("TEMPERATURE", T_field_ptr));
scalar_field_ptrs->insert(std::make_pair(
"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()));
scalar_field_ptrs->insert(std::make_pair(
"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()));
vector_field_ptrs->insert(std::make_pair("U", u_field_ptr));
vector_field_ptrs->insert(std::make_pair("FLUX", FLUX_field_ptr));
sym_tensor_field_ptrs->insert(std::make_pair(
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()));
sym_tensor_field_ptrs->insert(std::make_pair(
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()));
sym_tensor_field_ptrs->insert(std::make_pair(
"HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()));
sym_tensor_field_ptrs->insert(
std::make_pair("HENCKY_STRAIN", common_hencky_ptr->getMatLogC()));
tensor_field_ptrs->insert(
std::make_pair("GRAD", common_hencky_ptr->matGradPtr));
tensor_field_ptrs->insert(std::make_pair(
"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()));
} else {
scalar_field_ptrs->insert(std::make_pair("TEMPERATURE", T_field_ptr));
scalar_field_ptrs->insert(std::make_pair(
"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()));
scalar_field_ptrs->insert(std::make_pair(
"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()));
vector_field_ptrs->insert(std::make_pair("U", u_field_ptr));
vector_field_ptrs->insert(std::make_pair("FLUX", FLUX_field_ptr));
sym_tensor_field_ptrs->insert(
std::make_pair("STRAIN", common_plastic_ptr->mStrainPtr));
sym_tensor_field_ptrs->insert(
std::make_pair("STRESS", common_plastic_ptr->mStressPtr));
sym_tensor_field_ptrs->insert(std::make_pair(
"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()));
sym_tensor_field_ptrs->insert(std::make_pair(
"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()));
}
}
field_eval_data =
field_eval_data,
simple->getDomainFEName());
field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
auto no_rule = [](
int,
int,
int) {
return -1; };
auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV});
return 1.;
};
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
createCommonThermoPlasticOps<SPACE_DIM, IT, DomainEleOp>(
mField, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U", "EP",
"TAU", "T", 1, unity_2_args, unity_2_args, unity_3_args,
Sev::inform);
auto dispGradPtr = common_hencky_ptr->matGradPtr;
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto coeff_expansion_ptr = common_thermoelastic_ptr->getCoeffExpansionPtr();
auto ref_temp_ptr = common_thermoelastic_ptr->getRefTempPtr();
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
dispGradPtr));
plasticMultiplierFieldPtr = common_plastic_ptr->getPlasticTauPtr();
plasticStrainFieldPtr = common_plastic_ptr->getPlasticStrainPtr();
field_eval_fe_ptr->getOpPtrVector().push_back(
new typename H::template OpCalculateHenckyThermoPlasticStress<
SPACE_DIM,
"U", tempFieldPtr, common_hencky_ptr, coeff_expansion_ptr,
ref_temp_ptr));
field_eval_fe_ptr->getOpPtrVector().push_back(
common_hencky_ptr->getMatFirstPiolaStress(), stressFieldPtr));
field_eval_fe_ptr->getOpPtrVector().push_back(
} else {
field_eval_fe_ptr->getOpPtrVector().push_back(
new typename H::template OpCalculateLogC<SPACE_DIM, IT>(
"U", common_hencky_ptr));
stressFieldPtr = common_hencky_ptr->getMatFirstPiolaStress();
strainFieldPtr = common_hencky_ptr->getMatLogC();
};
}
auto set_time_monitor = [&](auto dm, auto solver) {
dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
} else {
test_monitor_ptr->preProcessHook = [&]() {
->evalFEAtThePoint<SPACE_DIM>(
field_eval_coords.data(), 1e-12,
simple->getProblemName(),
simple->getDomainFEName(), field_eval_data,
CHKERR VecZeroEntries(eval_num_vec);
if (tempFieldPtr->size()) {
CHKERR VecSetValue(eval_num_vec, 0, 1, ADD_VALUES);
}
CHKERR VecAssemblyBegin(eval_num_vec);
CHKERR VecAssemblyEnd(eval_num_vec);
double eval_num;
CHKERR VecSum(eval_num_vec, &eval_num);
if (!(int)eval_num) {
"atom test %d failed: did not find elements to evaluate "
"the field, check the coordinates",
}
if (tempFieldPtr->size()) {
auto t_temp = getFTensor0FromVec(*tempFieldPtr);
<< "Eval point T magnitude: " << t_temp;
if (
atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test <= 3 && fabs(t_temp - 554.48) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
if (
atom_test == 4 && fabs(t_temp - 250) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
if (
atom_test == 5 && fabs(t_temp - 1) > 1e-2) {
"atom test %d failed: wrong temperature value",
}
}
if (
atom_test == 8 && fabs(t_temp - 0.5) > 1e-12) {
"atom test %d failed: wrong temperature value",
atom_test);
}
}
if (fluxFieldPtr->size1()) {
auto t_flux = getFTensor1FromMat<SPACE_DIM>(*fluxFieldPtr);
auto flux_mag = sqrt(t_flux(
i) * t_flux(
i));
<< "Eval point FLUX magnitude: " << flux_mag;
if (
atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test <= 3 && fabs(flux_mag - 27008.0) > 2e1) {
"atom test %d failed: wrong flux value",
atom_test);
}
if (
atom_test == 4 && fabs(flux_mag - 150e3) > 1e-1) {
"atom test %d failed: wrong flux value",
atom_test);
}
if (
atom_test == 5 && fabs(flux_mag) > 1e-6) {
"atom test %d failed: wrong flux value",
atom_test);
}
}
}
if (dispFieldPtr->size1()) {
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*dispFieldPtr);
auto disp_mag = sqrt(t_disp(
i) * t_disp(
i));
<< "Eval point U magnitude: " << disp_mag;
if (
atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(disp_mag - 0.00345) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
fabs(disp_mag - 0.00265) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
fabs(t_disp(0) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5 &&
fabs(t_disp(1) - (std::sqrt(std::exp(0.2)) - 1)) > 1e-5) {
"atom test %d failed: wrong displacement value",
}
}
}
if (strainFieldPtr->size1()) {
auto t_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*strainFieldPtr);
auto t_strain_trace = t_strain(
i,
i);
if (
atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(t_strain_trace - 0.00679) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
fabs(t_strain_trace - 0.00522) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
if ((
atom_test == 5) && fabs(t_strain_trace - 0.2) > 1e-5) {
"atom test %d failed: wrong strain value",
atom_test);
}
}
}
if (stressFieldPtr->size1()) {
double von_mises_stress;
auto getVonMisesStress = [&](auto t_stress) {
von_mises_stress = std::sqrt(
0.5 *
((t_stress(0, 0) - t_stress(1, 1)) *
(t_stress(0, 0) - t_stress(1, 1)) +
(
SPACE_DIM == 3 ? (t_stress(1, 1) - t_stress(2, 2)) *
(t_stress(1, 1) - t_stress(2, 2))
: 0) +
(
SPACE_DIM == 3 ? (t_stress(2, 2) - t_stress(0, 0)) *
(t_stress(2, 2) - t_stress(0, 0))
: 0) +
6 * (t_stress(0, 1) * t_stress(0, 1) +
(
SPACE_DIM == 3 ? t_stress(1, 2) * t_stress(1, 2) : 0) +
(
SPACE_DIM == 3 ? t_stress(2, 0) * t_stress(2, 0) : 0))));
};
auto t_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*stressFieldPtr);
CHKERR getVonMisesStress(t_stress);
} else {
auto t_stress =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*stressFieldPtr);
CHKERR getVonMisesStress(t_stress);
}
<< "Eval point von Mises Stress: " << von_mises_stress;
if (
atom_test && fabs(test_monitor_ptr->ts_t - 10) < 1e-12) {
if (
atom_test == 1 && fabs(von_mises_stress - 523.0) > 5e-1) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 2 && fabs(von_mises_stress - 16.3) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 3 && fabs(von_mises_stress - 14.9) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
if (
atom_test == 5 && fabs(von_mises_stress) > 5e-2) {
"atom test %d failed: wrong von Mises stress value",
}
}
}
if (plasticMultiplierFieldPtr->size()) {
auto t_plastic_multiplier =
getFTensor0FromVec(*plasticMultiplierFieldPtr);
<< "Eval point TAU magnitude: " << t_plastic_multiplier;
if (
atom_test == 8 && fabs(t_plastic_multiplier - 1.5) > 1e-12) {
"atom test %d failed: wrong plastic multiplier value",
}
}
if (plasticStrainFieldPtr->size1()) {
auto t_plastic_strain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*plasticStrainFieldPtr);
<< "Eval point EP(0,0) = " << t_plastic_strain(0, 0)
<< ", EP(0,1) = " << t_plastic_strain(0, 1)
<< ", EP(1,1) = " << t_plastic_strain(1, 1);
if (fabs(t_plastic_strain(0, 0) - 0.005) > 1e-12) {
"atom test %d failed: wrong EP(0,0) value",
atom_test);
}
if (fabs(t_plastic_strain(0, 1) - 0.025) > 1e-12) {
"atom test %d failed: wrong EP(0,1) value",
atom_test);
}
if (fabs(t_plastic_strain(1, 1) - 0.015) > 1e-12) {
"atom test %d failed: wrong EP(1,1) value",
atom_test);
}
}
}
};
auto null = boost::shared_ptr<FEMethod>();
test_monitor_ptr, null);
};
};
auto set_schur_pc = [&](auto solver,
boost::shared_ptr<SetUpSchur> &schur_ptr) {
auto name_prb =
simple->getProblemName();
for (auto f : {"U", "FLUX"}) {
}
};
#ifdef ADD_CONTACT
for (auto f : {"SIGMA", "EP", "TAU", "T"}) {
}
#else
for (auto f : {"EP", "TAU", "T"}) {
}
#endif
};
if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
#ifdef ADD_CONTACT
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
{
{simple->getDomainFEName(),
{{"U", "U"},
{"SIGMA", "SIGMA"},
{"U", "SIGMA"},
{"SIGMA", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "EP"},
{"TAU", "U"},
{"T", "T"},
{"U", "T"},
{"T", "U"},
{"EP", "T"},
{"T", "EP"},
{"TAU", "T"},
{"T", "TAU"}
}},
{simple->getBoundaryFEName(),
{{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
}}
}
);
{dm_schur, dm_block}, block_mat_data,
{"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
true
);
};
#else
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
auto block_mat_data =
{{simple->getDomainFEName(),
{{"U", "U"},
{"EP", "EP"},
{"TAU", "TAU"},
{"U", "EP"},
{"EP", "U"},
{"EP", "TAU"},
{"TAU", "U"},
{"TAU", "EP"},
{"T", "T"},
{"T", "U"},
{"T", "FLUX"},
{"T", "EP"},
{"FLUX", "T"},
{"FLUX", "U"},
{"U", "T"},
{"EP", "T"},
{"TAU", "T"}
}}}
);
{dm_schur, dm_block}, block_mat_data,
{"EP", "TAU", "T"}, {nullptr, nullptr, nullptr}, false
);
};
#endif
auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
schur_ptr =
CHKERR schur_ptr->setUp(solver);
}
};
uXScatter = scatter_create(
D, 0);
uYScatter = scatter_create(
D, 1);
uZScatter = scatter_create(
D, 2);
auto create_solver = [pip_mng]() {
return pip_mng->createTSIM();
else
return pip_mng->createTSIM2();
};
PetscViewer viewer;
&viewer);
CHKERR VecLoad(solution_vector, viewer);
CHKERR PetscViewerDestroy(&viewer);
SCATTER_REVERSE);
}
auto solver = create_solver();
auto active_pre_lhs = []() {
};
auto active_post_lhs = [&]() {
auto get_iter = [&]() {
SNES snes;
int iter;
"Can not get iter");
return iter;
};
auto iter = get_iter();
if (iter >= 0) {
std::array<int, 5> activity_data;
std::fill(activity_data.begin(), activity_data.end(), 0);
activity_data.data(), activity_data.size(), MPI_INT,
int &active_points = activity_data[0];
int &avtive_full_elems = activity_data[1];
int &avtive_elems = activity_data[2];
int &nb_points = activity_data[3];
int &nb_elements = activity_data[4];
if (nb_points) {
double proc_nb_points =
100 * static_cast<double>(active_points) / nb_points;
double proc_nb_active =
100 * static_cast<double>(avtive_elems) / nb_elements;
double proc_nb_full_active = 100;
if (avtive_elems)
proc_nb_full_active =
100 * static_cast<double>(avtive_full_elems) / avtive_elems;
"Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
"elements %d "
"(%3.3f\%) nb full active elems %d (%3.3f\%)",
iter, nb_points, active_points, proc_nb_points,
avtive_elems, proc_nb_active, avtive_full_elems,
proc_nb_full_active, iter);
}
}
};
auto add_active_dofs_elem = [&](auto dm) {
auto fe_pre_proc = boost::make_shared<FEMethod>();
fe_pre_proc->preProcessHook = active_pre_lhs;
auto fe_post_proc = boost::make_shared<FEMethod>();
fe_post_proc->postProcessHook = active_post_lhs;
ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
};
auto set_essential_bc = [&](auto dm, auto solver) {
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
auto disp_time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
{disp_time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
mField, post_proc_rhs_ptr, 1.)();
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
};
} else {
CHKERR TS2SetSolution(solver,
D, DD);
}
auto set_up_adapt = [&](auto solver) {
TSAdapt adapt;
CHKERR TSGetAdapt(solver, &adapt);
};
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR add_active_dofs_elem(dm);
boost::shared_ptr<SetUpSchur> schur_ptr;
CHKERR set_schur_pc(solver, schur_ptr);
CHKERR set_essential_bc(dm, solver);
auto my_ctx = boost::make_shared<MyTsCtx>();
SNES snes;
CHKERR TSGetSNES(solver, &snes);
auto my_rhs = [](TS ts, PetscReal
t,
Vec u,
Vec u_t,
Vec F,
void *ctx) {
SNES snes;
double time_increment;
if (time_increment <
min_dt) {
"Minimum time increment exceeded");
}
int error = system("rm ./out_snes_residual_iter_*.h5m > /dev/null 2>&1");
}
double ts_dt;
"Cannot get timestep from TS object");
"Minimum timestep exceeded");
}
auto post_proc = [&](auto dm, auto f_res, auto sol, auto sol_dot,
auto out_name) {
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
auto &pip = post_proc_fe->getOpPtrVector();
"GEOMETRY");
auto disp_res_mat = boost::make_shared<MatrixDouble>();
auto tau_res_vec = boost::make_shared<VectorDouble>();
auto plastic_strain_res_mat = boost::make_shared<MatrixDouble>();
auto flux_res_mat = boost::make_shared<MatrixDouble>();
auto temp_res_vec = boost::make_shared<VectorDouble>();
"U", disp_res_mat, smart_f_res));
pip.push_back(
"EP", plastic_strain_res_mat, smart_f_res));
pip.push_back(
auto disp_mat = boost::make_shared<MatrixDouble>();
auto tau_vec = boost::make_shared<VectorDouble>();
auto plastic_strain_mat = boost::make_shared<MatrixDouble>();
auto flux_mat = boost::make_shared<MatrixDouble>();
auto temp_vec = boost::make_shared<VectorDouble>();
"EP", plastic_strain_mat));
auto disp_dot_mat = boost::make_shared<MatrixDouble>();
auto tau_dot_vec = boost::make_shared<VectorDouble>();
auto plastic_strain_dot_mat = boost::make_shared<MatrixDouble>();
auto flux_dot_mat = boost::make_shared<MatrixDouble>();
auto temp_dot_vec = boost::make_shared<VectorDouble>();
"U", disp_dot_mat, smart_sol_dot));
pip.push_back(
"EP", plastic_strain_dot_mat, smart_sol_dot));
pip.push_back(
auto make_d_mat = [&]() {
};
auto push_vol_ops = [&](auto &pip) {
return 1.;
};
auto [common_plastic_ptr, common_hencky_ptr, common_thermal_ptr,
common_thermoelastic_ptr, common_thermoplastic_ptr] =
m_field, "MAT_PLASTIC", "MAT_THERMAL", "MAT_THERMOELASTIC",
"MAT_THERMOPLASTIC", pip, "U", "EP", "TAU", "T", 1.,
unity_2_args, unity_2_args, unity_3_args, Sev::inform);
if (common_hencky_ptr) {
if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
"Wrong pointer for grad");
}
return std::make_tuple(common_plastic_ptr, common_hencky_ptr,
common_thermoplastic_ptr);
};
auto push_vol_post_proc_ops = [&](auto &pp_fe, auto &&p) {
auto &pip = pp_fe->getOpPtrVector();
auto [common_plastic_ptr, common_hencky_ptr, common_thermoplastic_ptr] =
p;
auto mat_strain_ptr = boost::make_shared<MatrixDouble>();
auto mat_stress_ptr = boost::make_shared<MatrixDouble>();
auto vec_temp_ptr = boost::make_shared<VectorDouble>();
auto mat_flux_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
auto x_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"RESIDUAL_TAU", tau_res_vec},
{"RESIDUAL_T", temp_res_vec},
{"TAU", tau_vec},
{"T", temp_vec},
{"RATE_TAU", tau_dot_vec},
{"RATE_T", temp_dot_vec},
{"ISOTROPIC_HARDENING",
common_plastic_ptr->getPlasticIsoHardeningPtr()},
{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()},
{"YIELD_FUNCTION",
common_plastic_ptr->getPlasticSurfacePtr()},
{"T", common_thermoplastic_ptr->getTempPtr()}},
{{"RESIDUAL_U", disp_res_mat},
{"RATE_U", disp_dot_mat},
{"U", u_ptr},
{"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()},
{"GEOMETRY", x_ptr}},
{{"GRAD", common_hencky_ptr->matGradPtr},
{"FIRST_PIOLA",
common_hencky_ptr->getMatFirstPiolaStress()}},
{{"RESIDUAL_EP", plastic_strain_res_mat},
{"RATE_EP", plastic_strain_dot_mat},
{"PLASTIC_STRAIN",
common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()},
{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
{"HENCKY_STRESS", common_hencky_ptr->getMatHenckyStress()}}
)
);
} else {
pip.push_back(
pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
{{"RESIDUAL_TAU", tau_res_vec},
{"RESIDUAL_T", temp_res_vec},
{"TAU", tau_vec},
{"T", temp_vec},
{"RATE_TAU", tau_dot_vec},
{"RATE_T", temp_dot_vec},
{"ISOTROPIC_HARDENING",
common_plastic_ptr->getPlasticIsoHardeningPtr()},
{"PLASTIC_SURFACE",
common_plastic_ptr->getPlasticSurfacePtr()},
{"PLASTIC_MULTIPLIER",
common_plastic_ptr->getPlasticTauPtr()},
{"YIELD_FUNCTION",
common_plastic_ptr->getPlasticSurfacePtr()},
{"T", common_thermoplastic_ptr->getTempPtr()}},
{{"RESIDUAL_U", disp_res_mat},
{"U", disp_mat},
{"RATE_U", disp_dot_mat},
{"U", u_ptr},
{"GEOMETRY", x_ptr},
{"FLUX", common_thermoplastic_ptr->getHeatFluxPtr()}},
{},
{{"RESIDUAL_PLASTIC_STRAIN", plastic_strain_res_mat},
{"PLASTIC_STRAIN", plastic_strain_mat},
{"RATE_PLASTIC_STRAIN", plastic_strain_dot_mat},
{"STRAIN", common_plastic_ptr->mStrainPtr},
{"STRESS", common_plastic_ptr->mStressPtr},
{"PLASTIC_STRAIN",
common_plastic_ptr->getPlasticStrainPtr()},
{"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
)
);
}
};
CHK_MOAB_THROW(push_vol_post_proc_ops(post_proc_fe, push_vol_ops(pip)),
"push_vol_post_proc_ops");
post_proc_fe);
post_proc_fe->writeFile(out_name);
};
"out_snes_residual_iter_" +
return set_RHS;
};
PetscBool is_output_residual_fields = PETSC_FALSE;
&is_output_residual_fields, PETSC_NULLPTR);
if (is_output_residual_fields == PETSC_TRUE) {
CHKERR TSSetIFunction(solver, PETSC_NULLPTR, my_rhs, my_ctx.get());
}
auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
ptr->ExRawPtr = this;
Range no_all_old_els, no_new_els, no_flipped_els;
no_new_els, no_flipped_els};
PetscBool rollback =
PETSC_TRUE;
ctx);
CHKERR TSSetApplicationContext(solver, (
void *)ctx);
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp";
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSetUp <= done";
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve";
CHKERR TSSolve(solver, PETSC_NULLPTR);
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"TSSolve <= done";
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
#ifdef ADD_CONTACT
#ifdef ENABLE_PYTHON_BINDING
Py_Initialize();
np::initialize();
#endif
#endif
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PLASTICITY"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "THERMAL"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "THERMOPLASTICITY"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "REMESHING"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
LogManager::setLog("PLASTICITY");
LogManager::setLog("THERMAL");
LogManager::setLog("THERMOPLASTICITY");
LogManager::setLog("REMESHING");
#ifdef ADD_CONTACT
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
#endif
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
} else {
char meshFileName[255];
meshFileName, 255, PETSC_NULLPTR);
}
}
#ifdef ADD_CONTACT
#ifdef ENABLE_PYTHON_BINDING
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif
return 0;
}
"Is expected that schur matrix is not "
"allocated. This is "
"possible only is PC is set up twice");
}
}
private:
};
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
"Is expected that schur matrix is not "
"allocated. This is "
"possible only is PC is set up twice");
}
DM solver_dm;
CHKERR TSGetDM(solver, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t, PetscReal
a,
Mat
A, Mat
B,
void *ctx) {
};
CHKERR TSSetIJacobian(solver,
A, P, swap_assemble, ts_ctx_ptr.get());
} else {
auto swap_assemble = [](TS ts, PetscReal
t,
Vec u,
Vec u_t,
Vec utt,
PetscReal
a, PetscReal aa, Mat
A, Mat
B,
void *ctx) {
};
CHKERR TSSetI2Jacobian(solver,
A, P, swap_assemble, ts_ctx_ptr.get());
}
auto set_ops = [&]() {
#ifndef ADD_CONTACT
pip_mng->getOpBoundaryLhsPipeline().push_front(
{
"EP",
"TAU",
"T"}, {
nullptr,
nullptr,
nullptr},
aoSchur,
S,
false,
false
));
pip_mng->getOpDomainLhsPipeline().push_front(
{
"EP",
"TAU",
"T"}, {
nullptr,
nullptr,
nullptr},
aoSchur,
S,
false,
false
));
#else
double eps_stab = 1e-4;
PETSC_NULLPTR);
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA", [eps_stab](double, double, double) {
return eps_stab;
}));
{"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
));
pip_mng->getOpDomainLhsPipeline().push_front(
{"SIGMA", "EP", "TAU", "T"}, {nullptr, nullptr, nullptr, nullptr},
));
#endif
};
auto set_assemble_elems = [&]() {
auto schur_asmb_pre_proc = boost::make_shared<FEMethod>();
schur_asmb_pre_proc->preProcessHook = [this]() {
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble Begin";
};
auto schur_asmb_post_proc = boost::make_shared<FEMethod>();
schur_asmb_post_proc->postProcessHook = [this, schur_asmb_post_proc]() {
MOFEM_LOG(
"TIMER", Sev::verbose) <<
"Lhs Assemble End";
CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
};
ts_ctx_ptr->getPreProcessIJacobian().push_front(schur_asmb_pre_proc);
ts_ctx_ptr->getPostProcessIJacobian().push_front(schur_asmb_post_proc);
};
auto set_pc = [&]() {
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
};
auto set_diagonal_pc = [&]() {
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULLPTR, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
true);
};
};
} else {
pip_mng->getOpBoundaryLhsPipeline().push_front(
pip_mng->getOpBoundaryLhsPipeline().push_back(
pip_mng->getOpDomainLhsPipeline().push_back(
}
}
boost::shared_ptr<SetUpSchur>
return boost::shared_ptr<SetUpSchur>(
}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr int SPACE_DIM
[Define dimension]
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define BITREFLEVEL_SIZE
max number of refinements
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
Order displacement.
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
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)
Get smart vector from DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMSubDMSetUp_MoFEM(DM subdm)
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
PetscErrorCode DMSetUp_MoFEM(DM dm)
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
virtual MoFEMErrorCode modify_problem_mask_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set dof mask ref level for problem
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
MoFEMErrorCode pushMarkDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true)
Mark DOFs on block entities for boundary conditions.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
auto createKSP(MPI_Comm comm)
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.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
auto createSNES(MPI_Comm comm)
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMSubData(DM dm)
Get sub problem data structure.
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
auto createDMNestSchurMat(DM dm)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpSchurAssembleBase * createOpSchurAssembleBegin()
auto createDMBlockMat(DM dm)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
FTensor::Index< 'J', 3 > J
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< SPACE_DIM > OpHdivT
Integrate Lhs div of base of flux times base of temperature (FLUX x T) and transpose of it,...
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBaseDotT
Integrate Rhs base of temperature time heat capacity times heat rate (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 1, SPACE_DIM > OpHDivTemp
Integrate Rhs div flux base times temperature (T)
MoFEMErrorCode addMatThermalBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermalParameters > blockedParamsPtr, double default_heat_conductivity, double default_heat_capacity, double default_thermal_conductivity_scale, double default_thermal_capacity_scale, Sev sev, ScalerFunTwoArgs thermal_conductivity_scaling_func, ScalerFunTwoArgs heat_capacity_scaling_func)
MoFEMErrorCode addMatThermoElasticBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string block_name, boost::shared_ptr< BlockedThermoElasticParameters > blockedParamsPtr, double default_coeff_expansion, double default_ref_temp, Sev sev)
MoFEMErrorCode addMatThermoPlasticBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string thermoplastic_block_name, boost::shared_ptr< ThermoPlasticBlockedParameters > blockedParamsPtr, boost::shared_ptr< ThermoElasticOps::BlockedThermalParameters > &blockedThermalParamsPtr, Sev sev, ScalerFunThreeArgs inelastic_heat_fraction_scaling)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr IntegrationType I
constexpr AssemblyType A
[Define dimension]
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr double t
plate stiffness
constexpr auto field_name
#define QUAD_2D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_2D_TABLE[]
OpBaseImpl< PETSC, EdgeEleOp > OpBase
constexpr double heat_conductivity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
OpBaseDotH OpBaseDivFlux
Integrate Rhs base of temperature times divergent of flux (T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 3, 3, 1 > OpHdivFlux
Integrating Rhs flux base (1/k) flux (FLUX)
BoundaryNaturalBC::OpFlux< NaturalTemperatureMeshsets, 3, SPACE_DIM > OpTemperatureBC
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpCapacity
Integrate Lhs base of temperature times (heat capacity) times base of temperature (T x T)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
void temp(int x, int y=10)
Rhs for testing EP mapping with initial conditions.
double getScale(const double time)
Get scaling at given time.
MoFEMErrorCode thermalBC()
[Mechanical boundary conditions]
boost::shared_ptr< VectorDouble > tempFieldPtr
MoFEMErrorCode opThermoPlasticFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, std::string thermal_block_name, std::string thermoelastic_block_name, std::string thermoplastic_block_name, Pip &pip, std::string u, std::string ep, std::string tau, std::string temperature)
boost::shared_ptr< DomainEle > reactionFe
boost::shared_ptr< MatrixDouble > strainFieldPtr
boost::shared_ptr< MatrixDouble > dispGradPtr
FieldApproximationBase base
Choice of finite element basis functions
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode mechanicalBC()
[Initial conditions]
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< VectorDouble > plasticMultiplierFieldPtr
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Run problem]
boost::shared_ptr< MatrixDouble > stressFieldPtr
MoFEM::Interface & mField
Reference to MoFEM interface.
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
boost::shared_ptr< MatrixDouble > fluxFieldPtr
MoFEMErrorCode setupProblem()
[Run problem]
boost::shared_ptr< MatrixDouble > dispFieldPtr
boost::shared_ptr< MatrixDouble > plasticStrainFieldPtr
MoFEMErrorCode opThermoPlasticFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, std::string thermal_block_name, std::string thermoelastic_block_name, std::string thermoplastic_block_name, Pip &pip, std::string u, std::string ep, std::string tau, std::string temperature)
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
MoFEMErrorCode initialConditions()
[Create common data]
auto createCommonThermoPlasticOps(MoFEM::Interface &m_field, std::string plastic_block_name, std::string thermal_block_name, std::string thermoelastic_block_name, std::string thermoplastic_block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, std::string temperature, double scale, ScalerFunTwoArgs thermal_conductivity_scaling, ScalerFunTwoArgs heat_capacity_scaling, ScalerFunThreeArgs inelastic_heat_fraction_scaling, Sev sev, bool with_rates=true)
Boundary condition manager for finite element problem setup.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to calculate residual side diagonal.
Class (Function) to enforce essential constrains.
Field evaluator interface.
Definition of the heat flux bc data structure.
Section manager is used to create indexes and sections.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Natural boundary conditions.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision scalar field values calculation.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Operator for symmetrizing tensor fields.
Unset indices on entities on finite element.
Template struct for dimension-specific finite element types.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Get domain left-hand side finite element.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
MoFEMErrorCode getOptions()
get options
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
MoFEM::Interface & mField
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Volume finite element base.
[Push operators to pipeline]
boost::shared_ptr< MoFEM::TsCtx > dmts_ctx
PetscInt snes_iter_counter
static std::array< int, 5 > activityData
MoFEM::Interface * m_field
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
static std::map< long int, MatrixDouble > mapRefCoords
SetEnrichedIntegration(boost::shared_ptr< Range > new_els)
SmartPetscObj< DM > subDM
field split sub dm
virtual ~SetUpSchurImpl()=default
SmartPetscObj< AO > aoSchur
SmartPetscObj< IS > fieldSplitIS
IS for split Schur block.
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
MoFEMErrorCode postProc()
MoFEM::Interface & mField
[Push operators to pipeline]
virtual MoFEMErrorCode setUp(SmartPetscObj< KSP >)=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
MoFEMErrorCode reSetUp(std::vector< std::string > fIelds, std::vector< int > oRders, BitRefLevel bIt)
MoFEMErrorCode mapFields(SmartPetscObj< DM > &intermediateDM, Range elsToRemove=Range(), Range elsToAdd=Range(), Range elsToFlip=Range(), PetscInt numVecs=1, Vec vecsToMapFrom[]=PETSC_NULLPTR, Vec vecsToMapTo[]=PETSC_NULLPTR)
Set of functions called by PETSc solver used to refine and update mesh.
virtual ~TSPrePostProc()=default
static MoFEMErrorCode tsPreProc(TS ts)
Pre process time step.
static MoFEMErrorCode tsPostProc(TS ts)
Post process time step.
static MoFEMErrorCode tsPreStage(TS ts)
Pre-stage processing for time stepping.
MoFEMErrorCode tsSetUp(TS ts)
Used to setup TS solver.
constexpr AssemblyType AT
constexpr IntegrationType IT
double iso_hardening_dtemp(double tau, double H, double omega_0, double Qinf, double omega_h, double b_iso, double sigmaY, double temp_0, double temp)
double default_heat_capacity_scale
constexpr AssemblyType AT
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
int post_processing_counter
double exp_softening(double temp_hat)
constexpr IntegrationType IT
auto init_T
Initialisation function for temperature field.
double default_thermal_conductivity_scale
ScalerFunTwoArgs heat_capacity_scaling
double dH_thermal_dtemp(double H, double omega_h)
ScalerFunThreeArgs inelastic_heat_fraction_scaling
PetscBool is_distributed_mesh
double dy_0_thermal_dtemp(double sigmaY, double omega_0)
boost::function< double(const double, const double, const double)> ScalerFunThreeArgs
boost::function< double(const double, const double)> ScalerFunTwoArgs
Example * thermoplastic_raw_ptr
PetscErrorCode MyTSResizeSetup(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
static std::unordered_map< TS, ResizeCtx * > ts_to_rctx
double Qinf_thermal(double Qinf, double omega_h, double temp_0, double temp)
const char *const ICTypes[]
double temp_hat(double temp)
auto postProcessHere(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, std::string output_name, int &counter= *(new int(0)))
int flux_order
Order of tau field.
MoFEMErrorCode TSIterationPreStage(TS ts, PetscReal stagetime)
int T_order
Order of ep field.
ElementsAndOps< SPACE_DIM >::DomainEleOnRange DomainEleOnRange
double y_0_thermal(double sigmaY, double omega_0, double temp_0, double temp)
double width
Width of Gaussian distribution.
auto uniform_distribution
auto printOnAllCores(MoFEM::Interface &m_field, const std::string &out_put_string, const auto &out_put_quantity)
MoFEMErrorCode SNESIterationMonitor(SNES snes, PetscInt its, PetscReal norm, void *ctx)
double omega_0
flow stress softening
PetscErrorCode MyTSResizeTransfer(TS ts, PetscInt nv, Vec ts_vecsin[], Vec ts_vecsout[], void *ctx)
auto postProcessPETScHere(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, Vec sol, std::string output_name)
double H_thermal(double H, double omega_h, double temp_0, double temp)
auto Gaussian_distribution
int geom_order
Order if fixed.
ScalerFunTwoArgs thermal_conductivity_scaling
double inelastic_heat_fraction
fraction of plastic dissipation converted to heat
double temp_0
reference temperature for thermal softening
static std::unordered_map< TS, MyTsCtx * > ts_to_ctx
double dQinf_thermal_dtemp(double Qinf, double omega_h)
double omega_h
hardening softening
NaturalBC< BoundaryEleOp >::Assembly< AT >::LinearForm< IT > BoundaryRhsBCs
double young_modulus
Young modulus.
double C1_k
Kinematic hardening.
BoundaryLhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
BoundaryRhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
double Qinf
Saturation yield stress.
FieldEvaluatorInterface::SetPtsData SetPtsData
ElementsAndOps< SPACE_DIM >::SideEle SideEle
#define EXECUTABLE_DIMENSION
PetscBool do_eval_field
Evaluate field.
PetscBool is_quasi_static
[Operators used for contact]
double visH
Viscous hardening.
double poisson_ratio
Poisson ratio.
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
PetscBool set_timer
Set timer.
DomainRhsBCs::OpFlux< PlasticOps::DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
double zeta
Viscous hardening.
NaturalBC< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT > BoundaryLhsBCs
int tau_order
Order of tau files.
double iso_hardening_exp(double tau, double b_iso)
double b_iso
Saturation exponent.
PetscBool is_large_strains
Large strains.
int geom_order
Order if fixed.
double sigmaY
Yield stress.
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
auto kinematic_hardening_dplastic_strain(double C1_k)
NaturalBC< DomainEleOp >::Assembly< AT >::LinearForm< IT > DomainRhsBCs
int ep_order
Order of ep files.
PostProcBrokenMeshInMoab< BoundaryEle > SkinPostProcEle
constexpr FieldSpace CONTACT_SPACE
[Specialisation for assembly]
#define FINITE_DEFORMATION_FLAG
NaturalBC< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS > DomainNaturalBCRhs
[Thermal problem]
BoundaryNaturalBC::OpFlux< NaturalForceMeshsets, 1, SPACE_DIM > OpForce
constexpr bool IS_LARGE_STRAINS
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, 1 > OpHeatSource
double default_heat_capacity
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpEssentialLhs< HeatFluxCubitBcData, 3, SPACE_DIM > OpEssentialFluxLhs
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce
DomainNaturalBCLhs::OpFlux< SetTargetTemperature, 1, 1 > OpSetTemperatureLhs
NaturalBC< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS > DomainNaturalBCLhs
DomainNaturalBCRhs::OpFlux< SetTargetTemperature, 1, 1 > OpSetTemperatureRhs
double default_coeff_expansion
EssentialBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpEssentialRhs< HeatFluxCubitBcData, 3, SPACE_DIM > OpEssentialFluxRhs
[Natural boundary conditions]
double default_heat_conductivity
NaturalBC< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS > BoundaryNaturalBC
[Body and heat source]
constexpr int SPACE_DIM
[Define dimension]
constexpr AssemblyType A
[Define dimension]
