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;
999 #endif //__ELECPHYSOPERATORS_HPP__