78 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
79 const char *list_bases[] = {
"ainsworth",
"demkowicz"};
81 PetscInt choice_base_value = AINSWORTH;
83 LASBASETOP, &choice_base_value, &flg);
84 if (flg != PETSC_TRUE)
87 if (choice_base_value == AINSWORTH)
89 else if (choice_base_value == DEMKOWICZ)
95 DMType dm_name =
"DMMOFEM";
99 auto core_log = logging::core::get();
101 LogManager::createSink(LogManager::getStrmWorld(),
"ATOM"));
102 LogManager::setLog(
"ATOM");
127 auto get_skin = [&]() {
133 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
137 auto filter_true_skin = [&](
auto skin) {
139 ParallelComm *pcomm =
141 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
142 PSTATUS_NOT, -1, &boundary_ents);
143 return boundary_ents;
146 auto boundary_ents = filter_true_skin(get_skin());
154 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM",
156 CHKERR bc_mng->removeBlockDOFsOnEntities(
159 auto project_ho_geometry = [&]() {
161 return m_field.
loop_dofs(
"GEOMETRY", ent_method);
163 CHKERR project_ho_geometry();
171 pip_mng->getOpDomainLhsPipeline().clear();
172 pip_mng->getOpDomainRhsPipeline().clear();
175 auto rule = [&](
int,
int,
int p) {
return 2 * p + 2; };
176 CHKERR pip_mng->setDomainRhsIntegrationRule(rule);
177 CHKERR pip_mng->setDomainLhsIntegrationRule(rule);
178 CHKERR pip_mng->setBoundaryRhsIntegrationRule(rule);
181 {{
"U", {-1, 1}}, {
"SIGMA", {-1, 1}}});
187 auto jacobian = [&](
double r) {
201 auto post_proc = [&](
auto dm,
auto f_res,
auto out_name) {
204 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
208 auto sigma_ptr = boost::make_shared<MatrixDouble>();
209 post_proc_fe->getOpPtrVector().push_back(
212 auto u_ptr = boost::make_shared<MatrixDouble>();
213 post_proc_fe->getOpPtrVector().push_back(
216 post_proc_fe->getOpPtrVector().push_back(
220 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
221 {}, {{
"U", u_ptr}}, {{
"SIGMA", sigma_ptr}}, {})
229 post_proc_fe->writeFile(out_name);
233 auto test_consistency_of_domain_and_bdy_integrals = [&]() {
237 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM, COORD_TYPE>;
253 auto ops_rhs_interior = [&](
auto &pip) {
257 auto u_ptr = boost::make_shared<MatrixDouble>();
258 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
263 pip.push_back(
new OpMixDivURhs(
"SIGMA", u_ptr, beta_domain));
265 new OpMixLambdaGradURhs(
"SIGMA", grad_u_ptr, beta_domain));
267 auto sigma_ptr = boost::make_shared<MatrixDouble>();
268 auto sigma_div_ptr = boost::make_shared<MatrixDouble>();
271 "SIGMA", sigma_div_ptr));
273 "SIGMA", sigma_ptr));
275 new OpMixUTimesDivLambdaRhs(
"U", sigma_div_ptr, beta_domain));
276 pip.push_back(
new OpMixUTimesLambdaRhs(
"U", sigma_ptr, beta_domain));
281 auto ops_rhs_boundary = [&](
auto &pip) {
285 auto u_ptr = boost::make_shared<MatrixDouble>();
287 auto traction_ptr = boost::make_shared<MatrixDouble>();
289 "SIGMA", traction_ptr));
294 pip.push_back(
new OpMixNormalLambdaURhs(
"SIGMA", u_ptr, beta_bdy));
295 pip.push_back(
new OpUTimeTractionRhs(
"U", traction_ptr, beta_bdy));
300 CHKERR ops_rhs_interior(pip_mng->getOpDomainRhsPipeline());
301 CHKERR ops_rhs_boundary(pip_mng->getOpBoundaryRhsPipeline());
304 pip_mng->getDomainRhsFE()->f =
f;
305 pip_mng->getBoundaryRhsFE()->f =
f;
310 simple->getDomainFEName(),
311 pip_mng->getDomainRhsFE());
313 simple->getBoundaryFEName(),
314 pip_mng->getBoundaryRhsFE());
320 CHKERR VecNorm(
f, NORM_2, &f_nrm2);
322 MOFEM_LOG(
"ATOM", Sev::inform) <<
"f_norm2 = " << f_nrm2;
323 if (std::fabs(f_nrm2) > 1e-10) {
325 "tensor_divergence_operator_res_vec.h5m");
334 auto test_lhs_ops = [&]() {
341 auto op_lhs_domain = [&](
auto &pip) {
346 auto unity = []() {
return 1; };
348 new OpMixDivULhs(
"SIGMA",
"U", unity, beta_domain,
true,
false));
350 new OpLambdaGraULhs(
"SIGMA",
"U", unity, beta_domain,
true,
false));
354 CHKERR op_lhs_domain(pip_mng->getOpDomainLhsPipeline());
356 auto diff_x = opt->setRandomFields(
simple->getDM(),
357 {{
"U", {-1, 1}}, {
"SIGMA", {-1, 1}}});
358 constexpr
double eps = 1e-5;
359 auto diff_res = opt->checkCentralFiniteDifference(
360 simple->getDM(),
simple->getDomainFEName(), pip_mng->getDomainRhsFE(),
367 CHKERR VecNorm(diff_res, NORM_2, &fnorm_res);
368 MOFEM_LOG_C(
"ATOM", Sev::inform,
"Test Lhs OPs %3.4e", fnorm_res);
369 if (std::fabs(fnorm_res) > 1e-8)
375 CHKERR test_consistency_of_domain_and_bdy_integrals();