v0.14.0
UnsaturatedFlow.hpp
/** \file UnsaturatedFlow.hpp
* \brief Mix implementation of transport element
* \example UnsaturatedFlow.hpp
*
* \ingroup mofem_mix_transport_elem
*
*/
#ifndef __UNSATURATD_FLOW_HPP__
#define __UNSATURATD_FLOW_HPP__
namespace MixTransport {
using PostProcEle = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
/**
* \brief Generic material model for unsaturated water transport
*
* \note This is abstact class, no physical material model is implemented
* here.
*/
struct GenericMaterial {
virtual ~GenericMaterial() {}
static double ePsilon0; ///< Regularization parameter
static double ePsilon1; ///< Regularization parameter
static double scaleZ; ///< Scale z direction
double sCale; ///< Scale time dependent eq.
double h; ///< hydraulic head
double h_t; ///< rate of hydraulic head
// double h_flux_residual; ///< residual at point
double K; ///< Hydraulic conductivity [L/s]
double diffK; ///< Derivative of hydraulic conductivity [L/s * L^2/F]
double C; ///< Capacity [S^2/L^2]
double diffC; ///< Derivative of capacity [S^2/L^2 * L^2/F ]
double tHeta; ///< Water content
double Se; ///< Effective saturation
Range tEts; ///< Elements with this material
double x, y, z; ///< in meters (L)
/**
* \brief Initialize head
* @return value of head
*/
virtual double initialPcEval() const = 0;
virtual void printMatParameters(const int id,
const std::string &prefix) const = 0;
virtual MoFEMErrorCode calK() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate hydraulic conductivity");
}
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate derivative of hydraulic "
"conductivity");
}
virtual MoFEMErrorCode calC() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
virtual MoFEMErrorCode calSe() {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
"Not implemented how to calculate capacity");
}
};
/**
* \brief Implementation of operators, problem and finite elements for
* unsaturated flow
*/
struct UnsaturatedFlowElement : public MixTransportElement {
DM dM; ///< Discrete manager for unsaturated flow problem
: MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
if (dM != PETSC_NULL) {
CHKERR DMDestroy(&dM);
CHKERRABORT(PETSC_COMM_WORLD, ierr);
}
}
typedef std::map<int, boost::shared_ptr<GenericMaterial>> MaterialsDoubleMap;
MaterialsDoubleMap dMatMap; ///< materials database
/**
* \brief For given element handle get material block Id
* @param ent finite element entity handle
* @param block_id reference to returned block id
* @return error code
*/
int &block_id) const {
for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
mit != dMatMap.end(); mit++) {
if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
block_id = mit->first;
}
}
"Element not found, no material data");
}
/**
* \brief Class storing information about boundary condition
*/
struct BcData {
double fixValue;
boost::function<double(const double x, const double y, const double z)>
BcData() : hookFun(NULL) {}
};
typedef map<int, boost::shared_ptr<BcData>> BcMap;
BcMap bcValueMap; ///< Store boundary condition for head capillary pressure
/**
* \brief Get value on boundary
* @param ent entity handle
* @param gg number of integration point
* @param x x-coordinate
* @param y y-coordinate
* @param z z-coordinate
* @param value returned value
* @return error code
*/
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
const double x, const double y, const double z,
double &value) {
int block_id = -1;
if (lastEvalBcValEnt == ent) {
} else {
for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
if (bcValueMap.at(block_id)->hookFun) {
value = bcValueMap.at(block_id)->hookFun(x, y, z);
} else {
value = bcValueMap.at(block_id)->fixValue;
}
} else {
value = 0;
}
}
/**
* \brief essential (Neumann) boundary condition (set fluxes)
* @param ent handle to finite element entity
* @param x coord
* @param y coord
* @param z coord
* @param flux reference to flux which is set by function
* @return [description]
*/
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
const double y, const double z, double &flux) {
int block_id = -1;
if (lastEvalBcFluxEnt == ent) {
} else {
for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
it++) {
if (it->second->eNts.find(ent) != it->second->eNts.end()) {
block_id = it->first;
}
}
}
if (block_id >= 0) {
if (bcFluxMap.at(block_id)->hookFun) {
flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
} else {
flux = bcFluxMap.at(block_id)->fixValue;
}
} else {
flux = 0;
}
}
/**
* \brief Evaluate boundary condition at the boundary
*/
struct OpRhsBcOnValues
boost::shared_ptr<MethodForForceScaling> valueScale;
/**
* \brief Constructor
*/
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
boost::shared_ptr<MethodForForceScaling> &value_scale)
fluxes_name, UserDataOperator::OPROW),
cTx(ctx), valueScale(value_scale) {}
VectorDouble nF; ///< Vector of residuals
/**
* \brief Integrate boundary condition
* @param side local index of entity
* @param type type of entity
* @param data data on entity
* @return error code
*/
MoFEMErrorCode doWork(int side, EntityType type,
if (data.getFieldData().size() == 0)
// Get EntityHandle of the finite element
// Resize and clear vector
nF.resize(data.getIndices().size());
nF.clear();
// Get number of integration points
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;
// get value of boundary condition
CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
const double w = getGaussPts()(2, gg) * 0.5;
const double beta = w * (value - z);
noalias(nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
}
// Scale vector if history evaluating method is given
if (valueScale) {
}
// Assemble vector
CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
&nF[0], ADD_VALUES);
}
};
/**
* \brief Assemble flux residual
*/
struct OpResidualFlux
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
nF.resize(nb_dofs, false);
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// block_data->printMatParameters(block_id,"Read material");
// Get base function
auto t_n_hdiv = data.getFTensor1N<3>();
// Get pressure
// Get flux
auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get integration weight
// Get volume
double vol = getVolume();
auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
divVec.resize(nb_dofs, false);
// Get material parameters
int nb_gauss_pts = data.getN().size1();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// Get divergence
for (auto &v : divVec) {
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);
CHKERR block_data->calK();
const double K = block_data->K;
const double scaleZ = block_data->scaleZ;
const double z = t_coords(2); /// z-coordinate at Gauss pt
// Calculate pressure gradient
noalias(nF) -= alpha * (t_h - z * scaleZ) * divVec;
// Calculate presure gradient from flux
FTensor::Tensor0<double *> t_nf(&*nF.begin());
for (int rr = 0; rr != nb_dofs; rr++) {
t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
++t_n_hdiv; // move to next base function
++t_nf; // move to next element in vector
}
++t_h; // move to next integration point
++t_flux;
++t_coords;
++t_w;
}
// Assemble residual
CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
&*data.getIndices().begin(), &*nF.begin(),
ADD_VALUES);
}
};
/**
* Assemble mass residual
*/
struct OpResidualMass
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
const int nb_dofs = data.getIndices().size();
if (nb_dofs == 0)
// Resize local element vector
nF.resize(nb_dofs, false);
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
// Get pressure rate
// Flux divergence
// Get integration weight
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Scale eq.
const double scale = block_data->sCale;
// Get volume
const double vol = getVolume();
// Get number of integration points
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);
CHKERR block_data->calC();
const double C = block_data->C;
// Calculate flux conservation
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;
}
// Assemble local vector
CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
ADD_VALUES);
}
};
struct OpTauDotSigma_HdivHdiv
const std::string flux_name)
flux_name, flux_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
/**
* \brief Assemble matrix
* @param row_side local index of row entity on element
* @param col_side local index of col entity on element
* @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param row_data data for row
* @param col_data data for col
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
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);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get base functions
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
// Get integration weight
// Get volume
const double vol = getVolume();
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->calK();
const double K = block_data->K;
// get integration weight and multiply by element volume
const double alpha = t_w * vol;
const double beta = alpha * (1 / K);
FTensor::Tensor0<double *> t_a(&*nN.data().begin());
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;
}
Mat a = getFEMethod()->ts_B;
CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
&*col_data.getIndices().begin(), &*nN.data().begin(),
ADD_VALUES);
// matrix is symmetric, assemble other part
if (row_side != col_side || row_type != col_type) {
transNN.resize(nb_col, nb_row);
noalias(transNN) = trans(nN);
CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
&*row_data.getIndices().begin(),
&*transNN.data().begin(), ADD_VALUES);
}
}
};
struct OpVU_L2L2
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
value_name, value_name, UserDataOperator::OPROWCOL),
cTx(ctx) {
sYmm = true;
}
/**
* \brief Assemble matrix
* @param row_side local index of row entity on element
* @param col_side local index of col entity on element
* @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
* @param row_data data for row
* @param col_data data for col
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
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);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
// // Get pressure
// auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
// Get pressure rate
// Get integration weight
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Scale eq.
const double scale = block_data->sCale;
// Time step factor
double ts_a = getFEMethod()->ts_a;
// get volume
const double vol = getVolume();
int nb_gauss_pts = row_data.getN().size1();
// get base functions on rows
auto t_n_row = row_data.getFTensor0N();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
// get integration weight and multiply by element volume
double alpha = t_w * vol * scale;
// evaluate material model at integration points
// to calculate capacity and tangent of capacity term
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->calC();
CHKERR block_data->calDiffC();
const double C = block_data->C;
const double diffC = block_data->diffC;
// assemble local entity tangent matrix
&*nN.data().begin());
// iterate base functions on rows
for (int kk = 0; kk != nb_row; kk++) {
// get first base function on column at integration point gg
auto t_n_col = col_data.getFTensor0N(gg, 0);
// iterate base functions on columns
for (int ll = 0; ll != nb_col; ll++) {
// assemble elements of local matrix
t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
++t_n_col; // move to next base function on column
++t_a; // move to next element in local tangent matrix
}
++t_n_row; // move to next base function on row
}
++t_w; // move to next integration weight
++t_coords; // move to next coordinate at integration point
++t_h; // move to next capillary head at integration point
// ++t_flux_residual;
++t_h_t; // move to next capillary head rate at integration point
}
Mat a = getFEMethod()->ts_B;
CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &*nN.data().begin(),
ADD_VALUES);
}
};
struct OpVDivSigma_L2Hdiv
/**
* \brief Constructor
*/
const std::string &val_name_row,
const std::string &flux_name_col)
val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
/**
* \brief Do calculations
* @param row_side local index of entity on row
* @param col_side local index of entity on column
* @param row_type type of row entity
* @param col_type type of col entity
* @param row_data row data structure carrying information about base
* functions, DOFs indices, etc.
* @param col_data column data structure carrying information about base
* functions, DOFs indices, etc.
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
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)
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
nN.resize(nb_row, nb_col, false);
divVec.resize(nb_col, false);
nN.clear();
// take derivatives of base functions
auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
// Scale eq.
const double scale = block_data->sCale;
int nb_gauss_pts = row_data.getN().size1();
for (int gg = 0; gg < nb_gauss_pts; gg++) {
double alpha = getGaussPts()(3, gg) * getVolume() * scale;
for (auto &v : divVec) {
v = t_base_diff_hdiv(i, i);
++t_base_diff_hdiv;
}
noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
}
CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
}
};
struct OpDivTauU_HdivL2
/**
* \brief Constructor
*/
const std::string &flux_name_col,
const std::string &val_name_row)
flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
cTx(ctx) {}
/**
* \brief Do calculations
* @param row_side local index of entity on row
* @param col_side local index of entity on column
* @param row_type type of row entity
* @param col_type type of col entity
* @param row_data row data structure carrying information about base
* functions, DOFs indices, etc.
* @param col_data column data structure carrying information about base
* functions, DOFs indices, etc.
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
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);
divVec.resize(nb_row, false);
nN.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Get pressure
// Get flux
auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// Get integration weight
// Get base function
auto t_n_hdiv_row = row_data.getFTensor1N<3>();
// Get derivatives of base functions
auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
// Get volume
double vol = getVolume();
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->calK();
CHKERR block_data->calDiffK();
const double K = block_data->K;
// const double z = t_coords(2);
const double KK = K * K;
const double diffK = block_data->diffK;
double alpha = t_w * vol;
for (auto &v : divVec) {
v = t_base_diff_hdiv(i, i);
++t_base_diff_hdiv;
}
noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
FTensor::Tensor0<double *> t_a(&*nN.data().begin());
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;
}
CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
&row_data.getIndices()[0], nb_col,
&col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
}
};
struct OpEvaluateInitiallHead
const std::string &val_name)
cTx(ctx) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (data.getFieldData().size() == 0)
auto nb_dofs = data.getFieldData().size();
if (nb_dofs != data.getN().size2()) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"wrong number of dofs");
}
nN.resize(nb_dofs, nb_dofs, false);
nF.resize(nb_dofs, false);
nN.clear();
nF.clear();
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// loop over integration points
auto nb_gauss_pts = getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; gg++) {
// get coordinates at integration point
block_data->x = getCoordsAtGaussPts()(gg, 0);
block_data->y = getCoordsAtGaussPts()(gg, 1);
block_data->z = getCoordsAtGaussPts()(gg, 2);
// get weight for integration rule
double alpha = getGaussPts()(2, gg);
nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
nF += alpha * block_data->initialPcEval() * data.getN(gg);
}
// factor matrix
// solve local problem
cholesky_solve(nN, nF, ublas::lower());
// set solution to vector
CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
&*nF.begin(), INSERT_VALUES);
}
};
struct OpIntegrateFluxes
/**
* \brief Constructor
*/
const std::string fluxes_name)
fluxes_name, UserDataOperator::OPROW),
cTx(ctx) {}
/**
* \brief Integrate boundary condition
* @param side local index of entity
* @param type type of entity
* @param data data on entity
* @return error code
*/
MoFEMErrorCode doWork(int side, EntityType type,
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
// Get base function
auto t_n_hdiv = data.getFTensor1N<3>();
// get normal of face
auto t_normal = getFTensor1Normal();
// Integration weight
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;
}
CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
}
};
/**
* Operator used to post-process results for unsaturated infiltration problem.
* Operator should with element for post-processing results, i.e.
* PostProcVolumeOnRefinedMesh
*/
struct OpPostProcMaterial
std::vector<EntityHandle> &mapGaussPts;
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
const std::string field_name)
field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
MoFEMErrorCode doWork(int side, EntityType type,
int nb_dofs = data.getFieldData().size();
if (nb_dofs == 0)
// Get EntityHandle of the finite element
// Get material block id
int block_id;
CHKERR cTx.getMaterial(fe_ent, block_id);
// Get material block
boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
// Set bloc Id
Tag th_id;
int def_block_id = -1;
CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
MB_TAG_CREAT | MB_TAG_SPARSE,
&def_block_id);
// Create mesh tag. Tags are created on post-processing mesh and
// visable in post-processor, e.g. Paraview
double zero = 0;
Tag th_theta;
CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_se;
CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
// Tag th_ks;
// CHKERR postProcMesh.tag_get_handle(
// "Ks",1,MB_TYPE_DOUBLE,th_ks,
// MB_TAG_CREAT|MB_TAG_SPARSE,&zero
// );
// CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
Tag th_k;
CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
Tag th_c;
CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
// Get pressure at integration points
// Coords at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
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);
// Calculate theta (water content) and save it on mesh tags
CHKERR block_data->calTheta();
double theta = block_data->tHeta;
CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
CHKERR block_data->calSe();
// Calculate Se (effective saturation and save it on the mesh tags)
double Se = block_data->Se;
CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
// Calculate K (hydraulic conductivity) and save it on the mesh tags
CHKERR block_data->calK();
double K = block_data->K;
CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
// Calculate water capacity and save it on the mesh tags
CHKERR block_data->calC();
double C = block_data->C;
CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
// Set block iD
CHKERR postProcMesh.tag_set_data(th_id, &mapGaussPts[gg], 1, &block_id);
++t_h;
++t_coords;
}
}
};
/**
* Finite element implementation called by TS monitor. Element calls other
* finite elements to evaluate material properties and save results on the
* mesh.
*
* \note Element overloaded only FEMethod::postProcess methods where other
* elements are called.
*/
struct MonitorPostProc : public MoFEM::FEMethod {
boost::shared_ptr<PostProcEle> postProc;
boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
const int fRequency;
boost::shared_ptr<PostProcEle> &post_proc,
boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
const int frequency)
: cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
fRequency(frequency) {}
}
}
// Get time step
int step;
CHKERR TSGetTimeStepNumber(ts, &step);
if ((step) % fRequency == 0) {
// Post-process results and save in the file
PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
CHKERR postProc->writeFile(
string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
}
// Integrate fluxes on faces where pressure head is applied
CHKERR VecZeroEntries(cTx.ghostFlux);
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
// Run finite element to integrate fluxes
CHKERR VecAssemblyBegin(cTx.ghostFlux);
CHKERR VecAssemblyEnd(cTx.ghostFlux);
// accumulate errors from processors
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
// scatter errors to all processors
CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
double *ghost_flux;
CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
ghost_flux[0]);
CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
}
};
/// \brief add fields
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
const int order) {
// Fields
// CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
// meshset consisting all entities in mesh
EntityHandle root_set = mField.get_moab().get_root_set();
// add entities to field
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, fluxes);
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, values);
CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
MBTET, values + "_t");
// CHKERR mField.add_ents_to_field_by_type(
// dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
// );
}
CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
CHKERR mField.set_field_order(root_set, MBTET, values, order);
CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
// CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
}
/// \brief add finite elements
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
const std::string &values_name) {
// Define element "MIX". Note that this element will work with fluxes_name
// and values_name. This reflect bilinear form for the problem
values_name + "_t");
// CHKERR
// mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
if (it->getName().compare(0, 4, "SOIL") != 0)
continue;
dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
}
// Define element to integrate natural boundary conditions, i.e. set values.
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "HEAD") != 0)
continue;
bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
}
// Define element to apply essential boundary conditions.
fluxes_name);
fluxes_name);
fluxes_name);
values_name);
if (it->getName().compare(0, 4, "FLUX") != 0)
continue;
bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
}
}
/**
* \brief Build problem
* @param ref_level mesh refinement on which mesh problem you like to built.
* @return error code
*/
BitRefLevel ref_level = BitRefLevel().set(0)) {
// Build fields
// Build finite elements
// Build adjacencies of degrees of freedom and elements
// create DM instance
CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
// setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
CHKERR DMSetType(dM, "DMMOFEM");
// mesh is portioned, each process keeps only part of problem
// creates problem in DM
CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
// discretised problem creates square matrix (that makes some optimizations)
// set DM options from command line
CHKERR DMSetFromOptions(dM);
// add finite elements
CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
// constructor data structures
CHKERR DMSetUp(dM);
// remove zero flux dofs
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
"MIX", "FLUXES", zero_flux_ents);
PetscSection section;
CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
CHKERR DMSetSection(dM, section);
// CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
CHKERR PetscSectionDestroy(&section);
}
boost::shared_ptr<ForcesAndSourcesCore>
feFaceBc; ///< Elemnet to calculate essential bc
boost::shared_ptr<ForcesAndSourcesCore>
feFaceRhs; ///< Face element apply natural bc
boost::shared_ptr<ForcesAndSourcesCore>
feVolInitialPc; ///< Calculate inital boundary conditions
boost::shared_ptr<ForcesAndSourcesCore>
feVolRhs; ///< Assemble residual vector
boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
boost::shared_ptr<MethodForForceScaling>
scaleMethodFlux; ///< Method scaling fluxes
boost::shared_ptr<MethodForForceScaling>
scaleMethodValue; ///< Method scaling values
boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
///< postprocess results at time step
boost::shared_ptr<VectorDouble>
headRateAtGaussPts; ///< Vector keeps head rate
// boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
/**
* \brief Set integration rule to volume elements
*
*/
struct VolRule {
int operator()(int, int, int p_data) const { return 3 * p_data; }
};
/**
* \brief Set integration rule to boundary elements
*
*/
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const {
return 2 * p_data;
}
};
std::vector<int> bcVecIds;
/**
* \brief Pre-peprocessing
* Set head pressure rate and get inital essential boundary conditions
*/
struct preProcessVol {
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
// Update pressure rates
CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
switch (fePtr->ts_ctx) {
case TSMethod::CTX_TSSETIFUNCTION:
CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
if (!cTx.bcIndices.empty()) {
double scale;
CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
cTx.bcIndices.end());
cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
}
CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(),
&*cTx.vecValsOnBc.begin());
VectorDouble::iterator vit = cTx.bcVecVals.begin();
const NumeredDofEntity *dof_ptr;
for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
it != cTx.bcVecIds.end(); it++, vit++) {
if (auto dof_ptr =
fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
.lock())
dof_ptr->getFieldData() = *vit;
}
} else {
cTx.bcVecIds.resize(0);
cTx.bcVecVals.resize(0);
cTx.vecValsOnBc.resize(0);
}
break;
default:
// don nothing
break;
}
}
};
/**
* \brief Post proces method for volume element
* Assemble vectors and matrices and apply essential boundary conditions
*/
struct postProcessVol {
boost::shared_ptr<ForcesAndSourcesCore> fePtr;
// std::string mArk;
boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
)
: cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
switch (fePtr->ts_ctx) {
case TSMethod::CTX_TSSETIJACOBIAN: {
CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
// MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
// std::string wait;
// std::cin >> wait;
CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
&*cTx.bcVecIds.begin(), 1, PETSC_NULL,
PETSC_NULL);
CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
// MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
// std::string wait;
// std::cin >> wait;
} break;
case TSMethod::CTX_TSSETIFUNCTION: {
CHKERR VecAssemblyBegin(fePtr->ts_F);
CHKERR VecAssemblyEnd(fePtr->ts_F);
if (!cTx.bcVecIds.empty()) {
&*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
INSERT_VALUES);
}
CHKERR VecAssemblyBegin(fePtr->ts_F);
CHKERR VecAssemblyEnd(fePtr->ts_F);
} break;
default:
// don nothing
break;
}
}
};
/**
* \brief Create finite element instances
* @param vol_rule integration rule for volume element
* @param face_rule integration rule for boundary element
* @return error code
*/
setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
// create finite element instances
feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
// set integration rule to elements
feFaceBc->getRuleHook = face_rule;
feFaceRhs->getRuleHook = face_rule;
feVolInitialPc->getRuleHook = vol_rule;
feVolLhs->getRuleHook = vol_rule;
feVolRhs->getRuleHook = vol_rule;
// set function hook for finite element postprocessing stage
feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
// Set Piola transform
CHKERR AddHOOps<2, 3, 3>::add(feFaceBc->getOpPtrVector(), {HDIV});
CHKERR AddHOOps<2, 3, 3>::add(feFaceRhs->getOpPtrVector(), {HDIV});
// create method for setting history for fluxes on boundary
scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
new TimeForceScale("-flux_history", false));
// create method for setting history for presure heads on boundary
scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
new TimeForceScale("-head_history", false));
// Set operator to calculate essential boundary conditions
feFaceBc->getOpPtrVector().push_back(
new OpEvaluateBcOnFluxes(*this, "FLUXES"));
// Set operator to calculate initial capillary pressure
CHKERR AddHOOps<3, 3, 3>::add(feVolInitialPc->getOpPtrVector(), {HDIV, L2});
feVolInitialPc->getOpPtrVector().push_back(
new OpEvaluateInitiallHead(*this, "VALUES"));
// set residual face from Neumann terms, i.e. applied pressure
feFaceRhs->getOpPtrVector().push_back(
new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
// set residual finite element operators
headRateAtGaussPts = boost::make_shared<VectorDouble>();
// resAtGaussPts = boost::make_shared<VectorDouble>();
CHKERR AddHOOps<3, 3, 3>::add(feVolRhs->getOpPtrVector(), {HDIV, L2});
feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
string("VALUES") + "_t", headRateAtGaussPts, MBTET));
feVolRhs->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
feVolRhs->getOpPtrVector().push_back(
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;
// set tangent matrix finite element operators
CHKERR AddHOOps<3, 3, 3>::add(feVolLhs->getOpPtrVector(), {HDIV, L2});
feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
string("VALUES") + "_t", headRateAtGaussPts, MBTET));
// feVolLhs->getOpPtrVector().push_back(
// new
// OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
// );
feVolLhs->getOpPtrVector().push_back(
new OpValuesAtGaussPts(*this, "VALUES"));
feVolLhs->getOpPtrVector().push_back(
new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
feVolLhs->getOpPtrVector().push_back(
new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
feVolLhs->getOpPtrVector().push_back(
new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
// Adding finite elements to DM, time solver will ask for it to assemble
// tangent matrices and residuals
boost::shared_ptr<FEMethod> null;
CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
// setting up post-processing
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"));
// Integrate fluxes on boundary
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");
CHKERR PetscOptionsInt(
"-how_often_output",
"frequency how often results are dumped on hard disk", "", frequency,
&frequency, NULL);
ierr = PetscOptionsEnd();
tsMonitor = boost::shared_ptr<FEMethod>(
new MonitorPostProc(*this, post_process, flux_integrate, frequency));
TsCtx *ts_ctx;
}
Vec D1; ///< Vector with inital head capillary pressure
Vec ghostFlux; ///< Ghost Vector of integrated fluxes
/**
* \brief Create vectors and matrices
* @return Error code
*/
MoFEMErrorCode createMatrices() {
CHKERR DMCreateMatrix(dM, &Aij);
CHKERR DMCreateGlobalVector(dM, &D0);
CHKERR VecDuplicate(D0, &D1);
CHKERR VecDuplicate(D0, &D);
CHKERR VecDuplicate(D0, &F);
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);
}
/**
* \brief Delete matrices and vector when no longer needed
* @return error code
*/
MoFEMErrorCode destroyMatrices() {
CHKERR MatDestroy(&Aij);
CHKERR VecDestroy(&D);
CHKERR VecDestroy(&D0);
CHKERR VecDestroy(&D1);
CHKERR VecDestroy(&F);
CHKERR VecDestroy(&ghostFlux);
}
/**
* \brief Calculate boundary conditions for fluxes
* @return Error code
*/
MoFEMErrorCode calculateEssentialBc() {
// clear vectors
CHKERR VecZeroEntries(D0);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
// clear essential bc indices, it could have dofs from other mesh refinement
bcIndices.clear();
// set operator to calculate essential boundary conditions
CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
CHKERR VecAssemblyBegin(D0);
CHKERR VecAssemblyEnd(D0);
CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
double norm2D0;
CHKERR VecNorm(D0, NORM_2, &norm2D0);
// CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
}
/**
* \brief Calculate inital pressure head distribution
* @return Error code
*/
MoFEMErrorCode calculateInitialPc() {
// clear vectors
CHKERR VecZeroEntries(D1);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
// Calculate initial pressure head on each element
CHKERR DMoFEMLoopFiniteElements(dM, "MIX", feVolInitialPc);
// Assemble vector
CHKERR VecAssemblyBegin(D1);
CHKERR VecAssemblyEnd(D1);
CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
// Calculate norm
double norm2D1;
CHKERR VecNorm(D1, NORM_2, &norm2D1);
// CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
}
/**
* \brief solve problem
* @return error code
*/
MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
if (set_initial_pc) {
// Set initial head
CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
}
// Initiate vector from data on the mesh
CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
// Create time solver
TS ts;
CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
// Use backward Euler method
CHKERR TSSetType(ts, TSBEULER);
// Set final time
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
// Setup solver from commabd line
CHKERR TSSetFromOptions(ts);
// Set DM to TS
CHKERR TSSetDM(ts, dM);
#if PETSC_VERSION_GE(3, 7, 0)
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
#endif
// Set-up monitor
TsCtx *ts_ctx;
CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
// This add SNES monitor, to show error by fields. It is dirty trick
// to add monitor, so code is hidden from doxygen
CHKERR TSSetSolution(ts, D);
CHKERR TSSetUp(ts);
SNES snes;
CHKERR TSGetSNES(ts, &snes);
#if PETSC_VERSION_GE(3, 7, 0)
{
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
CHKERR SNESMonitorSet(
snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
}
#else
{
CHKERR SNESMonitorSet(snes,
(MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
0, 0);
}
#endif
CHKERR TSSolve(ts, D);
// Get statistic form TS and print it
CHKERR TSGetTime(ts, &ftime);
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);
}
};
} // namespace MixTransport
#endif // __UNSATURATD_FLOW_HPP__
MixTransport::GenericMaterial::ePsilon0
static double ePsilon0
Regularization parameter.
Definition: UnsaturatedFlow.hpp:26
MixTransport::GenericMaterial::diffK
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
Definition: UnsaturatedFlow.hpp:36
MixTransport::GenericMaterial::calK
virtual MoFEMErrorCode calK()
Definition: UnsaturatedFlow.hpp:54
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
MixTransport::UnsaturatedFlowElement::OpRhsBcOnValues::nF
VectorDouble nF
Vector of residuals.
Definition: UnsaturatedFlow.hpp:256
MixTransport::UnsaturatedFlowElement::MonitorPostProc::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: UnsaturatedFlow.hpp:1073
MixTransport::UnsaturatedFlowElement::MonitorPostProc::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:1072
MixTransport::UnsaturatedFlowElement::OpVDivSigma_L2Hdiv::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:669
MixTransport::UnsaturatedFlowElement::BcData::BcData
BcData()
Definition: UnsaturatedFlow.hpp:150
MixTransport::MixTransportElement::MixTransportElement
MixTransportElement(MoFEM::Interface &m_field)
construction of this data structure
Definition: MixTransportElement.hpp:70
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
MoFEM::TSMethod::ts_a
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Definition: LoopMethods.hpp:160
MixTransport::UnsaturatedFlowElement::OpResidualFlux::divVec
VectorDouble divVec
Definition: UnsaturatedFlow.hpp:314
MixTransport::UnsaturatedFlowElement::OpResidualFlux::i
FTensor::Index< 'i', 3 > i
Definition: UnsaturatedFlow.hpp:315
MixTransport::UnsaturatedFlowElement::OpResidualFlux::nF
VectorDouble nF
Definition: UnsaturatedFlow.hpp:314
EntityHandle
MixTransport::UnsaturatedFlowElement::OpVU_L2L2::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:562
MoFEM::TSMethod::ts_F
Vec & ts_F
residual vector
Definition: LoopMethods.hpp:168
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MixTransport::GenericMaterial::diffC
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
Definition: UnsaturatedFlow.hpp:38
MixTransport::UnsaturatedFlowElement::dM
DM dM
Discrete manager for unsaturated flow problem.
Definition: UnsaturatedFlow.hpp:104
MixTransport::UnsaturatedFlowElement::tsMonitor
boost::shared_ptr< FEMethod > tsMonitor
Definition: UnsaturatedFlow.hpp:1299
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MixTransport::UnsaturatedFlowElement::lastEvalBcBlockValId
int lastEvalBcBlockValId
Definition: UnsaturatedFlow.hpp:156
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MixTransport::UnsaturatedFlowElement::OpVDivSigma_L2Hdiv::nN
MatrixDouble nN
Definition: UnsaturatedFlow.hpp:681
MixTransport::UnsaturatedFlowElement::OpRhsBcOnValues::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:244
MoFEM::TsCtx::getPostProcessMonitor
BasicMethodsSequence & getPostProcessMonitor()
Get the postProcess to do Monitor object.
Definition: TsCtx.hpp:144
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::OpTauDotSigma_HdivHdiv
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
Definition: UnsaturatedFlow.hpp:466
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MixTransport::GenericMaterial::h_t
double h_t
rate of hydraulic head
Definition: UnsaturatedFlow.hpp:33
MixTransport::GenericMaterial::x
double x
Definition: UnsaturatedFlow.hpp:44
MixTransport::GenericMaterial::printMatParameters
virtual void printMatParameters(const int id, const std::string &prefix) const =0
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MixTransport::UnsaturatedFlowElement::OpResidualMass::nF
VectorDouble nF
Definition: UnsaturatedFlow.hpp:403
MixTransport::GenericMaterial::calDiffK
virtual MoFEMErrorCode calDiffK()
Definition: UnsaturatedFlow.hpp:61
MixTransport::UnsaturatedFlowElement::postProcessVol::postProcessVol
postProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
Definition: UnsaturatedFlow.hpp:1396
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::transNN
MatrixDouble transNN
Definition: UnsaturatedFlow.hpp:474
MixTransport::UnsaturatedFlowElement::OpIntegrateFluxes::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:911
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::OpDivTauU_HdivL2
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx, const std::string &flux_name_col, const std::string &val_name_row)
Constructor.
Definition: UnsaturatedFlow.hpp:746
MixTransport::UnsaturatedFlowElement::OpResidualFlux::OpResidualFlux
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
Definition: UnsaturatedFlow.hpp:309
MixTransport::UnsaturatedFlowElement::OpRhsBcOnValues::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
Definition: UnsaturatedFlow.hpp:265
MixTransport::UnsaturatedFlowElement::OpPostProcMaterial::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:967
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:464
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
MixTransport::UnsaturatedFlowElement::OpPostProcMaterial::OpPostProcMaterial
OpPostProcMaterial(UnsaturatedFlowElement &ctx, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
Definition: UnsaturatedFlow.hpp:971
MixTransport::UnsaturatedFlowElement::OpEvaluateInitiallHead::OpEvaluateInitiallHead
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
Definition: UnsaturatedFlow.hpp:849
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::j
FTensor::Index< 'j', 3 > j
Definition: UnsaturatedFlow.hpp:476
MixTransport::UnsaturatedFlowElement::OpEvaluateInitiallHead::nF
VectorDouble nF
Definition: UnsaturatedFlow.hpp:856
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MixTransport::MixTransportElement::fluxesAtGaussPts
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
Definition: MixTransportElement.hpp:83
MixTransport::UnsaturatedFlowElement::BcData::fixValue
double fixValue
Definition: UnsaturatedFlow.hpp:147
MixTransport::UnsaturatedFlowElement::postProcessVol::operator()
MoFEMErrorCode operator()()
Definition: UnsaturatedFlow.hpp:1401
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
MixTransport::UnsaturatedFlowElement::postProcessVol::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:1393
MixTransport::UnsaturatedFlowElement::vecValsOnBc
VectorDouble vecValsOnBc
Definition: UnsaturatedFlow.hpp:1324
scale
double scale
Definition: plastic.cpp:119
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::i
FTensor::Index< 'i', 3 > i
Definition: UnsaturatedFlow.hpp:755
MixTransport::UnsaturatedFlowElement::lastEvalBcFluxEnt
EntityHandle lastEvalBcFluxEnt
Definition: UnsaturatedFlow.hpp:198
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
MixTransport::UnsaturatedFlowElement::OpVDivSigma_L2Hdiv::OpVDivSigma_L2Hdiv
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
Definition: UnsaturatedFlow.hpp:674
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MixTransport::GenericMaterial::C
double C
Capacity [S^2/L^2].
Definition: UnsaturatedFlow.hpp:37
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1240
MixTransport::UnsaturatedFlowElement::feVolLhs
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
Definition: UnsaturatedFlow.hpp:1294
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MixTransport::UnsaturatedFlowElement::MonitorPostProc::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: UnsaturatedFlow.hpp:1085
ROW
@ ROW
Definition: definitions.h:136
MixTransport::GenericMaterial::z
double z
in meters (L)
Definition: UnsaturatedFlow.hpp:44
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
MixTransport::UnsaturatedFlowElement::headRateAtGaussPts
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
Definition: UnsaturatedFlow.hpp:1303
MixTransport::GenericMaterial::calC
virtual MoFEMErrorCode calC()
Definition: UnsaturatedFlow.hpp:69
MixTransport::UnsaturatedFlowElement::MonitorPostProc::fRequency
const int fRequency
Definition: UnsaturatedFlow.hpp:1076
MixTransport::UnsaturatedFlowElement::MaterialsDoubleMap
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
Definition: UnsaturatedFlow.hpp:118
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
MixTransport::GenericMaterial::~GenericMaterial
virtual ~GenericMaterial()
Definition: UnsaturatedFlow.hpp:24
MixTransport::GenericMaterial::sCale
double sCale
Scale time dependent eq.
Definition: UnsaturatedFlow.hpp:30
MixTransport::UnsaturatedFlowElement::feFaceBc
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
Definition: UnsaturatedFlow.hpp:1287
MixTransport::GenericMaterial::y
double y
Definition: UnsaturatedFlow.hpp:44
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MixTransport::UnsaturatedFlowElement::OpEvaluateInitiallHead::nN
MatrixDouble nN
Definition: UnsaturatedFlow.hpp:855
MixTransport::UnsaturatedFlowElement::preProcessVol::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:1331
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MixTransport::GenericMaterial::scaleZ
static double scaleZ
Scale z direction.
Definition: UnsaturatedFlow.hpp:28
MixTransport::UnsaturatedFlowElement::MonitorPostProc::MonitorPostProc
MonitorPostProc(UnsaturatedFlowElement &ctx, boost::shared_ptr< PostProcEle > &post_proc, boost::shared_ptr< ForcesAndSourcesCore > flux_Integrate, const int frequency)
Definition: UnsaturatedFlow.hpp:1078
MixTransport::UnsaturatedFlowElement::OpVDivSigma_L2Hdiv::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
Definition: UnsaturatedFlow.hpp:696
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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
Definition: DMMoFEM.cpp:853
MixTransport::UnsaturatedFlowElement::OpPostProcMaterial::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: UnsaturatedFlow.hpp:969
MixTransport::UnsaturatedFlowElement::BcData::hookFun
boost::function< double(const double x, const double y, const double z)> hookFun
Definition: UnsaturatedFlow.hpp:149
MixTransport::UnsaturatedFlowElement::OpVU_L2L2::OpVU_L2L2
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
Definition: UnsaturatedFlow.hpp:564
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MixTransport::GenericMaterial::K
double K
Hydraulic conductivity [L/s].
Definition: UnsaturatedFlow.hpp:35
a
constexpr double a
Definition: approx_sphere.cpp:30
MixTransport::GenericMaterial::calSe
virtual MoFEMErrorCode calSe()
Definition: UnsaturatedFlow.hpp:90
MixTransport::UnsaturatedFlowElement::OpPostProcMaterial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: UnsaturatedFlow.hpp:979
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:741
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getFTensor1Normal
auto getFTensor1Normal()
get normal as tensor
Definition: FaceElementForcesAndSourcesCore.hpp:255
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MixTransport::UnsaturatedFlowElement::OpVU_L2L2::nN
MatrixDouble nN
Definition: UnsaturatedFlow.hpp:571
MixTransport::UnsaturatedFlowElement::OpIntegrateFluxes::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
Definition: UnsaturatedFlow.hpp:931
MixTransport::UnsaturatedFlowElement::VolRule::operator()
int operator()(int, int, int p_data) const
Definition: UnsaturatedFlow.hpp:1311
MoFEM::TSMethod::ts
TS ts
time solver
Definition: LoopMethods.hpp:155
MixTransport::UnsaturatedFlowElement::MonitorPostProc::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: UnsaturatedFlow.hpp:1090
double
convert.type
type
Definition: convert.py:64
MixTransport::UnsaturatedFlowElement::dMatMap
MaterialsDoubleMap dMatMap
materials database
Definition: UnsaturatedFlow.hpp:119
MixTransport::UnsaturatedFlowElement::postProcessVol::fePtr
boost::shared_ptr< ForcesAndSourcesCore > fePtr
Definition: UnsaturatedFlow.hpp:1394
MixTransport::UnsaturatedFlowElement::feVolRhs
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
Definition: UnsaturatedFlow.hpp:1293
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getNormal
VectorDouble & getNormal()
get triangle normal
Definition: FaceElementForcesAndSourcesCore.hpp:243
MixTransport::GenericMaterial::ePsilon1
static double ePsilon1
Regularization parameter.
Definition: UnsaturatedFlow.hpp:27
MixTransport::UnsaturatedFlowElement::OpResidualFlux::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:307
MixTransport::UnsaturatedFlowElement::~UnsaturatedFlowElement
~UnsaturatedFlowElement()
Definition: UnsaturatedFlow.hpp:111
MixTransport::GenericMaterial::calTheta
virtual MoFEMErrorCode calTheta()
Definition: UnsaturatedFlow.hpp:83
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MixTransport::UnsaturatedFlowElement::OpVDivSigma_L2Hdiv::divVec
VectorDouble divVec
Definition: UnsaturatedFlow.hpp:682
MixTransport::GenericMaterial::h
double h
hydraulic head
Definition: UnsaturatedFlow.hpp:32
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
MoFEM::ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Definition: ForcesAndSourcesCore.hpp:1265
MixTransport::UnsaturatedFlowElement::BcMap
map< int, boost::shared_ptr< BcData > > BcMap
Definition: UnsaturatedFlow.hpp:152
MixTransport::UnsaturatedFlowElement::preProcessVol::preProcessVol
preProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
Definition: UnsaturatedFlow.hpp:1335
MixTransport::UnsaturatedFlowElement::lastEvalBcValEnt
EntityHandle lastEvalBcValEnt
Definition: UnsaturatedFlow.hpp:155
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
MixTransport::MixTransportElement::divergenceAtGaussPts
VectorDouble divergenceAtGaussPts
divergence at integration points on element
Definition: MixTransportElement.hpp:82
MixTransport::UnsaturatedFlowElement::addFields
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
Definition: UnsaturatedFlow.hpp:1136
MixTransport::UnsaturatedFlowElement::preProcessVol::operator()
MoFEMErrorCode operator()()
Definition: UnsaturatedFlow.hpp:1340
MixTransport::UnsaturatedFlowElement::OpResidualMass::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: UnsaturatedFlow.hpp:405
MixTransport::UnsaturatedFlowElement::MonitorPostProc::fluxIntegrate
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
Definition: UnsaturatedFlow.hpp:1074
MoFEM::DMMoFEMCreateMoFEM
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.
Definition: DMMoFEM.cpp:114
MixTransport::UnsaturatedFlowElement::OpIntegrateFluxes::OpIntegrateFluxes
OpIntegrateFluxes(UnsaturatedFlowElement &ctx, const std::string fluxes_name)
Constructor.
Definition: UnsaturatedFlow.hpp:916
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MixTransport::UnsaturatedFlowElement::scaleMethodFlux
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
Definition: UnsaturatedFlow.hpp:1296
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MixTransport::UnsaturatedFlowElement::OpVU_L2L2::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
Definition: UnsaturatedFlow.hpp:583
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::divVec
VectorDouble divVec
Definition: UnsaturatedFlow.hpp:754
MixTransport::GenericMaterial::tEts
Range tEts
Elements with this material.
Definition: UnsaturatedFlow.hpp:42
MixTransport::UnsaturatedFlowElement::OpRhsBcOnValues::valueScale
boost::shared_ptr< MethodForForceScaling > valueScale
Definition: UnsaturatedFlow.hpp:245
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
MixTransport::UnsaturatedFlowElement::UnsaturatedFlowElement
UnsaturatedFlowElement(MoFEM::Interface &m_field)
Definition: UnsaturatedFlow.hpp:106
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
MixTransport::MixTransportElement::D0
Vec D0
Definition: MixTransportElement.hpp:475
MixTransport::UnsaturatedFlowElement::D1
Vec D1
Vector with inital head capillary pressure.
Definition: UnsaturatedFlow.hpp:1609
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MixTransport::MixTransportElement::bcIndices
set< int > bcIndices
Definition: MixTransportElement.hpp:85
MixTransport::GenericMaterial::calDiffC
virtual MoFEMErrorCode calDiffC()
Definition: UnsaturatedFlow.hpp:76
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:162
MixTransport::UnsaturatedFlowElement::bcVecIds
std::vector< int > bcVecIds
Definition: UnsaturatedFlow.hpp:1323
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
Definition: ForcesAndSourcesCore.hpp:1004
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MixTransport::GenericMaterial::tHeta
double tHeta
Water content.
Definition: UnsaturatedFlow.hpp:39
MixTransport
Definition: MaterialUnsaturatedFlow.hpp:11
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
Definition: UnsaturatedFlow.hpp:769
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MixTransport::UnsaturatedFlowElement::OpDivTauU_HdivL2::nN
MatrixDouble nN
Definition: UnsaturatedFlow.hpp:753
MixTransport::UnsaturatedFlowElement::OpEvaluateInitiallHead::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: UnsaturatedFlow.hpp:858
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MixTransport::UnsaturatedFlowElement::OpRhsBcOnValues::OpRhsBcOnValues
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name, boost::shared_ptr< MethodForForceScaling > &value_scale)
Constructor.
Definition: UnsaturatedFlow.hpp:250
MixTransport::UnsaturatedFlowElement::ghostFlux
Vec ghostFlux
Ghost Vector of integrated fluxes.
Definition: UnsaturatedFlow.hpp:1610
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
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
Definition: DMMoFEM.cpp:800
MixTransport::UnsaturatedFlowElement::bcVecVals
VectorDouble bcVecVals
Definition: UnsaturatedFlow.hpp:1324
MixTransport::UnsaturatedFlowElement::MonitorPostProc::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: UnsaturatedFlow.hpp:1095
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MixTransport::UnsaturatedFlowElement::preProcessVol::fePtr
boost::shared_ptr< ForcesAndSourcesCore > fePtr
Definition: UnsaturatedFlow.hpp:1332
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
cholesky_decompose
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
Definition: UnsaturatedFlow.hpp:488
cholesky_solve
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MixTransport::UnsaturatedFlowElement::feVolInitialPc
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
Definition: UnsaturatedFlow.hpp:1291
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MixTransport::UnsaturatedFlowElement::OpTauDotSigma_HdivHdiv::nN
MatrixDouble nN
Definition: UnsaturatedFlow.hpp:474
MixTransport::UnsaturatedFlowElement::BcData::eNts
Range eNts
Definition: UnsaturatedFlow.hpp:146
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
MixTransport::UnsaturatedFlowElement::OpResidualMass::OpResidualMass
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
Definition: UnsaturatedFlow.hpp:398
MixTransport::UnsaturatedFlowElement::OpResidualMass::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:396
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MixTransport::UnsaturatedFlowElement::OpPostProcMaterial::postProcMesh
moab::Interface & postProcMesh
Definition: UnsaturatedFlow.hpp:968
MixTransport::UnsaturatedFlowElement::OpIntegrateFluxes::i
FTensor::Index< 'i', 3 > i
Definition: UnsaturatedFlow.hpp:922
MixTransport::UnsaturatedFlowElement::getMaterial
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
Definition: UnsaturatedFlow.hpp:127
MoFEM::CoreInterface::set_field_order
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.
MixTransport::UnsaturatedFlowElement::buildProblem
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
Definition: UnsaturatedFlow.hpp:1241
MonitorPostProc
Definition: nonlinear_dynamics.cpp:30
MoFEM::TSMethod::ts_B
Mat & ts_B
Preconditioner for ts_A.
Definition: LoopMethods.hpp:172
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
MixTransport::MixTransportElement::mField
MoFEM::Interface & mField
Definition: MixTransportElement.hpp:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MixTransport::UnsaturatedFlowElement::OpResidualFlux::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: UnsaturatedFlow.hpp:317
MixTransport::GenericMaterial::Se
double Se
Effective saturation.
Definition: UnsaturatedFlow.hpp:40
MixTransport::UnsaturatedFlowElement::FaceRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: UnsaturatedFlow.hpp:1318
MixTransport::MixTransportElement::valuesAtGaussPts
VectorDouble valuesAtGaussPts
values at integration points on element
Definition: MixTransportElement.hpp:78
MixTransport::UnsaturatedFlowElement::bcFluxMap
BcMap bcFluxMap
Definition: UnsaturatedFlow.hpp:197
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MixTransport::UnsaturatedFlowElement::feFaceRhs
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
Definition: UnsaturatedFlow.hpp:1289
MoFEM::CoreInterface::add_field
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.
MixTransport::UnsaturatedFlowElement::scaleMethodValue
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
Definition: UnsaturatedFlow.hpp:1298
MixTransport::UnsaturatedFlowElement::getBcOnValues
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
Definition: UnsaturatedFlow.hpp:168
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MixTransport::GenericMaterial::initialPcEval
virtual double initialPcEval() const =0
Initialize head.
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
Definition: ForcesAndSourcesCore.hpp:1269
MixTransport::UnsaturatedFlowElement::addFiniteElements
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
Definition: UnsaturatedFlow.hpp:1172
MixTransport::UnsaturatedFlowElement::OpEvaluateInitiallHead::cTx
UnsaturatedFlowElement & cTx
Definition: UnsaturatedFlow.hpp:848
MixTransport::UnsaturatedFlowElement::setFiniteElements
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
Definition: UnsaturatedFlow.hpp:1446
MixTransport::UnsaturatedFlowElement::bcValueMap
BcMap bcValueMap
Store boundary condition for head capillary pressure.
Definition: UnsaturatedFlow.hpp:153
MixTransport::UnsaturatedFlowElement::getBcOnFluxes
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
Definition: UnsaturatedFlow.hpp:210
F
@ F
Definition: free_surface.cpp:394
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
MixTransport::UnsaturatedFlowElement::lastEvalBcBlockFluxId
int lastEvalBcBlockFluxId
Definition: UnsaturatedFlow.hpp:199
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698