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 main(
int argc,
char *argv[]) {
"-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";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "ELASTIC"));
LogManager::setLog("ELASTIC");
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)");
}
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
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 HookeElement::setBlocks(m_field, block_sets_ptr);
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(
CHKERR HookeElement::addElasticElement(m_field, block_sets_ptr,
"ELASTIC",
"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) {
CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
"DISPLACEMENT", "MESH_NODE_POSITIONS",
false, true, MBTET, data_at_pts);
}
if (mesh_has_prisms) {
CHKERR HookeElement::setOperators(
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;
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);
dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
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()) {
auto t_temp = getFTensor0FromVec(*scalar_field_ptr);
<< "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 = [&]() {
CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr,
"DISPLACEMENT",
"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;
}
const std::string default_options
Implementation of linear elastic material.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, PCMGSetUpViaApproxOrdersCtx *ctx, int verb)
Function build MG structure.
static PetscErrorCode ierr
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#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 CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ BODYFORCESSET
block name is "BODY_FORCES"
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
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 DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag 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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
virtual bool check_series(const std::string &name) const
check if series is in database
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
const FTensor::Tensor2< T, Dim, Dim > Vec
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 createKSP(MPI_Comm comm)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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)
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
data for calculation inertia forces
MoFEMErrorCode setBlocks()
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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.
FatPrismElementForcesAndSourcesCore(Interface &m_field)
Field evaluator interface.
Log manager is used to build and partition problems.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Get value at integration points for scalar field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Transform local reference derivatives of shape functions to global derivatives.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
Set data structures of MG pre-conditioner via approximation orders.
Operator post-procesing stresses for Hook isotropic material.
int getRuleThroughThickness(int order)
int getRuleTrianglesOnly(int order)
Implentation of thermal stress element.
int operator()(int, int, int) const