36#include <boost/numeric/ublas/vector_proxy.hpp>
37#include <boost/numeric/ublas/matrix.hpp>
38#include <boost/numeric/ublas/matrix_proxy.hpp>
39#include <boost/numeric/ublas/vector.hpp>
40#include <adolc/adolc.h>
45#include <boost/program_options.hpp>
47namespace po = boost::program_options;
62#include <boost/iostreams/tee.hpp>
63#include <boost/iostreams/stream.hpp>
66namespace bio = boost::iostreams;
72static char help[] =
"...\n\n";
81int main(
int argc,
char *argv[]) {
87 moab::Core mb_instance;
88 moab::Interface& moab = mb_instance;
90 PetscBool flg_gel_config,flg_file;
92 char gel_config_file[255];
94 PetscBool is_partitioned = PETSC_FALSE;
95 PetscBool is_atom_test = PETSC_FALSE;
99 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none"); CHKERRQ(
ierr);
100 ierr = PetscOptionsString(
108 ierr = PetscOptionsInt(
110 "default approximation order",
117 ierr = PetscOptionsBool(
118 "-my_is_partitioned",
119 "set if mesh is partitioned (this result that each process keeps only part of the mes",
"",
120 PETSC_FALSE,&is_partitioned,PETSC_NULL
122 ierr = PetscOptionsString(
124 "gel configuration file name",
"",
125 "gel_config.in",gel_config_file,255,&flg_gel_config
128 ierr = PetscOptionsBool(
130 "is used with testing, exit with error when diverged",
"",
131 PETSC_FALSE,&is_atom_test,PETSC_NULL
134 ierr = PetscOptionsEnd(); CHKERRQ(
ierr);
137 if(flg_file != PETSC_TRUE) {
138 SETERRQ(PETSC_COMM_SELF,1,
"*** ERROR -my_file (MESH FILE NEEDED)");
141 if(is_partitioned == PETSC_TRUE) {
145 "PARALLEL=BCAST_DELETE;"
146 "PARALLEL_RESOLVE_SHARED_ENTS;"
147 "PARTITION=PARALLEL_PARTITION;";
154 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
155 if(pcomm == NULL) pcomm =
new ParallelComm(&moab,PETSC_COMM_WORLD);
169 map<int,ThermalElement::FluxData> set_of_solvent_fluxes;
176 bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
236 map<int,BlockData> block_data;
239 ifstream ini_file(gel_config_file);
240 po::variables_map vm;
241 po::options_description config_file_options;
245 rval = moab.tag_get_handle(
246 "BLOCK_ID",1,MB_TYPE_INTEGER,th_block_id,MB_TAG_CREAT|MB_TAG_SPARSE,&def_block
251 const string &name = it->getName();
252 if(name.compare(0,3,
"GEL")!=0)
continue;
254 int block_id = it->getMeshsetId();
256 ostringstream str_order;
257 str_order <<
"block_" << block_id <<
".oRder";
258 config_file_options.add_options()
259 (str_order.str().c_str(),po::value<int>(
260 &block_data[block_id].oRder
261 )->default_value(
order));
263 rval = moab.get_entities_by_handle(
264 it->getMeshset(),block_ents,true
266 material_data[block_id].tEts = block_ents.subset_by_type(MBTET);
267 std::vector<int> block_id_vec(material_data[block_id].tEts.size(),it->getMeshsetId());
269 rval = moab.tag_set_data(
270 th_block_id,material_data[block_id].tEts,&*block_id_vec.begin()
273 ostringstream str_g_alpha;
274 str_g_alpha <<
"block_" << block_id <<
".gAlpha";
275 config_file_options.add_options()
276 (str_g_alpha.str().c_str(),po::value<double>(&material_data[block_id].gAlpha)->default_value(1));
277 ostringstream str_v_alpha;
278 str_v_alpha <<
"block_" << block_id <<
".vAlpha";
279 config_file_options.add_options()
280 (str_v_alpha.str().c_str(),po::value<double>(&material_data[block_id].vAlpha)->default_value(0));
281 ostringstream str_g_beta;
282 str_g_beta <<
"block_" << block_id <<
".gBeta";
283 config_file_options.add_options()
284 (str_g_beta.str().c_str(),po::value<double>(&material_data[block_id].gBeta)->default_value(1));
285 ostringstream str_v_beta;
286 str_v_beta <<
"block_" << block_id <<
".vBeta";
287 config_file_options.add_options()
288 (str_v_beta.str().c_str(),po::value<double>(&material_data[block_id].vBeta)->default_value(0));
289 ostringstream str_g_beta_hat;
290 str_g_beta_hat <<
"block_" << block_id <<
".gBetaHat";
291 config_file_options.add_options()
292 (str_g_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].gBetaHat)->default_value(1));
293 ostringstream str_v_beta_hat;
294 str_v_beta_hat <<
"block_" << block_id <<
".vBetaHat";
295 config_file_options.add_options()
296 (str_v_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].vBetaHat)->default_value(0));
297 ostringstream str_omega;
298 str_omega <<
"block_" << block_id <<
".oMega";
299 config_file_options.add_options()
300 (str_omega.str().c_str(),po::value<double>(&material_data[block_id].oMega)->default_value(1));
301 ostringstream str_viscosity;
302 str_viscosity <<
"block_" << block_id <<
".vIscosity";
303 config_file_options.add_options()
304 (str_viscosity.str().c_str(),po::value<double>(&material_data[block_id].vIscosity)->default_value(1));
305 ostringstream str_permeability;
306 str_permeability <<
"block_" << block_id <<
".pErmeability";
307 config_file_options.add_options()
308 (str_permeability.str().c_str(),po::value<double>(&material_data[block_id].pErmeability)->default_value(1));
309 ostringstream str_mu0;
310 str_mu0 <<
"block_" << block_id <<
".mU0";
311 config_file_options.add_options()
312 (str_mu0.str().c_str(),po::value<double>(&material_data[block_id].mU0)->default_value(1));
316 po::parsed_options parsed = parse_config_file(ini_file,config_file_options,
true);
321 const string &name = it->getName();
322 if(name.compare(0,3,
"GEL")!=0)
continue;
324 if(block_data[it->getMeshsetId()].oRder == -1)
continue;
325 if(block_data[it->getMeshsetId()].oRder ==
order)
continue;
326 PetscPrintf(PETSC_COMM_WORLD,
"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
329 Range ents_to_set_order;
330 rval = moab.get_adjacencies(block_ents,3,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
331 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
332 rval = moab.get_adjacencies(block_ents,2,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
333 rval = moab.get_adjacencies(block_ents,1,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
335 ierr = m_field.
set_field_order(ents_to_set_order,
"SPATIAL_POSITION",block_data[it->getMeshsetId()].oRder); CHKERRQ(
ierr);
336 ierr = m_field.
set_field_order(ents_to_set_order,
"CHEMICAL_LOAD",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
337 ierr = m_field.
set_field_order(ents_to_set_order,
"HAT_EPS",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
338 ierr = m_field.
set_field_order(ents_to_set_order,
"SPATIAL_POSITION_DOT",block_data[it->getMeshsetId()].oRder); CHKERRQ(
ierr);
339 ierr = m_field.
set_field_order(ents_to_set_order,
"CHEMICAL_LOAD_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
340 ierr = m_field.
set_field_order(ents_to_set_order,
"HAT_EPS_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
342 vector<string> additional_parameters;
343 additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
344 for(vector<string>::iterator vit = additional_parameters.begin();
345 vit!=additional_parameters.end();vit++) {
346 ierr = PetscPrintf(PETSC_COMM_WORLD,
"** WARNING Unrecognised option %s\n",vit->c_str()); CHKERRQ(
ierr);
348 }
catch (
const std::exception& ex) {
350 ss << ex.what() << endl;
358 ierr = PetscPrintf(PETSC_COMM_WORLD,
"** Block Id %d\n",mit->first); CHKERRQ(
ierr);
359 PetscPrintf(PETSC_COMM_WORLD,
"vAlpha = %5.4g\n",mit->second.vAlpha);
360 PetscPrintf(PETSC_COMM_WORLD,
"gAlpha = %5.4g\n",mit->second.gAlpha);
361 PetscPrintf(PETSC_COMM_WORLD,
"vBeta = %5.4g\n",mit->second.vBeta);
362 PetscPrintf(PETSC_COMM_WORLD,
"gBeta = %5.4g\n",mit->second.gBeta);
363 PetscPrintf(PETSC_COMM_WORLD,
"vBetaHat = %5.4g\n",mit->second.vBetaHat);
364 PetscPrintf(PETSC_COMM_WORLD,
"gBetaHat = %5.4g\n",mit->second.gBetaHat);
365 PetscPrintf(PETSC_COMM_WORLD,
"pErmeability = %5.4g\n",mit->second.pErmeability);
366 PetscPrintf(PETSC_COMM_WORLD,
"vIscosity = %5.4g\n",mit->second.vIscosity);
367 PetscPrintf(PETSC_COMM_WORLD,
"oMega = %5.4g\n",mit->second.oMega);
368 PetscPrintf(PETSC_COMM_WORLD,
"mU0 = %5.4g\n",mit->second.mU0);
375 if (!check_if_spatial_field_exist) {
377 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material); CHKERRQ(
ierr);
383 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
384 mit!=material_data.end();mit++
389 ierr = m_field.set_field(mit->second.mU0,MBVERTEX,vertices,
"CHEMICAL_LOAD"); CHKERRQ(
ierr);
415 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
416 mit!=material_data.end();
423 ierr = MetaNeummanForces::addNeumannBCElements(m_field,
"SPATIAL_POSITION"); CHKERRQ(
ierr);
439 if(it->getName().compare(0,18,
"FLUX_CHEMICAL_LOAD") == 0) {
441 ierr = it->getAttributes(data); CHKERRQ(
ierr);
443 SETERRQ(PETSC_COMM_SELF,1,
"Data inconsistency");
448 strcpy(set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.name,
"HeatFlu");
449 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
450 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
453 it->meshset,MBTRI,set_of_solvent_fluxes[it->getMeshsetId()].tRis,
true
456 set_of_solvent_fluxes[it->getMeshsetId()].tRis,MBTRI,
"CHEMICAL_LOAD_FLUX_FE"
488 common_data.spatialPositionNameDot =
"SPATIAL_POSITION_DOT";
489 common_data.muName =
"CHEMICAL_LOAD";
490 common_data.muNameDot =
"CHEMICAL_LOAD_DOT";
491 common_data.strainHatName =
"HAT_EPS";
492 common_data.strainHatNameDot =
"HAT_EPS_DOT";
497 for(
int ss = 0;ss<2;ss++) {
498 ierr =
addHOOpsVol(
"MESH_NODE_POSITIONS", *fe_ptr[ss],
true,
true,
false,
570 DMType dm_name =
"DMGEL";
573 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(
ierr);
574 ierr = DMSetType(dm,dm_name);CHKERRQ(
ierr);
576 ierr = DMSetFromOptions(dm); CHKERRQ(
ierr);
596 m_field,
"SPATIAL_POSITION",
A,
T,
F
598 spatial_position_bc.methodsOp.push_back(
new TimeForceScale(
"-my_displacements_history",
false));
600 m_field,
"CHEMICAL_LOAD",
"CHEMICAL_LOAD",
A,
T,
F
602 concentration_bc.methodsOp.push_back(
new TimeForceScale(
"-my_chemical_load_history",
false));
605 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
606 boost::ptr_map<string,NodalForce> nodal_forces;
607 boost::ptr_map<string,EdgeForce> edge_forces;
610 ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,
"SPATIAL_POSITION"); CHKERRQ(
ierr);
616 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
617 mit!=neumann_forces.end();mit++
619 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
622 boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
623 mit!=nodal_forces.end();mit++
625 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
628 boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
629 mit!=edge_forces.end();mit++
631 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
639 map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
640 for(;sit!=set_of_solvent_fluxes.end();sit++) {
656 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
657 for(;mit!=neumann_forces.end();mit++) {
662 boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
663 for(;fit!=nodal_forces.end();fit++) {
668 boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
669 for(;fit!=edge_forces.end();fit++) {
689 ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(
ierr);
690 ierr = TSSetType(ts,TSBEULER); CHKERRQ(
ierr);
703 ierr = TSSetIFunction(ts,
F,PETSC_NULL,PETSC_NULL); CHKERRQ(
ierr);
704 ierr = TSSetIJacobian(ts,
A,
A,PETSC_NULL,PETSC_NULL); CHKERRQ(
ierr);
713 ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
714 ierr = TSMonitorSet(ts,f_TSMonitorSet,
ts_ctx,PETSC_NULL); CHKERRQ(
ierr);
718 ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(
ierr);
719 ierr = TSSetFromOptions(ts); CHKERRQ(
ierr);
720 ierr = TSSetDM(ts,dm); CHKERRQ(
ierr);
721 #if PETSC_VERSION_GE(3,7,0)
722 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(
ierr);
725 ierr = TSGetTime(ts,&ftime); CHKERRQ(
ierr);
726 PetscInt steps,snesfails,rejects,nonlinits,linits;
727 ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(
ierr);
728 ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(
ierr);
729 ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(
ierr);
730 ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(
ierr);
731 ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(
ierr);
732 PetscPrintf(PETSC_COMM_WORLD,
733 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
734 steps,rejects,snesfails,ftime,nonlinits,linits
741 ierr = PetscPrintf(PETSC_COMM_WORLD,
"sum = %9.8e\n",sum); CHKERRQ(
ierr);
743 ierr = VecNorm(
T,NORM_2,&fnorm); CHKERRQ(
ierr);
744 ierr = PetscPrintf(PETSC_COMM_WORLD,
"fnorm = %9.8e\n",fnorm); CHKERRQ(
ierr);
746 if (fabs(sum - 6.82084809e+01) > 1e-6) {
749 if (fabs(fnorm - 3.66184607e+01) > 1e-6) {
756 ierr = DMDestroy(&dm); CHKERRQ(
ierr);
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProc
DEPRECATED typedef DirichletSpatialPositionsBc SpatialPositionsBCFEMethodPreAndPostProc
Implementation of Gel finite element.
Post-process fields on refined mesh.
Operators and data structures for thermal analysis.
Implementation of Gel finite element.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRQ_MOAB(a)
check error code of MoAB function
@ MOFEM_STD_EXCEPTION_THROW
@ MOFEM_ATOM_TEST_INVALID
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
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 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 DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
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 add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
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 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
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_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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
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
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Common data for gel model.
string spatialPositionName
definition of volume element
Calculating right hand side.
Calculate internal forces for solvent flux.
Assemble internal force vector.
Implementation of Gel constitutive model.
map< int, BlockMaterialData > blockMaterialData
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > constitutiveEquationPtr
User (hackable) Gel model.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
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.
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for Time Stepping (TS) solver.
operator for calculate heat flux and assemble to right hand side
Force scale operator for reading two columns.