v0.14.0
ContactPrismElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file ContactPrismElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of contact prism element
4 
5 */
6 
7 
8 
9 namespace MoFEM {
10 
13  : ForcesAndSourcesCore(m_field),
14  dataOnMaster{
15 
16  nullptr,
17  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOFIELD
18  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // H1
19  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HCURL
20  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HDIV
21  boost::make_shared<EntitiesFieldData>(MBENTITYSET) // L2
22 
23  },
24  dataOnSlave{
25 
26  nullptr,
27  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOFIELD
28  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // H1
29  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HCURL
30  boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HDIV
31  boost::make_shared<EntitiesFieldData>(MBENTITYSET) // L2
32 
33  },
34  derivedDataOnMaster{
35 
36  nullptr,
37  boost::make_shared<DerivedEntitiesFieldData>(
38  dataOnMaster[NOFIELD]),
39  boost::make_shared<DerivedEntitiesFieldData>(dataOnMaster[H1]),
40  boost::make_shared<DerivedEntitiesFieldData>(
41  dataOnMaster[HCURL]),
42  boost::make_shared<DerivedEntitiesFieldData>(
43  dataOnMaster[HDIV]),
44  boost::make_shared<DerivedEntitiesFieldData>(dataOnMaster[L2])
45 
46  },
47  derivedDataOnSlave{
48 
49  nullptr,
50  boost::make_shared<DerivedEntitiesFieldData>(
51  dataOnSlave[NOFIELD]),
52  boost::make_shared<DerivedEntitiesFieldData>(dataOnSlave[H1]),
53  boost::make_shared<DerivedEntitiesFieldData>(
54  dataOnSlave[HCURL]),
55  boost::make_shared<DerivedEntitiesFieldData>(
56  dataOnSlave[HDIV]),
57  boost::make_shared<DerivedEntitiesFieldData>(dataOnSlave[L2])
58 
59  },
60  dataH1Master(*dataOnMaster[H1].get()),
61  dataH1Slave(*dataOnSlave[H1].get()),
62  dataNoFieldSlave(*dataOnSlave[NOFIELD].get()),
63  dataNoFieldMaster(*dataOnMaster[NOFIELD].get()),
64  dataHcurlMaster(*dataOnMaster[HCURL].get()),
65  dataHcurlSlave(*dataOnSlave[HCURL].get()),
66  dataHdivMaster(*dataOnMaster[HDIV].get()),
67  dataL2Master(*dataOnMaster[L2].get()),
68  dataHdivSlave(*dataOnSlave[HDIV].get()),
69  dataL2Slave(*dataOnSlave[L2].get()), opContravariantTransform() {
70 
71  getUserPolynomialBase() =
72  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
73 
74  // Data on elements for proper spaces
75  dataOnMaster[H1]->setElementType(MBTRI);
76  derivedDataOnMaster[H1]->setElementType(MBTRI);
77  dataOnSlave[H1]->setElementType(MBTRI);
78  derivedDataOnSlave[H1]->setElementType(MBTRI);
79 
80  dataOnMaster[HDIV]->setElementType(MBTRI);
81  derivedDataOnMaster[HDIV]->setElementType(MBTRI);
82  dataOnSlave[HDIV]->setElementType(MBTRI);
83  derivedDataOnSlave[HDIV]->setElementType(MBTRI);
84 
85  dataOnMaster[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
87  dataOnSlave[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
89 
90  derivedDataOnMaster[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
92  derivedDataOnSlave[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
94 
95  CHK_THROW_MESSAGE(createDataOnElement(MBPRISM),
96  "Problem with creation data on element");
97 
98 }
99 
103 
104  if (rule < QUAD_2D_TABLE_SIZE) {
105  if (QUAD_2D_TABLE[rule]->dim != 2) {
106  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
107  }
108  if (QUAD_2D_TABLE[rule]->order < rule) {
109  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
110  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
111  }
113  // For master and slave
114  gaussPtsMaster.resize(3, nbGaussPts, false);
115  gaussPtsSlave.resize(3, nbGaussPts, false);
116 
117  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[1], 3,
118  &gaussPtsMaster(0, 0), 1);
119  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[2], 3,
120  &gaussPtsMaster(1, 0), 1);
121  cblas_dcopy(nbGaussPts, QUAD_2D_TABLE[rule]->weights, 1,
122  &gaussPtsMaster(2, 0), 1);
123 
125 
126  dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
127  false);
128 
129  dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
130  false);
131 
132  double *shape_ptr_master =
133  &*dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
134  cblas_dcopy(3 * nbGaussPts, QUAD_2D_TABLE[rule]->points, 1,
135  shape_ptr_master, 1);
136  double *shape_ptr_slave =
137  &*dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
138  cblas_dcopy(3 * nbGaussPts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr_slave,
139  1);
140  } else {
141  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
142  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
143  nbGaussPts = 0;
144  }
145 
147 }
148 
151 
152  if (numeredEntFiniteElementPtr->getEntType() != MBPRISM)
154 
155  EntitiesFieldData &data_div = *dataOnElement[HDIV];
156  EntitiesFieldData &data_curl = *dataOnElement[HCURL];
158 
160  auto get_coord_and_normal = [&]() {
162  int num_nodes;
163  const EntityHandle *conn;
164  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
165  coords.resize(num_nodes * 3, false);
166  CHKERR mField.get_moab().get_coords(conn, num_nodes,
167  &*coords.data().begin());
168  normal.resize(6, false);
169 
171  &tangentSlaveTwo}) {
172  v->resize(3);
173  v->clear();
174  }
175 
178 
179  auto get_vec_ptr = [](VectorDouble &vec_double, int r = 0) {
181  &vec_double(r + 0), &vec_double(r + 1), &vec_double(r + 2));
182  };
183 
184  auto t_coords_master = get_vec_ptr(coords);
185  auto t_coords_slave = get_vec_ptr(coords, 9);
186  auto t_normal_master = get_vec_ptr(normal);
187  auto t_normal_slave = get_vec_ptr(normal, 3);
188 
189  auto t_t1_master = get_vec_ptr(tangentMasterOne);
190  auto t_t2_master = get_vec_ptr(tangentMasterTwo);
191  auto t_t1_slave = get_vec_ptr(tangentSlaveOne);
192  auto t_t2_slave = get_vec_ptr(tangentSlaveTwo);
193 
194  const double *diff_ptr = Tools::diffShapeFunMBTRI.data();
196  &diff_ptr[0], &diff_ptr[1]);
197 
201 
203 
204  for (int nn = 0; nn != 3; ++nn) {
205  t_t1_master(i) += t_coords_master(i) * t_diff(N0);
206  t_t1_slave(i) += t_coords_slave(i) * t_diff(N0);
207  ++t_coords_master;
208  ++t_coords_slave;
209  ++t_diff;
210  }
211 
212  aRea[0] = sqrt(t_normal_master(i) * t_normal_master(i)) / 2.;
213  aRea[1] = sqrt(t_normal_slave(i) * t_normal_slave(i)) / 2.;
214 
215  t_t2_master(j) =
216  FTensor::levi_civita(i, j, k) * t_normal_master(k) * t_t1_master(i);
217  t_t2_slave(j) =
218  FTensor::levi_civita(i, j, k) * t_normal_slave(k) * t_t1_slave(i);
219 
221  };
222  CHKERR get_coord_and_normal();
223 
225 
226  // H1
227  if ((dataH1.spacesOnEntities[MBVERTEX]).test(H1)) {
228  CHKERR getEntitySense<MBEDGE>(dataH1);
229  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
230  CHKERR getEntitySense<MBTRI>(dataH1);
231  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
232  }
233 
234  // Hcurl
235  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
236  CHKERR getEntitySense<MBEDGE>(data_curl);
237  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
238  CHKERR getEntitySense<MBTRI>(data_curl);
239  CHKERR getEntityDataOrder<MBTRI>(data_curl, HCURL);
240  }
241 
242  // Hdiv
243  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
244  CHKERR getEntitySense<MBTRI>(data_div);
245  CHKERR getEntityDataOrder<MBTRI>(data_div, HDIV);
246  dataHdivSlave.spacesOnEntities[MBTRI].set(HDIV);
248  data_div.spacesOnEntities[MBTRI].set(HDIV);
249  }
250 
251  // L2
252  if ((dataH1.spacesOnEntities[MBTRI]).test(L2)) {
253  CHKERR getEntitySense<MBTRI>(data_l2);
254  CHKERR getEntityDataOrder<MBTRI>(data_l2, L2);
255  }
256 
257  auto clean_data = [](EntitiesFieldData &data) {
259  data.sPace.reset();
260  data.bAse.reset();
261  for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
262  data.spacesOnEntities[t].reset();
263  data.basesOnEntities[t].reset();
264  }
265  for (int s = 0; s != LASTSPACE; ++s)
266  data.basesOnSpaces[s].reset();
267 
269  };
270 
271  auto copy_data = [](EntitiesFieldData &data,
272  EntitiesFieldData &copy_data, const int shift) {
274 
275  if (shift != 0 && shift != 6) {
276  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
277  "Wrong shift for contact prism element");
278  }
279 
280  data.sPace = copy_data.sPace;
281  data.bAse = copy_data.bAse;
282  for (auto t : {MBVERTEX, MBEDGE, MBTRI}) {
283  data.spacesOnEntities[t] = copy_data.spacesOnEntities[t];
284  data.basesOnEntities[t] = copy_data.basesOnEntities[t];
285  data.basesOnSpaces[t] = copy_data.basesOnSpaces[t];
286  }
287 
288  for (int ii = 0; ii != 3; ++ii) {
289  data.dataOnEntities[MBEDGE][ii].getSense() =
290  copy_data.dataOnEntities[MBEDGE][ii + shift].getSense();
291  data.dataOnEntities[MBEDGE][ii].getOrder() =
292  copy_data.dataOnEntities[MBEDGE][ii + shift].getOrder();
293  }
294 
295  if (shift == 0) {
296  data.dataOnEntities[MBTRI][0].getSense() =
297  copy_data.dataOnEntities[MBTRI][3].getSense();
298  data.dataOnEntities[MBTRI][0].getOrder() =
299  copy_data.dataOnEntities[MBTRI][3].getOrder();
300  } else {
301  data.dataOnEntities[MBTRI][0].getSense() =
302  copy_data.dataOnEntities[MBTRI][4].getSense();
303  data.dataOnEntities[MBTRI][0].getOrder() =
304  copy_data.dataOnEntities[MBTRI][4].getOrder();
305  }
306 
308  };
309 
310  auto copy_data_hdiv = [](EntitiesFieldData &data,
311  EntitiesFieldData &copy_data,
312  const int shift) {
314 
315  if (shift != 3 && shift != 4) {
316  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
317  "Wrong shift for contact prism element");
318  }
319 
320  data.sPace = copy_data.sPace;
321  data.bAse = copy_data.bAse;
322 
323  for (auto t : {MBVERTEX, MBTRI}) {
324  data.spacesOnEntities[t] = copy_data.spacesOnEntities[MBVERTEX];
325  data.basesOnEntities[t] = copy_data.basesOnEntities[MBVERTEX];
326  data.basesOnSpaces[t] = copy_data.basesOnSpaces[MBVERTEX];
327  }
328 
329  auto &cpy_ent_dat = copy_data.dataOnEntities[MBTRI][shift];
330  auto &ent_dat = data.dataOnEntities[MBTRI][0];
331  ent_dat.getBase() = cpy_ent_dat.getBase();
332  ent_dat.getSpace() = cpy_ent_dat.getSpace();
333  ent_dat.getSense() = ent_dat.getSense();
334  ent_dat.getOrder() = cpy_ent_dat.getOrder();
335 
337  };
338 
339  CHKERR clean_data(dataH1Slave);
340  CHKERR copy_data(dataH1Slave, dataH1, 6);
341  CHKERR clean_data(dataH1Master);
342  CHKERR copy_data(dataH1Master, dataH1, 0);
343 
344  int order_data = getMaxDataOrder();
345  int order_row = getMaxRowOrder();
346  int order_col = getMaxColOrder(); // maybe two different rules?
347  int rule = getRule(order_row, order_col, order_data);
348 
349  if (rule >= 0) {
350 
352 
353  } else {
354 
355  // Master-Slave
356  if (gaussPtsMaster.size2() != gaussPtsSlave.size2())
357  SETERRQ(
358  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
359  "Number of Gauss Points at Master triangle is different than slave");
360 
361  CHKERR setGaussPts(order_row, order_col, order_data);
362  nbGaussPts = gaussPtsMaster.size2();
363  dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
364  false);
365  dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
366  false);
367 
368  if (nbGaussPts) {
369  CHKERR Tools::shapeFunMBTRI<1>(&*dataH1Master.dataOnEntities[MBVERTEX][0]
370  .getN(NOBASE)
371  .data()
372  .begin(),
373  &gaussPtsMaster(0, 0),
374  &gaussPtsMaster(1, 0), nbGaussPts);
375 
376  CHKERR Tools::shapeFunMBTRI<1>(
377  &*dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
378  &gaussPtsSlave(0, 0), &gaussPtsSlave(1, 0), nbGaussPts);
379  }
380  }
381 
382  if (nbGaussPts == 0)
384 
385  // Get coordinates on slave and master
386  {
387  coordsAtGaussPtsMaster.resize(nbGaussPts, 3, false);
388  coordsAtGaussPtsSlave.resize(nbGaussPts, 3, false);
389  for (int gg = 0; gg < nbGaussPts; gg++) {
390  for (int dd = 0; dd < 3; dd++) {
391  coordsAtGaussPtsMaster(gg, dd) = cblas_ddot(
392  3, &dataH1Master.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
393  &coords[dd], 3);
394  coordsAtGaussPtsSlave(gg, dd) = cblas_ddot(
395  3, &dataH1Slave.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0), 1,
396  &coords[9 + dd], 3);
397  }
398  }
399  }
400 
401  for (int space = HCURL; space != LASTSPACE; ++space)
402  if (dataOnElement[space]) {
403 
404  dataH1Master.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
405  dataOnMaster[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
406  dataH1Slave.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
407  dataOnSlave[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
408 
409  if (space == HDIV) {
410 
411  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
412 
413  CHKERR clean_data(dataHdivSlave);
414  CHKERR copy_data_hdiv(dataHdivSlave, data_div, 4);
415  CHKERR clean_data(dataHdivMaster);
416  CHKERR copy_data_hdiv(dataHdivMaster, data_div, 3);
417 
418  dataHdivMaster.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
419  dataOnMaster[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(
420  NOBASE);
421  dataHdivSlave.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
422  dataOnSlave[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(
423  NOBASE);
424  }
425  }
426  }
427 
428  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
429  if (dataH1.bAse.test(b)) {
430  switch (static_cast<FieldApproximationBase>(b)) {
433  if (dataH1.spacesOnEntities[MBVERTEX].test(H1)) {
434  CHKERR getUserPolynomialBase()->getValue(
436  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
437  dataH1Master, H1, static_cast<FieldApproximationBase>(b),
438  NOBASE)));
439 
440  CHKERR getUserPolynomialBase()->getValue(
442  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
443  dataH1Slave, H1, static_cast<FieldApproximationBase>(b),
444  NOBASE)));
445  }
446  break;
448  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
449 
450  CHKERR getUserPolynomialBase()->getValue(
452  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
453  dataHdivMaster, HDIV, static_cast<FieldApproximationBase>(b),
454  NOBASE)));
455  CHKERR getUserPolynomialBase()->getValue(
457  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
458  dataHdivSlave, HDIV, static_cast<FieldApproximationBase>(b),
459  NOBASE)));
460 
466 
467  }
468 
469  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
470  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
471  "Not yet implemented");
472  }
473  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
474  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
475  "Not yet implemented");
476  }
477  break;
478  default:
479  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
480  "Not yet implemented");
481  }
482  }
483  }
484 
485  // Iterate over operators
487 
489 }
490 
493 
494  constexpr std::array<UserDataOperator::OpType, 2> types{
496  std::array<std::string, 2> last_eval_field_name{std::string(), std::string()};
497 
498  auto oit = opPtrVector.begin();
499  auto hi_oit = opPtrVector.end();
500 
501  for (; oit != hi_oit; oit++) {
502 
503  oit->setPtrFE(this);
504 
505  if (oit->opType == UserDataOperator::OPSPACE) {
506 
507  // Set field
508  switch (oit->sPace) {
509  case NOSPACE:
510  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Unknown space");
511  case NOFIELD:
512  case H1:
513  case HCURL:
514  case HDIV:
515  case L2:
516  break;
517  default:
518  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
519  "Not implemented for this space", oit->sPace);
520  }
521 
522  // Reseat all data which all field dependent
523  dataOnMaster[oit->sPace]->resetFieldDependentData();
524  dataOnSlave[oit->sPace]->resetFieldDependentData();
525 
526  // Run operator
527  try {
528  CHKERR oit->opRhs(*dataOnMaster[oit->sPace], false);
529  CHKERR oit->opRhs(*dataOnSlave[oit->sPace], false);
530  }
531  CATCH_OP_ERRORS(*oit);
532 
533  } else {
534  boost::shared_ptr<EntitiesFieldData> op_master_data[2];
535  boost::shared_ptr<EntitiesFieldData> op_slave_data[2];
536 
537  for (int ss = 0; ss != 2; ss++) {
538 
539  const std::string field_name =
540  !ss ? oit->rowFieldName : oit->colFieldName;
541  const Field *field_struture = mField.get_field_structure(field_name);
542  const BitFieldId data_id = field_struture->getId();
543  const FieldSpace space = field_struture->getSpace();
544  op_master_data[ss] =
545  !ss ? dataOnMaster[space] : derivedDataOnMaster[space];
546  op_slave_data[ss] =
547  !ss ? dataOnSlave[space] : derivedDataOnSlave[space];
548 
549  if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
550  data_id)
551  .none())
552  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
553  "no data field < %s > on finite element < %s >",
554  field_name.c_str(), getFEName().c_str());
555 
556  if (oit->getOpType() & types[ss] ||
557  oit->getOpType() & UserDataOperator::OPROWCOL) {
558 
559  switch (space) {
560  case NOSPACE:
561  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown space");
562  break;
563  case NOFIELD:
564  case H1:
565  case HCURL:
566  case HDIV:
567  case L2:
568  break;
569  default:
570  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
571  "Not implemented for this space", space);
572  }
573 
574  if (last_eval_field_name[ss] != field_name) {
575  CHKERR getEntityFieldData(*op_master_data[ss], *op_slave_data[ss],
576  field_name, MBEDGE);
577 
578  if (!ss) {
579  struct Extractor {
580  boost::weak_ptr<EntityCacheNumeredDofs>
581  operator()(boost::shared_ptr<FieldEntity> &e) {
582  return e->entityCacheRowDofs;
583  }
584  };
585 
586  CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
587  field_name, getRowFieldEnts(), MBEDGE,
588  MBPRISM, Extractor());
589  } else {
590  struct Extractor {
591  boost::weak_ptr<EntityCacheNumeredDofs>
592  operator()(boost::shared_ptr<FieldEntity> &e) {
593  return e->entityCacheColDofs;
594  }
595  };
596  CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
597  field_name, getRowFieldEnts(), MBEDGE,
598  MBPRISM, Extractor());
599  }
600 
601  switch (space) {
602  case NOSPACE:
603  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
604  "unknown space");
605  break;
606  case H1: {
607 
608  auto get_indices = [&](auto &master, auto &slave, auto &ents,
609  auto &&ex) {
610  return getNodesIndices(
611  field_name, ents,
612  master.dataOnEntities[MBVERTEX][0].getIndices(),
613  master.dataOnEntities[MBVERTEX][0].getLocalIndices(),
614  slave.dataOnEntities[MBVERTEX][0].getIndices(),
615  slave.dataOnEntities[MBVERTEX][0].getLocalIndices(), ex);
616  };
617 
618  auto get_data = [&](EntitiesFieldData &master_data,
619  EntitiesFieldData &slave_data) {
620  return getNodesFieldData(
621  field_name,
622  master_data.dataOnEntities[MBVERTEX][0].getFieldData(),
623  slave_data.dataOnEntities[MBVERTEX][0].getFieldData(),
624  master_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
625  slave_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
626  master_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
627  slave_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
628  master_data.dataOnEntities[MBVERTEX][0].getSpace(),
629  slave_data.dataOnEntities[MBVERTEX][0].getSpace(),
630  master_data.dataOnEntities[MBVERTEX][0].getBase(),
631  slave_data.dataOnEntities[MBVERTEX][0].getBase());
632  };
633 
634  if (!ss) {
635 
636  struct Extractor {
637  boost::weak_ptr<EntityCacheNumeredDofs>
638  operator()(boost::shared_ptr<FieldEntity> &e) {
639  return e->entityCacheRowDofs;
640  }
641  };
642 
643  CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
644  getRowFieldEnts(), Extractor());
645  } else {
646  struct Extractor {
647  boost::weak_ptr<EntityCacheNumeredDofs>
648  operator()(boost::shared_ptr<FieldEntity> &e) {
649  return e->entityCacheColDofs;
650  }
651  };
652 
653  CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
654  getColFieldEnts(), Extractor());
655  }
656 
657  CHKERR get_data(*op_master_data[ss], *op_slave_data[ss]);
658 
659  } break;
660  case HCURL:
661  case HDIV:
662  case L2:
663  break;
664  case NOFIELD:
665  if (!getNinTheLoop()) {
666  // NOFIELD data are the same for each element, can be
667  // retrieved only once
668  const auto bit_number = mField.get_field_bit_number(field_name);
669  if (!ss) {
670  CHKERR getNoFieldRowIndices(*op_master_data[ss], bit_number);
671  CHKERR getNoFieldRowIndices(*op_slave_data[ss], bit_number);
672  } else {
673  CHKERR getNoFieldColIndices(*op_master_data[ss], bit_number);
674  CHKERR getNoFieldColIndices(*op_slave_data[ss], bit_number);
675  }
676  CHKERR getNoFieldFieldData(*op_master_data[ss], bit_number);
677  CHKERR getNoFieldFieldData(*op_slave_data[ss], bit_number);
678  }
679  break;
680  default:
681  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
682  "Not implemented for this space", space);
683  }
684  last_eval_field_name[ss] = field_name;
685  }
686  }
687  }
688 
689  int type;
690 
691  if (UserDataOperator *cast_oit =
692  dynamic_cast<UserDataOperator *>(&*oit)) {
693  type = cast_oit->getFaceType();
694  if (((oit->getOpType() & UserDataOperator::OPROW) ||
695  (oit->getOpType() & UserDataOperator::OPCOL)) &&
700  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
701  "Wrong combination of FaceType and OpType, OPROW or OPCOL "
702  "combined with face-face OpType");
703  }
704 
705  if ((oit->getOpType() & UserDataOperator::OPROWCOL) &&
708  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
709  "Wrong combination of FaceType and OpType, OPROWCOL "
710  "combined with face-face OpType");
711  }
712 
713  if (!type) {
714  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
715  "Face type is not set");
716  }
717  } else {
723  }
724 
725  if (oit->getOpType() & UserDataOperator::OPROW &&
727  try {
728  CHKERR oit->opRhs(*op_master_data[0], false);
729  }
730  CATCH_OP_ERRORS(*oit);
731  }
732 
733  if (oit->getOpType() & UserDataOperator::OPROW &&
735  try {
736  CHKERR oit->opRhs(*op_slave_data[0], false);
737  }
738  CATCH_OP_ERRORS(*oit);
739  }
740 
741  if (oit->getOpType() & UserDataOperator::OPCOL &&
743  try {
744  CHKERR oit->opRhs(*op_master_data[1], false);
745  }
746  CATCH_OP_ERRORS(*oit);
747  }
748 
749  if (oit->getOpType() & UserDataOperator::OPCOL &&
751  try {
752  CHKERR oit->opRhs(*op_slave_data[1], false);
753  }
754  CATCH_OP_ERRORS(*oit);
755  }
756 
757  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
759  try {
760  CHKERR oit->opLhs(*op_master_data[0], *op_master_data[1]);
761  }
762  CATCH_OP_ERRORS(*oit);
763  }
764 
765  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
767  try {
768  CHKERR oit->opLhs(*op_master_data[0], *op_slave_data[1]);
769  }
770  CATCH_OP_ERRORS(*oit);
771  }
772 
773  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
775  try {
776  CHKERR oit->opLhs(*op_slave_data[0], *op_master_data[1]);
777  }
778  CATCH_OP_ERRORS(*oit);
779  }
780 
781  if (oit->getOpType() & UserDataOperator::OPROWCOL &&
783  try {
784  CHKERR oit->opLhs(*op_slave_data[0], *op_slave_data[1]);
785  }
786  CATCH_OP_ERRORS(*oit);
787  }
788  }
789  }
791 }
792 
794  EntitiesFieldData &master_data, EntitiesFieldData &slave_data,
795  const std::string &field_name, const EntityType type_lo,
796  const EntityType type_hi) const {
798 
799  auto reset_data = [type_lo, type_hi](auto &data) {
800  for (EntityType t = type_lo; t != type_hi; ++t) {
801  for (auto &dat : data.dataOnEntities[t]) {
802  dat.getOrder() = 0;
803  dat.getBase() = NOBASE;
804  dat.getSpace() = NOSPACE;
805  dat.getFieldData().resize(0, false);
806  dat.getFieldDofs().resize(0, false);
807  dat.getFieldEntities().resize(0, false);
808  }
809  }
810  };
811  reset_data(master_data);
812  reset_data(slave_data);
813 
814  auto &field_ents = getDataFieldEnts();
815  auto bit_number = mField.get_field_bit_number(field_name);
816  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
817  bit_number, get_id_for_min_type(type_lo));
818  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
819  cmp_uid_lo);
820  if (lo != field_ents.end()) {
821  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
822  bit_number, get_id_for_max_type(type_hi));
823  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
824  if (lo != hi) {
825  for (auto it = lo; it != hi; ++it)
826  if (auto e = it->lock()) {
827 
828  auto get_data = [&](auto &data, auto type, auto side) {
829  auto &dat = data.dataOnEntities[type][side];
830  auto &ent_field_dofs = dat.getFieldDofs();
831  auto &ent_field_data = dat.getFieldData();
832  dat.getFieldEntities().resize(1, false);
833  dat.getFieldEntities()[0] = e.get();
834  dat.getBase() = e->getApproxBase();
835  dat.getSpace() = e->getSpace();
836  const int ent_order = e->getMaxOrder();
837  dat.getOrder() =
838  dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
839  const auto dof_ent_field_data = e->getEntFieldData();
840  const int nb_dofs_on_ent = e->getNbDofsOnEnt();
841  ent_field_data.resize(nb_dofs_on_ent, false);
842  noalias(ent_field_data) = e->getEntFieldData();
843  ent_field_dofs.resize(nb_dofs_on_ent, false);
844  std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
845  if (auto cache = e->entityCacheDataDofs.lock()) {
846  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
847  ent_field_dofs[(*dit)->getEntDofIdx()] =
848  reinterpret_cast<FEDofEntity *>((*dit).get());
849  }
850  }
851  };
852 
853  const EntityType type = e->getEntType();
854  const int side = e->getSideNumberPtr()->side_number;
855 
856  switch (type) {
857  case MBEDGE:
858 
859  if (side < 3)
860  get_data(master_data, type, side);
861  else if (side > 5)
862  get_data(slave_data, type, side - 6);
863 
864  break;
865  case MBTRI:
866 
867  if (side == 3)
868  get_data(master_data, type, 0);
869  if (side == 4)
870  get_data(slave_data, type, 0);
871 
872  break;
873  default:
874  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
875  "Entity type not implemented (FIXME)");
876  };
877 
878  const int brother_side = e->getSideNumberPtr()->brother_side_number;
879  if (brother_side != -1)
880  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
881  "Case with brother side not implemented (FIXME)");
882  }
883  }
884  }
885 
887 }
888 
890  const std::string field_name, VectorDouble &master_nodes_data,
891  VectorDouble &slave_nodes_data, VectorDofs &master_nodes_dofs,
892  VectorDofs &slave_nodes_dofs,
893 
894  VectorFieldEntities &master_field_entities, VectorFieldEntities &slave_field_entities,
895 
896  FieldSpace &master_space,
897  FieldSpace &slave_space, FieldApproximationBase &master_base,
898  FieldApproximationBase &slave_base) const {
900 
901  auto set_zero = [](auto &nodes_data, auto &nodes_dofs, auto &field_entities) {
902  nodes_data.resize(0, false);
903  nodes_dofs.resize(0, false);
904  field_entities.resize(0, false);
905  };
906  set_zero(master_nodes_data, master_nodes_dofs, master_field_entities);
907  set_zero(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
908 
909  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
910  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
911 
912  auto bit_number = (*field_it)->getBitNumber();
913  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
914  master_space = slave_space = (*field_it)->getSpace();
915  master_base = slave_base = (*field_it)->getApproxBase();
916 
917  auto &field_ents = getDataFieldEnts();
918  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
919  bit_number, get_id_for_min_type<MBVERTEX>());
920  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
921  cmp_uid_lo);
922  if (lo != field_ents.end()) {
923  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
924  bit_number, get_id_for_max_type<MBVERTEX>());
925  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
926  if (lo != hi) {
927 
928  int nb_dofs = 0;
929  for (auto it = lo; it != hi; ++it) {
930  if (auto e = it->lock()) {
931  nb_dofs += e->getEntFieldData().size();
932  }
933  }
934 
935  if (nb_dofs) {
936 
937  auto init_set = [&](auto &nodes_data, auto &nodes_dofs,
938  auto &field_entities) {
939  constexpr int num_nodes = 3;
940  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
941  nodes_data.resize(max_nb_dofs, false);
942  nodes_dofs.resize(max_nb_dofs, false);
943  field_entities.resize(num_nodes, false);
944  nodes_data.clear();
945  fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
946  fill(field_entities.begin(), field_entities.end(), nullptr);
947  };
948 
949  init_set(master_nodes_data, master_nodes_dofs, master_field_entities);
950  init_set(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
951 
952  for (auto it = lo; it != hi; ++it) {
953  if (auto e = it->lock()) {
954 
955  const auto &sn = e->getSideNumberPtr();
956  int side = sn->side_number;
957 
958  auto set_data = [&](auto &nodes_data, auto &nodes_dofs,
959  auto &field_entities, int side, int pos) {
960  field_entities[side] = e.get();
961  if (auto cache = e->entityCacheDataDofs.lock()) {
962  for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
963  ++dit) {
964  const auto dof_idx = (*dit)->getEntDofIdx();
965  nodes_data[pos + dof_idx] = (*dit)->getFieldData();
966  nodes_dofs[pos + dof_idx] =
967  reinterpret_cast<FEDofEntity *>((*dit).get());
968  }
969  }
970  };
971 
972  if (side < 3)
973  set_data(master_nodes_data, master_nodes_dofs,
974  master_field_entities, side, side * nb_dofs_on_vert);
975  else
976  set_data(slave_nodes_data, slave_nodes_dofs,
977  slave_field_entities, (side - 3),
978  (side - 3) * nb_dofs_on_vert);
979 
980  const int brother_side = sn->brother_side_number;
981  if (brother_side != -1)
982  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
983  "Not implemented (FIXME please)");
984  }
985  }
986  }
987  }
988  }
989  }
990 
992 }
993 
994 template <typename EXTRACTOR>
996  EntitiesFieldData &master_data, EntitiesFieldData &slave_data,
997  const std::string &field_name, FieldEntity_vector_view &ents_field,
998  const EntityType type_lo, const EntityType type_hi,
999  EXTRACTOR &&extractor) const {
1001 
1002  auto clear_data = [type_lo, type_hi](auto &data) {
1003  for (EntityType t = type_lo; t != type_hi; ++t) {
1004  for (auto &d : data.dataOnEntities[t]) {
1005  d.getIndices().resize(0, false);
1006  d.getLocalIndices().resize(0, false);
1007  }
1008  }
1009  };
1010 
1011  clear_data(master_data);
1012  clear_data(slave_data);
1013 
1014  auto bit_number = mField.get_field_bit_number(field_name);
1015  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1016  bit_number, get_id_for_min_type(type_lo));
1017  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1018  cmp_uid_lo);
1019  if (lo != ents_field.end()) {
1020  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1021  bit_number, get_id_for_max_type(type_hi));
1022  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1023  if (lo != hi) {
1024 
1025  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1026 
1027  for (auto it = lo; it != hi; ++it)
1028  if (auto e = it->lock()) {
1029 
1030  const EntityType type = e->getEntType();
1031  const int side = e->getSideNumberPtr()->side_number;
1032  const int brother_side = e->getSideNumberPtr()->brother_side_number;
1033  if (brother_side != -1)
1034  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1035  "Not implemented case");
1036 
1037  auto get_indices = [&](auto &data, const auto type, const auto side) {
1038  if (auto cache = extractor(e).lock()) {
1039  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
1040  auto &dof = (**dit);
1041  auto &dat = data.dataOnEntities[type][side];
1042  auto &ent_field_indices = dat.getIndices();
1043  auto &ent_field_local_indices = dat.getLocalIndices();
1044  if (ent_field_indices.empty()) {
1045  const int nb_dofs_on_ent = dof.getNbDofsOnEnt();
1046  ent_field_indices.resize(nb_dofs_on_ent, false);
1047  ent_field_local_indices.resize(nb_dofs_on_ent, false);
1048  std::fill(ent_field_indices.data().begin(),
1049  ent_field_indices.data().end(), -1);
1050  std::fill(ent_field_local_indices.data().begin(),
1051  ent_field_local_indices.data().end(), -1);
1052  }
1053  const int idx = dof.getEntDofIdx();
1054  ent_field_indices[idx] = dof.getPetscGlobalDofIdx();
1055  ent_field_local_indices[idx] = dof.getPetscLocalDofIdx();
1056  }
1057  }
1058  };
1059 
1060  switch (type) {
1061  case MBEDGE:
1062 
1063  if (side < 3)
1064  get_indices(master_data, type, side);
1065  else if (side > 5)
1066  get_indices(slave_data, type, side - 6);
1067 
1068  break;
1069  case MBTRI:
1070 
1071  if (side == 3)
1072  get_indices(master_data, type, 0);
1073  if (side == 4)
1074  get_indices(slave_data, type, 0);
1075 
1076  break;
1077  default:
1078  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1079  "Entity type not implemented");
1080  }
1081  }
1082  }
1083  }
1084 
1086 }
1087 
1088 template <typename EXTRACTOR>
1090  const std::string field_name, FieldEntity_vector_view &ents_field,
1091  VectorInt &master_nodes_indices, VectorInt &master_local_nodes_indices,
1092  VectorInt &slave_nodes_indices, VectorInt &slave_local_nodes_indices,
1093  EXTRACTOR &&extractor) const {
1095 
1096  master_nodes_indices.resize(0, false);
1097  master_local_nodes_indices.resize(0, false);
1098  slave_nodes_indices.resize(0, false);
1099  slave_local_nodes_indices.resize(0, false);
1100 
1101  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
1102  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
1103 
1104  auto bit_number = (*field_it)->getBitNumber();
1105  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
1106 
1107  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1108  bit_number, get_id_for_min_type<MBVERTEX>());
1109  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1110  cmp_uid_lo);
1111  if (lo != ents_field.end()) {
1112  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1113  bit_number, get_id_for_max_type<MBVERTEX>());
1114  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1115  if (lo != hi) {
1116 
1117  int nb_dofs = 0;
1118  for (auto it = lo; it != hi; ++it) {
1119  if (auto e = it->lock()) {
1120  if (auto cache = extractor(e).lock()) {
1121  nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
1122  }
1123  }
1124  }
1125 
1126  if (nb_dofs) {
1127 
1128  constexpr int num_nodes = 3;
1129  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
1130 
1131  auto set_vec_size = [&](auto &nodes_indices,
1132  auto &local_nodes_indices) {
1133  nodes_indices.resize(max_nb_dofs, false);
1134  local_nodes_indices.resize(max_nb_dofs, false);
1135  std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
1136  std::fill(local_nodes_indices.begin(), local_nodes_indices.end(),
1137  -1);
1138  };
1139 
1140  set_vec_size(master_nodes_indices, master_local_nodes_indices);
1141  set_vec_size(slave_nodes_indices, slave_local_nodes_indices);
1142 
1143  for (auto it = lo; it != hi; ++it) {
1144  if (auto e = it->lock()) {
1145 
1146  const int side = e->getSideNumberPtr()->side_number;
1147  const int brother_side =
1148  e->getSideNumberPtr()->brother_side_number;
1149  if (brother_side != -1)
1150  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1151  "Not implemented case");
1152 
1153  auto get_indices = [&](auto &nodes_indices,
1154  auto &local_nodes_indices,
1155  const auto side) {
1156  if (auto cache = extractor(e).lock()) {
1157  for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
1158  ++dit) {
1159  const int idx = (*dit)->getPetscGlobalDofIdx();
1160  const int local_idx = (*dit)->getPetscLocalDofIdx();
1161  const int pos =
1162  side * nb_dofs_on_vert + (*dit)->getDofCoeffIdx();
1163  nodes_indices[pos] = idx;
1164  local_nodes_indices[pos] = local_idx;
1165  }
1166  }
1167  };
1168 
1169  if (side < 3)
1170  get_indices(master_nodes_indices, master_local_nodes_indices,
1171  side);
1172  else if (side > 2)
1173  get_indices(slave_nodes_indices, slave_local_nodes_indices,
1174  side - 3);
1175  else
1176  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1177  "Impossible case");
1178  }
1179  }
1180  }
1181  }
1182  }
1183  }
1184 
1186 }
1187 
1190  const string fe_name,
1192  const int side_type, const EntityHandle ent_for_side) {
1193  return loopSide(fe_name, &fe_method, side_type, ent_for_side);
1194 }
1195 
1196 } // namespace MoFEM
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::FEMethod::getRowFieldEnts
auto & getRowFieldEnts() const
Definition: LoopMethods.hpp:423
MoFEM::BasicMethod::getNinTheLoop
int getNinTheLoop() const
get number of evaluated element in the loop
Definition: LoopMethods.hpp:207
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentSlaveOne
VectorDouble tangentSlaveOne
Definition: ContactPrismElementForcesAndSourcesCore.hpp:68
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentMasterOne
VectorDouble tangentMasterOne
Definition: ContactPrismElementForcesAndSourcesCore.hpp:69
MoFEM::ContactPrismElementForcesAndSourcesCore::derivedDataOnSlave
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:92
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTERSLAVE
@ FACEMASTERSLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:255
CATCH_OP_ERRORS
#define CATCH_OP_ERRORS(OP)
Definition: ForcesAndSourcesCore.cpp:1323
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVESLAVE
@ FACESLAVESLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:257
MoFEM::ForcesAndSourcesCore::getNoFieldFieldData
MoFEMErrorCode getNoFieldFieldData(const int bit_number, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
Definition: ForcesAndSourcesCore.cpp:855
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:140
LASTBASE
@ LASTBASE
Definition: definitions.h:69
MoFEM::ForcesAndSourcesCore::getMaxRowOrder
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Definition: ForcesAndSourcesCore.cpp:131
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::ContactPrismElementForcesAndSourcesCore::opContravariantTransform
OpSetContravariantPiolaTransformOnFace opContravariantTransform
Definition: ContactPrismElementForcesAndSourcesCore.hpp:70
EntityHandle
MoFEM::ContactPrismElementForcesAndSourcesCore::derivedDataOnMaster
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnMaster
Entity data on element entity columns fields.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:90
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPSPACE
@ OPSPACE
operator do Work is execute on space data
Definition: ForcesAndSourcesCore.hpp:570
MoFEM::EntitiesFieldData::basesOnEntities
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
Definition: EntitiesFieldData.hpp:52
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::ContactPrismElementForcesAndSourcesCore::ContactPrismElementForcesAndSourcesCore
ContactPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: ContactPrismElementForcesAndSourcesCore.cpp:12
MoFEM::cmp_uid_hi
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
Definition: ForcesAndSourcesCore.cpp:29
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpSetContravariantPiolaTransformOnFace::normalRawPtr
const VectorDouble * normalRawPtr
Definition: DataOperators.hpp:419
MoFEM::Field::getSpace
FieldSpace getSpace() const
Get field approximation space.
Definition: FieldMultiIndices.hpp:150
MoFEM::FieldEntity_vector_view
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Definition: FieldEntsMultiIndices.hpp:478
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1926
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentMasterTwo
VectorDouble tangentMasterTwo
Definition: ContactPrismElementForcesAndSourcesCore.hpp:69
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::FEMethod::getDataFieldEnts
const FieldEntity_vector_view & getDataFieldEnts() const
Definition: LoopMethods.hpp:412
MoFEM::ForcesAndSourcesCore::getMaxDataOrder
int getMaxDataOrder() const
Get max order of approximation for data fields.
Definition: ForcesAndSourcesCore.cpp:120
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
MoFEM::ContactPrismElementForcesAndSourcesCore::dataHdivMaster
EntitiesFieldData & dataHdivMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:101
MoFEM::ContactPrismElementForcesAndSourcesCore::getEntityFieldData
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &master_data, EntitiesFieldData &slave_data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
function that gets entity field data.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:793
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:67
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
MoFEM::ForcesAndSourcesCore::getMaxColOrder
int getMaxColOrder() const
Get max order of approximation for field in columns.
Definition: ForcesAndSourcesCore.cpp:135
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
MoFEM::ForcesAndSourcesCore::getNoFieldRowIndices
MoFEMErrorCode getNoFieldRowIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:459
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:50
sdf.r
int r
Definition: sdf.py:8
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::basesOnSpaces
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
Definition: EntitiesFieldData.hpp:54
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:252
MoFEM::ContactPrismElementForcesAndSourcesCore::dataOnSlave
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:81
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
MoFEM::Field::getId
const BitFieldId & getId() const
Get unique field id.
Definition: FieldMultiIndices.hpp:128
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnContactPrismSide &fe_method, const int side_type, const EntityHandle ent_for_side)
Definition: ContactPrismElementForcesAndSourcesCore.cpp:1189
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::ContactPrismElementForcesAndSourcesCore::getNodesIndices
MoFEMErrorCode getNodesIndices(const std::string field_name, FieldEntity_vector_view &ents_field, VectorInt &master_nodes_indices, VectorInt &master_local_nodes_indices, VectorInt &slave_nodes_indices, VectorInt &slave_local_nodes_indices, EXTRACTOR &&extractor) const
function that gets nodes indices.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:1089
MoFEM::ForcesAndSourcesCore::getNoFieldColIndices
MoFEMErrorCode getNoFieldColIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:474
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVEMASTER
@ FACESLAVEMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:256
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::ContactPrismElementForcesAndSourcesCore::setDefaultGaussPts
MoFEMErrorCode setDefaultGaussPts(const int rule)
Definition: ContactPrismElementForcesAndSourcesCore.cpp:101
MoFEM::ForcesAndSourcesCore::getUserPolynomialBase
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:97
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::VolumeElementForcesAndSourcesCoreOnContactPrismSide
Base volume element used to integrate on contact surface (could be extended to other volume elements)
Definition: VolumeElementForcesAndSourcesCoreOnContactPrismSide.hpp:23
MoFEM::ContactPrismElementForcesAndSourcesCore::getNodesFieldData
MoFEMErrorCode getNodesFieldData(const std::string field_name, VectorDouble &master_nodes_data, VectorDouble &slave_nodes_data, VectorDofs &master_nodes_dofs, VectorDofs &slave_nodes_dofs, VectorFieldEntities &master_field_entities, VectorFieldEntities &slave_field_entities, FieldSpace &master_space, FieldSpace &slave_space, FieldApproximationBase &master_base, FieldApproximationBase &slave_base) const
function that gets nodes field data.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:889
MoFEM::ContactPrismElementForcesAndSourcesCore::normal
VectorDouble normal
vector storing vector normal to master or slave element
Definition: ContactPrismElementForcesAndSourcesCore.hpp:56
convert.type
type
Definition: convert.py:64
MoFEM::FieldEntity::getLocalUniqueIdCalculate
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Definition: FieldEntsMultiIndices.hpp:136
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPCOL
@ OPCOL
operator doWork function is executed on FE columns
Definition: ForcesAndSourcesCore.hpp:568
MoFEM::ForcesAndSourcesCore::opPtrVector
boost::ptr_deque< UserDataOperator > opPtrVector
Vector of finite element users data operators.
Definition: ForcesAndSourcesCore.hpp:480
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:461
MoFEM::ContactPrismElementForcesAndSourcesCore::tangentSlaveTwo
VectorDouble tangentSlaveTwo
Definition: ContactPrismElementForcesAndSourcesCore.hpp:68
MoFEM::cmp_uid_lo
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
Definition: ForcesAndSourcesCore.cpp:18
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
MoFEM::EntitiesFieldData::sPace
std::bitset< LASTSPACE > sPace
spaces on element
Definition: EntitiesFieldData.hpp:42
MoFEM::ContactPrismElementForcesAndSourcesCore::gaussPtsMaster
MatrixDouble gaussPtsMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:63
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator
default operator for Contact Prism element
Definition: ContactPrismElementForcesAndSourcesCore.hpp:208
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::ContactPrismElementForcesAndSourcesCore::coordsAtGaussPtsSlave
MatrixDouble coordsAtGaussPtsSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:60
LASTSPACE
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:89
MoFEM::ForcesAndSourcesCore::getRule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Definition: ForcesAndSourcesCore.cpp:1920
MoFEM::get_id_for_max_type
EntityHandle get_id_for_max_type()
Definition: RefEntsMultiIndices.hpp:13
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::FEMethod::getFEName
auto getFEName() const
get finite element name
Definition: LoopMethods.hpp:391
MoFEM::ContactPrismElementForcesAndSourcesCore::dataH1Slave
EntitiesFieldData & dataH1Slave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:95
MoFEM::ContactPrismElementForcesAndSourcesCore::getEntityIndices
MoFEMErrorCode getEntityIndices(EntitiesFieldData &master_data, EntitiesFieldData &slave_data, const std::string &field_name, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
function that gets entity indices.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:995
MoFEM::ContactPrismElementForcesAndSourcesCore::coords
VectorDouble coords
Definition: ContactPrismElementForcesAndSourcesCore.hpp:57
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEM::EntitiesFieldData::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: EntitiesFieldData.hpp:45
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:990
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::ForcesAndSourcesCore::dataH1
EntitiesFieldData & dataH1
Definition: ForcesAndSourcesCore.hpp:471
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
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
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::ContactPrismElementForcesAndSourcesCore::nbGaussPts
int nbGaussPts
Definition: ContactPrismElementForcesAndSourcesCore.hpp:202
MoFEM::ForcesAndSourcesCore::UserDataOperator::loopSide
MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy, AdjCache *adj_cache=nullptr)
User calls this function to loop over elements on the side of face. This function calls finite elemen...
Definition: ForcesAndSourcesCore.cpp:1668
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::ContactPrismElementForcesAndSourcesCore::dataH1Master
EntitiesFieldData & dataH1Master
Definition: ContactPrismElementForcesAndSourcesCore.hpp:94
MoFEM::ContactPrismElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: ContactPrismElementForcesAndSourcesCore.cpp:149
MoFEM::ContactPrismElementForcesAndSourcesCore::coordsAtGaussPtsMaster
MatrixDouble coordsAtGaussPtsMaster
Definition: ContactPrismElementForcesAndSourcesCore.hpp:58
MoFEM::ContactPrismElementForcesAndSourcesCore::dataOnMaster
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnMaster
Entity data on element entity rows fields.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:80
MoFEM::OpSetContravariantPiolaTransformOnFace::normalShift
int normalShift
Shift in vector for linear geometry.
Definition: DataOperators.hpp:431
UBlasVector< double >
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
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::FEMethod::getColFieldEnts
auto & getColFieldEnts() const
Definition: LoopMethods.hpp:431
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::ContactPrismElementForcesAndSourcesCore::dataHdivSlave
EntitiesFieldData & dataHdivSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:102
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
MoFEM::VectorFieldEntities
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
Definition: EntitiesFieldData.hpp:31
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:383
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTERMASTER
@ FACEMASTERMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:254
MoFEM::field_it
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
Definition: Projection10NodeCoordsOnField.cpp:123
MoFEM::FEDofEntity
keeps information about indexed dofs for the finite element
Definition: DofsMultiIndices.hpp:267
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::ContactPrismElementForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ContactPrismElementForcesAndSourcesCore.cpp:491
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::get_id_for_min_type
EntityHandle get_id_for_min_type()
Definition: RefEntsMultiIndices.hpp:18
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::ContactPrismElementForcesAndSourcesCore::aRea
std::array< double, 2 > aRea
Array storing master and slave faces areas.
Definition: ContactPrismElementForcesAndSourcesCore.hpp:53
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:353
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
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::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACESLAVE
@ FACESLAVE
Definition: ContactPrismElementForcesAndSourcesCore.hpp:253
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
MoFEM::ContactPrismElementForcesAndSourcesCore::gaussPtsSlave
MatrixDouble gaussPtsSlave
Definition: ContactPrismElementForcesAndSourcesCore.hpp:65
MoFEM::ContactPrismElementForcesAndSourcesCore::UserDataOperator::FACEMASTER
@ FACEMASTER
Definition: ContactPrismElementForcesAndSourcesCore.hpp:252