v0.14.0
elasticity.cpp

The example shows how to solve the linear elastic problem. An example can read file with temperature field, then thermal stresses are included.

What example can do:

See example how code can be used [67],

Example what you can do with this

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;

/** \file elasticity.cpp
* \ingroup nonlinear_elastic_elem
* \example elasticity.cpp
The example shows how to solve the linear elastic problem. An example can read
file with temperature field, then thermal stresses are included.
What example can do:
- take into account temperature field, i.e. calculate thermal stresses and
deformation
- stationary and time depend field is considered
- take into account gravitational body forces
- take in account fluid pressure
- can work with higher order geometry definition
- works on distributed meshes
- multi-grid solver where each grid level is approximation order level
- each mesh block can have different material parameters and approximation
order
See example how code can be used \cite jordi:2017,
\image html SquelaDamExampleByJordi.png "Example what you can do with this
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:
- read mesh
- set-up problem
- run finite elements on the problem
- assemble matrices and vectors
- solve the linear system of equations
- save results
If you like to see how to implement finite elements, material, are other parts
of the code, look here;
- Hooke material, see \ref Hooke
- Thermal-stress assembly, see \ref ThermalElement
- Body forces element, see \ref BodyForceConstantField
- Fluid pressure element, see \ref FluidPressure
- The general implementation of an element for arbitrary Lagrangian-Eulerian
formulation for a nonlinear elastic problem is here \ref
NonlinearElasticElement. Here we limit ourselves to Hooke equation and fix
mesh, so the problem becomes linear. Not that elastic element is implemented
with automatic differentiation.
*/
using namespace MoFEM;
#include <Hooke.hpp>
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;
Range tRis;
BlockOptionData() : oRder(-1), yOung(-1), pOisson(-2), initTemp(0) {}
};
/// Set integration rule
struct VolRule {
int operator()(int, int, int order) const { return 2 * order; }
};
int getRuleTrianglesOnly(int order);
int getRuleThroughThickness(int order);
};
int PrismFE::getRuleTrianglesOnly(int order) { return 2 * order; };
int PrismFE::getRuleThroughThickness(int order) { return 2 * 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();
}
}
MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::setLog("ELASTIC");
MOFEM_LOG_TAG("ELASTIC", "elasticity")
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "ELASTIC_SYNC"));
LogManager::setLog("ELASTIC_SYNC");
MOFEM_LOG_TAG("ELASTIC_SYNC", "elastic_sync");
try {
PetscBool flg_block_config, flg_file;
char mesh_file_name[255];
char block_config_file[255];
PetscInt test_nb = 0;
PetscInt order = 2;
PetscBool is_partitioned = PETSC_FALSE;
PetscBool is_calculating_frequency = PETSC_FALSE;
PetscBool is_post_proc_volume = PETSC_TRUE;
// Select base
enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, JACOBI, LASBASETOP };
const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier",
"jacobi"};
PetscInt choice_base_value = LOBATTO;
// Read options from command line
ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_file);
CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
order, &order, PETSC_NULL);
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);
CHKERR PetscOptionsBool(
"-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();
// Throw error if file with mesh is not provided
if (flg_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
// Create mesh database
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
// Create moab communicator
// Create separate MOAB communicator, it is duplicate of PETSc communicator.
// NOTE That this should eliminate potential communication problems between
// MOAB and PETSC functions.
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);
// Read whole mesh or part of it if partitioned
if (is_partitioned == PETSC_TRUE) {
// This is a case of distributed mesh and algebra. In this case each
// processor keeps only one part of the problem.
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR moab.load_file(mesh_file_name, 0, option);
} else {
// If that case we have distributed algebra, i.e. assembly of vectors and
// matrices is in parallel, but whole mesh is stored on all processors.
// Solver and matrix scales well, however problem set-up of problem is
// not fully parallel.
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
}
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
// Print boundary conditions and material parameters
MeshsetsManager *meshsets_mng_ptr;
CHKERR m_field.getInterface(meshsets_mng_ptr);
CHKERR meshsets_mng_ptr->printDisplacementSet();
CHKERR meshsets_mng_ptr->printForceSet();
CHKERR meshsets_mng_ptr->printMaterialsSet();
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;
// Set bit refinement level to all entities (we do not refine mesh in
// this example so all entities have the same bit refinement level)
BitRefLevel bit_level0;
bit_level0.set(0);
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 3, bit_level0);
// CHECK IF EDGE BLOCKSET EXIST AND IF IT IS ADD ALL ENTITIES FROM IT
// CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
// MESHSET_OF_EDGE_BLOCKSET, 1, bit_level0);
if (bit->getName().compare(0, 3, "ROD") == 0) {
CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
0, 1, bit_level0);
}
}
// Declare approximation fields
switch (choice_base_value) {
case LEGENDRE:
break;
case LOBATTO:
break;
case BERNSTEIN_BEZIER:
break;
case JACOBI:
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Base not implemented");
};
CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_DENSE,
// We can use higher oder geometry to define body
CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, base, 3, MB_TAG_DENSE,
// Declare problem
// Add entities (by tets) to the field (all entities in the mesh, root_set
// = 0 )
CHKERR m_field.add_ents_to_field_by_dim(0, 3, "DISPLACEMENT");
CHKERR m_field.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
// Get all edges in the mesh
Range all_edges;
CHKERR m_field.get_moab().get_entities_by_type(0, MBEDGE, all_edges, true);
// Get edges associated with simple rod
Range edges_in_simple_rod;
if (bit->getName().compare(0, 3, "ROD") == 0) {
Range edges;
CHKERR m_field.get_moab().get_entities_by_type(bit->getMeshset(),
MBEDGE, edges, true);
edges_in_simple_rod.merge(edges);
}
}
CHKERR m_field.add_ents_to_field_by_type(edges_in_simple_rod, MBEDGE,
"DISPLACEMENT");
// Set order of edge in rod to be 1
CHKERR m_field.set_field_order(edges_in_simple_rod, "DISPLACEMENT", 1);
// Get remaining edges (not associated with simple rod) to set order
Range edges_to_set_order;
edges_to_set_order = subtract(all_edges, edges_in_simple_rod);
// Set approximation order.
// See Hierarchical Finite Element Bases on Unstructured Tetrahedral
// Meshes.
CHKERR m_field.set_field_order(0, MBPRISM, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(0, MBHEX, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(0, MBQUAD, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(edges_to_set_order, "DISPLACEMENT", order);
CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", order);
else
CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
// Set order of approximation of geometry.
// Apply 2nd order only on skin (or in whole body)
auto setting_second_order_geometry = [&m_field]() {
// Setting geometry order everywhere
Range tets, edges;
CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, tets);
CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, edges,
moab::Interface::UNION);
// Setting 2nd geometry order only on skin
// Range tets, faces, edges;
// Skinner skin(&m_field.get_moab());
// CHKERR skin.find_skin(0,tets,false,faces);
// CHKERR m_field.get_moab().get_adjacencies(
// faces,1,false,edges,moab::Interface::UNION
// );
// CHKERR
// m_field.getInterface<CommInterface>()->synchroniseEntities(edges);
CHKERR m_field.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
};
CHKERR setting_second_order_geometry();
// Configure blocks by parsing config file. It allows setting
// approximation order for each block independently.
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)
->default_value(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);
Range block_ents;
CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
true);
Range ents_to_set_order;
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);
CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
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());
}
}
// Update material parameters. Set material parameters block by
// block.
m_field, BLOCKSET | MAT_ELASTICSET, it)) {
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, "Block %d", id);
MOFEM_LOG_C("ELASTIC", Sev::inform, "\tYoung modulus %3.4g", bd.E);
MOFEM_LOG_C("ELASTIC", Sev::inform, "\tPoisson ratio %3.4g",
bd.PoissonRatio);
}
};
// Add elastic element
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>>();
CHKERR ConvectiveMassElement::setBlocks(m_field, mass_block_sets_ptr);
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();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_lhs_ptr, true, false, false,
false);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fe_rhs_ptr, true, false, false,
false);
boost::shared_ptr<ForcesAndSourcesCore> prism_fe_lhs_ptr(
new PrismFE(m_field));
boost::shared_ptr<ForcesAndSourcesCore> prism_fe_rhs_ptr(
new PrismFE(m_field));
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);
Skinner skin(&m_field.get_moab());
Range skin_faces; // skin faces from 3d ents
CHKERR skin.find_skin(0, elastic_element_ents, false, skin_faces);
Range proc_skin;
if (is_partitioned) {
CHKERR pcomm->filter_pstatus(skin_faces,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &proc_skin);
} else {
proc_skin = skin_faces;
}
CHKERR m_field.add_finite_element("POST_PROC_SKIN");
CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
"DISPLACEMENT");
CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
"DISPLACEMENT");
CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
"DISPLACEMENT");
if (m_field.check_field("TEMP")) {
// CHKERR m_field.modify_finite_element_add_field_row("POST_PROC_SKIN",
// "TEMP");
// CHKERR m_field.modify_finite_element_add_field_col("POST_PROC_SKIN",
// "TEMP");
CHKERR m_field.modify_finite_element_add_field_data("POST_PROC_SKIN",
"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) {
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));
}
}
// Add spring boundary condition applied on surfaces.
// This is only declaration not implementation.
CHKERR MetaSpringBC::addSpringElements(m_field, "DISPLACEMENT",
"MESH_NODE_POSITIONS");
// Implementation of spring element
// Create new instances of face elements for springs
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
CHKERR MetaSpringBC::setSpringOperators(m_field, fe_spring_lhs_ptr,
fe_spring_rhs_ptr, "DISPLACEMENT",
"MESH_NODE_POSITIONS");
// Add Simple Rod elements
// This is only declaration not implementation.
"MESH_NODE_POSITIONS");
// CHKERR m_field.add_ents_to_finite_element_by_type(edges_in_simple_rod,
// MBEDGE, "SIMPLE_ROD");
// Implementation of Simple Rod element
// Create new instances of edge elements for Simple Rod
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");
// Add body force element. This is only declaration of element. not its
// implementation.
CHKERR m_field.add_finite_element("BODY_FORCE");
"DISPLACEMENT");
"DISPLACEMENT");
"DISPLACEMENT");
"MESH_NODE_POSITIONS");
m_field, BLOCKSET | BODYFORCESSET, it)) {
Range tets;
CHKERR m_field.get_moab().get_entities_by_dimension(it->meshset, 3, tets,
true);
CHKERR m_field.add_ents_to_finite_element_by_dim(tets, 3, "BODY_FORCE");
}
CHKERR m_field.build_finite_elements("BODY_FORCE");
// Add Neumann forces, i.e. pressure or traction forces applied on body
// surface. This is only declaration not implementation.
CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
// Add fluid pressure finite elements. This is special pressure on the
// surface from fluid, i.e. pressure which linearly change with the depth.
FluidPressure fluid_pressure_fe(m_field);
// This function only declare element. Element is implemented by operators
// in class FluidPressure.
fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
// Add elements for thermo elasticity if temperature field is defined.
ThermalStressElement thermal_stress_elem(m_field);
// Check if TEMP field exist, and then add element.
if (!m_field.check_field("TEMP")) {
bool add_temp_field = false;
if (block_data[it->getMeshsetId()].initTemp != 0) {
add_temp_field = true;
break;
}
}
if (add_temp_field) {
MB_TAG_SPARSE, MF_ZERO);
CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "TEMP");
CHKERR m_field.set_field_order(0, MBVERTEX, "TEMP", 1);
}
}
if (m_field.check_field("TEMP")) {
CHKERR thermal_stress_elem.addThermalStressElement(
"ELASTIC", "DISPLACEMENT", "TEMP");
}
// All is declared, at this point build fields first,
CHKERR m_field.build_fields();
// If 10-node test are on the mesh, use mid nodes to set HO-geometry. Class
// Projection10NodeCoordsOnField
// do the trick.
Projection10NodeCoordsOnField ent_method_material(m_field,
"MESH_NODE_POSITIONS");
CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
if (m_field.check_field("TEMP")) {
if (block_data[it->getMeshsetId()].initTemp != 0) {
MOFEM_LOG_C("ELASTIC", Sev::inform,
"Set block %d temperature to %3.2g\n", it->getMeshsetId(),
block_data[it->getMeshsetId()].initTemp);
Range block_ents;
CHKERR moab.get_entities_by_handle(it->meshset, block_ents, true);
Range vertices;
CHKERR moab.get_connectivity(block_ents, vertices, true);
CHKERR m_field.getInterface<FieldBlas>()->setField(
block_data[it->getMeshsetId()].initTemp, MBVERTEX, vertices,
"TEMP");
}
}
}
// Build database for elements. Actual implementation of element is not need
// here, only elements has to be declared.
// Build adjacencies between elements and field entities
CHKERR m_field.build_adjacencies(bit_level0);
// Register MOFEM DM implementation in PETSc
// Create DM manager
auto dm = createDM(PETSC_COMM_WORLD, "MOFEM");
CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "ELASTIC_PROB", bit_level0);
CHKERR DMSetFromOptions(dm);
CHKERR DMMoFEMSetIsPartitioned(dm, is_partitioned);
// Add elements to DM manager
CHKERR DMMoFEMAddElement(dm, "ELASTIC");
CHKERR DMMoFEMAddElement(dm, "SPRING");
CHKERR DMMoFEMAddElement(dm, "SIMPLE_ROD");
CHKERR DMMoFEMAddElement(dm, "BODY_FORCE");
CHKERR DMMoFEMAddElement(dm, "FLUID_PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
CHKERR DMMoFEMAddElement(dm, "POST_PROC_SKIN");
CHKERR DMSetUp(dm);
// Create matrices & vectors. Note that native PETSc DM interface is used,
// but under the PETSc interface MoFEM implementation is running.
auto D = vectorDuplicate(F);
auto D0 = vectorDuplicate(F);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
// Initialise mass matrix
if (is_calculating_frequency == PETSC_TRUE) {
Mij = matDuplicate(Aij, MAT_DO_NOT_COPY_VALUES);
CHKERR MatSetOption(Mij, MAT_SPD, PETSC_TRUE);
// MatView(Mij, PETSC_VIEWER_STDOUT_SELF);
}
// Assign global matrix/vector contributed by springs
fe_spring_lhs_ptr->ksp_B = Aij;
fe_spring_rhs_ptr->ksp_f = F;
// Assign global matrix/vector contributed by Simple Rod
fe_simple_rod_lhs_ptr->ksp_B = Aij;
fe_simple_rod_rhs_ptr->ksp_f = F;
// Zero vectors and matrices
CHKERR VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecZeroEntries(D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
CHKERR MatZeroEntries(Aij);
// Below particular implementations of finite elements are used to assemble
// problem matrixes and vectors. Implementation of element does not change
// how element is declared.
// Assemble Aij and F. Define Dirichlet bc element, which sets constrains
// to MatrixDouble and the right hand side vector.
// if normally defined boundary conditions are not found,
// DirichletDisplacementBc will try to use DISPLACEMENT blockset. Two
// implementations are available, depending how BC is defined on mesh file.
auto dirichlet_bc_ptr = boost::make_shared<DirichletDisplacementBc>(
m_field, "DISPLACEMENT", Aij, D0, F);
// That sets Dirichlet bc objects that problem is linear, i.e. no newton
// (SNES) solver is run for this problem.
dirichlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
// D0 vector will store initial displacements
CHKERR VecZeroEntries(D0);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
// Run dirichlet_bc, from that on the mesh set values in vector D0. Run
// implementation
// of particular dirichlet_bc.
CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
// Set values from D0 on the field (on the mesh)
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D0, INSERT_VALUES, SCATTER_REVERSE);
// Calculate residual forces as result of applied kinematic constrains. Run
// implementation
// of particular finite element implementation. Look how
// NonlinearElasticElement is implemented,
// in that case. We run NonlinearElasticElement with hook material.
// Calculate right hand side vector
fe_rhs_ptr->snes_f = F;
prism_fe_rhs_ptr->snes_f = F;
MOFEM_LOG("ELASTIC", Sev::inform) << "Assemble external force vector ...";
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_rhs_ptr);
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_rhs_ptr);
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
// Assemble matrix
fe_lhs_ptr->snes_B = Aij;
prism_fe_lhs_ptr->snes_B = Aij;
MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate stiffness matrix ...";
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_lhs_ptr);
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_lhs_ptr);
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
// Assemble springs
CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_lhs_ptr);
CHKERR DMoFEMLoopFiniteElements(dm, "SPRING", fe_spring_rhs_ptr);
// Assemble Simple Rod
CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_lhs_ptr);
CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", fe_simple_rod_rhs_ptr);
if (is_calculating_frequency == PETSC_TRUE) {
// Assemble mass matrix
fe_mass_ptr->snes_B = Mij;
MOFEM_LOG("ELASTIC", Sev::inform) << "Calculate mass matrix ...";
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_mass_ptr);
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
}
// MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
// Assemble pressure and traction forces. Run particular implemented for do
// this, see
// MetaNeumannForces how this is implemented.
boost::ptr_map<std::string, NeumannForcesSurface> neumann_forces;
F, "DISPLACEMENT");
{
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumann_forces.begin();
for (; mit != neumann_forces.end(); mit++) {
CHKERR DMoFEMLoopFiniteElements(dm, mit->first.c_str(),
&mit->second->getLoopFe());
}
}
// Assemble forces applied to nodes, see implementation in NodalForce
boost::ptr_map<std::string, NodalForce> nodal_forces;
CHKERR MetaNodalForces::setOperators(m_field, nodal_forces, F,
"DISPLACEMENT");
{
boost::ptr_map<std::string, NodalForce>::iterator fit =
nodal_forces.begin();
for (; fit != nodal_forces.end(); fit++) {
CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(),
&fit->second->getLoopFe());
}
}
// Assemble edge forces
boost::ptr_map<std::string, EdgeForce> edge_forces;
CHKERR MetaEdgeForces::setOperators(m_field, edge_forces, F,
"DISPLACEMENT");
{
auto fit = edge_forces.begin();
for (; fit != edge_forces.end(); fit++) {
auto &fe = fit->second->getLoopFe();
CHKERR DMoFEMLoopFiniteElements(dm, fit->first.c_str(), &fe);
}
}
// Assemble body forces, implemented in BodyForceConstantField
BodyForceConstantField body_forces_methods(m_field);
m_field, BLOCKSET | BODYFORCESSET, it)) {
CHKERR body_forces_methods.addBlock("DISPLACEMENT", F,
it->getMeshsetId());
}
CHKERR DMoFEMLoopFiniteElements(dm, "BODY_FORCE",
&body_forces_methods.getLoopFe());
// Assemble fluid pressure forces
CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fluid_pressure_fe.getLoopFe(),
false, false);
CHKERR fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
"DISPLACEMENT", F, false, true);
CHKERR DMoFEMLoopFiniteElements(dm, "FLUID_PRESSURE_FE",
&fluid_pressure_fe.getLoopFe());
// Apply kinematic constrain to right hand side vector and matrix
CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
// Matrix View
PetscViewerPushFormat(
PETSC_VIEWER_STDOUT_SELF,
PETSC_VIEWER_ASCII_MATLAB); /// PETSC_VIEWER_ASCII_DENSE,
/// PETSC_VIEWER_ASCII_MATLAB
// MatView(Aij, PETSC_VIEWER_STDOUT_SELF);
// MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
// std::string wait;
// std::cin >> wait;
if (is_calculating_frequency == PETSC_TRUE) {
CHKERR MatAssemblyBegin(Mij, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(Mij, MAT_FINAL_ASSEMBLY);
}
// Set matrix positive defined and symmetric for Cholesky and icc
// pre-conditioner
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
CHKERR VecScale(F, -1);
// Create solver
auto solver = createKSP(PETSC_COMM_WORLD);
CHKERR KSPSetDM(solver, dm);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetOperators(solver, Aij, Aij);
// Setup multi-grid pre-conditioner if set from command line
{
// from PETSc example ex42.c
PetscBool same = PETSC_FALSE;
PC pc;
CHKERR KSPGetPC(solver, &pc);
PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
if (same) {
pc, createPCMGSetUpViaApproxOrdersCtx(dm, Aij, true));
CHKERR PCSetFromOptions(pc);
} else {
// Operators are already set, do not use DM for doing that
CHKERR KSPSetDMActive(solver, PETSC_FALSE);
}
}
CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
CHKERR KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
// Set up solver
CHKERR KSPSetUp(solver);
auto set_post_proc_skin = [&](auto &post_proc_skin) {
CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", post_proc_skin, false,
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);
if (m_field.check_field("TEMP")) {
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) {
// Add operators to the elements, starting with some generic operators
CHKERR post_proc.generateReferenceElementMesh();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", post_proc, true, false, false,
false);
CHKERR post_proc.addFieldValuesPostProc("DISPLACEMENT");
CHKERR post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
if (m_field.check_field("TEMP")) {
CHKERR post_proc.addFieldValuesPostProc("TEMP");
CHKERR post_proc.addFieldValuesGradientPostProc("TEMP");
}
// Add problem specific operator on element to post-process stresses
post_proc.getOpPtrVector().push_back(new PostProcHookStress(
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(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
prism_post_proc.getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
CHKERR prism_post_proc.addFieldValuesPostProc("DISPLACEMENT");
CHKERR prism_post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR prism_post_proc.addFieldValuesGradientPostProc("DISPLACEMENT");
prism_post_proc.getOpPtrVector().push_back(new PostProcHookStress(
m_field, prism_post_proc.postProcMesh, prism_post_proc.mapGaussPts,
"DISPLACEMENT", prism_post_proc.commonData, block_sets_ptr.get()));
};
PostProcFaceOnRefinedMesh post_proc_skin(m_field);
PostProcFatPrismOnRefinedMesh prism_post_proc(m_field);
PostProcEdgeOnRefinedMesh post_proc_edge(m_field);
PostProcVolumeOnRefinedMesh post_proc(m_field);
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;
CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
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) {
field_eval_data = m_field.getInterface<FieldEvaluatorInterface>()
->getData<VolumeElementForcesAndSourcesCore>();
CHKERR m_field.getInterface<FieldEvaluatorInterface>()->buildTree3D(
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;
if (m_field.check_field("TEMP")) {
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("TEMP", scalar_field_ptr));
}
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>("DISPLACEMENT", vector_field_ptr));
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>("DISPLACEMENT", tensor_field_ptr));
}
// Temperature field is defined on the mesh
if (m_field.check_field("TEMP")) {
// Create thermal vector
Vec F_thermal;
CHKERR VecDuplicate(F, &F_thermal);
// Set up implementation for calculation of thermal stress vector. Look
// how thermal stresses and vector is assembled in ThermalStressElement.
CHKERR thermal_stress_elem.setThermalStressRhsOperators(
"DISPLACEMENT", "TEMP", F_thermal);
SeriesRecorder *recorder_ptr;
CHKERR m_field.getInterface(recorder_ptr);
// Read time series and do thermo-elastic analysis, this is when time
// dependent
// temperature problem was run before on the mesh. It means that before
// non-stationary
// problem was solved for temperature and filed "TEMP" is stored for
// subsequent time
// steps in the recorder.
if (recorder_ptr->check_series("THEMP_SERIES")) {
// This is time dependent case, so loop of data series stored by tape
// recorder.
// Loop over time steps
for (_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(recorder_ptr, "THEMP_SERIES",
sit)) {
MOFEM_LOG_C("ELASTIC", Sev::inform, "Process step %d",
sit->get_step_number());
// Load field data for this time step
CHKERR recorder_ptr->load_series_data("THEMP_SERIES",
sit->get_step_number());
CHKERR VecZeroEntries(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
// Calculate the right-hand side vector as result of thermal stresses.
// It MetaNodalForces
// that on "ELASTIC" element data structure the element implementation
// in thermal_stress_elem
// is executed.
dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
// Assemble vector
CHKERR VecAssemblyBegin(F_thermal);
CHKERR VecAssemblyEnd(F_thermal);
// Accumulate ghost dofs
CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
// Calculate norm of vector and print values
PetscReal nrm_F;
CHKERR VecNorm(F, NORM_2, &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);
// check this !!!
CHKERR VecAXPY(F_thermal, 1, F);
// Set dirichlet boundary to thermal stresses vector
dirichlet_bc_ptr->snes_x = D;
dirichlet_bc_ptr->snes_f = F_thermal;
CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
// Solve problem
CHKERR KSPSolve(solver, F_thermal, D);
// Add boundary conditions vector
CHKERR VecAXPY(D, 1., D0);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
// Save data on the mesh
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
// Save data on mesh
CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
if (field_eval_flag) {
->evalFEAtThePoint3D(
field_eval_coords.data(), 1e-12, "ELASTIC_PROB", "ELASTIC",
field_eval_data, m_field.get_comm_rank(),
m_field.get_comm_rank(), nullptr, MF_EXIST, QUIET);
if (scalar_field_ptr->size()) {
auto t_temp = getFTensor0FromVec(*scalar_field_ptr);
MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
<< "Eval point TEMP: " << t_temp;
}
if (vector_field_ptr->size1()) {
auto t_disp = getFTensor1FromMat<3>(*vector_field_ptr);
MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
<< "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);
MOFEM_LOG("ELASTIC_SYNC", Sev::inform)
<< "Eval point DISPLACEMENT_GRAD trace: " << t_disp_grad(i, i);
}
}
// Post-process results
if (is_post_proc_volume == PETSC_TRUE) {
MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
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 ...";
CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_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 {
// This is a case when stationary problem for temperature was solved.
CHKERR VecZeroEntries(F_thermal);
CHKERR VecGhostUpdateBegin(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F_thermal, INSERT_VALUES, SCATTER_FORWARD);
// Calculate the right-hand side vector with thermal stresses
dm, "ELASTIC", &thermal_stress_elem.getLoopThermalStressRhs());
// Assemble vector
CHKERR VecAssemblyBegin(F_thermal);
CHKERR VecAssemblyEnd(F_thermal);
// Accumulate ghost dofs
CHKERR VecGhostUpdateBegin(F_thermal, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F_thermal, ADD_VALUES, SCATTER_REVERSE);
// Calculate norm
PetscReal nrm_F;
CHKERR VecNorm(F, NORM_2, &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);
// Add thermal stress vector and other forces vector
CHKERR VecScale(F_thermal, -1);
CHKERR VecAXPY(F_thermal, 1, F);
// Apply kinetic boundary conditions
dirichlet_bc_ptr->snes_x = D;
dirichlet_bc_ptr->snes_f = F_thermal;
CHKERR DMoFEMPostProcessFiniteElements(dm, dirichlet_bc_ptr.get());
// Solve problem
CHKERR KSPSolve(solver, F_thermal, D);
CHKERR VecAXPY(D, 1., D0);
// Update ghost values for solution vector
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
// Save data on mesh
if (is_post_proc_volume == PETSC_TRUE) {
MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file ...";
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
// Save results to file
if (!test_nb)
CHKERR post_proc.writeFile("out.h5m");
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
}
MOFEM_LOG("ELASTIC", Sev::inform) << "Write output file skin ...";
CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
if (!test_nb)
CHKERR post_proc_skin.writeFile("out_skin.h5m");
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
}
// Destroy vector, no needed any more
CHKERR VecDestroy(&F_thermal);
} else {
// Elastic analysis no temperature field
// VecView(F, PETSC_VIEWER_STDOUT_WORLD);
// Solve for vector D
CHKERR KSPSolve(solver, F, D);
// VecView(D, PETSC_VIEWER_STDOUT_WORLD);
// cerr << F;
// Add kinetic boundary conditions
CHKERR VecAXPY(D, 1., D0);
// Update ghost values
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
// Save data from vector on mesh
CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
// Post-process results
MOFEM_LOG("ELASTIC", Sev::inform) << "Post-process start ...";
if (is_post_proc_volume == PETSC_TRUE) {
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &post_proc);
}
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", &prism_post_proc);
CHKERR DMoFEMLoopFiniteElements(dm, "SIMPLE_ROD", &post_proc_edge);
CHKERR DMoFEMLoopFiniteElements(dm, "POST_PROC_SKIN", &post_proc_skin);
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
// Write mesh in parallel (using h5m MOAB format, writing is in parallel)
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");
MOFEM_LOG("ELASTIC", Sev::inform) << "done";
}
if (is_calculating_frequency == PETSC_TRUE) {
// Calculate mode mass, m = u^T * M * u
Vec u1;
VecDuplicate(D, &u1);
CHKERR MatMult(Mij, D, u1);
double mode_mass;
CHKERR VecDot(u1, D, &mode_mass);
MOFEM_LOG_C("ELASTIC", Sev::inform, "Mode mass %6.4e\n", mode_mass);
Vec v1;
VecDuplicate(D, &v1);
CHKERR MatMult(Aij, D, v1);
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);
}
// Calculate elastic energy
auto calculate_strain_energy = [&]() {
CHKERR HookeElement::calculateEnergy(dm, block_sets_ptr, "DISPLACEMENT",
"MESH_NODE_POSITIONS", false, true,
v_energy);
// Print elastic 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)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged!");
break;
case 2:
if (fabs(energy - 5.6475e-03) > 1e-4)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged!");
break;
case 3:
if (fabs(energy - 7.4679e-03) > 1e-4)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged!");
break;
case 4:
if (fabs(energy - 2.4992e+00) > 1e-3)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged!");
break;
// FIXME: Here are missing regersion tests
case 8: {
double min;
CHKERR VecMin(D, PETSC_NULL, &min);
constexpr double expected_val = 0.10001;
if (fabs(min + expected_val) > 1e-10)
SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged! %3.4e != %3.4e", min, expected_val);
} break;
case 9: {
if (fabs(energy - 4.7416e-04) > 1e-8)
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"atom test diverged!");
}
default:
break;
}
};
CHKERR calculate_strain_energy();
MPI_Comm_free(&moab_comm_world);
}
return 0;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MetaNeumannForces::addNeumannBCElements
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MetaNodalForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:674
MoFEM::SeriesRecorder::check_series
virtual bool check_series(const std::string &name) const
check if series is in database
Definition: SeriesRecorder.cpp:301
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::createPCMGSetUpViaApproxOrdersCtx
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
Definition: PCMGSetUpViaApproxOrders.cpp:630
MoFEM::SeriesRecorder
Definition: SeriesRecorder.hpp:25
PrismFE::getRuleTrianglesOnly
int getRuleTrianglesOnly(int order)
Definition: prism_elements_from_surface.cpp:436
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
PostProcFatPrismOnRefinedMesh
Postprocess on prism.
Definition: PostProcOnRefMesh.hpp:973
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
ClipperLib::pi
static const double pi
Definition: clipper.cpp:89
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
MoFEM::CoreInterface::add_ents_to_field_by_dim
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3145
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MetaNodalForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore
FatPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: FatPrismElementForcesAndSourcesCore.cpp:11
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:304
MoFEM::CoreInterface::add_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
MetaSimpleRodElement::setSimpleRodOperators
static MoFEMErrorCode setSimpleRodOperators(MoFEM::Interface &m_field, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_lhs_ptr, boost::shared_ptr< EdgeElementForcesAndSourcesCore > fe_simple_rod_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Implementation of SimpleRod element. Set operators to calculate LHS and RHS.
Definition: SimpleRodElement.cpp:298
BodyForceConstantField
Body forces elements.
Definition: BodyForce.hpp:12
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
PrismFE::getRuleThroughThickness
int getRuleThroughThickness(int order)
Definition: prism_elements_from_surface.cpp:437
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3172
MetaNeumannForces::setMomentumFluxOperators
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Definition: SurfacePressure.cpp:2069
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
EdgeEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Definition: elasticity.cpp:94
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
Hooke.hpp
Implementation of linear elastic material.
BODYFORCESSET
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:175
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition: PCMGSetUpViaApproxOrders.cpp:302
FluidPressure
Fluid pressure forces.
Definition: FluidPressure.hpp:20
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
MoFEM::SeriesRecorder::load_series_data
virtual MoFEMErrorCode load_series_data(const std::string &serie_name, const int step_number)
Definition: SeriesRecorder.cpp:309
_IT_SERIES_STEPS_BY_NAME_FOR_LOOP_
#define _IT_SERIES_STEPS_BY_NAME_FOR_LOOP_(RECORDER, NAME, IT)
loop over recorded series step
Definition: SeriesRecorder.hpp:205
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PrismFE
Definition: prism_elements_from_surface.cpp:62
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
PostProcEdgeOnRefinedMesh
Postprocess on edge.
Definition: PostProcOnRefMesh.hpp:1148
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:234
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
MoFEM::DMMoFEMCreateMoFEM
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.
Definition: DMMoFEM.cpp:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:325
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FTensor::Index< 'i', 3 >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
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
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
MetaEdgeForces::addElement
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
Range
ConvectiveMassElement::setBlocks
MoFEMErrorCode setBlocks()
Definition: ConvectiveMassElement.cpp:1761
MetaSpringBC::setSpringOperators
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
Definition: SpringElement.cpp:1178
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:290
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MetaSimpleRodElement::addSimpleRodElements
static MoFEMErrorCode addSimpleRodElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare SimpleRod element.
Definition: SimpleRodElement.cpp:270
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
addElasticElement
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
Definition: HookeElement.cpp:533
main
int main(int argc, char *argv[])
Definition: elasticity.cpp:96
setOperators
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
BlockOptionData
Definition: elasticity.cpp:62
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PCMGSetUpViaApproxOrders
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
Definition: PCMGSetUpViaApproxOrders.cpp:634
MoFEM::CoreInterface::set_field_order
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.
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Vec >
help
static char help[]
Definition: elasticity.cpp:57
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
ThermalStressElement
Implentation of thermal stress element.
Definition: ThermalStressElement.hpp:15
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::LogManager
Log manager is used to build and partition problems.
Definition: LogManager.hpp:26
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::CoreInterface::add_field
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.
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: elasticity.cpp:74
F
@ F
Definition: free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MetaEdgeForces::setOperators
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97