1#ifndef __UFOPERATORS2D_HPP__
2#define __UFOPERATORS2D_HPP__
15using EntData = DataForcesAndSourcesCore::EntData;
22constexpr double eps = 1e-16;
45 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Material options",
"none");
46 CHKERR PetscOptionsScalar(
"-K_s",
"K_s",
"",
K_s, &
K_s, PETSC_NULL);
47 CHKERR PetscOptionsScalar(
"-saturated_conductivity",
"K_s",
"",
K_s, &
K_s,
49 CHKERR PetscOptionsScalar(
"-l",
"l",
"",
ll, &
ll, PETSC_NULL);
56 CHKERR PetscOptionsScalar(
"-n",
"n",
"",
nn, &
nn, PETSC_NULL);
58 CHKERR PetscOptionsScalar(
"-h_s",
"h_s",
"",
h_s, &
h_s, PETSC_NULL);
59 ierr = PetscOptionsEnd();
73 if (head_norm <
h_s) {
85 return pow(S_e,
ll) * pow(((1 - F_e) / (1 - F_1)), 2);
88 template <
typename T>
T get_Fe(T eff_satu) {
89 const T m = 1 - 1.0 /
nn;
91 return pow(1 - pow(S_eStar, 1.0 /
m),
m);
103 const T m = 1 - 1.0 /
nn;
105 if (head_norm <
h_s) {
108 pow(1 + pow(-
alpha * head,
nn),
m + 1);
116 std::complex<double> cx_head = head +
eps * 1
i;
118 return capacity.imag() /
eps;
122 std::complex<double> cx_head = head +
eps * 1
i;
124 return conductivity.imag() /
eps;
144 const int nb_dofs = data.getIndices().size();
146 vecF.resize(nb_dofs,
false);
148 const int nb_integration_pts =
getGaussPts().size2();
152 auto t_grad = getFTensor1FromMat<2>(
commonData->grads);
154 auto t_base = data.getFTensor0N();
155 auto t_diff_base = data.getFTensor1DiffN<2>();
162 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
163 const double a = vol * t_w;
168 for (
int rr = 0; rr != nb_dofs; ++rr) {
170 a * (t_base * C_h * t_dot_val + K_h * t_diff_base(
i) * t_grad(
i));
205 const int nb_row_dofs = row_data.getIndices().size();
206 const int nb_col_dofs = col_data.getIndices().size();
209 if (nb_row_dofs && nb_col_dofs) {
210 mat.resize(nb_row_dofs, nb_col_dofs,
false);
212 const int nb_integration_pts =
getGaussPts().size2();
213 auto t_row_base = row_data.getFTensor0N();
217 auto t_grad = getFTensor1FromMat<2>(
commonData->grads);
219 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
226 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
227 const double a = vol * t_w;
233 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
234 auto t_col_base = col_data.getFTensor0N(gg, 0);
235 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
236 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
239 a * (t_row_base * (DC_h * t_dot_val + C_h * ts_a) * t_col_base +
240 DK_h * t_grad(
i) * t_row_diff_base(
i) * t_col_base +
241 K_h * t_row_diff_base(
i) * t_col_diff_base(
i));
279 const int nb_dofs = data.getIndices().size();
287 vecF.resize(nb_dofs,
false);
289 const int nb_integration_pts =
getGaussPts().size2();
290 auto t_row_base = data.getFTensor0N();
295 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
296 const double a = vol * t_w;
299 for (
int rr = 0; rr != nb_dofs; ++rr) {
300 vecF[rr] += t_row_base *
h *
a;
323 boost::shared_ptr<PostProc> &post_proc,
double &err)
335 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
#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.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr double natural_bc_values
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
constexpr double essential_bc_values
bool sYmm
If true assume that matrix is symmetric structure.
default operator for EDGE element
structure for User Loop Methods on finite elements
default operator for TRI element
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
double get_diffCapacity(double head)
T get_conductivity(T head, double head_norm)
T get_waterContent(T head)
T get_capacity(T head, double head_norm)
T get_effSaturation(T head)
T get_eff_satuStar(T eff_satu)
T get_relativeConductivity(T head)
MoFEMErrorCode setOptions()
double get_diffConductivity(double head)
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc, double &err)
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< PostProc > postProc
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr< PreviousData > &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhs(std::string field, boost::shared_ptr< PreviousData > &data)
MatrixDouble grads
Gradients of field "u" at integration points.
VectorDouble values
Values of field "u" at integration points.
VectorDouble dot_values
Rate of values of field "u" at integration points.