v0.9.0
forces_and_sources_h1_continuity_check.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #include <MoFEM.hpp>
16 
17 namespace bio = boost::iostreams;
18 using bio::stream;
19 using bio::tee_device;
20 
21 using namespace MoFEM;
22 
23 static char help[] = "...\n\n";
24 
25 static const double face_coords[4][9] = {{0, 0, 0, 1, 0, 0, 0, 0, 1},
26  {1, 0, 0, 0, 1, 0, 0, 0, 1},
27  {0, 0, 0, 0, 1, 0, 0, 0, 1},
28  {0, 0, 0, 1, 0, 0, 0, 1, 0}};
29 
30 static const double edge_coords[6][6] = {
31  {0, 0, 0, 1, 0, 0}, {1, 0, 0, 0, 1, 0}, {0, 1, 0, 0, 0, 0},
32  {0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 1}, {0, 1, 0, 0, 0, 1}};
33 
34 int main(int argc, char *argv[]) {
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  ierr = PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
49  255, &flg);
50  CHKERRG(ierr);
51 #else
52  ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
53  mesh_file_name, 255, &flg);
54  CHKERRG(ierr);
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 
64  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
65  if (pcomm == NULL)
66  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
67 
68  // Create MoFEM (Joseph) database
69  MoFEM::Core core(moab);
70  MoFEM::Interface &m_field = core;
71 
72  // set entitities bit level
73  BitRefLevel bit_level0;
74  bit_level0.set(0);
75  EntityHandle meshset_level0;
76  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
77 
78  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
79  0, 3, bit_level0);
80 
81 
82  // Fields
83  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
84  3);
85 
86  CHKERR m_field.add_field("H1", H1, AINSWORTH_LEGENDRE_BASE, 1);
87 
88 
89  // FE
90  CHKERR m_field.add_finite_element("TET_FE");
91  CHKERR m_field.add_finite_element("TRI_FE");
92  CHKERR m_field.add_finite_element("SKIN_FE");
93  CHKERR m_field.add_finite_element("EDGE_FE");
94 
95 
96  // Define rows/cols and element data
97  CHKERR m_field.modify_finite_element_add_field_row("TET_FE", "H1");
98  CHKERR m_field.modify_finite_element_add_field_col("TET_FE", "H1");
99  CHKERR m_field.modify_finite_element_add_field_data("TET_FE", "H1");
101  "MESH_NODE_POSITIONS");
102 
103 
104  CHKERR m_field.modify_finite_element_add_field_row("SKIN_FE", "H1");
105  CHKERR m_field.modify_finite_element_add_field_col("SKIN_FE", "H1");
106  CHKERR m_field.modify_finite_element_add_field_data("SKIN_FE", "H1");
108  "MESH_NODE_POSITIONS");
109 
110 
111  CHKERR m_field.modify_finite_element_add_field_row("TRI_FE", "H1");
112  CHKERR m_field.modify_finite_element_add_field_col("TRI_FE", "H1");
113  CHKERR m_field.modify_finite_element_add_field_data("TRI_FE", "H1");
115  "MESH_NODE_POSITIONS");
116 
117 
118  CHKERR m_field.modify_finite_element_add_field_row("EDGE_FE", "H1");
119  CHKERR m_field.modify_finite_element_add_field_col("EDGE_FE", "H1");
120  CHKERR m_field.modify_finite_element_add_field_data("EDGE_FE", "H1");
122  "MESH_NODE_POSITIONS");
123 
124  // Problem
125  CHKERR m_field.add_problem("TEST_PROBLEM");
126 
127 
128  // set finite elements for problem
129  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TET_FE");
130  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "SKIN_FE");
131  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TRI_FE");
132  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE_FE");
133 
134  // set refinement level for problem
135  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
136 
137 
138  // meshset consisting all entities in mesh
139  EntityHandle root_set = moab.get_root_set();
140  // add entities to field
141  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "H1");
142 
143 
144  // add entities to finite element
145  CHKERR
146  m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "TET_FE");
147 
148 
149  Range tets;
150  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
151  BitRefLevel().set(0), BitRefLevel().set(), MBTET, tets);
152 
153  Skinner skin(&moab);
154  Range skin_faces; // skin faces from 3d ents
155  CHKERR skin.find_skin(0, tets, false, skin_faces);
156  CHKERR m_field.add_ents_to_finite_element_by_type(skin_faces, MBTRI,
157  "SKIN_FE");
158 
159  Range faces;
160  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
161  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, faces);
162 
163  faces = subtract(faces, skin_faces);
164  CHKERR m_field.add_ents_to_finite_element_by_type(faces, MBTRI, "TRI_FE");
165 
166  Range edges;
167  CHKERR moab.get_adjacencies(faces, 1, false, edges, moab::Interface::UNION);
168  CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE, "EDGE_FE");
169 
170  Range skin_edges;
171  CHKERR moab.get_adjacencies(skin_faces, 1, false, skin_edges,
172  moab::Interface::UNION);
173  CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE, "EDGE_FE");
174 
175 
176  // set app. order
177  int order = 4;
178  CHKERR m_field.set_field_order(root_set, MBTET, "H1", order);
179  CHKERR m_field.set_field_order(root_set, MBTRI, "H1", order);
180  CHKERR m_field.set_field_order(root_set, MBEDGE, "H1", order);
181 
182  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
183  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
184  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
185  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
186  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
187 
188  /****/
189  // build database
190  // build field
191  CHKERR m_field.build_fields();
192 
193  // build finite elemnts
194  CHKERR m_field.build_finite_elements();
195 
196  // build adjacencies
197  CHKERR m_field.build_adjacencies(bit_level0);
198 
199 
200  ProblemsManager *prb_mng_ptr;
201  CHKERR m_field.getInterface(prb_mng_ptr);
202 
203  // build problem
204  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
205 
206 
207  // project geometry form 10 node tets on higher order approx. functions
208  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
209  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
210 
211  /****/
212  // mesh partitioning
213  // partition
214  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
215  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
216 
217  // what are ghost nodes, see Petsc Manual
218  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
219 
220 
221  Vec v;
222  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("TEST_PROBLEM",
223  ROW, &v);
224  CHKERR VecSetRandom(v, PETSC_NULL);
225 
226  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
227  "TEST_PROBLEM", ROW, v, INSERT_VALUES, SCATTER_REVERSE);
228 
229  CHKERR VecDestroy(&v);
230 
231 
232  typedef tee_device<std::ostream, std::ofstream> TeeDevice;
233  typedef stream<TeeDevice> TeeStream;
234  std::ofstream ofs("forces_and_sources_hdiv_continuity_check.txt");
235  TeeDevice my_tee(std::cout, ofs);
236  TeeStream my_split(my_tee);
237 
238  struct OpTetFluxes
240 
241  MoFEM::Interface &mField;
242  Tag tH;
243 
244  OpTetFluxes(MoFEM::Interface &m_field, Tag th)
246  "H1", UserDataOperator::OPROW),
247  mField(m_field), tH(th) {}
248 
249  MoFEMErrorCode doWork(int side, EntityType type,
252 
253  if (data.getFieldData().size() == 0)
255 
256  if (type == MBTRI) {
257 
258  boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
259  getNumeredEntFiniteElementPtr();
260  SideNumber_multiIndex &side_table = mofem_fe->getSideNumberTable();
261  EntityHandle face = side_table.get<1>()
262  .find(boost::make_tuple(type, side))
263  ->get()
264  ->ent;
265  int sense = side_table.get<1>()
266  .find(boost::make_tuple(type, side))
267  ->get()
268  ->sense;
269 
270  double t = 0;
271  int nb_dofs = data.getN().size2();
272  for (int dd = 0; dd < nb_dofs; dd++) {
273  t += data.getN(side)(dd) * data.getFieldData()[dd];
274  }
275 
276  double *t_ptr;
277  CHKERR mField.get_moab().tag_get_by_ptr(tH, &face, 1,
278  (const void **)&t_ptr);
279 
280  *t_ptr += sense * t;
281  }
282 
283  if (type == MBEDGE) {
284 
285  boost::shared_ptr<const NumeredEntFiniteElement> mofem_fe =
286  getNumeredEntFiniteElementPtr();
287  SideNumber_multiIndex &side_table = mofem_fe->getSideNumberTable();
288  EntityHandle edge = side_table.get<1>()
289  .find(boost::make_tuple(type, side))
290  ->get()
291  ->ent;
292  Range adj_tets;
293  CHKERR mField.get_moab().get_adjacencies(&edge, 1, 3, false, adj_tets,
294  moab::Interface::UNION);
295  const int nb_adj_tets = adj_tets.size();
296 
297  double t = 0;
298  int nb_dofs = data.getN().size2();
299  for (int dd = 0; dd < nb_dofs; dd++) {
300  t += data.getN(4 + side)(dd) * data.getFieldData()[dd];
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  *t_ptr += t / nb_adj_tets;
308  }
309 
311  }
312  };
313 
314  struct MyTetFE : public VolumeElementForcesAndSourcesCore {
315 
316  MyTetFE(MoFEM::Interface &m_field)
318  int getRule(int order) { return -1; };
319 
320  MatrixDouble N_tri;
321  MoFEMErrorCode setGaussPts(int order) {
322 
324 
325  try {
326 
327  N_tri.resize(1, 3);
328  CHKERR ShapeMBTRI(&N_tri(0, 0), G_TRI_X1, G_TRI_Y1, 1);
329 
330 
331  gaussPts.resize(4, 4 + 6);
332  int ff = 0;
333  for (; ff < 4; ff++) {
334  int dd = 0;
335  for (; dd < 3; dd++) {
336  gaussPts(dd, ff) =
337  cblas_ddot(3, &N_tri(0, 0), 1, &face_coords[ff][dd], 3);
338  }
339  gaussPts(3, ff) = G_TRI_W1[0];
340  }
341 
342  int ee = 0;
343  for (; ee < 6; ee++) {
344  int dd = 0;
345  for (; dd < 3; dd++) {
346  gaussPts(dd, 4 + ee) =
347  (edge_coords[ee][0 + dd] + edge_coords[ee][3 + dd]) * 0.5;
348  }
349  gaussPts(3, 4 + ee) = 1;
350  }
351 
352  // std::cerr << gaussPts << std::endl;
353 
354  } catch (std::exception &ex) {
355  std::ostringstream ss;
356  ss << "thorw in method: " << ex.what() << " at line " << __LINE__
357  << " in file " << __FILE__;
358  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
359  }
360 
362  }
363  };
364 
365  struct OpFacesSkinFluxes
367 
368  MoFEM::Interface &mField;
369  Tag tH1, tH2;
370  TeeStream &mySplit;
371 
372  OpFacesSkinFluxes(MoFEM::Interface &m_field, Tag th1, Tag th2,
373  TeeStream &my_split)
375  "H1", UserDataOperator::OPROW),
376  mField(m_field), tH1(th1), tH2(th2), mySplit(my_split) {}
377 
378  MoFEMErrorCode doWork(int side, EntityType type,
381 
382  if (type != MBTRI)
384  EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
385 
386  double *t_ptr;
387  CHKERR mField.get_moab().tag_get_by_ptr(tH1, &face, 1,
388  (const void **)&t_ptr);
389 
390  double *tn_ptr;
391  CHKERR mField.get_moab().tag_get_by_ptr(tH2, &face, 1,
392  (const void **)&tn_ptr);
393 
394 
395  *tn_ptr = *t_ptr;
396 
397  int dd = 0;
398  int nb_dofs = data.getN().size2();
399  for (; dd < nb_dofs; dd++) {
400  double val = data.getFieldData()[dd];
401  *tn_ptr += -data.getN()(0, dd) * val;
402  }
403 
404  const double eps = 1e-8;
405  if (fabs(*tn_ptr) > eps) {
406  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
407  "H1 continuity failed %6.4e", *tn_ptr);
408  }
409 
410  mySplit.precision(5);
411  mySplit << face << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
412 
414  }
415  };
416 
417  struct OpFacesFluxes
419 
420  MoFEM::Interface &mField;
421  Tag tH1, tH2;
422  TeeStream &mySplit;
423 
424  OpFacesFluxes(MoFEM::Interface &m_field, Tag _th1, Tag _th2,
425  TeeStream &my_split)
427  "H1", UserDataOperator::OPROW),
428  mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
429 
430  MoFEMErrorCode doWork(int side, EntityType type,
433 
434  if (type != MBTRI)
436 
437  EntityHandle face = getNumeredEntFiniteElementPtr()->getEnt();
438 
439  double *t_ptr;
440  CHKERR mField.get_moab().tag_get_by_ptr(tH1, &face, 1,
441  (const void **)&t_ptr);
442 
443  double *tn_ptr;
444  CHKERR mField.get_moab().tag_get_by_ptr(tH2, &face, 1,
445  (const void **)&tn_ptr);
446 
447 
448  *tn_ptr = *t_ptr;
449 
450  const double eps = 1e-8;
451  if (fabs(*tn_ptr) > eps) {
452  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
453  "H1 continuity failed %6.4e", *tn_ptr);
454  }
455 
456  mySplit.precision(5);
457 
458  mySplit << face << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
459 
461  }
462  };
463 
464  struct MyTriFE : public FaceElementForcesAndSourcesCore {
465 
466  MyTriFE(MoFEM::Interface &m_field)
467  : FaceElementForcesAndSourcesCore(m_field) {}
468  int getRule(int order) { return -1; };
469 
470  MoFEMErrorCode setGaussPts(int order) {
472 
473  gaussPts.resize(3, 1);
474  gaussPts(0, 0) = G_TRI_X1[0];
475  gaussPts(1, 0) = G_TRI_Y1[0];
476  gaussPts(2, 0) = G_TRI_W1[0];
477 
479  }
480  };
481 
482  struct OpEdgesFluxes
484 
485  MoFEM::Interface &mField;
486  Tag tH1, tH2;
487  TeeStream &mySplit;
488 
489  OpEdgesFluxes(MoFEM::Interface &m_field, Tag _th1, Tag _th2,
490  TeeStream &my_split)
492  "H1", UserDataOperator::OPROW),
493  mField(m_field), tH1(_th1), tH2(_th2), mySplit(my_split) {}
494 
495  MoFEMErrorCode doWork(int side, EntityType type,
498 
499  if (type != MBEDGE)
501 
502  EntityHandle edge = getNumeredEntFiniteElementPtr()->getEnt();
503 
504  double *t_ptr;
505  CHKERR mField.get_moab().tag_get_by_ptr(tH1, &edge, 1,
506  (const void **)&t_ptr);
507 
508  double *tn_ptr;
509  CHKERR mField.get_moab().tag_get_by_ptr(tH2, &edge, 1,
510  (const void **)&tn_ptr);
511 
512 
513  *tn_ptr = *t_ptr;
514 
515  double tn = 0;
516  int nb_dofs = data.getN().size2();
517  int dd = 0;
518  for (; dd < nb_dofs; dd++) {
519  double val = data.getFieldData()[dd];
520  tn += data.getN()(0, dd) * val;
521  }
522 
523  // mySplit << *tn_ptr << " " << tn << " " << getLength() << endl;
524  *tn_ptr -= tn;
525 
526  // mySplit << getTangetAtGaussPts() << " " << getDirection() << endl;
527 
528  // cerr << t_ptr[0] << " " << t_ptr[1] << " " << t_ptr[2] << endl;
529 
530  // const double eps = 1e-8;
531  // if(fabs(*tn_ptr)>eps) {
532  // SETERRQ1(
533  // PETSC_COMM_SELF,
534  // MOFEM_ATOM_TEST_INVALID,
535  // "H1 continuity failed %6.4e",
536  // *tn_ptr
537  // );
538  // }
539 
540  mySplit.precision(5);
541 
542  mySplit << edge << " " << /*std::fixed <<*/ fabs(*tn_ptr) << std::endl;
543 
545  }
546  };
547 
548  struct MyEdgeFE : public EdgeElementForcesAndSourcesCore {
549 
550  MyEdgeFE(MoFEM::Interface &m_field)
551  : EdgeElementForcesAndSourcesCore(m_field) {}
552  int getRule(int order) { return -1; };
553 
554  MoFEMErrorCode setGaussPts(int order) {
556 
557  gaussPts.resize(2, 1);
558  gaussPts(0, 0) = 0.5;
559  gaussPts(1, 0) = 1;
560 
562  }
563  };
564 
565  MyTetFE tet_fe(m_field);
566  MyTriFE tri_fe(m_field);
567  MyTriFE skin_fe(m_field);
568  MyEdgeFE edge_fe(m_field);
569 
570  Tag th1;
571  double def_val[] = {0, 0, 0};
572  CHKERR moab.tag_get_handle("T", 3, MB_TYPE_DOUBLE, th1,
573  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
574 
575  tet_fe.getOpPtrVector().push_back(new OpTetFluxes(m_field, th1));
576 
577  Tag th2;
578  CHKERR moab.tag_get_handle("TN", 1, MB_TYPE_DOUBLE, th2,
579  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
580 
581  tri_fe.getOpPtrVector().push_back(
582  new OpFacesFluxes(m_field, th1, th2, my_split));
583  skin_fe.getOpPtrVector().push_back(
584  new OpFacesSkinFluxes(m_field, th1, th2, my_split));
585  edge_fe.getOpPtrVector().push_back(
586  new OpEdgesFluxes(m_field, th1, th2, my_split));
587 
588  for (Range::iterator fit = faces.begin(); fit != faces.end(); fit++) {
589  CHKERR moab.tag_set_data(th1, &*fit, 1, &def_val);
590 
591  CHKERR moab.tag_set_data(th2, &*fit, 1, &def_val);
592 
593  }
594 
595  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TET_FE", tet_fe);
596  my_split << "internal\n";
597  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TRI_FE", tri_fe);
598  my_split << "skin\n";
599  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "SKIN_FE", skin_fe);
600  my_split << "edges\n";
601  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE_FE", edge_fe);
602 
603 
604  EntityHandle meshset;
605  CHKERR moab.create_meshset(MESHSET_SET, meshset);
606  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
607  BitRefLevel().set(0), BitRefLevel().set(), MBTRI, meshset);
608 
609  CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
610 
611  }
612  CATCH_ERRORS;
613 
615 }
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
int main(int argc, char *argv[])
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
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
virtual moab::Interface & get_moab()=0
static const double G_TRI_X1[]
Definition: fem_tools.h:247
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
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
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
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 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...
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
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
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
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
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
static const double edge_coords[6][6]
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:36
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
static const double face_coords[4][9]
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
#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
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
const VectorDouble & getFieldData() const
get dofs values
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)