The example shows how to solve the linear elastic problem.
using namespace boost::numeric;
static char help[] =
"-order approximation order\n"
"\n";
double yOung;
double pOisson;
symm) {
pOisson = 0.1;
yOung = 10;
const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
tD(0, 0, 0, 0) = 1 - pOisson;
tD(1, 1, 1, 1) = 1 - pOisson;
tD(2, 2, 2, 2) = 1 - pOisson;
tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
tD(0, 0, 1, 1) = pOisson;
tD(1, 1, 0, 0) = pOisson;
tD(0, 0, 2, 2) = pOisson;
tD(2, 2, 0, 0) = pOisson;
tD(1, 1, 2, 2) = pOisson;
tD(2, 2, 1, 1) = pOisson;
tD(
i,
j,
k,
l) *= coefficient;
}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
if (!nbRows)
if (!nbCols)
K.resize(nbRows, nbCols,
false);
nbIntegrationPts = getGaussPts().size2();
if (row_side == col_side && row_type == col_type) {
isDiag = true;
} else {
isDiag = false;
}
CHKERR iNtegrate(row_data, col_data);
CHKERR aSsemble(row_data, col_data);
}
protected:
int nbRows;
int nbCols;
int nbIntegrationPts;
bool isDiag;
&
m(
r + 0,
c + 0), &
m(
r + 0,
c + 1), &
m(
r + 0,
c + 2),
&
m(
r + 1,
c + 0), &
m(
r + 1,
c + 1), &
m(
r + 1,
c + 2),
&
m(
r + 2,
c + 0), &
m(
r + 2,
c + 1), &
m(
r + 2,
c + 2));
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
const double a = t_w * vol;
for (int rr = 0; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
for (int cc = 0; cc != nbCols / 3; ++cc) {
a * (tD(
i,
j,
k,
l) * (t_row_diff_base(
j) * t_col_diff_base(
l)));
++t_col_diff_base;
++t_m;
}
++t_row_diff_base;
}
++t_w;
}
}
const int *row_indices = &*row_data.
getIndices().data().begin();
const int *col_indices = &*col_data.
getIndices().data().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
&*
K.data().begin(), ADD_VALUES);
if (!isDiag && sYmm) {
&*
K.data().begin(), ADD_VALUES);
}
}
};
double pressureVal;
pressureVal(pressure_val) {}
if (nb_dofs == 0)
nF.resize(nb_dofs, false);
nF.clear();
const int nb_gauss_pts = data.
getN().size1();
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
&nF[2]);
for (int bb = 0; bb != nb_dofs / 3; ++bb) {
t_nf(
i) += (
w * pressureVal * t_base) * t_normal(
i);
++t_nf;
++t_base;
}
++t_w;
}
&nF[0], ADD_VALUES);
}
};
Range fixFaces, fixNodes, fixSecondNode;
const Range &fix_second_node)
fixSecondNode(fix_second_node) {
}
std::set<int> set_fix_dofs;
if (dit->get()->getDofCoeffIdx() == 2) {
if (fixFaces.find(dit->get()->getEnt()) != fixFaces.end()) {
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
if (fixSecondNode.find(dit->get()->getEnt()) != fixSecondNode.end()) {
if (dit->get()->getDofCoeffIdx() == 1) {
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
if (fixNodes.find(dit->get()->getEnt()) != fixNodes.end()) {
set_fix_dofs.insert(dit->get()->getPetscGlobalDofIdx());
}
}
std::vector<int> fix_dofs(set_fix_dofs.size());
std::copy(set_fix_dofs.begin(), set_fix_dofs.end(), fix_dofs.begin());
CHKERR MatAssemblyBegin(ksp_B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(ksp_B, MAT_FINAL_ASSEMBLY);
CHKERR VecAssemblyBegin(ksp_f);
CHKERR VecDuplicate(ksp_f, &x);
CHKERR MatZeroRowsColumns(ksp_B, fix_dofs.size(), &fix_dofs[0], 1, x,
ksp_f);
}
};
int operator()(int, int, int p) const { return 2 * (p - 1); }
};
int operator()(int, int, int p) const { return p; }
};
int main(
int argc,
char *argv[]) {
const string default_options = "-ksp_type fgmres \n"
"-pc_type lu \n"
"-pc_factor_mat_solver_type mumps \n"
"-mat_mumps_icntl_20 0 \n"
"-ksp_monitor\n";
string param_file = "param_file.petsc";
if (!static_cast<bool>(ifstream(param_file))) {
std::ofstream file(param_file.c_str(), std::ios::ate);
if (file.is_open()) {
file << default_options;
file.close();
}
}
try {
PetscBool flg_test = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"SimpleElasticProblem",
"none");
PETSC_NULL);
CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
Range fix_faces, pressure_faces, fix_nodes, fix_second_node;
enum MyBcTypes {
FIX_BRICK_FACES = 1,
FIX_NODES = 2,
BRICK_PRESSURE_FACES = 3,
FIX_NODES_Y_DIR = 4
};
int id =
bit->getMeshsetId();
if (id == FIX_BRICK_FACES) {
fix_faces, true);
CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 0,
false, adj_ents,
moab::Interface::UNION);
CHKERR m_field.
get_moab().get_adjacencies(fix_faces, 1,
false, adj_ents,
moab::Interface::UNION);
fix_faces.merge(adj_ents);
} else if (id == FIX_NODES) {
fix_nodes, true);
} else if (id == BRICK_PRESSURE_FACES) {
meshset, 2, pressure_faces, true);
} else if (id ==
FIX_NODES_Y_DIR) {
meshset, 0, fix_second_node, true);
} else {
}
}
3);
DM dm;
"PRESSURE");
boost::shared_ptr<VolumeElementForcesAndSourcesCore> elastic_fe(
boost::shared_ptr<FaceElementForcesAndSourcesCore> pressure_fe(
elastic_fe->getRuleHook =
VolRule();
elastic_fe->getOpPtrVector().push_back(
new OpK());
pressure_fe->getOpPtrVector().push_back(
new OpPressure());
boost::shared_ptr<FEMethod> fix_dofs_fe(
boost::shared_ptr<FEMethod> null_fe;
null_fe);
CHKERR DMCreateGlobalVector(dm, &x);
elastic_fe);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
VecView(x, PETSC_VIEWER_STDOUT_WORLD);
auto get_post_proc_ele = [&]() {
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto post_proc_ele = boost::make_shared<
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getOpPtrVector().push_back(
grad_ptr));
post_proc_ele->getOpPtrVector().push_back(
post_proc_ele->getPostProcMesh(), post_proc_ele->getMapGaussPts(),
{},
{{"U", u_ptr}},
{{"GRAD", grad_ptr}},
{})
);
return post_proc_ele;
};
if (auto post_proc = get_post_proc_ele()) {
post_proc);
CHKERR post_proc->writeFile(
"out.h5m");
}
{
if (flg_test == PETSC_TRUE) {
const double x_vec_norm_const = 0.4;
double norm_check;
CHKERR VecNorm(x, NORM_INFINITY, &norm_check);
if (std::abs(norm_check - x_vec_norm_const) / x_vec_norm_const >
1.e-10) {
"test failed (nrm 2 %6.4e)", norm_check);
}
}
}
}
return 0;
}