v0.14.0
Loading...
Searching...
No Matches
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
9namespace 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
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
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
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);
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(
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{
495 UserDataOperator::OPROW, UserDataOperator::OPCOL};
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(
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)) &&
696 ((type & UserDataOperator::FACEMASTERMASTER) ||
697 (type & UserDataOperator::FACEMASTERSLAVE) ||
698 (type & UserDataOperator::FACESLAVEMASTER) ||
699 (type & UserDataOperator::FACESLAVESLAVE))) {
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) &&
706 ((type & UserDataOperator::FACEMASTER) ||
707 (type & UserDataOperator::FACESLAVE))) {
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 {
718 type = UserDataOperator::FACEMASTER | UserDataOperator::FACESLAVE |
719 UserDataOperator::FACEMASTERMASTER |
720 UserDataOperator::FACEMASTERSLAVE |
721 UserDataOperator::FACESLAVEMASTER |
722 UserDataOperator::FACESLAVESLAVE;
723 }
724
725 if (oit->getOpType() & UserDataOperator::OPROW &&
726 (type & UserDataOperator::FACEMASTER)) {
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 &&
734 (type & UserDataOperator::FACESLAVE)) {
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 &&
742 (type & UserDataOperator::FACEMASTER)) {
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 &&
750 (type & UserDataOperator::FACESLAVE)) {
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 &&
758 (type & UserDataOperator::FACEMASTERMASTER)) {
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 &&
766 (type & UserDataOperator::FACEMASTERSLAVE)) {
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 &&
774 (type & UserDataOperator::FACESLAVEMASTER)) {
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 &&
782 (type & UserDataOperator::FACESLAVESLAVE)) {
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
994template <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
1088template <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
1189ContactPrismElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes(
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
static Number< 0 > N0
#define CATCH_OP_ERRORS(OP)
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ NOBASE
Definition: definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:89
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
const int dim
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
BlockParamData * cache
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
EntityHandle get_id_for_max_type()
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
EntityHandle get_id_for_min_type()
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
const Field_multiIndex * fieldsPtr
raw pointer to fields container
int getNinTheLoop() const
get number of evaluated element in the loop
MoFEMErrorCode operator()()
function is run for every finite element
VectorDouble normal
vector storing vector normal to master or slave element
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.
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.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnSlave
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnMaster
Entity data on element entity rows fields.
std::array< double, 2 > aRea
Array storing master and slave faces areas.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnMaster
Entity data on element entity columns fields.
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.
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.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnSlave
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Deprecated interface functions.
this class derive data form other data structure
Class used to pass element data to calculate base functions on tet,triangle,edge.
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::bitset< LASTSPACE > sPace
spaces on element
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
keeps information about indexed dofs for the finite element
auto & getRowFieldEnts() const
auto getFEName() const
get finite element name
const FieldEntity_vector_view & getDataFieldEnts() const
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
auto & getColFieldEnts() const
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Provide data structure for (tensor) field approximation.
FieldSpace getSpace() const
Get field approximation space.
const BitFieldId & getId() const
Get unique field id.
MultiIndex Tag for field name.
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...
structure to get information form mofem into EntitiesFieldData
int getMaxRowOrder() const
Get max order of approximation for field in rows.
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getNoFieldRowIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
boost::ptr_deque< UserDataOperator > opPtrVector
Vector of finite element users data operators.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
int getMaxColOrder() const
Get max order of approximation for field in columns.
MoFEMErrorCode getNoFieldColIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode getNoFieldFieldData(const int bit_number, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
int normalShift
Shift in vector for linear geometry.
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:355
Calculate base functions on triangle.
Base volume element used to integrate on contact surface (could be extended to other volume elements)...
int npoints
Definition: quad.h:29