64 {
65
67
68 try {
69
70 moab::Core mb_instance;
71 moab::Interface &moab = mb_instance;
72
73
76
77
78 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
79 const char *list_bases[] = {"ainsworth", "demkowicz"};
80 PetscBool flg;
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)
93
94
95 DMType dm_name = "DMMOFEM";
97
98
99 auto core_log = logging::core::get();
100 core_log->add_sink(
104
105
107
108
110
112
119
121
122
124
126
127 auto get_skin = [&]() {
130 body_ents);
133 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
134 return skin_ents;
135 };
136
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;
144 };
145
146 auto boundary_ents = filter_true_skin(get_skin());
149
150
152
154 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"SYM",
156 CHKERR bc_mng->removeBlockDOFsOnEntities(
158
159 auto project_ho_geometry = [&]() {
161 return m_field.
loop_dofs(
"GEOMETRY", ent_method);
162 };
163 CHKERR project_ho_geometry();
164
165
167
169
170
171 pip_mng->getOpDomainLhsPipeline().clear();
172 pip_mng->getOpDomainRhsPipeline().clear();
173
174
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);
179
181 {{"U", {-1, 1}}, {"SIGMA", {-1, 1}}});
183 SCATTER_REVERSE);
184
185
186
187 auto jacobian = [&](
double r) {
189 return 2 * M_PI * r;
190 else
191 return 1.;
192 };
193
195 return jacobian(r);
196 };
198 return -jacobian(r);
199 };
200
201 auto post_proc = [&](auto dm, auto f_res, auto out_name) {
203 auto post_proc_fe =
204 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
205
207
208 auto sigma_ptr = boost::make_shared<MatrixDouble>();
209 post_proc_fe->getOpPtrVector().push_back(
211 sigma_ptr));
212 auto u_ptr = boost::make_shared<MatrixDouble>();
213 post_proc_fe->getOpPtrVector().push_back(
215
216 post_proc_fe->getOpPtrVector().push_back(
217
219
220 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
221 {}, {{"U", u_ptr}}, {{"SIGMA", sigma_ptr}}, {})
222
223 );
224
226 SCATTER_REVERSE);
228 post_proc_fe);
229 post_proc_fe->writeFile(out_name);
231 };
232
233 auto test_consistency_of_domain_and_bdy_integrals = [&]() {
237 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM, COORD_TYPE>;
244
247
252
253 auto ops_rhs_interior = [&](auto &pip) {
256 "GEOMETRY");
257 auto u_ptr = boost::make_shared<MatrixDouble>();
258 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
261 "U", grad_u_ptr));
262
263 pip.push_back(
new OpMixDivURhs(
"SIGMA", u_ptr, beta_domain));
264 pip.push_back(
265 new OpMixLambdaGradURhs("SIGMA", grad_u_ptr, beta_domain));
266
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));
274 pip.push_back(
277
279 };
280
281 auto ops_rhs_boundary = [&](auto &pip) {
284 "GEOMETRY");
285 auto u_ptr = boost::make_shared<MatrixDouble>();
287 auto traction_ptr = boost::make_shared<MatrixDouble>();
289 "SIGMA", traction_ptr));
290
291
292
294 pip.push_back(new OpMixNormalLambdaURhs("SIGMA", u_ptr, beta_bdy));
295 pip.push_back(new OpUTimeTractionRhs("U", traction_ptr, beta_bdy));
296
298 };
299
300 CHKERR ops_rhs_interior(pip_mng->getOpDomainRhsPipeline());
301 CHKERR ops_rhs_boundary(pip_mng->getOpBoundaryRhsPipeline());
302
304 pip_mng->getDomainRhsFE()->f =
f;
305 pip_mng->getBoundaryRhsFE()->f =
f;
306
308
310 simple->getDomainFEName(),
311 pip_mng->getDomainRhsFE());
313 simple->getBoundaryFEName(),
314 pip_mng->getBoundaryRhsFE());
315
316 CHKERR VecAssemblyBegin(f);
318
319 double f_nrm2;
320 CHKERR VecNorm(f, NORM_2, &f_nrm2);
321
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");
327 }
328
330 };
331
332
333
334 auto test_lhs_ops = [&]() {
340
341 auto op_lhs_domain = [&](auto &pip) {
344 "GEOMETRY");
345
346 auto unity = []() { return 1; };
347 pip.push_back(
348 new OpMixDivULhs(
"SIGMA",
"U", unity, beta_domain,
true,
false));
349 pip.push_back(
352 };
353
354 CHKERR op_lhs_domain(pip_mng->getOpDomainLhsPipeline());
355
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(),
363
364
365
366 double fnorm_res;
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)
371
373 };
374
375 CHKERR test_consistency_of_domain_and_bdy_integrals();
377}
379
381}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#define CATCH_ERRORS
Catch errors.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr)
Generate random fileds.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
constexpr FieldSpace HDIV_SPACE
constexpr CoordinateTypes COORD_TYPE