v0.9.1
prism_elements_from_surface.cpp
Go to the documentation of this file.
1 /** \file prism_elements_from_surface.cpp
2  \example prism_elements_from_surface.cpp
3  \brief Adding prims on the surface and checking conformity between quads,
4  triangles and prism
5 
6 */
7 
8 /* This file is part of MoFEM.
9  * MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13  *
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21 
22 #include <MoFEM.hpp>
23 
24 using namespace MoFEM;
25 
26 static char help[] = "...\n\n";
27 
28 static constexpr int precision_exponent = 5;
29 static constexpr int number_of_prisms_layers = 18;
30 static constexpr double delta =
31  1. / static_cast<double>(number_of_prisms_layers);
32 static constexpr std::array<double, 3> d3 = {0, 0, 0};
33 static constexpr std::array<double, 3> d4 = {0, 0, delta};
34 
36 
37  inline static double getArg(double x) {
38  return std::round(x * pow(10., precision_exponent - 1));
39  };
40 
41  int x, y, z;
43  CoordsAndHandle(const double *coords, EntityHandle v)
44  : x(getArg(coords[0])), y(getArg(coords[1])), z(getArg(coords[2])),
45  node(v) {}
46 };
47 
48 typedef multi_index_container<
50  indexed_by<
51 
52  hashed_unique<composite_key<
53  CoordsAndHandle, member<CoordsAndHandle, int, &CoordsAndHandle::x>,
54  member<CoordsAndHandle, int, &CoordsAndHandle::y>,
55  member<CoordsAndHandle, int, &CoordsAndHandle::z>>>
56 
57  >>
59 
61 
62  PrismOp(moab::Interface &post_proc, MapCoords &map_coords);
63  MoFEMErrorCode doWork(int side, EntityType type,
65 
66 public:
69  std::vector<EntityHandle> nodeHandles;
70 };
71 
73 
74  PrismFE(MoFEM::Interface &m_field, MatrixDouble &tri_coords);
75 
76  int getRuleTrianglesOnly(int order);
77  int getRuleThroughThickness(int order);
78 
79  MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only);
80 
81  MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness);
82 
83 private:
86 };
87 
88 template <typename OP> struct Op : public OP {
89 
90  Op(moab::Interface &post_proc, MapCoords &map_coords, EntityHandle prism);
91  MoFEMErrorCode doWork(int side, EntityType type,
93 
94 public:
98  std::vector<EntityHandle> nodeHandles;
99 };
100 
102 
103  TriFE(MoFEM::Interface &m_field, MatrixDouble &tri_coords,
104  EntityHandle prims);
105 
106  int getRule(int order_row, int order_col, int order_data);
107 
108  MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data);
109 
110 private:
113 };
114 
116 
117  QuadFE(MoFEM::Interface &m_field, std::array<Range, 3> &edges_blocks,
118  EntityHandle prims);
119 
120  int getRule(int order_row, int order_col, int order_data);
121 
122  MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data);
123 
124 private:
125  std::array<Range, 3> &edgeBlocks;
127 };
128 
130 
131  EdgeFE(MoFEM::Interface &m_field, std::array<Range, 3> &edges_blocks,
132  EntityHandle prims);
133 
134  int getRule(int order_row, int order_col, int order_data);
135 
136  MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data);
137 
138 private:
139  std::array<Range, 3> &edgeBlocks;
141 };
142 
143 int main(int argc, char *argv[]) {
144 
145  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
146 
147  try {
148 
149  moab::Core mb_instance;
150  moab::Interface &moab = mb_instance;
151  int rank;
152  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
153 
154  // Read parameters from line command
155  PetscBool flg = PETSC_TRUE;
156  char mesh_file_name[255];
157  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
158  255, &flg);
159  if (flg != PETSC_TRUE)
160  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
161  "error -my_file (MESH FILE NEEDED)");
162 
163  // Read mesh to MOAB
164  const char *option;
165  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
166  CHKERR moab.load_file(mesh_file_name, 0, option);
167  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
168  if (pcomm == NULL)
169  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
170 
171  Range tris;
172  CHKERR moab.get_entities_by_type(0, MBTRI, tris, false);
173  Range tris_verts;
174  CHKERR moab.get_connectivity(tris, tris_verts);
175  MatrixDouble tri_coords(tris_verts.size(), 3);
176  CHKERR moab.get_coords(tris_verts, &tri_coords(0, 0));
177 
178  MoFEM::Core core(moab);
179  MoFEM::Interface &m_field = core;
180 
181  std::array<Range, 3> edge_block;
182  for (auto b : {1, 2, 3})
184  b, BLOCKSET, 1, edge_block[b - 1]);
185 
186  PrismsFromSurfaceInterface *prisms_from_surface_interface;
187  CHKERR m_field.getInterface(prisms_from_surface_interface);
188 
189  Range prisms;
190  CHKERR prisms_from_surface_interface->createPrisms(tris, prisms);
191  prisms_from_surface_interface->setThickness(prisms, d3.data(), d4.data());
192  Range add_prims_layer;
193  Range extrude_prisms = prisms;
194 
195  for (int ll = 1; ll != number_of_prisms_layers; ++ll) {
196  prisms_from_surface_interface->createdVertices.clear();
197  CHKERR prisms_from_surface_interface->createPrismsFromPrisms(
198  extrude_prisms, false, add_prims_layer);
199  prisms_from_surface_interface->setThickness(add_prims_layer, d3.data(),
200  d4.data());
201  extrude_prisms = add_prims_layer;
202  prisms.merge(add_prims_layer);
203  add_prims_layer.clear();
204  }
205 
206  EntityHandle meshset;
207  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER, meshset);
208  CHKERR moab.add_entities(meshset, prisms);
209 
210  MapCoords map_coords;
211  Range verts;
212  CHKERR moab.get_connectivity(prisms, verts);
213  MatrixDouble coords(verts.size(), 3);
214  CHKERR moab.get_coords(verts, &coords(0, 0));
215 
216  for (size_t v = 0; v != verts.size(); ++v)
217  map_coords.insert(CoordsAndHandle(&coords(v, 0), verts[v]));
218 
219  EntityHandle one_prism_meshset;
220  CHKERR moab.create_meshset(MESHSET_SET | MESHSET_TRACK_OWNER,
221  one_prism_meshset);
222 
223  std::array<double, 18> one_prism_coords = {0, 0, 0, 1, 0, 0, 0, 1, 0,
224  0, 0, 1, 1, 0, 1, 0, 1, 1};
225  std::array<EntityHandle, 6> one_prism_nodes;
226  for (int n = 0; n != 6; ++n)
227  CHKERR moab.create_vertex(&one_prism_coords[3 * n], one_prism_nodes[n]);
228  EntityHandle one_prism;
229  CHKERR m_field.get_moab().create_element(MBPRISM, one_prism_nodes.data(), 6,
230  one_prism);
231  Range one_prism_range;
232  one_prism_range.insert(one_prism);
233  CHKERR moab.add_entities(one_prism_meshset, one_prism_range);
234  Range one_prism_verts;
235  CHKERR moab.get_connectivity(one_prism_range, one_prism_verts);
236  CHKERR moab.add_entities(one_prism_meshset, one_prism_verts);
237  Range one_prism_adj_ents;
238  for (int d = 1; d != 3; ++d)
239  CHKERR moab.get_adjacencies(one_prism_range, d, true, one_prism_adj_ents,
240  moab::Interface::UNION);
241  CHKERR moab.add_entities(one_prism_meshset, one_prism_adj_ents);
242 
243  BitRefLevel bit_level0;
244  bit_level0.set(0);
245  CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
246  one_prism_range, bit_level0);
247  CHKERR prisms_from_surface_interface->seedPrismsEntities(one_prism_range,
248  bit_level0);
249 
250  // Fields
251  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
252  CHKERR m_field.add_ents_to_field_by_type(one_prism_meshset, MBPRISM,
253  "FIELD1", VERY_NOISY);
254 
255  CHKERR m_field.set_field_order(one_prism_meshset, MBVERTEX, "FIELD1", 1);
256  CHKERR m_field.set_field_order(one_prism_meshset, MBEDGE, "FIELD1", 5,
257  VERY_NOISY);
258  CHKERR m_field.set_field_order(one_prism_meshset, MBTRI, "FIELD1", 5);
259  CHKERR m_field.set_field_order(one_prism_meshset, MBQUAD, "FIELD1", 5,
260  VERY_NOISY);
261  CHKERR m_field.set_field_order(one_prism_meshset, MBPRISM, "FIELD1", 7,
262  VERY_NOISY);
263  CHKERR m_field.build_fields(VERY_NOISY);
264 
265  // FE
266  CHKERR m_field.add_finite_element("PRISM");
267  CHKERR m_field.add_finite_element("EDGE");
268  CHKERR m_field.add_finite_element("TRI");
269  CHKERR m_field.add_finite_element("QUAD");
270 
271  // Define rows/cols and element data
272  CHKERR m_field.modify_finite_element_add_field_row("PRISM", "FIELD1");
273  CHKERR m_field.modify_finite_element_add_field_col("PRISM", "FIELD1");
274  CHKERR m_field.modify_finite_element_add_field_data("PRISM", "FIELD1");
275  // Define rows/cols and element data
276  CHKERR m_field.modify_finite_element_add_field_row("EDGE", "FIELD1");
277  CHKERR m_field.modify_finite_element_add_field_col("EDGE", "FIELD1");
278  CHKERR m_field.modify_finite_element_add_field_data("EDGE", "FIELD1");
279  // Define rows/cols and element data
280  CHKERR m_field.modify_finite_element_add_field_row("TRI", "FIELD1");
281  CHKERR m_field.modify_finite_element_add_field_col("TRI", "FIELD1");
282  CHKERR m_field.modify_finite_element_add_field_data("TRI", "FIELD1");
283  // Define rows/cols and element data
284  CHKERR m_field.modify_finite_element_add_field_row("QUAD", "FIELD1");
285  CHKERR m_field.modify_finite_element_add_field_col("QUAD", "FIELD1");
286  CHKERR m_field.modify_finite_element_add_field_data("QUAD", "FIELD1");
287 
288  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBEDGE,
289  "EDGE");
290  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBTRI,
291  "TRI");
292  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset, MBQUAD,
293  "QUAD");
294  CHKERR m_field.add_ents_to_finite_element_by_type(one_prism_meshset,
295  MBPRISM, "PRISM");
296 
297  // build finite elemnts
298  CHKERR m_field.build_finite_elements();
299  // //build adjacencies
300  CHKERR m_field.build_adjacencies(bit_level0);
301 
302  // Problem
303  CHKERR m_field.add_problem("TEST_PROBLEM");
304 
305  // set finite elements for problem
306  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "PRISM");
307  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE");
308  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TRI");
309  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "QUAD");
310  // set refinement level for problem
311  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
312 
313  // build problem
314  ProblemsManager *prb_mng_ptr;
315  CHKERR m_field.getInterface(prb_mng_ptr);
316  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
317  // partition
318  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
319  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
320  // what are ghost nodes, see Petsc Manual
321  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
322 
323  PrismFE fe_prism(m_field, tri_coords);
324  fe_prism.getOpPtrVector().push_back(new PrismOp(moab, map_coords));
325  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "PRISM", fe_prism);
326 
327  EdgeFE fe_edge(m_field, edge_block, one_prism);
328  fe_edge.getOpPtrVector().push_back(
330  moab, map_coords, one_prism));
331  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE", fe_edge);
332 
333  TriFE fe_tri(m_field, tri_coords, one_prism);
334  fe_tri.getOpPtrVector().push_back(
336  moab, map_coords, one_prism));
337  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TRI", fe_tri);
338 
339  QuadFE fe_quad(m_field, edge_block, one_prism);
340  fe_quad.getOpPtrVector().push_back(
342  moab, map_coords, one_prism));
343  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "QUAD", fe_quad);
344 
345  CHKERR moab.write_file("prism_mesh.vtk", "VTK", "", &meshset, 1);
346  CHKERR moab.write_file("one_prism_mesh.vtk", "VTK", "", &one_prism_meshset,
347  1);
348  }
349  CATCH_ERRORS;
350 
352  return 0;
353 }
354 
357  "FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
358  postProc(post_proc), mapCoords(map_coords) {}
359 
360 MoFEMErrorCode PrismOp::doWork(int side, EntityType type,
362  constexpr double def_val[] = {0, 0, 0};
364  switch (type) {
365  case MBVERTEX:
366  case MBEDGE:
367  case MBTRI:
368  case MBQUAD:
369  case MBPRISM:
370  break;
371  default:
373  }
374  if (type == MBTRI && (side != 3 && side != 4))
376  if (type == MBQUAD && (side == 3 || side == 4))
378 
379  const int nb_dofs = data.getIndices().size();
380  for (int dd = 0; dd != nb_dofs; ++dd)
381  data.getFieldDofs()[dd]->getFieldData() = data.getIndices()[dd];
382 
383  if (type == MBVERTEX) {
384  const size_t nb_gauss_pts = getGaussPts().size2();
385  auto &coords_at_pts = getGaussPts();
386  nodeHandles.reserve(nb_gauss_pts);
387  nodeHandles.clear();
388  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
389  auto t = boost::make_tuple(CoordsAndHandle::getArg(coords_at_pts(0, gg)),
390  CoordsAndHandle::getArg(coords_at_pts(1, gg)),
391  CoordsAndHandle::getArg(coords_at_pts(2, gg)));
392 
393  auto it = mapCoords.find(t);
394  if (it != mapCoords.end())
395  nodeHandles.emplace_back(it->node);
396  else
397  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vertex not found");
398  }
399 
400  MatrixDouble node_coords(nb_gauss_pts, 3);
401  CHKERR postProc.get_coords(&nodeHandles[0], nodeHandles.size(),
402  &node_coords(0, 0));
403  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
404  constexpr double eps = 1e-12;
405  for (auto d : {0, 1, 2})
406  if (std::abs(node_coords(gg, d) - getGaussPts()(d, gg)) > eps) {
407  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
408  "Inconsistency between node coords and integration "
409  "points");
410  }
411  for (auto d : {0, 1, 2})
412  if (std::abs(node_coords(gg, d) - getCoordsAtGaussPts()(gg, d)) > eps)
413  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
414  "Inconsistency between node coords and integration "
415  "points");
416  }
417 
418  Tag th;
419  CHKERR postProc.tag_get_handle("Coords", 3, MB_TYPE_DOUBLE, th,
420  MB_TAG_CREAT | MB_TAG_DENSE, def_val);
421  CHKERR postProc.tag_set_data(th, &nodeHandles[0], nodeHandles.size(),
422  &getCoordsAtGaussPts()(0, 0));
423  }
424 
425  auto to_str = [](auto i) { return boost::lexical_cast<std::string>(i); };
426  std::string tag_name_base =
427  "PrismType" + to_str(type) + "Side" + to_str(side);
428  std::cout << tag_name_base << endl;
429 
430  MatrixDouble trans_base = trans(data.getN());
431  if (trans_base.size2() != nodeHandles.size())
432  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "wrong size %d != %d",
433  trans_base.size2(), nodeHandles.size());
434  for (size_t rr = 0; rr != trans_base.size1(); ++rr) {
435  auto tag_name = tag_name_base + "Base" + to_str(rr);
436  Tag th;
437  CHKERR postProc.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, th,
438  MB_TAG_CREAT | MB_TAG_DENSE, def_val);
439  CHKERR postProc.tag_set_data(th, &nodeHandles[0], nodeHandles.size(),
440  &trans_base(rr, 0));
441  }
442 
444 }
445 
447  : FatPrismElementForcesAndSourcesCore(m_field), triCoords(tri_coords) {}
448 
449 int PrismFE::getRuleTrianglesOnly(int order) { return -1; };
450 int PrismFE::getRuleThroughThickness(int order) { return -1; };
451 
454  gaussPtsTrianglesOnly.resize(3, triCoords.size1(), false);
455  gaussPtsTrianglesOnly.clear();
456  for (int gg = 0; gg != triCoords.size1(); ++gg)
457  for (int dd = 0; dd != 2; ++dd)
458  gaussPtsTrianglesOnly(dd, gg) = triCoords(gg, dd);
460 }
461 
465  gaussPtsThroughThickness.clear();
466  for (int gg = 0; gg != number_of_prisms_layers + 1; ++gg)
467  gaussPtsThroughThickness(0, gg) = delta * gg;
468 
470 }
471 
473  EntityHandle prism)
474  : FaceElementForcesAndSourcesCore(m_field), triCoords(tri_coords),
475  prism(prism) {}
476 
477 int TriFE::getRule(int order_row, int order_col, int order_data) { return -1; }
478 
479 MoFEMErrorCode TriFE::setGaussPts(int order_row, int order_col,
480  int order_data) {
482 
483  const EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
484  int side, sense, offset;
485  CHKERR mField.get_moab().side_number(prism, ent, side, sense, offset);
486  std::array<int, 2> swap = {0, 1};
487  if (side == 3)
488  swap = std::array<int, 2>{1, 0};
489 
490  gaussPts.resize(3, triCoords.size1(), false);
491  gaussPts.clear();
492  for (int gg = 0; gg != triCoords.size1(); ++gg)
493  for (int dd = 0; dd != 2; ++dd)
494  gaussPts(dd, gg) = triCoords(gg, swap[dd]);
495 
497 }
498 
499 QuadFE::QuadFE(MoFEM::Interface &m_field, std::array<Range, 3> &edge_blocks,
500  EntityHandle prism)
501  : FaceElementForcesAndSourcesCore(m_field), edgeBlocks(edge_blocks),
502  prism(prism) {}
503 
504 int QuadFE::getRule(int order_row, int order_col, int order_data) { return -1; }
505 
506 MoFEMErrorCode QuadFE::setGaussPts(int order_row, int order_col,
507  int order_data) {
509 
510  const EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
511  int side, sense, offset;
512  CHKERR mField.get_moab().side_number(prism, ent, side, sense, offset);
513 
514  Range edge_verts;
515  CHKERR mField.get_moab().get_connectivity(edgeBlocks[side], edge_verts);
516  MatrixDouble edge_coords(edge_verts.size(), 3);
517  CHKERR mField.get_moab().get_coords(edge_verts, &edge_coords(0, 0));
518 
519  constexpr double normal[3][2] = {{-1, 0}, {-0.5, 0.5}, {0, 1}};
520  constexpr double origin[3][2] = {{1, 0}, {1, 0}, {0, 0}};
521  constexpr int swap[3][2] = {{0, 1}, {0, 1}, {1, 0}};
522  gaussPts.resize(3, edge_verts.size() * (number_of_prisms_layers + 1), false);
523  int gg = 0;
524  for (size_t rr = 0; rr != edge_verts.size(); ++rr) {
525  const double x = edge_coords(rr, 0) - origin[side][0];
526  const double y = edge_coords(rr, 1) - origin[side][1];
527  const double edge_dist = x * normal[side][0] + y * normal[side][1];
528  for (size_t cc = 0; cc != number_of_prisms_layers + 1; ++cc, ++gg) {
529  gaussPts(swap[side][0], gg) = edge_dist;
530  gaussPts(swap[side][1], gg) = delta * cc;
531  }
532  }
533 
535 }
536 
537 EdgeFE::EdgeFE(MoFEM::Interface &m_field, std::array<Range, 3> &edge_blocks,
538  EntityHandle prism)
539  : EdgeElementForcesAndSourcesCore(m_field), edgeBlocks(edge_blocks),
540  prism(prism) {}
541 
542 int EdgeFE::getRule(int order_row, int order_col, int order_data) { return -1; }
543 
544 MoFEMErrorCode EdgeFE::setGaussPts(int order_row, int order_col,
545  int order_data) {
547 
548  const EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
549  int side, sense, offset;
550  CHKERR mField.get_moab().side_number(prism, ent, side, sense, offset);
551 
552  if (side >= 3 && side <= 5) {
553 
554  gaussPts.resize(2, number_of_prisms_layers + 1, false);
555  for (size_t gg = 0; gg != number_of_prisms_layers + 1; ++gg)
556  gaussPts(0, gg) = delta * gg;
557 
558  // cerr << gaussPts << endl;
559 
560  } else {
561 
562  constexpr double normal[3][2] = {{1, 0}, {-0.5, 0.5}, {0, 1}};
563  constexpr double origin[3][2] = {{0, 1}, {1, 0}, {0, 0}};
564  constexpr int side_map[9] = {0, 1, 2, -1, -1, -1, 0, 1, 2};
565 
566  int num_nodes;
567  const EntityHandle *conn;
568  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
569  MatrixDouble coords(num_nodes, 3);
570  CHKERR mField.get_moab().get_coords(conn, num_nodes,
571  &*coords.data().begin());
572  side = side_map[side];
573 
574  Range edge_verts;
575  CHKERR mField.get_moab().get_connectivity(edgeBlocks[side], edge_verts);
576  MatrixDouble edge_coords(edge_verts.size(), 3);
577  CHKERR mField.get_moab().get_coords(edge_verts, &edge_coords(0, 0));
578  gaussPts.resize(2, edge_verts.size(), false);
579 
580  for (size_t gg = 0; gg != edge_verts.size(); ++gg) {
581  const double x = edge_coords(gg, 0) - origin[side][0];
582  const double y = edge_coords(gg, 1) - origin[side][1];
583  const double edge_dist = x * normal[side][0] + y * normal[side][1];
584  gaussPts(0, gg) = edge_dist;
585  }
586  }
587 
589 }
590 
591 template <typename OP>
592 Op<OP>::Op(moab::Interface &post_proc, MapCoords &map_coords,
593  EntityHandle prism)
594  : OP("FIELD1", "FIELD1", ForcesAndSourcesCore::UserDataOperator::OPROW),
595  postProc(post_proc), mapCoords(map_coords), prism(prism) {}
596 
597 template <typename OP>
598 MoFEMErrorCode Op<OP>::doWork(int side, EntityType type,
600  constexpr double def_val[] = {0, 0, 0};
602  switch (type) {
603  case MBVERTEX:
604  case MBEDGE:
605  case MBTRI:
606  case MBQUAD:
607  break;
608  default:
610  }
611 
612  const int nb_dofs = data.getIndices().size();
613  for (int dd = 0; dd != nb_dofs; ++dd)
614  if (data.getFieldData()[dd] != data.getIndices()[dd]) {
615  std::cerr << "Indices: " << data.getIndices() << std::endl;
616  std::cerr << "Local indices: " << data.getLocalIndices() << std::endl;
617  std::cerr << "Data: " << data.getFieldData() << std::endl;
618  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
619  "Indicices/data inconsistency %3.1f != %d",
620  data.getFieldData()[dd], data.getIndices()[dd]);
621  }
622 
623  const EntityHandle fe_ent = OP::getFEEntityHandle();
624  const EntityHandle ent = OP::getSideEntity(side, type);
625  int side_prism, sense, offset;
626  if (type == MBVERTEX) {
627  CHKERR postProc.side_number(prism, fe_ent, side_prism, sense, offset);
628  } else
629  CHKERR postProc.side_number(prism, ent, side_prism, sense, offset);
630 
631  if (type == MBVERTEX) {
632  auto &coords_at_pts = OP::getCoordsAtGaussPts();
633  const size_t nb_gauss_pts = coords_at_pts.size1();
634 
635  nodeHandles.reserve(nb_gauss_pts);
636  nodeHandles.clear();
637  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
638  auto t = boost::make_tuple(CoordsAndHandle::getArg(coords_at_pts(gg, 0)),
639  CoordsAndHandle::getArg(coords_at_pts(gg, 1)),
640  CoordsAndHandle::getArg(coords_at_pts(gg, 2)));
641 
642  auto it = mapCoords.find(t);
643  if (it != mapCoords.end())
644  nodeHandles.emplace_back(it->node);
645  else
646  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vertex not found");
647  }
648  }
649 
650  auto to_str = [](auto i) { return boost::lexical_cast<std::string>(i); };
651  std::string tag_name_base =
652  "FEType" + to_str(OP::getNumeredEntFiniteElementPtr()->getEntType()) +
653  "Type" + to_str(type) + "Side" + to_str(side_prism);
654 
655  std::string tag_prism_name_base =
656  "PrismType" + to_str(type) + "Side" + to_str(side_prism);
657 
658  MatrixDouble trans_base = trans(data.getN());
659  MatrixDouble prism_base(trans_base.size1(), trans_base.size2());
660  if (trans_base.size2() != nodeHandles.size())
661  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "wrong size %d != %d",
662  trans_base.size2(), nodeHandles.size());
663 
664  for (size_t rr = 0; rr != trans_base.size1(); ++rr) {
665 
666  std::string tag_name = tag_name_base + "Base";
667  if (type == MBVERTEX) {
668  EntityHandle node = data.getFieldDofs()[rr]->getEnt();
669  int prism_mode_side;
670  CHKERR postProc.side_number(prism, node, prism_mode_side, sense, offset);
671  tag_name += to_str(prism_mode_side);
672  } else {
673  tag_name += to_str(rr);
674  }
675 
676  std::cout << tag_name << endl;
677 
678  Tag th;
679  CHKERR postProc.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, th,
680  MB_TAG_CREAT | MB_TAG_DENSE, def_val);
681  CHKERR postProc.tag_set_data(th, &nodeHandles[0], nodeHandles.size(),
682  &trans_base(rr, 0));
683 
684  if (type != MBVERTEX) {
685  auto tag_prism_name = tag_prism_name_base + "Base" + to_str(rr);
686  Tag th_prism;
687  CHKERR postProc.tag_get_handle(tag_prism_name.c_str(), th_prism);
688  CHKERR postProc.tag_get_data(th_prism, &nodeHandles[0],
689  nodeHandles.size(), &prism_base(rr, 0));
690  }
691  }
692 
693  auto sum_matrix = [](MatrixDouble &m) {
694  double s = 0;
695  for (unsigned int ii = 0; ii < m.size1(); ii++) {
696  for (unsigned int jj = 0; jj < m.size2(); jj++) {
697  s += std::abs(m(ii, jj));
698  }
699  }
700  return s;
701  };
702 
703  if (type != MBVERTEX) {
704  prism_base -= trans_base;
705  double sum = sum_matrix(prism_base);
706  constexpr double eps = 1e-6;
707 
708  if (std::abs(sum) > eps)
709  cout << "Inconsistent base " << tag_prism_name_base << " "
710  << tag_name_base << " sum " << sum << endl;
711 
712  if (std::abs(sum) > eps)
713  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
714  "Inconsistent base %s sum %6.4e", tag_prism_name_base.c_str(),
715  sum);
716  }
717 
719 }
std::array< Range, 3 > & edgeBlocks
MatrixDouble & triCoords
int getRule(int order_row, int order_col, int order_data)
another variant of getRule
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
static double sum_matrix(MatrixDouble &m)
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
Deprecated interface functions.
Problem manager is used to build and partition problems.
QuadFE(MoFEM::Interface &m_field, std::array< Range, 3 > &edges_blocks, EntityHandle prims)
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MatrixDouble gaussPts
Matrix of integration points.
static double getArg(double x)
static constexpr std::array< double, 3 > d4
multi_index_container< CoordsAndHandle, indexed_by< hashed_unique< composite_key< CoordsAndHandle, member< CoordsAndHandle, int, &CoordsAndHandle::x >, member< CoordsAndHandle, int, &CoordsAndHandle::y >, member< CoordsAndHandle, int, &CoordsAndHandle::z > > > > > MapCoords
virtual moab::Interface & get_moab()=0
moab::Interface & postProc
static constexpr std::array< double, 3 > d3
MoFEMErrorCode setThickness(const Range &prisms, const double director3[], const double director4[])
MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimensionNodeset can contain nodes,...
MoFEMErrorCode seedPrismsEntities(Range &prisms, const BitRefLevel &bit, int verb=-1)
Seed prism entities by bit level.
int getRuleThroughThickness(int order)
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
static constexpr double delta
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
std::vector< EntityHandle > nodeHandles
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
MatrixDouble & triCoords
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.
int getRule(int order_row, int order_col, int order_data)
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
Interface for managing meshsets containing materials and boundary conditions.
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:513
int getRuleTrianglesOnly(int order)
PrismFE(MoFEM::Interface &m_field, MatrixDouble &tri_coords)
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]
static constexpr int number_of_prisms_layers
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 > numeredEntFiniteElementPtr
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
static constexpr int precision_exponent
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
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
constexpr int order
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
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.
MoFEMErrorCode createPrisms(const Range &ents, Range &prisms, int verb=-1)
Make prisms from triangles.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
std::map< EntityHandle, EntityHandle > createdVertices
MoFEMErrorCode createPrismsFromPrisms(const Range &prisms, bool from_down, Range &out_prisms, int verb=-1)
Make prisms by extruding top or bottom prisms.
MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
EdgeFE(MoFEM::Interface &m_field, std::array< Range, 3 > &edges_blocks, EntityHandle prims)
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]
static char help[]
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MapCoords & mapCoords
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:601
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 doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
ForcesAndSourcesCore::UserDataOperator UserDataOperator
moab::Interface & postProc
std::vector< EntityHandle > nodeHandles
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:290
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
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
int main(int argc, char *argv[])
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
Op(moab::Interface &post_proc, MapCoords &map_coords, EntityHandle prism)
TriFE(MoFEM::Interface &m_field, MatrixDouble &tri_coords, EntityHandle prims)
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
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....
PrismOp(moab::Interface &post_proc, MapCoords &map_coords)
structure to get information form mofem into DataForcesAndSourcesCore
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
FatPrism finite elementUser is implementing own operator at Gauss points level, by own object derived...
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
continuous field
Definition: definitions.h:176
CoordsAndHandle(const double *coords, EntityHandle v)
EntityHandle prism
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....
const VectorDouble & getFieldData() const
get dofs values
MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
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...
std::array< Range, 3 > & edgeBlocks
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)
int getRule(int order_row, int order_col, int order_data)
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26