v0.14.0
hcurl_continuity_check.cpp
Go to the documentation of this file.
1 
2 
3 #include <MoFEM.hpp>
4 
5 namespace bio = boost::iostreams;
6 using bio::stream;
7 using bio::tee_device;
8 
9 using namespace MoFEM;
10 
11 static char help[] = "...\n\n";
12 
13 static const double face_coords[4][9] = {{0, 0, 0, 1, 0, 0, 0, 0, 1},
14  {1, 0, 0, 0, 1, 0, 0, 0, 1},
15  {0, 0, 0, 0, 1, 0, 0, 0, 1},
16  {0, 0, 0, 1, 0, 0, 0, 1, 0}};
17 
18 static const double edge_coords[6][6] = {
19  {0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 0},
20  {0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 1}};
21 
22 int main(int argc, char *argv[]) {
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_IMPOSSIBLE_CASE, "base not set");
62  }
63 
64  FieldApproximationBase base = AINSWORTH_LEGENDRE_BASE;
65  if (choice_base_value == AINSWORTH) {
66  base = AINSWORTH_LEGENDRE_BASE;
67  } else if (choice_base_value == DEMKOWICZ) {
68  base = DEMKOWICZ_JACOBI_BASE;
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");
100  CHKERR m_field.modify_finite_element_add_field_data("TET_FE",
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");
106  CHKERR m_field.modify_finite_element_add_field_data("SKIN_FE",
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");
112  CHKERR m_field.modify_finite_element_add_field_data("TRI_FE",
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");
118  CHKERR m_field.modify_finite_element_add_field_data("EDGE_FE",
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
189  CHKERR m_field.build_finite_elements();
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
229  : public VolumeElementForcesAndSourcesCore::UserDataOperator {
230 
231  MoFEM::Interface &mField;
232  Tag tH;
233 
234  OpTetFluxes(MoFEM::Interface &m_field, Tag th)
235  : VolumeElementForcesAndSourcesCore::UserDataOperator(
236  "HCURL", UserDataOperator::OPROW),
237  mField(m_field), tH(th) {}
238 
239  MoFEMErrorCode doWork(int side, EntityType type,
240  EntitiesFieldData::EntData &data) {
241  MoFEMFunctionBegin;
242 
243  if (data.getFieldData().size() == 0)
244  MoFEMFunctionReturnHot(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 
312  MoFEMFunctionReturn(0);
313  }
314  };
315 
316  struct MyTetFE : public VolumeElementForcesAndSourcesCore {
317 
318  MyTetFE(MoFEM::Interface &m_field)
319  : VolumeElementForcesAndSourcesCore(m_field) {}
320  int getRule(int order) { return -1; };
321 
322  MatrixDouble N_tri;
323  MoFEMErrorCode setGaussPts(int order) {
324 
325  MoFEMFunctionBegin;
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 
362  MoFEMFunctionReturn(0);
363  }
364  };
365 
366  struct OpFacesSkinFluxes
367  : public FaceElementForcesAndSourcesCore::UserDataOperator {
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)
375  : FaceElementForcesAndSourcesCore::UserDataOperator(
376  "HCURL", UserDataOperator::OPROW),
377  mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
378 
379  MoFEMErrorCode doWork(int side, EntityType type,
380  EntitiesFieldData::EntData &data) {
381  MoFEMFunctionBegin;
382 
383  if (type != MBTRI)
384  MoFEMFunctionReturnHot(0);
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 
426  MoFEMFunctionReturn(0);
427  }
428  };
429 
430  struct OpFacesFluxes
431  : public FaceElementForcesAndSourcesCore::UserDataOperator {
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)
439  : FaceElementForcesAndSourcesCore::UserDataOperator(
440  "HCURL", UserDataOperator::OPROW),
441  mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
442 
443  MoFEMErrorCode doWork(int side, EntityType type,
444  EntitiesFieldData::EntData &data) {
445  MoFEMFunctionBegin;
446 
447  if (type != MBTRI)
448  MoFEMFunctionReturnHot(0);
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 
477  MoFEMFunctionReturn(0);
478  }
479  };
480 
481  struct MyTriFE : public FaceElementForcesAndSourcesCore {
482 
483  MyTriFE(MoFEM::Interface &m_field)
484  : FaceElementForcesAndSourcesCore(m_field) {}
485  int getRule(int order) { return -1; };
486 
487  MoFEMErrorCode setGaussPts(int order) {
488  MoFEMFunctionBeginHot;
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 
495  MoFEMFunctionReturnHot(0);
496  }
497  };
498 
499  struct OpEdgesFluxes
500  : public EdgeElementForcesAndSourcesCore::UserDataOperator {
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)
508  : EdgeElementForcesAndSourcesCore::UserDataOperator(
509  "HCURL", UserDataOperator::OPROW),
510  mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
511 
512  MoFEMErrorCode doWork(int side, EntityType type,
513  EntitiesFieldData::EntData &data) {
514  MoFEMFunctionBegin;
515 
516  if (type != MBEDGE)
517  MoFEMFunctionReturnHot(0);
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 
569  MoFEMFunctionReturn(0);
570  }
571  };
572 
573  struct MyEdgeFE : public EdgeElementForcesAndSourcesCore {
574 
575  MyEdgeFE(MoFEM::Interface &m_field)
576  : EdgeElementForcesAndSourcesCore(m_field) {}
577  int getRule(int order) { return -1; };
578 
579  MoFEMErrorCode setGaussPts(int order) {
580  MoFEMFunctionBeginHot;
581 
582  gaussPts.resize(2, 1);
583  gaussPts(0, 0) = 0.5;
584  gaussPts(1, 0) = 1;
585 
586  MoFEMFunctionReturnHot(0);
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(
625  new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
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(
632  new OpHOSetCovariantPiolaTransformOnFace3D(HCURL));
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(
639  new OpHOSetContravariantPiolaTransformOnEdge3D(HCURL));
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  }
662  CATCH_ERRORS;
663 
664  CHKERR MoFEM::Core::Finalize();
665 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
TeeStream
stream< TeeDevice > TeeStream
Definition: forces_and_sources_testing_contact_prism_element.cpp:10
G_TRI_X1
static const double G_TRI_X1[]
Definition: fem_tools.h:312
MoFEM::CoreInterface::loop_finite_elements
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.
MoFEM::OpGetHOTangentsOnEdge
Calculate tangent vector on edge form HO geometry approximation.
Definition: HODataOperators.hpp:369
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
EntityHandle
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
OpFacesFluxes
Definition: hdiv_divergence_operator.cpp:31
MoFEM.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MoFEM::OpCalculateHOCoords
Calculate HO coordinates at gauss points.
Definition: HODataOperators.hpp:41
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
TeeDevice
tee_device< std::ostream, std::ofstream > TeeDevice
Definition: forces_and_sources_testing_contact_prism_element.cpp:9
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::OpSetHOCovariantPiolaTransform
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
Definition: HODataOperators.hpp:206
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::CoreInterface::add_field
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.
convert.type
type
Definition: convert.py:64
SideNumber_multiIndex
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.
Definition: RefEntsMultiIndices.hpp:101
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
G_TRI_Y1
static const double G_TRI_Y1[]
Definition: fem_tools.h:313
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpGetHONormalsOnFace
Calculate normals at Gauss points of triangle element.
Definition: HODataOperators.hpp:281
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
MoFEM::ProblemsManager::partitionFiniteElements
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
Definition: ProblemsManager.cpp:2167
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
MoFEM::EntitiesFieldData::EntData::getVectorN
const MatrixAdaptor getVectorN(const FieldApproximationBase base, const int gg)
get Hdiv of base functions at Gauss pts
Definition: EntitiesFieldData.hpp:1434
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
MoFEM::CoreTmp< 0 >::Initialize
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
FTensor::dd
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
MoFEM::OpHOSetCovariantPiolaTransformOnFace3D
transform Hcurl base fluxes from reference element to physical triangle
Definition: HODataOperators.hpp:336
face_coords
static const double face_coords[4][9]
Definition: hcurl_continuity_check.cpp:13
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
G_TRI_W1
static const double G_TRI_W1[]
Definition: fem_tools.h:314
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
MoFEM::ProblemsManager::partitionSimpleProblem
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
Definition: ProblemsManager.cpp:1537
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
help
static char help[]
Definition: hcurl_continuity_check.cpp:11
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:159
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
main
int main(int argc, char *argv[])
Definition: hcurl_continuity_check.cpp:22
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::OpSetHOInvJacVectorBase
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Definition: HODataOperators.hpp:91
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::CoreInterface::set_field_order
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.
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
edge_coords
static const double edge_coords[6][6]
Definition: hcurl_continuity_check.cpp:18
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::OpHOSetContravariantPiolaTransformOnEdge3D
transform Hcurl base fluxes from reference element to physical edge
Definition: HODataOperators.hpp:317
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow