v0.13.1
Functions | Variables
hcurl_continuity_check.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
static const double face_coords [4][9]
 
static const double edge_coords [6][6]
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 22 of file hcurl_continuity_check.cpp.

22 {
23
24 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
25
26 try {
27
28 moab::Core mb_instance;
29 moab::Interface &moab = mb_instance;
30 int rank;
31 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
32
33 PetscBool flg = PETSC_TRUE;
34 char mesh_file_name[255];
35#if PETSC_VERSION_GE(3, 6, 4)
36 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
37 255, &flg);
38
39#else
40 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
41 mesh_file_name, 255, &flg);
42
43#endif
44 if (flg != PETSC_TRUE) {
45 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
46 }
47 const char *option;
48 option = "";
49 CHKERR moab.load_file(mesh_file_name, 0, option);
50
51 // Select base
52 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
53
54 const char *list_bases[] = {"ainsworth", "demkowicz"};
55
56 PetscInt choice_base_value = AINSWORTH;
57 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
58 LASBASETOP, &choice_base_value, &flg);
59
60 if (flg != PETSC_TRUE) {
61 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
62 }
63
65 if (choice_base_value == AINSWORTH) {
67 } else if (choice_base_value == DEMKOWICZ) {
69 }
70
71 // Create MoFEM (Joseph) database
72 MoFEM::Core core(moab);
73 MoFEM::Interface &m_field = core;
74
75 // set entities bit level
76 BitRefLevel bit_level0;
77 bit_level0.set(0);
78 EntityHandle meshset_level0;
79 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
80
81 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
82 0, 3, bit_level0);
83
84 // Fields
85 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
86 3);
87
88 CHKERR m_field.add_field("HCURL", HCURL, base, 1);
89
90 // FE
91 CHKERR m_field.add_finite_element("TET_FE");
92 CHKERR m_field.add_finite_element("TRI_FE");
93 CHKERR m_field.add_finite_element("SKIN_FE");
94 CHKERR m_field.add_finite_element("EDGE_FE");
95
96 // Define rows/cols and element data
97 CHKERR m_field.modify_finite_element_add_field_row("TET_FE", "HCURL");
98 CHKERR m_field.modify_finite_element_add_field_col("TET_FE", "HCURL");
99 CHKERR m_field.modify_finite_element_add_field_data("TET_FE", "HCURL");
101 "MESH_NODE_POSITIONS");
102
103 CHKERR m_field.modify_finite_element_add_field_row("SKIN_FE", "HCURL");
104 CHKERR m_field.modify_finite_element_add_field_col("SKIN_FE", "HCURL");
105 CHKERR m_field.modify_finite_element_add_field_data("SKIN_FE", "HCURL");
107 "MESH_NODE_POSITIONS");
108
109 CHKERR m_field.modify_finite_element_add_field_row("TRI_FE", "HCURL");
110 CHKERR m_field.modify_finite_element_add_field_col("TRI_FE", "HCURL");
111 CHKERR m_field.modify_finite_element_add_field_data("TRI_FE", "HCURL");
113 "MESH_NODE_POSITIONS");
114
115 CHKERR m_field.modify_finite_element_add_field_row("EDGE_FE", "HCURL");
116 CHKERR m_field.modify_finite_element_add_field_col("EDGE_FE", "HCURL");
117 CHKERR m_field.modify_finite_element_add_field_data("EDGE_FE", "HCURL");
119 "MESH_NODE_POSITIONS");
120
121 // Problem
122 CHKERR m_field.add_problem("TEST_PROBLEM");
123
124 // set finite elements for problem
125 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TET_FE");
126 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "SKIN_FE");
127 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TRI_FE");
128 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE_FE");
129
130 // set refinement level for problem
131 CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
132
133 // meshset consisting all entities in mesh
134 EntityHandle root_set = moab.get_root_set();
135 // add entities to field
136 CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "HCURL");
137
138 // add entities to finite element
139 CHKERR
140 m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "TET_FE");
141
142 Range tets;
143 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
144 BitRefLevel().set(0), BitRefLevel().set(), MBTET, tets);
145
146 Skinner skin(&moab);
147 Range skin_faces; // skin faces from 3d ents
148 CHKERR skin.find_skin(0, tets, false, skin_faces);
149
150 CHKERR m_field.add_ents_to_finite_element_by_type(skin_faces, MBTRI,
151 "SKIN_FE");
152
153 Range faces;
154 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
155 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, faces);
156
157 faces = subtract(faces, skin_faces);
158 CHKERR m_field.add_ents_to_finite_element_by_type(faces, MBTRI, "TRI_FE");
159
160 Range edges;
161 CHKERR moab.get_adjacencies(faces, 1, false, edges, moab::Interface::UNION);
162
163 CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE, "EDGE_FE");
164
165 Range skin_edges;
166 CHKERR moab.get_adjacencies(skin_faces, 1, false, skin_edges,
167 moab::Interface::UNION);
168
169 CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE, "EDGE_FE");
170
171 // set app. order
172 int order = 4;
173 CHKERR m_field.set_field_order(root_set, MBTET, "HCURL", order);
174 CHKERR m_field.set_field_order(root_set, MBTRI, "HCURL", order);
175 CHKERR m_field.set_field_order(root_set, MBEDGE, "HCURL", order);
176
177 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
178 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
179 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
180 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
181 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
182
183 /****/
184 // build database
185 // build field
186 CHKERR m_field.build_fields();
187
188 // build finite elemnts
190
191 // build adjacencies
192 CHKERR m_field.build_adjacencies(bit_level0);
193
194 // build problem
195 ProblemsManager *prb_mng_ptr;
196 CHKERR m_field.getInterface(prb_mng_ptr);
197 CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
198
199 // project geometry form 10 node tets on higher order approx. functions
200 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
201 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
202
203 /****/
204 // mesh partitioning
205 // partition
206 CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
207 CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
208
209 // what are ghost nodes, see Petsc Manual
210 CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
211
212 Vec v;
213 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
214 ROW, &v);
215 CHKERR VecSetRandom(v, PETSC_NULL);
216
217 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
218 "TEST_PROBLEM", ROW, v, INSERT_VALUES, SCATTER_REVERSE);
219
220 CHKERR VecDestroy(&v);
221
222 typedef tee_device<std::ostream, std::ofstream> TeeDevice;
223 typedef stream<TeeDevice> TeeStream;
224 std::ofstream ofs("forces_and_sources_hdiv_continuity_check.txt");
225 TeeDevice my_tee(std::cout, ofs);
226 TeeStream my_split(my_tee);
227
228 struct OpTetFluxes
230
231 MoFEM::Interface &mField;
232 Tag tH;
233
234 OpTetFluxes(MoFEM::Interface &m_field, Tag th)
236 "HCURL", UserDataOperator::OPROW),
237 mField(m_field), tH(th) {}
238
239 MoFEMErrorCode doWork(int side, EntityType type,
242
243 if (data.getFieldData().size() == 0)
245
246 if (type == MBTRI) {
247
248 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
249 getNumeredEntFiniteElementPtr();
250 SideNumber_multiIndex &side_table = mofem_fe->getSideNumberTable();
251 EntityHandle face = side_table.get<1>()
252 .find(boost::make_tuple(type, side))
253 ->get()
254 ->ent;
255 int sense = side_table.get<1>()
256 .find(boost::make_tuple(type, side))
257 ->get()
258 ->sense;
259
260 // cerr << data.getVectorN() << endl;
261
262 VectorDouble t(3, 0);
263 int nb_dofs = data.getN().size2() / 3;
264 for (int dd = 0; dd < nb_dofs; dd++) {
265 for (int ddd = 0; ddd < 3; ddd++) {
266 t(ddd) +=
267 data.getVectorN<3>(side)(dd, ddd) * data.getFieldData()[dd];
268 }
269 }
270
271 double *t_ptr;
272 CHKERR mField.get_moab().tag_get_by_ptr(tH, &face, 1,
273 (const void **)&t_ptr);
274
275 for (int dd = 0; dd < 3; dd++) {
276 t_ptr[dd] += sense * t[dd];
277 }
278 }
279
280 if (type == MBEDGE) {
281
282 boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
283 getNumeredEntFiniteElementPtr();
284 SideNumber_multiIndex &side_table = mofem_fe->getSideNumberTable();
285 EntityHandle edge = side_table.get<1>()
286 .find(boost::make_tuple(type, side))
287 ->get()
288 ->ent;
289 Range adj_tets;
290 CHKERR mField.get_moab().get_adjacencies(&edge, 1, 3, false, adj_tets,
291 moab::Interface::UNION);
292 const int nb_adj_tets = adj_tets.size();
293
294 VectorDouble t(3, 0);
295 int nb_dofs = data.getN().size2() / 3;
296 for (int dd = 0; dd < nb_dofs; dd++) {
297 for (int ddd = 0; ddd < 3; ddd++) {
298 t(ddd) += data.getVectorN<3>(4 + side)(dd, ddd) *
299 data.getFieldData()[dd];
300 }
301 }
302
303 double *t_ptr;
304 CHKERR mField.get_moab().tag_get_by_ptr(tH, &edge, 1,
305 (const void **)&t_ptr);
306
307 for (int dd = 0; dd < 3; dd++) {
308 t_ptr[dd] += t[dd] / nb_adj_tets;
309 }
310 }
311
313 }
314 };
315
316 struct MyTetFE : public VolumeElementForcesAndSourcesCore {
317
318 MyTetFE(MoFEM::Interface &m_field)
320 int getRule(int order) { return -1; };
321
322 MatrixDouble N_tri;
324
326
327 try {
328
329 N_tri.resize(1, 3);
330 CHKERR ShapeMBTRI(&N_tri(0, 0), G_TRI_X1, G_TRI_Y1, 1);
331
332 gaussPts.resize(4, 4 + 6);
333 int ff = 0;
334 for (; ff < 4; ff++) {
335 int dd = 0;
336 for (; dd < 3; dd++) {
337 gaussPts(dd, ff) =
338 cblas_ddot(3, &N_tri(0, 0), 1, &face_coords[ff][dd], 3);
339 }
340 gaussPts(3, ff) = G_TRI_W1[0];
341 }
342
343 int ee = 0;
344 for (; ee < 6; ee++) {
345 int dd = 0;
346 for (; dd < 3; dd++) {
347 gaussPts(dd, 4 + ee) =
348 (edge_coords[ee][0 + dd] + edge_coords[ee][3 + dd]) * 0.5;
349 }
350 gaussPts(3, 4 + ee) = 1;
351 }
352
353 // std::cerr << gaussPts << std::endl;
354
355 } catch (std::exception &ex) {
356 std::ostringstream ss;
357 ss << "throw in method: " << ex.what() << " at line " << __LINE__
358 << " in file " << __FILE__;
359 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
360 }
361
363 }
364 };
365
366 struct OpFacesSkinFluxes
368
369 MoFEM::Interface &mField;
370 Tag tH1, tH2;
371 TeeStream &mySplit;
372
373 OpFacesSkinFluxes(MoFEM::Interface &m_field, Tag th1, Tag th2,
374 TeeStream &my_split)
376 "HCURL", UserDataOperator::OPROW),
377 mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
378
379 MoFEMErrorCode doWork(int side, EntityType type,
382
383 if (type != MBTRI)
385
386 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
387
388 double *t_ptr;
389 CHKERR mField.get_moab().tag_get_by_ptr(tH1, &face, 1,
390 (const void **)&t_ptr);
391
392 double *tn_ptr;
393 CHKERR mField.get_moab().tag_get_by_ptr(tH2, &face, 1,
394 (const void **)&tn_ptr);
395
396 *tn_ptr = getTangent1AtGaussPts()(0, 0) * t_ptr[0] +
397 getTangent1AtGaussPts()(0, 1) * t_ptr[1] +
398 getTangent1AtGaussPts()(0, 2) * t_ptr[2] +
399 getTangent2AtGaussPts()(0, 0) * t_ptr[0] +
400 getTangent2AtGaussPts()(0, 1) * t_ptr[1] +
401 getTangent2AtGaussPts()(0, 2) * t_ptr[2];
402
403 int nb_dofs = data.getN().size2() / 3;
404 int dd = 0;
405 for (; dd < nb_dofs; dd++) {
406 double val = data.getFieldData()[dd];
407 *tn_ptr +=
408 -getTangent1AtGaussPts()(0, 0) * data.getN()(0, 3 * dd + 0) *
409 val -
410 getTangent1AtGaussPts()(0, 1) * data.getN()(0, 3 * dd + 1) * val -
411 getTangent1AtGaussPts()(0, 2) * data.getN()(0, 3 * dd + 2) * val -
412 getTangent2AtGaussPts()(0, 0) * data.getN()(0, 3 * dd + 0) * val -
413 getTangent2AtGaussPts()(0, 1) * data.getN()(0, 3 * dd + 1) * val -
414 getTangent2AtGaussPts()(0, 2) * data.getN()(0, 3 * dd + 2) * val;
415 }
416
417 const double eps = 1e-8;
418 if (fabs(*tn_ptr) > eps) {
419 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
420 "HCurl continuity failed %6.4e", *tn_ptr);
421 }
422
423 mySplit.precision(5);
424 mySplit << face << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
425
427 }
428 };
429
430 struct OpFacesFluxes
432
433 MoFEM::Interface &mField;
434 Tag tH1, tH2;
435 TeeStream &mySplit;
436
437 OpFacesFluxes(MoFEM::Interface &m_field, Tag _th1, Tag _th2,
438 TeeStream &my_split)
440 "HCURL", UserDataOperator::OPROW),
441 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
442
446
447 if (type != MBTRI)
449
450 EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
451
452 double *t_ptr;
453 CHKERR mField.get_moab().tag_get_by_ptr(tH1, &face, 1,
454 (const void **)&t_ptr);
455
456 double *tn_ptr;
457 CHKERR mField.get_moab().tag_get_by_ptr(tH2, &face, 1,
458 (const void **)&tn_ptr);
459
460 *tn_ptr = getTangent1AtGaussPts()(0, 0) * t_ptr[0] +
461 getTangent1AtGaussPts()(0, 1) * t_ptr[1] +
462 getTangent1AtGaussPts()(0, 2) * t_ptr[2] +
463 getTangent2AtGaussPts()(0, 0) * t_ptr[0] +
464 getTangent2AtGaussPts()(0, 1) * t_ptr[1] +
465 getTangent2AtGaussPts()(0, 2) * t_ptr[2];
466
467 const double eps = 1e-8;
468 if (fabs(*tn_ptr) > eps) {
469 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
470 "HCurl continuity failed %6.4e", *tn_ptr);
471 }
472
473 mySplit.precision(5);
474
475 mySplit << face << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
476
478 }
479 };
480
481 struct MyTriFE : public FaceElementForcesAndSourcesCore {
482
483 MyTriFE(MoFEM::Interface &m_field)
485 int getRule(int order) { return -1; };
486
489
490 gaussPts.resize(3, 1);
491 gaussPts(0, 0) = G_TRI_X1[0];
492 gaussPts(1, 0) = G_TRI_Y1[0];
493 gaussPts(2, 0) = G_TRI_W1[0];
494
496 }
497 };
498
499 struct OpEdgesFluxes
501
502 MoFEM::Interface &mField;
503 Tag tH1, tH2;
504 TeeStream &mySplit;
505
506 OpEdgesFluxes(MoFEM::Interface &m_field, Tag _th1, Tag _th2,
507 TeeStream &my_split)
509 "HCURL", UserDataOperator::OPROW),
510 mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
511
512 MoFEMErrorCode doWork(int side, EntityType type,
515
516 if (type != MBEDGE)
518
519 EntityHandle edge = getNumeredEntFiniteElementPtr()->getEnt();
520
521 double *t_ptr;
522 CHKERR mField.get_moab().tag_get_by_ptr(tH1, &edge, 1,
523 (const void **)&t_ptr);
524
525 double *tn_ptr;
526 CHKERR mField.get_moab().tag_get_by_ptr(tH2, &edge, 1,
527 (const void **)&tn_ptr);
528
529 *tn_ptr = getTangentAtGaussPts()(0, 0) * t_ptr[0] +
530 getTangentAtGaussPts()(0, 1) * t_ptr[1] +
531 getTangentAtGaussPts()(0, 2) * t_ptr[2];
532
533 double tn = 0;
534 unsigned int nb_dofs = data.getN().size2() / 3;
535 if (nb_dofs != data.getFieldData().size()) {
536 SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
537 "Number of dofs on edge and number of base functions not "
538 "equla %d != %d",
539 nb_dofs, data.getFieldData().size());
540 }
541
542 for (unsigned int dd = 0; dd != nb_dofs; ++dd) {
543 double val = data.getFieldData()[dd];
544 tn += getTangentAtGaussPts()(0, 0) * data.getN()(0, 3 * dd + 0) *
545 val +
546 getTangentAtGaussPts()(0, 1) * data.getN()(0, 3 * dd + 1) *
547 val +
548 getTangentAtGaussPts()(0, 2) * data.getN()(0, 3 * dd + 2) *
549 val;
550 }
551
552 // mySplit << *tn_ptr << " " << tn << " " << getLength() << endl;
553 *tn_ptr -= tn;
554
555 // mySplit << getTangentAtGaussPts() << " " << getDirection() << endl;
556
557 // cerr << t_ptr[0] << " " << t_ptr[1] << " " << t_ptr[2] << endl;
558
559 const double eps = 1e-8;
560 if (fabs(*tn_ptr) > eps) {
561 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
562 "HCurl continuity failed %6.4e", *tn_ptr);
563 }
564
565 mySplit.precision(5);
566
567 mySplit << edge << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
568
570 }
571 };
572
573 struct MyEdgeFE : public EdgeElementForcesAndSourcesCore {
574
575 MyEdgeFE(MoFEM::Interface &m_field)
577 int getRule(int order) { return -1; };
578
581
582 gaussPts.resize(2, 1);
583 gaussPts(0, 0) = 0.5;
584 gaussPts(1, 0) = 1;
585
587 }
588 };
589
590 MyTetFE tet_fe(m_field);
591 MyTriFE tri_fe(m_field);
592 MyTriFE skin_fe(m_field);
593 MyEdgeFE edge_fe(m_field);
594
595 Tag th1;
596 double def_val[] = {0, 0, 0};
597 CHKERR moab.tag_get_handle("T", 3, MB_TYPE_DOUBLE, th1,
598 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
599
600 auto material_grad_mat = boost::make_shared<MatrixDouble>();
601 auto material_det_vec = boost::make_shared<VectorDouble>();
602 auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
603
604 tet_fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
605 "MESH_NODE_POSITIONS", material_grad_mat));
606 tet_fe.getOpPtrVector().push_back(new OpInvertMatrix<3>(
607 material_grad_mat, material_det_vec, material_inv_grad_mat));
608 tet_fe.getOpPtrVector().push_back(new OpSetHOWeights(material_det_vec));
609 tet_fe.getOpPtrVector().push_back(
610 new OpSetHOCovariantPiolaTransform(HCURL, material_inv_grad_mat));
611 tet_fe.getOpPtrVector().push_back(
612 new OpSetHOInvJacVectorBase(HCURL, material_inv_grad_mat));
613 tet_fe.getOpPtrVector().push_back(new OpTetFluxes(m_field, th1));
614
615 Tag th2;
616 CHKERR moab.tag_get_handle("TN", 1, MB_TYPE_DOUBLE, th2,
617 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
618
619 tri_fe.getOpPtrVector().push_back(
620 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
621 tri_fe.getOpPtrVector().push_back(
622 new OpCalculateHOCoords("MESH_NODE_POSITIONS"));
623
624 tri_fe.getOpPtrVector().push_back(
626 tri_fe.getOpPtrVector().push_back(
627 new OpFacesFluxes(m_field, th1, th2, my_split));
628
629 skin_fe.getOpPtrVector().push_back(
630 new OpGetHONormalsOnFace("MESH_NODE_POSITIONS"));
631 skin_fe.getOpPtrVector().push_back(
633 skin_fe.getOpPtrVector().push_back(
634 new OpFacesSkinFluxes(m_field, th1, th2, my_split));
635
636 edge_fe.getOpPtrVector().push_back(
637 new OpGetHOTangentsOnEdge("MESH_NODE_POSITIONS"));
638 edge_fe.getOpPtrVector().push_back(
640 edge_fe.getOpPtrVector().push_back(
641 new OpEdgesFluxes(m_field, th1, th2, my_split));
642
643 for (Range::iterator fit = faces.begin(); fit != faces.end(); fit++) {
644 CHKERR moab.tag_set_data(th1, &*fit, 1, &def_val);
645 CHKERR moab.tag_set_data(th2, &*fit, 1, &def_val);
646 }
647
648 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TET_FE", tet_fe);
649 my_split << "internal\n";
650 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TRI_FE", tri_fe);
651 my_split << "skin\n";
652 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "SKIN_FE", skin_fe);
653 my_split << "edges\n";
654 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE_FE", edge_fe);
655
656 EntityHandle meshset;
657 CHKERR moab.create_meshset(MESHSET_SET, meshset);
658 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
659 BitRefLevel().set(0), BitRefLevel().set(), MBTRI, meshset);
660 CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
661 }
663
665}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static const double eps
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
static const double G_TRI_Y1[]
Definition: fem_tools.h:313
static const double G_TRI_W1[]
Definition: fem_tools.h:314
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
static const double G_TRI_X1[]
Definition: fem_tools.h:312
tee_device< std::ostream, std::ofstream > TeeDevice
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
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.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
static char help[]
static const double face_coords[4][9]
static const double edge_coords[6][6]
const double v
phase velocity of light in medium (cm/ns)
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
CoreTmp< 0 > Core
Definition: Core.hpp:1086
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
constexpr double t
plate stiffness
Definition: plate.cpp:59
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
const MatrixAdaptor getVectorN(const FieldApproximationBase base, const int gg)
get Hdiv of base functions at Gauss pts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MatrixDouble gaussPts
Matrix of integration points.
Calculate HO coordinates at gauss points.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Calculate normals at Gauss points of triangle element.
Calculate tangent vector on edge form HO geometry approximation.
transform Hcurl base fluxes from reference element to physical edge
transform Hcurl base fluxes from reference element to physical triangle
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)

Variable Documentation

◆ edge_coords

const double edge_coords[6][6]
static
Initial value:
= {
{0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 1}}

Definition at line 18 of file hcurl_continuity_check.cpp.

◆ face_coords

const double face_coords[4][9]
static
Initial value:
= {{0, 0, 0, 1, 0, 0, 0, 0, 1},
{1, 0, 0, 0, 1, 0, 0, 0, 1},
{0, 0, 0, 0, 1, 0, 0, 0, 1},
{0, 0, 0, 1, 0, 0, 0, 1, 0}}

Definition at line 13 of file hcurl_continuity_check.cpp.

◆ help

char help[] = "...\n\n"
static

Definition at line 11 of file hcurl_continuity_check.cpp.