v0.9.0
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 34 of file hcurl_continuity_check.cpp.

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

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 30 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 25 of file hcurl_continuity_check.cpp.

◆ help

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

Definition at line 23 of file hcurl_continuity_check.cpp.