#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
struct GenericMaterial {
const std::string &prefix) const = 0;
virtual MoFEMErrorCode
calK() {
"Not implemented how to calculate hydraulic conductivity");
}
"Not implemented how to calculate derivative of hydraulic "
"conductivity");
}
virtual MoFEMErrorCode
calC() {
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
"Not implemented how to calculate capacity");
}
virtual MoFEMErrorCode
calSe() {
"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) {
flux =
bcFluxMap.at(block_id)->hookFun(x, y, z);
} else {
}
} else {
flux = 0;
}
}
struct OpRhsBcOnValues
UnsaturatedFlowElement &
cTx;
boost::shared_ptr<MethodForForceScaling>
valueScale;
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
boost::shared_ptr<MethodForForceScaling> &value_scale)
fluxes_name, UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
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);
noalias(
nF) += beta * prod(data.getVectorN<3>(gg),
getNormal());
}
}
CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
}
};
struct OpResidualFlux
UnsaturatedFlowElement &
cTx;
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
flux_name, UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
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);
noalias(
nF) -= alpha * (t_h - z * scaleZ) *
divVec;
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
UnsaturatedFlowElement &
cTx;
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
val_name, UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
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++) {
const double alpha = t_w * vol *
scale;
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;
}
CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*
nF.begin(),
ADD_VALUES);
}
};
struct OpTauDotSigma_HdivHdiv
UnsaturatedFlowElement &
cTx;
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx,
const std::string flux_name)
flux_name, flux_name, UserDataOperator::OPROWCOL),
}
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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;
const double beta = alpha * (1 /
K);
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;
}
CHKERR MatSetValues(
a, nb_row, &*row_data.getIndices().begin(), nb_col,
&*col_data.getIndices().begin(), &*
nN.data().begin(),
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
CHKERR MatSetValues(
a, nb_col, &*col_data.getIndices().begin(), nb_row,
&*row_data.getIndices().begin(),
&*
transNN.data().begin(), ADD_VALUES);
}
}
};
struct OpVU_L2L2
UnsaturatedFlowElement &
cTx;
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
value_name, value_name, UserDataOperator::OPROWCOL),
}
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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++) {
double alpha = t_w * vol *
scale;
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;
}
CHKERR MatSetValues(
a, nb_row, &row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &*
nN.data().begin(),
ADD_VALUES);
}
};
struct OpVDivSigma_L2Hdiv
UnsaturatedFlowElement &
cTx;
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx,
const std::string &val_name_row,
const std::string &flux_name_col)
val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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;
}
noalias(
nN) += alpha * outer_prod(row_data.getN(gg),
divVec);
}
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &
nN(0, 0), ADD_VALUES);
}
};
struct OpDivTauU_HdivL2
UnsaturatedFlowElement &
cTx;
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx,
const std::string &flux_name_col,
const std::string &val_name_row)
flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
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;
}
noalias(
nN) -= alpha * outer_prod(
divVec, col_data.getN(gg));
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_a += beta * t_n_col;
++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
UnsaturatedFlowElement &
cTx;
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx,
const std::string &val_name)
val_name, UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
if (data.getFieldData().size() == 0)
auto nb_dofs = data.getFieldData().size();
if (nb_dofs != 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 (auto 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);
}
CHKERR VecSetValues(
cTx.
D1, nb_dofs, &*data.getIndices().begin(),
&*
nF.begin(), INSERT_VALUES);
}
};
struct OpIntegrateFluxes
UnsaturatedFlowElement &
cTx;
OpIntegrateFluxes(UnsaturatedFlowElement &ctx,
const std::string fluxes_name)
fluxes_name, UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
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
UnsaturatedFlowElement &
cTx;
OpPostProcMaterial(UnsaturatedFlowElement &ctx,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
EntitiesFieldData::EntData &data) {
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;
}
}
};
UnsaturatedFlowElement &
cTx;
boost::shared_ptr<PostProcEle>
postProc;
boost::shared_ptr<PostProcEle> &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]);
}
};
MoFEMErrorCode
addFields(
const std::string &values,
const std::string &fluxes,
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");
}
}
BitRefLevel ref_level = BitRefLevel().set(0)) {
CHKERR DMMoFEMSetIsPartitioned(
dM, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(
dM, PETSC_TRUE);
CHKERR DMMoFEMAddElement(
dM,
"MIX_BCFLUX");
CHKERR DMMoFEMAddElement(
dM,
"MIX_BCVALUE");
"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 3 * p_data; }
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data;
}
};
struct preProcessVol {
UnsaturatedFlowElement &
cTx;
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
preProcessVol(
UnsaturatedFlowElement &ctx,
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 {
UnsaturatedFlowElement &
cTx;
boost::shared_ptr<ForcesAndSourcesCore>
fePtr;
postProcessVol(
UnsaturatedFlowElement &ctx,
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;
}
}
};
MoFEMErrorCode
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 OpEvaluateBcOnFluxes(*this, "FLUXES"));
new OpEvaluateInitiallHead(*this, "VALUES"));
CHKERR AddHOOps<3, 3, 3>::add(
feVolRhs->getOpPtrVector(), {HDIV, L2});
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;
CHKERR AddHOOps<3, 3, 3>::add(
feVolLhs->getOpPtrVector(), {HDIV, L2});
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;
auto get_post_process_ele = [&]() {
auto post_process = boost::make_shared<PostProcEle>(
mField);
CHKERR AddHOOps<3, 3, 3>::add(post_process->getOpPtrVector(), {HDIV, L2});
auto values_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto flux_ptr = boost::make_shared<MatrixDouble>();
post_process->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("VALUES", values_ptr));
post_process->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("VALUES", grad_ptr));
post_process->getOpPtrVector().push_back(
new OpCalculateHVecVectorField<3>("FLUXES", flux_ptr));
using OpPPMap = OpPostProcMapInMoab<3, 3>;
post_process->getOpPtrVector().push_back(
new OpPPMap(post_process->getPostProcMesh(),
post_process->getMapGaussPts(),
{{"VALUES", values_ptr}},
{{"GRAD", grad_ptr}, {"FLUXES", flux_ptr}},
{}, {}
));
return post_process;
};
auto post_process = get_post_process_ele();
post_process->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
post_process->getOpPtrVector().push_back(
new OpPostProcMaterial(*this, post_process->getPostProcMesh(),
post_process->getMapGaussPts(), "VALUES"));
boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
CHKERR AddHOOps<2, 3, 3>::add(flux_integrate->getOpPtrVector(), {HDIV});
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();
}
Vec D1;
Vec ghostFlux;
MoFEMErrorCode createMatrices() {
CHKERR DMCreateMatrix(dM, &Aij);
CHKERR DMCreateGlobalVector(dM, &D0);
int ghosts[] = {0};
int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
&ghostFlux);
}
MoFEMErrorCode destroyMatrices() {
CHKERR VecDestroy(&ghostFlux);
}
MoFEMErrorCode calculateEssentialBc() {
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
bcIndices.clear();
CHKERR DMoFEMLoopFiniteElements(dM,
"MIX_BCFLUX", feFaceBc);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
double norm2D0;
CHKERR VecNorm(D0, NORM_2, &norm2D0);
PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
}
MoFEMErrorCode calculateInitialPc() {
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMLoopFiniteElements(dM,
"MIX", feVolInitialPc);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
double norm2D1;
CHKERR VecNorm(D1, NORM_2, &norm2D1);
PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
}
MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
if (set_initial_pc) {
CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
}
CHKERR DMoFEMMeshToLocalVector(dM,
D, INSERT_VALUES, SCATTER_FORWARD);
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
SNES snes;
#if PETSC_VERSION_GE(3, 7, 0)
{
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
}
#else
{
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
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
static PetscErrorCode ierr
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.
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.
#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
const double v
phase velocity of light in medium (cm/ns)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
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
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< PostProcEle > postProc
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
MoFEMErrorCode preProcess()
function is run at the beginning of loop
UnsaturatedFlowElement & cTx
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
UnsaturatedFlowElement & cTx
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
FTensor::Index< 'i', 3 > i
UnsaturatedFlowElement & cTx
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
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
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
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
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()()
UnsaturatedFlowElement & cTx
boost::shared_ptr< ForcesAndSourcesCore > fePtr
boost::shared_ptr< ForcesAndSourcesCore > fePtr
MoFEMErrorCode operator()()
UnsaturatedFlowElement & cTx
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
boost::shared_ptr< FEMethod > tsMonitor
BcMap bcValueMap
Store boundary condition for head capillary pressure.
int lastEvalBcBlockFluxId
~UnsaturatedFlowElement()
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
DM dM
Discrete manager for unsaturated flow problem.
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.
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
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.
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
double getVolume() const
element volume (linear geometry)
Force scale operator for reading two columns.