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