1 #ifndef __ELECPHYSOPERATORS2D_HPP__
2 #define __ELECPHYSOPERATORS2D_HPP__
5 #include <BasicFiniteElements.hpp>
10 FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY |
11 FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
12 FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
15 EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY |
16 EdgeElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
37 const double alpha = 0.01;
38 const double gma = 0.002;
39 const double b = 0.15;
40 const double c = 8.00;
41 const double mu1 = 0.20;
42 const double mu2 = 0.30;
72 jac.resize(2, 2,
false);
83 int nb_dofs = data.getIndices().size();
90 int size2 = data.getN().size2();
91 if (3 * nb_dofs !=
static_cast<int>(data.getN().size2()))
93 "wrong number of dofs");
94 nN.resize(nb_dofs, nb_dofs,
false);
95 nF.resize(nb_dofs,
false);
99 auto t_row_tau = data.getFTensor1N<3>();
101 auto dir = getDirection();
102 double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
110 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
111 const double a = t_w * vol;
112 for (
int rr = 0; rr != nb_dofs; rr++) {
113 auto t_col_tau = data.getFTensor1N<3>(gg, 0);
115 for (
int cc = 0; cc != nb_dofs; cc++) {
116 nN(rr, cc) += a * (t_row_tau(
i) * t_normal(
i)) *
117 (t_col_tau(
i) * t_normal(
i));
128 for (
auto &dof : data.getFieldDofs()) {
129 dof->getFieldData() =
nF[dof->getEntDofIdx()];
147 struct OpInitialMass :
public OpFaceEle {
148 OpInitialMass(
const std::string &mass_field, Range &inner_surface,
double &init_val)
156 int nb_dofs = data.getFieldData().size();
162 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
164 "wrong number of dofs");
165 nN.resize(nb_dofs, nb_dofs,
false);
166 nF.resize(nb_dofs,
false);
170 auto t_row_mass = data.getFTensor0N();
174 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
175 const double a = t_w * vol;
177 for (
int rr = 0; rr != nb_dofs; rr++) {
178 auto t_col_mass = data.getFTensor0N(gg, 0);
179 nF[rr] += a *
r * t_row_mass;
180 for (
int cc = 0; cc != nb_dofs; cc++) {
181 nN(rr, cc) += a * t_row_mass * t_col_mass;
192 for (
auto &dof : data.getFieldDofs()) {
193 dof->getFieldData() =
nF[dof->getEntDofIdx()];
204 struct OpSolveRecovery :
public OpFaceEle {
205 typedef boost::function<double(
const double,
const double,
const double)>
208 boost::shared_ptr<PreviousData> &data_u,
209 boost::shared_ptr<PreviousData> &data_v,
Method runge_kutta4)
212 boost::shared_ptr<PreviousData>
dataU;
213 boost::shared_ptr<PreviousData>
dataV;
221 int nb_dofs = data.getFieldData().size();
224 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
226 "wrong number of dofs");
227 nN.resize(nb_dofs, nb_dofs,
false);
228 nF.resize(nb_dofs,
false);
238 auto t_row_mass = data.getFTensor0N();
242 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
243 const double a = t_w * vol;
245 for (
int rr = 0; rr != nb_dofs; rr++) {
246 auto t_col_mass = data.getFTensor0N(gg, 0);
247 nF[rr] += a * vn * t_row_mass;
248 for (
int cc = 0; cc != nb_dofs; cc++) {
249 nN(rr, cc) += a * t_row_mass * t_col_mass;
262 for (
auto &dof : data.getFieldDofs()) {
263 dof->getFieldData() =
nF[dof->getEntDofIdx()];
279 typedef boost::function<double(
const double,
const double)>
282 boost::shared_ptr<PreviousData> &common_datau,
283 boost::shared_ptr<PreviousData> &common_datav,
FUVal rhs_u)
290 const int nb_dofs = data.getIndices().size();
293 if (nb_dofs !=
static_cast<int>(data.getN().size2()))
295 "wrong number of dofs");
296 vecF.resize(nb_dofs,
false);
297 mat.resize(nb_dofs, nb_dofs,
false);
300 const int nb_integration_pts =
getGaussPts().size2();
303 auto t_row_v_base = data.getFTensor0N();
309 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
310 const double a = vol * t_w;
316 double f = t_u_value * (1.0 - t_u_value);
317 for (
int rr = 0; rr != nb_dofs; ++rr) {
318 double rhs =
rhsU(t_u_value, t_v_value);
319 auto t_col_v_base = data.getFTensor0N(gg, 0);
320 vecF[rr] += a * rhs * t_row_v_base;
322 for (
int cc = 0; cc != nb_dofs; ++cc) {
323 mat(rr, cc) += a * t_row_v_base * t_col_v_base;
414 boost::shared_ptr<PreviousData> &data,
422 const int nb_dofs = data.getIndices().size();
441 vecF.resize(nb_dofs,
false);
444 const int nb_integration_pts =
getGaussPts().size2();
445 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
447 auto t_tau_base = data.getFTensor1N<3>();
449 auto t_tau_grad = data.getFTensor2DiffN<3, 2>();
454 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
457 const double K_inv = 1. /
K;
458 const double a = vol * t_w;
459 for (
int rr = 0; rr < nb_dofs; ++rr) {
460 double div_base = t_tau_grad(0, 0) + t_tau_grad(1, 1);
461 vecF[rr] += (K_inv * t_tau_base(
i) * t_flux_value(
i) -
462 div_base * t_mass_value) *
488 typedef boost::function<double(
const double,
const double)>
491 boost::shared_ptr<PreviousData> &datau,
492 boost::shared_ptr<PreviousData> &datav,
FUval rhs_u,
493 std::map<int, BlockData> &
block_map, Range &stim_region)
502 const int nb_dofs = data.getIndices().size();
522 vecF.resize(nb_dofs,
false);
524 const int nb_integration_pts =
getGaussPts().size2();
530 auto t_row_v_base = data.getFTensor0N();
543 double duration = 5.0;
545 if (
T -
dt < c_time && c_time <=
T + duration) {
553 for (
int gg = 0; gg < nb_integration_pts; ++gg) {
554 const double a = vol * t_w;
559 double rhsu =
rhsU(t_u_value, t_v_value);
560 for (
int rr = 0; rr < nb_dofs; ++rr) {
561 vecF[rr] += (t_row_v_base * (t_mass_dot + t_flux_div - 0*rhsu -
factor * stim)) * a;
608 EntityType col_type,
EntData &row_data,
611 const int nb_row_dofs = row_data.getIndices().size();
612 const int nb_col_dofs = col_data.getIndices().size();
614 if (nb_row_dofs && nb_col_dofs) {
632 mat.resize(nb_row_dofs, nb_col_dofs,
false);
634 const int nb_integration_pts =
getGaussPts().size2();
637 auto t_row_tau_base = row_data.getFTensor1N<3>();
642 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
643 const double a = vol * t_w;
645 const double K_inv = 1. /
K;
646 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
647 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
648 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
649 mat(rr, cc) += (K_inv * t_row_tau_base(
i) * t_col_tau_base(
i)) * a;
659 if (row_side != col_side || row_type != col_type) {
660 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
681 boost::shared_ptr<PreviousData> &data,
689 EntityType col_type,
EntData &row_data,
692 const int nb_row_dofs = row_data.getIndices().size();
693 const int nb_col_dofs = col_data.getIndices().size();
695 if (nb_row_dofs && nb_col_dofs) {
712 mat.resize(nb_row_dofs, nb_col_dofs,
false);
714 const int nb_integration_pts =
getGaussPts().size2();
716 auto t_row_tau_base = row_data.getFTensor1N<3>();
718 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 2>();
720 auto t_flux_value = getFTensor1FromMat<3>(
commonData->flux_values);
723 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
724 const double a = vol * t_w;
726 const double K_inv = 1. /
K;
727 const double K_diff =
B;
729 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
730 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
731 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
732 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1);
733 mat(rr, cc) += (-(t_row_tau_base(
i) * t_flux_value(
i) * K_inv *
734 K_inv * K_diff * t_col_v_base) -
735 (div_row_base * t_col_v_base)) *
767 EntityType col_type,
EntData &row_data,
770 const int nb_row_dofs = row_data.getIndices().size();
771 const int nb_col_dofs = col_data.getIndices().size();
773 if (nb_row_dofs && nb_col_dofs) {
774 mat.resize(nb_row_dofs, nb_col_dofs,
false);
776 const int nb_integration_pts =
getGaussPts().size2();
778 auto t_row_v_base = row_data.getFTensor0N();
781 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
782 const double a = vol * t_w;
783 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
784 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 2>(gg, 0);
785 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
786 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1);
787 mat(rr, cc) += (t_row_v_base * div_col_base) * a;
807 typedef boost::function<double(
const double,
const double)>
FUval;
809 boost::shared_ptr<PreviousData> &datau,
810 boost::shared_ptr<PreviousData> &datav,
819 boost::shared_ptr<PreviousData> &
dataU;
820 boost::shared_ptr<PreviousData> &
dataV;
825 EntityType col_type,
EntData &row_data,
829 const int nb_row_dofs = row_data.getIndices().size();
830 const int nb_col_dofs = col_data.getIndices().size();
831 if (nb_row_dofs && nb_col_dofs) {
833 mat.resize(nb_row_dofs, nb_col_dofs,
false);
835 const int nb_integration_pts =
getGaussPts().size2();
837 auto t_row_v_base = row_data.getFTensor0N();
845 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
846 const double a = vol * t_w;
847 double dfu =
DRhs_u(t_u_value, t_v_value);
848 for (
int rr = 0; rr != nb_row_dofs; ++rr) {
849 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
850 for (
int cc = 0; cc != nb_col_dofs; ++cc) {
851 mat(rr, cc) += ((ts_a - 0*dfu) * t_row_v_base * t_col_v_base) * a;
864 if (row_side != col_side || row_type != col_type) {
865 transMat.resize(nb_col_dofs, nb_row_dofs,
false);
973 Monitor(MPI_Comm &comm,
const int &rank, SmartPetscObj<DM> &dm,
974 boost::shared_ptr<PostProc> &post_proc)
983 "out_level_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
990 SmartPetscObj<DM>
dM;
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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)
Executes FEMethod for finite elements in DM.
DataForcesAndSourcesCore::EntData EntData
const int save_every_nth_step
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::vector< double, DoubleAllocator > VectorDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
const double r
rate factor
DataForcesAndSourcesCore::EntData EntData
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
MoFEMErrorCode operator()()
MoFEMErrorCode postProcess()
MoFEMErrorCode preProcess()
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.
default operator for EDGE element
default operator for TRI element
double getMeasure()
get measure of element
auto getFTensor1CoordsAtGaussPts()
get coordinates at Gauss pts.
Face finite element switched.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
friend class UserDataOperator
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
[Push operators to pipeline]