#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) {
flux =
bcFluxMap.at(block_id)->hookFun(x, y, z);
} else {
}
} else {
flux = 0;
}
}
struct OpRhsBcOnValues
boost::shared_ptr<MethodForForceScaling>
valueScale;
boost::shared_ptr<MethodForForceScaling> &value_scale)
fluxes_name, UserDataOperator::OPROW),
if (data.getFieldData().size() == 0)
EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
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;
x = getCoordsAtGaussPts()(gg, 0);
y = getCoordsAtGaussPts()(gg, 1);
z = getCoordsAtGaussPts()(gg, 2);
double value;
const double w = getGaussPts()(2, gg) * 0.5;
const double beta =
w * (value - z);
noalias(
nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
}
Vec f = getFEMethod()->ts_F;
}
}
};
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>();
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
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),
}
EntityType col_type,
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;
}
&*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),
}
EntityType col_type,
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),
EntityType col_type,
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);
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
}
&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),
EntityType col_type,
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>();
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;
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
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>();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
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,
const std::string field_name)
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 DMSetDefaultSection(
dM, section);
CHKERR DMSetDefaultGlobalSection(
dM, 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>(
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
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->get_postProcess_to_do_Monitor().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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
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 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.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#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.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
DataForcesAndSourcesCore::EntData EntData
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, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::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
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::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, DataForcesAndSourcesCore::EntData &data)
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
UnsaturatedFlowElement & cTx
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
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, DataForcesAndSourcesCore::EntData &data)
Integrate boundary condition.
VectorDouble nF
Vector of residuals.
FTensor::Index< 'j', 3 > j
UnsaturatedFlowElement & cTx
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble matrix.
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations.
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::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 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.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
auto getFTensor0IntegrationWeight()
Get integration weights.
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(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
default operator for TET element
MoFEMErrorCode getDivergenceOfHDivBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, VectorDouble &div)
Get divergence of base functions at integration point.
double getVolume() const
element volume (linear geometry)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Force scale operator for reading two columns.