v0.13.2
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 auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
452 &m(0, 0), &m(0, 1), &m(0, 2));
453 };
454
455 auto get_tensor_vec = [](VectorDouble &n, const int r = 0) {
456 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1),
457 &n(r + 2));
458 };
459
460 CHKERR getUserPolynomialBase()->getValue(
462 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
464 NOBASE)));
465 CHKERR getUserPolynomialBase()->getValue(
467 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
468 dataHdivSlave, HDIV, static_cast<FieldApproximationBase>(b),
469 NOBASE)));
470
476
477 }
478
479 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
480 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
481 "Not yet implemented");
482 }
483 if (dataH1.spacesOnEntities[MBTET].test(L2)) {
484 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
485 "Not yet implemented");
486 }
487 break;
488 default:
489 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
490 "Not yet implemented");
491 }
492 }
493 }
494
495 // Iterate over operators
497
499}
500
503
504 const EntityType type = numeredEntFiniteElementPtr->getEntType();
505 constexpr std::array<UserDataOperator::OpType, 2> types{
506 UserDataOperator::OPROW, UserDataOperator::OPCOL};
507 std::array<std::string, 2> last_eval_field_name{std::string(), std::string()};
508
509 auto oit = opPtrVector.begin();
510 auto hi_oit = opPtrVector.end();
511
512 for (; oit != hi_oit; oit++) {
513
514 oit->setPtrFE(this);
515
516 if (oit->opType == UserDataOperator::OPSPACE) {
517
518 // Set field
519 switch (oit->sPace) {
520 case NOSPACE:
521 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Unknown space");
522 case NOFIELD:
523 case H1:
524 case HCURL:
525 case HDIV:
526 case L2:
527 break;
528 default:
529 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
530 "Not implemented for this space", oit->sPace);
531 }
532
533 // Reseat all data which all field dependent
534 dataOnMaster[oit->sPace]->resetFieldDependentData();
535 dataOnSlave[oit->sPace]->resetFieldDependentData();
536
537 // Run operator
538 try {
539 CHKERR oit->opRhs(*dataOnMaster[oit->sPace], false);
540 CHKERR oit->opRhs(*dataOnSlave[oit->sPace], false);
541 }
542 CATCH_OP_ERRORS(*oit);
543
544 } else {
545 boost::shared_ptr<EntitiesFieldData> op_master_data[2];
546 boost::shared_ptr<EntitiesFieldData> op_slave_data[2];
547
548 for (int ss = 0; ss != 2; ss++) {
549
550 const std::string field_name =
551 !ss ? oit->rowFieldName : oit->colFieldName;
552 const Field *field_struture = mField.get_field_structure(field_name);
553 const BitFieldId data_id = field_struture->getId();
554 const FieldSpace space = field_struture->getSpace();
555 op_master_data[ss] =
556 !ss ? dataOnMaster[space] : derivedDataOnMaster[space];
557 op_slave_data[ss] =
558 !ss ? dataOnSlave[space] : derivedDataOnSlave[space];
559
560 if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
561 data_id)
562 .none())
563 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
564 "no data field < %s > on finite element < %s >",
565 field_name.c_str(), getFEName().c_str());
566
567 if (oit->getOpType() & types[ss] ||
568 oit->getOpType() & UserDataOperator::OPROWCOL) {
569
570 switch (space) {
571 case NOSPACE:
572 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown space");
573 break;
574 case NOFIELD:
575 case H1:
576 case HCURL:
577 case HDIV:
578 case L2:
579 break;
580 default:
581 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
582 "Not implemented for this space", space);
583 }
584
585 if (last_eval_field_name[ss] != field_name) {
586 CHKERR getEntityFieldData(*op_master_data[ss], *op_slave_data[ss],
587 field_name, MBEDGE);
588
589 if (!ss) {
590 struct Extractor {
591 boost::weak_ptr<EntityCacheNumeredDofs>
592 operator()(boost::shared_ptr<FieldEntity> &e) {
593 return e->entityCacheRowDofs;
594 }
595 };
596
597 CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
598 field_name, getRowFieldEnts(), MBEDGE,
599 MBPRISM, Extractor());
600 } else {
601 struct Extractor {
602 boost::weak_ptr<EntityCacheNumeredDofs>
603 operator()(boost::shared_ptr<FieldEntity> &e) {
604 return e->entityCacheColDofs;
605 }
606 };
607 CHKERR getEntityIndices(*op_master_data[ss], *op_slave_data[ss],
608 field_name, getRowFieldEnts(), MBEDGE,
609 MBPRISM, Extractor());
610 }
611
612 switch (space) {
613 case NOSPACE:
614 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
615 "unknown space");
616 break;
617 case H1: {
618
619 auto get_indices = [&](auto &master, auto &slave, auto &ents,
620 auto &&ex) {
621 return getNodesIndices(
622 field_name, ents,
623 master.dataOnEntities[MBVERTEX][0].getIndices(),
624 master.dataOnEntities[MBVERTEX][0].getLocalIndices(),
625 slave.dataOnEntities[MBVERTEX][0].getIndices(),
626 slave.dataOnEntities[MBVERTEX][0].getLocalIndices(), ex);
627 };
628
629 auto get_data = [&](EntitiesFieldData &master_data,
630 EntitiesFieldData &slave_data) {
631 return getNodesFieldData(
633 master_data.dataOnEntities[MBVERTEX][0].getFieldData(),
634 slave_data.dataOnEntities[MBVERTEX][0].getFieldData(),
635 master_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
636 slave_data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
637 master_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
638 slave_data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
639 master_data.dataOnEntities[MBVERTEX][0].getSpace(),
640 slave_data.dataOnEntities[MBVERTEX][0].getSpace(),
641 master_data.dataOnEntities[MBVERTEX][0].getBase(),
642 slave_data.dataOnEntities[MBVERTEX][0].getBase());
643 };
644
645 if (!ss) {
646
647 struct Extractor {
648 boost::weak_ptr<EntityCacheNumeredDofs>
649 operator()(boost::shared_ptr<FieldEntity> &e) {
650 return e->entityCacheRowDofs;
651 }
652 };
653
654 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
655 getRowFieldEnts(), Extractor());
656 } else {
657 struct Extractor {
658 boost::weak_ptr<EntityCacheNumeredDofs>
659 operator()(boost::shared_ptr<FieldEntity> &e) {
660 return e->entityCacheColDofs;
661 }
662 };
663
664 CHKERR get_indices(*op_master_data[ss], *op_slave_data[ss],
665 getColFieldEnts(), Extractor());
666 }
667
668 CHKERR get_data(*op_master_data[ss], *op_slave_data[ss]);
669
670 } break;
671 case HCURL:
672 case HDIV:
673 case L2:
674 break;
675 case NOFIELD:
676 if (!getNinTheLoop()) {
677 // NOFIELD data are the same for each element, can be
678 // retrieved only once
679 const auto bit_number = mField.get_field_bit_number(field_name);
680 if (!ss) {
681 CHKERR getNoFieldRowIndices(*op_master_data[ss], bit_number);
682 CHKERR getNoFieldRowIndices(*op_slave_data[ss], bit_number);
683 } else {
684 CHKERR getNoFieldColIndices(*op_master_data[ss], bit_number);
685 CHKERR getNoFieldColIndices(*op_slave_data[ss], bit_number);
686 }
687 CHKERR getNoFieldFieldData(*op_master_data[ss], bit_number);
688 CHKERR getNoFieldFieldData(*op_slave_data[ss], bit_number);
689 }
690 break;
691 default:
692 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
693 "Not implemented for this space", space);
694 }
695 last_eval_field_name[ss] = field_name;
696 }
697 }
698 }
699
700 int type;
701
702 if (UserDataOperator *cast_oit =
703 dynamic_cast<UserDataOperator *>(&*oit)) {
704 type = cast_oit->getFaceType();
705 if (((oit->getOpType() & UserDataOperator::OPROW) ||
706 (oit->getOpType() & UserDataOperator::OPCOL)) &&
707 ((type & UserDataOperator::FACEMASTERMASTER) ||
708 (type & UserDataOperator::FACEMASTERSLAVE) ||
709 (type & UserDataOperator::FACESLAVEMASTER) ||
710 (type & UserDataOperator::FACESLAVESLAVE))) {
711 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
712 "Wrong combination of FaceType and OpType, OPROW or OPCOL "
713 "combined with face-face OpType");
714 }
715
716 if ((oit->getOpType() & UserDataOperator::OPROWCOL) &&
717 ((type & UserDataOperator::FACEMASTER) ||
718 (type & UserDataOperator::FACESLAVE))) {
719 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
720 "Wrong combination of FaceType and OpType, OPROWCOL "
721 "combined with face-face OpType");
722 }
723
724 if (!type) {
725 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
726 "Face type is not set");
727 }
728 } else {
729 type = UserDataOperator::FACEMASTER | UserDataOperator::FACESLAVE |
730 UserDataOperator::FACEMASTERMASTER |
731 UserDataOperator::FACEMASTERSLAVE |
732 UserDataOperator::FACESLAVEMASTER |
733 UserDataOperator::FACESLAVESLAVE;
734 }
735
736 if (oit->getOpType() & UserDataOperator::OPROW &&
737 (type & UserDataOperator::FACEMASTER)) {
738 try {
739 CHKERR oit->opRhs(*op_master_data[0], false);
740 }
741 CATCH_OP_ERRORS(*oit);
742 }
743
744 if (oit->getOpType() & UserDataOperator::OPROW &&
745 (type & UserDataOperator::FACESLAVE)) {
746 try {
747 CHKERR oit->opRhs(*op_slave_data[0], false);
748 }
749 CATCH_OP_ERRORS(*oit);
750 }
751
752 if (oit->getOpType() & UserDataOperator::OPCOL &&
753 (type & UserDataOperator::FACEMASTER)) {
754 try {
755 CHKERR oit->opRhs(*op_master_data[1], false);
756 }
757 CATCH_OP_ERRORS(*oit);
758 }
759
760 if (oit->getOpType() & UserDataOperator::OPCOL &&
761 (type & UserDataOperator::FACESLAVE)) {
762 try {
763 CHKERR oit->opRhs(*op_slave_data[1], false);
764 }
765 CATCH_OP_ERRORS(*oit);
766 }
767
768 if (oit->getOpType() & UserDataOperator::OPROWCOL &&
769 (type & UserDataOperator::FACEMASTERMASTER)) {
770 try {
771 CHKERR oit->opLhs(*op_master_data[0], *op_master_data[1]);
772 }
773 CATCH_OP_ERRORS(*oit);
774 }
775
776 if (oit->getOpType() & UserDataOperator::OPROWCOL &&
777 (type & UserDataOperator::FACEMASTERSLAVE)) {
778 try {
779 CHKERR oit->opLhs(*op_master_data[0], *op_slave_data[1]);
780 }
781 CATCH_OP_ERRORS(*oit);
782 }
783
784 if (oit->getOpType() & UserDataOperator::OPROWCOL &&
785 (type & UserDataOperator::FACESLAVEMASTER)) {
786 try {
787 CHKERR oit->opLhs(*op_slave_data[0], *op_master_data[1]);
788 }
789 CATCH_OP_ERRORS(*oit);
790 }
791
792 if (oit->getOpType() & UserDataOperator::OPROWCOL &&
793 (type & UserDataOperator::FACESLAVESLAVE)) {
794 try {
795 CHKERR oit->opLhs(*op_slave_data[0], *op_slave_data[1]);
796 }
797 CATCH_OP_ERRORS(*oit);
798 }
799 }
800 }
802}
803
805 EntitiesFieldData &master_data, EntitiesFieldData &slave_data,
806 const std::string &field_name, const EntityType type_lo,
807 const EntityType type_hi) const {
809
810 auto reset_data = [type_lo, type_hi](auto &data) {
811 for (EntityType t = type_lo; t != type_hi; ++t) {
812 for (auto &dat : data.dataOnEntities[t]) {
813 dat.getOrder() = 0;
814 dat.getBase() = NOBASE;
815 dat.getSpace() = NOSPACE;
816 dat.getFieldData().resize(0, false);
817 dat.getFieldDofs().resize(0, false);
818 dat.getFieldEntities().resize(0, false);
819 }
820 }
821 };
822 reset_data(master_data);
823 reset_data(slave_data);
824
825 auto &field_ents = getDataFieldEnts();
826 auto bit_number = mField.get_field_bit_number(field_name);
827 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
828 bit_number, get_id_for_min_type(type_lo));
829 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
830 cmp_uid_lo);
831 if (lo != field_ents.end()) {
832 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
833 bit_number, get_id_for_max_type(type_hi));
834 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
835 if (lo != hi) {
836 for (auto it = lo; it != hi; ++it)
837 if (auto e = it->lock()) {
838
839 auto get_data = [&](auto &data, auto type, auto side) {
840 auto &dat = data.dataOnEntities[type][side];
841 auto &ent_field_dofs = dat.getFieldDofs();
842 auto &ent_field_data = dat.getFieldData();
843 dat.getFieldEntities().resize(1, false);
844 dat.getFieldEntities()[0] = e.get();
845 dat.getBase() = e->getApproxBase();
846 dat.getSpace() = e->getSpace();
847 const int ent_order = e->getMaxOrder();
848 dat.getOrder() =
849 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
850 const auto dof_ent_field_data = e->getEntFieldData();
851 const int nb_dofs_on_ent = e->getNbDofsOnEnt();
852 ent_field_data.resize(nb_dofs_on_ent, false);
853 noalias(ent_field_data) = e->getEntFieldData();
854 ent_field_dofs.resize(nb_dofs_on_ent, false);
855 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
856 if (auto cache = e->entityCacheDataDofs.lock()) {
857 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
858 ent_field_dofs[(*dit)->getEntDofIdx()] =
859 reinterpret_cast<FEDofEntity *>((*dit).get());
860 }
861 }
862 };
863
864 const EntityType type = e->getEntType();
865 const int side = e->getSideNumberPtr()->side_number;
866
867 switch (type) {
868 case MBEDGE:
869
870 if (side < 3)
871 get_data(master_data, type, side);
872 else if (side > 5)
873 get_data(slave_data, type, side - 6);
874
875 break;
876 case MBTRI:
877
878 if (side == 3)
879 get_data(master_data, type, 0);
880 if (side == 4)
881 get_data(slave_data, type, 0);
882
883 break;
884 default:
885 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
886 "Entity type not implemented (FIXME)");
887 };
888
889 const int brother_side = e->getSideNumberPtr()->brother_side_number;
890 if (brother_side != -1)
891 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
892 "Case with brother side not implemented (FIXME)");
893 }
894 }
895 }
896
898}
899
901 const std::string field_name, VectorDouble &master_nodes_data,
902 VectorDouble &slave_nodes_data, VectorDofs &master_nodes_dofs,
903 VectorDofs &slave_nodes_dofs,
904
905 VectorFieldEntities &master_field_entities, VectorFieldEntities &slave_field_entities,
906
907 FieldSpace &master_space,
908 FieldSpace &slave_space, FieldApproximationBase &master_base,
909 FieldApproximationBase &slave_base) const {
911
912 auto set_zero = [](auto &nodes_data, auto &nodes_dofs, auto &field_entities) {
913 nodes_data.resize(0, false);
914 nodes_dofs.resize(0, false);
915 field_entities.resize(0, false);
916 };
917 set_zero(master_nodes_data, master_nodes_dofs, master_field_entities);
918 set_zero(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
919
920 auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
921 if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
922
923 auto bit_number = (*field_it)->getBitNumber();
924 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
925 master_space = slave_space = (*field_it)->getSpace();
926 master_base = slave_base = (*field_it)->getApproxBase();
927
928 auto &field_ents = getDataFieldEnts();
929 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
930 bit_number, get_id_for_min_type<MBVERTEX>());
931 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
932 cmp_uid_lo);
933 if (lo != field_ents.end()) {
934 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
935 bit_number, get_id_for_max_type<MBVERTEX>());
936 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
937 if (lo != hi) {
938
939 int nb_dofs = 0;
940 for (auto it = lo; it != hi; ++it) {
941 if (auto e = it->lock()) {
942 nb_dofs += e->getEntFieldData().size();
943 }
944 }
945
946 if (nb_dofs) {
947
948 auto init_set = [&](auto &nodes_data, auto &nodes_dofs,
949 auto &field_entities) {
950 constexpr int num_nodes = 3;
951 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
952 nodes_data.resize(max_nb_dofs, false);
953 nodes_dofs.resize(max_nb_dofs, false);
954 field_entities.resize(num_nodes, false);
955 nodes_data.clear();
956 fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
957 fill(field_entities.begin(), field_entities.end(), nullptr);
958 };
959
960 init_set(master_nodes_data, master_nodes_dofs, master_field_entities);
961 init_set(slave_nodes_data, slave_nodes_dofs, slave_field_entities);
962
963 for (auto it = lo; it != hi; ++it) {
964 if (auto e = it->lock()) {
965
966 const auto &sn = e->getSideNumberPtr();
967 int side = sn->side_number;
968
969 auto set_data = [&](auto &nodes_data, auto &nodes_dofs,
970 auto &field_entities, int side, int pos) {
971 field_entities[side] = e.get();
972 if (auto cache = e->entityCacheDataDofs.lock()) {
973 for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
974 ++dit) {
975 const auto dof_idx = (*dit)->getEntDofIdx();
976 nodes_data[pos + dof_idx] = (*dit)->getFieldData();
977 nodes_dofs[pos + dof_idx] =
978 reinterpret_cast<FEDofEntity *>((*dit).get());
979 }
980 }
981 };
982
983 if (side < 3)
984 set_data(master_nodes_data, master_nodes_dofs,
985 master_field_entities, side, side * nb_dofs_on_vert);
986 else
987 set_data(slave_nodes_data, slave_nodes_dofs,
988 slave_field_entities, (side - 3),
989 (side - 3) * nb_dofs_on_vert);
990
991 const int brother_side = sn->brother_side_number;
992 if (brother_side != -1)
993 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
994 "Not implemented (FIXME please)");
995 }
996 }
997 }
998 }
999 }
1000 }
1001
1003}
1004
1005template <typename EXTRACTOR>
1007 EntitiesFieldData &master_data, EntitiesFieldData &slave_data,
1008 const std::string &field_name, FieldEntity_vector_view &ents_field,
1009 const EntityType type_lo, const EntityType type_hi,
1010 EXTRACTOR &&extractor) const {
1012
1013 auto clear_data = [type_lo, type_hi](auto &data) {
1014 for (EntityType t = type_lo; t != type_hi; ++t) {
1015 for (auto &d : data.dataOnEntities[t]) {
1016 d.getIndices().resize(0, false);
1017 d.getLocalIndices().resize(0, false);
1018 }
1019 }
1020 };
1021
1022 clear_data(master_data);
1023 clear_data(slave_data);
1024
1025 auto bit_number = mField.get_field_bit_number(field_name);
1026 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1027 bit_number, get_id_for_min_type(type_lo));
1028 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1029 cmp_uid_lo);
1030 if (lo != ents_field.end()) {
1031 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1032 bit_number, get_id_for_max_type(type_hi));
1033 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1034 if (lo != hi) {
1035
1036 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1037
1038 for (auto it = lo; it != hi; ++it)
1039 if (auto e = it->lock()) {
1040
1041 const EntityType type = e->getEntType();
1042 const int side = e->getSideNumberPtr()->side_number;
1043 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1044 if (brother_side != -1)
1045 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1046 "Not implemented case");
1047
1048 auto get_indices = [&](auto &data, const auto type, const auto side) {
1049 if (auto cache = extractor(e).lock()) {
1050 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
1051 auto &dof = (**dit);
1052 auto &dat = data.dataOnEntities[type][side];
1053 auto &ent_field_indices = dat.getIndices();
1054 auto &ent_field_local_indices = dat.getLocalIndices();
1055 if (ent_field_indices.empty()) {
1056 const int nb_dofs_on_ent = dof.getNbDofsOnEnt();
1057 ent_field_indices.resize(nb_dofs_on_ent, false);
1058 ent_field_local_indices.resize(nb_dofs_on_ent, false);
1059 std::fill(ent_field_indices.data().begin(),
1060 ent_field_indices.data().end(), -1);
1061 std::fill(ent_field_local_indices.data().begin(),
1062 ent_field_local_indices.data().end(), -1);
1063 }
1064 const int idx = dof.getEntDofIdx();
1065 ent_field_indices[idx] = dof.getPetscGlobalDofIdx();
1066 ent_field_local_indices[idx] = dof.getPetscLocalDofIdx();
1067 }
1068 }
1069 };
1070
1071 switch (type) {
1072 case MBEDGE:
1073
1074 if (side < 3)
1075 get_indices(master_data, type, side);
1076 else if (side > 5)
1077 get_indices(slave_data, type, side - 6);
1078
1079 break;
1080 case MBTRI:
1081
1082 if (side == 3)
1083 get_indices(master_data, type, 0);
1084 if (side == 4)
1085 get_indices(slave_data, type, 0);
1086
1087 break;
1088 default:
1089 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1090 "Entity type not implemented");
1091 }
1092 }
1093 }
1094 }
1095
1097}
1098
1099template <typename EXTRACTOR>
1101 const std::string field_name, FieldEntity_vector_view &ents_field,
1102 VectorInt &master_nodes_indices, VectorInt &master_local_nodes_indices,
1103 VectorInt &slave_nodes_indices, VectorInt &slave_local_nodes_indices,
1104 EXTRACTOR &&extractor) const {
1106
1107 master_nodes_indices.resize(0, false);
1108 master_local_nodes_indices.resize(0, false);
1109 slave_nodes_indices.resize(0, false);
1110 slave_local_nodes_indices.resize(0, false);
1111
1112 auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
1113 if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
1114
1115 auto bit_number = (*field_it)->getBitNumber();
1116 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
1117
1118 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1119 bit_number, get_id_for_min_type<MBVERTEX>());
1120 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
1121 cmp_uid_lo);
1122 if (lo != ents_field.end()) {
1123 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1124 bit_number, get_id_for_max_type<MBVERTEX>());
1125 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
1126 if (lo != hi) {
1127
1128 int nb_dofs = 0;
1129 for (auto it = lo; it != hi; ++it) {
1130 if (auto e = it->lock()) {
1131 if (auto cache = extractor(e).lock()) {
1132 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
1133 }
1134 }
1135 }
1136
1137 if (nb_dofs) {
1138
1139 constexpr int num_nodes = 3;
1140 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
1141
1142 auto set_vec_size = [&](auto &nodes_indices,
1143 auto &local_nodes_indices) {
1144 nodes_indices.resize(max_nb_dofs, false);
1145 local_nodes_indices.resize(max_nb_dofs, false);
1146 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
1147 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(),
1148 -1);
1149 };
1150
1151 set_vec_size(master_nodes_indices, master_local_nodes_indices);
1152 set_vec_size(slave_nodes_indices, slave_local_nodes_indices);
1153
1154 for (auto it = lo; it != hi; ++it) {
1155 if (auto e = it->lock()) {
1156
1157 const int side = e->getSideNumberPtr()->side_number;
1158 const int brother_side =
1159 e->getSideNumberPtr()->brother_side_number;
1160 if (brother_side != -1)
1161 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1162 "Not implemented case");
1163
1164 auto get_indices = [&](auto &nodes_indices,
1165 auto &local_nodes_indices,
1166 const auto side) {
1167 if (auto cache = extractor(e).lock()) {
1168 for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
1169 ++dit) {
1170 const int idx = (*dit)->getPetscGlobalDofIdx();
1171 const int local_idx = (*dit)->getPetscLocalDofIdx();
1172 const int pos =
1173 side * nb_dofs_on_vert + (*dit)->getDofCoeffIdx();
1174 nodes_indices[pos] = idx;
1175 local_nodes_indices[pos] = local_idx;
1176 }
1177 }
1178 };
1179
1180 if (side < 3)
1181 get_indices(master_nodes_indices, master_local_nodes_indices,
1182 side);
1183 else if (side > 2)
1184 get_indices(slave_nodes_indices, slave_local_nodes_indices,
1185 side - 3);
1186 else
1187 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1188 "Impossible case");
1189 }
1190 }
1191 }
1192 }
1193 }
1194 }
1195
1197}
1198
1200ContactPrismElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes(
1201 const string fe_name,
1203 const int side_type, const EntityHandle ent_for_side) {
1204 return loopSide(fe_name, &fe_method, side_type, ent_for_side);
1205}
1206
1207} // 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
const int dim
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
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)
User call this function to loop over elements on the side of face. This function calls finite element...
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)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:352
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