#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();
}
}
CHKERR DMRegister_MoFEM(
"DMMOFEM");
PetscBool test_mat = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL,
"",
"-test_mat", &test_mat,
PETSC_NULL);
if (test_mat == PETSC_TRUE) {
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 {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
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;";
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());
}
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>(
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>(
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) {
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.
static PetscErrorCode ierr
#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.
#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.
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
void printMatParameters(const int id, const std::string &prefix) const
void printTheta(const double b, const double e, double s, const std::string &prefix)
void printC(const double b, const double e, double s, const std::string &prefix)
void printKappa(const double b, const double e, double s, const std::string &prefix)
Class storing information about boundary condition.
Implementation of operators, problem and finite elements for unsaturated flow.
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.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.