1#ifndef __ELECPHYSOPERATORS2D_HPP__
2#define __ELECPHYSOPERATORS2D_HPP__
9using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch<0>;
11using BoundaryEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch<0>;
16using EntData = DataForcesAndSourcesCore::EntData;
32const double alpha = 0.01;
33const double gma = 0.002;
36const double mu1 = 0.20;
37const double mu2 = 0.30;
67 jac.resize(2, 2,
false);
78 int nb_dofs = data.getIndices().size();
85 int size2 = data.getN().size2();
86 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2()))
88 "wrong number of dofs");
89 nN.resize(nb_dofs, nb_dofs,
false);
90 nF.resize(nb_dofs,
false);
94 auto t_row_tau = data.getFTensor1N<3>();
96 auto dir = getDirection();
97 double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
105 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
106 const double a = t_w * vol;
107 for (
int rr = 0; rr != nb_dofs; rr++) {
108 auto t_col_tau = data.getFTensor1N<3>(gg, 0);
110 for (
int cc = 0; cc != nb_dofs; cc++) {
111 nN(rr, cc) +=
a * (t_row_tau(
i) * t_normal(
i)) *
112 (t_col_tau(
i) * t_normal(
i));
123 for (
auto &dof : data.getFieldDofs()) {
124 dof->getFieldData() =
nF[dof->getEntDofIdx()];
151 int nb_dofs = data.getFieldData().size();
157 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
159 "wrong number of dofs");
160 nN.resize(nb_dofs, nb_dofs,
false);
161 nF.resize(nb_dofs,
false);
165 auto t_row_mass = data.getFTensor0N();
169 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
170 const double a = t_w * vol;
172 for (
int rr = 0; rr != nb_dofs; rr++) {
173 auto t_col_mass = data.getFTensor0N(gg, 0);
174 nF[rr] +=
a * r * t_row_mass;
175 for (
int cc = 0; cc != nb_dofs; cc++) {
176 nN(rr, cc) +=
a * t_row_mass * t_col_mass;
187 for (
auto &dof : data.getFieldDofs()) {
188 dof->getFieldData() =
nF[dof->getEntDofIdx()];
199struct OpSolveRecovery :
public OpFaceEle {
200 typedef boost::function<
double(
const double,
const double,
const double)>
203 boost::shared_ptr<PreviousData> &data_u,
204 boost::shared_ptr<PreviousData> &data_v,
Method runge_kutta4)
207 boost::shared_ptr<PreviousData>
dataU;
208 boost::shared_ptr<PreviousData>
dataV;
216 int nb_dofs = data.getFieldData().size();
219 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
221 "wrong number of dofs");
222 nN.resize(nb_dofs, nb_dofs,
false);
223 nF.resize(nb_dofs,
false);
233 auto t_row_mass = data.getFTensor0N();
237 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
238 const double a = t_w * vol;
240 for (
int rr = 0; rr != nb_dofs; rr++) {
241 auto t_col_mass = data.getFTensor0N(gg, 0);
242 nF[rr] +=
a * vn * t_row_mass;
243 for (
int cc = 0; cc != nb_dofs; cc++) {
244 nN(rr, cc) +=
a * t_row_mass * t_col_mass;
257 for (
auto &dof : data.getFieldDofs()) {
258 dof->getFieldData() =
nF[dof->getEntDofIdx()];
274 typedef boost::function<
double(
const double,
const double)>
277 boost::shared_ptr<PreviousData> &common_datau,
278 boost::shared_ptr<PreviousData> &common_datav,
FUVal rhs_u)
285 const int nb_dofs = data.getIndices().size();
288 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
290 "wrong number of dofs");
291 vecF.resize(nb_dofs,
false);
292 mat.resize(nb_dofs, nb_dofs,
false);
295 const int nb_integration_pts =
getGaussPts().size2();
298 auto t_row_v_base = data.getFTensor0N();
304 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
305 const double a = vol * t_w;
311 double f = t_u_value * (1.0 - t_u_value);
312 for (
int rr = 0; rr != nb_dofs; ++rr) {
313 double rhs =
rhsU(t_u_value, t_v_value);
314 auto t_col_v_base = data.getFTensor0N(gg, 0);
315 vecF[rr] +=
a * rhs * t_row_v_base;
317 for (
int cc = 0; cc != nb_dofs; ++cc) {
318 mat(rr, cc) +=
a * t_row_v_base * t_col_v_base;
409 boost::shared_ptr<PreviousData> &data,
410 std::map<int, BlockData> &block_map)
417 const int nb_dofs = data.getIndices().size();
436 vecF.resize(nb_dofs,
false);
439 const int nb_integration_pts =
getGaussPts().size2();
440 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
442 auto t_tau_base = data.getFTensor1N<3>();
444 auto t_tau_grad = data.getFTensor2DiffN<3, 2>();
449 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
452 const double K_inv = 1. /
K;
453 const double a = vol * t_w;
454 for (
int rr = 0; rr < nb_dofs; ++rr) {
455 double div_base = t_tau_grad(0, 0) + t_tau_grad(1, 1);
456 vecF[rr] += (K_inv * t_tau_base(
i) * t_flux_value(
i) -
457 div_base * t_mass_value) *
483 typedef boost::function<
double(
const double,
const double)>
486 boost::shared_ptr<PreviousData> &datau,
487 boost::shared_ptr<PreviousData> &datav,
FUval rhs_u,
488 std::map<int, BlockData> &block_map,
Range &stim_region)
497 const int nb_dofs = data.getIndices().size();
517 vecF.resize(nb_dofs,
false);
519 const int nb_integration_pts =
getGaussPts().size2();
525 auto t_row_v_base = data.getFTensor0N();
538 double duration = 5.0;
540 if (
T -
dt < c_time && c_time <=
T + duration) {
548 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
549 const double a = vol * t_w;
554 double rhsu =
rhsU(t_u_value, t_v_value);
555 for (
int rr = 0; rr < nb_dofs; ++rr) {
556 vecF[rr] += (t_row_v_base * (t_mass_dot + t_flux_div - 0*rhsu -
factor * stim)) *
a;
596 std::map<int, BlockData> &block_map)
606 const int nb_row_dofs = row_data.getIndices().size();
607 const int nb_col_dofs = col_data.getIndices().size();
609 if (nb_row_dofs && nb_col_dofs) {
627 mat.resize(nb_row_dofs, nb_col_dofs,
false);
629 const int nb_integration_pts =
getGaussPts().size2();
632 auto t_row_tau_base = row_data.getFTensor1N<3>();
637 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
638 const double a = vol * t_w;
640 const double K_inv = 1. /
K;
641 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
642 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
643 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
644 mat(rr, cc) += (K_inv * t_row_tau_base(
i) * t_col_tau_base(
i)) *
a;
654 if (row_side != col_side || row_type != col_type) {
655 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
676 boost::shared_ptr<PreviousData> &data,
677 std::map<int, BlockData> &block_map)
687 const int nb_row_dofs = row_data.getIndices().size();
688 const int nb_col_dofs = col_data.getIndices().size();
690 if (nb_row_dofs && nb_col_dofs) {
707 mat.resize(nb_row_dofs, nb_col_dofs,
false);
709 const int nb_integration_pts =
getGaussPts().size2();
711 auto t_row_tau_base = row_data.getFTensor1N<3>();
713 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 2>();
715 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
718 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
719 const double a = vol * t_w;
721 const double K_inv = 1. /
K;
722 const double K_diff =
B;
724 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
725 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
726 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
727 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1);
728 mat(rr, cc) += (-(t_row_tau_base(
i) * t_flux_value(
i) * K_inv *
729 K_inv * K_diff * t_col_v_base) -
730 (div_row_base * t_col_v_base)) *
765 const int nb_row_dofs = row_data.getIndices().size();
766 const int nb_col_dofs = col_data.getIndices().size();
768 if (nb_row_dofs && nb_col_dofs) {
769 mat.resize(nb_row_dofs, nb_col_dofs,
false);
771 const int nb_integration_pts =
getGaussPts().size2();
773 auto t_row_v_base = row_data.getFTensor0N();
776 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
777 const double a = vol * t_w;
778 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
779 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 2>(gg, 0);
780 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
781 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1);
782 mat(rr, cc) += (t_row_v_base * div_col_base) *
a;
802 typedef boost::function<
double(
const double,
const double)>
FUval;
804 boost::shared_ptr<PreviousData> &datau,
805 boost::shared_ptr<PreviousData> &datav,
814 boost::shared_ptr<PreviousData> &
dataU;
815 boost::shared_ptr<PreviousData> &
dataV;
824 const int nb_row_dofs = row_data.getIndices().size();
825 const int nb_col_dofs = col_data.getIndices().size();
826 if (nb_row_dofs && nb_col_dofs) {
828 mat.resize(nb_row_dofs, nb_col_dofs,
false);
830 const int nb_integration_pts =
getGaussPts().size2();
832 auto t_row_v_base = row_data.getFTensor0N();
840 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
841 const double a = vol * t_w;
842 double dfu =
DRhs_u(t_u_value, t_v_value);
843 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
844 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
845 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
846 mat(rr, cc) += ((ts_a - 0*dfu) * t_row_v_base * t_col_v_base) *
a;
859 if (row_side != col_side || row_type != col_type) {
860 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
969 boost::shared_ptr<PostProc> &post_proc)
978 "out_level_" + boost::lexical_cast<std::string>(
ts_step) +
".h5m");
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.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ 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 DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
const int save_every_nth_step
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< 0 > BoundaryEle
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
BoundaryEle::UserDataOperator OpBoundaryEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
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.
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< PostProc > postProc
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsTauTau(std::string flux_field, boost::shared_ptr< PreviousData > &commonData, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauV(std::string flux_field, std::string mass_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field, boost::shared_ptr< PreviousData > &datau, boost::shared_ptr< PreviousData > &datav, FUval Drhs_u)
boost::function< double(const double, const double)> FUval
boost::shared_ptr< PreviousData > & dataV
boost::shared_ptr< PreviousData > & dataU
boost::shared_ptr< PreviousData > commonDatav
boost::function< double(const double, const double)> FUVal
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonDatau
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &common_datau, boost::shared_ptr< PreviousData > &common_datav, FUVal rhs_u)
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
boost::function< double(const double, const double)> FUval
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonDatau
boost::shared_ptr< PreviousData > commonDatav
std::map< int, BlockData > setOfBlock
OpAssembleStiffRhsV(std::string flux_field, boost::shared_ptr< PreviousData > &datau, boost::shared_ptr< PreviousData > &datav, FUval rhs_u, std::map< int, BlockData > &block_map, Range &stim_region)
Range & essential_bd_ents
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface, double &init_val)
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> Method
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
bool sYmm
If true assume that matrix is symmetric structure.
structure for User Loop Methods on finite elements
default operator for TRI element
double getMeasure()
get measure of element
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ 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