17 #error "MoFEM need to be compiled with ADOL-C"
20 using namespace MoFEM;
22 static char help[] =
"...\n\n";
24 double GenericMaterial::ePsilon0 = 0;
25 double GenericMaterial::ePsilon1 = 0;
26 double GenericMaterial::scaleZ = 1;
28 map<std::string, CommonMaterialData::RegisterHook>
29 RegisterMaterials::mapOfRegistredMaterials;
31 int main(
int argc,
char *argv[]) {
33 const string default_options =
"-ksp_type fgmres\n"
35 "-pc_factor_mat_solver_type mumps \n"
36 "-mat_mumps_icntl_20 0 \n";
38 string param_file =
"param_file.petsc";
39 if (!
static_cast<bool>(ifstream(param_file))) {
40 std::ofstream file(param_file.c_str(), std::ios::ate);
42 file << default_options;
54 PetscBool test_mat = PETSC_FALSE;
58 if (test_mat == PETSC_TRUE) {
63 van_genuchten.
printTheta(-1e-16, -1e4, -1e-12,
"hTheta");
64 van_genuchten.
printKappa(-1e-16, -1, -1e-12,
"hK");
65 van_genuchten.
printC(-1e-16, -1, -1e-12,
"hC");
76 PetscBool flg_mesh_file_name = PETSC_TRUE;
78 PetscBool flg_conf_file_name = PETSC_TRUE;
79 char conf_file_name[255];
83 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Unsaturated flow options",
86 ierr = PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
89 ierr = PetscOptionsString(
90 "-configure",
"material and bc configuration file name",
"",
91 "unsaturated.cfg", conf_file_name, 255, &flg_conf_file_name);
93 ierr = PetscOptionsInt(
"-my_order",
"default approximation order",
"",
96 ierr = PetscOptionsEnd();
99 if (flg_mesh_file_name != PETSC_TRUE) {
101 "File name not set (e.g. -my_name mesh.h5m)");
103 if (flg_conf_file_name != PETSC_TRUE) {
105 "File name not set (e.g. -config unsaturated.cfg)");
110 auto moab_comm_wrap =
111 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
114 new ParallelComm(&moab, moab_comm_wrap->get_comm());
117 option =
"PARALLEL=READ_PART;"
118 "PARALLEL_RESOLVE_SHARED_ENTS;"
119 "PARTITION=PARALLEL_PARTITION;";
131 PetscPrintf(PETSC_COMM_WORLD,
132 "Read meshsets and added meshsets from bc.cfg\n");
134 PetscPrintf(PETSC_COMM_WORLD,
"%s",
135 static_cast<std::ostringstream &
>(
136 std::ostringstream().seekp(0) << *it << endl)
147 ifstream ini_file(conf_file_name);
148 po::variables_map vm;
149 po::options_description config_file_options;
150 Range domain_ents, bc_boundary_ents;
152 map<int, CommonMaterialData> material_blocks;
153 map<int, double> head_blocks;
154 map<int, double> flux_blocks;
157 if (it->getName().compare(0, 4,
"SOIL") != 0)
160 const int block_id = it->getMeshsetId();
161 std::string block_name =
162 "mat_block_" + boost::lexical_cast<std::string>(block_id);
163 material_blocks[block_id].blockId = block_id;
164 material_blocks[block_id].addOptions(config_file_options, block_name);
167 if (it->getName().compare(0, 4,
"HEAD") != 0)
170 const int block_id = it->getMeshsetId();
171 std::string block_name =
172 "head_block_" + boost::lexical_cast<std::string>(block_id);
173 config_file_options.add_options()(
174 (block_name +
".head").c_str(),
175 po::value<double>(&head_blocks[block_id])->default_value(0));
178 if (it->getName().compare(0, 4,
"FLUX") != 0)
181 const int block_id = it->getMeshsetId();
182 std::string block_name =
183 "flux_block_" + boost::lexical_cast<std::string>(block_id);
184 config_file_options.add_options()(
185 (block_name +
".flux").c_str(),
186 po::value<double>(&head_blocks[block_id])->default_value(0));
188 po::parsed_options parsed =
189 parse_config_file(ini_file, config_file_options,
true);
193 if (it->getName().compare(0, 4,
"SOIL") != 0)
195 const int block_id = it->getMeshsetId();
196 uf.dMatMap[block_id] = RegisterMaterials::mapOfRegistredMaterials.at(
197 material_blocks[block_id].matName)(material_blocks.at(block_id));
198 if (!uf.dMatMap.at(block_id)) {
200 "Material block not set");
204 it->meshset, MBTET, uf.dMatMap.at(block_id)->tEts,
true);
205 domain_ents.merge(uf.dMatMap.at(block_id)->tEts);
206 uf.dMatMap.at(block_id)->printMatParameters(block_id,
"Read material");
208 std::vector<std::string> additional_parameters;
209 additional_parameters =
210 collect_unrecognized(parsed.options, po::include_positional);
211 for (std::vector<std::string>::iterator vit = additional_parameters.begin();
212 vit != additional_parameters.end(); vit++) {
214 "** WARNING Unrecognized option %s\n", vit->c_str());
218 if (it->getName().compare(0, 4,
"HEAD") != 0)
221 const int block_id = it->getMeshsetId();
223 uf.bcValueMap[block_id] =
224 boost::shared_ptr<UnsaturatedFlowElement::BcData>(
227 std::vector<double> attributes;
228 CHKERR it->getAttributes(attributes);
229 if (attributes.size() != 1)
231 "Wrong number of head parameters %d", attributes.size());
232 uf.bcValueMap[block_id]->fixValue = attributes[0];
233 std::string block_name =
234 "head_block_" + boost::lexical_cast<std::string>(block_id);
235 if (vm.count((block_name) +
".head")) {
236 uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
242 it->meshset, MBTRI, uf.bcValueMap[block_id]->eNts,
true);
243 bc_boundary_ents.merge(uf.bcValueMap[block_id]->eNts);
249 if (it->getName().compare(0, 4,
"FLUX") != 0)
252 const int block_id = it->getMeshsetId();
254 uf.bcFluxMap[block_id] =
255 boost::shared_ptr<UnsaturatedFlowElement::BcData>(
258 std::vector<double> attributes;
259 CHKERR it->getAttributes(attributes);
261 uf.bcFluxMap[block_id]->fixValue = attributes[0];
262 std::string block_name =
263 "flux_block_" + boost::lexical_cast<std::string>(block_id);
264 if (vm.count((block_name) +
".flux")) {
265 uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
269 it->meshset, MBTRI, uf.bcFluxMap[block_id]->eNts,
true);
270 bc_boundary_ents.merge(uf.bcFluxMap[block_id]->eNts);
271 max_flux_id = max_flux_id > block_id ? max_flux_id : block_id + 1;
275 auto get_zero_flux_entities = [&m_field,
276 &pcomm](
const auto &domain_ents,
277 const auto &bc_boundary_ents) {
278 Range zero_flux_ents;
281 CHKERR skin.find_skin(0, domain_ents,
false, domain_skin);
283 CHKERR pcomm->filter_pstatus(domain_skin,
284 PSTATUS_SHARED | PSTATUS_MULTISHARED,
285 PSTATUS_NOT, -1, &zero_flux_ents);
286 zero_flux_ents = subtract(zero_flux_ents, bc_boundary_ents);
287 return zero_flux_ents;
291 CHKERR uf.addFiniteElements(
"FLUXES",
"VALUES");
293 get_zero_flux_entities(domain_ents, bc_boundary_ents));
294 CHKERR uf.createMatrices();
295 CHKERR uf.setFiniteElements();
296 CHKERR uf.calculateEssentialBc();
297 CHKERR uf.calculateInitialPc();
299 CHKERR uf.destroyMatrices();