#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
struct GenericMaterial {
const std::string &prefix) const = 0;
"Not implemented how to calculate hydraulic conductivity");
}
"Not implemented how to calculate derivative of hydraulic "
"conductivity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
};
struct UnsaturatedFlowElement : public MixTransportElement {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
}
int &block_id) const {
for (MaterialsDoubleMap::const_iterator mit =
dMatMap.begin();
if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
block_id = mit->first;
}
}
"Element not found, no material data");
}
struct BcData {
boost::function<
double(
const double x,
const double y,
const double z)>
};
typedef map<int, boost::shared_ptr<BcData>>
BcMap;
const double x, const double y, const double z,
double &value) {
int block_id = -1;
} else {
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
value =
bcValueMap.at(block_id)->hookFun(x, y, z);
} else {
}
} else {
value = 0;
}
}
const double y,
const double z,
double &
flux) {
int block_id = -1;
} else {
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
} else {
}
} else {
}
}
struct OpRhsBcOnValues
boost::shared_ptr<MethodForForceScaling>
valueScale;
boost::shared_ptr<MethodForForceScaling> &value_scale)
fluxes_name, UserDataOperator::OPROW),
if (data.getFieldData().size() == 0)
nF.resize(data.getIndices().size());
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double x, y, z;
double value;
const double beta =
w * (value - z);
}
}
}
};
struct OpResidualFlux
flux_name, UserDataOperator::OPROW),
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv = data.getFTensor1N<3>();
auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
divVec.resize(nb_dofs,
false);
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
v = t_base_diff_hdiv(
i,
i);
++t_base_diff_hdiv;
}
const double alpha = t_w * vol;
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
const double K = block_data->K;
const double scaleZ = block_data->scaleZ;
const double z = t_coords(2);
for (int rr = 0; rr != nb_dofs; rr++) {
t_nf +=
alpha * (1 /
K) * (t_n_hdiv(
i) * t_flux(
i));
++t_n_hdiv;
++t_nf;
}
++t_h;
++t_flux;
++t_coords;
++t_w;
}
&*data.getIndices().begin(), &*
nF.begin(),
ADD_VALUES);
}
};
struct OpResidualMass
val_name, UserDataOperator::OPROW),
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
const double scale = block_data->sCale;
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
const double C = block_data->C;
noalias(
nF) += (
alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
++t_h;
++t_h_t;
++t_div_flux;
++t_coords;
++t_w;
}
ADD_VALUES);
}
};
struct OpTauDotSigma_HdivHdiv
const std::string flux_name)
flux_name, flux_name, UserDataOperator::OPROWCOL),
}
const int nb_row = row_data.getIndices().size();
const int nb_col = col_data.getIndices().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
const double K = block_data->K;
const double alpha = t_w * vol;
for (int kk = 0; kk != nb_row; kk++) {
auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
for (int ll = 0; ll != nb_col; ll++) {
t_a +=
beta * (t_n_hdiv_row(
j) * t_n_hdiv_col(
j));
++t_n_hdiv_col;
++t_a;
}
++t_n_hdiv_row;
}
++t_coords;
++t_h;
++t_w;
}
&*col_data.getIndices().begin(), &*
nN.data().begin(),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
&*row_data.getIndices().begin(),
&*
transNN.data().begin(), ADD_VALUES);
}
}
};
struct OpVU_L2L2
value_name, value_name, UserDataOperator::OPROWCOL),
}
int nb_row = row_data.getIndices().size();
int nb_col = col_data.getIndices().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
auto t_n_row = row_data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->h_t = t_h_t;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calDiffC();
const double C = block_data->C;
const double diffC = block_data->diffC;
for (int kk = 0; kk != nb_row; kk++) {
auto t_n_col = col_data.getFTensor0N(gg, 0);
for (int ll = 0; ll != nb_col; ll++) {
t_a += (
alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
++t_n_col;
++t_a;
}
++t_n_row;
}
++t_w;
++t_coords;
++t_h;
++t_h_t;
}
&col_data.getIndices()[0], &*
nN.data().begin(),
ADD_VALUES);
}
};
struct OpVDivSigma_L2Hdiv
const std::string &val_name_row,
const std::string &flux_name_col)
val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
int nb_row = row_data.getFieldData().size();
int nb_col = col_data.getFieldData().size();
if (nb_row == 0)
if (nb_col == 0)
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
nN.resize(nb_row, nb_col,
false);
auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
v = t_base_diff_hdiv(
i,
i);
++t_base_diff_hdiv;
}
}
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &
nN(0, 0), ADD_VALUES);
}
};
struct OpDivTauU_HdivL2
const std::string &flux_name_col,
const std::string &val_name_row)
flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
int nb_row = row_data.getFieldData().size();
int nb_col = col_data.getFieldData().size();
if (nb_row == 0)
if (nb_col == 0)
nN.resize(nb_row, nb_col,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calDiffK();
const double K = block_data->K;
const double diffK = block_data->diffK;
double alpha = t_w * vol;
v = t_base_diff_hdiv(
i,
i);
++t_base_diff_hdiv;
}
for (int rr = 0; rr != nb_row; rr++) {
double beta =
alpha * (-diffK / KK) * (t_n_hdiv_row(
i) * t_flux(
i));
auto t_n_col = col_data.getFTensor0N(gg, 0);
for (int cc = 0; cc != nb_col; cc++) {
++t_n_col;
++t_a;
}
++t_n_hdiv_row;
}
++t_w;
++t_coords;
++t_h;
++t_flux;
}
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &
nN(0, 0), ADD_VALUES);
}
};
struct OpEvaluateInitiallHead
const std::string &val_name)
val_name, UserDataOperator::OPROW),
if (data.getFieldData().size() == 0)
int nb_dofs = data.getFieldData().size();
int nb_gauss_pts = data.getN().size1();
if (nb_dofs != static_cast<int>(data.getN().size2())) {
"wrong number of dofs");
}
nN.resize(nb_dofs, nb_dofs,
false);
nF.resize(nb_dofs,
false);
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
for (int gg = 0; gg < nb_gauss_pts; gg++) {
nN +=
alpha * outer_prod(data.getN(gg), data.getN(gg));
nF +=
alpha * block_data->initialPcEval() * data.getN(gg);
}
&*
nF.begin(), INSERT_VALUES);
}
};
struct OpIntegrateFluxes
const std::string fluxes_name)
fluxes_name, UserDataOperator::OPROW),
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
auto t_n_hdiv = data.getFTensor1N<3>();
double flux_on_entity = 0;
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
auto t_data = data.getFTensor0FieldData();
for (int rr = 0; rr != nb_dofs; rr++) {
flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(
i) * t_normal(
i));
++t_n_hdiv;
++t_data;
}
++t_w;
}
}
};
struct OpPostProcMaterial
std::vector<EntityHandle> &map_gauss_pts,
field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
int block_id;
boost::shared_ptr<GenericMaterial> &block_data =
cTx.
dMatMap.at(block_id);
Tag th_id;
int def_block_id = -1;
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_block_id);
double zero = 0;
Tag th_theta;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_se;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_k;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_c;
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
block_data->h = t_h;
block_data->x = t_coords(0);
block_data->y = t_coords(1);
block_data->z = t_coords(2);
CHKERR block_data->calTheta();
double theta = block_data->tHeta;
double Se = block_data->Se;
double K = block_data->K;
double C = block_data->C;
++t_h;
++t_coords;
}
}
};
boost::shared_ptr<PostProcVolumeOnRefinedMesh>
postProc;
boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc,
boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
const int frequency)
}
}
int step;
PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
}
double *ghost_flux;
CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Flux at time %6.4g %6.4g\n",
ts_t,
ghost_flux[0]);
}
};
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
MBTET, fluxes);
MBTET, values);
MBTET, values + "_t");
}
}
const std::string &values_name) {
values_name + "_t");
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
dMatMap[it->getMeshsetId()]->tEts, MBTET,
"MIX");
}
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
bcValueMap[it->getMeshsetId()]->eNts, MBTRI,
"MIX_BCVALUE");
}
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
bcFluxMap[it->getMeshsetId()]->eNts, MBTRI,
"MIX_BCFLUX");
}
}
"MIX", "FLUXES", zero_flux_ents);
PetscSection section;
CHKERR PetscSectionDestroy(§ion);
}
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
boost::shared_ptr<ForcesAndSourcesCore>
feVolLhs;
boost::shared_ptr<MethodForForceScaling>
boost::shared_ptr<MethodForForceScaling>
boost::shared_ptr<VectorDouble>
int operator()(
int,
int,
int p_data)
const {
return 2 * p_data + p_data; }
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data;
}
};
struct preProcessVol {
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr
)
CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
fePtr->problemPtr,
"VALUES",
string(
"VALUES") +
"_t",
ROW,
fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
case TSMethod::CTX_TSSETIFUNCTION:
CHKERR VecSetOption(
fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
}
const NumeredDofEntity *dof_ptr;
for (std::vector<int>::iterator it =
cTx.
bcVecIds.begin();
if (auto dof_ptr =
fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
.lock())
dof_ptr->getFieldData() = *vit;
}
} else {
}
break;
default:
break;
}
}
};
struct postProcessVol {
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr
)
case TSMethod::CTX_TSSETIJACOBIAN: {
CHKERR MatAssemblyBegin(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
PETSC_NULL);
CHKERR MatAssemblyBegin(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
fePtr->ts_B, MAT_FINAL_ASSEMBLY);
} break;
case TSMethod::CTX_TSSETIFUNCTION: {
INSERT_VALUES);
}
} break;
default:
break;
}
}
};
ForcesAndSourcesCore::RuleHookFun face_rule =
FaceRule()) {
feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(
mField));
new OpHOSetContravariantPiolaTransformOnFace3D(
HDIV));
new OpHOSetContravariantPiolaTransformOnFace3D(
HDIV));
new OpEvaluateBcOnFluxes(*this, "FLUXES"));
new OpEvaluateInitiallHead(*this, "VALUES"));
feVolRhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(
new OpValuesAtGaussPts(*this, "VALUES"));
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
feVolRhs->getOpPtrVector().push_back(
new OpResidualFlux(*
this,
"FLUXES"));
feVolRhs->getOpPtrVector().push_back(
new OpResidualMass(*
this,
"VALUES"));
feVolRhs->getOpPtrVector().back().opType =
ForcesAndSourcesCore::UserDataOperator::OPROW;
feVolLhs->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(
new OpValuesAtGaussPts(*this, "VALUES"));
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpVU_L2L2(*
this,
"VALUES"));
new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
boost::shared_ptr<FEMethod> null;
boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
CHKERR post_process->generateReferenceElementMesh();
CHKERR post_process->addFieldValuesPostProc(
"VALUES");
CHKERR post_process->addFieldValuesPostProc(
"VALUES_t");
CHKERR post_process->addFieldValuesPostProc(
"FLUXES");
post_process->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
post_process->getOpPtrVector().push_back(
new OpPostProcMaterial(*this, post_process->postProcMesh,
post_process->mapGaussPts, "VALUES"));
boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
flux_integrate->getOpPtrVector().push_back(
new OpIntegrateFluxes(*this, "FLUXES"));
int frequency = 1;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Monitor post proc",
"none");
"-how_often_output",
"frequency how often results are dumped on hard disk", "", frequency,
&frequency, NULL);
ierr = PetscOptionsEnd();
TsCtx *ts_ctx;
ts_ctx->getPostProcessMonitor().push_back(
tsMonitor);
}
int ghosts[] = {0};
CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
}
}
CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D0, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
D0, INSERT_VALUES, SCATTER_REVERSE);
double norm2D0;
PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
}
CHKERR VecGhostUpdateBegin(
D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(
D1, INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(
D1, INSERT_VALUES, SCATTER_REVERSE);
double norm2D1;
PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
}
if (set_initial_pc) {
}
TS ts;
CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
CHKERR TSSetType(ts, TSBEULER);
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
TsCtx *ts_ctx;
SNES snes;
#if PETSC_VERSION_GE(3, 7, 0)
{
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
void *))SNESMonitorFields,
}
#else
{
void *))SNESMonitorFields,
0, 0);
}
#endif
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetTimeStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
PetscPrintf(PETSC_COMM_WORLD,
"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
"%D, linits %D\n",
steps, rejects, snesfails, ftime, nonlinits, linits);
}
};
}
#endif
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
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 DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
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
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 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 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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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.
#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.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
constexpr auto VecSetValues
constexpr auto MatSetValues
double flux
impulse magnitude
constexpr auto field_name
Set integration rule to boundary elements.
double tHeta
Water content.
virtual MoFEMErrorCode calC()
virtual MoFEMErrorCode calSe()
virtual void printMatParameters(const int id, const std::string &prefix) const =0
virtual MoFEMErrorCode calTheta()
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static double ePsilon0
Regularization parameter.
virtual MoFEMErrorCode calK()
virtual MoFEMErrorCode calDiffC()
double sCale
Scale time dependent eq.
double K
Hydraulic conductivity [L/s].
Range tEts
Elements with this material.
double C
Capacity [S^2/L^2].
double Se
Effective saturation.
virtual ~GenericMaterial()
static double scaleZ
Scale z direction.
virtual double initialPcEval() const =0
Initialize head.
double h_t
rate of hydraulic head
static double ePsilon1
Regularization parameter.
virtual MoFEMErrorCode calDiffK()
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
MoFEM::Interface & mField
VectorDouble valuesAtGaussPts
values at integration points on element
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
VectorDouble divergenceAtGaussPts
divergence at integration points on element
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
boost::function< double(const double x, const double y, const double z)> hookFun
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
UnsaturatedFlowElement & cTx
MonitorPostProc(UnsaturatedFlowElement &ctx, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc, boost::shared_ptr< ForcesAndSourcesCore > flux_Integrate, const int frequency)
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
FTensor::Index< 'i', 3 > i
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx, const std::string &flux_name_col, const std::string &val_name_row)
Constructor.
UnsaturatedFlowElement & cTx
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
OpIntegrateFluxes(UnsaturatedFlowElement &ctx, const std::string fluxes_name)
Constructor.
FTensor::Index< 'i', 3 > i
UnsaturatedFlowElement & cTx
OpPostProcMaterial(UnsaturatedFlowElement &ctx, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
FTensor::Index< 'i', 3 > i
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name, boost::shared_ptr< MethodForForceScaling > &value_scale)
Constructor.
UnsaturatedFlowElement & cTx
boost::shared_ptr< MethodForForceScaling > valueScale
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
VectorDouble nF
Vector of residuals.
FTensor::Index< 'j', 3 > j
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
UnsaturatedFlowElement & cTx
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
UnsaturatedFlowElement & cTx
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
UnsaturatedFlowElement & cTx
int operator()(int, int, int p_data) const
MoFEMErrorCode operator()()
postProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
UnsaturatedFlowElement & cTx
boost::shared_ptr< ForcesAndSourcesCore > fePtr
boost::shared_ptr< ForcesAndSourcesCore > fePtr
MoFEMErrorCode operator()()
preProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
UnsaturatedFlowElement & cTx
MoFEMErrorCode calculateEssentialBc()
Calculate boundary conditions for fluxes.
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
Vec D1
Vector with inital head capillary pressure.
map< int, boost::shared_ptr< BcData > > BcMap
MaterialsDoubleMap dMatMap
materials database
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
Vec ghostFlux
Ghost Vector of integrated fluxes.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
std::vector< int > bcVecIds
EntityHandle lastEvalBcValEnt
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
EntityHandle lastEvalBcFluxEnt
UnsaturatedFlowElement(MoFEM::Interface &m_field)
boost::shared_ptr< FEMethod > tsMonitor
BcMap bcValueMap
Store boundary condition for head capillary pressure.
MoFEMErrorCode solveProblem(bool set_initial_pc=true)
solve problem
int lastEvalBcBlockFluxId
~UnsaturatedFlowElement()
MoFEMErrorCode destroyMatrices()
Delete matrices and vector when no longer needed.
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
MoFEMErrorCode calculateInitialPc()
Calculate inital pressure head distribution.
DM dM
Discrete manager for unsaturated flow problem.
MoFEMErrorCode createMatrices()
Create vectors and matrices.
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
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.
virtual int get_comm_rank() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
structure for User Loop Methods on finite elements
default operator for TRI element
VectorDouble & getNormal()
get triangle normal
auto getFTensor1Normal()
get normal as tensor
friend class UserDataOperator
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Vec & ts_F
residual vector
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double getVolume() const
element volume (linear geometry)
friend class UserDataOperator
Force scale operator for reading two columns.