static char help[] =
"...\n\n";
};
struct OpAssembleMass :
OpEle {
OpAssembleMass(boost::shared_ptr<CommonData> &data)
sYmm = true;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
const int nb_integration_pts = getGaussPts().size2();
mat.resize(nb_row_dofs, nb_col_dofs, false);
mat.clear();
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = t_w * vol;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; ++cc) {
mat(rr, cc) +=
a * t_row_base * t_col_base;
++t_col_base;
}
++t_row_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transMat.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transMat) = trans(mat);
ADD_VALUES);
}
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
struct OpAssembleSlowRhs :
OpEle {
OpAssembleSlowRhs(boost::shared_ptr<CommonData> &data)
if (nb_dofs) {
vecF.resize(nb_dofs, false);
vecF.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
const double f =
a *
r * t_val * (1 - t_val /
k);
for (int rr = 0; rr != nb_dofs; ++rr) {
const double b =
f * t_base;
vecF[rr] += b;
++t_base;
}
++t_val;
++t_w;
}
CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
ADD_VALUES);
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
template <
int DIM>
struct OpAssembleStiffRhs :
OpEle {
OpAssembleStiffRhs(boost::shared_ptr<CommonData> &data)
if (nb_dofs) {
vecF.resize(nb_dofs, false);
vecF.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_grad = getFTensor1FromMat<DIM>(commonData->grad);
auto t_w = getFTensor0IntegrationWeight();
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
for (int rr = 0; rr != nb_dofs; ++rr) {
vecF[rr] +=
a * (t_base * t_dot_val +
D * t_diff_base(
i) * t_grad(
i));
++t_diff_base;
++t_base;
}
++t_dot_val;
++t_grad;
++t_w;
}
CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
PETSC_TRUE);
ADD_VALUES);
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
template <
int DIM>
struct OpAssembleStiffLhs :
OpEle {
OpAssembleStiffLhs(boost::shared_ptr<CommonData> &data)
sYmm = true;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntData &row_data,
const int nb_row_dofs = row_data.
getIndices().size();
const int nb_col_dofs = col_data.
getIndices().size();
if (nb_row_dofs && nb_col_dofs) {
mat.resize(nb_row_dofs, nb_col_dofs, false);
mat.clear();
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
const double ts_a = getFEMethod()->ts_a;
const double vol = getMeasure();
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double a = vol * t_w;
for (int rr = 0; rr != nb_row_dofs; ++rr) {
for (int cc = 0; cc != nb_col_dofs; ++cc) {
mat(rr, cc) +=
a * (t_row_base * t_col_base * ts_a +
D * t_row_diff_base(
i) * t_col_diff_base(
i));
++t_col_base;
++t_col_diff_base;
}
++t_row_base;
++t_row_diff_base;
}
++t_w;
}
ADD_VALUES);
if (row_side != col_side || row_type != col_type) {
transMat.resize(nb_col_dofs, nb_row_dofs, false);
noalias(transMat) = trans(mat);
&transMat(0, 0), ADD_VALUES);
}
}
}
private:
boost::shared_ptr<CommonData> commonData;
};
: dM(dm), postProc(post_proc){};
"out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
};
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
try {
DMType dm_name = "DMMOFEM";
1);
boost::shared_ptr<CommonData> data(
new CommonData());
auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
boost::shared_ptr<Ele> vol_ele_slow_rhs(
new Ele(m_field));
boost::shared_ptr<Ele> vol_ele_stiff_rhs(
new Ele(m_field));
boost::shared_ptr<Ele> vol_ele_stiff_lhs(
new Ele(m_field));
vol_ele_slow_rhs->getOpPtrVector().push_back(
vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
auto solve_for_g = [&]() {
if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
SCATTER_REVERSE);
CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
*vol_ele_slow_rhs->vecAssembleSwitch = false;
}
CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
vol_ele_slow_rhs->ts_F);
};
vol_ele_slow_rhs->postProcessHook = solve_for_g;
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
vol_ele_stiff_rhs->getOpPtrVector().push_back(
vol_ele_stiff_rhs->getOpPtrVector().push_back(
vol_ele_stiff_rhs->getOpPtrVector().push_back(
vol_ele_stiff_rhs->getOpPtrVector().push_back(
vol_ele_stiff_rhs->getOpPtrVector().push_back(
vol_ele_stiff_rhs->getOpPtrVector().push_back(
new OpAssembleStiffRhs<2>(data));
vol_ele_stiff_lhs->getOpPtrVector().push_back(
vol_ele_stiff_lhs->getOpPtrVector().push_back(
vol_ele_stiff_lhs->getOpPtrVector().push_back(
vol_ele_stiff_lhs->getOpPtrVector().push_back(
new OpAssembleStiffLhs<2>(data));
auto vol_rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
vol_ele_slow_rhs->getRuleHook = vol_rule;
vol_ele_stiff_rhs->getRuleHook = vol_rule;
vol_ele_stiff_lhs->getRuleHook = vol_rule;
auto post_proc = boost::make_shared<PostProcEle>(m_field);
boost::shared_ptr<ForcesAndSourcesCore> null;
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc->getOpPtrVector().push_back(
post_proc->getOpPtrVector().push_back(
post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
{{"u", u_ptr}},
{},
{},
{}
)
);
auto dm = simple_interface->
getDM();
if (!inner_surface.empty()) {
Range inner_surface_verts;
CHKERR moab.get_connectivity(inner_surface, inner_surface_verts,
false);
u0, MBVERTEX, inner_surface_verts,
"u");
}
}
ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &edges_part);
CHKERR moab.get_connectivity(edges_part, edges_verts,
false);
unite(edges_verts, edges_part));
CHKERR MatZeroEntries(data->M);
boost::shared_ptr<Ele> vol_mass_ele(
new Ele(m_field));
vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
vol_mass_ele);
CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
CHKERR KSPSetOperators(data->ksp, data->M, data->M);
CHKERR KSPSetFromOptions(data->ksp);
CHKERR TSSetType(ts, TSARKIMEX);
CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
vol_ele_stiff_lhs, null, null);
vol_ele_stiff_rhs, null, null);
vol_ele_slow_rhs, null, null);
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, post_proc));
monitor_ptr, null, null);
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
}
return 0;
}