1 #ifndef __UFOPERATORS2D_HPP__
2 #define __UFOPERATORS2D_HPP__
5 #include <BasicFiniteElements.hpp>
22 constexpr
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);
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));
201 EntityType col_type,
EntData &row_data,
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;
322 Monitor(MPI_Comm &comm,
const int &rank, SmartPetscObj<DM> &dm,
323 boost::shared_ptr<PostProc> &post_proc,
double &err)
335 "out_level_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
342 SmartPetscObj<DM>
dM;
#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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr auto VecSetValues
constexpr auto MatSetValues
constexpr double natural_bc_values
MoFEM::FaceElementForcesAndSourcesCore VolEle
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
constexpr double essential_bc_values
double h
convective heat coefficient
FTensor::Index< 'm', 3 > m
bool sYmm
If true assume that matrix is symmetric structure.
default operator for TRI element
double getMeasure()
get measure of element
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
friend class UserDataOperator
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
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()
MoFEMErrorCode operator()()
boost::shared_ptr< PostProc > postProc
MoFEMErrorCode preProcess()
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.