The example shows how to solve the linear elastic problem. An example can read file with temperature field, then thermal stresses are included.
code. Analysis of the arch dam of Susqueda, located in Catalonia (Spain)" width=800px
This is an example of application code; it does not show how elements are implemented. Example presents how to:
If you like to see how to implement finite elements, material, are other parts of the code, look here;
using namespace boost::numeric;
static char help[] =
"-my_block_config set block data\n"
"-my_order approximation order\n"
"-my_is_partitioned set if mesh is partitioned\n"
"\n";
int oRder;
int iD;
double yOung;
double pOisson;
double initTemp;
};
int operator()(
int,
int,
int order)
const {
return 2 *
order; }
};
int getRuleTrianglesOnly(
int order);
int getRuleThroughThickness(
int order);
};
int main(
int argc,
char *argv[]) {
const string default_options = "-ksp_type gmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-ksp_monitor \n"
"-snes_type newtonls \n"
"-snes_linesearch_type basic \n"
"-snes_atol 1e-8 \n"
"-snes_rtol 1e-8 \n"
"-snes_monitor \n"
"-ts_monitor \n"
"-ts_type beuler \n";
string param_file = "param_file.petsc";
if (!static_cast<bool>(ifstream(param_file))) {
std::ofstream file(param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file << default_options;
file.close();
}
}
auto core_log = logging::core::get();
core_log->add_sink(
core_log->add_sink(
try {
PetscBool flg_block_config, flg_file;
char block_config_file[255];
PetscInt test_nb = 0;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_calculating_frequency = PETSC_FALSE;
PetscBool is_post_proc_volume = PETSC_TRUE;
enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
"jacobi"};
PetscInt choice_base_value = LOBATTO;
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
CHKERR PetscOptionsInt(
"-my_order",
"default approximation order",
"",
CHKERR PetscOptionsEList(
"-base",
"approximation base",
"", list_bases,
LASBASETOP, list_bases[choice_base_value],
&choice_base_value, PETSC_NULL);
CHKERR PetscOptionsInt(
"-is_atom_test",
"ctest number",
"", test_nb,
&test_nb, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_partitioned",
"set if mesh is partitioned (this result that each "
"process keeps only one part of the mesh)",
"", is_partitioned, &is_partitioned, PETSC_NULL);
CHKERR PetscOptionsString(
"-my_block_config",
"elastic configure file name",
"", "block_conf.in", block_config_file, 255,
&flg_block_config);
"-my_is_calculating_frequency", "set if frequency will be calculated",
"", is_calculating_frequency, &is_calculating_frequency, PETSC_NULL);
CHKERR PetscOptionsBool(
"-my_is_post_proc_volume",
"if true post proc volume", "", is_post_proc_volume,
&is_post_proc_volume, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
MPI_Comm moab_comm_world;
MPI_Comm_dup(PETSC_COMM_WORLD, &moab_comm_world);
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
if (pcomm == NULL)
pcomm = new ParallelComm(&moab, moab_comm_world);
if (is_partitioned == PETSC_TRUE) {
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
} else {
const char *option;
option = "";
}
bool mesh_has_tets = false;
bool mesh_has_prisms = false;
int nb_tets = 0;
int nb_hexs = 0;
int nb_prisms = 0;
CHKERR moab.get_number_entities_by_type(0, MBTET, nb_tets,
true);
CHKERR moab.get_number_entities_by_type(0, MBHEX, nb_hexs,
true);
CHKERR moab.get_number_entities_by_type(0, MBPRISM, nb_prisms,
true);
mesh_has_tets = (nb_tets + nb_hexs) > 0;
mesh_has_prisms = nb_prisms > 0;
bit_level0.set(0);
0, 3, bit_level0);
if (
bit->getName().compare(0, 3,
"ROD") == 0) {
0, 1, bit_level0);
}
}
switch (choice_base_value) {
case LEGENDRE:
break;
case LOBATTO:
break;
case BERNSTEIN_BEZIER:
break;
case JACOBI:
break;
default:
};
CHKERR m_field.
get_moab().get_entities_by_type(0, MBEDGE, all_edges,
true);
Range edges_in_simple_rod;
if (
bit->getName().compare(0, 3,
"ROD") == 0) {
MBEDGE, edges, true);
edges_in_simple_rod.merge(edges);
}
}
"DISPLACEMENT");
Range edges_to_set_order;
edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
else
auto setting_second_order_geometry = [&m_field]() {
moab::Interface::UNION);
};
CHKERR setting_second_order_geometry();
std::map<int, BlockOptionData> block_data;
auto setting_blocks_data_and_order_from_config_file =
[&m_field, &moab, &block_data, flg_block_config, block_config_file,
order](boost::shared_ptr<std::map<int, BlockData>> &block_sets_ptr) {
if (flg_block_config) {
ifstream ini_file(block_config_file);
po::variables_map vm;
po::options_description config_file_options;
it)) {
std::ostringstream str_order;
str_order << "block_" << it->getMeshsetId()
<< ".displacement_order";
config_file_options.add_options()(
str_order.str().c_str(),
po::value<int>(&block_data[it->getMeshsetId()].oRder)
std::ostringstream str_cond;
str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
config_file_options.add_options()(
str_cond.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].yOung)
->default_value(-1));
std::ostringstream str_capa;
str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
config_file_options.add_options()(
str_capa.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].pOisson)
->default_value(-2));
std::ostringstream str_init_temp;
str_init_temp << "block_" << it->getMeshsetId()
<< ".initial_temperature";
config_file_options.add_options()(
str_init_temp.str().c_str(),
po::value<double>(&block_data[it->getMeshsetId()].initTemp)
->default_value(0));
}
po::parsed_options parsed =
parse_config_file(ini_file, config_file_options, true);
store(parsed, vm);
po::notify(vm);
it)) {
if (block_data[it->getMeshsetId()].oRder == -1)
continue;
if (block_data[it->getMeshsetId()].oRder ==
order)
continue;
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Set block %d order to %d",
it->getMeshsetId(),
block_data[it->getMeshsetId()].oRder);
CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
CHKERR moab.get_adjacencies(block_ents, 3,
false,
ents_to_set_order,
moab::Interface::UNION);
ents_to_set_order = ents_to_set_order.subset_by_dimension(3);
CHKERR moab.get_adjacencies(block_ents, 2,
false,
ents_to_set_order,
moab::Interface::UNION);
CHKERR moab.get_adjacencies(block_ents, 1,
false,
ents_to_set_order,
moab::Interface::UNION);
ents_to_set_order);
ents_to_set_order, "DISPLACEMENT",
block_data[it->getMeshsetId()].oRder);
}
std::vector<std::string> additional_parameters;
additional_parameters =
collect_unrecognized(parsed.options, po::include_positional);
for (std::vector<std::string>::iterator vit =
additional_parameters.begin();
vit != additional_parameters.end(); vit++) {
MOFEM_LOG_C(
"ELASTIC", Sev::warning,
"Unrecognized option %s",
vit->c_str());
}
}
const int id = it->getMeshsetId();
auto &bd = (*block_sets_ptr)[id];
if (block_data[id].yOung > 0)
bd.E = block_data[id].yOung;
if (block_data[id].pOisson >= -1)
bd.PoissonRatio = block_data[id].pOisson;
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tYoung modulus %3.4g", bd.E);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"\tPoisson ratio %3.4g",
bd.PoissonRatio);
}
};
boost::shared_ptr<std::map<int, HookeElement::BlockData>> block_sets_ptr =
boost::make_shared<std::map<int, HookeElement::BlockData>>();
CHKERR setting_blocks_data_and_order_from_config_file(block_sets_ptr);
boost::shared_ptr<std::map<int, MassBlockData>> mass_block_sets_ptr =
boost::make_shared<std::map<int, MassBlockData>>();
auto fe_lhs_ptr =
boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
auto fe_rhs_ptr =
boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
fe_lhs_ptr->getRuleHook =
VolRule();
fe_rhs_ptr->getRuleHook =
VolRule();
false);
false);
boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
"DISPLACEMENT",
"MESH_NODE_POSITIONS", false);
auto add_skin_element_for_post_processing = [&]() {
Range elastic_element_ents;
"ELASTIC", 3, elastic_element_ents);
CHKERR skin.find_skin(0, elastic_element_ents,
false, skin_faces);
if (is_partitioned) {
CHKERR pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
} else {
proc_skin = skin_faces;
}
"DISPLACEMENT");
"DISPLACEMENT");
"DISPLACEMENT");
"TEMP");
}
"POST_PROC_SKIN", "MESH_NODE_POSITIONS");
"POST_PROC_SKIN");
};
CHKERR add_skin_element_for_post_processing();
auto data_at_pts = boost::make_shared<HookeElement::DataAtIntegrationPts>();
if (mesh_has_tets) {
"DISPLACEMENT", "MESH_NODE_POSITIONS",
false, true, MBTET, data_at_pts);
}
if (mesh_has_prisms) {
prism_fe_lhs_ptr, prism_fe_rhs_ptr, block_sets_ptr, "DISPLACEMENT",
"MESH_NODE_POSITIONS", false, true, MBPRISM, data_at_pts);
}
if (test_nb == 4) {
auto thermal_strain =
constexpr double alpha = 1;
t_thermal_strain(
i,
j) = alpha * t_coords(2) *
t_kd(
i,
j);
return t_thermal_strain;
};
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAnalyticalInternalStrain_dx<0>(
"DISPLACEMENT", data_at_pts, thermal_strain));
}
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_mass_ptr(
for (auto &sit : *block_sets_ptr) {
for (auto &mit : *mass_block_sets_ptr) {
fe_mass_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateMassMatrix("DISPLACEMENT",
"DISPLACEMENT", sit.second,
mit.second, data_at_pts));
}
}
"MESH_NODE_POSITIONS");
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
fe_spring_rhs_ptr, "DISPLACEMENT",
"MESH_NODE_POSITIONS");
"MESH_NODE_POSITIONS");
boost::shared_ptr<EdgeEle> fe_simple_rod_lhs_ptr(
new EdgeEle(m_field));
boost::shared_ptr<EdgeEle> fe_simple_rod_rhs_ptr(
new EdgeEle(m_field));
m_field, fe_simple_rod_lhs_ptr, fe_simple_rod_rhs_ptr, "DISPLACEMENT",
"MESH_NODE_POSITIONS");
"DISPLACEMENT");
"DISPLACEMENT");
"DISPLACEMENT");
"MESH_NODE_POSITIONS");
CHKERR m_field.
get_moab().get_entities_by_dimension(it->meshset, 3, tets,
true);
}
fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
bool add_temp_field = false;
if (block_data[it->getMeshsetId()].initTemp != 0) {
add_temp_field = true;
break;
}
}
if (add_temp_field) {
}
}
CHKERR thermal_stress_elem.addThermalStressElement(
"ELASTIC", "DISPLACEMENT", "TEMP");
}
"MESH_NODE_POSITIONS");
if (block_data[it->getMeshsetId()].initTemp != 0) {
"Set block %d temperature to %3.2g\n", it->getMeshsetId(),
block_data[it->getMeshsetId()].initTemp);
CHKERR moab.get_entities_by_handle(it->meshset, block_ents,
true);
CHKERR moab.get_connectivity(block_ents, vertices,
true);
block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
"TEMP");
}
}
}
auto dm =
createDM(PETSC_COMM_WORLD,
"MOFEM");
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
if (is_calculating_frequency == PETSC_TRUE) {
CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
}
fe_spring_lhs_ptr->ksp_B = Aij;
fe_spring_rhs_ptr->ksp_f =
F;
fe_simple_rod_lhs_ptr->ksp_B = Aij;
fe_simple_rod_rhs_ptr->ksp_f =
F;
CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
m_field,
"DISPLACEMENT", Aij, D0,
F);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
prism_fe_rhs_ptr->snes_f =
F;
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Assemble external force vector ...";
fe_lhs_ptr->snes_B = Aij;
prism_fe_lhs_ptr->snes_B = Aij;
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate stiffness matrix ...";
if (is_calculating_frequency == PETSC_TRUE) {
fe_mass_ptr->snes_B = Mij;
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Calculate mass matrix ...";
}
boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
{
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumann_forces.begin();
for (; mit != neumann_forces.end(); mit++) {
&mit->second->getLoopFe());
}
}
boost::ptr_map<std::string, NodalForce> nodal_forces;
"DISPLACEMENT");
{
boost::ptr_map<std::string, NodalForce>::iterator fit =
nodal_forces.begin();
for (; fit != nodal_forces.end(); fit++) {
&fit->second->getLoopFe());
}
}
boost::ptr_map<std::string, EdgeForce> edge_forces;
"DISPLACEMENT");
{
auto fit = edge_forces.begin();
for (; fit != edge_forces.end(); fit++) {
auto &fe = fit->second->getLoopFe();
}
}
CHKERR body_forces_methods.addBlock(
"DISPLACEMENT",
F,
it->getMeshsetId());
}
&body_forces_methods.getLoopFe());
false, false);
CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
"DISPLACEMENT",
F,
false,
true);
&fluid_pressure_fe.getLoopFe());
PetscViewerPushFormat(
PETSC_VIEWER_STDOUT_SELF,
PETSC_VIEWER_ASCII_MATLAB);
if (is_calculating_frequency == PETSC_TRUE) {
CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
}
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
CHKERR VecGhostUpdateBegin(
F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
F, ADD_VALUES, SCATTER_REVERSE);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetOperators(solver, Aij, Aij);
{
PetscBool same = PETSC_FALSE;
PC pc;
PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
if (same) {
} else {
CHKERR KSPSetDMActive(solver, PETSC_FALSE);
}
}
CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
auto set_post_proc_skin = [&](auto &post_proc_skin) {
false);
CHKERR post_proc_skin.generateReferenceElementMesh();
CHKERR post_proc_skin.addFieldValuesPostProc(
"DISPLACEMENT");
CHKERR post_proc_skin.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
"DISPLACEMENT", "ELASTIC", data_at_pts->hMat, true);
CHKERR post_proc_skin.addFieldValuesGradientPostProcOnSkin(
"MESH_NODE_POSITIONS", "ELASTIC", data_at_pts->HMat, false);
CHKERR post_proc_skin.addFieldValuesPostProc(
"TEMP");
CHKERR post_proc_skin.addFieldValuesGradientPostProc(
"TEMP");
}
post_proc_skin.getOpPtrVector().push_back(
new HookeElement::OpPostProcHookeElement<
"DISPLACEMENT", data_at_pts, *block_sets_ptr,
post_proc_skin.postProcMesh, post_proc_skin.mapGaussPts, true,
true));
};
auto set_post_proc_tets = [&](auto &post_proc) {
CHKERR post_proc.generateReferenceElementMesh();
false);
CHKERR post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
CHKERR post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
CHKERR post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
CHKERR post_proc.addFieldValuesPostProc(
"TEMP");
CHKERR post_proc.addFieldValuesGradientPostProc(
"TEMP");
}
m_field, post_proc.postProcMesh, post_proc.mapGaussPts,
"DISPLACEMENT", post_proc.commonData, block_sets_ptr.get()));
};
auto set_post_proc_edge = [&](auto &post_proc_edge) {
CHKERR post_proc_edge.generateReferenceElementMesh();
CHKERR post_proc_edge.addFieldValuesPostProc(
"DISPLACEMENT");
};
auto set_post_proc_prisms = [&](auto &prism_post_proc) {
CHKERR prism_post_proc.generateReferenceElementMesh();
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
prism_post_proc.getOpPtrVector().push_back(
prism_post_proc.getOpPtrVector().push_back(
CHKERR prism_post_proc.addFieldValuesPostProc(
"DISPLACEMENT");
CHKERR prism_post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
CHKERR prism_post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT");
m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
"DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
};
CHKERR set_post_proc_skin(post_proc_skin);
CHKERR set_post_proc_tets(post_proc);
CHKERR set_post_proc_prisms(prism_post_proc);
CHKERR set_post_proc_edge(post_proc_edge);
PetscBool field_eval_flag = PETSC_FALSE;
std::array<double, 3> field_eval_coords;
boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> field_eval_data;
PetscInt coords_dim = 3;
field_eval_coords.data(), &coords_dim,
&field_eval_flag);
auto scalar_field_ptr = boost::make_shared<VectorDouble>();
auto vector_field_ptr = boost::make_shared<MatrixDouble>();
auto tensor_field_ptr = boost::make_shared<MatrixDouble>();
if (field_eval_flag) {
->getData<VolumeElementForcesAndSourcesCore>();
field_eval_data, "ELASTIC");
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().push_back(
}
field_eval_fe_ptr->getOpPtrVector().push_back(
field_eval_fe_ptr->getOpPtrVector().push_back(
}
CHKERR VecDuplicate(
F, &F_thermal);
CHKERR thermal_stress_elem.setThermalStressRhsOperators(
"DISPLACEMENT", "TEMP", F_thermal);
sit)) {
sit->get_step_number());
sit->get_step_number());
CHKERR VecZeroEntries(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
CHKERR VecAssemblyBegin(F_thermal);
CHKERR VecAssemblyEnd(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
PetscReal nrm_F;
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
PetscReal nrm_F_thermal;
CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
nrm_F_thermal);
CHKERR VecScale(F_thermal, -1);
dirichlet_bc_ptr->snes_x =
D;
dirichlet_bc_ptr->snes_f = F_thermal;
CHKERR KSPSolve(solver, F_thermal,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
if (field_eval_flag) {
->evalFEAtThePoint3D(
field_eval_coords.data(), 1e-12, "ELASTIC_PROB", "ELASTIC",
if (scalar_field_ptr->size()) {
<< "Eval point TEMP: " << t_temp;
}
if (vector_field_ptr->size1()) {
auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
<< "Eval point DISPLACEMENT magnitude: "
<< sqrt(t_disp(
i) * t_disp(
i));
}
if (tensor_field_ptr->size1()) {
auto t_disp_grad = getFTensor2FromMat<3, 3>(*tensor_field_ptr);
<<
"Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(
i,
i);
}
}
if (is_post_proc_volume == PETSC_TRUE) {
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
std::ostringstream o1;
o1 << "out_" << sit->step_number << ".h5m";
if (!test_nb)
CHKERR post_proc.writeFile(o1.str().c_str());
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
}
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
&post_proc_skin);
std::ostringstream o1_skin;
o1_skin << "out_skin_" << sit->step_number << ".h5m";
if (!test_nb)
CHKERR post_proc_skin.writeFile(o1_skin.str().c_str());
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"done ...";
}
} else {
CHKERR VecZeroEntries(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
CHKERR VecAssemblyBegin(F_thermal);
CHKERR VecAssemblyEnd(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
PetscReal nrm_F;
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F = %6.4e", nrm_F);
PetscReal nrm_F_thermal;
CHKERR VecNorm(F_thermal, NORM_2, &nrm_F_thermal);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"norm2 F_thermal = %6.4e",
nrm_F_thermal);
CHKERR VecScale(F_thermal, -1);
dirichlet_bc_ptr->snes_x =
D;
dirichlet_bc_ptr->snes_f = F_thermal;
CHKERR KSPSolve(solver, F_thermal,
D);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
if (is_post_proc_volume == PETSC_TRUE) {
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
if (!test_nb)
CHKERR post_proc.writeFile(
"out.h5m");
}
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file skin ...";
if (!test_nb)
CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
}
CHKERR VecDestroy(&F_thermal);
} else {
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Post-process start ...";
if (is_post_proc_volume == PETSC_TRUE) {
}
MOFEM_LOG(
"ELASTIC", Sev::inform) <<
"Write output file ...";
if (mesh_has_tets) {
if (is_post_proc_volume == PETSC_TRUE) {
if (!test_nb)
CHKERR post_proc.writeFile(
"out.h5m");
}
if (!test_nb)
CHKERR post_proc_skin.writeFile(
"out_skin.h5m");
}
if (mesh_has_prisms) {
if (!test_nb)
CHKERR prism_post_proc.writeFile(
"prism_out.h5m");
}
if (!edges_in_simple_rod.empty())
if (!test_nb)
CHKERR post_proc_edge.writeFile(
"out_edge.h5m");
}
if (is_calculating_frequency == PETSC_TRUE) {
double mode_mass;
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode mass %6.4e\n", mode_mass);
double mode_stiffness;
CHKERR VecDot(v1,
D, &mode_stiffness);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Mode stiffness %6.4e\n",
mode_stiffness);
double frequency;
double pi = 3.14159265359;
frequency = std::sqrt(mode_stiffness / mode_mass) / (2 *
pi);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Frequency %6.4e", frequency);
}
auto calculate_strain_energy = [&]() {
"MESH_NODE_POSITIONS", false, true,
v_energy);
double energy;
CHKERR VecSum(v_energy, &energy);
MOFEM_LOG_C(
"ELASTIC", Sev::inform,
"Elastic energy %6.4e", energy);
switch (test_nb) {
case 1:
if (fabs(energy - 17.129) > 1e-3)
"atom test diverged!");
break;
case 2:
if (fabs(energy - 5.6475e-03) > 1e-4)
"atom test diverged!");
break;
case 3:
if (fabs(energy - 7.4679e-03) > 1e-4)
"atom test diverged!");
break;
case 4:
if (fabs(energy - 2.4992e+00) > 1e-3)
"atom test diverged!");
break;
case 8: {
double min;
CHKERR VecMin(
D, PETSC_NULL, &min);
constexpr double expected_val = 0.10001;
if (fabs(min + expected_val) > 1e-10)
"atom test diverged! %3.4e != %3.4e", min, expected_val);
} break;
case 9: {
if (fabs(energy - 4.7416e-04) > 1e-8)
"atom test diverged!");
}
default:
break;
}
};
CHKERR calculate_strain_energy();
MPI_Comm_free(&moab_comm_world);
}
return 0;
}