#include <BasicFiniteElements.hpp>
#ifndef WITH_ADOL_C
#error "MoFEM need to be compiled with ADOL-C"
#endif
static char help[] =
"...\n\n";
double GenericMaterial::ePsilon0 = 0;
double GenericMaterial::ePsilon1 = 0;
double GenericMaterial::scaleZ = 1;
map<std::string, CommonMaterialData::RegisterHook>
RegisterMaterials::mapOfRegistredMaterials;
int main(
int argc,
char *argv[]) {
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n";
std::ofstream file(
param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file.close();
}
}
PetscBool test_mat = PETSC_FALSE;
PETSC_NULL);
if (test_mat == PETSC_TRUE) {
CommonMaterialData data;
MaterialVanGenuchten van_genuchten(data);
van_genuchten.printMatParameters(-1, "testing");
van_genuchten.printTheta(-1e-16, -1e4, -1e-12, "hTheta");
van_genuchten.printKappa(-1e-16, -1, -1e-12, "hK");
van_genuchten.printC(-1e-16, -1, -1e-12, "hC");
return 0;
}
try {
PetscBool flg_mesh_file_name = PETSC_TRUE;
PetscBool flg_conf_file_name = PETSC_TRUE;
char conf_file_name[255];
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Unsaturated flow options",
"none");
ierr = PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
ierr = PetscOptionsString(
"-configure", "material and bc configuration file name", "",
"unsaturated.cfg", conf_file_name, 255, &flg_conf_file_name);
ierr = PetscOptionsInt(
"-my_order",
"default approximation order",
"",
ierr = PetscOptionsEnd();
if (flg_mesh_file_name != PETSC_TRUE) {
"File name not set (e.g. -my_name mesh.h5m)");
}
if (flg_conf_file_name != PETSC_TRUE) {
"File name not set (e.g. -config unsaturated.cfg)");
}
auto moab_comm_wrap =
boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
if (pcomm == NULL)
pcomm =
new ParallelComm(&moab, moab_comm_wrap->get_comm());
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
MeshsetsManager *meshsets_manager_ptr;
CHKERR meshsets_manager_ptr->setMeshsetFromFile();
PetscPrintf(PETSC_COMM_WORLD,
"Read meshsets and added meshsets from bc.cfg\n");
PetscPrintf(PETSC_COMM_WORLD, "%s",
static_cast<std::ostringstream &>(
std::ostringstream().seekp(0) << *it << endl)
.str()
.c_str());
}
UnsaturatedFlowElement uf(m_field);
ifstream ini_file(conf_file_name);
po::variables_map vm;
po::options_description config_file_options;
Range domain_ents, bc_boundary_ents;
map<int, CommonMaterialData> material_blocks;
map<int, double> head_blocks;
map<int, double> flux_blocks;
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
const int block_id = it->getMeshsetId();
std::string block_name =
"mat_block_" + boost::lexical_cast<std::string>(block_id);
material_blocks[block_id].blockId = block_id;
material_blocks[block_id].addOptions(config_file_options, block_name);
}
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
const int block_id = it->getMeshsetId();
std::string block_name =
"head_block_" + boost::lexical_cast<std::string>(block_id);
config_file_options.add_options()(
(block_name + ".head").c_str(),
po::value<double>(&head_blocks[block_id])->default_value(0));
}
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
const int block_id = it->getMeshsetId();
std::string block_name =
"flux_block_" + boost::lexical_cast<std::string>(block_id);
config_file_options.add_options()(
(block_name + ".flux").c_str(),
po::value<double>(&head_blocks[block_id])->default_value(0));
}
po::parsed_options parsed =
parse_config_file(ini_file, config_file_options, true);
store(parsed, vm);
po::notify(vm);
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
const int block_id = it->getMeshsetId();
uf.dMatMap[block_id] = RegisterMaterials::mapOfRegistredMaterials.at(
material_blocks[block_id].matName)(material_blocks.at(block_id));
if (!uf.dMatMap.at(block_id)) {
"Material block not set");
}
it->meshset, MBTET, uf.dMatMap.at(block_id)->tEts, true);
domain_ents.merge(uf.dMatMap.at(block_id)->tEts);
uf.dMatMap.at(block_id)->printMatParameters(block_id, "Read material");
}
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++) {
"** WARNING Unrecognized option %s\n", vit->c_str());
}
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
const int block_id = it->getMeshsetId();
uf.bcValueMap[block_id] =
boost::shared_ptr<UnsaturatedFlowElement::BcData>(
new UnsaturatedFlowElement::BcData());
std::vector<double> attributes;
CHKERR it->getAttributes(attributes);
if (attributes.size() != 1)
"Wrong number of head parameters %d", attributes.size());
uf.bcValueMap[block_id]->fixValue = attributes[0];
std::string block_name =
"head_block_" + boost::lexical_cast<std::string>(block_id);
if (vm.count((block_name) + ".head")) {
uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
}
it->meshset, MBTRI, uf.bcValueMap[block_id]->eNts, true);
bc_boundary_ents.merge(uf.bcValueMap[block_id]->eNts);
}
int max_flux_id = 0;
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
const int block_id = it->getMeshsetId();
uf.bcFluxMap[block_id] =
boost::shared_ptr<UnsaturatedFlowElement::BcData>(
new UnsaturatedFlowElement::BcData());
std::vector<double> attributes;
CHKERR it->getAttributes(attributes);
uf.bcFluxMap[block_id]->fixValue = attributes[0];
std::string block_name =
"flux_block_" + boost::lexical_cast<std::string>(block_id);
if (vm.count((block_name) + ".flux")) {
uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
}
it->meshset, MBTRI, uf.bcFluxMap[block_id]->eNts, true);
bc_boundary_ents.merge(uf.bcFluxMap[block_id]->eNts);
max_flux_id = max_flux_id > block_id ? max_flux_id : block_id + 1;
}
auto get_zero_flux_entities = [&m_field,
&pcomm](const auto &domain_ents,
const auto &bc_boundary_ents) {
Range zero_flux_ents;
Range domain_skin;
CHKERR skin.find_skin(0, domain_ents,
false, domain_skin);
CHKERR pcomm->filter_pstatus(domain_skin,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &zero_flux_ents);
zero_flux_ents = subtract(zero_flux_ents, bc_boundary_ents);
return zero_flux_ents;
};
CHKERR uf.addFiniteElements(
"FLUXES",
"VALUES");
get_zero_flux_entities(domain_ents, bc_boundary_ents));
CHKERR uf.setFiniteElements();
CHKERR uf.calculateEssentialBc();
CHKERR uf.calculateInitialPc();
}
return 0;
}
const std::string default_options
Mix implementation of transport element.
Mix implementation of transport element.
Mix implementation of transport element.
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets 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.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
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
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
DeprecatedCoreInterface Interface
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() 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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
int main(int argc, char *argv[])